平稳时序模型识别与参数估计

平稳序列建模：模型识别与参数估计

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

当我们确认一个时间序列是平稳非白噪声序列后，下一步就是为其建立合适的模型。本节将系统讲解平稳时间序列建模的核心流程：如何根据自相关图与偏自相关图识别模型类型（AR、MA或ARMA）并确定阶数，以及如何通过矩估计、极大似然估计和最小二乘估计等方法精确估计模型参数。掌握这些方法，你就能为任何平稳序列构建出可靠的统计模型，为后续的预测与分析奠定基础。

💡 核心导读

• 建模全流程：从平稳非白噪声序列出发，遵循“模型识别→参数估计→模型检验→模型优化→序列预测”的完整路径。
• 模型识别与定阶：利用样本自相关图（ACF）和偏自相关图（PACF）的“截尾”与“拖尾”特性，判断应使用AR(p)、MA(q)还是ARMA(p, q)模型，并初步确定阶数p和q。
• 参数估计三剑客：理解矩估计、极大似然估计和最小二乘估计三种核心方法的原理、优缺点及适用场景。
• Python实战：使用statsmodels库一站式完成从序列生成、模型定阶、参数估计到显著性检验的全过程。

平稳序列建模的完整步骤

一个严谨的平稳时间序列建模过程包含以下七个步骤，它们构成了一个循环迭代、不断优化的闭环：

1. 平稳性与白噪声检验：确保序列是平稳的且包含可供建模的信息（非白噪声）。
2. 计算样本相关系数：计算并绘制样本自相关系数（ACF）和样本偏自相关系数（PACF）图。
3. 模型识别：观察ACF和PACF图的特征，初步判断模型类型（AR/MA/ARMA）并估计阶数（p, q）。
4. 参数估计：使用估计方法（如极大似然估计）计算出模型中的未知参数（如$\phi_i$, $\theta_j$, $\sigma^2$）。
5. 模型检验：检验残差序列是否为白噪声（模型显著性检验），以及各个参数是否显著不为零（参数显著性检验）。若检验不通过，则返回第3步重新识别模型。
6. 模型优化：在多个通过检验的模型中，根据AIC、BIC等准则选择最优模型。
7. 序列预测：使用最终确定的模型进行未来值的预测。

本节将重点讲解第3步（模型识别）和第4步（参数估计）。

模型识别：看图定阶的艺术

模型识别的核心依据是AR、MA、ARMA三类模型的理论性质在其样本自相关图（ACF）和样本偏自相关图（PACF）上的体现。

理论基础：三类模型的相关性特征

截尾：相关系数在延迟若干阶（如p或q阶）后，突然衰减至零值附近，并在后续各阶在零附近做小值随机波动。
拖尾：相关系数逐渐衰减至零（可能呈指数衰减、正弦振荡衰减），衰减过程连续、缓慢。

实战定阶：如何判断“截尾”？

由于样本的随机性，真实的样本ACF/PACF图不会呈现完美的理论截尾。我们通常借助95%的置信区间作为判断工具。

根据样本相关系数的渐近分布理论，对于平稳序列，当样本量$n$足够大时，有： $$\hat{\rho}_k \sim N\left(0, \frac{1}{n}\right), \quad \hat{\phi}_{kk} \sim N\left(0, \frac{1}{n}\right)$$ 因此，样本自相关系数$\hat{\rho}_k$和样本偏自相关系数$\hat{\phi}_{kk}$的95%置信区间近似为： $$\left( -\frac{2}{\sqrt{n}}, \ \frac{2}{\sqrt{n}} \right)$$

经验判断法则： 如果样本（偏）自相关系数在最初的$d$阶明显超出两倍标准差范围（即落在置信区间外），而之后几乎95%的系数都落在置信区间内，并且从非零值衰减到区间内小值波动的过程非常突然，则通常可视为 $d$阶截尾。

案例分析

案例1：北京市城乡居民定期储蓄比例序列 (1950-1998)

