时间：2026-01-09 分类：数学
论文《矩阵值时间序列自回归模型的稳健估计》发表在《应用数学学报》，版权归《应用数学学报》所有。本文来自网络平台，仅供参考。
摘要
为解决矩阵自回归模型(MAR)的最小二乘估计易受异常值或厚尾分布误差影响而偏离真值的问题，本文对损失函数进行改进，基于投影法和迭代最小二乘法提出RE-PROJ和RE-ILS方法求解MAR模型的稳健估计。两种方法在构造估计量的过程中均借鉴了M估计的思想，模拟数据表明上述两种方法能够更好地抵御异常值对参数估计的影响。本文进一步讨论了MAR(p)模型阶数选择的BIC准则以及矩阵观测值存在相关性时带有减秩假定的稳健估计方法。实际数据分析结果表明，在异常值较多或尖峰厚尾的矩阵值观测数据中，本文提出的稳健估计相比于最小二乘估计在拟合和预测方面更具优势。
关键词
矩阵自回归模型;稳健估计;阶数选择;减秩
1 引言
随着数据科学的快速发展，矩阵值时间序列在金融、经济、工程等多个领域中日益常见，例如3D视频数据、多变量传感器监测数据、高维金融市场数据等。与传统的向量时间序列不同，矩阵值时间序列的行和列均具有实际意义，是结构化的数据形式[1-3]。如何对这类数据进行有效的统计建模与分析，提取其中的关键信息，成为当前研究的热点问题。
早期的研究中，Hannan[1]、Engle和Granger[2]、Tsay[3]等学者围绕向量自回归(VAR)模型展开了大量研究，为多变量时间序列分析提供了基础框架。然而，VAR模型存在维度灾难问题：当矩阵值时间序列被向量化后，即便是最简单的VAR(1)模型，其系数矩阵维度也会达到(mn×mn)(其中(m)和(n)分别为矩阵的行和列维度)，导致参数估计难度大、解释性差[4-6]。
为解决这一问题，Zhou和Li[4]、Zhao和Leng[5]、Wang[6]等学者尝试通过正则化等方法简化模型，但效果有限。Chen等[7]提出了矩阵自回归(MAR)模型，该模型通过Kronecker积结构约束系数矩阵，大幅减少了待估参数数量，为矩阵值时间序列的建模提供了新的思路。Chen等[7]针对MAR(1)模型提出了投影法(PROJ)和迭代最小二乘法(ILS)两种估计方法，通过目标函数求导给出了系数矩阵的迭代求解公式。但这两种方法均基于(L_2)损失函数，在存在异常值或误差服从厚尾分布的场景下，估计结果易受干扰而偏离真值[8-10]。
稳健估计方法能够有效抵御异常值和厚尾误差的影响，常见的包括Huber估计[10]、Tukey估计[12]、(L_p)估计[13]、Cauchy估计[14]、Welsch估计[15]等。本文借鉴M估计的思想，将(L_2)损失函数替换为稳健损失函数，提出适用于MAR模型的稳健估计方法，以提升模型在复杂数据场景下的拟合与预测性能。
本文的主要贡献如下：
1. 基于投影法和迭代最小二乘法，提出RE-PROJ和RE-ILS两种稳健估计方法，通过构造稳健损失函数，降低异常值对参数估计的影响;
2. 讨论了MAR(p)模型阶数选择的BIC准则，为模型阶数的确定提供依据;
3. 针对矩阵观测值存在相关性的情况，提出带有减秩假定的稳健估计方法(RR-RE-PROJ和RR-RE-ILS);
4. 通过模拟数据和实际数据分析，验证了所提方法的有效性和优越性。
2 MAR模型及传统估计方法
2.1 MAR模型定义
给定矩阵值时间序列({X_t}(t=1,cdots,T))，每个时点(t)的观测个体(X_t)为(m×n)维矩阵。MAR(1)模型的定义如下：
[X_t = A X_{t-1} B^T + E_t] (1)
其中：
- (X_t)经过中心化处理，模型不包含常数项;
- (A)和(B)分别为(m×m)维和(n×n)维的系数矩阵，为保证识别性，要求(|A|_F=1)且(A(1,1))的符号为正((|cdot|_F)表示矩阵的Frobenius范数);
- (E_t)是(m×n)维白噪声矩阵，时间维度上与(E_{t-1})不相关，截面维度上元素之间可能存在相关性。
2.2 传统估计方法
Chen等[7]针对MAR(1)模型提出了两种传统估计方法：
2.2.1 投影法(PROJ)
将(X_t)向量化为( ilde{X}_t = vec(X_t))((vec(cdot))表示按列拉直算子)，( ilde{E}_t = vec(E_t))，则MAR(1)模型转化为VAR(1)模型：
[ ilde{X}_t = Phi ilde{X}_{t-1} + ilde{E}_t] (2)
其中(Phi = B otimes A)((otimes)表示Kronecker积)。(Phi)的最小二乘估计为：
[hat{Phi} = left(sum_{t=2}^T ilde{X}_t ilde{X}_{t-1}^T ight)left(sum_{t=2}^T ilde{X}_{t-1} ilde{X}_{t-1}^T ight)^{-1}] (3)
通过Kronecker积分解(hat{Phi})，并调整(A)和(B)的范数与符号，得到系数矩阵估计((hat{A}, hat{B}))。
2.2.2 迭代最小二乘法(ILS)
构造目标函数：
[Q(A,B) = sum_{t=2}^T ext{trace}left((X_t - A X_{t-1} B^T)^T (X_t - A X_{t-1} B^T) ight)] (4)
令目标函数对(A)和(B)的偏导数为零，得到迭代更新公式：
[A^{(k+1)} = left(sum_{t=2}^T X_t (B^{(k)})^T X_{t-1}^T ight)left(sum_{t=2}^T A^{(k)} X_{t-1} (B^{(k)})^T (B^{(k)}) X_{t-1}^T ight)^{-1}] (5)
[B^{(k+1)} = left(sum_{t=2}^T X_{t-1}^T (A^{(k+1)})^T X_t ight)left(sum_{t=2}^T X_{t-1}^T (A^{(k+1)})^T A^{(k+1)} X_{t-1} ight)^{-1}] (6)
给定初始值((A^{(0)}, B^{(0)}))，迭代更新直至收敛，得到估计((hat{A}, hat{B}))。
3 MAR模型的稳健估计方法
3.1 稳健损失函数构造
传统估计方法基于(L_2)损失函数，易受异常值影响。本文借鉴M估计思想，将损失函数替换为稳健损失函数( ho(cdot))，满足：
- 非负且对称;
- 在零点取得唯一最小值;
- 增长速度慢于平方函数。
定义矩阵算子( ho(cdot))：对于任意矩阵(Y)，( ho(Y))与(Y)同规模，且( ho(Y)_{ij} = ho(Y_{ij}))。
3.2 稳健估计方法
3.2.1 基于投影法的稳健估计(RE-PROJ)
将MAR(1)模型向量化后，构造稳健目标函数：
[Q(Phi) = sum_{t=2}^T 1_m^T ho( ilde{X}_t - Phi ilde{X}_{t-1}) 1_n] (7)
其中(1_m)和(1_n)分别为(m)维和(n)维全1向量。设( ho)的导函数为(varphi(cdot))(矩阵算子定义同( ho(cdot)))，令目标函数对(Phi)的导数为零，得到：
[sum_{t=2}^T W_t ( ilde{X}_t - Phi ilde{X}_{t-1}) ilde{X}_{t-1}^T = 0] (8)
其中(W_t = varphi( ilde{X}_t - Phi^{(k)} ilde{X}_{t-1}) circ gamma( ilde{X}_t - Phi^{(k)} ilde{X}_{t-1}))((circ)表示Hadamard积，(gamma(Y)_{ij}=1/Y_{ij})为倒数算子)。
迭代求解(Phi)，收敛后通过Kronecker积分解得到((hat{A}, hat{B}))，算法步骤如下：
算法1：RE-PROJ估计步骤
1. 初始化(Phi^{(0)})(可采用PROJ方法的初始估计);
2. 计算残差( ilde{E}_t^{(k)} = ilde{X}_t - Phi^{(k)} ilde{X}_{t-1})，构造权重矩阵(W_t^{(k)});
3. 求解(Phi^{(k+1)} = left(sum_{t=2}^T W_t^{(k)} ilde{X}_t ilde{X}_{t-1}^T ight)left(sum_{t=2}^T W_t^{(k)} ilde{X}_{t-1} ilde{X}_{t-1}^T ight)^{-1});
4. 重复步骤2-3，直至(|Phi^{(k+1)} - Phi^{(k)}|_F < epsilon)((epsilon)为收敛阈值);
5. 对(hat{Phi})进行Kronecker积分解，调整(A)和(B)的范数与符号，得到((hat{A}, hat{B}))。
3.2.2 基于迭代最小二乘法的稳健估计(RE-ILS)
构造稳健目标函数：
[Q(A,B) = sum_{t=2}^T 1_m^T ho(X_t - A X_{t-1} B^T) 1_n] (9)
令目标函数对(A)和(B)的偏导数为零，得到迭代更新公式：
[vec(A^{(k+1)}) = left(sum_{t=2}^T G_t^{(k)} ilde{Z}_t^{(k)} G_t^{(k)T} ight)^{-1} sum_{t=2}^T G_t^{(k)} ilde{Z}_t^{(k)} vec(X_t)] (10)
[vec((B^{(k+1)})^T) = left(sum_{t=2}^T H_{t-1}^{(k)} check{Z}_t^{(k)} H_{t-1}^{(k)T} ight)^{-1} sum_{t=2}^T H_{t-1}^{(k)} check{Z}_t^{(k)} vec(X_t^T)] (11)
其中：
- (Z_t^{(k)} = varphi(X_t - A^{(k)} X_{t-1} (B^{(k)})^T) circ gamma(X_t - A^{(k)} X_{t-1} (B^{(k)})^T));
- ( ilde{Z}_t^{(k)} = diag[vec(Z_t^{(k)})])，(check{Z}_t^{(k)} = diag[vec((Z_t^{(k)})^T)]);
- (G_t^{(k)} = [(X_{t-1} (B^{(k)})^T) otimes I_m])，(H_{t-1}^{(k)} = [(X_{t-1}^T (A^{(k)})^T) otimes I_n])。
算法2：RE-ILS估计步骤
1. 初始化((A^{(0)}, B^{(0)}))(可采用ILS方法的初始估计);
2. 计算残差(E_t^{(k)} = X_t - A^{(k)} X_{t-1} (B^{(k)})^T)，构造权重矩阵(Z_t^{(k)});
3. 按公式(10)-(11)更新(A^{(k+1)})和(B^{(k+1)})，并调整范数与符号;
4. 重复步骤2-3，直至(|A^{(k+1)} - A^{(k)}|_F + |B^{(k+1)} - B^{(k)}|_F < epsilon);
5. 收敛后得到((hat{A}, hat{B}))。
3.3 常用稳健损失函数
1. (L_1)损失函数：( ho(x) = |x|)，对异常值敏感度低;
2. Huber损失函数[10]：
[
ho(x) =
egin{cases}
x^2/2, & |x| leq k
k(|x| - k/2), & |x| > k
end{cases}
]
其中(k)通常取1.345(hat{sigma})，(hat{sigma} = ext{MAR}/0.6745)(MAR为残差绝对值的中位数);
3. Tukey损失函数[12]：
[
ho(x) =
egin{cases}
k^2/6 [1 - (1 - (x/k)^2)^3], & |x| leq k
k^2/6, & |x| > k
end{cases}
]
其中(k)通常取4.685(hat{sigma})，对大残差赋予零权重。
4 模拟数据分析
4.1 模拟数据生成
模拟数据基于MAR(1)模型(X_t = A X_{t-1} B^T + E_t)生成，设置如下：
- 样本量(T)：100、200、400、5000;
- 矩阵规模((m,n))：(3,2)、(6,4)、(9,6);
- 平稳性条件：( ho(A) cdot ho(B) = 0.5)(( ho(cdot))表示谱半径)，(|A|_F=1)且(A(1,1)>0);
- 误差项(E_t)的8种分布设置：
1. 联合正态分布，协方差矩阵为单位阵;
2. 联合正态分布，协方差矩阵(sum_{ij}=2×0.5^{|i-j|});
3. 联合T(3)分布，独立;
4. 联合T(3)分布，相关系数矩阵(C_{ij}=0.5^{|i-j|});
5. 联合T(1)分布，独立;
6. 联合T(1)分布，相关系数矩阵(C_{ij}=0.5^{|i-j|});
7. 联合La(0,1)分布，独立;
8. 联合La(0,1)分布，协方差矩阵(sum_{ij}=0.5^{|i-j|})。
每种设置重复模拟100次，采用PROJ、ILS、RE-PROJ、RE-ILS四种方法估计参数，以(log(|hat{B} otimes hat{A} - B otimes A|_F^2))作为估计误差指标。
4.2 模拟结果分析
模拟结果如图1-8所示，主要结论如下：
1. 样本量(T)增大时，四种方法的估计误差均减小，且稳健估计的优势更明显;
2. 矩阵规模((m,n))较小时，正态分布下四种方法差异不大;规模较大时，RE-ILS和ILS表现更优;
3. 误差服从厚尾分布(T(3)、T(1)、La(0,1))时，RE-PROJ和RE-ILS显著优于PROJ和ILS，其中RE-ILS表现最佳;
4. 误差项是否相关对估计结果影响较小;
5. 样本量较小时，RE-PROJ表现可能不如ILS，但随着样本量增加，其优势逐渐显现。
(图1-8 不同误差分布下四种方法的估计误差对比图，原文中图表编号及对应分布如下：
图1：正态分布(独立);图2：正态分布(相关);图3：T(3)分布(独立);图4：T(3)分布(相关);
图5：T(1)分布(独立);图6：T(1)分布(相关);图7：La(0,1)分布(独立);图8：La(0,1)分布(相关))
5 扩展讨论
5.1 MAR(p)模型阶数选择的BIC准则
对于MAR(p)模型：
[X_t = sum_{i=1}^p A_i X_{t-i} B_i^T + E_t] (12)
定义BIC准则为：
[BIC(p) = log|hat{sum}| + frac{(m^2 + n^2)p log T}{T}] (13)
其中(hat{sum})为残差项(vec(E_t))的协方差矩阵估计，选择使BIC(p)最小的(p)作为最优阶数。
通过模拟验证(表1)，BIC准则能够有效选择MAR(p)模型的阶数，稳健高维协方差估计下的准确性略高于普通协方差估计。
表1 BIC准则对阶数p选择的概率
| 阶数 | 普通协方差矩阵估计 | | | | 稳健高维协方差估计 | | | |
||--|--|--|--|--|--|--|--|
| | 设置1(T(1)独立) | 设置2(T(1)相关) | 设置1(T(1)独立) | 设置2(T(1)相关) | 设置1(T(1)独立) | 设置2(T(1)相关) | 设置1(T(1)独立) | 设置2(T(1)相关) |
| | RE-PROJ | RE-ILS | RE-PROJ | RE-ILS | RE-PROJ | RE-ILS | RE-PROJ | RE-ILS |
| p=1 | 1% | 1% | 0% | 0% | 4% | 7% | 7% | 3% |
| p=2 | 10% | 8% | 9% | 1% | 16% | 18% | 14% | 12% |
| p=3 | 80% | 80% | 85% | 89% | 73% | 75% | 78% | 83% |
| p=4 | 4% | 4% | 3% | 6% | 5% | 0% | 1% | 1% |
| p=5 | 5% | 7% | 3% | 4% | 2% | 0% | 0% | 1% |
| 准确性 | 80% | 80% | 85% | 89% | 73% | 75% | 78% | 83% |
5.2 带有减秩假定的稳健估计方法
当系数矩阵(A)和(B)不满秩时，提出减秩稳健估计方法：
5.2.1 减秩投影法稳健估计(RR-RE-PROJ)
在RE-PROJ迭代过程中，对每次得到的(A^{(k)})和(B^{(k)})进行SVD分解，取前( ilde{m})(( ilde{m}
5.2.2 减秩迭代最小二乘法稳健估计(RR-RE-ILS)
类似RR-RE-PROJ，在RE-ILS迭代过程中，对(A^{(k+1)})和(B^{(k+1)})进行SVD分解并降秩，调整范数与符号后继续迭代。
算法3：RR-RE-PROJ估计步骤
1. 按RE-PROJ步骤得到(A^{(k)})和(B^{(k)});
2. 对(A^{(k)})进行SVD分解：(A^{(k)} = U sum V^T)，取前( ilde{m})列得到(A_R^{(k)} = U_1 sum_1 V_1^T)((rank(A_R^{(k)})= ilde{m}));
3. 同理得到(B_R^{(k)});
4. 重复迭代，直至收敛。
算法4：RR-RE-ILS估计步骤
1. 按RE-ILS步骤得到(A^{(k+1)})和(B^{(k+1)});
2. 对(A^{(k+1)})和(B^{(k+1)})进行SVD分解并降秩，得到(A_R^{(k+1)})和(B_R^{(k+1)});
3. 调整范数与符号，重复迭代直至收敛。
模拟结果(图9-10)表明：
1. 样本量增大时，估计误差减小;
2. 数据量充足((m^2 + n^2 < T))时，满秩稳健估计表现更优;
3. 低维数据((m^2 + n^2 > T))时，减秩稳健估计(RR-RE-PROJ、RR-RE-ILS)显著优于满秩估计。
(图9：系数矩阵满秩时四种方法的估计误差;图10：系数矩阵不满秩时四种方法的估计误差)
6 实际数据分析
6.1 数据描述
采用两组实际数据验证所提方法：
- 数据1：宏观经济数据(2000年1月-2021年5月，(T=257))，包含4个指标(CPI、RPI、Export、进口)和3个地区(北京、天津、河北)，为(4×3)维矩阵;
- 数据2：行业股票数据(2000年1月-2021年5月，(T=257))，包含5个行业(轻工制造、医药生物、建筑装饰、电子、食品饮料)和4个指标(收盘价、成交量、成交额、收益率)，为(5×4)维矩阵。
两组数据均存在异常值和尖峰厚尾特征(表2)，峰度均超过6，异常值占比大于0.27%(标准正态分布(3sigma)以外的概率)。
表2 实际数据的异常值特征与拟合误差
| 数据 | 峰度 | 异常值比例 | 原始数据离差平方和 | 拟合误差(SSE) | | | |
||||--|||||
| | | | | PROJ | ILS | RE-PROJ | RE-ILS |
| 数据1 | 9.11 | 1.79% | 3081 | 399.05 | 367.70 | 398.20 | 368.82 |
| 数据2 | 12.00 | 2.10% | 5156 | 5125.94 | 2461.88 | 5048.34 | 1751.27 |
6.2 实证结果分析
计算四种方法的拟合误差(SSE)和预测误差((FE_1)为平均预测误差，(FE_2)为中位数预测误差，分别基于(L_2)和(L_infty)范数)，结果如表3所示：
表3 实际数据的预测误差
| 数据 | 误差类型 | PROJ | ILS | RE-PROJ | RE-ILS |
||-||--||--|
| 数据1 | (FE_1(L_2)) | 0.33 | 0.21 | 0.26 | 0.17 |
| | (FE_2(L_2)) | 0.28 | 0.19 | 0.22 | 0.15 |
| | (FE_1(L_infty)) | 0.34 | 0.23 | 0.27 | 0.18 |
| | (FE_2(L_infty)) | 0.28 | 0.20 | 0.22 | 0.16 |
| 数据2 | (FE_1(L_2)) | 6.25 | 1.54 | 5.89 | 1.46 |
| | (FE_2(L_2)) | 5.76 | 1.39 | 4.94 | 1.32 |
| | (FE_1(L_infty)) | 8.83 | 1.57 | 8.40 | 1.52 |
| | (FE_2(L_infty)) | 6.81 | 1.28 | 6.35 | 1.25 |
实证结果表明：
1. 拟合误差方面，ILS和RE-ILS显著优于PROJ和RE-PROJ，其中RE-ILS表现最佳;
2. 预测误差方面，无论是(L_2)还是(L_infty)范数，RE-ILS的预测误差均最小，ILS次之，PROJ和RE-PROJ表现较差;
3. 数据2的异常值和峰度更显著，RE-ILS的优势更突出，验证了稳健估计在复杂数据场景下的优越性。
7 结论与展望
本文针对矩阵自回归(MAR)模型的最小二乘估计易受异常值和厚尾误差影响的问题，提出了RE-PROJ和RE-ILS两种稳健估计方法，通过构造稳健损失函数，有效降低了异常值对参数估计的干扰。进一步讨论了MAR(p)模型阶数选择的BIC准则和带有减秩假定的稳健估计方法，丰富了MAR模型的理论体系。
模拟数据和实际数据分析均表明，所提稳健估计方法在异常值较多或尖峰厚尾的数据中，拟合与预测性能显著优于传统最小二乘估计，尤其是RE-ILS方法表现最佳。
未来研究方向：
1. 探索其他稳健损失函数(如分位数损失函数[23])在MAR模型中的应用;
2. 结合正则化方法[24]，提出带有惩罚项的稳健估计，进一步提升模型的解释性;
3. 将MAR模型的稳健估计拓展到高维数据场景，解决更大规模矩阵值时间序列的建模问题。
8 定理证明
定理1(RE-PROJ估计的渐近正态性)
假设：
1. (E_1,cdots,E_T)独立同分布，均值为0，协方差矩阵为(sum)，且满足Lebesgue密度函数连续;
2. 满足平稳性条件( ho(A) cdot ho(B) < 1)，且矩阵(A)、(B)和(sum)非奇异;
3. 稳健损失函数的导函数(varphi(cdot))满足正则条件。
则当(T o infty)时，有：
[
sqrt{T} egin{pmatrix} vec(hat{A}_1 - A) vec(hat{B}_1^T - B^T) end{pmatrix} Rightarrow N(0, V_0 Xi_1 V_0^T)
]
其中：
- (hat{A}_1)、(hat{B}_1)为RE-PROJ估计;
- (V_0 = [eta otimes I_m, I_n otimes alpha])，(alpha = vec(A))，(eta = vec(B));
- (Xi_1 = G(Lambda))，(Lambda = (M Gamma_0 M^T)^{-1} otimes (M sum M^T))，(M = diag(vec(W_t^s circ W_t^s)))，(Gamma_0 = Cov(vec(X_t)));
- (W_t^s = varphi(E_t))，(G(cdot))为矩阵重新排布算子。
定理2(RE-ILS估计的渐近正态性)
在定理1的假设条件下，当(T o infty)时，有：
[
sqrt{T} egin{pmatrix} vec(hat{A}_2 - A) vec(hat{B}_2^T - B^T) end{pmatrix} Rightarrow N(0, Xi_2)
]
其中：
- (hat{A}_2)、(hat{B}_2)为RE-ILS估计;
- (Xi_2 = H^{-1} mathbb{E}(G_t M sum M G_t^T) H^{-1});
- (H = mathbb{E}(G_t M G_t^T) + gamma gamma^T)，(gamma = (alpha^T, 0^T)^T in mathbb{R}^{m^2 + n^2});
- (G_t^T = [(B X_t^T) otimes I_m, I_n otimes (A X_t)])。
(详细证明过程见原文附录)
参考文献
[1] Hannan E J. Multiple time series[M]. John Wiley & Sons, 1970.
[2] Engle R F, Granger C W J. Co-integration and error correction: representation, estimation, and testing[J]. Econometrica: Journal of the Econometric Society, 1987: 251–276.
[3] Tsay R S. Testing and modeling multivariate threshold models[J]. Journal of the American Statistical Association, 1998, 93(443): 1188–1202.
[4] Zhou H, Li L. Regularized matrix regression[J]. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2014, 76(2): 463–483.
[5] Zhao J, Leng C. Structured lasso for regression with matrix covariates[J]. Statistica Sinica, 2014: 799–814.
[6] Wang D, Liu X, Chen R. Factor models for matrix-valued high-dimensional time series[J]. Journal of Econometrics, 2019, 208(1): 231–248.
[7] Chen R, Xiao H, Yang D. Autoregressive models for matrix-valued time series[J]. Journal of Econometrics, 2021, 222(1): 539–560.
[8] Kamm J, Nagy J G. Optimal Kronecker product approximation of block Toeplitz matrices[J]. SIAM Journal on Matrix Analysis and Applications, 2000, 22(1): 155–172.
[9] Machado J A F. Robust model selection and M-estimation[J]. Econometric Theory, 1993: 478–493.
[10] Huber P J. Robust Estimation of a Location Parameter[J]. The Annals of Mathematical Statistics, 1964: 73–101.
[11] Sheskin D J. Handbook of parametric and nonparametric statistical procedures[M]. CRC Press, 2020.
[12] Sinova B, Van Aelst S. Advantages of M-estimators of location for fuzzy numbers based on Tukey’s biweight loss function[J]. International Journal of Approximate Reasoning, 2018, 93: 219–237.
[13] Hofmann M. (L_p) estimation of the diffusion coefficient[J]. Bernoulli, 1999, 5(3): 447–481.
[14] McCullagh P. M¨obius transformation and Cauchy parameter estimation[J]. Annals of Statistics, 1996, 24(2): 787–808.
[15] Krasker W S, Welsch R E. Efficient bounded-influence regression estimation[J]. Journal of the American Statistical Association, 1982, 77(379): 595–604.
[16] 陈丙振, 孔令臣, 尚蓬勃. 稳健矩阵回归模型及算法研究[J]. 数值计算与计算机应用, 2018, 40(4): 402–417.
[17] Kamm J, Nagy J G. Optimal Kronecker product approximation of block Toeplitz matrices[J]. SIAM Journal on Matrix Analysis and Applications, 2000, 22(1): 155–172.
[18] Goes J, Lerman G, Nadler B. Robust sparse covariance estimation by thresholding Tyler’s M-estimator[J]. The Annals of Statistics, 2020, 48(1): 86–110.
[19] L¨utkepohl H. Forecasting contemporaneously aggregated vector ARMA processes[J]. Journal of Business & Economic Statistics, 1984, 2(3): 201–214.
[20] Stock J H, Watson M W. Forecasting inflation[J]. Journal of Monetary Economics, 1999, 44(2): 293–335.
[21] Stock J H, Watson M W. Macroeconomic forecasting using diffusion indexes[J]. Journal of Business & Economic Statistics, 2002, 20(2): 147–162.
[22] Lam C, Yao Q. Factor modeling for high-dimensional time series: inference for the number of factors[J]. The Annals of Statistics, 2012, 40(2): 694–726.
[23] Koenker R, Hallock K F. Quantile regression[J]. Journal of Economic Perspectives, 2001, 15(4): 143–156.
[24] Hoerl A E, Kennard R W. Ridge regression: applications to nonorthogonal problems[J]. Technometrics, 1970, 12(1): 69–82.
[25] Dunsmuir W, Hannan E J. Vector linear time series models[J]. Advances in Applied Probability, 1976, 8(2): 339–364.
[26] Muler N, Yohai V J. Robust estimation for vector autoregressive models[J]. Computational Statistics & Data Analysis, 2013, 65: 68–79.
[27] Hall P, Heyde C C. Martingale Limit Theory and Its Applications[M]. Academic Press, New York, 1980.
