自回归模型（AR 模型）

小象实战讲义 · 时间序列分析

在上一讲中，我们掌握了如何检验一个时间序列是否平稳、是否具有非随机的内在规律。本讲将聚焦于平稳且非白噪声序列的核心建模方法——自回归模型（AR）。我们将深入理解 AR 模型的定义、平稳性条件、统计性质及其识别特征。学完本节，你将能够判断一个平稳序列是否适合用 AR 模型拟合，并理解其自相关与偏自相关图所蕴含的模型信息。

💡 核心导读

• AR 模型定义：理解 p 阶自回归模型 AR(p) 的数学结构，掌握其中心化变换与延迟算子表示。
• 平稳性判别：掌握 AR 模型平稳的充要条件——特征根在单位圆内，并了解平稳域的概念。
• 统计性质：推导平稳 AR 模型的均值、方差、协方差，并理解其自相关系数（拖尾）与偏自相关系数（p 阶截尾）的核心特征。
• 模型识别：学会通过观察样本自相关图（ACF）与偏自相关图（PACF）来初步识别 AR 模型。
• Python 实战：使用 statsmodels 库模拟、检验并分析一个平稳的 AR(2) 过程。

方法性工具回顾

在正式进入 AR 模型前，我们先回顾三个关键的数学工具，它们将极大简化后续模型的表达与推导。

差分运算

差分是处理非平稳序列（如趋势、季节）的基础运算。

• 一阶差分：$\nabla X_t = X_t - X_{t-1}$。
• p 阶差分：$\nabla^p X_t = \nabla^{p-1} X_t - \nabla^{p-1} X_{t-1}$，即连续差分 p 次。
• k 步差分：$\nabla_k X_t = X_t - X_{t-k}$，即间隔 k 期的观测值相减。

延迟算子

延迟算子 $B$ 是一个“时间指针”，其作用是将序列值的时间向前拨动。

• 定义：$B X_t = X_{t-1}$，进而有 $B^k X_t = X_{t-k}$。
• 性质：延迟算子满足与常数类似的运算律，如结合律、分配律等，这使得我们可以像处理多项式一样处理它。
• 表示差分：利用延迟算子，差分运算可以简洁地表示为：
  • p 阶差分：$\nabla^p X_t = (1 - B)^p X_t$
  • k 步差分：$\nabla_k X_t = (1 - B^k) X_t$

线性差分方程

AR 模型的求解与线性差分方程理论紧密相关。

• 齐次方程：形如 $a_0 z_t + a_1 z_{t-1} + \dots + a_p z_{t-p} = 0$ 的方程。
• 求解：通过构建特征方程 $a_0 \lambda^p + a_1 \lambda^{p-1} + \dots + a_p = 0$，求解特征根 $\lambda_1, \dots, \lambda_p$，进而得到通解。方程的稳定性（解不发散）要求特征根在单位圆内（$|\lambda_i| < 1$），这与 AR 模型的平稳性条件一致。

AR 模型的定义与结构

自回归模型的核心思想是：序列当前值 $X_t$ 可以表示为自身过去若干期值 $X_{t-1}, X_{t-2}, \dots, X_{t-p}$ 的线性组合，再加上一个随机扰动。

标准定义

具有如下结构的模型称为 p 阶自回归模型，简记为 $AR(p)$：

$$\begin{cases} X_t = \phi_0 + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \dots + \phi_p X_{t-p} + \varepsilon_t \\ \phi_p \neq 0 \\ \varepsilon_t \sim WN(0, \sigma_{\varepsilon}^2) \\ E(\varepsilon_t X_s) = 0, \quad \forall s < t \end{cases}$$

其中：

• $\phi_0$ 为常数项（截距）。
• $\phi_1, \phi_2, \dots, \phi_p$ 为自回归系数，且最高阶系数 $\phi_p \neq 0$。
• $\varepsilon_t$ 是均值为 0、方差为 $\sigma_{\varepsilon}^2$ 的白噪声序列，且与过去的所有观测值 $X_s (s < t)$ 不相关。

中心化 AR 模型

若序列均值 $\mu = E(X_t)$ 不为零，我们可以通过中心化变换，得到一个均值为零的序列，从而简化模型。

