平稳时间序列建模代码实战

平稳时间序列建模实战：从理论到 Python 实现

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

在掌握了平稳时间序列分析的理论基础后，本节将带领大家将理论知识转化为实践能力。我们将学习如何使用 Python 语言实现 AR、MA、ARMA 等经典时序模型的识别、估计、检验与预测全流程，并通过真实数据案例演示如何从时序图、自相关图出发，一步步构建出最优的预测模型。

💡 核心导读

• 模型识别原理：掌握如何通过自相关图（ACF）与偏自相关图（PACF）的”截尾”与”拖尾”特征，初步判断 AR、MA 或 ARMA 模型类型及其阶数。
• Python 建模全流程：学习使用 statsmodels 库完成平稳序列的模型拟合、参数估计、模型显著性检验（残差白噪声检验）与参数显著性检验。
• 模型优化与选择：理解 AIC 与 BIC 准则在模型选择中的作用，学会在多个有效模型中选择相对最优的模型。
• 序列预测实战：掌握点预测与区间预测的实现方法，理解 AR、MA 模型预测行为的本质差异（如 MA 模型的有限记忆性）。
• 稀疏系数模型：了解如何通过约束特定滞后项的系数为零，构建更精简的模型结构。

平稳时间序列建模流程回顾

平稳时间序列建模遵循一套严谨的流程，确保模型的可靠性与预测的有效性。完整的建模步骤如下图所示：

flowchart TD
    A[平稳非白噪声序列] --> B[计算样本自相关<br>与偏自相关系数]
    B --> C{模型识别}
    C --> D[AR(p)模型]
    C --> E[MA(q)模型]
    C --> F[ARMA(p,q)模型]
    D --> G[参数估计]
    E --> G
    F --> G
    G --> H{模型检验}
    H -- 通过 --> I[模型优化]
    H -- 未通过 --> C
    I --> J[序列预测]

建模始于对序列平稳性与纯随机性的检验。只有平稳且非白噪声的序列才值得进一步建模。随后，我们通过分析样本自相关系数（ACF）和偏自相关系数（PACF）的特征来识别可能的模型类型与阶数。

模型识别：ACF 与 PACF 的”语言”

模型识别的核心依据是 AR(p)、MA(q) 和 ARMA(p,q) 模型自相关与偏自相关函数的理论性质：

截尾：相关系数在延迟若干阶后突然衰减至零值附近，并保持在置信区间内波动。 拖尾：相关系数按负指数或正弦振荡等方式逐渐衰减至零，不会在有限阶后截断。

理论回顾：ARMA 模型定义

ARMA(p,q) 模型具有以下结构：

$$\begin{cases} X_t = \phi_0 + \phi_1 X_{t-1} + \phi_2 X_{t-2} + \cdots + \phi_p X_{t-p} + \varepsilon_t - \theta_1 \varepsilon_{t-1} - \theta_2 \varepsilon_{t-2} - \cdots - \theta_q \varepsilon_{t-q} \\ \phi_p \neq 0, \quad \theta_q \neq 0 \\ E(\varepsilon_t) = 0, \quad Var(\varepsilon_t) = \sigma^2, \quad E(\varepsilon_t \varepsilon_s) = 0, s \neq t \end{cases}$$

其中 $\varepsilon_t \sim WN(0, \sigma^2)$ 为白噪声序列。特别当 $\phi_0 = 0$ 时，称为中心化 ARMA(p,q) 模型。

引入延迟算子 $B$，令 $\Phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p$，$\Theta(B) = 1 - \theta_1 B - \theta_2 B^2 - \cdots - \theta_q B^q$，则中心化模型可简写为：

$$\Phi(B)X_t = \Theta(B)\varepsilon_t$$

AR(p) 和 MA(q) 模型是 ARMA(p,q) 模型的特例：

• 当 $q = 0$ 时，ARMA(p,0) 即为 AR(p) 模型
• 当 $p = 0$ 时，ARMA(0,q) 即为 MA(q) 模型

