ARIMA模型进阶——具有季节性的SARIMAX模型实现

SARIMAX模型是一种强大的时间序列分析和预测模型。它是对ARIMA模型的扩展，增加了对季节性成分和外部因素的考虑

1. SARIMAX模型的主要组成部分

• AR： SARIMAX模型中的自回归部分表示当前观测值与过去观测值的线性关系， 这反映了时间序列数据中存在的自相关性，即过去的值对当前值有影响
• I： 差分操作，用于使时间序列平稳化， 平稳化是为了消除数据中的趋势和季节性，使模型更容易捕捉数据的本质规律（差分操作主要用于使时间序列平稳，它并不一定能够完全消除季节性）
• MA： 移动平均部分表示当前观测值与过去白噪声（随机误差）的线性关系。 它帮助模型捕捉数据中的非自相关性
• 季节性成分： SARIMAX模型引入了季节性成分，允许模型捕捉数据中周期性的季节性变化。 这对于处理具有明显季节性趋势的时间序列数据非常重要
• eXogenous factors： 外部因素部分允许引入额外的外部变量，这些变量可能影响时间序列数据。 这使得模型更加灵活，可以考虑除时间序列自身之外的其他因素

2. 代码实列

2.1 数据展示

          
import pandas as pd
          
df = pd.read_excel('data.xlsx',index_col=0, parse_dates=['Month'])
          
df.info()
          
df.head()
      

该数据集为1949-1960某国际航空公司乘客数据

          
import matplotlib.pyplot as plt
          
plt.rcParams['font.sans-serif'] = 'SimHei' # 设置中文显示
          
plt.rcParams['axes.unicode_minus'] = False
          
# 绘制原始数据和滤波后的数据
          
plt.figure(figsize=(15, 5))
          
plt.plot(df['data'], label='原始数据', color='r',  alpha=0.3)
          
plt.title('时序图')
          
plt.grid(True)
          
plt.legend()
          
plt.show()
      

2.2 数据集分割

          
train_data = df['1949-01-01':'1959-12-01']
          
test_data = df['1960-01-01':'1960-12-01']
      

将数据集分割为训练集和测试集，其目的是评估模型的性能和泛化能力， 训练集包含从1949年1月到1959年12月的时间段的数据。这部分数据用于训练模型，即拟合 SARIMAX 模型的参数， 测试集 包含1960年1月到1960年12月的时间段的数据。这部分数据用于评估模型的性能，验证模型在未来时间点上的预测准确度

2.3 SARIMAX参数选择

          
from statsmodels.tsa.statespace.sarimax import SARIMAX
          
import warnings
          

          
# 定义迭代范围
          
max_p = 3  # p（自回归阶数）
          
max_d = 2  # d（差分阶数）数据以及平稳不需要差分
          
max_q = 1  # q（移动平均阶数）
          

          
max_P = 3  # 季节性自回归阶数
          
max_D = 1  # 季节性差分阶数
          
max_Q = 1  # 季节性移动平均阶数
          
seasonal_period = 12 # 季节周期为12 一年一周期12个月
          

          
best_aic = float('inf')  # 初始化最小 AIC 值
          
best_params = None  # 初始化最优参数组合
          

          
# 迭代搜索最优参数组合
          
for p in range(max_p + 1):
          
    for d in range(max_d + 1):
          
        for q in range(max_q + 1):
          
            for P in range(max_P + 1):
          
                for D in range(max_D + 1):
          
                    for Q in range(max_Q + 1):
          
                        if p + d + q > 0:
          
                            try:
          
                                # 使用 SARIMAX 模型拟合数据
          
                                with warnings.catch_warnings():
          
                                    warnings.filterwarnings("ignore", category=UserWarning)
          
                                    model = SARIMAX(train_data['data'], order=(p, d, q), 
          
                                                    seasonal_order=(P, D, Q, seasonal_period),
          
                                                   enforce_stationarity=False,# 是否强制模型保持平稳性
          
                                                   enforce_invertibility=False )# 可逆性的参数
          
                                    results = model.fit()
          

          
                                # 获取拟合的 AIC 值
          
                                aic = results.aic
          

          
                                # 更新最小 AIC 值和对应的参数组合
          
                                if aic < best_aic:
          
                                    best_aic = aic
          
                                    best_params = (p, d, q, P, D, Q)
          

          
                            except:
          
                                continue
          

          