1. 计算序列均值：对于平稳 $AR(p)$ 模型，其均值 $\mu = \frac{\phi_0}{1 - \phi_1 - \dots - \phi_p}$。
2. 定义中心化序列：$Y_t = X_t - \mu$。
3. 得到中心化 $AR(p)$ 模型： $$Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \dots + \phi_p Y_{t-p} + \varepsilon_t$$ 此时，模型不再包含常数项 $\phi_0$。

延迟算子表示

引入延迟算子 $B$，中心化 $AR(p)$ 模型可以写成极其紧凑的形式：

$$\Phi(B) Y_t = \varepsilon_t$$

其中，$\Phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \dots - \phi_p B^p$ 称为 自回归系数多项式。这个表示法将模型的核心结构清晰地封装在一个多项式里。

AR 模型的平稳性判别

并非所有形式上符合 $AR(p)$ 定义的模型都对应平稳序列。模型的平稳性完全由其自回归系数 $\phi_1, \dots, \phi_p$ 决定。

判别原因与示例

考虑以下四个模型：

1. $X_t = 0.8 X_{t-1} + \varepsilon_t$
2. $X_t = -1.1 X_{t-1} + \varepsilon_t$
3. $X_t = X_{t-1} - 0.5 X_{t-2} + \varepsilon_t$
4. $X_t = X_{t-1} + 0.5 X_{t-2} + \varepsilon_t$

通过模拟这些模型生成的时序图可以发现，模型 (1) 和 (3) 的序列波动平稳，而模型 (2) 的方差随时间爆炸式增长，模型 (4) 则呈现出明显的趋势。因此，(2) 和 (4) 是不平稳的 AR 模型。

平稳性判别方法

1. 特征根判别法（常用）：

   • 条件：$AR(p)$ 模型平稳的 充要条件 是其 p 个特征根都在单位圆内（即模长 $|\lambda_i| < 1$）。
   • 特征根：是自回归系数多项式 $\Phi(B)=0$ 的根的倒数。等价地，要求 $\Phi(B)=0$ 的根都在 单位圆外（模长 $>1$）。

2. 平稳域判别法：

   • 对于低阶模型，可以直接给出系数 $\phi_i$ 的取值范围（平稳域）。
   • $AR(1)$：平稳域为 $|\phi_1| < 1$。
   • $AR(2)$：平稳域为 $\phi_2 + \phi_1 < 1, \quad \phi_2 - \phi_1 < 1, \quad |\phi_2| < 1$。

应用上述判别法可知，示例中：

• 模型(1)：$\phi_1=0.8$，满足 $|\phi_1|<1$，平稳。
• 模型(2)：$\phi_1=-1.1$，不满足 $|\phi_1|<1$，不平稳。
• 模型(3)：$\phi_1=1, \phi_2=-0.5$，代入 $AR(2)$ 平稳域验证，平稳。
• 模型(4)：$\phi_1=1, \phi_2=0.5$，不满足平稳域条件，不平稳。

平稳 AR 模型的统计性质

对于满足平稳性条件的 $AR(p)$ 模型，我们可以推导出其一系列重要的统计性质，这些性质是模型识别和参数估计的基础。

均值与方差

• 均值：对于中心化模型，$E(Y_t) = 0$。对于非中心化模型，$E(X_t) = \mu = \frac{\phi_0}{1 - \sum_{i=1}^p \phi_i}$。
• 方差：平稳 AR 模型的方差是一个常数。求解需要借助 格林函数 (Green Function)。模型可以表示为无限阶的移动平均形式（传递形式）： $$Y_t = \sum_{j=0}^{\infty} G_j \varepsilon_{t-j}$$ 其中 $G_j$ 是格林函数，满足特定递推公式。方差则为： $$\gamma_0 = Var(Y_t) = \sigma_{\varepsilon}^2 \sum_{j=0}^{\infty} G_j^2$$ 以 $AR(1)$ 模型 $Y_t = \phi_1 Y_{t-1} + \varepsilon_t$ 为例，其方差为 $\gamma_0 = \frac{\sigma_{\varepsilon}^2}{1 - \phi_1^2}$。

自协方差与自相关系数

