当前位置：首页 > news >正文

进阶篇第 4 篇：驾驭季节性波动 - SARIMA 模型实战

news 来源：原创 2025/7/20 18:16:03

进阶篇第 4 篇：驾驭季节性波动 - SARIMA 模型实战

Wavy Patterns or Seasonal Landscape
(图片来源: Pixabay on Pexels)

在上一篇中，我们深入探索了经典的 ARIMA(p,d,q) 模型。它通过整合自回归 (AR)、差分 (I) 和移动平均 (MA) 提供了一个强大的框架来对（处理后）平稳的时间序列进行建模。然而，我们也观察到，对于像 CO2 浓度这样具有明显季节性波动的数据，仅仅使用 ARIMA 可能不足以完全捕捉其动态，模型的残差中可能仍会留下明显的季节性模式。

那么，如何让 ARIMA 模型直接拥抱并处理这种周期性的起伏呢？答案就是它的扩展版本——SARIMA 模型，全称季节性差分整合移动平均自回归模型 (Seasonal Autoregressive Integrated Moving Average)。

本篇，我们将一起：

解构 SARIMA(p,d,q)(P,D,Q)s： 揭开季节性参数的神秘面纱。
升级建模流程： 学习如何在 Box-Jenkins 方法中融入季节性考量。
实战演练： 使用 Python (statsmodels) 对具有季节性的数据构建、诊断和预测 SARIMA 模型。
对比与提升： 看看 SARIMA 相对于简单 ARIMA 在处理季节性数据上的优势。

准备好让你的时间序列模型“感知”季节的变化了吗？

解构 SARIMA(p,d,q)(P,D,Q)s 模型

SARIMA 模型在 ARIMA(p,d,q) 的基础上，增加了一组季节性参数 (P,D,Q)s，专门用来描述数据在季节性周期上的相关性结构。

(p, d, q): 这部分与 ARIMA 模型完全相同，描述的是序列的短期、非季节性相关性。
- p: 非季节性 AR 阶数
- d: 非季节性差分次数
- q: 非季节性 MA 阶数
(P, D, Q)s: 这部分是新增的，描述的是序列的长期、季节性相关性。
- s (Seasonality): 季节性周期长度。这是最关键的参数之一，必须根据数据特点指定（例如，对于月度数据年度季节性，s=12；对于季度数据年度季节性，s=4；对于日数据周季节性，s=7）。
- P (Seasonal AR Order): 季节性自回归阶数。表示当前值与过去 P 个季节周期的对应值（例如 y(t-s), y(t-2s), …）之间的关系。
- D (Seasonal Difference Order): 季节性差分次数。表示为了消除季节性而需要进行的差分次数。通常 D 取 0 或 1。季节性差分计算为 y'(t) = y(t) - y(t-s)。
- Q (Seasonal MA Order): 季节性移动平均阶数。表示当前值与过去 Q 个季节周期的对应预测误差 (ε(t-s), ε(t-2s), …) 之间的关系。

直观理解： SARIMA 模型认为，一个时间点的值不仅与它紧邻的前几个点有关（由 (p,d,q) 描述），也与去年（或上个周期）同一时间点附近的值和误差有关（由 (P,D,Q)s 描述）。它通过两个 ARIMA 类的结构（一个作用于短期，一个作用于季节性滞后）巧妙地结合起来。

升级建模流程：融入季节性考量

SARIMA 的建模流程仍然遵循 Box-Jenkins 的思想（识别 -> 估计 -> 诊断 -> 预测），但在识别阶段需要额外考虑季节性：

步骤 1：模型识别 (Identification - 升级版)

目标： 确定非季节性差分 d、季节性差分 D，以及初步判断非季节性阶数 (p, q) 和季节性阶数 (P, Q)。
操作：
- a. (非季节性) 平稳性检验与差分 (定 d): 与 ARIMA 相同，通过 ADF 检验等确定需要的普通差分次数 d。得到差分 d 次后的序列 y'。
- b. 季节性检验与差分 (定 D):
  - 观察 y' 的 ACF/PACF 图在季节性滞后位置 (s, 2s, 3s, …) 是否有显著的尖峰或缓慢衰减。如果 ACF 在这些位置缓慢衰减，通常暗示需要进行季节性差分。
  - 进行一次季节性差分 (D=1)，得到 y'' = y'(t) - y'(t-s)。通常 D=1 就足够了。
  - 再次检验 y'' 的平稳性，并观察其 ACF/PACF 图，看季节性模式是否已消除。
