时间序列分析：AR(1)与AR(2)模型从原理到Python实战

做时间序列分析的人，迟早会在某个深夜盯着一条缓慢衰减的自相关图发呆：这条线到底该用一个滞后项还是两个？我当年学 AR 模型的时候，教材把 AR(1)、AR(2) 单独拎出来讲了一整节，心里还嫌啰嗦——都是一样的递推结构，直接上 AR(p) 不就完了。后来越做越多才明白，这两阶模型几乎把时间序列里所有容易踩的坑都演了一遍：平稳性边界怎么划、参数为什么要加约束、什么时候会出现那种看着有周期性其实是虚假周期的波动、定阶错了会有什么后果，全在里面。这篇就把 AR(1) 和 AR(2) 掰开揉碎讲一遍，从 Yule-Walker 方程的推导，到用 Python 完整的模拟、定阶、估计、诊断、预测流程，再到几个只有真正跑过数据才会遇到的问题。前一百字里我把关键词都说清楚了：时间序列分析、AR 模型、AR(1)、AR(2)。

内容偏向实操，但推导部分我尽量写得能让人看懂，不跳步。适合刚学完平稳性、自相关函数、白噪声这些概念，准备动手建第一个模型的人；也适合已经在用statsmodels跑AutoReg、但说不清里面参数含义的人。

1. 为什么 AR(1) 和 AR(2) 值得单独拎出来讲

1.1 从"今天的值跟昨天的值有关"说起

AR 模型的核心思想朴素得不像个统计模型：今天的数据，是过去若干个数据的加权和，再加上一点当下无法解释的扰动。写成公式就是 X_t = φ₁X_{t-1} + φ₂X_{t-2} + … + φ_pX_{t-p} + ε_t，其中的 p 就是阶数。AR(1) 是说今天只跟昨天有关，AR(2) 是说今天跟昨天和前天都有关。

这个假设放在很多真实场景里是成立的。比如某只股票每天的波动率、某个城市一天的用电负荷、某种传染病每日新增、某个电商平台日活跃用户数——这些量都有一个共同特点：相邻时点的值高度相关，而且这种相关性会随着时间间隔拉长而衰减。AR 模型要做的就是把这层"记忆"用几个参数刻画出来。

那为什么偏偏是 1 阶和 2 阶被反复拿出来讲？因为它们是唯一两个能用纸笔算出全部解析解的阶数。p≥3 之后，Yule-Walker 方程变成一个需要数值求解的线性方程组，平稳域也变成高维空间里一个难以画出来的区域，你很难再对参数形成直观感觉。1 阶和 2 阶则是可以画在二维平面上、可以手推自相关函数、可以肉眼判断平稳性的。把这两阶吃透，后面遇到高阶模型，心里是有底的。

1.2 两阶模型是理解高阶模型的跳板

我在实际项目里几乎没遇到过纯 AR(1) 就够了的数据，但每一次建模的前几步，都会先看一眼一阶滞后系数。原因很简单：φ₁ 的大小直接决定这条序列的"记忆长度"，也就是自相关衰减的快慢。一个 φ₁=0.95 的序列，衰减极慢，看起来像随机游走，差分之后可能就平稳了；一个 φ₁=0.2 的序列，今天的信息两三天就忘干净了，用 AR(1) 描述基本够。

AR(2) 则多了另一层信息：二阶项 φ₂ 的符号。φ₂ 为正，说明序列存在同向的惯性叠加，自相关曲线是单调衰减或者缓慢波动；φ₂ 为负，往往对应那种"涨一天跌一天"式的交替振荡，这在日频数据里非常常见。很多人看到自相关图一波一波地荡，第一反应是"这数据有周期性，得做季节性差分"，其实只要加一个负的二阶项就能解释掉，贸然差分反而把信号毁了。

所以这两阶模型不是教学上的简化，而是实际建模中两把最常用的尺子。

1.3 需要的预备知识清单

在往下走之前，这几件事得先有概念，不然推导部分会卡住：