• 自协方差函数 $\gamma_k$：度量 $Y_t$ 与 $Y_{t-k}$ 的协方差。对于平稳 $AR(p)$ 模型，它满足 Yule-Walker 方程： $$\gamma_k = \phi_1 \gamma_{k-1} + \phi_2 \gamma_{k-2} + \dots + \phi_p \gamma_{k-p}, \quad k \ge 1$$
• 自相关系数 $\rho_k$：$\rho_k = \frac{\gamma_k}{\gamma_0}$。将 Yule-Walker 方程两边同除以 $\gamma_0$，即得到自相关系数的递推公式： $$\rho_k = \phi_1 \rho_{k-1} + \phi_2 \rho_{k-2} + \dots + \phi_p \rho_{k-p}, \quad k \ge 1$$

AR 模型自相关系数的核心性质是拖尾性：$\rho_k$ 不会在有限阶后突然截断为 0，而是随着 $k$ 增大，按负指数或正弦衰减的方式逐渐趋向于 0。衰减的模式（单调、振荡、伪周期）由特征根决定。

偏自相关系数

偏自相关系数 $\phi_{kk}$ 度量的是在给定中间 $k-1$ 个随机变量 $Y_{t-1}, \dots, Y_{t-k+1}$ 的条件下，$Y_t$ 与 $Y_{t-k}$ 之间的相关性。

AR 模型偏自相关系数的核心性质是 p 阶截尾性：对于 $AR(p)$ 模型，有 $$\phi_{kk} = 0, \quad \forall k > p$$ 这意味着，在滞后阶数 $k$ 大于模型阶数 $p$ 之后，偏自相关系数在理论上应为零。这一性质是我们在实际中识别 AR 模型阶数 $p$ 的关键依据。

Python 实战：AR 过程模拟与性质验证

下面，我们使用 Python 模拟一个平稳的 $AR(2)$ 过程，并验证其平稳性条件及自相关特征。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller

# 设置随机种子，保证结果可复现
np.random.seed(42)
n = 300  # 序列长度

# 1. 定义并模拟一个平稳的 AR(2) 过程
# 模型：X_t = 0.5*X_{t-1} + 0.3*X_{t-2} + e_t
# 滞后算子多项式系数：[1, -0.5, -0.3] 代表 1 - 0.5B - 0.3B^2
ar_params = np.array([1, -0.5, -0.3])  # 注意：statsmodels 使用此格式
ma_params = np.array([1])              # 移动平均部分为 1 (即只有白噪声)
ar_process = ArmaProcess(ar_params, ma_params)

# 检验平稳性 (内部使用特征根判别)
is_stationary = ar_process.isstationary
roots = ar_process.arroots
root_magnitudes = np.abs(roots)

# 模拟生成序列数据
simulated_data = ar_process.generate_sample(nsample=n)

print("=== AR(2) 过程平稳性检验 ===")
print(f"模型平稳性判据 (isstationary): {is_stationary}")
print(f"特征多项式根: {np.round(roots, 4).tolist()}")
print(f"特征根模长 (|root| > 1 则平稳): {np.round(root_magnitudes, 4).tolist()}")
print(f"\n=== 模拟序列描述性统计 ===")
print(f"前 5 项观测值: {np.round(simulated_data[:5], 4).tolist()}")
print(f"序列样本均值: {simulated_data.mean():.4f}")
print(f"序列样本方差: {simulated_data.var():.4f}")

# 2. ADF 单位根检验 (另一种平稳性检验)
adf_result = adfuller(simulated_data)
print(f"\n=== ADF 单位根检验结果 ===")
print(f"ADF 统计量: {adf_result[0]:.4f}")
print(f"P-value: {adf_result[1]:.4f}")
print(f"结论: {'序列平稳 (P<0.05)' if adf_result[1] < 0.05 else '序列可能非平稳'}")

# 3. 绘制时序图、ACF 和 PACF
fig, axes = plt.subplots(3, 1, figsize=(10, 8))

# 时序图
axes[0].plot(simulated_data)
axes[0].set_title('模拟 AR(2) 序列时序图')
axes[0].set_xlabel('时间 t')
axes[0].set_ylabel('X_t')
axes[0].grid(True, alpha=0.3)

# 自相关图 (ACF)
plot_acf(simulated_data, lags=40, ax=axes[1], title='样本自相关图 (ACF) - 呈现拖尾性')
axes[1].grid(True, alpha=0.3)

