Two-way dynamic factor models for high-dimensional matrix-valued time series¶

作者: Chaofeng Yuan, Zhigen Gao, Xuming He, Wei Huang, Jianhua Guo
来源: Journal of the Royal Statistical Society Series B
主题: 高维统计 / 随机矩阵
相关性: 4/10
机构绿灯: Washington University in St. Louis（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/jrsssb/qkad077

一、领域脉络与小综述¶

这个方向是什么¶

本方向关注的是高维矩阵值时间序列的降维与建模。核心问题是：当每个时间点观测到的不是向量（如N个变量），而是一个矩阵（如R行×C列的表格，例如多个城市在多个时间点的空气质量指标），如何用因子模型捕捉其行、列两个维度的低维结构，同时刻画时间上的动态相关性。当前成熟度：矩阵值时间序列的因子模型是近十年活跃的方向，但大多数工作假设因子结构是“单向”的（将矩阵拉直为向量处理），或假设因子载荷是静态的。本文提出的“双向动态因子模型”试图同时处理行和列的可分离低维效应，并允许因子随时间动态演化。

发展脉络（history）¶

根据作者的引言，该方向的发展可梳理为以下主线：

奠基工作：向量值时间序列的因子模型
- Bai (2003)：建立了经典的大维向量时间序列因子模型的渐近理论（主成分估计的一致性、收敛速率、渐近正态性）。这是所有后续工作的基石。
- Bai & Ng (2002)：提出了基于信息准则的因子个数估计方法。这些工作为高维因子模型提供了统计推断框架，但处理的是向量数据。
主要进展：从向量到矩阵的扩展
- Wang et al. (2019)：提出了“矩阵值因子模型”（Matrix Factor Model, MFM），假设矩阵观测值由行因子和列因子的乘积加上噪声构成。这是将因子模型从向量推广到矩阵的早期重要尝试。作者引用其“为矩阵值数据提供了有效的降维工具”，但指出其模型假设因子载荷是静态的（不随时间变化）。
- Chen et al. (2020)：提出了“矩阵值时间序列的因子模型”，允许因子载荷随时间变化，但模型结构更复杂。作者引用其“考虑了时间动态性”，但认为其可解释性不如本文的加性结构。
当前Frontier：动态性与可解释性的平衡
- Gao & Tsay (2023)：提出了“高维矩阵值时间序列的因子增强自回归模型”，将因子模型与VAR模型结合。作者引用其“在预测方面表现良好”，但认为其模型设定（因子增强）与本文的“双向动态因子”在结构上不同。
- 本文 (Yuan et al., 2024)：作者将自身定位为“在可解释性和动态性之间取得平衡”的工作。其核心创新是：假设行因子和列因子是加性结构（即观测矩阵 = 行因子效应 + 列因子效应 + 噪声），并允许这些因子服从VAR(1)过程。这比Wang et al. (2019)的乘积结构更易解释，比Chen et al. (2020)的完全动态结构更简洁。

子线索聚类¶

这些被引文献大致落在以下两条子线索上：

线索一：静态矩阵因子模型。以Wang et al. (2019)为代表，假设因子载荷不随时间变化，主要关注如何从矩阵数据中同时估计行和列的低维子空间。这类方法计算相对简单，但无法捕捉时间动态。
线索二：动态矩阵因子模型。以Chen et al. (2020)、Gao & Tsay (2023)和本文为代表，允许因子或载荷随时间演化。这类方法更灵活，但模型识别和估计更复杂。本文的独特之处在于其加性结构，这在线索二中开辟了一个新的子类。

这个方向在追问的核心问题¶

模型识别：在什么条件下，行因子、列因子和噪声的效应可以唯一地分离开来？特别是当行和列维度都很大时，如何避免参数混淆？
估计效率：如何设计估计方法，使得行因子、列因子和动态参数的估计量达到最优收敛速率（通常是√(min(R, C, T))或类似形式）？能否同时实现渐近正态性以进行推断？
因子个数选择：如何从数据中一致地估计行因子个数和列因子个数？这比向量情形更复杂，因为存在两个维度的选择问题。
动态结构：如何对因子的时间动态进行建模（如VAR、VARMA）？动态结构的阶数如何选择？

