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_nameNH.TsdSST, header1, index_col0) 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:]) # 排除第一个NaN3. 模型定阶与参数选择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, lags36, axax1) plot_pacf(diff_12, lags36, axax2) 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()诊断图包括标准化残差时序图 - 应无明显趋势残差直方图 - 应接近正态分布正态Q-Q图 - 点应大致在直线上残差自相关图 - 应无显著自相关4.2 白噪声检验使用Ljung-Box检验验证残差是否为白噪声from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(results.resid, lags12) print(p值:, lb_test[1])如果所有p值都大于0.05说明残差是白噪声模型拟合良好。5. 预测与评估5.1 时间序列预测使用训练好的模型进行预测# 预测未来12个月 forecast results.get_forecast(steps12) 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预测值, colorred) plt.fill_between(range(len(data), len(data)12), conf_int.iloc[:,0], conf_int.iloc[:,1], colorpink, alpha0.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(fMSE: {mse:.4f}) print(fMAE: {mae:.4f}) print(fR²: {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_col0) 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(steps12) 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预测值, colorred) plt.fill_between(range(len(ts_data), len(ts_data)12), conf_int.iloc[:,0], conf_int.iloc[:,1], colorpink, alpha0.3) plt.legend() plt.show()在实现过程中有几个常见坑需要注意数据量不足会导致模型波动大建议至少使用3-4个完整周期的数据差分过度会使数据失去原有特征ADF检验可以帮助判断参数搜索范围不宜过大否则计算成本会急剧上升残差检验不通过时需要回到数据预处理阶段重新检查