SARIMA模型实战：Python实现季节性时间序列预测

1. 项目缘起为什么季节性时序预测绕不开SARIMA做数据分析或者业务预测的朋友十有八九都跟时间序列打过交道。无论是预测下个月的销售额还是预估明天的网站流量时间序列模型都是工具箱里的常客。在众多模型中ARIMA自回归积分滑动平均模型因其强大的理论基础和广泛的适用性几乎成了入门必学。但很多新手包括当年的我在兴冲冲地用ARIMA跑完数据后经常会遇到一个尴尬的局面模型在训练集上表现还行一到预测未来尤其是预测带有明显周期性波动的数据时结果就惨不忍睹。比如预测电商的“双十一”销量或者预测夏季的空调耗电量ARIMA模型给出的预测线往往是一条平滑的直线完全忽略了那些周而复始的峰谷。这就是ARIMA模型的局限性它擅长捕捉趋势和随机性但对季节性Seasonality这种固定周期的重复模式显得力不从心。而现实世界的数据尤其是商业、气象、交通等领域的数据季节性几乎是标配。于是SARIMA季节性自回归积分滑动平均模型应运而生。你可以把它理解为ARIMA的“威力加强版”专门为处理季节性时间序列而生。今天我就结合自己踩过的坑和实战经验带你彻底搞懂SARIMA模型并用Python手把手实现一个完整的预测流程。无论你是想快速上手解决一个具体的业务预测问题还是想夯实时间序列分析的基础这篇文章都能给你提供一条清晰的路径。2. SARIMA模型核心原理不只是ARIMA的简单叠加在深入代码之前我们必须先理解SARIMA到底在做什么。很多教程一上来就扔公式SARIMA(p, d, q)(P, D, Q)s让人一头雾水。我们换个方式把它拆解成几个核心部分来理解。2.1 基石ARIMA模型的三驾马车SARIMA建立在ARIMA之上所以我们先快速回顾ARIMA的三大组件AR (p) - 自回归部分当前时刻的值是过去p个时刻值的线性组合。简单说就是“历史会重演”。p表示用过去多少期的数据来预测现在。比如p3就是用前3天的数据来预测第4天。I (d) - 差分部分这是为了让时间序列变得“平稳”。很多数据有长期上升或下降的趋势非平稳直接建模会出问题。d表示需要进行几次差分运算来消除趋势。一次差分就是用今天的数据减去昨天的数据得到的是“变化量”这个变化量序列往往更平稳。MA (q) - 移动平均部分当前时刻的值是过去q个时刻的预测误差白噪声的线性组合。简单说就是模型会“记住”自己过去犯了多大的错并用来修正现在的预测。q表示考虑过去多少期的预测误差。ARIMA(p,d,q) 就是把这三部分组合起来对非季节性部分进行建模。2.2 核心扩展引入季节性分量 (P, D, Q)sSARIMA 在 ARIMA 的基础上额外增加了一套完全相同的机制专门用来刻画季节性模式。这就是括号里的(P, D, Q)s。季节性周期 s这是最关键的一个参数。它定义了季节性的长度。比如月度数据有明显的年度规律那么s12。日度数据有明显的周规律那么s7。小时级数据有明显的日规律那么s24。季度数据那么s4。确定s是应用SARIMA的第一步通常基于业务常识或数据可视化如周期图来判断。季节性AR(P)不是用相邻时刻的数据而是用间隔一个季节性周期的历史数据来做回归。例如对于月度数据 (s12)P1意味着用去年同月的值来帮助预测今年这个月的值。季节性差分(D)为了消除季节性趋势比如每年销售额都在增长。操作是当前值 - 去年同期的值。D表示进行几次这样的季节性差分。通常D1就够了。季节性MA(Q)与季节性AR类似它考虑的是间隔一个季节性周期的历史预测误差。所以一个SARIMA(1,1,1)(1,1,1,12)模型可以解读为对非季节性部分用一阶差分 (d1) 消除趋势并用1阶自回归 (p1) 和1阶移动平均 (q1) 建模。对季节性部分周期为12用一阶季节性差分 (D1) 消除季节性趋势并用1阶季节性自回归 (P1) 和1阶季节性移动平均 (Q1) 建模。注意季节性差分 (D) 和普通差分 (d) 的执行顺序很重要。通常先做季节性差分再做普通差分。因为季节性波动可能掩盖长期趋势先消除季节性能让长期趋势更明显。2.3 模型定阶如何确定 (p,d,q)(P,D,Q)这是SARIMA建模中最具挑战性的一步。我们无法预先知道这7个参数的最佳值。传统统计学方法依赖ACF自相关函数和PACF偏自相关函数图来人工判断但这需要大量经验且对于季节性数据ACF/PACF图会在季节性周期倍数处出现截尾或拖尾判断起来更加复杂。现代更通用的做法是“网格搜索”我们设定一个参数范围让计算机自动尝试所有可能的参数组合然后选择一个评价指标最优的模型。最常用的评价指标是AIC赤池信息准则或BIC贝叶斯信息准则。它们的核心思想是权衡模型的拟合优度与复杂度值越小越好意味着模型用更少的参数达到了更好的拟合效果。虽然网格搜索计算量大但在Python的pmdarima库后文会介绍的帮助下这个过程可以自动化。我们只需要理解其原理算法会遍历我们设定的参数空间拟合每一个SARIMA模型计算其AIC最后选出AIC最小的那一组(p,d,q)(P,D,Q)作为最优模型。3. 环境准备与数据实战从0到1构建SARIMA预测模型理论说得再多不如一行代码。我们用一个经典的、自带强季节性的数据集来演示——航空乘客数据。这个数据集记录了1949年至1960年每月的航空乘客数量有明显的增长趋势和年度季节性。3.1 环境搭建与库导入首先确保你的Python环境安装了必要的库。我强烈建议使用Anaconda创建独立环境避免包冲突。# 在终端或Anaconda Prompt中创建环境可选 conda create -n timeseries python3.9 conda activate timeseries # 安装核心库 pip install pandas numpy matplotlib statsmodels pmdarima scikit-learn关键库说明pandas,numpy: 数据处理基石。matplotlib: 数据可视化。statsmodels: 包含完整的统计模型包括SARIMA的实现 (SARIMAX)。pmdarima: 神器封装了自动定阶ARIMA/SARIMA的功能能极大简化模型选择流程。scikit-learn: 用于后续的模型评估如计算误差指标。接下来在Jupyter Notebook或Python脚本中导入库并加载数据。import pandas as pd import numpy as np import matplotlib.pyplot as plt from statsmodels.tsa.statespace.sarimax import SARIMAX from statsmodels.tsa.seasonal import seasonal_decompose from pmdarima import auto_arima from sklearn.metrics import mean_absolute_error, mean_squared_error import warnings warnings.filterwarnings(ignore) # 忽略一些不影响运行的警告 # 加载航空乘客数据 url https://raw.githubusercontent.com/jbrownlee/Datasets/master/airline-passengers.csv df pd.read_csv(url) df[Month] pd.to_datetime(df[Month]) df.set_index(Month, inplaceTrue) series df[Passengers] print(series.head()) print(f\n数据形状: {series.shape}) print(f时间范围: {series.index.min()} 到 {series.index.max()})3.2 数据探索与可视化看见趋势和季节建模前必须“看清”你的数据。可视化是最直观的方式。plt.figure(figsize(12, 6)) plt.plot(series) plt.title(Monthly Airline Passengers (1949-1960)) plt.xlabel(Year) plt.ylabel(Passengers (Thousands)) plt.grid(True) plt.show()运行后你会看到一条明显向上增长且每年冬季波谷、夏季波峰的曲线。这直观地告诉我们数据有上升趋势和周期为12个月年度的季节性。为了更精确地分解我们可以使用seasonal_decompose函数。# 加法模型分解 (如果趋势和季节性幅度随时间基本不变) result_add seasonal_decompose(series, modeladditive, period12) # 乘法模型分解 (如果季节性幅度随趋势增长而增大本例更适用) result_mul seasonal_decompose(series, modelmultiplicative, period12) fig, axes plt.subplots(4, 2, figsize(15, 10)) result_add.observed.plot(axaxes[0, 0], titleObserved (Additive)) result_add.trend.plot(axaxes[1, 0], titleTrend (Additive)) result_add.seasonal.plot(axaxes[2, 0], titleSeasonal (Additive)) result_add.resid.plot(axaxes[3, 0], titleResidual (Additive)) result_mul.observed.plot(axaxes[0, 1], titleObserved (Multiplicative)) result_mul.trend.plot(axaxes[1, 1], titleTrend (Multiplicative)) result_mul.seasonal.plot(axaxes[2, 1], titleSeasonal (Multiplicative)) result_mul.resid.plot(axaxes[3, 1], titleResidual (Multiplicative)) plt.tight_layout() plt.show()观察分解图特别是残差图。乘法模型的残差看起来更随机、更平稳围绕0波动而加法模型的残差在后期方差似乎变大了。这印证了我们的观察航空乘客数量的季节性波动幅度是随着整体趋势增长而放大的。因此在本例中我们后续应该考虑使用乘法形式的SARIMA在statsmodels中可以通过对数据取对数来近似实现乘法模型。3.3 关键一步平稳性检验与数据变换SARIMA模型要求序列是“平稳”的即均值和方差不随时间变化。我们的原始数据显然不平稳有趋势。我们通过两种方式处理对数变换压缩数据尺度稳定方差同时将乘法关系转化为加法关系。np.log(series)差分消除趋势。包括普通差分和季节性差分。我们可以使用Augmented Dickey-Fuller (ADF) 检验来定量判断序列是否平稳。原假设是“序列非平稳”。p值小于显著性水平如0.05时我们拒绝原假设认为序列平稳。from statsmodels.tsa.stattools import adfuller def adf_test(timeseries): print(Results of Dickey-Fuller Test:) dftest adfuller(timeseries, autolagAIC) dfoutput pd.Series(dftest[0:4], index[Test Statistic, p-value, #Lags Used, Number of Observations Used]) for key, value in dftest[4].items(): dfoutput[fCritical Value ({key})] value print(dfoutput) if dfoutput[p-value] 0.05: print(\n结论序列是平稳的 (拒绝原假设)) else: print(\n结论序列是非平稳的 (无法拒绝原假设)) print(原始序列的ADF检验:) adf_test(series) print(\n *50 \n) print(对数变换后序列的ADF检验:) adf_test(np.log(series))通常对数变换后的序列仍然不平稳p值0.05因为它还有趋势。我们需要进一步差分。3.4 自动化模型定阶使用pmdarima解放双手手动看ACF/PACF图定阶对于季节性数据非常困难。这里我们祭出大杀器pmdarima.auto_arima。它会自动进行差分包括季节性差分确定d和D并搜索最佳的p, q, P, Q参数。# 对数据取对数以适应潜在的多重季节性乘法模型 log_series np.log(series) # 使用 auto_arima 自动寻找最佳 SARIMA 参数 # 设置 seasonalTrue, m12 表示我们处理的是周期为12的季节性数据 # traceTrue 会打印搜索过程 # suppress_warningsTrue 忽略一些提示性警告 # stepwiseTrue 使用逐步搜索更快False则进行更彻底的网格搜索更慢 auto_model auto_arima(log_series, start_p0, start_q0, max_p3, max_q3, # 非季节性ARMA最大阶数 start_P0, start_Q0, max_P2, max_Q2, # 季节性ARMA最大阶数通常不需要太高 m12, # 季节性周期 seasonalTrue, dNone, # 自动检测最优d DNone, # 自动检测最优D traceTrue, error_actionignore, suppress_warningsTrue, stepwiseTrue, information_criterionaic) # 使用AIC准则 print(auto_model.summary())运行后auto_arima会输出它找到的最佳模型参数例如可能是SARIMAX(0,1,1)(0,1,1,12)。同时summary()会给出模型的详细统计结果包括系数显著性、AIC/BIC值、残差诊断等。请务必关注系数coef的P值P|z|如果P值很大比如0.05说明该系数不显著对应的参数如某个AR或MA项可能没必要包含在模型中。不过auto_arima通常能给出一个不错的起点。3.5 手动建模、训练与诊断拿到auto_arima推荐的参数后我们可以用statsmodels的SARIMAX函数手动拟合模型以便获得更细致的控制和分析。# 假设 auto_arima 给出的最佳模型是 (0,1,1)(0,1,1,12) order (0, 1, 1) # (p, d, q) seasonal_order (0, 1, 1, 12) # (P, D, Q, s) # 使用对数序列进行拟合 model SARIMAX(log_series, orderorder, seasonal_orderseasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) model_fit model.fit(dispFalse) # dispFalse 不显示迭代信息 print(model_fit.summary())查看summary除了参数更重要的是进行残差诊断。一个好的模型其残差应该类似于白噪声均值为0方差恒定无自相关。# 绘制残差诊断图 model_fit.plot_diagnostics(figsize(12, 8)) plt.tight_layout() plt.show()诊断图包含四个子图标准化残差图残差应该随机分布在0附近没有明显模式。残差直方图 KDE密度估计最好与正态分布曲线红色重合表示残差近似正态分布。正态Q-Q图点应大致分布在红色直线上表示残差符合正态分布。残差自相关图ACF各阶滞后的自相关系数应迅速落入蓝色阴影区域置信区间表示无显著自相关。如果残差ACF图在某些滞后阶数特别是季节性滞后如1224上仍有显著 spikes说明模型未能完全捕捉季节性可能需要调整P或Q。3.6 进行预测并评估模型诊断通过后我们就可以进行预测了。我们需要将预测值从对数尺度转换回原始尺度。# 设定预测期数例如预测未来24个月2年 forecast_steps 24 # 使用 get_forecast 方法可以得到预测结果和置信区间 forecast_result model_fit.get_forecast(stepsforecast_steps) # 获取预测均值在对数尺度下 forecast_log forecast_result.predicted_mean # 获取置信区间在对数尺度下 confidence_interval_log forecast_result.conf_int() # 将对数尺度的预测值和置信区间转换回原始尺度 forecast np.exp(forecast_log) confidence_interval np.exp(confidence_interval_log) # 创建未来日期索引 last_date series.index[-1] future_dates pd.date_range(startlast_date pd.DateOffset(months1), periodsforecast_steps, freqMS) # 将预测结果组装成DataFrame forecast_series pd.Series(forecast, indexfuture_dates) lower_series pd.Series(confidence_interval.iloc[:, 0], indexfuture_dates) upper_series pd.Series(confidence_interval.iloc[:, 1], indexfuture_dates) # 绘制历史数据与预测结果 plt.figure(figsize(14, 7)) plt.plot(series, labelHistorical Data) plt.plot(forecast_series, labelForecast, colorred) plt.fill_between(future_dates, lower_series, upper_series, colorpink, alpha0.3, label95% Confidence Interval) plt.title(SARIMA Forecast for Airline Passengers) plt.xlabel(Year) plt.ylabel(Passengers) plt.legend() plt.grid(True) plt.show()为了评估模型在历史数据上的拟合效果我们可以计算一些常见的误差指标如MAE平均绝对误差、RMSE均方根误差和MAPE平均绝对百分比误差。MAPE尤其具有解释性它表示平均预测误差的百分比。# 获取模型对历史数据的拟合值in-sample forecast fitted_values np.exp(model_fit.fittedvalues) # 注意转换回原始尺度 # 由于差分前面部分数据没有拟合值需要对齐 # 对于 (0,1,1)(0,1,1,12) 模型由于进行了一阶普通差分和一阶季节性差分会丢失前13个数据点 start_idx 13 # 具体数值取决于 (dD*s) aligned_actual series.iloc[start_idx:] aligned_fitted fitted_values.iloc[start_idx:] # 计算误差指标 mae mean_absolute_error(aligned_actual, aligned_fitted) rmse np.sqrt(mean_squared_error(aligned_actual, aligned_fitted)) mape np.mean(np.abs((aligned_actual - aligned_fitted) / aligned_actual)) * 100 print(f模型在历史数据上的表现) print(fMAE: {mae:.2f}) print(fRMSE: {rmse:.2f}) print(fMAPE: {mape:.2f}%)一个MAPE小于10%的模型通常被认为是不错的小于5%则非常优秀。但这高度依赖于具体行业和数据波动性。4. 实战避坑指南与高阶技巧纸上得来终觉浅绝知此事要躬行。在实际项目中应用SARIMA你会遇到比教科书例子复杂得多的情况。下面分享几个我踩过坑后总结的关键点。4.1 如何判断该用加法模型还是乘法模型这是建模的第一步选错了模型形式效果会大打折扣。加法模型适用于季节性波动的幅度不随时间序列水平趋势变化的情况。即趋势和季节性相互独立地叠加。在图上表现为季节性的“波峰波谷”的宽度振幅在整个时间范围内大致恒定。乘法模型适用于季节性波动的幅度随着时间序列水平的上升/下降而同比增大的情况。在图上表现为序列值越大季节性波动就越剧烈。本文的航空乘客数据就是典型例子。判断方法可视化观察画出时间序列图。如果随着趋势上升季节性的“锯齿”明显变大倾向于乘法模型。分解观察使用seasonal_decompose分别用加法和乘法模型分解。观察残差图。残差应该看起来是随机的、没有模式的、方差大致恒定的。哪个模型的残差图更符合这个特征就选哪个。通常乘法模型的残差更稳定。业务理解很多经济、销售数据都是百分比增长的其季节性波动也往往是比例性的更适合乘法模型。在SARIMA中的实现statsmodels的SARIMAX默认是加法模型。要实现乘法模型通常的做法是先对原始数据取自然对数将乘法关系转化为加法关系因为 log(a*b) log(a) log(b)。然后对取对数后的序列拟合加法SARIMA模型最后将预测结果用指数函数np.exp()转换回来。我们上面的示例正是采用了这种方法。4.2 处理更复杂的季节性多重季节性与傅里叶项现实中的数据可能包含多个季节性周期。例如小时级的电力负荷数据同时具有日季节性(s124)每天用电高峰在傍晚。周季节性(s2168)工作日和周末的用电模式不同。年季节性(s38760)夏季和冬季用电量不同。标准的SARIMA(P,D,Q,s)只能处理一个季节性周期。对于多重季节性有几种应对策略SARIMA结合外部回归量将其他季节性周期通过傅里叶级数正弦余弦项或虚拟变量如“是否为周末”作为外生变量加入模型。SARIMAX中的X就支持外生变量。使用更高级的模型如 ProphetFacebook开源或 TBATS它们原生支持多重季节性。数据重采样如果主要关心某一个季节性可以重采样数据。例如将小时数据按天聚合来重点分析周季节性。当季节性周期非常长如s365对于日数据时标准的季节性(P,D,Q)参数会使得模型非常庞大难以估计。此时使用傅里叶项作为外生变量来近似季节性是一种常用且高效的技巧。pmdarima也支持通过seasonal参数设置复杂的季节性结构。4.3 模型评估与调优不要迷信AICauto_arima基于AIC选模是一个很好的起点但绝不是终点。你必须进行样本外预测Out-of-Sample Forecast验证。正确做法滚动预测验证将数据按时间顺序分为训练集和测试集例如用前80%的数据训练预测后20%。在训练集上拟合模型。预测测试集的第一步记录误差。将真实的测试集第一个值纳入训练集重新拟合模型或更新模型状态预测下一步。重复步骤3-4直到遍历整个测试集。这个过程模拟了真实世界中利用最新数据不断更新模型进行预测的场景。计算在整个测试集上的平均误差如MAPE。这个指标比模型在训练集上的拟合优度如AIC更能反映模型的真实预测能力。# 滚动预测验证示例框架 def rolling_forecast_validation(data, train_size, order, seasonal_order): train, test data[:train_size], data[train_size:] history list(train) predictions [] for t in range(len(test)): model SARIMAX(history, orderorder, seasonal_orderseasonal_order) model_fit model.fit(dispFalse) yhat model_fit.forecast()[0] # 预测下一步 predictions.append(yhat) history.append(test[t]) # 将真实值加入历史模拟实时更新 # 注意实际中频繁重新拟合整个模型计算成本高可考虑使用 update 方法 # 计算误差 mape np.mean(np.abs((np.array(test) - np.array(predictions)) / np.array(test))) * 100 return predictions, mape # 注意此示例使用全部历史数据重新拟合效率低。对于SARIMAX可以使用 model_fit.append 或 model_fit.apply 来更新模型效率更高。如果滚动预测的误差很大可能需要回到auto_arima阶段调整搜索范围如增大max_p,max_q或者尝试不同的变换如 Box-Cox 变换。4.4 当预测结果不理想时你的检查清单数据真的平稳了吗再次检查ADF检验结果。确保经过正确的差分d和D后序列是平稳的。可以绘制差分后的序列图观察。季节性周期s确定对了吗这是根本。用业务知识结合周期图、季节性分解图反复确认。残差是白噪声吗仔细查看诊断图中的ACF图。如果在滞后1, 12, 24等处仍有显著的自相关说明模型没有捕捉到某些模式需要增加p, q, P, Q的阶数。有没有异常值或结构性突变比如促销活动、政策变化导致的数据骤变。SARIMA无法处理这种“断层”可能需要先检测并处理异常值或使用包含干预分析Intervention Analysis的模型。趋势是线性的吗SARIMA通过差分来处理趋势这隐含地假设趋势是随机游走或多项式形式。如果趋势是指数型的对数变换可能先于差分。如果趋势复杂可能需要先进行去趋势化处理或使用其他模型。参数是否显著检查model_fit.summary()中系数的P值。如果某些项的P值很大0.1可以考虑移除该参数重新拟合一个更简洁的模型。SARIMA是一个强大但复杂的模型。它假设数据生成过程是线性的、参数是固定的。对于非线性、波动聚集如金融数据或具有复杂外部影响的数据可能需要考虑神经网络如LSTM、Prophet或GARCH等模型。没有万能的模型只有最适合具体场景的模型。理解SARIMA的原理和局限能让你在时间序列预测的武器库中又多了一件得心应手的利器。