1. 从“平稳”说起为什么时序预测要先“稳”下来聊到时间序列预测很多刚入门的朋友会一头扎进各种复杂的模型里比如LSTM、Prophet恨不得马上跑出个漂亮的预测曲线。但干了这么多年数据分析我越来越觉得磨刀不误砍柴工在动模型之前先把数据“看明白”比什么都重要。而“看明白”的第一步就是理解“平稳性”这个概念。你可以把它想象成你要在一片土地上盖房子如果这片地本身就在缓慢下沉或者高低不平非平稳那你无论在上面盖多漂亮的房子用多复杂的模型最终都可能因为地基问题而倒塌。平稳时间序列就是那块相对坚实、平整的地基。那么到底什么是平稳时间序列学术定义有点绕简单来说一个平稳的时间序列它的统计特性——主要是均值、方差和自协方差——不随时间推移而改变。这意味着序列没有明显的趋势比如持续上涨或下降也没有明显的季节性剧烈波动整个序列看起来像是在一个固定的水平线上下随机波动。这种“稳定性”是很多经典时序预测模型比如我们今天要重点讲的AR模型能够成立的前提假设。如果数据不平稳模型的基本假设就错了预测结果自然不可靠。为什么平稳性这么重要因为很多经典的统计预测方法其数学原理都建立在“未来与过去具有相同的统计规律”这个基础上。只有序列平稳了我们才能用历史数据中总结出的规律比如昨天和今天的关系去可靠地推断未来。否则你可能会用三年前的“牛市”规律去预测今天的“震荡市”结果可想而知。在实际项目中我见过太多因为忽略平稳性检验直接上ARIMA导致预测完全失真的案例。所以我的经验是拿到时序数据先别急着调参画个图做个检验看看它“稳不稳”。2. 平稳性的“体检报告”如何诊断你的时间序列知道了平稳性的重要性接下来就是实操怎么判断我的数据是不是平稳的这里我分享两个最常用、也最核心的方法肉眼观察图示法和量化检验统计检验。两者结合判断更准。2.1 图示法第一眼的直觉判断图示法是最快、最直观的方法。通常我们会看两个图时序图和自相关图。时序图就是把数据按时间顺序画出来。你盯着图看问自己几个问题数据是不是围绕着一个固定的水平线上下波动有没有明显的长期上升或下降趋势有没有周期固定、幅度变化的波动比如销售额每年“双十一”爆增如果答案是“围绕水平线随机波动”那初步看起来是平稳的。如果有一条明显的向上或向下的“线”或者有规律性的“波浪”那很可能就是非平稳的包含了趋势或季节性成分。自相关图则是更专业的工具。它展示的是时间序列在不同时间间隔滞后阶数上的自相关性。对于平稳序列它的自相关系数通常会快速衰减到零附近比如在滞后2、3期之后就接近0呈现出“截尾”或“拖尾但迅速衰减”的特征。而对于非平稳序列自相关系数会衰减得非常慢或者长期维持在较高水平因为历史值对当前值的影响持久而强烈。在Python里用statsmodels.graphics.tsaplots.plot_acf可以轻松画出这个图。注意图示法虽然直观但主观性强。尤其是当趋势或季节性不那么明显时不同的人可能会有不同的判断。所以它通常作为初步筛查我们还需要更客观的统计检验来下结论。2.2 统计检验ADF检验的原理与解读最常用的统计检验是Augmented Dickey-Fuller (ADF) 检验。它的原假设是时间序列存在单位根即序列是非平稳的。备择假设是序列不存在单位根是平稳的。这个检验会输出一个ADF统计量和对应的p值。我们决策的关键就看p值如果p值小于显著性水平通常取0.05我们就拒绝原假设认为序列是平稳的。如果p值大于0.05我们无法拒绝原假设倾向于认为序列是非平稳的。在Python中用statsmodels.tsa.stattools.adfuller函数可以一键完成检验。但这里有个关键点ADF检验的结果对检验中包含的项如常数项、趋势项很敏感。函数有几个参数需要注意regression: 这个参数决定了检验方程的形式。‘c’表示只包含常数项‘ct’表示包含常数项和线性趋势项‘ctt’表示包含常数项、线性及二次趋势项‘n’表示都不包含。如果你从时序图已经看出有明显趋势就该用‘ct’。如果不确定一个保守的做法是都试试或者参考一些自动确定阶数的包。autolag: 选择最佳滞后阶数的准则比如‘AIC’或‘BIC’。通常让函数自动选择即可。我常用的代码片段和解读逻辑是这样的from statsmodels.tsa.stattools import adfuller result adfuller(series, autolagAIC) # series是你的时间序列数据 print(ADF Statistic: %f % result[0]) print(p-value: %f % result[1]) print(Critical Values:) for key, value in result[4].items(): print(\t%s: %.3f % (key, value)) # 解读 if result[1] 0.05: print(序列是平稳的) else: print(序列是非平稳的)不仅要看p值也建议对比一下ADF统计量和不同置信度下的临界值。如果ADF统计量比临界值更负也是平稳的信号。2.3 处理非平稳数据差分与变换如果你的数据不幸被诊断为“非平稳”别慌这是常态。大多数真实世界的数据比如股票价格、月度销售额、气温记录都是非平稳的。我们有标准的处理方法让它“平稳化”。差分是最常用、最有效的方法。原理很简单计算当前时刻的值与前一时刻值的差值。一阶差分可以消除线性趋势二阶差分可以消除二次曲线趋势。在Python中用pandas的.diff()方法就能轻松实现。# 一阶差分 diff_1 series.diff().dropna() # 二阶差分 diff_2 series.diff().diff().dropna()差分之后一定要重新对差分后的序列做平稳性检验直到通过ADF检验为止。但也要注意过度差分会导致序列方差增大并可能引入不必要的相关性一般差分1-2次就够了。对数变换是另一种常用方法特别是当序列具有指数增长趋势或方差随时间增大的情况这称为异方差。取对数可以压缩数据的尺度使指数趋势变为线性趋势同时稳定方差。通常可以先取对数再对取对数后的序列进行差分这种方法在金融时间序列如股价分析中非常普遍。import numpy as np log_series np.log(series) # 然后再对 log_series 进行差分和平稳性检验3. 自回归模型用“昨天的自己”预测“明天的自己”当我们的时间序列变得平稳之后就可以请出今天的主角——自回归模型了。AR模型的全称是AutoRegressive Model它的核心思想非常直观且强大一个时间序列当前时刻的值可以用它过去若干个时刻值的线性组合再加上一个随机扰动白噪声来解释。用公式表示就是X_t c φ_1 * X_{t-1} φ_2 * X_{t-2} ... φ_p * X_{t-p} ε_t其中X_t是当前时刻的值。c是常数项。φ_1, φ_2, ..., φ_p是模型参数代表了过去各时刻对当前时刻的影响权重。p是模型的阶数意思是我们要用过去多少期的数据。ε_t是均值为0、方差恒定的白噪声代表那些无法用历史数据解释的随机波动。你可以把它想象成一个“自我对话”的过程。比如预测明天的气温AR模型认为明天的气温主要取决于今天、昨天、前天的气温分别乘以不同的权重再加上一些无法预料的随机变化比如突然来的冷空气。模型要学习的就是这些权重φ到底是多少。3.1 模型阶数p的选择AIC与BIC准则模型摆在那里第一个实际问题就是这个p我到底该选几用过去1期5期还是10期这就是模型定阶。选得太小比如p1模型可能太简单无法捕捉数据中完整的自相关结构这称为“欠拟合”。选得太大比如p20模型会变得非常复杂不仅计算量大还可能把噪声也当成规律学进来导致在新的数据上表现很差这称为“过拟合”。那怎么科学地选呢统计学给了我们两个非常实用的信息准则AIC和BIC。它们的思想都是在“模型拟合优度”和“模型复杂度”之间找一个平衡。AIC赤池信息准则。它鼓励数据拟合得更好但对模型复杂度的惩罚相对温和一些。BIC贝叶斯信息准则。它对模型复杂度的惩罚更严厉随着样本量增大惩罚力度会更强因此BIC倾向于选择更简单的模型。在实战中我通常的做法是设定一个候选阶数范围比如从0到15。用每个阶数p去拟合AR模型并计算对应的AIC和BIC值。画出AIC和BIC随p变化的折线图。通常AIC和BIC的值会随着p增加先下降后上升。选择使AIC或BIC值最小的那个p作为模型阶数。如果AIC和BIC选出的p不同在样本量不大时我可能更倾向于用BIC选的因为它更保守模型更简洁稳健。在Python的statsmodels库中拟合AR模型后可以直接获取AIC和BIC值。我们也可以通过循环来计算import statsmodels.api as sm # 假设 ts_data 是已经处理好的平稳时间序列 best_aic np.inf best_bic np.inf best_p_aic 0 best_p_bic 0 for p in range(0, 16): # 尝试0到15阶 try: model sm.tsa.AutoReg(ts_data, lagsp, old_namesFalse).fit() if model.aic best_aic: best_aic model.aic best_p_aic p if model.bic best_bic: best_bic model.bic best_p_bic p except: continue print(f根据AIC最佳阶数 p {best_p_aic}) print(f根据BIC最佳阶数 p {best_p_bic})3.2 模型拟合与参数估计OLS方法确定了阶数p下一步就是估计模型里的那些参数c, φ_1, ..., φ_p。对于AR模型最常用的方法是普通最小二乘法。它的目标很直观找到一组参数使得模型预测值c φ_1 * X_{t-1} ... φ_p * X_{t-p}与实际观测值X_t之间的差距的平方和最小。statsmodels库的AutoReg类在调用.fit()方法时默认使用的就是OLS。拟合完成后我们可以通过model.params查看所有估计出的参数值。这些参数有明确的统计意义例如φ_1显著为正且接近1说明上一期对本期有很强的正向影响序列的惯性很大如果φ_1为负则可能意味着一种均值回复的特性。3.3 模型检验残差分析的重要性模型参数估计好了是不是就能直接用来预测了别急还有关键一步模型检验。一个合格的AR模型要求其残差序列ε_t是一个白噪声序列。什么是白噪声就是均值为零、方差恒定、且各时刻之间完全无关的随机序列。如果残差不是白噪声说明模型还没有把数据中的规律提取干净还有信息藏在残差里那么这个模型就是不充分的。怎么检验残差是不是白噪声呢主要看两点自相关性检验最常用的是Ljung-Box检验。它的原假设是残差序列在检验的滞后阶数内没有显著的自相关性。我们期望看到的结果是p值大于0.05这样我们就无法拒绝原假设认为残差是白噪声。正态性检验虽然AR模型本身不要求残差严格服从正态分布但如果残差近似正态模型的很多统计推断会更稳健。可以用QQ图或Shapiro-Wilk检验来观察。在Python中可以方便地进行这些检验from statsmodels.stats.diagnostic import acorr_ljungbox import scipy.stats as stats # 获取模型残差 residuals model.resid # Ljung-Box检验检验前10阶的自相关 lb_test acorr_ljungbox(residuals, lags10, return_dfTrue) print(lb_test) # 我们希望所有滞后阶数的p值都大于0.05 # 正态性检验Shapiro-Wilk适用于小样本 shapiro_test stats.shapiro(residuals) print(fShapiro-Wilk test statistic: {shapiro_test[0]}, p-value: {shapiro_test[1]}) # p值大于0.05则不能拒绝正态性原假设只有通过了残差白噪声检验我们才能比较有信心地说这个AR模型已经较好地捕捉了数据中的线性自相关结构可以用于预测了。4. Python全流程实战从数据到预测理论说了这么多现在我们用一个完整的例子串起来。假设我们有一组某产品过去100天的日销量数据我们要用AR模型预测未来5天的销量。4.1 数据准备与平稳性处理首先我们生成一组模拟数据并引入一个轻微的趋势来模拟真实情况。import numpy as np import pandas as pd import matplotlib.pyplot as plt import statsmodels.api as sm from statsmodels.tsa.stattools import adfuller from statsmodels.graphics.tsaplots import plot_acf, plot_pacf # 设置随机种子保证可复现 np.random.seed(42) # 生成模拟数据基础水平 线性趋势 自回归成分 噪声 n 100 time np.arange(n) trend 0.05 * time # 轻微线性趋势 # 生成一个AR(2)过程的基础序列 ar_component np.zeros(n) ar_component[0] np.random.randn() ar_component[1] np.random.randn() for t in range(2, n): ar_component[t] 0.6 * ar_component[t-1] - 0.2 * ar_component[t-2] np.random.randn() * 0.5 # 组合成最终序列 original_series 10 trend ar_component # 转换为pandas Series并创建时间索引 dates pd.date_range(start2023-01-01, periodsn, freqD) ts pd.Series(original_series, indexdates) # 绘制原始序列 plt.figure(figsize(12, 6)) plt.plot(ts, labelOriginal Series) plt.title(Original Time Series (with Trend)) plt.xlabel(Date) plt.ylabel(Value) plt.legend() plt.grid(True) plt.show()从图上我们能看出一个明显的上升趋势这显然是非平稳的。接下来我们进行ADF检验确认并进行差分处理。# 1. 对原始序列进行ADF检验 result_original adfuller(ts, autolagAIC) print(f原始序列 ADF p-value: {result_original[1]:.4f}) if result_original[1] 0.05: print(- 原始序列是非平稳的需要进行差分。) # 2. 进行一阶差分 ts_diff ts.diff().dropna() # 绘制差分后序列 plt.figure(figsize(12, 6)) plt.plot(ts_diff, labelDifferenced Series (1st Order)) plt.title(Time Series After First-Order Differencing) plt.xlabel(Date) plt.ylabel(Differenced Value) plt.legend() plt.grid(True) plt.show() # 3. 对差分后序列进行ADF检验 result_diff adfuller(ts_diff, autolagAIC) print(f一阶差分后序列 ADF p-value: {result_diff[1]:.4f}) if result_diff[1] 0.05: print(- 一阶差分后序列是平稳的。) else: print(- 一阶差分后序列仍非平稳可能需要二阶差分。)在这个模拟例子中一阶差分后序列应该能通过平稳性检验。我们后续的建模将基于这个平稳的差分序列ts_diff进行。但请注意最终预测结果需要“还原”到原始尺度。4.2 模型识别、定阶与拟合序列平稳后我们通过自相关图和偏自相关图来初步判断AR模型的阶数并用AIC/BIC准则精确确定。# 绘制自相关图(ACF)和偏自相关图(PACF) fig, axes plt.subplots(1, 2, figsize(15, 4)) plot_acf(ts_diff, lags20, axaxes[0]) plot_pacf(ts_diff, lags20, axaxes[1], methodywm) # 使用ywm方法计算PACF axes[0].set_title(Autocorrelation Function (ACF)) axes[1].set_title(Partial Autocorrelation Function (PACF)) plt.show()对于AR模型我们主要看偏自相关图。如果PACF在滞后p阶后突然截尾即之后的系数在置信区间内不显著那么p可能就是AR模型的阶数。从图上我们可以观察截尾点。然后我们用AIC/BIC准则来精确选择p。# 使用AIC/BIC准则自动选择最佳阶数p (在差分后的平稳序列上) best_aic np.inf best_bic np.inf best_p_aic 0 best_p_bic 0 max_lag 15 # 最大尝试阶数 for p in range(0, max_lag1): try: # 使用AutoReg注意这里拟合的是平稳的差分序列 ts_diff model_temp sm.tsa.AutoReg(ts_diff, lagsp, old_namesFalse).fit() if model_temp.aic best_aic: best_aic model_temp.aic best_p_aic p if model_temp.bic best_bic: best_bic model_temp.bic best_p_bic p except Exception as e: print(f拟合AR({p})时出错: {e}) continue print(fAIC推荐的最佳阶数 p {best_p_aic} (AIC{best_aic:.2f})) print(fBIC推荐的最佳阶数 p {best_p_bic} (BIC{best_bic:.2f})) # 通常选择BIC推荐的阶数因为它更倾向于简洁模型防止过拟合 selected_p best_p_bic print(f\n最终选定模型阶数 p {selected_p})假设我们选定了p2。接下来我们用这个阶数来正式拟合模型。# 使用选定的阶数p拟合AR模型 model_ar sm.tsa.AutoReg(ts_diff, lagsselected_p, old_namesFalse).fit() print(model_ar.summary())summary()会输出非常详细的报告包括每个参数的估计值、标准误、t统计量和p值。我们需要关注参数的显著性p值小于0.05通常认为显著以及模型的整体拟合优度R-squared等。4.3 模型诊断与预测拟合好模型后必须进行残差诊断。# 获取残差 residuals model_ar.resid # 1. 绘制残差序列图 plt.figure(figsize(12, 6)) plt.plot(residuals) plt.axhline(y0, colorr, linestyle--) plt.title(Residuals of AR Model) plt.xlabel(Date) plt.ylabel(Residual) plt.grid(True) plt.show() # 2. 残差自相关检验 (Ljung-Box) from statsmodels.stats.diagnostic import acorr_ljungbox lb_test acorr_ljungbox(residuals, lags[10], return_dfTrue) # 检验前10阶 print(Ljung-Box Test for Residuals:) print(lb_test) # 我们希望p值 0.05 # 3. 绘制残差分布直方图和QQ图 fig, axes plt.subplots(1, 2, figsize(12, 4)) axes[0].hist(residuals, bins20, edgecolorblack, alpha0.7) axes[0].set_title(Histogram of Residuals) axes[0].set_xlabel(Residual) axes[0].set_ylabel(Frequency) import scipy.stats as stats stats.probplot(residuals, distnorm, plotaxes[1]) axes[1].set_title(Q-Q Plot of Residuals) plt.tight_layout() plt.show()如果残差图看起来是随机的Ljung-Box检验p值大于0.05Q-Q图上的点大致在一条直线上那么模型诊断通过。最后我们进行预测。这里有个关键点我们是在平稳的差分序列ts_diff上建模的所以预测结果也是差分序列的未来值。我们需要将其还原到原始序列的尺度上。# 预测未来5步在差分序列上 forecast_steps 5 forecast_diff model_ar.forecast(stepsforecast_steps) print(f预测的未来 {forecast_steps} 期差分值:\n{forecast_diff}) # 将差分预测值还原为原始序列的预测值 # 因为 ts_diff ts[t] - ts[t-1]所以 ts[t] ts[t-1] ts_diff[t] # 我们需要原始序列的最后一个值作为起点 last_original_value ts.iloc[-1] forecast_original [] current_value last_original_value for diff_val in forecast_diff: next_value current_value diff_val forecast_original.append(next_value) current_value next_value forecast_original pd.Series(forecast_original, indexpd.date_range(startts.index[-1] pd.Timedelta(days1), periodsforecast_steps, freqD)) print(f\n还原后的未来 {forecast_steps} 期原始序列预测值:\n{forecast_original}) # 可视化 plt.figure(figsize(12, 6)) plt.plot(ts, labelHistorical Data, colorblue) plt.plot(forecast_original, labelAR Model Forecast, colorred, markero) plt.fill_between(forecast_original.index, forecast_original - 1.96 * np.std(residuals), # 近似95%置信区间 forecast_original 1.96 * np.std(residuals), colorred, alpha0.2, label95% Confidence Interval) plt.title(AR Model Forecast vs Historical Data) plt.xlabel(Date) plt.ylabel(Value) plt.legend() plt.grid(True) plt.show()这个还原过程是ARIMA模型中“积分”部分的逆运算。我们得到了未来5天的点预测值并且用残差的标准差构造了一个简单的预测区间这能让我们对预测的不确定性有个直观的认识。5. AR模型的局限与实战避坑指南AR模型是时序预测的基石但它并非万能。理解它的局限性能帮助你在正确的地方使用它并避免很多常见的坑。局限性仅适用于线性关系AR模型本质是线性模型。它假设当前值与过去值呈线性关系。如果真实数据中存在复杂的非线性依赖如周期性突变、状态切换AR模型可能力不从心。对异常值敏感由于基于最小二乘法异常值会显著影响参数估计导致模型偏离。只能捕捉自身历史信息AR模型是“自”回归只利用序列自身的历史信息。如果预测目标强烈依赖于其他外部变量如促销活动、天气对销量的影响纯AR模型会遗漏这些重要信息此时需要考虑带外生变量的ARX模型或向量自回归模型。要求序列平稳这是最重要的前提。对非平稳序列直接使用AR模型预测结果往往是无意义的。实战避坑指南平稳性检验不是一次性的对于长期预测项目数据特征可能随时间变化。建议定期如每月或每季度重新检验序列的平稳性必要时重新差分或变换。差分不是越多越好我见过有人为了追求ADF检验的p值足够小连续做三阶甚至四阶差分。这通常会导致序列失去经济或物理意义并放大噪声。一般一阶或二阶差分足矣。差分后序列的均值应在零附近波动。小心“伪回归”如果你用AR模型去拟合一个纯随机游走非平稳序列可能会得到一个拟合优度R²很高的模型但这完全是虚假的。这就是为什么必须先做平稳性检验。样本量要充足AR模型需要估计p1个参数p个自回归系数加一个常数。经验法则是样本量至少是参数数量的10-20倍。如果只有100个数据点却去拟合一个AR(10)模型结果很可能不稳定。预测区间会迅速变宽AR模型做多步预测时是将上一步的预测值作为下一步的输入。这种迭代方式会导致预测误差不断累积和放大。因此AR模型通常只适合做短期预测比如未来1-5期。从上面代码绘制的预测区间图也能看出越往后的预测不确定性区间越宽。模型需要更新世界在变数据的生成过程也可能在变。一个基于去年数据训练的AR模型今年可能就不适用了。对于在线预测系统需要设计模型重训练或在线学习的机制。在我自己的项目中AR模型常常作为一个基准模型。我会先用它跑出一个结果然后再尝试更复杂的模型如ARIMA、SARIMA、机器学习模型。如果复杂模型的提升不大那么简洁明了的AR模型往往是更优的选择因为它更容易解释和维护。记住在时间序列预测里往往不是模型越复杂越好而是越合适越好。从平稳性检验到AR模型这套经典流程为你提供了一个坚实可靠的起点。