- c. 非季节性阶数识别 (定 p, q):
  - 观察最终差分后的序列 y'' 的 ACF/PACF 图的前几个滞后 (例如 lag 1, 2, 3…)。
  - 根据截尾/拖尾模式初步判断非季节性的 p 和 q。
- d. 季节性阶数识别 (定 P, Q):
  - 观察最终差分后的序列 y'' 的 ACF/PACF 图在季节性滞后位置 (s, 2s, …) 是否有显著的尖峰。
  - 如果 ACF 在 lag s 处截尾，PACF 在 s, 2s, … 处拖尾，可能暗示需要 季节性 MA(1)，即 Q=1, P=0。
  - 如果 PACF 在 lag s 处截尾，ACF 在 s, 2s, … 处拖尾，可能暗示需要 季节性 AR(1)，即 P=1, Q=0。
  - 如果都在 s 处不显著，可能 P=0, Q=0。
  - 季节性阶数 P, Q 通常不会很大（常见 0 或 1）。
- 记录： 记录下你确定的 d, D, s 以及候选的 (p,q) 和 (P,Q) 组合。

步骤 2：参数估计 (Estimation)

使用 statsmodels.tsa.statespace.SARIMAX 模块来拟合模型，需要同时提供 order=(p,d,q) 和 seasonal_order=(P,D,Q,s) 参数。

步骤 3：模型诊断 (Diagnostic Checking)

与 ARIMA 完全相同：检查最终模型的残差是否为白噪声（残差图、ACF/PACF 图、Ljung-Box 检验）。如果残差仍有结构（尤其是在季节性滞后处），需要返回步骤 1 调整阶数。

步骤 4：预测 (Forecasting)

使用拟合好的 SARIMAX 模型进行预测。

Python 实战：构建 SARIMA 模型

我们再次使用月度 CO2 数据，它具有明显的上升趋势和年度季节性 (s=12)。