Python 实战：从模拟到真实数据建模

1. 模拟验证：理论与实践的桥梁

我们首先通过模拟数据验证理论性质。使用 Python 可以方便地生成符合特定 ARMA 过程的序列，并观察其 ACF/PACF 图是否与理论一致。

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.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings('ignore')

# 设置随机种子确保可重复性
np.random.seed(42)

# 1. 模拟 MA(2) 过程：ε_t ~ N(0,1), θ1=0.5, θ2=0.2
print("=== 模拟 MA(2) 过程验证 ===")
ma_coef = np.array([0.5, 0.2])  # MA系数
ar_coef = np.array([1])         # AR部分为1（无AR项）
process_ma = ArmaProcess(ar_coef, np.r_[1, ma_coef])
ma_data = pd.Series(process_ma.generate_sample(nsample=5000))

# 拟合 MA(2) 模型
model_ma = ARIMA(ma_data, order=(0, 0, 2))
result_ma = model_ma.fit()
print(f"MA(2) 拟合参数: ma.L1 = {result_ma.params['ma.L1']:.4f}, "
      f"ma.L2 = {result_ma.params['ma.L2']:.4f}")
print(f"理论参数: θ1 = 0.5000, θ2 = 0.2000")
print(f"残差方差估计: {result_ma.params['sigma2']:.4f} (理论值: 1.0000)")
print(f"AIC值: {result_ma.aic:.4f}")

# 2. 模拟 AR(2) 过程：φ1=0.5, φ2=0.2
print("\n=== 模拟 AR(2) 过程验证 ===")
ar_coef = np.array([0.5, 0.2])
process_ar = ArmaProcess(np.r_[1, -ar_coef], [1])
ar_data = pd.Series(process_ar.generate_sample(nsample=5000))

# 拟合 AR(2) 模型
model_ar = ARIMA(ar_data, order=(2, 0, 0))
result_ar = model_ar.fit()
print(f"AR(2) 拟合参数: ar.L1 = {result_ar.params['ar.L1']:.4f}, "
      f"ar.L2 = {result_ar.params['ar.L2']:.4f}")
print(f"理论参数: φ1 = 0.5000, φ2 = 0.2000")
print(f"残差方差估计: {result_ar.params['sigma2']:.4f} (理论值: 1.0000)")
print(f"AIC值: {result_ar.aic:.4f}")

# 3. 模拟 ARMA(2,2) 过程
print("\n=== 模拟 ARMA(2,2) 过程验证 ===")
ar_coef = np.array([0.5, 0.2])
ma_coef = np.array([0.2, 0.3])
process_arma = ArmaProcess(np.r_[1, -ar_coef], np.r_[1, ma_coef])
arma_data = pd.Series(process_arma.generate_sample(nsample=5000))

# 拟合 ARMA(2,2) 模型
model_arma = ARIMA(arma_data, order=(2, 0, 2))
result_arma = model_arma.fit()
print(f"ARMA(2,2) 拟合参数:")
print(f"  AR部分: ar.L1 = {result_arma.params['ar.L1']:.4f}, ar.L2 = {result_arma.params['ar.L2']:.4f}")
print(f"  MA部分: ma.L1 = {result_arma.params['ma.L1']:.4f}, ma.L2 = {result_arma.params['ma.L2']:.4f}")
print(f"理论参数: φ1=0.5, φ2=0.2, θ1=0.2, θ2=0.3")
print(f"残差方差估计: {result_arma.params['sigma2']:.4f} (理论值: 1.0000)")

运行结果分析： 从输出可以看到，即使使用 5000 个样本点，拟合参数与理论真值仍存在微小偏差，这是参数估计的固有特性。但偏差很小，说明估计方法是有效的。样本量越大，估计通常越精确。

2. ACF/PACF 图可视化验证