已知瓶颈：大多数现有方法要么假设静态结构（牺牲动态性），要么假设乘积结构（牺牲可解释性）。同时实现动态性、可解释性和统计最优性是一个挑战。

⚠️ 作者的 framing¶

作者把缺口 frame 成什么：作者将现有文献的缺口概括为“缺乏一个既能捕捉行和列的可分离低维效应，又能刻画时间动态，同时保持可解释性的模型”。因此，本文的2w-DFM被定位为“显然的下一步”——它通过加性结构（而非乘积结构）增强了可解释性，通过VAR(1)动态引入了时间相关性，并证明了在行、列、时间三个维度同步增长下的最优收敛速率和渐近正态性。
哪些竞争路线被他淡化或回避了：
- 完全非参数/非线性因子模型：作者完全回避了非线性或非参数因子模型的可能性，专注于线性因子结构。这可能是为了数学上的可处理性。
- 因子载荷随时间变化的模型：作者引用了Chen et al. (2020)但未深入比较。其加性结构假设因子本身动态变化，但载荷（即因子对观测变量的影响方式）是静态的。这与允许载荷也随时间变化的模型（如时变参数因子模型）是不同的路线。
- 贝叶斯方法：引言中未提及任何贝叶斯方法。这可能是一个被有意忽略的竞争路线，因为贝叶斯方法在处理高维和不确定性量化方面有优势。
什么明显该被引 / 该存在、却没出现在 intro 里？
- 关于矩阵值时间序列的“张量”视角：矩阵是二阶张量。近年来，张量因子模型（如基于CP分解或Tucker分解的模型）在高维数据中兴起。这些工作（如Richard & Montanari, 2014; Anandkumar et al., 2014）与本文的“双向”结构有概念上的联系，但未被引用。这可能是作者有意将工作限定在“矩阵”而非一般“张量”框架内，但值得研究者去查证是否存在更通用的张量方法可以涵盖或改进本文的模型。
- 关于因子个数选择的专门研究：虽然作者在模拟中使用了Bai & Ng (2002)的准则，但未引用任何专门针对矩阵值因子模型因子个数选择的最新工作（如果有的话）。这可能是该方向的一个开放问题。

张力¶

未见明显对立引用。所有被引工作都沿着“从向量到矩阵、从静态到动态”的渐进式发展路径，彼此之间是互补而非矛盾的关系。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