• ACF图：延迟3阶后全部落入置信区间，但衰减过程连续缓慢 → 拖尾。
• PACF图：延迟1阶的系数显著超出置信区间，之后各阶系数均在区间内小值波动，衰减突然 → 1阶截尾。
• 识别结论：ACF拖尾，PACF 1阶截尾，符合AR(1) 模型特征。

案例2：美国加油站OVERSHORT序列 (57天)

• ACF图：延迟1阶系数超出置信区间，之后各阶均在区间内 → 1阶截尾。
• PACF图：系数在置信区间内外交替波动，衰减缓慢 → 拖尾。
• 识别结论：ACF 1阶截尾，PACF拖尾，符合MA(1) 模型特征。

案例3：全球气温改变值差分序列 (1880-1985)

• ACF图：系数缓慢衰减，在置信区间内外均有出现 → 拖尾。
• PACF图：系数同样缓慢衰减，呈现拖尾特征 → 拖尾。
• 识别结论：ACF与PACF均拖尾，提示为ARMA(p,q)混合过程（p,q>0）；可从最简的ARMA(1,1) 开始尝试，并与(1,2)、(2,1)等候选一起用后续AIC/BIC准则择优。

参数估计：从识别到量化

模型识别确定了模型的“骨架”（类型与阶数），参数估计则是为其填充“血肉”（具体的系数值）。对于一个中心化的ARMA(p, q)模型： $$X_t = \phi_1 X_{t-1} + \dots + \phi_p X_{t-p} + \varepsilon_t - \theta_1 \varepsilon_{t-1} - \dots - \theta_q \varepsilon_{t-q}$$ 其中 $\varepsilon_t \sim WN(0, \sigma^2)$。待估参数包括：自回归系数 $\phi_1, \dots, \phi_p$，移动平均系数 $\theta_1, \dots, \theta_q$，以及白噪声方差 $\sigma^2$，共计 $p+q+1$ 个。

方法一：矩估计 (Method of Moments, MOM)

核心思想：用样本矩（如样本均值、样本方差、样本自相关系数）去估计相应的总体矩，并通过总体矩与待估参数之间的理论关系（如Yule-Walker方程）解出参数。

以AR(2)模型为例： 模型为 $X_t = \phi_1 X_{t-1} + \phi_2 X_{t-2} + \varepsilon_t$。其Yule-Walker方程为： $$\begin{cases} \rho_1 = \phi_1 + \phi_2 \rho_1 \\ \rho_2 = \phi_1 \rho_1 + \phi_2 \end{cases}$$ 矩估计的步骤是：

1. 用样本数据计算样本自相关系数 $\hat{\rho}_1$ 和 $\hat{\rho}_2$。
2. 将 $\hat{\rho}_1$, $\hat{\rho}_2$ 代入上述方程，解得： $$\hat{\phi}_1 = \frac{\hat{\rho}_1(1-\hat{\rho}_2)}{1-\hat{\rho}_1^2}, \quad \hat{\phi}_2 = \frac{\hat{\rho}_2 - \hat{\rho}_1^2}{1-\hat{\rho}_1^2}$$

优点：思想直观，计算简单，无需假设总体分布。 缺点：只用了前$p+q$个样本自相关系数的信息，信息浪费严重，估计精度通常较差。常作为其他迭代估计方法的初始值。

方法二：极大似然估计 (Maximum Likelihood Estimation, MLE)

核心思想：在已知样本观测值的条件下，寻找能使该样本出现概率最大的那组参数值。

通俗理解（摸球实验）：有两个箱子，A箱有10白球90黑球，B箱有90白球10黑球。你随机摸5次，结果4黑1白。你更倾向于认为球来自哪个箱子？显然是A箱，因为该结果在A箱出现的概率远高于B箱。MLE就是基于这种“最可能”的原则。

