用Python实现时间序列模型实战——Day 8: 季节性ARIMA模型 (SARIMA)

一、学习内容

1. SARIMA 模型的定义与公式推导

SARIMA 模型：

SARIMA 模型是扩展了 ARIMA 模型的一种方法，全称为季节性自回归积分滑动平均模型（Seasonal AutoRegressive Integrated Moving Average）。它结合了 ARIMA 模型的非季节性部分和季节性成分，用于处理具有季节性模式的时间序列数据。
SARIMA 模型可以表示为

$\Phi_P(B^m) \phi_p(B) \nabla^d \nabla_m^D y_t = \mu + \Theta_Q(B^m) \theta_q(B) \epsilon_t$

其中：

$y_t$ 是时间 $t$ 的观察值。
$\phi_p(B)$ 是非季节性自回归多项式，阶数为 $p$ 。
$\Phi_P(B^m)$ 是季节性自回归多项式，阶数为 $P$ 。
$\theta_q(B)$ 是非季节性移动平均多项式，阶数为 $q$ 。
$\Theta_Q(B^m)$ 是季节性移动平均多项式，阶数为 $Q$ 。
$\nabla^d$ 表示非季节性差分操作，次数为 $d$ 。
$\nabla_m^D$ 表示季节性差分操作，次数为 $D$ 。
$\epsilon_t$ 是白噪声项。
$m$ 是季节性周期。

SARIMA 模型的公式：

SARIMA 模型的数学表达式为：

$(1 - \Phi_1 B^{12})(1 - \phi_1 B)(1 - B)(1 - B^{12})y_t = \mu + (1 + \Theta_1 B^{12})(1 + \theta_1 B) \epsilon_t$

其中：

$\phi_1$ 是非季节性自回归项的系数。
$\Phi_1$ 是季节性自回归项的系数。
$\theta_1$ 是非季节性移动平均项的系数。
$\Theta_1$ 是季节性移动平均项的系数。
$B$ 是滞后算子，表示滞后 1 期。
$\epsilon_t$ 是白噪声项。

2. SARIMA 模型的参数识别与选择 (P, D, Q)

非季节性参数识别：

$p$ （非季节性自回归项）：通过偏自相关函数（PACF）图确定，如果 PACF 在某个滞后期后截尾，截尾点即为 $p$ 的值。
$d$ （非季节性差分次数）：通过差分使序列平稳，差分次数 $d$ 通常根据时间序列的平稳性检验或观察差分后的序列图来选择。
$q$ （非季节性移动平均项）：通过自相关函数（ACF）图确定，如果 ACF 在某个滞后期后截尾，截尾点即为 $q$ 的值。

季节性参数识别：

$P$ （季节性自回归项）：通过季节性 PACF 图确定。
$D$ （季节性差分次数）：通过季节性差分使序列平稳，季节性差分次数 $D$ 根据季节性周期和数据的季节性特征确定。
$Q$ （季节性移动平均项）：通过季节性 ACF 图确定。
$m$ （季节性周期）：季节性周期 $m$ 通常是已知的，如每年或每季度。

3. SARIMA 模型的季节性成分处理

季节性成分处理：

在 SARIMA 模型中，季节性成分通过季节性自回归（SAR）、季节性差分（SD）、季节性移动平均（SMA）项进行处理。这些季节性成分使模型能够捕捉周期性模式，并提升对季节性数据的预测能力。

二、实战案例

下面我们将使用 statsmodels 库对航空乘客数据进行 SARIMA 模型的拟合与预测。

1. 数据加载与原始数据可视化

import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf# 加载时间序列数据集
url = "https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv"
data = pd.read_csv(url, parse_dates=['Month'], index_col='Month')# 绘制原始数据的时间序列图
plt.figure(figsize=(12, 6))
plt.plot(data['Passengers'], label='Original')
plt.title('Original Time Series')
plt.xlabel('Date')
plt.ylabel('Number of Passengers')
plt.legend()
plt.grid(True)
plt.show()

程序解释：

载入航空乘客数据集，并绘制原始时间序列图，观察数据的趋势和季节性波动。

结果输出：

原始数据的时间序列图：图表显示了航空乘客数量随时间的变化，表现出明显的上升趋势和季节性波动。

2. 一阶差分与季节性差分