• 弱平稳：均值恒为常数，方差恒为常数，协方差只跟时间间隔有关，与绝对时点无关。
• 白噪声：均值为零、方差恒定、彼此不相关的随机序列，记作 ε_t ~ WN(0, σ²)。它是 AR 模型的"原料"。
• 自协方差 γ_k 与自相关 ρ_k：γ_k = Cov(X_t, X_{t-k})，ρ_k = γ_k / γ_0。ρ_0 恒等于 1。
• 滞后算子 L：LX_t = X_{t-1}，它的作用就是"往回拨一格"，写成算子形式能让方程看起来干净很多。

工具层面，我后面统一用 Python，主力是numpy、pandas、statsmodels和matplotlib。statsmodels里有两套 AR 接口，statsmodels.tsa.ar_model.AutoReg是较新的那套，还有个老的AR类，我在第 5 节会说清楚该用哪个、为什么。

2. AR(1) 模型：结构、平稳条件与统计性质

2.1 模型表达式与"记忆"的数学形式

AR(1) 写成两种等价形式：

不带常数项：X_t = φX_{t-1} + ε_t

带常数项：X_t = c + φX_{t-1} + ε_t

两者只差一个均值。带常数项那一版，均值 μ = c/(1-φ)，你可以自己代进去验证：两边同时取期望，μ = c + φμ，移项就得到 μ = c/(1-φ)。这意味着常数项不是"平均水平"，而是控制均值的参数，它的实际含义要跟 φ 一起看。很多人第一次拟合带常数项的 AR(1)，看到截距估计值是 0.03 就以为序列均值是 0.03，其实真正的均值是 0.03/(1-0.8)=0.15，差了五倍。这个坑我在早期报告里踩过。

另外注意，ε_t 必须是白噪声，这条假设有两个含义：一是它跟所有历史 X 都不相关，二是它自身没有自相关。如果残差里还有结构，说明模型没提干净，要么阶数不够，要么该上 ARMA。

2.2 平稳条件为什么一定是 |φ| < 1

这一点值得认真推一遍，因为它是整个 AR 模型的地基。

用滞后算子把 AR(1) 写成 (1 - φL)X_t = ε_t，然后形式地解出：

X_t = (1 - φL)^{-1} ε_t = ε_t + φε_{t-1} + φ²ε_{t-2} + φ³ε_{t-3} + …

这是把 X_t 展开成过去所有扰动的一个加权和（一般称为 MA(∞) 表示）。要让它有意义，也就是让这个无穷级数收敛，需要 Σ|φ|^k < ∞，这等价于 |φ| < 1。

从另一个角度看更直观。反复迭代 X_t = φX_{t-1}+ε_t，能得到 X_t = φ^t X_0 + Σ_{i=0}^{t-1} φ^i ε_{t-i}。初始值的影响以 φ^t 的速度衰减，只有 |φ|<1 时这个影响才会消失，序列才会"忘掉"起点、收敛到同一个分布。φ=1 时序列就是随机游走，方差随时间线性增长，根本不平稳；|φ|>1 时初始扰动被无限放大，序列爆炸。

注意：平稳性条件是针对方程的参数，不是针对样本数据的表现。你完全可能生成一条看起来波动范围稳定的序列，但它的 φ 估计值是 1.02，那它长期是发散的，只是样本太短没看出来。

2.3 均值、方差与自相关函数的完整推导

在平稳前提下，AR(1) 的统计性质可以一行一行推出来，建议自己动手写一遍。

均值：两边取期望，μ = φμ，只要 φ≠1 就有 μ = 0（不带常数项）。这也是为什么很多教材干脆不提均值，默认序列已中心化。

自协方差：两边同乘 X_{t-k}（k≥1）再取期望：

E[X_t X_{t-k}] = φE[X_{t-1}X_{t-k}] + E[ε_t X_{t-k}]

左边是 γ_k，右边第一项是 φγ_{k-1}，第二项因为 ε_t 与过去不相关所以为 0。于是得到递归式：

γ_k = φγ_{k-1}，k ≥ 1

这个差分方程的解是 γ_k = φ^k γ_0。

方差：k=0 时特殊一点，两边平方取期望：γ_0 = φ²γ_0 + σ²，解得 γ_0 = σ²/(1-φ²)。注意这个式子又一次说明 |φ|<1 的必要性，φ→1 时方差趋于无穷。

自相关函数：把两式相除，得到极其干净的结果：

ρ_k = φ^k，k ≥ 0