# 打印最小 AIC 值和对应的参数组合
          
print("Best Parameters (p, d, q, P, D, Q):", best_params)
          
print("Best AIC:", best_aic)
      

代码通过迭代搜索不同的SARIMAX模型参数组合，从而找到在训练集上拟合效果最优的模型。具体而言，它循环遍历了指定范围内的自回归阶数（p）、差分阶数（d）、移动平均阶数（q）、季节性自回归阶数（P）、季节性差分阶数（D）、季节性移动平均阶数（Q）等参数。在每一组参数的组合下，使用 SARIMAX 模型拟合训练集数据，并计算模型的AIC（赤池信息准则）值。 最终，输出最小AIC值和对应的最优参数组合，这表示在给定的参数范围内，选择这组参数可以得到对训练集拟合效果最好的 SARIMAX 模型

2.4 模型训练

          
model = SARIMAX(train_data['data'], 
          
              order=(3, 2, 1),
          
              seasonal_order=(3,1,0,12),
          
              enforce_stationarity=False,
          
              enforce_invertibility=False)
          
results = model.fit()
          

          
results.plot_diagnostics(figsize=(20, 14))
          
plt.show()
      

在上面的图表中，可以观察到残差是不相关的（右下图），并且不表现出明显的季节性（左上图），此外，残差大致呈零均值的正态分布（右上图），左下角的 qq-plot 表明残差的有序分布（蓝点）大致遵循标准正态分布 N(0, 1) 的样本的线性趋势

在SARIMAX模型中，这些诊断结果的意义在于验证模型的合适性，确保模型在拟合数据时没有遗漏结构，而且残差是符合模型假设的，如果残差不相关、无季节性、零均值正态分布，模型通常被认为是有效的

2.5 模型预测

          
plt.figure(figsize=(15, 5))
          
plt.plot(train_data['data'], label='训练集', color='r',  alpha=0.7)
          
plt.plot(train_data.index[24::], results.fittedvalues.dropna().iloc[24::], label='模型在训练集上预测', color='k',  alpha=0.6)
          
plt.plot(test_data, label='测试集', color='y',  alpha=0.8)
          
plt.plot(results.forecast(steps=12), label='模型在测试集上预测', color='c',  alpha=0.4)
          
plt.plot(results.forecast(steps=36).iloc[12::], label='预测未来两年', color='b',  alpha=0.6)
          
plt.grid(True)
          
plt.legend()
          
plt.show()
      

从图中可以看出，该模型在时间序列建模方面做得很好，在训练集上灰色和红色的线非常接近，重点关注的测试集也存在高度吻合，蓝色的线为未来预测，认为这样的一个简单SARIMA模型能够进行预测未来数据

2.6 测试集 MAPE计算

          
# 获取预测值和实际值
          
forecast_values = results.forecast(steps=12)
          
actual_values = test_data['data']
          

          
# 计算 MAPE
          
def calculate_mape(actual, forecast):
          
    return np.mean(np.abs((actual - forecast) / actual)) * 100
          

          
mape = calculate_mape(actual_values, forecast_values)
          

          
print(f"MAPE: {mape:.2f}%")
      

在测试集上计算平均绝对百分比误差（MAPE）是评估时间序列预测模型性能的一种方法，MAPE度量的是模型的预测值与实际值之间的百分比误差的平均值，通过SARIMAX模型对时间序列数据的建模与预测，得到了一个在测试集上表现良好的模型，其平均绝对百分比误差（MAPE）为4.50%。这意味着模型对未来的预测相对准确，具备一定的可信度

结论

通过迭代搜索不同SARIMAX模型参数组合，从训练集中选择最优模型，然后在测试集上评估其性能，最终输出预测结果和模型评估指标，从而避免了ARIMA模型在处理具有季节性的时间序列数据时的局限性，提高了对季节性变化的更灵活和准确的建模能力，当然其中有些参数不需要进行迭代确认如差分阶数可根据ADF平稳性检验进行确定，参考以下链接时间序列——平稳性检验，更多详情不一一赘述，请参考往期文章

如果你对类似于这样的文章感兴趣。

欢迎关注、点赞、转发~