# 偏自相关图 (PACF)
plot_pacf(simulated_data, lags=40, ax=axes[2], title='样本偏自相关图 (PACF) - 2阶后截尾')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

运行上述代码，你将得到类似以下的输出和图表：

=== AR(2) 过程平稳性检验 ===
模型平稳性判据 (isstationary): True
特征多项式根: [-2.8403, 1.1736]
特征根模长 (|root| > 1 则平稳): [2.8403, 1.1736]

=== 模拟序列描述性统计 ===
前 5 项观测值: [0.4967, 0.1101, 0.8517, 1.9819, 1.0123]
序列样本均值: -0.0525
序列样本方差: 1.8365

=== ADF 单位根检验结果 ===
ADF 统计量: -5.6814
P-value: 0.0000
结论: 序列平稳 (P<0.05)

结果解读：

1. 平稳性：isstationary 为 True，且两个特征根的模长均大于 1（等价于 $\Phi(B)=0$ 的根在单位圆外），ADF 检验 p 值远小于 0.05，三者均确认该 AR(2) 过程是平稳的。
2. 统计性质：模拟序列的样本均值和方差在固定值附近波动，符合平稳序列特征。
3. 图形特征：
   • 时序图：序列围绕均值上下波动，无明显趋势或周期性。
   • ACF 图：自相关系数随着滞后阶数增加逐渐衰减至 0 附近，呈现典型的“拖尾”特征。
   • PACF 图：偏自相关系数在滞后 2 阶之后，其值几乎全部落在置信区间内（即与 0 无显著差异），呈现出明显的“2 阶截尾”特征。这正是 $AR(2)$ 模型的标志性特征。

📝 动手练一练

1. 平稳性判断：对于 $AR(1)$ 模型 $X_t = \phi X_{t-1} + \varepsilon_t$，若 $\phi = 1.05$，该模型是否平稳？请分别用特征根判别法和直观理解（如对方差的影响）解释。
2. 模型识别：假设你拿到一个平稳序列，绘制其 ACF 和 PACF 图后，发现 ACF 拖尾，而 PACF 在滞后 3 阶后截尾。你认为最适合拟合该序列的 ARMA 族模型是什么？阶数如何？

参考答案：

1. 不平稳。特征根判别法：$\phi = 1.05$，其模 $|\phi| > 1$，不满足 $AR(1)$ 平稳条件 $|\phi| < 1$。直观上，将模型递归展开，$X_t$ 的方差会包含 $(1.05^{2t})$ 项，导致方差随时间指数增长，序列发散。
2. 最适合的模型是 $AR(3)$。因为 PACF 的 p 阶截尾性是 AR(p) 模型的独有特征，ACF 拖尾也符合 AR 模型的性质。PACF 在 3 阶后截尾，提示自回归阶数 $p=3$。

本章小结

本节我们深入探讨了自回归模型（AR）的核心内容：

• 定义与表示：理解了 $AR(p)$ 模型用当期与过去 p 期值的线性关系来解释序列动态，并掌握了其中心化形式和简洁的延迟算子表示 $\Phi(B)X_t = \varepsilon_t$。
• 平稳性：认识到并非所有 AR 模型都平稳，其平稳的充要条件是自回归系数多项式的根在单位圆外（或特征根在单位圆内）。我们学习了特征根判别法和平稳域判别法。
• 统计性质：掌握了平稳 AR 模型的均值、方差公式，并重点理解了其 自相关系数 (ACF) 拖尾 和 偏自相关系数 (PACF) p 阶截尾 这两个最关键的识别特征。
• Python 实战：通过代码模拟了一个平稳的 $AR(2)$ 过程，验证了其平稳性，并直观观察到了 ACF 拖尾和 PACF 在 2 阶后截尾的典型图形特征。

行动清单

1. 概念自查：能否不看书，独立写出 $AR(p)$ 模型的定义、平稳性条件、ACF 和 PACF 的核心性质？
2. 代码复现：运行讲义中的 Python 代码，尝试修改 ar_params 为不平稳的参数（如 [1, -1.2, -0.3]），重新运行并观察时序图、ACF/PACF 图发生了何种变化。
3. 图形识别：在网上或教科书中寻找一些真实时间序列的 ACF/PACF 图，尝试判断哪些可能适合用 AR 模型拟合，并估计其可能的阶数 p。

— 小象教研组