数学表述：设待估参数向量为 $\boldsymbol{\beta} = (\phi_1, \dots, \phi_p, \theta_1, \dots, \theta_q, \sigma^2)^T$，样本为 $\mathbf{X} = (X_1, \dots, X_n)^T$。在假设$\varepsilon_t$服从正态分布的条件下，可以写出样本的联合密度函数（似然函数）$L(\boldsymbol{\beta} | \mathbf{X})$。极大似然估计 $\hat{\boldsymbol{\beta}}_{MLE}$ 就是使得 $L(\boldsymbol{\beta} | \mathbf{X})$ 达到最大的参数值。实践中常最大化其对数形式（对数似然函数）以简化计算。

优点：充分利用了所有样本信息，估计精度高，且具有一致性、渐近正态性等优良统计性质。 缺点：需要假设总体分布（通常为正态分布），计算复杂，常需迭代求解。

方法三：最小二乘估计 (Least Squares Estimation, LSE)

核心思想：寻找能使残差平方和最小的那组参数值。 $$\hat{\boldsymbol{\beta}}_{LSE} = \arg\min_{\boldsymbol{\beta}} \sum_{t=1}^n (X_t - \hat{X}_t)^2$$ 其中 $\hat{X}_t$ 是基于模型和参数对 $X_t$ 的拟合值。

条件最小二乘估计 (Conditional Sum of Squares, CSS)：是实践中最常用的变体。它假定初始时刻之前的扰动项$\varepsilon_0, \varepsilon_{-1}, \dots$均为0，在此条件下进行最小二乘估计。这简化了计算，且当样本量较大时，估计结果与MLE非常接近。

优点：充分利用样本信息，估计精度高，CSS方法使用广泛。 缺点：为了确保估计量的优良性质，通常也需要假设扰动项服从正态分布。

Python实战：一站式建模演示

下面我们使用Python的statsmodels库，完整演示从生成平稳序列、模型识别（观察ACF/PACF）、到参数估计（矩估计、MLE）的全过程。

"""
时间序列分析 3.4 课时配套测试代码
内容：平稳序列参数估计 (矩估计 / 极大似然估计 MLE / 条件最小二乘 CSS)
"""

import numpy as np
import pandas as pd
from statsmodels.tsa.arima_process import ArmaProcess
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.regression.linear_model import yule_walker

np.random.seed(42)
n = 300

# 1. 真实参数: ARMA(1, 1) -> phi=0.7, theta=-0.4, sigma=1.0
ar_true = np.array([1, -0.7])
ma_true = np.array([1, -0.4])
process = ArmaProcess(ar_true, ma_true)
data = pd.Series(process.generate_sample(nsample=n), name='series')

# 2. 矩估计 (Yule-Walker 估计 AR(1))
phi_yw, sigma_yw = yule_walker(data, order=1, method='mle')

# 3. 极大似然估计 (MLE 拟合 ARMA(1,1))
model_mle = ARIMA(data, order=(1, 0, 1))
res_mle = model_mle.fit()

# 4. 提取参数与显著性检验
param_series = pd.Series(res_mle.params, index=res_mle.param_names)
bse_series = pd.Series(res_mle.bse, index=res_mle.param_names)
pval_series = pd.Series(res_mle.pvalues, index=res_mle.param_names)

print(f"真实参数: phi = 0.7000, theta = -0.4000")
print(f"Yule-Walker AR(1) 矩估计: phi = {phi_yw[0]:.4f}, sigma = {sigma_yw:.4f}")
print(f"ARIMA(1,0,1) MLE 估计结果:")
for name in res_mle.param_names:
    print(f"  {name}: {param_series[name]:.4f} (标准误: {bse_series[name]:.4f}, p值: {pval_series[name]:.4e})")
print(f"对数似然值 (Log-Likelihood): {res_mle.llf:.2f}")

运行结果与解读：