这是 AR(1) 最标志性的性质：ACF 以几何速度衰减（"拖尾"）。φ>0 时单调衰减，φ<0 时正负交替衰减。

给你一个能记住的换算：衰减到一半需要的期数约为 ln(0.5)/ln(φ)。φ=0.8 时约 3.1 期，φ=0.9 时约 6.6 期，φ=0.95 时约 13.5 期。做业务预测时这个数字很有用——它告诉你"这个模型的有效预测视野大概有多远"。

2.4 偏自相关：为什么它在 1 阶之后就断了

ACF 是"X_t 和 X_{t-k} 的总相关"，这个总相关里包含了一部分"绕路"来的相关性。比如 ρ₂ 不为零，可能只是因为 X_t 跟 X_{t-1} 相关、X_{t-1} 又跟 X_{t-2} 相关，并不是 X_t 真的跟 X_{t-2} 有直接联系。

偏自相关函数（PACF）就是把这层中介效应剔除之后的"净相关"。计算 φ_kk 的方法是：用 X_t 对 X_{t-1},…,X_{t-k} 做回归，取最后一个自变量的系数。对 AR(1) 来说：

• k=1：回归只有一项，φ₁₁ = ρ₁ = φ
• k=2：用 X_t 对 X_{t-1}、X_{t-2} 做回归。因为真实过程里 X_{t-2} 的系数是 0，样本估计出来的 φ₂₂ 应该接近 0

k≥2 之后全部为 0，这叫"截尾"。所以 AR(1) 的识别特征很明确：ACF 拖尾，PACF 一阶截尾。

但"接近 0"到底多接近算 0？标准做法是看 2/√n 这条置信带。n=500 时带宽约 ±0.089。如果是复根导致的伪周期，PACF 可能会在几个滞后上看起来超出这个带宽，这时候别急着加阶，先看整体形状——拖着一个小尾巴跳来跳去，通常只是噪声。

3. AR(2) 模型：二阶滞后带来的新变化

3.1 模型形式与特征方程

把阶数加到 2：

X_t = φ₁X_{t-1} + φ₂X_{t-2} + ε_t

用滞后算子写成 (1 - φ₁L - φ₂L²)X_t = ε_t。这个二次多项式的根决定了序列的全部动态性质。

令 z 为特征方程 1 - φ₁z - φ₂z² = 0 的根，平稳性要求所有根的模大于 1（也就是落在单位圆外）。等价地，用 λ = 1/z 表示为 λ² - φ₁λ - φ₂ = 0，要求两根 λ₁、λ₂ 都在单位圆内（|λ|<1）。

这里我建议不要死记公式，而是理解结构：AR(1) 时特征方程是 1-φz=0，根是 1/φ，要求 |1/φ|>1 即 |φ|<1，跟前面的结论一致。AR(2) 只是把这个条件扩展到了二次情形。

3.2 平稳域是个三角形，怎么划出来的

把"两根都在单位圆内"翻译成对 φ₁、φ₂ 的约束，经过代数整理可以得到三条不等式：

• φ₁ + φ₂ < 1
• φ₂ - φ₁ < 1
• |φ₂| < 1

在 (φ₁, φ₂) 平面上，这三条线围出一个三角形，三个顶点分别是 (2, -1)、(-2, -1)、(0, 1)。这个三角区域就是 AR(2) 的平稳域。我每次做 AR(2) 建模，估计完参数之后都会做一件事：把点画到这张图上，看看有没有落在三角形外面。落到外面意味着模型理论上是发散的，哪怕残差检验看起来还过得去，也不能拿它做长期预测。

顺便说一句，这个检查非常便宜，比跑任何诊断检验都快。我在一次电力负荷的项目里就靠它发现了一个问题：数据做了错误的差分，导致 φ₂ 的估计值落在 -1.1 附近，稳稳地在三角形外面，说明差分过头了。

3.3 ACF 的两种形态与"伪周期"识别

AR(2) 的自相关同样满足 Yule-Walker 递归：

ρ_k = φ₁ρ_{k-1} + φ₂ρ_{k-2}，k ≥ 3

起始两个值是：

• ρ₁ = φ₁/(1-φ₂)
• ρ₂ = φ₁ρ₁ + φ₂ = φ₁²/(1-φ₂) + φ₂