# 可视化对比 AR(2) 与 MA(2) 的 ACF/PACF 特征
fig, axes = plt.subplots(2, 2, figsize=(12, 8))

# AR(2) 的 ACF 和 PACF
plot_acf(ar_data, ax=axes[0, 0], lags=30, title='AR(2) 自相关图 (ACF)')
plot_pacf(ar_data, ax=axes[0, 1], lags=30, title='AR(2) 偏自相关图 (PACF)')

# MA(2) 的 ACF 和 PACF
plot_acf(ma_data, ax=axes[1, 0], lags=30, title='MA(2) 自相关图 (ACF)')
plot_pacf(ma_data, ax=axes[1, 1], lags=30, title='MA(2) 偏自相关图 (PACF)')

plt.tight_layout()
plt.show()

图形解读：

• AR(2) 序列：ACF 呈现拖尾特征（逐渐衰减），PACF 在 2 阶后明显截尾（落入蓝色置信带内）。
• MA(2) 序列：ACF 在 2 阶后明显截尾，PACF 呈现拖尾特征。

这与理论完全一致，验证了通过 ACF/PACF 识别模型类型的可行性。

3. 模型预测特性对比

AR 模型与 MA 模型的预测行为有本质区别：

# MA(2) 模型预测
forecast_ma = result_ma.get_forecast(steps=6)
pred_ma = forecast_ma.predicted_mean
print(f"MA(2) 未来6期预测值:\n{pred_ma}")

# AR(2) 模型预测  
forecast_ar = result_ar.get_forecast(steps=6)
pred_ar = forecast_ar.predicted_mean
print(f"\nAR(2) 未来6期预测值:\n{pred_ar}")

关键发现：

• MA(2) 预测：2 期之后的预测值收敛到序列均值（本例中为 0），这是因为 MA 模型只有有限记忆性，超过其阶数 q 后，预测不再依赖历史信息。
• AR(2) 预测：所有期数的预测值都不同，AR 模型具有无限记忆性，通过自回归结构将历史信息不断传递到未来预测中。

真实数据建模实战

1. 数据准备与初步分析

# 读取真实数据（这里用模拟数据替代，实际应用时替换为真实数据文件）
# 假设我们有一个化学过程测量序列
np.random.seed(123)
n_points = 80
# 生成一个具有 AR(1) 特性的序列作为示例数据
real_data = np.zeros(n_points)
real_data[0] = np.random.randn()
for t in range(1, n_points):
    real_data[t] = 0.7 * real_data[t-1] + np.random.randn()
real_data = real_data + 50  # 添加常数均值

real_series = pd.Series(real_data, name='chemical_process')

# 时序图观察
plt.figure(figsize=(10, 4))
plt.plot(real_series, 'b-', linewidth=1)
plt.title('化学过程测量值时序图')
plt.xlabel('时间')
plt.ylabel('测量值')
plt.grid(True, alpha=0.3)
plt.show()

2. 平稳性与纯随机性检验

from statsmodels.tsa.stattools import adfuller

# ADF 平稳性检验
adf_result = adfuller(real_series)
print(f"ADF 检验统计量: {adf_result[0]:.4f}")
print(f"p-value: {adf_result[1]:.4f}")
if adf_result[1] < 0.05:
    print("结论: 序列平稳 (拒绝原假设)")
else:
    print("结论: 序列非平稳 (不能拒绝原假设)")

# Ljung-Box 纯随机性检验
lb_test = acorr_ljungbox(real_series, lags=12, return_df=True)
print(f"\nLjung-Box 检验 (滞后12阶):")
print(f"统计量: {lb_test['lb_stat'].iloc[-1]:.4f}, p-value: {lb_test['lb_pvalue'].iloc[-1]:.4f}")
if lb_test['lb_pvalue'].iloc[-1] < 0.05:
    print("结论: 序列非白噪声 (拒绝原假设，适合建模)")
else:
    print("结论: 序列为白噪声 (无需进一步建模)")

3. 模型识别与定阶

