用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()

程序解释：

使用拟合的 SARIMA 模型进行未来 24 个月的乘客数量预测。
绘制预测结果图，并包含置信区间（confidence interval），展示模型预测的范围。

结果输出：

模型预测结果：预测结果图展示了未来 24 个月的航空乘客数量预测值，并在图中包含了置信区间，表示预测的不确定性范围。预测结果与观测数据的趋势和季节性波动相一致，表明 SARIMA 模型能够有效捕捉季节性时间序列的数据结构。

三、运行结果分析

1. 季节性处理

SARIMA 模型通过季节性自回归、季节性差分和季节性移动平均的组合，成功捕捉了时间序列中的季节性成分。这种处理方式特别适用于具有周期性模式的时间序列。

2. ACF 和 PACF 图分析

通过季节性差分后的 ACF 和 PACF 图，确定了 SARIMA 模型的参数。这些图表帮助识别了数据中的自相关性和季节性成分，为模型参数的选择提供了依据。

3. SARIMA 模型的拟合与预测

SARIMA 模型拟合后的结果表明，模型能够有效地解释和预测时间序列数据。AIC/BIC 指标可以用来评估模型的优劣。
预测结果显示了未来的乘客数量，并且预测结果与观测数据的季节性波动相一致，表明模型的预测效果良好。

通过这次学习，您掌握了如何构建和应用 SARIMA 模型来处理具有季节性模式的时间序列数据。SARIMA 模型在实际应用中非常有用，特别是在预测那些具有周期性模式的数据时，如销售数据、气象数据等。

四、补充更新——非季节性时间序列的案例分析

以下是一个非季节性时间序列的案例分析，它展示了如何使用 ARIMA 模型来分析和预测数据。这个案例使用的时间序列数据没有明显的季节性成分，主要体现趋势和随机性。

1. 案例：某公司月度销售额预测

假设我们有一家公司的月度销售额数据，该数据呈现出一定的趋势性，但没有显著的季节性波动。我们将使用 ARIMA 模型来分析和预测未来的销售额。

2. 数据生成与可视化

首先，我们生成一个模拟的非季节性时间序列数据，并进行可视化：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA

# 生成非季节性时间序列数据
np.random.seed(42)
time = np.arange(1, 101)
sales = 50 + 2 * time + np.random.normal(loc=0.0, scale=10, size=len(time))  # 线性趋势 + 随机噪声

# 创建数据框
data = pd.DataFrame({'Month': pd.date_range(start='2020-01', periods=len(time), freq='M'), 'Sales': sales})
data.set_index('Month', inplace=True)

# 绘制原始数据的时间序列图
plt.figure(figsize=(10, 6))
plt.plot(data, label='Monthly Sales')
plt.title('Company Monthly Sales Time Series')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.legend()
plt.grid(True)
plt.show()

数据说明：

在这个示例中，数据展示了一个月度销售额的时间序列。我们假设这个数据没有明显的季节性变化，但呈现出一定的上升趋势和随机波动。

结果输出：

3. ARIMA 模型的拟合

接下来，我们将对该时间序列数据进行差分处理，以消除趋势，确保数据平稳，然后使用 ARIMA 模型进行拟合。

# 差分以使数据平稳
data_diff = data.diff().dropna()

# 使用 ARIMA 模型进行拟合
model = ARIMA(data, order=(1, 1, 1))
results = model.fit()

# 输出模型摘要
print(results.summary())

结果输出：

 SARIMAX Results                                
==============================================================================
Dep. Variable:                  Sales   No. Observations:                  100
Model:                 ARIMA(1, 1, 1)   Log Likelihood                -379.852
Date:                Mon, 02 Sep 2024   AIC                            765.704
Time:                        21:36:34   BIC                            773.490
Sample:                    01-31-2020   HQIC                           768.854
                         - 04-30-2028                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1         -0.1250      0.199     -0.628      0.530      -0.515       0.265
ma.L1         -0.4721      0.192     -2.461      0.014      -0.848      -0.096
sigma2       125.4660     18.589      6.750      0.000      89.033     161.899
===================================================================================
Ljung-Box (L1) (Q):                   2.96   Jarque-Bera (JB):                 0.73
Prob(Q):                              0.09   Prob(JB):                         0.70
Heteroskedasticity (H):               1.13   Skew:                            -0.19
Prob(H) (two-sided):                  0.73   Kurtosis:                         3.17
===================================================================================

# 进行未来12个月的预测
forecast = results.get_forecast(steps=12)
forecast_ci = forecast.conf_int()

# 绘制预测结果
plt.figure(figsize=(10, 6))
plt.plot(data, label='Observed Sales')
plt.plot(forecast.predicted_mean.index, forecast.predicted_mean, color='red', label='Forecasted Sales')
plt.fill_between(forecast_ci.index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='red', alpha=0.3)
plt.title('ARIMA Model Sales Forecast')
plt.xlabel('Date')
plt.ylabel('Sales')
plt.legend()
plt.grid(True)
plt.show()

结果输出：