拿个具体例子算一遍。φ₁=0.6、φ₂=0.3：

• ρ₁ = 0.6/0.7 ≈ 0.857
• ρ₂ = 0.6×0.857 + 0.3 ≈ 0.814
• ρ₃ = 0.6×0.814 + 0.3×0.857 ≈ 0.745
• ρ₄ ≈ 0.6×0.745 + 0.3×0.814 ≈ 0.691

单调衰减，形态跟 AR(1) 很像，只是衰减慢一点。此时判别式 φ₁²+4φ₂ = 0.36+1.2 = 1.56 > 0，特征根是实根。

换一组：φ₁=0.5、φ₂=-0.7。先验证平稳性：φ₁+φ₂ = -0.2 < 1 成立，φ₂-φ₁ = -1.2 < 1 成立，|φ₂|=0.7 < 1 成立，落在三角形内。再看判别式：φ₁²+4φ₂ = 0.25-2.8 = -2.55 < 0，复根。这时序列会呈现明显的振荡，周期约为

T = 2π / arccos( φ₁ / (2√(-φ₂)) )

代入数值：2√0.7 ≈ 1.673，φ₁/1.673 ≈ 0.299，arccos(0.299) ≈ 1.267 弧度，T ≈ 2π/1.267 ≈ 4.96，也就是大约 5 期一个循环。

这个公式在实际工作中很值钱。当你看到一条日频序列大概每周波动一次，有些人的本能是加一个 7 阶的季节项，但如果 φ₂ 为负且由它产生的伪周期恰好是 5 到 7 期，用 AR(2) 就能解释，模型参数量少得多，也更容易解释。当然反过来也成立：如果不算这个周期就瞎拟合，可能会把阶数定得过高。

提示：复根带来的振荡是"随机相位"的，波峰波谷的位置每次都不太一样，跟真正由外部因素驱动的严格季节性完全是两回事。判断方法是看这张 ACF 图的波峰是否每次都能对齐到同一个滞后位置，对不齐的基本就是 AR 内生振荡。

3.4 PACF 两阶截尾

AR(2) 的 PACF 在 k=1 和 k=2 处不为 0，k≥3 之后理论上全部为 0。具体表达式：

• φ₁₁ = ρ₁
• φ₂₂ = (ρ₂ - ρ₁²) / (1 - ρ₁²)

用上面那组 φ₁=0.6、φ₂=0.3 的数据算：ρ₁=0.857，ρ₂=0.814，代入得 φ₂₂ = (0.814 - 0.734)/(1 - 0.734) = 0.08/0.266 ≈ 0.30。有意思的是这个值恰好等于 φ₂，这不是巧合——AR(2) 的最后一个偏自相关理论上就等于 φ₂。你可以用这个性质做快速验算：估计出来的 φ₂ 和 PACF 第二阶的值应该对得上。

所以识别表长这样：

模型	ACF 形态	PACF 形态
AR(1)	几何衰减（拖尾）	1 阶后截尾
AR(2)	衰减/振荡（拖尾）	2 阶后截尾
AR(p)	拖尾	p 阶后截尾
MA(q)	q 阶后截尾	拖尾
ARMA(p,q)	拖尾	拖尾

这张表建议背下来，实际建模时第一眼的判断就靠它。但也要知道它的局限：当 p 或 q 大于 2、或者两个过程的特征互相抵消时，肉眼很难分清"截尾"和"拖尾"，最终还是要靠信息准则和残差检验落地。

4. 参数估计与定阶：从 Yule-Walker 到最小二乘

4.1 Yule-Walker 方程：把手推结果直接变成估计方法

把 AR(2) 的自相关递归式在 k=1、k=2 处写全：

• ρ₁ = φ₁ + φ₂ρ₁
• ρ₂ = φ₁ρ₁ + φ₂

这其实就是 Yule-Walker 方程的一个特例。一般形式用矩阵写出来是 Γφ = γ，其中 Γ 是用 ρ₁…ρ_{p-1} 构成的 Toeplitz 矩阵，γ 是 ρ₁…ρ_p 组成的向量。解出来：

φ = Γ⁻¹γ

对 AR(2) 有闭式解：

