企业官网建设流程全解析

做时间序列分析的人，几乎都会遇到同一个问题：手里一串数据，今天和昨天相关，昨天和前天相关，但这种相关性到底是怎么衰减的？是一天比一天弱，还是每隔几期又反弹回来？要回答这个问题，最直观的方式就是看它的自相关矩阵。而在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 ) 再取期望。以 ( k=0 )、( k=1 )、( k=2 ) 代进去，会得到三组方程，最后可以解出：

[ \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, monthly

tsset后面跟时间变量和频率。如果不确定频率，用describe date先看看变量格式。声明完时间变量之后，先画个时序图：

tsline y

这一步别跳过。我习惯先肉眼看一遍：有没有明显趋势？有没有周期性？大概在哪个水平附近波动？ARMA模型要求序列是平稳的，如果你看到明显的上升趋势，后面就要考虑差分，或者退一步说，至少要意识到均值的估计会受影响。

接下来用相关图和单位根检验做个定量判断：

corrgram y, lags(20) dfuller y

corrgram会同时给出ACF、PACF和Q统计量；dfuller是ADF检验。这两条命令配合使用，基本能判断序列是否平稳。

如果序列不平稳，常见处理是先取对数消除方差扩张，再一阶差分消除趋势；如果季节性明显，可能还要做季节差分。检验到平稳之后再进入识别环节。

3.2 用ACF和PACF判断AR还是MA

识别ARMA阶数最常用的工具就是corrgram输出的ACF和PACF。原则并不复杂，我整理成表格：

举个例子。假设数据由一个AR(2)过程生成，样本量300，corrgram输出大致会是这样：

LAG AC PAC Q Prob>Q 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 aroots：

estat aroots

它会画出AR和MA多项式的特征根倒数，并且用图形判断是否都落在单位圆内。如果所有点都在单位圆内，说明模型平稳可逆。文本输出里也会给出特征多项式的根的模长，我的判断标准是：所有根的模长严格小于1，如果有任何一个落在0.98以上，就要怀疑样本量不够或者模型设定有问题。

接下来是残差检验。模型估计完，把残差存下来：

predict res, resid wntestq res, lags(10)

wntestq是Ljung-Box Q检验，原假设是残差没有自相关。p值大于0.05说明残差基本是白噪声，模型拟合到位。我通常会再画一眼残差的ACF：

corrgram 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 L`i' = L`i'.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.15；ACF该拖尾，结果前几阶就正负交替。遇到这种情况，我的建议是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的输出找线索，这两样东西能回答大部分疑惑。