# 绘制 ACF 和 PACF 图
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
plot_acf(real_series, ax=axes[0], lags=30, title='自相关图 (ACF)')
plot_pacf(real_series, ax=axes[1], lags=30, title='偏自相关图 (PACF)')
plt.tight_layout()
plt.show()

# 基于图形特征初步判断：
# 1. ACF 拖尾衰减 -> 可能为 AR 或 ARMA 模型
# 2. PACF 1阶后截尾 -> 可能为 AR(1) 模型

4. 模型拟合与比较

# 尝试拟合不同模型
models = {}
results = {}

# AR(1) 模型
models['AR(1)'] = ARIMA(real_series, order=(1, 0, 0))
results['AR(1)'] = models['AR(1)'].fit()

# MA(2) 模型 (假设我们看到 ACF 2阶截尾)
models['MA(2)'] = ARIMA(real_series, order=(0, 0, 2))
results['MA(2)'] = models['MA(2)'].fit()

# ARMA(1,1) 模型
models['ARMA(1,1)'] = ARIMA(real_series, order=(1, 0, 1))
results['ARMA(1,1)'] = models['ARMA(1,1)'].fit()

# 比较模型效果
comparison = pd.DataFrame({
    'Model': list(results.keys()),
    'AIC': [results[m].aic for m in results],
    'BIC': [results[m].bic for m in results],
    'Log-Likelihood': [results[m].llf for m in results]
})
print("模型比较:")
print(comparison.sort_values('AIC'))

# 选择 AIC 最小的模型
best_model_name = comparison.loc[comparison['AIC'].idxmin(), 'Model']
print(f"\n根据 AIC 准则，最优模型为: {best_model_name}")
best_result = results[best_model_name]
print(f"\n{best_model_name} 模型参数:")
print(best_result.summary())

5. 模型检验

# 残差分析
residuals = best_result.resid

# 残差时序图
plt.figure(figsize=(10, 6))
plt.subplot(2, 2, 1)
plt.plot(residuals, 'b-', linewidth=0.8)
plt.axhline(y=0, color='r', linestyle='--', alpha=0.5)
plt.title('残差序列时序图')
plt.xlabel('时间')
plt.ylabel('残差')

# 残差 ACF 图
plt.subplot(2, 2, 2)
plot_acf(residuals, ax=plt.gca(), lags=20, title='残差自相关图')

# 残差 PACF 图
plt.subplot(2, 2, 3)
plot_pacf(residuals, ax=plt.gca(), lags=20, title='残差偏自相关图')

# 残差正态性检验 (Q-Q图)
plt.subplot(2, 2, 4)
from scipy import stats
stats.probplot(residuals, dist="norm", plot=plt)
plt.title('残差 Q-Q 图')

plt.tight_layout()
plt.show()

# Ljung-Box 检验残差是否为白噪声
lb_resid = acorr_ljungbox(residuals, lags=[10, 15, 20], return_df=True)
print("残差白噪声检验 (Ljung-Box):")
for lag in [10, 15, 20]:
    pval = lb_resid.loc[lb_resid['lb_stat'].index == lag, 'lb_pvalue'].values[0]
    print(f"  滞后{lag}阶: p-value = {pval:.4f}", 
          "-> 白噪声" if pval > 0.05 else "-> 非白噪声")

6. 序列预测与可视化

# 进行未来6期预测
forecast_steps = 6
forecast = best_result.get_forecast(steps=forecast_steps)
pred_mean = forecast.predicted_mean
pred_ci = forecast.conf_int(alpha=0.05)  # 95% 置信区间

# 可视化预测结果
plt.figure(figsize=(12, 6))

# 绘制历史数据
history_len = len(real_series)
plt.plot(range(history_len), real_series, 'b-', linewidth=1.5, label='历史数据')

# 绘制预测值
pred_index = range(history_len, history_len + forecast_steps)
plt.plot(pred_index, pred_mean, 'r-', linewidth=2, label='点预测')