# 进行一阶差分以去除趋势
data_diff = data.diff().dropna()# 进行季节性差分以去除季节性成分
data_seasonal_diff = data_diff.diff(12).dropna()# 绘制差分后的时间序列图
plt.figure(figsize=(12, 6))
plt.plot(data_seasonal_diff['Passengers'], label='Seasonally Differenced')
plt.title('Seasonally Differenced Time Series')
plt.xlabel('Date')
plt.ylabel('Differenced Passengers')
plt.legend()
plt.grid(True)
plt.show()

程序解释：

对时间序列进行一阶差分，以去除趋势。
对差分后的数据进行季节性差分，以去除季节性成分。季节性差分是针对周期为 12 的数据进行差分（即 diff(12)），消除一年中的季节性波动。

结果输出：

差分后的时间序列图：一阶差分后的数据消除了趋势，而季节性差分进一步消除了季节性成分，表现为更加平稳的时间序列。

3. ACF 和 PACF 图的绘制

# 绘制 ACF 和 PACF 图
plt.figure(figsize=(12, 6))
plt.subplot(121)
plot_acf(data_seasonal_diff, lags=40, ax=plt.gca())
plt.title('ACF of Seasonally Differenced Data')plt.subplot(122)
plot_pacf(data_seasonal_diff, lags=40, ax=plt.gca())
plt.title('PACF of Seasonally Differenced Data')plt.tight_layout()
plt.show()

程序解释：

绘制季节性差分后的 ACF 和 PACF 图，帮助确定 SARIMA 模型的参数。

结果输出：

ACF 和 PACF 图：ACF 和 PACF 图展示了季节性差分后的数据自相关结构，通过这些图可以确定 SARIMA 模型的参数。

4. SARIMA 模型拟合

# SARIMA 模型拟合 (p, d, q) x (P, D, Q, m) = (1, 1, 1) x (1, 1, 1, 12)
model = SARIMAX(data['Passengers'], order=(1, 1, 1), seasonal_order=(1, 1, 1, 12))
results = model.fit()# 输出模型摘要
print(results.summary())

程序解释：

使用 SARIMAX 函数拟合 SARIMA 模型。本例中选择 $(p, d, q) = (1, 1, 1)$ 和 $(P, D, Q, m) = (1, 1, 1, 12)$ 来拟合模型。
输出模型的摘要信息，包括 AIC、BIC 等指标。

结果输出：

SARIMAX Results                                      
==========================================================================================
Dep. Variable:                         Passengers   No. Observations:                  144
Model:             SARIMAX(1, 1, 1)x(1, 1, 1, 12)   Log Likelihood                -506.149
Date:                            Mon, 02 Sep 2024   AIC                           1022.299
Time:                                    12:59:16   BIC                           1036.675
Sample:                                01-01-1949   HQIC                          1028.140- 12-01-1960                                         
Covariance Type:                              opg                                         
==============================================================================coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1272      0.356     -0.357      0.721      -0.825       0.570
ma.L1         -0.2149      0.325     -0.660      0.509      -0.853       0.423
ar.S.L12      -0.9272      0.214     -4.341      0.000      -1.346      -0.509
ma.S.L12       0.8395      0.309      2.717      0.007       0.234       1.445
sigma2       130.7819     15.420      8.481      0.000     100.559     161.005
===================================================================================
Ljung-Box (L1) (Q):                   0.00   Jarque-Bera (JB):                 7.05
Prob(Q):                              0.99   Prob(JB):                         0.03
Heteroskedasticity (H):               2.65   Skew:                             0.13
Prob(H) (two-sided):                  0.00   Kurtosis:                         4.11
===================================================================================

SARIMA 模型的摘要：模型摘要中包含了各参数的估计值、标准误差、t 统计量，以及 AIC/BIC 等信息准则，用于评估模型的拟合效果。

5. 模型预测

# 进行预测
pred = results.get_forecast(steps=24)
pred_ci = pred.conf_int()# 绘制预测结果
plt.figure(figsize=(12, 6))
plt.plot(data.index, data['Passengers'], label='Observed')
plt.plot(pred.predicted_mean.index, pred.predicted_mean, color='red', label='Forecast')
plt.fill_between(pred_ci.index, pred_ci.iloc[:, 0], pred_ci.iloc[:, 1], color='pink', alpha=0.3)
plt.title('SARIMA Model Forecast')
plt.xlabel('Date')
plt.ylabel('Number of Passengers')
plt.legend()
plt.grid(True)
plt.show()

程序解释：