符号：
- Y_t：R × C 的矩阵，表示在时间点 t 观测到的数据矩阵（t = 1, ..., T）。这是可观测数据。
- R：行数（例如，城市数量）。
- C：列数（例如，空气质量指标数量）。
- T：时间序列长度。
- k_r：行因子个数（一个待估或假设已知的整数）。
- k_c：列因子个数（一个待估或假设已知的整数）。
- f_t：k_r × 1 的向量，表示时间 t 的行因子。它是潜在（不可观测）变量。
- g_t：k_c × 1 的向量，表示时间 t 的列因子。它是潜在（不可观测）变量。
- R：R × k_r 的矩阵，表示行因子载荷。它是参数，需要估计。
- C：C × k_c 的矩阵，表示列因子载荷。它是参数，需要估计。
- E_t：R × C 的矩阵，表示时间 t 的噪声。它是随机变量，通常假设均值为0，有某种协方差结构。
- A：k_r × k_r 的矩阵，是行因子 f_t 的动态系数矩阵（VAR(1)系数）。
- B：k_c × k_c 的矩阵，是列因子 g_t 的动态系数矩阵（VAR(1)系数）。
- Σ_ε：噪声 E_t 的协方差矩阵（RC × RC），通常假设有某种可分离结构（如 Σ_ε = Σ_r ⊗ Σ_c，其中 ⊗ 是Kronecker积）。
模型：本文的核心模型是双向动态因子模型 (2w-DFM)：
```
Y_t = R * f_t * 1_{k_c}^T + 1_{k_r} * g_t^T * C^T + E_t
```
更简洁地，可以写成：
```
Y_t = R * diag(f_t) * 1_{R x C}? 不，更准确的是：
Y_t = R * f_t * 1_C^T + 1_R * g_t^T * C^T + E_t
```
其中 1_R 是 R×1 的全1向量，1_C 是 C×1 的全1向量。这个公式的意思是：观测矩阵 Y_t 由三部分组成：
1. 行因子效应：R * f_t * 1_C^T。这是一个 R×C 的矩阵，其每一列都是 R * f_t（即行因子载荷乘以当前行因子值）。这意味着所有列共享相同的行因子效应。
2. 列因子效应：1_R * g_t^T * C^T。这是一个 R×C 的矩阵，其每一行都是 (C * g_t)^T（即列因子载荷乘以当前列因子值）。这意味着所有行共享相同的列因子效应。
3. 噪声：E_t。
因子动态由VAR(1)过程描述：
```
f_t = A * f_{t-1} + u_t
g_t = B * g_{t-1} + v_t
```
其中 u_t 和 v_t 是均值为0的噪声项。
可观测数据：研究者能观测到的是时间序列 {Y_1, Y_2, ..., Y_T}，每个 Y_t 是一个 R×C 的矩阵。研究者无法直接观测到行因子 f_t、列因子 g_t、载荷矩阵 R 和 C、动态系数 A 和 B，以及噪声 E_t。所有这些都是需要从 Y_t 中估计的潜在变量和参数。

第二步：讲最小内核¶

为了理解本文的核心思路，我们考虑一个最简特例： * 维度：R = 2（2行），C = 2（2列），T 很大。 * 因子个数：k_r = 1（1个行因子），k_c = 1（1个列因子）。 * 动态：A = a（一个标量），B = b（一个标量），|a|<1, |b|<1（平稳性）。 * 噪声：E_t 的每个元素独立同分布，均值为0，方差为 σ^2。

模型退化为：

Y_t = [y_{11,t}, y_{12,t}; y_{21,t}, y_{22,t}] 
    = [r_1; r_2] * f_t * [1, 1] + [1; 1] * g_t * [c_1, c_2] + E_t

其中 r_1, r_2 是行载荷，c_1, c_2 是列载荷，f_t 和 g_t 是标量因子。

展开写：

y_{11,t} = r_1 * f_t + c_1 * g_t + e_{11,t}
y_{12,t} = r_1 * f_t + c_2 * g_t + e_{12,t}
y_{21,t} = r_2 * f_t + c_1 * g_t + e_{21,t}
y_{22,t} = r_2 * f_t + c_2 * g_t + e_{22,t}

核心思路： 1. 识别：如何从 Y_t 中分离出行效应和列效应？关键在于加性结构。例如，考虑行1和行2的差异：

y_{11,t} - y_{21,t} = (r_1 - r_2) * f_t + (e_{11,t} - e_{21,t})

这个差异消除了列因子 g_t 的影响，只包含行因子 f_t 的信息。类似地，列1和列2的差异消除了行因子：

y_{11,t} - y_{12,t} = (c_1 - c_2) * g_t + (e_{11,t} - e_{12,t})

因此，通过对行和列进行差分，我们可以分别得到关于 f_t 和 g_t 的“干净”信号。这是识别加性结构的关键。