# 绘制置信区间
plt.fill_between(pred_index, 
                  pred_ci.iloc[:, 0], 
                  pred_ci.iloc[:, 1], 
                  color='r', alpha=0.2, label='95% 置信区间')

plt.axvline(x=history_len-0.5, color='k', linestyle='--', alpha=0.7)
plt.title(f'{best_model_name} 模型预测结果')
plt.xlabel('时间')
plt.ylabel('测量值')
plt.legend(loc='best')
plt.grid(True, alpha=0.3)
plt.xlim(0, history_len + forecast_steps)
plt.show()

print(f"未来{forecast_steps}期预测值:")
for i in range(forecast_steps):
    print(f"  第{i+1}期: {pred_mean.iloc[i]:.4f} "
          f"[{pred_ci.iloc[i, 0]:.4f}, {pred_ci.iloc[i, 1]:.4f}]")

7. 稀疏系数模型（约束参数）

在某些情况下，我们可能希望某些滞后项的系数为零，构建更简洁的模型：

# 示例：拟合 AR(3) 模型，但强制第二个滞后项系数为0
# 注意：statsmodels 的 ARIMA 不支持直接约束参数，但可以通过以下方式近似实现

# 方法：先拟合完整模型，然后检查系数显著性
model_ar3 = ARIMA(real_series, order=(3, 0, 0))
result_ar3 = model_ar3.fit()
print("完整 AR(3) 模型参数:")
print(result_ar3.params)
print(f"\n参数 t 检验 p-values:")
print(result_ar3.pvalues)

# 如果 ar.L2 不显著 (p-value > 0.05)，可考虑拟合 AR(1)+AR(3) 模型
# 实际中可通过重新设定 order 或使用更灵活的模型设定

实战：ARMA 模型预测全流程演示

以下代码完整演示了从数据生成、模型拟合到预测评估的全过程：

"""
时间序列分析 3.6 课时配套实战代码
内容：ARMA模型预测全流程 (点预测、条件方差递推与95%预测置信区间构建)
"""

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

np.random.seed(42)
n_train = 200
n_test = 5
n_total = n_train + n_test

# 1. 模拟生成完整序列
ar_true = np.array([1, -0.7])
ma_true = np.array([1, 0.3])
process = ArmaProcess(ar_true, ma_true)
all_data = process.generate_sample(nsample=n_total)

train_series = pd.Series(all_data[:n_train], name='train')
true_future = all_data[n_train:]

# 2. 训练模型并进行 5 步外推预测
model = ARIMA(train_series, order=(1, 0, 1))
res = model.fit()

forecast_res = res.get_forecast(steps=n_test)
pred_mean = forecast_res.predicted_mean
conf_int = forecast_res.conf_int(alpha=0.05)  # 95% 置信区间
se_mean = forecast_res.se_mean

print(f"历史训练集样本数: {n_train}, 外推预测步数: {n_test}")
print(f"拟合参数: ar.L1 = {res.params['ar.L1']:.4f}, ma.L1 = {res.params['ma.L1']:.4f}, sigma2 = {res.params['sigma2']:.4f}")
print("\n外推预测与 95% 置信区间:")
for i in range(n_test):
    step = i + 1
    p_val = pred_mean.iloc[i]
    lower = conf_int.iloc[i, 0]
    upper = conf_int.iloc[i, 1]
    se = se_mean.iloc[i]
    actual = true_future[i]
    in_interval = lower <= actual <= upper
    print(f"  Step {step}: 预测值 = {p_val:.4f}, 标准误 = {se:.4f}, 95%区间 = [{lower:.4f}, {upper:.4f}], 真实值 = {actual:.4f}, 命中 = {in_interval}")

# 计算预测均方根误差 (RMSE)
rmse = np.sqrt(np.mean((pred_mean.values - true_future) ** 2))
print(f"\n5 步外推 RMSE: {rmse:.4f}")

运行结果：