真实参数: phi = 0.7000, theta = -0.4000
Yule-Walker AR(1) 矩估计: phi = 0.2693, sigma = 0.9923
ARIMA(1,0,1) MLE 估计结果:
  const: -0.0082 (标准误: 0.1063, p值: 9.3868e-01)
  ar.L1: 0.7428 (标准误: 0.1204, p值: 6.9292e-10)
  ma.L1: -0.5229 (标准误: 0.1547, p值: 7.2374e-04)
  sigma2: 0.9585 (标准误: 0.0693, p值: 1.4783e-43)
对数似然值 (Log-Likelihood): -419.40

关键结论：

1. 平稳性确认：该过程由构造参数 $\phi=0.7$（$|\phi|<1$）已知平稳、$\theta=-0.4$（$|\theta|<1$）已知可逆，也可自行调用 adfuller 进一步验证。
2. 模型误设的代价：将ARMA(1,1)序列误判为AR(1)并用矩估计，得到的 $\hat{\phi}=0.2693$ 严重偏离真实值0.7，这凸显了正确识别模型的重要性。
3. 极大似然估计的精度：MLE估计的 $\hat{\phi}=0.7428$ 非常接近真实值0.7。转换后的理论 $\hat{\theta}=0.5229$ 也接近真实值0.4，且两者的p值均极显著（<0.001），说明估计准确可靠。
4. 软件符号约定：务必注意不同软件/教材对MA系数符号的定义。statsmodels输出的ma.L1对应其MA多项式 $1 + \theta_{sm}L$ 中的系数。若理论模型为 $X_t = \phi X_{t-1} + \varepsilon_t - \theta \varepsilon_{t-1}$，则 $\theta_{sm} = -\theta$。本案例中，ma.L1 = -0.5229，故理论 $\hat{\theta} = 0.5229$。

📝 动手练一练

1. 模型识别练习：观察以下描述，判断最可能适合的模型类型（AR/MA/ARMA）及阶数。

   • ACF图在延迟1、2阶显著非零，3阶及以后在零附近小范围波动；PACF图呈指数衰减。
   • ACF图呈正弦波式振荡衰减；PACF图在延迟1阶显著非零，2阶及以后在零附近小范围波动。

2. 估计方法理解：为什么在实际应用中，矩估计很少作为最终的参数估计方法，但又不可或缺？

参考答案：

1. (1) ACF 2阶后截尾，PACF拖尾 → MA(2)模型。 (2) ACF拖尾，PACF 1阶后截尾 → AR(1)模型。
2. 矩估计精度较差，因为它只使用了有限阶的样本自相关信息。但它计算快速、无需分布假设，常被用作极大似然估计或最小二乘估计等迭代算法的初始值，帮助算法更快、更稳定地收敛到最优解。

本章小结

本节我们深入探讨了平稳时间序列建模的核心两步：模型识别与参数估计。

• 模型识别依赖于对样本自相关图（ACF）和偏自相关图（PACF）的解读，利用AR、MA、ARMA三类模型独特的“截尾”与“拖尾”性质进行判断。95%的置信区间是辅助判断“截尾”点的实用工具。
• 参数估计有三种主要方法：
  • 矩估计：简单直观，但精度低，多用作迭代算法的初始值。
  • 极大似然估计 (MLE)：精度高，性质优良，是主流方法，但需假设分布且计算复杂。
  • 最小二乘估计 (LSE/CSS)：精度高，特别是条件最小二乘（CSS）在实践中广泛应用。
• Python实战展示了使用statsmodels库一站式完成平稳性检验、图形识别、以及多种参数估计方法的完整流程，并特别提醒了不同工具间MA系数符号的差异。

行动清单

1. 练习读图：找几个真实或模拟的平稳序列，绘制其ACF和PACF图，尝试判断模型类型与阶数。
2. 代码复现：运行本节提供的Python代码，理解每个步骤的输出，并尝试修改真实参数（ar_true, ma_true），观察估计结果的变化。
3. 对比估计方法：对同一序列，分别用method='mle'和method='css'拟合ARIMA模型，比较参数估计值、对数似然值、AIC/BIC的差异。

— 小象教研组