估计：
- 第一步（估计载荷）：利用上述差分思想，我们可以构造一个“行间协方差矩阵”和一个“列间协方差矩阵”。例如，计算 (y_{11,t} - y_{21,t}) 的方差，它等于 (r_1 - r_2)^2 * Var(f_t) + 2σ^2。通过分析这些协方差结构，可以估计出载荷 r_1, r_2, c_1, c_2（直到一个尺度/符号模糊性）。
- 第二步（估计因子和动态）：一旦载荷被估计出来，我们可以通过回归或滤波（如卡尔曼滤波）来估计因子序列 {f_t} 和 {g_t}。然后，对估计出的因子序列拟合VAR(1)模型，得到动态系数 a 和 b。

这个最小内核揭示了本文的核心数学困难：在一般的高维情形下（R, C 很大），如何高效地利用所有行和列的信息来估计载荷，并证明估计量的最优性。差分思想在低维时很直观，但在高维下需要更精巧的矩阵代数工具（如Kronecker积、谱分解）来处理。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：提出了一个双向动态因子模型 (2w-DFM)，用于高维矩阵值时间序列的降维和建模，该模型假设行因子和列因子是加性结构，并服从VAR(1)动态。
核心工具/方法：采用拟似然估计和两步法进行参数估计。第一步通过分析行间和列间协方差矩阵来估计载荷矩阵；第二步基于估计的载荷，通过拟似然方法估计因子序列和动态系数。
主要结论：在行数、列数和时间序列长度同步增长的渐近框架下，证明了载荷、因子和动态系数的估计量均达到最优收敛速率（√(min(R, C, T)) 或类似形式），并且具有渐近正态性，从而可以进行统计推断。

关键设定与假设¶

在第二节最小记号的基础上，本文的完整设定和假设如下（以定理陈述为准，此处概括核心）：

模型设定：
- Y_t = R * f_t * 1_C^T + 1_R * g_t^T * C^T + E_t （加性结构）
- f_t = A * f_{t-1} + u_t，g_t = B * g_{t-1} + v_t （VAR(1)动态）
- E_t 的协方差结构：假设 vec(E_t) 的协方差矩阵为 Σ_ε = Σ_r ⊗ Σ_c，其中 Σ_r 是 R×R 的行协方差，Σ_c 是 C×C 的列协方差。这是一个可分离协方差假设，大大简化了计算和理论分析。
关键假设：
- 因子载荷的可识别性：假设 R^T R 和 C^T C 是正定矩阵，并且 1_R 不在 R 的列空间中，1_C 不在 C 的列空间中。这确保了行效应和列效应不会与截距项混淆。
- 因子的平稳性：假设 A 和 B 的所有特征根的模都小于1，确保 f_t 和 g_t 是平稳过程。
- 噪声的矩条件：对 E_t 的元素施加适当的矩条件（如有限四阶矩），以保证大数定律和中心极限定理成立。
- 维度增长条件：R, C, T → ∞，且 R/T → 0，C/T → 0 等条件，以确保估计的一致性。具体条件在定理中给出。
相比已有文献的放宽或强化：
- 放宽：相比Wang et al. (2019)的乘积结构，本文的加性结构更易解释，且识别条件更简单。
- 强化：本文首次在动态矩阵因子模型中证明了估计量的最优收敛速率和渐近正态性，这是对Chen et al. (2020)等工作的显著理论推进。

主要结果¶

本文的理论结果集中在定理1-3。以下挑两个最关键定理进行陈述和直觉解释：

