SARIMA模型实战：从数据预处理到预测评估的完整Python实现

1. 时间序列分析与SARIMA模型基础

时间序列分析是数据科学领域的重要分支，它专门研究按时间顺序排列的数据点。想象你每天记录体重变化，这些数据点按日期排列就构成了一个典型的时间序列。在实际应用中，时间序列分析可以帮助我们预测股票价格、销售额变化、气象数据等具有时间依赖性的数据。

SARIMA（季节性自回归综合移动平均）模型是ARIMA模型的升级版，专门处理具有季节性特征的数据。比如空调销量每年夏天都会上升，这种周期性变化就是典型的季节性特征。SARIMA模型通过(p,d,q)×(P,D,Q,s)这组参数来描述数据的特征：

• p：自回归项阶数
• d：差分次数
• q：移动平均项阶数
• P：季节性自回归阶数
• D：季节性差分次数
• Q：季节性移动平均阶数
• s：季节周期长度

与传统ARIMA相比，SARIMA增加了处理季节性变化的能力。举个例子，如果分析月度销售数据，s通常设为12（一年12个月）；如果是季度数据，s则设为4。

2. 数据准备与预处理实战

2.1 数据获取与探索

我们从NASA戈达德空间研究所获取了北美陆地表面温度数据，时间跨度为1990年1月至2020年11月。原始数据格式如下：

import pandas as pd data = pd.read_excel('temperature.xlsx', sheet_name='NH.Ts+dSST', header=1, index_col=0) print(data.head())

数据探索是建模的第一步。我们先绘制原始数据的时序图：

import matplotlib.pyplot as plt plt.figure(figsize=(12,6)) plt.plot(data.values.ravel()[:-1]) # 排除最后一个空值 plt.title('北美陆地温度变化(1990-2020)') plt.xlabel('月份') plt.ylabel('温度(℃)') plt.show()

2.2 平稳化处理技巧

原始数据通常包含趋势和季节性成分。我们采用一阶12步差分来消除这些影响：

# 一阶差分消除趋势 diff_1 = data.diff(1).dropna() # 12步差分消除季节性 diff_12 = diff_1.diff(12).dropna()

平稳性检验使用ADF单位根检验：

from statsmodels.tsa.stattools import adfuller result = adfuller(diff_12) print('ADF统计量:', result[0]) print('p值:', result[1]) print('临界值:', result[4])

如果p值小于0.05，可以认为数据已平稳。对于不平稳的数据，可以尝试Box-Cox变换：

from scipy.stats import boxcox transformed, _ = boxcox(data.values.flatten()[1:]) # 排除第一个NaN

3. 模型定阶与参数选择

3.1 自相关与偏自相关分析

ACF和PACF图是定阶的重要工具：

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf fig, (ax1, ax2) = plt.subplots(2,1, figsize=(12,8)) plot_acf(diff_12, lags=36, ax=ax1) plot_pacf(diff_12, lags=36, ax=ax2) plt.show()

通过观察截尾和拖尾特征，可以初步判断p、q值。但这种方法比较主观，我们需要更客观的准则。

3.2 自动化参数搜索

使用AIC/BIC准则进行参数网格搜索：

from itertools import product from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings("ignore") # 参数范围设置 ps = range(0, 3) # AR阶数 ds = range(0, 2) # 差分次数 qs = range(0, 3) # MA阶数 Ps = range(0, 2) # 季节性AR阶数 Ds = range(0, 2) # 季节性差分 Qs = range(0, 2) # 季节性MA阶数 s = 12 # 季节周期 # 生成所有参数组合 params_list = list(product(ps, ds, qs, Ps, Ds, Qs)) def find_best_params(data, params_list): best_bic = float('inf') best_params = None for param in params_list: try: model = SARIMAX(data, order=(param[0], param[1], param[2]), seasonal_order=(param[3], param[4], param[5], s)) results = model.fit(disp=-1) if results.bic < best_bic: best_bic = results.bic best_params = param except: continue return best_params, best_bic best_params, best_bic = find_best_params(data, params_list) print(f'最优参数: SARIMA{best_params[:3]}×{best_params[3:]}{s}') print(f'最小BIC值: {best_bic}')

4. 模型训练与诊断

4.1 模型拟合与残差分析

使用最优参数训练模型：

