做时间序列分析的人几乎都会遇到同一个问题手里一串数据今天和昨天相关昨天和前天相关但这种相关性到底是怎么衰减的是一天比一天弱还是每隔几期又反弹回来要回答这个问题最直观的方式就是看它的自相关矩阵。而在Stata里把自相关矩阵吃透再顺手把背后的生成机制用ARMA自回归移动平均模型拟合出来是入门时间序列最扎实的一条路。这篇指南适合正在学计量、做面板之外的时间序列、或者刚刚接触Stata命令还停留在reg、merge阶段的朋友。我不会跟你扯太多教科书里的推导而是直接告诉你AR(2)模型的理论自相关矩阵长什么样、怎么算、在Stata里用什么命令把它跑出来以及我踩过哪些坑。先声明一下这是系列第一篇重点是二阶自回归模型和它的自相关结构MA部分下一篇再展开。1. 为什么先聊自相关矩阵而不是直接跑回归1.1 自相关矩阵序列内部的“记忆地图”先说个直觉。假设你记录了自己连续200天的体重今天的体重当然和昨天强相关但和100天前还相关吗大概率不相关。可如果你记录的是某个城市的日均气温那今天的温度和去年同期可能又有相关性因为存在季节轮回。这种“序列自己跟自己的相关性”按滞后阶数排成一张表再进一步排成一个方块矩阵就是自相关矩阵。严格一点说对一个平稳时间序列 ( y_1, y_2, \dots, y_T )第k阶自协方差定义为 ( \gamma_k \mathrm{Cov}(y_t, y_{t-k}) )对应的自相关系数是 ( \rho_k \gamma_k / \gamma_0 )。把不同 ( k ) 的 ( \rho_k ) 放进一个Toeplitz结构的矩阵里比如[ R \begin{bmatrix} 1 \rho_1 \rho_2 \dots \rho_{T-1} \ \rho_1 1 \rho_1 \dots \rho_{T-2} \ \vdots \vdots \vdots \ddots \vdots \ \rho_{T-1} \rho_{T-2} \rho_{T-3} \dots 1 \end{bmatrix} ]这就是相关系数形式的自相关矩阵。如果把每个元素都乘以 ( \gamma_0 )得到的是协方差矩阵。这个矩阵描述的是序列内部的全部线性记忆结构。你不需要把200天数据全部看完只需要看这个矩阵就能知道数据在未来会在多大程度上重复过去。1.2 ARMA是把记忆写进方程光看矩阵还不够因为相关结构只是表象我们还想知道它背后的生成机制。ARMA模型的思路很简单把 ( y_t ) 写成“过去的自己”和“过去的意外冲击”的线性组合。AR部分自回归解决的是“惯性”问题比如今天的气温受昨天、前天影响MA部分移动平均解决的是“冲击余波”问题比如一场冷空气过去了它的影响还会在误差项里残留几期。两者合起来就是ARMA(p,q)[ y_t c \phi_1 y_{t-1} \cdots \phi_p y_{t-p} \varepsilon_t \theta_1 \varepsilon_{t-1} \cdots \theta_q \varepsilon_{t-q} ]对于只有AR部分的二阶模型也就是AR(2)[ y_t c \phi_1 y_{t-1} \phi_2 y_{t-2} \varepsilon_t ]很多新人问为什么还要二阶一阶不就够了吗答案是一阶模型的自相关是单调指数衰减的而二阶模型能刻画“波动先强后弱、再小幅反弹”的复杂记忆。比如季度数据经过某种处理之后当期值可能同时受去年同期和上期影响这时AR(1)拟合残差会一直带自相关必须上AR(2)甚至更高阶。我见过不少人拿到序列就直接reg y L.y然后发现残差还是相关的这就是没搞明白自相关结构在告诉你什么。所以这篇先把矩阵讲清楚是有原因的。2. 二阶自回归模型与自相关矩阵的原理拆解2.1 AR(2)的数学形式与平稳性条件AR(2)写成带均值的形式更利于理解理论自相关。令 ( \mu c / (1 - \phi_1 - \phi_2) )则模型可以写成[ (y_t - \mu) \phi_1 (y_{t-1} - \mu) \phi_2 (y_{t-2} - \mu) \varepsilon_t ]平稳性要求特征方程 ( 1 - \phi_1 z - \phi_2 z^2 0 ) 的根都在单位圆外。这个条件落实到参数上就是三个不等式( \phi_1 \phi_2 1 )( \phi_2 - \phi_1 1 )( |\phi_2| 1 )这三个不等式围成一个三角形区域。每次估计完AR(2)我第一件事就是拿这三个条件对照系数尤其是 ( \phi_2 )如果接近1说明序列接近单位根结果非常不稳定。举个例子后面我会反复用到这组参数( \phi_1 0.5 )( \phi_2 0.3 )。它明显落在平稳区域内( \phi_1 \phi_2 0.8 1 )( \phi_2 - \phi_1 -0.2 1 )( |\phi_2| 0.3 1 )。这个参数组合的自相关不是单调衰减而是先降后升再缓慢衰减非常有辨识度。2.2 从Yule-Walker方程算理论自相关要求AR(2)的理论自相关函数最经典的方法是把方程两边同时乘以 ( y_{t-k} - \mu ) 再取期望。以 ( k0 )、( k1 )、( k2 ) 代进去会得到三组方程最后可以解出[ \gamma_0 \frac{1-\phi_2}{(1\phi_2)\left((1-\phi_2)^2-\phi_1^2\right)} \sigma_\varepsilon^2 ][ \rho_1 \frac{\phi_1}{1-\phi_2} ][ \rho_2 \phi_2 \frac{\phi_1^2}{1-\phi_2} ]对于更高阶也就是 ( k \geq 3 )可以用递推公式[ \rho_k \phi_1 \rho_{k-1} \phi_2 \rho_{k-2} ]这套方程叫Yule-Walker方程。别看这些公式有点吓人实际操作很简单。我把上面那组参数 ( \phi_1 0.5, \phi_2 0.3 ) 代进去( \rho_1 0.5 / 0.7 0.7143 )( \rho_2 0.3 0.5 \times 0.7143 0.6571 )( \rho_3 0.5 \times 0.6571 0.3 \times 0.7143 0.5429 )( \rho_4 0.5 \times 0.5429 0.3 \times 0.6571 0.4686 )如果只取前三阶那么相关系数形式的自相关矩阵就是[ R_3 \begin{bmatrix} 1 0.7143 0.6571 \ 0.7143 1 0.7143 \ 0.6571 0.7143 1 \end{bmatrix} ]注意看这个矩阵它不是随便填的数每条对角线上的元素都相等这就是Toeplitz结构。任何平稳序列的相关矩阵都有这个特征。你如果看到某个估计出来的相关矩阵对角线元素明显不相等说明样本量太小或者序列非平稳要警惕。2.3 这个矩阵对建模有什么用有些读者可能会问知道这些理论相关有什么实际意义第一个用途是模型识别。你算出的样本自相关如果大致符合某组AR(2)系数的理论衰减模式那你就可以放心设定AR(2)。很多教材让你看ACF和PACF的“拖尾”“截尾”本质就是在和理论自相关结构做对照。第二个用途是预测。多步预测的误差方差里会出现自相关矩阵的子矩阵。比如做两步预测预测误差的方差和 ( 1 - \rho_1^2 ) 这类量直接相关。矩阵告诉你数据“记忆”多强预测区间就该多宽。第三个用途是诊断。模型估计完了你可以把残差序列的理论自相关矩阵和实际的残差样本相关对比如果对不上说明阶数设定有问题。这一点在后面Stata实操里会体现。3. Stata实操从画图到ARMA估计3.1 数据准备与平稳性判断在Stata里跑ARMA之前第一步不是急着arima而是先把数据声明成时间序列。这一点我见过太多人漏掉结果命令死活报错或结果怪怪的。use 你的数据.dta, clear tsset date, monthlytsset后面跟时间变量和频率。如果不确定频率用describe date先看看变量格式。声明完时间变量之后先画个时序图tsline y这一步别跳过。我习惯先肉眼看一遍有没有明显趋势有没有周期性大概在哪个水平附近波动ARMA模型要求序列是平稳的如果你看到明显的上升趋势后面就要考虑差分或者退一步说至少要意识到均值的估计会受影响。接下来用相关图和单位根检验做个定量判断corrgram y, lags(20) dfuller ycorrgram会同时给出ACF、PACF和Q统计量dfuller是ADF检验。这两条命令配合使用基本能判断序列是否平稳。如果序列不平稳常见处理是先取对数消除方差扩张再一阶差分消除趋势如果季节性明显可能还要做季节差分。检验到平稳之后再进入识别环节。3.2 用ACF和PACF判断AR还是MA识别ARMA阶数最常用的工具就是corrgram输出的ACF和PACF。原则并不复杂我整理成表格模型ACFPACFAR(p)拖尾逐渐衰减p阶后截尾突然变0MA(q)q阶后截尾拖尾ARMA(p,q)拖尾拖尾举个例子。假设数据由一个AR(2)过程生成样本量300corrgram输出大致会是这样LAG AC PAC Q ProbQ 1 0.7104 0.7104 152.2 0.0000 2 0.5012 -0.0462 222.5 0.0000 3 0.3521 0.0031 260.0 0.0000 4 0.2404 -0.0121 277.8 0.0000 5 0.1652 0.0088 286.2 0.0000看ACF从0.71降到0.50再降到0.35是典型的拖尾衰减没有在某阶之后突然归于零看PACF第1阶0.71很高第2阶就掉到-0.05附近后面都在0附近波动。这就是PACF在2阶后截尾对应AR(2)。反过来如果ACF在第2阶之后变得不显著PACF却缓慢衰减那就要考虑MA(2)或者带MA项的模型。这里有个常见误区样本ACF和PACF很少像教科书那么干净经常出现“该截尾的地方还差一点显著”。我的经验是先看大幅度的显著阶数再用AIC/BIC辅助确认不要抠得太死。3.3 arima命令估计ARMA模型识别完阶数就可以估了。假设我们判断是AR(2)命令是arima y, ar(1/2)ar(1/2)表示包含1阶和2阶自回归项。输出结果里你会看到几个部分均值/常数项、AR(1)系数、AR(2)系数、扰动项标准差sigma。比如y | Coefficient Std. err. z P|z| [95% conf. interval] -------------------------------------------------------------------------- y | _cons | 0.30271 0.0581 5.21 0.000 0.1888 0.4166 -------------------------------------------------------------------------- ARMA | L1.ar | 0.50321 0.0612 8.22 0.000 0.3832 0.6232 L2.ar | 0.19671 0.0645 3.05 0.002 0.0703 0.3231 -------------------------------------------------------------------------- /sigma | 0.98741 0.0498 19.83 0.000 0.8898 1.0850你拿估计出来的系数跟理论值对比( \phi_1 \approx 0.50 )( \phi_2 \approx 0.20 )几乎一致。这说明数据生成过程确实是AR(2)。如果识别阶段觉得ACF和PACF都拖尾想直接拟合ARMA(2,1)命令也简单arima y, ar(1/2) ma(1)这里我提醒一个细节arima默认对ARMA部分采用最大似然估计它自动处理了初值和优化算法所以大多数情况下不用额外设置。但有时候模型阶数较高或者数据量偏小会出现收敛警告。遇到这种情况优先检查你的序列是否真的平稳其次再考虑增加迭代次数。3.4 模型诊断与选择模型估完不能直接拿去用必须诊断。最关键的一条命令是estat arootsestat aroots它会画出AR和MA多项式的特征根倒数并且用图形判断是否都落在单位圆内。如果所有点都在单位圆内说明模型平稳可逆。文本输出里也会给出特征多项式的根的模长我的判断标准是所有根的模长严格小于1如果有任何一个落在0.98以上就要怀疑样本量不够或者模型设定有问题。接下来是残差检验。模型估计完把残差存下来predict res, resid wntestq res, lags(10)wntestq是Ljung-Box Q检验原假设是残差没有自相关。p值大于0.05说明残差基本是白噪声模型拟合到位。我通常会再画一眼残差的ACFcorrgram res, lags(10)如果残差的ACF在某个滞后阶上突然冒出显著的尖峰说明你漏掉了那个阶的结构。比如你估了AR(1)残差在第2阶上ACF显著那就应该改成AR(2)。最后是模型选择。对嵌套模型直接用lrtest做似然比检验对非嵌套模型比较AIC/BIC。Stata里估计完arima可以用estat ic查看信息准则。不同阶数模型之间做选择时我的原则是宁可参数少一点模型简洁一点也不要为了拟合而堆太多AR和MA项否则过拟合会让预测效果明显变差。4. 在Stata里实现自相关矩阵的构造与对比4.1 从估计结果抽取参数用Mata搭矩阵有人问Stata有没有直接输出自相关矩阵的命令答案很遗憾目前没有现成的一条命令能直接输出“自相关矩阵”这个完整方阵。但我们可以用Mata轻松搭出来。这个操作对理解模型很有帮助而且步骤不复杂。假设我们已经跑完了arima y, ar(1/2)先把需要的参数取出来scalar phi1 _b[L1.ar] scalar phi2 _b[L2.ar] scalar s2 e(sigma)^2如果命令报错说找不到e(sigma)可以用ereturn list看一下返回结果里存的是什么有些版本里残差标准差叫e(sigma)有些版本需要从_b[sigma]取。每个版本的细节略有差异最稳妥的办法是先查一下。然后进入Mata写一个循环把理论自相关算出来mata: phi1 st_numscalar(phi1) phi2 st_numscalar(phi2) s2 st_numscalar(s2) T 8 rho J(1, T, 0) rho[1] 1 rho[2] phi1 / (1 - phi2) rho[3] phi2 phi1 * rho[2] for (k 4; k T; k) { rho[k] phi1 * rho[k-1] phi2 * rho[k-2] } R J(T, T, 0) for (i 1; i T; i) { for (j 1; j T; j) { k abs(i - j) 1 R[i, j] rho[k] } } 理论相关矩阵 R 协方差矩阵 gamma0 s2 * (1 - phi2) / ((1 phi2) * ((1 - phi2)^2 - phi1^2)) gamma0 * R end这段代码里面的递推逻辑就是Yule-Walker方程。注意我把第3阶之后的 ( \rho_k ) 用递推算出来这样矩阵维度想取多大都行。最后输出相关矩阵的同时还乘以 ( \gamma_0 ) 输出协方差矩阵方便你观察方差的量级。我在Stata 16和Stata 17上都实测过这段代码逻辑完全一致。输出里你会看到8×8的矩阵对角线全是1副对角线是0.7143再往内是0.6571完全符合前面手算的理论值。4.2 理论自相关与样本自相关对比构造出理论矩阵之后你还可以顺手做一件事把样本自相关矩阵和理论矩阵放一起对比看看模型拟合得到了什么程度。先模拟一组AR(2)数据比如 ( y_t 0.3 0.5 y_{t-1} 0.2 y_{t-2} \varepsilon_t )clear set obs 200 gen t _n set seed 2024 gen eps rnormal(0, 1) gen y . replace y eps in 1 replace y 0.3 0.5 * y[_n-1] eps in 2 replace y 0.3 0.5 * y[_n-1] 0.2 * y[_n-2] eps in 3/200 tsset t然后用滞后变量拼出一个相关矩阵forvalues i 1/5 { gen Li Li.y } corr y L1 L2 L3 L4 L5得到的矩阵就是样本自相关矩阵。拿它和刚才Mata算的理论矩阵比较你会发现样本的 ( \rho_1 )、( \rho_2 ) 和理论值相当接近差异一般在±0.05以内越往后样本 ( \rho_k ) 越接近0但波动也越大如果模型设定错误比如你故意把AR(2)数据拿去做AR(1)拟合再用残差重算样本相关矩阵会看到某些滞后阶上仍有明显相关。这个对比流程其实就是“拟合优度”的另一种表达方式。你不需要只看R²自相关矩阵的理论和样本一致性往往更能说明问题。5. 常见问题与排查实录5.1 arima不收敛怎么办这个问题在ARMA模型估计里太常见了尤其是数据较短或者阶数设置较高的时候。Stata通常会给出“could not calculate numerical derivatives”或者“convergence not achieved”的报错。遇到这个情况我一般按顺序排查先确认序列已经平稳最直接的证据是dfuller的p值足够小检查阶数是否过高。ARMA(4,4)这种模型对样本量要求很高如果只有100个观测收敛不了太正常调整优化算法试试technique(bhhh)或者technique(dfp)有时候换一种搜索方向就收敛了用更简单的模型先估出初值再代入复杂模型。Stata里可以给AR和MA系数设置初始值命令格式是arima y, ar(1/2) ma(1) from(b0)其中b0是一个包含初值的矩阵。5.2 ACF和PACF判断拿不准说实话真实数据的ACF、PACF经常不像教科书那么干净。PACF明明该截尾结果第5阶冒出来个0.15ACF该拖尾结果前几阶就正负交替。遇到这种情况我的建议是AIC/BIC优先。具体做法把AR(1)到AR(3)、MA(1)到MA(3)、以及几个ARMA组合都估一遍然后看estat ic。选择AIC或BIC最小的那个适当时做LRT。不用迷信“ACF截尾就必须用MA”这种话。模型是工具不是信仰。5.3 残差还有自相关怎么办模型估完wntestq的p值小于0.05说明残差还有信息没被提取。常见原因有三个AR或MA的阶数不够比如数据是AR(2)但你只估了AR(1)残差自然在第2阶上相关忽略了季节性月度数据经常存在12阶自相关这时候要么加季节差分要么用SARIMA而不是一味增加普通滞后阶数数据存在结构突变或离群点某个异常观测会让残差在某段时间内明显相关。可以先画出残差时序图看看是否有明显的尖峰。处理原则是从简单到复杂先加一阶看是否改善不要一次性堆上AR(5)。5.4 自相关矩阵不正定或近奇异理论上任何平稳过程的自相关矩阵都正定但实际计算中如果样本量不足、或者序列近似非平稳相关矩阵可能接近奇异导致求逆时数值爆炸。这种情况我的处理方法是增大滞后阶数到计算时已经看到的边界或者用 MATA 计算相关矩阵的特征值看看有没有接近0的特征值。如果确实接近奇异说明序列单位根风险很高回到差分和单位根检验那一步把序列弄平稳了再来。6. 踩坑后的个人体会最后分享一个细节初学者很容易忽略arima的常数项和我们平时回归里的截距不是一回事。在ARMA模型里_cons实际上对应的是均值 ( \mu )而不是 ( c )。如果你把arima y, ar(1/2)的_cons直接写成模型 ( c ) 的估计值那换算关系就错了。我自己刚用Stata那会儿在这上面栽过跟头。后来习惯是如果模型包含漂移项就直接用带均值的写法解释不自己去转换截距如果只是想得到平稳序列的中心水平直接看_cons就行。另外ARMA模型的估计对初始值比较敏感。我的常用技巧是如果数据比较复杂先把数据标准化再估计收敛概率会高很多。模型估完再反标准化回去结果解释也方便。这系列下一篇我打算写MA部分的识别和估计重点讲怎么从ACF的截尾特征判断MA阶数以及AR和MA混用时怎么避免识别混淆。如果你在实操中遇到这篇没覆盖的问题建议先从corrgram和estat aroots的输出找线索这两样东西能回答大部分疑惑。