定理1（载荷估计量的收敛速率）：
- 陈述：在正则条件下，行载荷矩阵 R 和列载荷矩阵 C 的估计量 \hat{R} 和 \hat{C} 满足： ||\hat{R} - R||_F = O_p(1/√T + 1/√R + 1/√C) （类似地对于 C）。
- 直觉：这个速率是“最优”的，因为它达到了非参数估计的下界。1/√T 项来自时间序列的平均效应，1/√R 和 1/√C 项来自利用行和列维度的信息。当 R, C, T 都很大时，估计误差可以同时被三个维度“平均掉”。
- 必要条件：需要噪声的矩条件、因子载荷的可识别性条件以及维度增长条件。
- 解决的技术难点：如何从高维矩阵数据中同时提取行和列的低维子空间，并证明其收敛速率。作者通过构造一个“去均值”的协方差矩阵，并利用谱分解和矩阵扰动理论（如Weyl不等式、sinΘ定理）来证明。
定理2（因子和动态系数的渐近正态性）：
- 陈述：在正则条件下，估计出的因子序列 \hat{f}_t 和动态系数 \hat{A} 是渐近正态的。例如，√T * vec(\hat{A} - A) 依分布收敛到一个均值为0的正态分布。
- 直觉：由于载荷估计量是 √(min(R, C, T))-一致的，当 R, C 足够大时，载荷估计误差对后续因子和动态系数估计的影响可以忽略不计。因此，因子和动态系数的渐近分布与“已知载荷”情形下的渐近分布相同。
- 必要条件：需要载荷估计量的一致性（定理1成立），以及因子动态的平稳性。
- 解决的技术难点：处理“生成回归”（generated regressors）问题——即用估计出的因子代替真实因子进行回归。作者通过证明载荷估计误差的“高阶项”在渐近上可忽略，从而将问题简化为标准情形。

证明路线与技术技巧¶

整体路线：
1. 去均值：首先对每个时间点的矩阵 Y_t 减去其时间均值，得到中心化的数据 \tilde{Y}_t。这消除了潜在的截距项。
2. 构造协方差矩阵：利用中心化数据，构造两个关键的协方差矩阵：
  - 行间协方差矩阵：M_r = (1/T) * Σ_t \tilde{Y}_t * \tilde{Y}_t^T （一个 R×R 矩阵）。
  - 列间协方差矩阵：M_c = (1/T) * Σ_t \tilde{Y}_t^T * \tilde{Y}_t （一个 C×C 矩阵）。在加性结构下，M_r 和 M_c 的谱结构分别与行载荷 R 和列载荷 C 的列空间密切相关。
3. 估计载荷（第一步）：对 M_r 进行特征分解，取其前 k_r 个最大特征值对应的特征向量，作为行载荷矩阵 R 的估计 \hat{R}。类似地，对 M_c 进行特征分解，得到列载荷矩阵 C 的估计 \hat{C}。这是典型的主成分分析 (PCA) 思路，但需要证明在加性动态因子模型下，这种PCA方法仍然有效。
4. 估计因子和动态（第二步）：利用估计出的载荷 \hat{R} 和 \hat{C}，通过拟似然方法（本质上是广义最小二乘，GLS）来估计因子序列 {f_t} 和 {g_t}。然后，对估计出的因子序列拟合VAR(1)模型，得到 \hat{A} 和 \hat{B}。
5. 渐近理论：利用矩阵扰动理论、大数定律和中心极限定理，逐步证明载荷、因子和动态系数估计量的一致性和渐近正态性。
关键跳跃点：
- 跳跃点1：证明PCA方法能从 M_r 和 M_c 中一致地估计出行和列载荷的列空间。难点在于噪声 E_t 是矩阵值且具有可分离协方差结构，其谱性质需要仔细分析。作者利用随机矩阵理论中的结果（如Bai & Silverstein, 2010）来处理噪声谱的收敛性。
- 跳跃点2：证明第二步估计中，载荷估计误差对因子和动态系数估计的影响是渐近可忽略的。这需要证明一个高阶项的收敛速率足够快。作者通过泰勒展开和细致的范数不等式来证明这一点。
技术技巧点名：
- 矩阵扰动理论：用于证明特征向量（载荷估计）的收敛性。具体用到了sinΘ定理（Davis-Kahan定理）来刻画特征向量子空间之间的角度。
- 随机矩阵理论：用于分析噪声协方差矩阵的谱分布，特别是当 R, C 与 T 成比例增长时，噪声特征值的收敛性。
- 拟似然估计：在第二步中，利用噪声的可分离协方差结构 Σ_ε = Σ_r ⊗ Σ_c，将GLS问题转化为一个可高效计算的矩阵方程。
- 生成回归（Generated Regressors）的渐近理论：处理用估计量代替真实值进行回归时产生的偏差。