model = SARIMAX(data, order=(best_params[0], best_params[1], best_params[2]), seasonal_order=(best_params[3], best_params[4], best_params[5], s)) results = model.fit(disp=-1) # 残差诊断 results.plot_diagnostics(figsize=(15,12)) plt.show()

诊断图包括：

1. 标准化残差时序图 - 应无明显趋势
2. 残差直方图 - 应接近正态分布
3. 正态Q-Q图 - 点应大致在直线上
4. 残差自相关图 - 应无显著自相关

4.2 白噪声检验

使用Ljung-Box检验验证残差是否为白噪声：

from statsmodels.stats.diagnostic import acorr_ljungbox lb_test = acorr_ljungbox(results.resid, lags=12) print('p值:', lb_test[1])

如果所有p值都大于0.05，说明残差是白噪声，模型拟合良好。

5. 预测与评估

5.1 时间序列预测

使用训练好的模型进行预测：

# 预测未来12个月 forecast = results.get_forecast(steps=12) forecast_mean = forecast.predicted_mean conf_int = forecast.conf_int() # 绘制预测结果 plt.figure(figsize=(12,6)) plt.plot(data.values.ravel()[:-1], label='观测值') plt.plot(range(len(data), len(data)+12), forecast_mean, label='预测值', color='red') plt.fill_between(range(len(data), len(data)+12), conf_int.iloc[:,0], conf_int.iloc[:,1], color='pink', alpha=0.3) plt.legend() plt.show()

5.2 模型评估指标

使用MSE、MAE和R²评估预测效果：

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score # 假设test_data是真实值 mse = mean_squared_error(test_data, forecast_mean) mae = mean_absolute_error(test_data, forecast_mean) r2 = r2_score(test_data, forecast_mean) print(f'MSE: {mse:.4f}') print(f'MAE: {mae:.4f}') print(f'R²: {r2:.4f}')

在实际项目中，我遇到过MSE很低但预测曲线明显偏离的情况。这时候需要结合业务知识判断，不能完全依赖指标。比如温度预测中，0.5℃的误差对日常穿衣影响不大，但对农作物生长可能很关键。

6. 完整代码实现

以下是整合后的完整代码示例：

# 导入所需库 import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.tsa.statespace.sarimax import SARIMAX from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from sklearn.metrics import mean_squared_error from itertools import product # 1. 数据加载与预处理 data = pd.read_excel('temperature.xlsx', index_col=0) ts_data = data.values.ravel()[:-1] # 2. 平稳性检验与处理 diff_1 = pd.Series(ts_data).diff(1).dropna() diff_12 = diff_1.diff(12).dropna() # 3. 参数搜索 def find_best_params(data): ps = range(0, 3) ds = range(0, 2) qs = range(0, 3) Ps = range(0, 2) Ds = range(0, 2) Qs = range(0, 2) s = 12 params_list = list(product(ps, ds, qs, Ps, Ds, Qs)) best_bic = np.inf best_params = None for param in params_list: try: model = SARIMAX(data, order=(param[0], param[1], param[2]), seasonal_order=(param[3], param[4], param[5], s)) results = model.fit(disp=-1) if results.bic < best_bic: best_bic = results.bic best_params = param except: continue return best_params, best_bic best_params, best_bic = find_best_params(ts_data) # 4. 模型训练 final_model = SARIMAX(ts_data, order=(best_params[0], best_params[1], best_params[2]), seasonal_order=(best_params[3], best_params[4], best_params[5], 12)) final_results = final_model.fit(disp=-1) # 5. 预测与评估 forecast = final_results.get_forecast(steps=12) forecast_mean = forecast.predicted_mean conf_int = forecast.conf_int() # 可视化 plt.figure(figsize=(12,6)) plt.plot(ts_data, label='观测值') plt.plot(range(len(ts_data), len(ts_data)+12), forecast_mean, label='预测值', color='red') plt.fill_between(range(len(ts_data), len(ts_data)+12), conf_int.iloc[:,0], conf_int.iloc[:,1], color='pink', alpha=0.3) plt.legend() plt.show()

在实现过程中，有几个常见坑需要注意：

1. 数据量不足会导致模型波动大，建议至少使用3-4个完整周期的数据
2. 差分过度会使数据失去原有特征，ADF检验可以帮助判断
3. 参数搜索范围不宜过大，否则计算成本会急剧上升
4. 残差检验不通过时，需要回到数据预处理阶段重新检查