• φ₁ = ρ₁(1-ρ₂)/(1-ρ₁²)
• φ₂ = (ρ₂-ρ₁²)/(1-ρ₁²)

噪声方差同样有闭式：σ² = γ₀(1 - φ₁ρ₁ - φ₂ρ₂)。

Yule-Walker 估计的做法是：用样本自相关 r_k 替换理论值 ρ_k，直接代入公式。它计算快、一定落在平稳域内（这是它的一个独特优势），但小样本时偏差比较大，尤其是 φ 接近 1 或者序列接近非平稳时。另外注意，Yule-Walker 估计默认不带常数项，实际用的时候要先中心化数据。

4.2 最小二乘与极大似然：我更常用的两种

条件最小二乘（CLS）的思路是把 AR 模型当成一个普通回归：把 X_t 当因变量，X_{t-1},…,X_{t-p} 当自变量，前 p 个观测拿来做初始条件。损失函数是残差平方和：

Σ_{t=p+1}^{n} (X_t - c - φ₁X_{t-1} - φ₂X_{t-2})²

最小化它对参数求导，得到一组线性方程组，解起来很快。极大似然（MLE）在 CLS 基础上多考虑了初始观测的联合密度，理论上更优，小样本时差别比较明显。

statsmodels.tsa.ar_model.AutoReg用的是条件极大似然，参数顺序在输出里是const、x.L1、x.L2。这一点特别容易搞混，我第一次用的时候把const当成 φ₁ 填进了公式里，算出来的长期均值离谱到没法看。

在实践中我通常这么做：先用 Yule-Walker 快速拿一组初始值和对模型的直觉，正式报告里的参数则用极大似然的结果。两者差异很大时，说明样本可能有非平稳趋势或者异常值，需要先回去处理数据。

4.3 定阶：ACF/PACF 初判加信息准则复核

定阶我一般走两步。

第一步是图形判断。画 ACF 和 PACF，看 PACF 在第几阶掉到置信带以内。如果 PACF 第一阶高出很多、第二阶贴着带子边缘，那 AR(1) 和 AR(2) 都值得试；如果第二阶明显超出，就上 AR(2)。

第二步是算信息准则。对 p=1 到 pmax 逐个拟合，比较 AIC 或 BIC：

• AIC = -2lnL + 2k
• BIC = -2lnL + k·ln(n)

其中 k 是参数个数（含常数项和方差），n 是有效样本量。AIC 倾向于选稍大的模型，BIC 惩罚更重、倾向于更简洁的模型。我的经验是：样本量小于 200 时听 BIC，样本量大时两者差别不大，但如果它们打架，我会顺着 AIC 选的那个再看一眼残差，如果残差干净就用小的那个，毕竟简洁性本身有价值。

pmax 一般取 min(10, n/10) 之类的上限，别一上来就试到 50 阶。

4.4 参数估计的不确定性：别只看点估计

AR(1) 的参数估计近似标准误有个很好用的公式：se(φ̂) ≈ √((1-φ²)/n)。n=500、φ=0.8 时，se ≈ √(0.36/500) ≈ 0.027，95% 置信区间大约是 0.8±0.053。这说明什么？说明你要区分 φ=0.8 和 φ=0.9 这两个假设，n=500 根本不够——它们的差距远大于估计误差，但如果你只有 100 个样本，se 会放大到 0.06，那就更分不清了。

这个数量级的感觉在做项目时非常关键。有人拿着 60 个观测的月度数据建 AR(2)，估出来 φ₂=0.15，p 值 0.4，然后纠结"要不要保留这一项"。我的做法很直接：先算一下这个样本量能分辨的最小效应大概是多少，如果 0.15 本来就在噪声水平之下，纠结没有意义，直接用更简洁的模型。

5. Python 实操：把 AR(1) 和 AR(2) 完整跑一遍

5.1 生成可控的模拟数据

自己造数据是最好的学习方式，因为真实参数已知，估出来对不对一眼就知道。下面这段代码生成一个 AR(2) 过程：

import numpy as np import pandas as pd import matplotlib.pyplot as plt from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.ar_model import AutoReg from statsmodels.stats.diagnostic import acorr_ljungbox np.random.seed(42) n = 500 phi1, phi2 = 0.6, 0.3 eps = np.random.normal(0, 1.0, n) x = np.zeros(n) x[0] = eps[0] x[1] = phi1 * x[0] + eps[1] for t in range(2, n): x[t] = phi1 * x[t-1] + phi2 * x[t-2] + eps[t] series = pd.Series(x, name="ar2_sim") print(series.describe())