import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.stats.diagnostic import acorr_ljungbox# --- 1. 数据准备与识别 ---
co2_data = sm.datasets.co2.load_pandas().data['co2']
co2_data = co2_data.interpolate().resample('M').mean()# a. 非季节性差分 (定 d)
# 从上一篇已知，d=1 可能使趋势平稳
co2_d1 = co2_data.diff().dropna()# b. 季节性检验与差分 (定 D, s=12)
s = 12
# 绘制差分一次后数据的 ACF/PACF，观察季节性滞后
# fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# plot_acf(co2_d1, ax=axes[0], lags=40, title='ACF after 1 non-seasonal diff')
# plot_pacf(co2_d1, ax=axes[1], lags=40, title='PACF after 1 non-seasonal diff', method='ywm')
# plt.suptitle('Check for Seasonality after d=1', fontsize=16)
# plt.show()
# 预期：ACF 在 lag 12, 24 处可能仍然很高且缓慢衰减，提示需要 D=1co2_d1_D1 = co2_d1.diff(s).dropna() # 进行一次季节性差分 D=1
d=1
D=1print(f"--- Data after d={d} non-seasonal and D={D} seasonal differencing (s={s}) ---")
# 检验最终差分后数据的平稳性 (可选)
adf_test = sm.tsa.stattools.adfuller(co2_d1_D1)
print(f'ADF p-value after differencing: {adf_test[1]:.3f}') # 应远小于 0.05# c & d. 识别 p, q, P, Q
# 绘制最终差分后数据的 ACF/PACF 图
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
plot_acf(co2_d1_D1, ax=axes[0], lags=40, title='ACF after d=1, D=1')
plot_pacf(co2_d1_D1, ax=axes[1], lags=40, title='PACF after d=1, D=1', method='ywm')
plt.suptitle(f'ACF/PACF after d={d}, D={D}, s={s}', fontsize=16)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plt.show()
# 根据图形判断:
# 短期滞后 (lag 1, 2...): 可能 ACF 在 lag 1 截尾，PACF 拖尾 -> 尝试 q=1, p=0 ?
# 季节性滞后 (lag 12, 24...): 可能 ACF 在 lag 12 截尾，PACF 拖尾 -> 尝试 Q=1, P=0 ?
# 结合起来，一个可能的候选模型是 SARIMA(0,1,1)(0,1,1,12)
# 这只是一个示例判断，真实分析可能需要尝试更多组合# --- 2. 参数估计 ---
# 选定一个候选模型，例如 SARIMA(1,1,1)(1,1,1,12)
p, q = 1, 1
P, Q = 1, 1# 注意：模型是在原始数据上拟合，内部处理差分
model = SARIMAX(co2_data,order=(p, d, q),seasonal_order=(P, D, Q, s),enforce_stationarity=False, # Let model handle differencingenforce_invertibility=False)
model_fit = model.fit(disp=False) # disp=False to turn off verbose outputprint(f"\n--- Fitting SARIMA({p},{d},{q})({P},{D},{Q},{s}) Model ---")
print(model_fit.summary())# --- 3. 模型诊断 ---
print("\n--- Model Diagnostics ---")
residuals = model_fit.resid
fig = model_fit.plot_diagnostics(figsize=(12, 8))
plt.tight_layout()
plt.show()lb_test = acorr_ljungbox(residuals, lags=[2*s], return_df=True) # 检验到 2 倍季节周期
print("\nLjung-Box Test for Residual Autocorrelation:")
print(lb_test)
# 检查 p-value 是否 > 0.05# --- 4. 预测 ---
n_forecast = 36 # 预测未来 3 年
forecast_result = model_fit.get_forecast(steps=n_forecast)
forecast_values = forecast_result.predicted_mean
confidence_intervals = forecast_result.conf_int()plt.figure(figsize=(12, 6))
plt.plot(co2_data.index, co2_data, label='Observed')
# plt.plot(co2_data['2000':].index, co2_data['2000':], label='Observed (Zoomed)') # 可选：放大看近期
plt.plot(forecast_values.index, forecast_values, label='SARIMA Forecast', color='red')
plt.fill_between(confidence_intervals.index,confidence_intervals.iloc[:, 0],confidence_intervals.iloc[:, 1],color='pink', alpha=0.5, label='95% Confidence Interval')
plt.title(f'SARIMA({p},{d},{q})({P},{D},{Q},{s}) Forecast for CO2 Data')
plt.xlabel('Date')
plt.ylabel('CO2 Concentration')
plt.legend()
plt.show()

在这里插入图片描述

解读与对比：

最终差分后的 ACF/PACF: 观察季节性差分 (D=1) 是否有效消除了季节性滞后处的强相关性。根据剩余的模式判断 p,q, P,Q。
模型摘要: 查看季节性参数 (ar.S, ma.S) 是否显著。比较 SARIMA 模型的 AIC/BIC 值与上一篇简单 ARIMA 模型的 AIC/BIC 值，SARIMA 的值通常会更低，表明模型拟合更好。
诊断图: 重点关注残差 ACF 图在季节性滞后位置是否还有显著尖峰。理想的 SARIMA 模型应该能消除这些季节性残留。检查 Ljung-Box p 值。
预测图: 观察 SARIMA 的预测是否能同时捕捉到趋势和季节性模式，这通常是它相对于简单 ARIMA 的主要优势。

(可选) 关于 SARIMAX 中的 ‘X’

statsmodels 实现的名字是 SARIMAX，那个 ‘X’ 代表外生变量 (eXogenous variables)。这意味着你可以在模型中加入其他可能影响目标时间序列的外部时间序列变量。

例如，在预测销售额时，你可以将广告支出、节假日标志等作为外生变量加入模型。

# 假设你有外生变量 exog_data (一个 Pandas DataFrame/Series，索引与主数据对齐)
# model = SARIMAX(co2_data, exog=exog_data, order=(p,d,q), seasonal_order=(P,D,Q,s))
# model_fit = model.fit()
# 预测时也需要提供未来的外生变量值:
# forecast = model_fit.get_forecast(steps=n_forecast, exog=future_exog_data)

注意： 使用外生变量进行预测时，你必须拥有（或能够准确预测）未来时期的外生变量值，这可能是一个挑战。

总结

本篇我们成功地将 ARIMA 模型升级，学会了使用 SARIMA(p,d,q)(P,D,Q)s 模型来处理具有明显季节性的时间序列数据：

理解了季节性参数 P, D, Q, s 的含义。
掌握了如何在模型识别阶段确定季节性差分 D 和季节性阶数 P, Q。
实践了使用 statsmodels.tsa.statespace.SARIMAX 进行完整的建模、诊断和预测流程。
认识到 SARIMA 在捕捉季节性模式上的优势，以及 SARIMAX 加入外生变量的能力。

SARIMA 是处理带有季节性趋势数据的强大经典工具。熟练掌握它，你将能解决更多样的时间序列预测问题。

下一篇预告：
虽然 SARIMA 很强大，但有时参数选择和模型诊断仍然繁琐。有没有更“自动化”、更易于使用的现代预测工具，尤其擅长处理复杂的季节性和节假日效应？下一篇，我们将探索由 Facebook 开源的流行预测库——Prophet！

敬请期待！

（尝试对 CO2 数据使用不同的 SARIMA 阶数组合，哪个组合的 AIC/BIC 最低且能通过诊断？改变季节性周期 s 会有什么影响？欢迎分享你的实验结果！）