历史训练集样本数: 200, 外推预测步数: 5
拟合参数: ar.L1 = 0.6454, ma.L1 = 0.3037, sigma2 = 0.8581

外推预测与 95% 置信区间:
  Step 1: 预测值 = -1.3314, 标准误 = 0.9263, 95%区间 = [-3.1471, 0.4842], 真实值 = -0.9992, 命中 = True
  Step 2: 预测值 = -0.9172, 标准误 = 1.2771, 95%区间 = [-3.4203, 1.5859], 真实值 = -0.0313, 命中 = True
  Step 3: 预测值 = -0.6498, 标准误 = 1.3975, 95%区间 = [-3.3888, 2.0892], 真实值 = 1.2294, 命中 = True
  Step 4: 预测值 = -0.4773, 标准误 = 1.4446, 95%区间 = [-3.3087, 2.3542], 真实值 = 2.2393, 命中 = True
  Step 5: 预测值 = -0.3659, 标准误 = 1.4638, 95%区间 = [-3.2350, 2.5031], 真实值 = 0.5060, 命中 = True

5 步外推 RMSE: 1.5853

结果解读：

1. 所有真实值都落在 95% 预测区间内，说明区间预测是有效的。
2. 随着预测步数增加，预测标准误逐渐增大，置信区间变宽，这符合预测不确定性随时间增大的规律。
3. RMSE 为 1.5853，提供了预测精度的整体度量。

📝 动手练一练

1. 模型识别练习：使用 statsmodels 生成一个 ARMA(1,1) 过程，其中 $\phi_1=0.6$, $\theta_1=-0.3$。绘制其 ACF 和 PACF 图，观察是否呈现理论上的”双拖尾”特征。然后尝试分别拟合 AR(1)、MA(1) 和 ARMA(1,1) 模型，比较它们的 AIC 值，验证哪个模型最优。

2. 预测对比实验：生成一个 MA(2) 序列（$\theta_1=0.5, \theta_2=0.2$）和一个 AR(2) 序列（$\phi_1=0.5, \phi_2=0.2$），分别进行 10 步预测。观察并解释：

   • MA(2) 模型在 2 期后的预测值有何特点？
   • AR(2) 模型的长期预测趋势如何？
   • 两种模型的预测置信区间变化规律有何不同？

参考答案：

1. ARMA(1,1) 的 ACF 和 PACF 都应呈现拖尾特征。理论上，ARMA(1,1) 模型的 AIC 值应该最小，因为它才是数据生成的真实模型。
2. MA(2) 在 2 期后预测值收敛到序列均值，置信区间宽度趋于稳定；AR(2) 预测值会逐渐收敛到长期均值，但每期预测值都不同，置信区间宽度持续增加。

本章小结

本节我们完成了从时间序列建模理论到 Python 实战的跨越，掌握了：

核心要点回顾：

1. 模型识别：通过 ACF/PACF 的截尾/拖尾特征初步判断模型类型（AR/MA/ARMA）及阶数。
2. 参数估计：使用 statsmodels.tsa.arima.model.ARIMA 进行模型拟合，支持极大似然估计等多种方法。
3. 模型检验：通过残差白噪声检验（Ljung-Box）验证模型提取信息是否充分；通过参数 t 检验判断模型是否最简。
4. 模型优化：使用 AIC/BIC 准则在多个有效模型中选择相对最优模型。
5. 序列预测：掌握点预测与区间预测方法，理解 AR 与 MA 模型预测行为的本质差异。

行动清单：

1. 对任意平稳非白噪声序列，尝试完整的建模流程：平稳性检验 → 纯随机性检验 → ACF/PACF 分析 → 模型拟合 → 模型检验 → 预测。
2. 使用 AIC/BIC 准则比较至少 3 个不同模型的拟合效果，选择最优模型并解释原因。
3. 对选定的最优模型进行 5-10 期预测，并绘制包含置信区间的预测图。

— 小象教研组