两点说明。一是前两个点的处理：严格来说需要从平稳分布抽样，我这里用了简单的初始条件，会让前几十个观测受初始值影响，所以后面做参数估计时我会把前 100 个点丢掉（burn-in）。二是噪声方差设成 1，理论方差应该是 σ²/(1-φ₁ρ₁-φ₂ρ₂)，代入 0.6、0.3 和前面算的 ρ₁=0.857、ρ₂=0.814，得到 1/(1-0.514-0.244) ≈ 4.13，标准差约 2.03。跑一下describe()，如果样本标准差在这个量级附近，说明生成逻辑没问题。

生成完先画图，这一步别省。折线图看趋势和异常值，直方图看分布形态。我见过不少人直接跳到拟合，结果数据里明明有个断崖式的跳变（设备故障导致的一小段异常），模型硬生生把它当成真实的动态吸收进参数里，后面所有结论都偏了。

5.2 看一眼 ACF 和 PACF

fig, axes = plt.subplots(2, 1, figsize=(10, 7)) plot_acf(series.iloc[100:], lags=30, ax=axes[0]) plot_pacf(series.iloc[100:], lags=30, ax=axes[1], method="ywm") plt.tight_layout() plt.show()

plot_pacf的method参数值得留意。老版本默认用的是 Yule-Walker 方法，新版本默认改成了ldbiased，在小样本下两者画出来的置信带不一样，容易让人误判。我一般显式写成method="ywm"，跟教材上的结论对得上。

你会看到 ACF 缓慢衰减且全程为正（因为两个系数都为正），PACF 在第一、二阶显著，第三阶开始落进置信带。跟第 3 节的推导吻合。

5.3 拟合、读参数、做残差检验

model = AutoReg(series.iloc[100:], lags=2, old_names=False) res = model.fit() print(res.params) print(res.bse) print(res.aic, res.bic)

输出里的参数顺序是const、ar.L1、ar.L2。真实值是 0.6 和 0.3，估计值一般会在 0.55~0.65 和 0.25~0.35 之间，样本量 400 的情况下这个波动是正常的。

残差检验我固定做三件事：

1. 时序图：残差应该围绕 0 波动，没有趋势、没有明显的方差聚集（波动聚集说明该上 GARCH 而不是加 AR 阶数）。
2. 残差的 ACF：理论上应该全部落在置信带内，也就是没有自相关残留。
3. Ljung-Box 检验：滞后取 10、20 阶，看 p 值是否大于 0.05。

lb = acorr_ljungbox(res.resid, lags=[10, 20], return_df=True) print(lb)

三项都通过，模型才算站得住。如果 Ljung-Box 在滞后 10 阶显著，说明还有结构没提干净，回去看 PACF 是不是隐约有第三个峰。

5.4 定阶的自动化流程

把 p=1 到 8 遍历一遍，用循环比 AIC/BIC：

rows = [] for p in range(1, 9): r = AutoReg(series.iloc[100:], lags=p, old_names=False).fit() rows.append({ "p": p, "aic": r.aic, "bic": r.bic, "lb_p_lag10": acorr_ljungbox(r.resid, lags=[10])["lb_pvalue"].values[0] }) df_cmp = pd.DataFrame(rows).sort_values("bic") print(df_cmp)

跑下来你会看到 AIC 和 BIC 大概率都指向 p=2，而且 p=2 的 Ljung-Box p 值最高。这个交叉验证的过程比单看一个准则可靠得多。这里要多说一句：AutoReg默认每个 p 都带常数项，如果你比较不同 p 的 AIC，参数个数必须算一致，包括常数项和方差参数，否则比较没有意义。

5.5 预测与预测区间

fc = res.get_prediction(start=len(series.iloc[100:]), end=len(series.iloc[100:])+19) fc_df = fc.summary_frame(alpha=0.05) print(fc_df.head())