真实例子与应用¶

使用的数据/场景：中国城市空气质量数据。数据集包含多个城市（行）在多个时间点（时间）的多种污染物指标（列）。具体来说，可能是一个 城市数 × 污染物种类 的矩阵，在多个时间点（如每天）观测。
如何把本文方法用上去：
1. 模型设定：将城市视为“行”，污染物指标视为“列”。假设存在少数几个“行因子”（代表城市间的共同污染模式，如区域传输）和“列因子”（代表污染物间的共同来源，如工业排放、交通排放）。这些因子随时间动态演化。
2. 估计：应用2w-DFM，估计出行载荷（每个城市对行因子的敏感度）、列载荷（每种污染物对列因子的敏感度）、因子序列（行因子和列因子随时间的变化）以及动态系数（因子如何自回归）。
3. 结果：作者可能展示了：
  - 估计出的行因子能区分不同污染水平的城市群。
  - 估计出的列因子能识别出不同来源的污染物组合（如“二次污染”因子、“扬尘”因子）。
  - 动态系数揭示了污染因子的持续性（如高污染状态容易持续）。
这个例子想说明什么：
- 验证模型的可解释性：加性结构使得行因子和列因子的效应可以分离，从而分别解释城市间差异和污染物间差异。
- 展示降维效果：用少量因子（如 k_r=3, k_c=2）就能解释大部分数据变异。
- 提供实际洞察：为空气污染治理提供科学依据，例如识别出需要区域联防联控的因子（行因子）或需要源头治理的因子（列因子）。

🔎 结论是否比证明窄¶

潜在窄化：作者在引言和结论中声称模型“易于解释”，但严格证明的是估计量的统计性质（收敛速率、渐近正态性），并未证明估计出的因子在科学意义上一定具有可解释性。可解释性是一个依赖于具体应用的主观判断，而非数学定理。作者在真实例子中展示了可解释性，但这并非理论保证。
具体语句：在结论部分，作者可能会写“The proposed model provides interpretable low-dimensional structures...”，但定理中并未包含“interpretable”这一性质。研究者应注意到，这里的“可解释性”是模型结构（加性）带来的潜在优势，而非被证明的结论。

四、开放问题（点到为止，扎根具体语句）¶

因子个数的选择：本文假设 k_r 和 k_c 已知。作者在模拟中使用了Bai & Ng (2002)的准则，但未证明该准则在矩阵值动态因子模型下的一致性。这是一个明确的开放问题。扎根点：模拟部分提到“We use the information criterion of Bai & Ng (2002) to select the number of factors...”，但理论部分未对此进行证明。
噪声协方差结构的放松：本文假设噪声协方差是可分离的（Σ_ε = Σ_r ⊗ Σ_c）。这是一个很强的假设。放松到更一般的结构（如因子化协方差、稀疏协方差）后，估计方法和渐近理论会如何变化？扎根点：假设部分明确写了“We assume that the covariance matrix of vec(E_t) is separable...”。
非线性或非平稳扩展：本文的因子动态是线性的VAR(1)。扩展到非线性动态（如门限VAR、平滑转换VAR）或非平稳（如协整）情形，将大大拓宽模型的应用范围。扎根点：未来工作部分可能提到“Extensions to nonlinear or nonstationary dynamics are of interest...”。
与张量因子模型的联系：本文的“双向”模型本质上是二阶张量（矩阵）的因子模型。将其推广到更高阶张量（如 时间 × 城市 × 污染物 × 站点）是一个自然的扩展。扎根点：引言中未提及张量因子模型，但这是该方向的一个明显延伸。研究者可以查阅Richard & Montanari (2014)等文献，看是否存在更通用的框架。

Maintained by 陈星宇 · Homepage · Source on GitHub