1. 项目概述从ARIMA到SARIMA的跃迁如果你做过时间序列预测大概率绕不开ARIMA模型。它经典、强大是很多数据分析师和建模初学者的首选。但当你处理的数据呈现出明显的周期性波动时——比如电商网站的日销售额在周末飙升电力负荷在夏季和冬季形成高峰或者某种商品的月度销量在每年春节前后规律性上涨——标准的ARIMA模型就会显得力不从心。它捕捉不到这种固定周期的重复模式预测的曲线会变得平滑丢失掉关键的季节性信息。这时你就需要请出它的升级版季节性自回归积分滑动平均模型也就是SARIMA。SARIMA模型可以看作是ARIMA模型在时间序列分析领域的“完全体”。它在ARIMA的三个核心参数自回归阶数p、差分阶数d、移动平均阶数q之外额外引入了一套专门处理季节性成分的参数P, D, Q, s。这套组合拳让模型具备了“双线作战”的能力一套参数p,d,q负责捕捉数据中非季节性的趋势和随机波动另一套参数P,D,Q,s则专门用来建模和预测那些以固定周期s比如s12代表月度数据s7代表周度数据s4代表季度数据重复出现的季节性规律。我最初在数模竞赛和商业分析项目中接触SARIMA时也被它复杂的参数搞晕过。但一旦理清思路你会发现它的预测能力提升是立竿见影的。这篇文章我就以一个从业者的角度带你彻底搞懂SARIMA模型的原理、核心参数的意义并手把手用Python完成从数据准备、模型识别、参数估计到预测评估的全流程。我们不止讲“怎么做”更重点剖析“为什么这么做”以及在实际操作中那些容易踩坑的细节。无论你是备战数模国赛还是处理实际的业务预测问题这篇内容都能给你一套可直接复现的可靠方案。2. 核心原理拆解SARIMA模型的“双重人格”要驾驭SARIMA模型必须理解它的数学表达和参数含义。很多人一看到公式就头疼我们不妨用更直观的方式来理解。2.1 模型数学表达与参数释义SARIMA模型的完整记法是 SARIMA(p, d, q)(P, D, Q)s。这看起来像一串密码我们把它拆开看(p, d, q)这是非季节性部分就是经典ARIMA模型的参数。p (自回归阶数)表示当前值受过去多少个时间点的自身值影响。比如p2意味着今天的销售额受到昨天和前天销售额的线性影响。它捕捉的是序列的“记忆”效应。d (差分阶数)为了使非平稳序列变得平稳而进行的差分次数。平稳性是时间序列建模的一个重要前提意味着序列的统计特性如均值、方差不随时间变化。一次差分就是计算相邻点的差值今天值-昨天值常用于消除线性趋势。q (移动平均阶数)表示当前值受过去多少个时间点的随机误差白噪声影响。它捕捉的是序列受到的突发性“冲击”效应。(P, D, Q)s这是季节性部分是SARIMA的精华所在。P (季节性自回归阶数)表示当前值受过去多少个季节性周期的自身值影响。例如对于月度数据(s12)P1意味着本月的销售额受到去年同月销售额的影响。D (季节性差分阶数)为了使序列的季节性成分变得平稳而进行的季节性差分次数。季节性差分是指计算当前值与上一个周期同期值的差值。例如对于月度数据一次季节性差分就是“本月值 - 去年同月值”这能有效消除固定的年度周期性趋势。Q (季节性移动平均阶数)表示当前值受过去多少个季节性周期的随机误差影响。s (季节周期长度)这是最关键的一个参数定义了季节性的周期。必须根据你的数据特性来指定。常见的有天数据s7每周周期月数据s12年度周期季度数据s4年度周期小时数据s24日周期SARIMA模型的数学公式是上述两部分算子的乘积作用于时间序列。简单理解模型认为任何一个时间点的值都是由“非季节性规律”和“季节性规律”共同作用再加上随机噪声的结果。建模的过程就是通过数据把这四组参数(p,d,q)和(P,D,Q)给估计出来。2.2 与ARIMA的核心区别与应用场景为了更清晰地看到SARIMA的升级之处我们把它和ARIMA放在一起对比特性维度ARIMA模型SARIMA模型核心能力处理趋势和非季节性自相关同时处理趋势、非季节性自相关和季节性自相关参数(p, d, q) 三组(p, d, q)(P, D, Q)s四组多一个季节周期s数据要求对周期性数据预测效果差必须用于具有明显、稳定周期性的数据模型复杂度相对简单参数少更复杂参数组合多计算量更大典型应用股价无固定周期、某些平滑的指标预测零售销售额周/月/年周期、能源负荷日/年周期、旅游客流节假日周期从应用场景来看当你面对的数据图表呈现出“波浪形”或“锯齿形”的规律重复时SARIMA就是你的首选工具。例如零售与电商预测未来几个月不同品类商品的周度销量必须考虑周末效应和年度促销节如双十一、黑色星期五的影响。能源行业预测未来一周的每小时电力负荷必须考虑工作日/休息日的用电模式差异以及日间的用电高峰低谷。交通物流预测机场未来数月的客流量必须考虑暑假、春运等季节性出行高峰。数模竞赛在国赛或美赛中遇到涉及经济指标、气象数据、传染病传播等有明显周期性的题目SARIMA是一个强有力的基准模型。注意SARIMA假设季节性模式是固定且稳定的。如果数据的季节性模式随着时间在缓慢变化例如由于消费习惯改变周末的销售峰值在逐年降低那么纯SARIMA的预测效果可能会打折扣可能需要结合其他模型或进行数据预处理。3. 实战前准备Python环境与数据理解理论说得再多不如动手跑一遍。在开始建模之前我们需要确保环境就绪并且真正理解手头的数据。3.1 工具链选型为什么是statsmodelsPython中进行时间序列分析有很多库比如pmdarima自动ARIMA、prophetFacebook开源等。但对于学习SARIMA的原理和进行精细控制我强烈推荐使用statsmodels库。原因有三权威与透明statsmodels是统计学领域的标准库之一其实现经过学术界和工业界的广泛检验。它输出的结果非常详细包括每个参数的估计值、标准误、t统计量和p值这对于我们理解模型、诊断问题至关重要。控制力强它允许你手动指定所有参数(p,d,q)(P,D,Q,s)而不是完全黑盒自动化。在数模竞赛或需要解释性的业务场景中这种控制力非常重要。生态完整它提供了完整的建模流程支持从平稳性检验ADF、白噪声检验Ljung-Box到模型诊断残差分析图一站式解决。当然pmdarima的auto_arima函数在快速寻找最优参数组合时非常方便可以作为辅助工具。但本文以理解为核心所以我们从statsmodels入手。安装非常简单如果你使用Anaconda它通常已经内置。如果没有一行命令搞定pip install statsmodels pandas numpy matplotlib seaborn这里一并安装了数据分析三件套pandas,numpy和绘图库matplotlib、seaborn。3.2 数据准备与探索性分析没有数据一切模型都是空中楼阁。我们使用一个经典的、具有明显季节性的数据集航空旅客数据AirPassengers。这个数据集记录了1949年至1960年每月航空旅客的数量是时间序列分析的“Hello World”。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf import warnings warnings.filterwarnings(ignore) # 忽略一些不影响运行的警告 sns.set_style(whitegrid) # 1. 加载数据 # 这里我们直接从statsmodels内置数据集获取 from statsmodels.datasets import get_rdataset data get_rdataset(AirPassengers) df data.data # 查看数据前几行 print(df.head()) print(f\n数据形状: {df.shape}) print(f时间范围: {df[time].min()} 到 {df[time].max()})运行后你会看到数据包含两列time年月和value旅客数。接下来是至关重要的可视化步骤这能帮你直观判断是否适合使用SARIMA。# 2. 将时间列转换为datetime格式并设为索引 df[date] pd.to_datetime(df[time].astype(str)) df.set_index(date, inplaceTrue) ts df[value] # 获取时间序列对象 # 3. 绘制原始序列图 plt.figure(figsize(14, 6)) plt.plot(ts, colorsteelblue, linewidth2) plt.title(AirPassengers Monthly Totals (1949-1960), fontsize15) plt.xlabel(Year) plt.ylabel(Number of Passengers (Thousands)) plt.grid(True, alpha0.3) plt.show()这张图会清晰地展示出四个关键特征明显上升趋势旅客总数随着年份增长而增加。稳定的年度季节性每年夏季年中出现一个高峰冬季年末年初是低谷周期为12个月。季节性振幅扩大随着趋势上升每年波峰和波谷的差值振幅似乎在变大。这提示数据可能具有乘性季节性季节波动的幅度与趋势水平成正比而非加性季节性波动幅度恒定。经典的SARIMA模型通常处理加性季节性对于乘性季节性需要对数据先取对数转换将其变为加性。非平稳性均值和方差随时间变化。为了确认季节性我们可以绘制季节分解图from statsmodels.tsa.seasonal import seasonal_decompose # 进行加法模型分解 result_add seasonal_decompose(ts, modeladditive, period12) # 进行乘法模型分解更可能适合本例 result_mul seasonal_decompose(ts, modelmultiplicative, period12) fig, axes plt.subplots(4, 2, figsize(16, 12)) 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()通过对比你会发现乘法模型的残差Residual更随机、更接近白噪声而加法模型的残差在序列后期仍有明显的周期性。这证实了我们的猜测原始数据具有乘性季节性。因此标准的建模流程是先对原始数据取自然对数将乘性关系转化为加性关系再用SARIMA建模最后将预测值指数变换回去。# 对原序列取对数以稳定方差将乘性季节性转为加性 ts_log np.log(ts) plt.figure(figsize(14, 6)) plt.subplot(1,2,1) plt.plot(ts, colorsteelblue) plt.title(Original Series) plt.subplot(1,2,2) plt.plot(ts_log, colorcoral) plt.title(Log-Transformed Series) plt.tight_layout() plt.show()观察取对数后的序列图你会发现上升趋势依然存在但季节性波动的幅度变得相对均匀了这正是我们想要的。4. 模型识别与定阶寻找最优参数组合这是SARIMA建模中最具技巧性的一步。我们需要确定(p,d,q)和(P,D,Q,s)这7个参数。s我们已经知道是12月度数据。剩下的6个参数我们通过“平稳性检验”确定d和D通过“自相关图(ACF)和偏自相关图(PACF)”初步判断p,q和P,Q。4.1 平稳性检验与差分阶数确定首先检验取对数后的序列ts_log是否平稳。使用Augmented Dickey-Fuller (ADF)检验。def adf_test(timeseries): 执行ADF检验并打印结果 print(Results of Dickey-Fuller Test:) dftest adfuller(timeseries, autolagAIC) # AIC准则自动选择滞后阶数 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[Critical Value (%s) % key] value print(dfoutput) # 判断 if dfoutput[p-value] 0.05: print(\n结论序列是平稳的 (p-value 0.05)) else: print(\n结论序列是非平稳的 (p-value 0.05)) print(对原始对数序列进行ADF检验) adf_test(ts_log)结果大概率显示p-value大于0.05说明ts_log非平稳。我们需要进行差分。先尝试一阶普通差分消除趋势ts_log_diff1 ts_log.diff().dropna() # 一阶差分并删除NaN print(\n对一阶差分后的序列进行ADF检验) adf_test(ts_log_diff1)如果一阶差分后序列平稳了那么d1。如果还不平稳可能需要二阶差分(d2)但实践中d通常为0,1,2。接下来处理季节性。我们对已经进行了一阶普通差分的序列ts_log_diff1再进行一次周期为12的季节性差分以消除季节性不平稳ts_log_diff1_seasonal ts_log_diff1.diff(12).dropna() # 周期s12的季节性差分 print(\n对一阶差分季节性差分后的序列进行ADF检验) adf_test(ts_log_diff1_seasonal)如果这个序列平稳了那么D1。对于这个数据集通常d1,D1就能达到平稳。实操心得差分不是越多越好。过差分over-differencing虽然能让序列变得平稳但会引入不必要的噪声和结构降低模型预测精度并且使ACF/PACF图难以解读。原则是用最少的差分次数使序列达到平稳。通常d和D都不会超过2。4.2 利用ACF/PACF图初步定阶在获得平稳序列例如ts_log_diff1_seasonal后我们绘制它的自相关函数ACF和偏自相关函数PACF图。这两个图是确定p, q, P, Q的主要工具。ACF图描述当前序列与自身滞后k期序列之间的相关性。PACF图在控制了中间滞后项1,2,...,k-1期的影响后当前序列与滞后k期序列之间的纯相关性。fig, axes plt.subplots(1, 2, figsize(16, 4)) # 绘制平稳序列的ACF和PACF图 plot_acf(ts_log_diff1_seasonal, lags40, axaxes[0]) # 看40个滞后 plot_pacf(ts_log_diff1_seasonal, lags40, axaxes[1], methodywm) # 推荐使用ywm方法 axes[0].set_title(ACF of Stationary Series (Diff Seasonal Diff)) axes[1].set_title(PACF of Stationary Series (Diff Seasonal Diff)) plt.tight_layout() plt.show()如何解读确定非季节性参数(p, q)p (AR阶数)看PACF图。PACF在滞后p阶后出现截尾即突然下降到置信区间内像被砍断一样那么p就可能是这个值。图中可能在滞后1阶或2阶后截尾。q (MA阶数)看ACF图。ACF在滞后q阶后出现截尾那么q就可能是这个值。图中可能在滞后1阶后截尾。确定季节性参数(P, Q)观察ACF和PACF在季节周期倍数处即滞后12, 24, 36...的峰值。如果ACF在滞后12阶处有一个显著的峰值并且缓慢衰减可能暗示季节性自回归成分(P)。如果PACF在滞后12阶处有一个显著的峰值并且缓慢衰减可能暗示季节性移动平均成分(Q)。更常见的是在滞后12阶处ACF和PACF都出现一个显著的“钉子”峰值这通常意味着需要一阶季节性差分(D1)并且可能还需要一个季节性自回归项(P1)或季节性移动平均项(Q1)。观察我们的图可能会发现ACF在滞后1阶后截尾提示q1。PACF在滞后1阶或2阶后截尾提示p1或p2。在滞后12阶、24阶处ACF和PACF都有显著峰值提示需要P1和/或Q1。这是一种基于经验的“看图说话”并不精确。因此我们通常基于此确定一个参数搜索范围然后通过网格搜索Grid Search配合信息准则如AIC来找到最优组合。5. 模型拟合、诊断与预测有了初步的参数范围我们就可以开始拟合模型并进行严格的诊断以确保模型是合适的。5.1 网格搜索与AIC准则选优我们使用statsmodels.tsa.statespace.sarimax.SARIMAX来构建模型。为了避免手动尝试所有组合我们写一个简单的网格搜索。这里我们搜索一个较小的空间作为示例。import itertools from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings(ignore) # 忽略拟合过程中的警告 # 定义参数搜索范围 p d q range(0, 2) # 假设p,d,q在0,1中选 P D Q range(0, 2) # 假设P,D,Q在0,1中选 s 12 # 季节周期固定为12 # 生成所有参数组合 pdq list(itertools.product(p, d, q)) seasonal_pdq list(itertools.product(P, D, Q, [s])) best_aic float(inf) best_order None best_seasonal_order None best_model None print(开始网格搜索...) for param in pdq: for param_seasonal in seasonal_pdq: try: # 注意这里我们用原始的对数序列 ts_log 进行拟合 # SARIMAX 内部会自己根据 order 和 seasonal_order 进行差分 model SARIMAX(ts_log, orderparam, seasonal_orderparam_seasonal, enforce_stationarityFalse, # 让模型自己处理平稳性 enforce_invertibilityFalse) # 让模型自己处理可逆性 results model.fit(dispFalse) # dispFalse不显示迭代信息 current_aic results.aic if current_aic best_aic: best_aic current_aic best_order param best_seasonal_order param_seasonal best_model results print(f当前最优: SARIMA{param}x{param_seasonal} - AIC:{current_aic:.2f}) except Exception as e: # 某些参数组合可能导致无法估计跳过 # print(f参数{param}x{param_seasonal}拟合失败: {e}) continue print(f\n 网格搜索结束 ) print(f最优参数组合: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC值: {best_aic:.2f})AICAkaike Information Criterion准则衡量的是模型的拟合优度和复杂度的平衡AIC值越小越好。运行这段代码需要一些时间因为它要拟合多个模型。对于我们的示例数据一个常见的最优结果是SARIMA(1,1,1)x(1,1,1,12)。5.2 模型拟合结果解读与诊断用最优参数拟合模型并查看详细的摘要报告。# 使用最优参数重新拟合模型为了获得完整输出 final_model SARIMAX(ts_log, orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_results final_model.fit(dispTrue) # dispTrue显示迭代过程 print(final_results.summary())摘要报告非常丰富重点关注以下几点系数表coef查看每个估计参数ar.L1, ma.L1, ar.S.L12等的值、标准误、z统计量和P|z|。P值小于0.05通常认为该系数显著不为零。如果某个系数的P值很大比如0.1可以考虑在简化模型时去掉它。信息准则AIC, BIC, HQIC的值用于模型比较。残差诊断更重要的步骤是可视化诊断。一个好的模型其残差应该类似于白噪声均值为0方差恒定无自相关。# 绘制诊断图 final_results.plot_diagnostics(figsize(15, 10)) plt.tight_layout() plt.show()诊断图包含四个子图标准化残差图残差应该围绕0随机波动没有明显的趋势或周期性。残差直方图核密度估计理想情况下应与正态分布曲线红色基本吻合检验残差是否服从正态分布。正态Q-Q图点应大致分布在红色直线上检验残差的正态性。残差自相关图Correlogram所有滞后阶数的自相关系数都应落在蓝色阴影区域置信区间内表明残差无自相关是白噪声。如果诊断图不理想例如残差自相关图在滞后1阶或季节滞后处仍有显著峰值说明模型未能完全捕捉数据中的信息可能需要调整参数。5.3 样本内拟合与样本外预测模型诊断通过后我们就可以进行预测了。首先看看模型在历史数据上的拟合效果样本内拟合。# 获取样本内预测值动态预测使用直到该时间点的所有信息 pred_dynamic final_results.get_prediction(startpd.to_datetime(1958-01-01), dynamicFalse, full_resultsTrue) pred_dynamic_ci pred_dynamic.conf_int() # 获取置信区间 # 因为我们对对数序列建模所以预测值也是对数需要指数变换回去 pred_dynamic_values np.exp(pred_dynamic.predicted_mean) actual_values ts[1958-01-01:] # 绘制拟合效果图 plt.figure(figsize(14, 7)) plt.plot(ts, labelObserved, colorblue, alpha0.6) plt.plot(pred_dynamic_values, labelDynamic Forecast (1958 onwards), colorred, alpha0.8, linestyle--) plt.fill_between(pred_dynamic_ci.index, np.exp(pred_dynamic_ci.iloc[:, 0]), np.exp(pred_dynamic_ci.iloc[:, 1]), colorr, alpha0.1) plt.title(SARIMA Model - In-Sample Forecast vs Actuals) plt.xlabel(Date) plt.ylabel(AirPassengers) plt.legend() plt.show()接下来进行真正的样本外预测预测未来一段时间例如未来24个月的值。# 预测未来24个月2年 forecast_steps 24 # 获取预测值及置信区间对数尺度 forecast_object final_results.get_forecast(stepsforecast_steps) forecast_values_log forecast_object.predicted_mean forecast_ci_log forecast_object.conf_int() # 指数变换回到原始尺度 forecast_values np.exp(forecast_values_log) forecast_ci np.exp(forecast_ci_log) # 生成未来日期索引 last_date ts.index[-1] future_dates pd.date_range(startlast_date pd.DateOffset(months1), periodsforecast_steps, freqMS) # 绘制最终预测图 plt.figure(figsize(14, 8)) plt.plot(ts, labelHistorical Data, colorblue) plt.plot(future_dates, forecast_values, labelForecast, colorred, linestyle--) plt.fill_between(future_dates, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorr, alpha0.15, label95% Confidence Interval) plt.title(fSARIMA{best_order}x{best_seasonal_order} - {forecast_steps} Months Forecast) plt.xlabel(Date) plt.ylabel(AirPassengers) plt.legend() plt.grid(True, alpha0.3) plt.show() # 也可以打印出具体的预测数值 forecast_df pd.DataFrame({ Date: future_dates, Forecast: forecast_values.values, Lower CI: forecast_ci.iloc[:, 0].values, Upper CI: forecast_ci.iloc[:, 1].values }) print(forecast_df.head(10))观察预测图你应该能看到模型不仅延续了上升趋势还完美地复现了年度季节性波动。置信区间给出了预测的不确定性范围。6. 避坑指南与高级技巧在实际项目中严格按照上述流程走一遍你大概率能得到一个可用的SARIMA模型。但要想让模型更稳健、预测更准还需要注意以下这些我踩过的坑和总结的技巧。6.1 参数选择中的常见陷阱过度依赖自动定阶工具像pmdarima.auto_arima这样的工具非常方便但它给出的“最优”模型有时在业务上解释性不强或者因为陷入局部最优而忽略了更简单的模型。最佳实践是用自动工具给出建议再结合ACF/PACF图和自己对业务的理解进行微调。一个AIC值稍高但参数更简洁、更易解释的模型往往比一个复杂模型更受业务方欢迎。忽略残差诊断拟合完模型看一眼AIC不错就直接预测这是大忌。必须做残差诊断。如果残差不是白噪声说明还有信息未被模型提取预测就会有偏差。常见的残差问题包括残差自相关在ACF图上有显著滞后。解决方案尝试增加p或q或P,Q。残差异方差残差方差随时间增大。解决方案检查是否需要对数据做更复杂的变换如Box-Cox变换而不仅仅是取对数。残差非正态这不一定致命但可能影响预测区间的准确性。对于点预测影响不大但对于需要概率预测的场景需要注意。季节性周期s判断错误这是最致命的错误。如果周期设错整个季节性部分就全错了。务必通过业务知识如周、月、季、年和数据可视化绘制多年数据叠加图双重确认s的值。对于多重季节性如同时存在周周期和年周期标准SARIMA无法处理需要考虑其他模型如TBATS或Prophet。6.2 处理外生变量与节假日效应标准的SARIMA是单变量模型只利用自身的历史数据进行预测。但在现实中很多序列受外部因素影响。例如销售额受促销活动影响客流量受天气影响。这时可以使用SARIMAX模型X代表外生变量。statsmodels的SARIMAX类本身就支持外生变量。你需要准备一个与外生变量对应的DataFrame然后在fit方法中传入exog参数。# 假设我们有一个外生变量‘promotion’表示是否有促销0/1 exog_data pd.DataFrame(...) # 需要与ts_log具有相同的索引和长度 model_with_exog SARIMAX(ts_log, order(1,1,1), seasonal_order(1,1,1,12), exogexog_data) # 传入外生变量 results_with_exog model_with_exog.fit()预测时也需要提供未来期的外生变量值。对于节假日效应一种常见做法是将其作为虚拟变量哑变量加入外生变量。例如为春节所在的月份创建一个变量值为1其他月份为0。6.3 模型评估与滚动预测在正式使用模型预测未来之前必须评估其预测性能。不能只看在全部历史数据上的拟合效果这会导致过拟合乐观估计。正确的方法是使用时间序列交叉验证或滚动预测。思路是在历史数据上模拟一个不断向前滚动的预测过程。用前N个月的数据训练一个SARIMA模型。预测未来M个月。将预测值与真实值比较计算误差如均方根误差RMSE平均绝对百分比误差MAPE。将训练窗口向后移动M个月或1个月重复步骤1-3。将所有滚动窗口的误差平均得到模型在“未知”数据上的性能估计。from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error def rolling_forecast(series, order, seasonal_order, train_size, test_size, s): 执行滚动预测评估 series: 时间序列对数变换后 train_size: 初始训练集大小 test_size: 每次预测的步长通常为季节周期s的倍数 history list(series[:train_size]) predictions [] actuals [] for t in range(train_size, len(series), test_size): # 用当前历史数据拟合模型 model SARIMAX(history, orderorder, seasonal_orderseasonal_order, enforce_stationarityFalse) model_fit model.fit(dispFalse) # 预测未来test_size步 forecast model_fit.forecast(stepstest_size) predictions.extend(forecast) actuals.extend(series[t:ttest_size].values) # 将真实值更新到历史中模拟时间推移 history.extend(series[t:ttest_size].values) # 计算误差注意预测值是对数需要转换回原始尺度再与原始尺度真实值比较 # 更合理的做法将对数尺度的预测值和真实值都指数变换回原始尺度再计算误差。 pred_original np.exp(predictions) actual_original np.exp(actuals) rmse np.sqrt(mean_squared_error(actual_original, pred_original)) mape mean_absolute_percentage_error(actual_original, pred_original) * 100 return pred_original, actual_original, rmse, mape # 示例使用最后3年36个月作为测试期每次滚动预测12个月 train_len len(ts_log) - 36 preds, actuals, rmse_val, mape_val rolling_forecast(ts_log, orderbest_order, seasonal_orderbest_seasonal_order, train_sizetrain_len, test_size12, # 每次预测一年 s12) print(f滚动预测RMSE: {rmse_val:.2f}) print(f滚动预测MAPE: {mape_val:.2f}%)MAPE平均绝对百分比误差是一个很好的业务指标比如MAPE5%意味着平均预测误差在5%左右。这个评估结果比样本内拟合的R²更有说服力。7. 在数模竞赛与业务中的实战要点最后结合我参加数模和做业务预测的经验分享几点SARIMA模型的实战心得。在数模竞赛中基础模型高级呈现SARIMA是一个很好的基准模型。在论文中你需要清晰地展示数据平稳化处理取对数、差分的过程图、ACF/PACF定阶分析图、模型诊断图、预测效果对比图。图表质量直接影响评分。结合其他模型不要只用一个SARIMA。可以将其与指数平滑ETS、线性回归考虑趋势和季节虚拟变量甚至简单的移动平均进行对比。用RMSE、MAPE等指标说明为什么SARIMA更优。强调假设检验在论文中写明你检查了残差的白噪声性Ljung-Box检验、正态性Jarque-Bera检验等证明模型是充分的。这体现了建模的严谨性。处理预测区间SARIMA可以给出预测的置信区间。在论文中画出这个区间并讨论其意义例如“我们有95%的把握认为未来销售额将落在X到Y之间”这比只给一个点预测值要高级得多。在业务分析中自动化与监控业务预测需要定期更新如每周/每月。可以编写自动化脚本定期用新数据重新训练模型并生成预测报告。同时监控模型的预测误差如果MAPE持续上升说明模型可能失效需要重新调整或更换。理解业务上限SARIMA是统计模型它基于历史模式外推。如果业务发生了根本性变化如新产品上市、政策巨变历史模式将不再适用模型的预测会失灵。这时需要结合业务判断或引入能够捕捉突变的外生变量。解释性优先在商业场景一个SARIMA(0,1,1)x(0,1,1,12)模型可能比一个SARIMA(2,1,2)x(1,1,1,12)模型更受欢迎即使后者AIC稍低因为前者更简洁更容易向非技术人员解释“模型认为下个月的值主要受上个月的误差和去年同月的值影响”。SARIMA模型是一个强大的工具但它不是银弹。它最适合于具有稳定趋势和固定季节性的中短期预测。对于长期预测不确定性会急剧增加对于模式快速变化的数据需要更灵活的模型。掌握它意味着你拥有了时间序列预测武器库中一件经典且可靠的重武器。从理解原理到代码实现再到避开实战中的各种坑这个过程本身就是数据分析能力的一次扎实提升。 SARIMA模型实战：从原理到Python代码实现季节性时间序列预测
1. 项目概述从ARIMA到SARIMA的跃迁如果你做过时间序列预测大概率绕不开ARIMA模型。它经典、强大是很多数据分析师和建模初学者的首选。但当你处理的数据呈现出明显的周期性波动时——比如电商网站的日销售额在周末飙升电力负荷在夏季和冬季形成高峰或者某种商品的月度销量在每年春节前后规律性上涨——标准的ARIMA模型就会显得力不从心。它捕捉不到这种固定周期的重复模式预测的曲线会变得平滑丢失掉关键的季节性信息。这时你就需要请出它的升级版季节性自回归积分滑动平均模型也就是SARIMA。SARIMA模型可以看作是ARIMA模型在时间序列分析领域的“完全体”。它在ARIMA的三个核心参数自回归阶数p、差分阶数d、移动平均阶数q之外额外引入了一套专门处理季节性成分的参数P, D, Q, s。这套组合拳让模型具备了“双线作战”的能力一套参数p,d,q负责捕捉数据中非季节性的趋势和随机波动另一套参数P,D,Q,s则专门用来建模和预测那些以固定周期s比如s12代表月度数据s7代表周度数据s4代表季度数据重复出现的季节性规律。我最初在数模竞赛和商业分析项目中接触SARIMA时也被它复杂的参数搞晕过。但一旦理清思路你会发现它的预测能力提升是立竿见影的。这篇文章我就以一个从业者的角度带你彻底搞懂SARIMA模型的原理、核心参数的意义并手把手用Python完成从数据准备、模型识别、参数估计到预测评估的全流程。我们不止讲“怎么做”更重点剖析“为什么这么做”以及在实际操作中那些容易踩坑的细节。无论你是备战数模国赛还是处理实际的业务预测问题这篇内容都能给你一套可直接复现的可靠方案。2. 核心原理拆解SARIMA模型的“双重人格”要驾驭SARIMA模型必须理解它的数学表达和参数含义。很多人一看到公式就头疼我们不妨用更直观的方式来理解。2.1 模型数学表达与参数释义SARIMA模型的完整记法是 SARIMA(p, d, q)(P, D, Q)s。这看起来像一串密码我们把它拆开看(p, d, q)这是非季节性部分就是经典ARIMA模型的参数。p (自回归阶数)表示当前值受过去多少个时间点的自身值影响。比如p2意味着今天的销售额受到昨天和前天销售额的线性影响。它捕捉的是序列的“记忆”效应。d (差分阶数)为了使非平稳序列变得平稳而进行的差分次数。平稳性是时间序列建模的一个重要前提意味着序列的统计特性如均值、方差不随时间变化。一次差分就是计算相邻点的差值今天值-昨天值常用于消除线性趋势。q (移动平均阶数)表示当前值受过去多少个时间点的随机误差白噪声影响。它捕捉的是序列受到的突发性“冲击”效应。(P, D, Q)s这是季节性部分是SARIMA的精华所在。P (季节性自回归阶数)表示当前值受过去多少个季节性周期的自身值影响。例如对于月度数据(s12)P1意味着本月的销售额受到去年同月销售额的影响。D (季节性差分阶数)为了使序列的季节性成分变得平稳而进行的季节性差分次数。季节性差分是指计算当前值与上一个周期同期值的差值。例如对于月度数据一次季节性差分就是“本月值 - 去年同月值”这能有效消除固定的年度周期性趋势。Q (季节性移动平均阶数)表示当前值受过去多少个季节性周期的随机误差影响。s (季节周期长度)这是最关键的一个参数定义了季节性的周期。必须根据你的数据特性来指定。常见的有天数据s7每周周期月数据s12年度周期季度数据s4年度周期小时数据s24日周期SARIMA模型的数学公式是上述两部分算子的乘积作用于时间序列。简单理解模型认为任何一个时间点的值都是由“非季节性规律”和“季节性规律”共同作用再加上随机噪声的结果。建模的过程就是通过数据把这四组参数(p,d,q)和(P,D,Q)给估计出来。2.2 与ARIMA的核心区别与应用场景为了更清晰地看到SARIMA的升级之处我们把它和ARIMA放在一起对比特性维度ARIMA模型SARIMA模型核心能力处理趋势和非季节性自相关同时处理趋势、非季节性自相关和季节性自相关参数(p, d, q) 三组(p, d, q)(P, D, Q)s四组多一个季节周期s数据要求对周期性数据预测效果差必须用于具有明显、稳定周期性的数据模型复杂度相对简单参数少更复杂参数组合多计算量更大典型应用股价无固定周期、某些平滑的指标预测零售销售额周/月/年周期、能源负荷日/年周期、旅游客流节假日周期从应用场景来看当你面对的数据图表呈现出“波浪形”或“锯齿形”的规律重复时SARIMA就是你的首选工具。例如零售与电商预测未来几个月不同品类商品的周度销量必须考虑周末效应和年度促销节如双十一、黑色星期五的影响。能源行业预测未来一周的每小时电力负荷必须考虑工作日/休息日的用电模式差异以及日间的用电高峰低谷。交通物流预测机场未来数月的客流量必须考虑暑假、春运等季节性出行高峰。数模竞赛在国赛或美赛中遇到涉及经济指标、气象数据、传染病传播等有明显周期性的题目SARIMA是一个强有力的基准模型。注意SARIMA假设季节性模式是固定且稳定的。如果数据的季节性模式随着时间在缓慢变化例如由于消费习惯改变周末的销售峰值在逐年降低那么纯SARIMA的预测效果可能会打折扣可能需要结合其他模型或进行数据预处理。3. 实战前准备Python环境与数据理解理论说得再多不如动手跑一遍。在开始建模之前我们需要确保环境就绪并且真正理解手头的数据。3.1 工具链选型为什么是statsmodelsPython中进行时间序列分析有很多库比如pmdarima自动ARIMA、prophetFacebook开源等。但对于学习SARIMA的原理和进行精细控制我强烈推荐使用statsmodels库。原因有三权威与透明statsmodels是统计学领域的标准库之一其实现经过学术界和工业界的广泛检验。它输出的结果非常详细包括每个参数的估计值、标准误、t统计量和p值这对于我们理解模型、诊断问题至关重要。控制力强它允许你手动指定所有参数(p,d,q)(P,D,Q,s)而不是完全黑盒自动化。在数模竞赛或需要解释性的业务场景中这种控制力非常重要。生态完整它提供了完整的建模流程支持从平稳性检验ADF、白噪声检验Ljung-Box到模型诊断残差分析图一站式解决。当然pmdarima的auto_arima函数在快速寻找最优参数组合时非常方便可以作为辅助工具。但本文以理解为核心所以我们从statsmodels入手。安装非常简单如果你使用Anaconda它通常已经内置。如果没有一行命令搞定pip install statsmodels pandas numpy matplotlib seaborn这里一并安装了数据分析三件套pandas,numpy和绘图库matplotlib、seaborn。3.2 数据准备与探索性分析没有数据一切模型都是空中楼阁。我们使用一个经典的、具有明显季节性的数据集航空旅客数据AirPassengers。这个数据集记录了1949年至1960年每月航空旅客的数量是时间序列分析的“Hello World”。import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf import warnings warnings.filterwarnings(ignore) # 忽略一些不影响运行的警告 sns.set_style(whitegrid) # 1. 加载数据 # 这里我们直接从statsmodels内置数据集获取 from statsmodels.datasets import get_rdataset data get_rdataset(AirPassengers) df data.data # 查看数据前几行 print(df.head()) print(f\n数据形状: {df.shape}) print(f时间范围: {df[time].min()} 到 {df[time].max()})运行后你会看到数据包含两列time年月和value旅客数。接下来是至关重要的可视化步骤这能帮你直观判断是否适合使用SARIMA。# 2. 将时间列转换为datetime格式并设为索引 df[date] pd.to_datetime(df[time].astype(str)) df.set_index(date, inplaceTrue) ts df[value] # 获取时间序列对象 # 3. 绘制原始序列图 plt.figure(figsize(14, 6)) plt.plot(ts, colorsteelblue, linewidth2) plt.title(AirPassengers Monthly Totals (1949-1960), fontsize15) plt.xlabel(Year) plt.ylabel(Number of Passengers (Thousands)) plt.grid(True, alpha0.3) plt.show()这张图会清晰地展示出四个关键特征明显上升趋势旅客总数随着年份增长而增加。稳定的年度季节性每年夏季年中出现一个高峰冬季年末年初是低谷周期为12个月。季节性振幅扩大随着趋势上升每年波峰和波谷的差值振幅似乎在变大。这提示数据可能具有乘性季节性季节波动的幅度与趋势水平成正比而非加性季节性波动幅度恒定。经典的SARIMA模型通常处理加性季节性对于乘性季节性需要对数据先取对数转换将其变为加性。非平稳性均值和方差随时间变化。为了确认季节性我们可以绘制季节分解图from statsmodels.tsa.seasonal import seasonal_decompose # 进行加法模型分解 result_add seasonal_decompose(ts, modeladditive, period12) # 进行乘法模型分解更可能适合本例 result_mul seasonal_decompose(ts, modelmultiplicative, period12) fig, axes plt.subplots(4, 2, figsize(16, 12)) 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()通过对比你会发现乘法模型的残差Residual更随机、更接近白噪声而加法模型的残差在序列后期仍有明显的周期性。这证实了我们的猜测原始数据具有乘性季节性。因此标准的建模流程是先对原始数据取自然对数将乘性关系转化为加性关系再用SARIMA建模最后将预测值指数变换回去。# 对原序列取对数以稳定方差将乘性季节性转为加性 ts_log np.log(ts) plt.figure(figsize(14, 6)) plt.subplot(1,2,1) plt.plot(ts, colorsteelblue) plt.title(Original Series) plt.subplot(1,2,2) plt.plot(ts_log, colorcoral) plt.title(Log-Transformed Series) plt.tight_layout() plt.show()观察取对数后的序列图你会发现上升趋势依然存在但季节性波动的幅度变得相对均匀了这正是我们想要的。4. 模型识别与定阶寻找最优参数组合这是SARIMA建模中最具技巧性的一步。我们需要确定(p,d,q)和(P,D,Q,s)这7个参数。s我们已经知道是12月度数据。剩下的6个参数我们通过“平稳性检验”确定d和D通过“自相关图(ACF)和偏自相关图(PACF)”初步判断p,q和P,Q。4.1 平稳性检验与差分阶数确定首先检验取对数后的序列ts_log是否平稳。使用Augmented Dickey-Fuller (ADF)检验。def adf_test(timeseries): 执行ADF检验并打印结果 print(Results of Dickey-Fuller Test:) dftest adfuller(timeseries, autolagAIC) # AIC准则自动选择滞后阶数 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[Critical Value (%s) % key] value print(dfoutput) # 判断 if dfoutput[p-value] 0.05: print(\n结论序列是平稳的 (p-value 0.05)) else: print(\n结论序列是非平稳的 (p-value 0.05)) print(对原始对数序列进行ADF检验) adf_test(ts_log)结果大概率显示p-value大于0.05说明ts_log非平稳。我们需要进行差分。先尝试一阶普通差分消除趋势ts_log_diff1 ts_log.diff().dropna() # 一阶差分并删除NaN print(\n对一阶差分后的序列进行ADF检验) adf_test(ts_log_diff1)如果一阶差分后序列平稳了那么d1。如果还不平稳可能需要二阶差分(d2)但实践中d通常为0,1,2。接下来处理季节性。我们对已经进行了一阶普通差分的序列ts_log_diff1再进行一次周期为12的季节性差分以消除季节性不平稳ts_log_diff1_seasonal ts_log_diff1.diff(12).dropna() # 周期s12的季节性差分 print(\n对一阶差分季节性差分后的序列进行ADF检验) adf_test(ts_log_diff1_seasonal)如果这个序列平稳了那么D1。对于这个数据集通常d1,D1就能达到平稳。实操心得差分不是越多越好。过差分over-differencing虽然能让序列变得平稳但会引入不必要的噪声和结构降低模型预测精度并且使ACF/PACF图难以解读。原则是用最少的差分次数使序列达到平稳。通常d和D都不会超过2。4.2 利用ACF/PACF图初步定阶在获得平稳序列例如ts_log_diff1_seasonal后我们绘制它的自相关函数ACF和偏自相关函数PACF图。这两个图是确定p, q, P, Q的主要工具。ACF图描述当前序列与自身滞后k期序列之间的相关性。PACF图在控制了中间滞后项1,2,...,k-1期的影响后当前序列与滞后k期序列之间的纯相关性。fig, axes plt.subplots(1, 2, figsize(16, 4)) # 绘制平稳序列的ACF和PACF图 plot_acf(ts_log_diff1_seasonal, lags40, axaxes[0]) # 看40个滞后 plot_pacf(ts_log_diff1_seasonal, lags40, axaxes[1], methodywm) # 推荐使用ywm方法 axes[0].set_title(ACF of Stationary Series (Diff Seasonal Diff)) axes[1].set_title(PACF of Stationary Series (Diff Seasonal Diff)) plt.tight_layout() plt.show()如何解读确定非季节性参数(p, q)p (AR阶数)看PACF图。PACF在滞后p阶后出现截尾即突然下降到置信区间内像被砍断一样那么p就可能是这个值。图中可能在滞后1阶或2阶后截尾。q (MA阶数)看ACF图。ACF在滞后q阶后出现截尾那么q就可能是这个值。图中可能在滞后1阶后截尾。确定季节性参数(P, Q)观察ACF和PACF在季节周期倍数处即滞后12, 24, 36...的峰值。如果ACF在滞后12阶处有一个显著的峰值并且缓慢衰减可能暗示季节性自回归成分(P)。如果PACF在滞后12阶处有一个显著的峰值并且缓慢衰减可能暗示季节性移动平均成分(Q)。更常见的是在滞后12阶处ACF和PACF都出现一个显著的“钉子”峰值这通常意味着需要一阶季节性差分(D1)并且可能还需要一个季节性自回归项(P1)或季节性移动平均项(Q1)。观察我们的图可能会发现ACF在滞后1阶后截尾提示q1。PACF在滞后1阶或2阶后截尾提示p1或p2。在滞后12阶、24阶处ACF和PACF都有显著峰值提示需要P1和/或Q1。这是一种基于经验的“看图说话”并不精确。因此我们通常基于此确定一个参数搜索范围然后通过网格搜索Grid Search配合信息准则如AIC来找到最优组合。5. 模型拟合、诊断与预测有了初步的参数范围我们就可以开始拟合模型并进行严格的诊断以确保模型是合适的。5.1 网格搜索与AIC准则选优我们使用statsmodels.tsa.statespace.sarimax.SARIMAX来构建模型。为了避免手动尝试所有组合我们写一个简单的网格搜索。这里我们搜索一个较小的空间作为示例。import itertools from statsmodels.tsa.statespace.sarimax import SARIMAX import warnings warnings.filterwarnings(ignore) # 忽略拟合过程中的警告 # 定义参数搜索范围 p d q range(0, 2) # 假设p,d,q在0,1中选 P D Q range(0, 2) # 假设P,D,Q在0,1中选 s 12 # 季节周期固定为12 # 生成所有参数组合 pdq list(itertools.product(p, d, q)) seasonal_pdq list(itertools.product(P, D, Q, [s])) best_aic float(inf) best_order None best_seasonal_order None best_model None print(开始网格搜索...) for param in pdq: for param_seasonal in seasonal_pdq: try: # 注意这里我们用原始的对数序列 ts_log 进行拟合 # SARIMAX 内部会自己根据 order 和 seasonal_order 进行差分 model SARIMAX(ts_log, orderparam, seasonal_orderparam_seasonal, enforce_stationarityFalse, # 让模型自己处理平稳性 enforce_invertibilityFalse) # 让模型自己处理可逆性 results model.fit(dispFalse) # dispFalse不显示迭代信息 current_aic results.aic if current_aic best_aic: best_aic current_aic best_order param best_seasonal_order param_seasonal best_model results print(f当前最优: SARIMA{param}x{param_seasonal} - AIC:{current_aic:.2f}) except Exception as e: # 某些参数组合可能导致无法估计跳过 # print(f参数{param}x{param_seasonal}拟合失败: {e}) continue print(f\n 网格搜索结束 ) print(f最优参数组合: SARIMA{best_order}x{best_seasonal_order}) print(f最优AIC值: {best_aic:.2f})AICAkaike Information Criterion准则衡量的是模型的拟合优度和复杂度的平衡AIC值越小越好。运行这段代码需要一些时间因为它要拟合多个模型。对于我们的示例数据一个常见的最优结果是SARIMA(1,1,1)x(1,1,1,12)。5.2 模型拟合结果解读与诊断用最优参数拟合模型并查看详细的摘要报告。# 使用最优参数重新拟合模型为了获得完整输出 final_model SARIMAX(ts_log, orderbest_order, seasonal_orderbest_seasonal_order, enforce_stationarityFalse, enforce_invertibilityFalse) final_results final_model.fit(dispTrue) # dispTrue显示迭代过程 print(final_results.summary())摘要报告非常丰富重点关注以下几点系数表coef查看每个估计参数ar.L1, ma.L1, ar.S.L12等的值、标准误、z统计量和P|z|。P值小于0.05通常认为该系数显著不为零。如果某个系数的P值很大比如0.1可以考虑在简化模型时去掉它。信息准则AIC, BIC, HQIC的值用于模型比较。残差诊断更重要的步骤是可视化诊断。一个好的模型其残差应该类似于白噪声均值为0方差恒定无自相关。# 绘制诊断图 final_results.plot_diagnostics(figsize(15, 10)) plt.tight_layout() plt.show()诊断图包含四个子图标准化残差图残差应该围绕0随机波动没有明显的趋势或周期性。残差直方图核密度估计理想情况下应与正态分布曲线红色基本吻合检验残差是否服从正态分布。正态Q-Q图点应大致分布在红色直线上检验残差的正态性。残差自相关图Correlogram所有滞后阶数的自相关系数都应落在蓝色阴影区域置信区间内表明残差无自相关是白噪声。如果诊断图不理想例如残差自相关图在滞后1阶或季节滞后处仍有显著峰值说明模型未能完全捕捉数据中的信息可能需要调整参数。5.3 样本内拟合与样本外预测模型诊断通过后我们就可以进行预测了。首先看看模型在历史数据上的拟合效果样本内拟合。# 获取样本内预测值动态预测使用直到该时间点的所有信息 pred_dynamic final_results.get_prediction(startpd.to_datetime(1958-01-01), dynamicFalse, full_resultsTrue) pred_dynamic_ci pred_dynamic.conf_int() # 获取置信区间 # 因为我们对对数序列建模所以预测值也是对数需要指数变换回去 pred_dynamic_values np.exp(pred_dynamic.predicted_mean) actual_values ts[1958-01-01:] # 绘制拟合效果图 plt.figure(figsize(14, 7)) plt.plot(ts, labelObserved, colorblue, alpha0.6) plt.plot(pred_dynamic_values, labelDynamic Forecast (1958 onwards), colorred, alpha0.8, linestyle--) plt.fill_between(pred_dynamic_ci.index, np.exp(pred_dynamic_ci.iloc[:, 0]), np.exp(pred_dynamic_ci.iloc[:, 1]), colorr, alpha0.1) plt.title(SARIMA Model - In-Sample Forecast vs Actuals) plt.xlabel(Date) plt.ylabel(AirPassengers) plt.legend() plt.show()接下来进行真正的样本外预测预测未来一段时间例如未来24个月的值。# 预测未来24个月2年 forecast_steps 24 # 获取预测值及置信区间对数尺度 forecast_object final_results.get_forecast(stepsforecast_steps) forecast_values_log forecast_object.predicted_mean forecast_ci_log forecast_object.conf_int() # 指数变换回到原始尺度 forecast_values np.exp(forecast_values_log) forecast_ci np.exp(forecast_ci_log) # 生成未来日期索引 last_date ts.index[-1] future_dates pd.date_range(startlast_date pd.DateOffset(months1), periodsforecast_steps, freqMS) # 绘制最终预测图 plt.figure(figsize(14, 8)) plt.plot(ts, labelHistorical Data, colorblue) plt.plot(future_dates, forecast_values, labelForecast, colorred, linestyle--) plt.fill_between(future_dates, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], colorr, alpha0.15, label95% Confidence Interval) plt.title(fSARIMA{best_order}x{best_seasonal_order} - {forecast_steps} Months Forecast) plt.xlabel(Date) plt.ylabel(AirPassengers) plt.legend() plt.grid(True, alpha0.3) plt.show() # 也可以打印出具体的预测数值 forecast_df pd.DataFrame({ Date: future_dates, Forecast: forecast_values.values, Lower CI: forecast_ci.iloc[:, 0].values, Upper CI: forecast_ci.iloc[:, 1].values }) print(forecast_df.head(10))观察预测图你应该能看到模型不仅延续了上升趋势还完美地复现了年度季节性波动。置信区间给出了预测的不确定性范围。6. 避坑指南与高级技巧在实际项目中严格按照上述流程走一遍你大概率能得到一个可用的SARIMA模型。但要想让模型更稳健、预测更准还需要注意以下这些我踩过的坑和总结的技巧。6.1 参数选择中的常见陷阱过度依赖自动定阶工具像pmdarima.auto_arima这样的工具非常方便但它给出的“最优”模型有时在业务上解释性不强或者因为陷入局部最优而忽略了更简单的模型。最佳实践是用自动工具给出建议再结合ACF/PACF图和自己对业务的理解进行微调。一个AIC值稍高但参数更简洁、更易解释的模型往往比一个复杂模型更受业务方欢迎。忽略残差诊断拟合完模型看一眼AIC不错就直接预测这是大忌。必须做残差诊断。如果残差不是白噪声说明还有信息未被模型提取预测就会有偏差。常见的残差问题包括残差自相关在ACF图上有显著滞后。解决方案尝试增加p或q或P,Q。残差异方差残差方差随时间增大。解决方案检查是否需要对数据做更复杂的变换如Box-Cox变换而不仅仅是取对数。残差非正态这不一定致命但可能影响预测区间的准确性。对于点预测影响不大但对于需要概率预测的场景需要注意。季节性周期s判断错误这是最致命的错误。如果周期设错整个季节性部分就全错了。务必通过业务知识如周、月、季、年和数据可视化绘制多年数据叠加图双重确认s的值。对于多重季节性如同时存在周周期和年周期标准SARIMA无法处理需要考虑其他模型如TBATS或Prophet。6.2 处理外生变量与节假日效应标准的SARIMA是单变量模型只利用自身的历史数据进行预测。但在现实中很多序列受外部因素影响。例如销售额受促销活动影响客流量受天气影响。这时可以使用SARIMAX模型X代表外生变量。statsmodels的SARIMAX类本身就支持外生变量。你需要准备一个与外生变量对应的DataFrame然后在fit方法中传入exog参数。# 假设我们有一个外生变量‘promotion’表示是否有促销0/1 exog_data pd.DataFrame(...) # 需要与ts_log具有相同的索引和长度 model_with_exog SARIMAX(ts_log, order(1,1,1), seasonal_order(1,1,1,12), exogexog_data) # 传入外生变量 results_with_exog model_with_exog.fit()预测时也需要提供未来期的外生变量值。对于节假日效应一种常见做法是将其作为虚拟变量哑变量加入外生变量。例如为春节所在的月份创建一个变量值为1其他月份为0。6.3 模型评估与滚动预测在正式使用模型预测未来之前必须评估其预测性能。不能只看在全部历史数据上的拟合效果这会导致过拟合乐观估计。正确的方法是使用时间序列交叉验证或滚动预测。思路是在历史数据上模拟一个不断向前滚动的预测过程。用前N个月的数据训练一个SARIMA模型。预测未来M个月。将预测值与真实值比较计算误差如均方根误差RMSE平均绝对百分比误差MAPE。将训练窗口向后移动M个月或1个月重复步骤1-3。将所有滚动窗口的误差平均得到模型在“未知”数据上的性能估计。from sklearn.metrics import mean_squared_error, mean_absolute_percentage_error def rolling_forecast(series, order, seasonal_order, train_size, test_size, s): 执行滚动预测评估 series: 时间序列对数变换后 train_size: 初始训练集大小 test_size: 每次预测的步长通常为季节周期s的倍数 history list(series[:train_size]) predictions [] actuals [] for t in range(train_size, len(series), test_size): # 用当前历史数据拟合模型 model SARIMAX(history, orderorder, seasonal_orderseasonal_order, enforce_stationarityFalse) model_fit model.fit(dispFalse) # 预测未来test_size步 forecast model_fit.forecast(stepstest_size) predictions.extend(forecast) actuals.extend(series[t:ttest_size].values) # 将真实值更新到历史中模拟时间推移 history.extend(series[t:ttest_size].values) # 计算误差注意预测值是对数需要转换回原始尺度再与原始尺度真实值比较 # 更合理的做法将对数尺度的预测值和真实值都指数变换回原始尺度再计算误差。 pred_original np.exp(predictions) actual_original np.exp(actuals) rmse np.sqrt(mean_squared_error(actual_original, pred_original)) mape mean_absolute_percentage_error(actual_original, pred_original) * 100 return pred_original, actual_original, rmse, mape # 示例使用最后3年36个月作为测试期每次滚动预测12个月 train_len len(ts_log) - 36 preds, actuals, rmse_val, mape_val rolling_forecast(ts_log, orderbest_order, seasonal_orderbest_seasonal_order, train_sizetrain_len, test_size12, # 每次预测一年 s12) print(f滚动预测RMSE: {rmse_val:.2f}) print(f滚动预测MAPE: {mape_val:.2f}%)MAPE平均绝对百分比误差是一个很好的业务指标比如MAPE5%意味着平均预测误差在5%左右。这个评估结果比样本内拟合的R²更有说服力。7. 在数模竞赛与业务中的实战要点最后结合我参加数模和做业务预测的经验分享几点SARIMA模型的实战心得。在数模竞赛中基础模型高级呈现SARIMA是一个很好的基准模型。在论文中你需要清晰地展示数据平稳化处理取对数、差分的过程图、ACF/PACF定阶分析图、模型诊断图、预测效果对比图。图表质量直接影响评分。结合其他模型不要只用一个SARIMA。可以将其与指数平滑ETS、线性回归考虑趋势和季节虚拟变量甚至简单的移动平均进行对比。用RMSE、MAPE等指标说明为什么SARIMA更优。强调假设检验在论文中写明你检查了残差的白噪声性Ljung-Box检验、正态性Jarque-Bera检验等证明模型是充分的。这体现了建模的严谨性。处理预测区间SARIMA可以给出预测的置信区间。在论文中画出这个区间并讨论其意义例如“我们有95%的把握认为未来销售额将落在X到Y之间”这比只给一个点预测值要高级得多。在业务分析中自动化与监控业务预测需要定期更新如每周/每月。可以编写自动化脚本定期用新数据重新训练模型并生成预测报告。同时监控模型的预测误差如果MAPE持续上升说明模型可能失效需要重新调整或更换。理解业务上限SARIMA是统计模型它基于历史模式外推。如果业务发生了根本性变化如新产品上市、政策巨变历史模式将不再适用模型的预测会失灵。这时需要结合业务判断或引入能够捕捉突变的外生变量。解释性优先在商业场景一个SARIMA(0,1,1)x(0,1,1,12)模型可能比一个SARIMA(2,1,2)x(1,1,1,12)模型更受欢迎即使后者AIC稍低因为前者更简洁更容易向非技术人员解释“模型认为下个月的值主要受上个月的误差和去年同月的值影响”。SARIMA模型是一个强大的工具但它不是银弹。它最适合于具有稳定趋势和固定季节性的中短期预测。对于长期预测不确定性会急剧增加对于模式快速变化的数据需要更灵活的模型。掌握它意味着你拥有了时间序列预测武器库中一件经典且可靠的重武器。从理解原理到代码实现再到避开实战中的各种坑这个过程本身就是数据分析能力的一次扎实提升。