AR(2) 的长期预测会收敛到均值（这里是 0），而且预测区间的宽度会逐步趋于 γ₀ 的开方，也就是序列自身的标准差。这个性质有个实际含义：如果业务方问"你预测三个月后是多少"，诚实一点的回答是"我只能说它大概落在均值附近这个范围里"——点预测本身的信息量随步长迅速衰减。我通常会把预测区间一并给出去，而不是只给一条线。

另外一个实用细节：用get_prediction而不是手动递推，前者会自动处理置信区间的计算，包括参数估计不确定性带来的额外方差。手动递推虽然简单，但容易漏掉这一块，让区间看着比实际窄。

6. 踩坑记录与常见问题排查

6.1 常见问题速查表

现象	可能原因	处理方式
估计出的 φ 落在平稳域外	数据未差分、趋势未去、异常值干扰	先做差分和趋势剔除，检查异常点
PACF 第 3、4 阶也冒出峰	阶数不够，或存在未建模的季节性	提高 pmax 重试，或改用 SARIMA
ACF 衰减极慢（几十阶都不落）	接近单位根	做单位根检验，考虑差分
残差 Ljung-Box 显著	阶数不足或模型选择错误	加阶、试 MA 项，检查是否该用 ARMA
参数估计值波动大、标准误高	样本量太小或严重共线性	缩短滞后阶数，延长样本区间
长期预测发散	参数在外但没检查	把 (φ₁,φ₂) 画到平稳三角形里核对

这张表是从我自己和同事的调试记录里整理出来的，实际遇到问题时照着走一遍，八成能定位到原因。

6.2 三个必须提醒的坑

第一个坑：用statsmodels老接口AR类时容易忽略的默认行为。

老接口statsmodels.tsa.ar_model.AR默认用 Yule-Walker 方法，且不带常数项，参数索引方式也跟新接口不一样。如果你的代码是从某篇 2018 年以前的博客抄来的，很可能踩这个坑，估出来的值跟理论差一截还找不到原因。建议统一用AutoReg，并且显式写上lags和old_names=False（后者只影响参数显示名称，不影响数值）。

第二个坑：数据预处理对 AR 参数的影响非常大。

AR 模型是建立在平稳假设上的，但真实数据的平稳往往需要先处理：去趋势、去均值、处理异常值、必要时做一阶差分。每一步都会改变 AR 参数的物理含义。举个例子，原始序列做一阶差分之后，AR(1) 的系数往往变成负的，这不是模型坏了，而是"差分"这个操作本身引入了负相关。你要是拿差分后的参数去解释原始序列的惯性，结论会完全反过来。

第三个坑：AR 阶数不能凭"业务上有几个影响因素"来决定。

我见过有人因为业务上说"影响销量的是上周价格和老爹上周的促销"，就直接上 AR(2)，结果残差一堆自相关。滞后阶数是统计意义上的记忆长度，跟业务因素个数没有对应关系。影响因素该放 AR 还是放外生变量（也就是 ARX 模型），得看它是滞后观测还是当期可得的解释变量。这个区分搞错了，模型的可解释性再漂亮也没用。

6.3 我个人的几条经验

第一条，永远先画图再拟合。折线图、ACF、PACF 三张图加起来不到十行代码，却能省掉后面几个小时的返工。

第二条，参数估计之后一定要跑平稳性检查。把估计出的 (φ₁, φ₂) 标到那张三角形图上，这个动作花不了 30 秒，但它能拦住一批"看起来拟合得挺好、实际会发散"的模型。

第三条，样本量决定了你能分辨多细的结构。前面算过，AR(1) 参数的标准误大概是 √((1-φ²)/n)，这个公式可以反过来用：想要把标准误压到 0.03 以内，在 φ=0.8 的情况下大概需要 400 个有效观测。做数据收集方案的时候，这个数字比拍脑袋估个"至少几百条"靠谱得多。

最后分享一个小技巧。当你不确定该用 AR(1) 还是 AR(2) 时，可以两边都拟合，然后比较残差方差。如果 AR(2) 相对 AR(1) 的残差方差下降不到 5%，我基本就不加那一阶了。参数量增加带来的过拟合风险和解释成本，通常超过那点微弱的拟合改善。这个 5% 不是硬标准，但它比单纯盯 p 值更贴近实际业务判断——因为业务方真正关心的是模型能不能用，而不是某个系数在统计上是不是显著。