ARIMAX多变量时间序列预测实战：数据预处理到滚动回测

简介：基于ARIMAX的多变量预测模型Python源码与配套数据集，面向统计学、数据科学及相关工科专业的毕业设计、课程设计和期末大作业场景，适合需要快速上手多变量时序预测项目、又担心代码无法运行的学习者。资源包共8个文件，包含2个Python脚本（数据预处理、模型构建与预测）、2个CSV数据集、2张可视化结果图、1份README说明文档及.gitignore配置，整体压缩包仅148KB，结构精炼清晰，便于按流程逐步复现。项目为经导师指导并通过的高分设计，代码完整可运行，能帮助读者理解ARIMAX模型处理外生变量的思路、数据归一化与预测对比的完整流程；自带数据与结果图可直接验证效果，也可替换自己的数据集进行扩展实验。资源已有55人学习/下载，对希望用最短时间完成高质量课程项目并展示模型能力的同学有较强参考价值。

1. 多变量预测模型为什么值得用 ARIMAX 而不是 ARIMA

ARIMAX 这个词乍听起来像是 ARIMA 加了个外生变量就完事，但真正落地时坑全在细节里。当你要预测的目标不仅受自身历史影响，还明显受外部因素——比如促销、天气、价格变动——驱动时，普通的 ARIMA 会丢信息。这份源码拆开来看，datapre.py 负责把多列数据对齐、做平稳性检验，arimax.py 承担定阶、拟合、预测和残差诊断，data.csv 是原始数据，datacf.csv 从命名看就是差分后的序列或另一份对照数据。它适合两类人：一是正在做毕业设计或期末大作业、需要把完整流程跑通并讲清楚每一步的同学，二是平时只用过单变量 ARIMA、想看看外生变量建模边界在哪的工程师。用流量预估或销量预测这类场景来理解，模型给出的不是单序列外推，而是把外部变量压缩成线性权重再叠加到自回归项上。

2. 数据准备：datapre.py 里的对齐、差分与平稳性检验

2.1 先把 data.csv 读成时间索引加数值列

拿到压缩包后第一件事不是直接跑 arimax.py，而是先看数据长什么样。data.csv 里通常是一列时间戳加一列目标值 y，再加若干列外生变量。读入的时候要注意，时间列必须解析成 pandas 的 DatetimeIndex，否则后面SARIMAX的exog对齐会非常麻烦。

import pandas as pd df = pd.read_csv("data.csv", parse_dates=True, index_col=0) print(df.shape) print(df.dtypes) print(df.isnull().sum())

这里parse_dates=True让 pandas 尝试把第一列解析成时间，index_col=0把时间列设置为索引。打印出来的三个信息分别回答：样本量多少、每列类型是什么、有没有缺失值。时序数据里缺失值不能像普通表格那样直接删行，常见做法是用ffill()前向填充，或者interpolate(method="linear")线性插值。如果某列缺失值超过 30%，这一列基本可以放弃，填充出来的长段常数值对模型系数估计是干扰。

datacf.csv 从命名看大概率是差分后的数据或者备用的预测对比数据。跑完 datapre.py 之后可以读出来和原数据对一下行数，确认预处理到底做了哪一步。我一般会先只输出两个文件的 shape，再决定后续用哪个作为建模输入。

2.2 用 ADF 检验确定差分阶数 d

ARIMAX 和 ARIMA 一样，要求目标序列至少是趋势平稳的。判断平稳性最直接的做法是画图，但图上看着平稳不代表统计上通过了检验。ADF 检验是这里最常用的工具，原假设是序列存在单位根，也就是不平稳。

from statsmodels.tsa.stattools import adfuller def adf_test(series, name): res = adfuller(series.dropna(), autolag="AIC") print(f"{name}: ADF={res[0]:.3f}, p={res[1]:.4f}") return res[1] < 0.05 print("y", adf_test(df["y"], "y")) for c in ["x1", "x2", "promo"]: if c in df.columns: print("exog", adf_test(df[c], c))

autolag="AIC"表示自动选择滞后阶数来消除自相关，比固定 lag 更稳。返回值里第一个是 ADF 统计量，第二个是 p 值。p 小于 0.05 时拒绝单位根假设，认为序列平稳；p 大于 0.05 时做一阶差分，然后重新检验。

实际项目中 d = 2 非常少见。如果一阶差分后还不平稳，多数情况是序列里有结构突变或者异常点，直接做二阶差分会把有效波动一起抹掉。拆这个项目的时候我在 datapre.py 里看到它对 y 和 x1 都做了 ADF 判断，这就是为什么模型在测试集上比直接对原始序列建 ARIMA 更稳。

2.3 外生变量的类型转换、对齐与差分口径

外生变量必须是数值型。如果 data.csv 里有促销标识这类 0/1 变量，直接读进来是整型没问题；但如果某列是金额含千分位逗号，pandas 会读成字符串，必须强制转换。

x_cols = ["x1", "x2", "promo"] for c in x_cols: df[c] = pd.to_numeric(df[c], errors="coerce").ffill()

errors="coerce"的意思是把无法解析的值变成 NaN，再配合ffill()用上一个有效值补齐。这里有一个我拆代码时特别注意的点：如果目标序列 y 需要 d=1 差分，外生变量也要考虑差分口径。statsmodels 的 SARIMAX 对endog做差分，但exog并不会自动做同样变换。也就是说，当你给order=(p, 1, q)时，模型做的事是回归Δy_t对exog_t，如果 x 本身带明显趋势，系数会被趋势主导。

d = 1 # 来自上一节的 ADF 检验 if d > 0: df["y_model"] = df["y"].diff(d).dropna() exog_model = df[x_cols].diff(d).dropna() else: df["y_model"] = df["y"] exog_model = df[x_cols]

这里对 y 做了 d 阶差分，同时对 x 列做了同阶差分，保证进入模型的外生变量和endog在一个量纲体系下。差分之后行数会少 d 行，所以后面训练测试切分要基于差分后的 DataFrame，而不是原始的 df。很多初学者在这里踩坑：差分后没有 dropna，结果 SARIMAX 拟合时把 NaN 带进去直接报错。

外生变量还有一个滞后对齐的问题。如果你要预测的是 t 期的 y，而 x 在 t 期是已知的（比如当天的计划促销活动），直接用当期 x 没问题。如果 x 只能拿到 t 期的上一期值，那就必须df["x1"] = df["x1"].shift(1)。这一步决定了预测口径，也直接决定后面get_forecast时外生变量怎么填。

3. arimax.py 实现解析：SARIMAX 封装、定阶与残差诊断

3.1 为什么用 statsmodels 的 SARIMAX 来跑 ARIMAX

ARIMAX 的数学形式可以写成：

y_t = c + Σ φ_i y_{t-i} + Σ θ_j ε_{t-j} + Σ β_k X_{k,t} + ε_t

区别在于多了一项外生变量 X 及其系数 β。一种常见错误做法是先用最小二乘回归把 y 对 X 回归，拿残差再建 ARIMA。这种做法把回归和时序估计拆成了两步，理论上会有偏差，因为第一步没有考虑残差里的自相关结构，β 的估计效率低。statsmodels 的SARIMAX是状态空间实现，把外生回归项和 ARIMA 误差项放在同一个似然函数里联合估计，这才是严格意义上的 ARIMAX。

从工程角度说，SARIMAX也省事：它能直接处理order，支持exog，输出里有 AIC、AICc 和残差序列，不用自己拼装。所以 arimax.py 里有很大概率直接 import 了SARIMAX，而不是再用普通 ARIMA 加手工回归。

3.2 从创建模型到预测的最小闭环

from statsmodels.tsa.statespace.sarimax import SARIMAX train = df_model.iloc[:-12] test = df_model.iloc[-12:] model = SARIMAX( endog=train["y_model"], exog=train[x_cols], order=(2, 1, 2), enforce_stationarity=False, enforce_invertibility=False, ) fit = model.fit(disp=False) forecast = fit.get_forecast(steps=12, exog=test[x_cols]) yhat = forecast.predicted_mean print(yhat.head())

endog是差分后的目标序列，exog是同步对齐好的外生变量矩阵，order=(2, 1, 2)里的三个数分别对应 AR 阶数 p、差分阶数 d、MA 阶数 q。这段代码里 p=2、d=1、q=2 只是一个初始值，真正上线时要靠网格搜索定。enforce_stationarity=False和enforce_invertibility=False是放开约束，避免在边界参数上直接抛异常，代价是可能拟合出非平稳的 AR 项，所以只能用于搜索阶段，最终模型要拿这两个变量保持默认值再验证一遍。

get_forecast里的exog必须覆盖预测区间的全部外生变量，行数和steps一致，列名和训练时一致。这里最容易出的错是训练集里 x 列有 3 列，预测时只给了 2 列，statsmodels 会报 KeyError。预测结果yhat的索引和test对齐，可以直接拿去算误差。

3.3 定阶：ACF/PACF 初筛加 AIC 网格搜索

理论上可以靠 ACF 和 PACF 图判断 p、q：ACF 拖尾、PACF 截尾看成 AR 项，反过来看成 MA 项。但实际数据很少给出干净的截尾图，尤其是多变量模型里外生变量还会吸收一部分自相关信息。所以我更习惯用网格搜索把 p、q 从 0 到 3 全扫一遍，用 AIC 和 AICc 做选择。

import itertools import warnings from statsmodels.tools.sm_exceptions import ConvergenceWarning warnings.filterwarnings("ignore", category=ConvergenceWarning) results = [] for p, q in itertools.product(range(4), range(4)): try: trial = SARIMAX( train["y_model"], exog=train[x_cols], order=(p, 1, q), enforce_stationarity=False, enforce_invertibility=False, ).fit(disp=False) results.append((trial.aic, trial.aicc, p, q)) except Exception: continue results.sort() for aic, aicc, p, q in results[:5]: print(f"p={p}, q={q}, AIC={aic:.2f}, AICc={aicc:.2f}")

range(4)意味着 p、q 都从 0 扫到 3，一共 16 个组合。代码里用 try-except 包住拟合适配器，是因为部分阶数组合（比如 p=0、q=0 加外生变量）可能收敛失败，直接跳过比中断整个搜索更好。排序时用 AIC 也可以，但样本量小于 100 时 AICc 更可靠，因为 AICc 在小样本下对参数个数有更强的惩罚。

注意 AIC 只能在同一份训练集上横向比较，不能拿这个项目的 AIC 去跟另一个数据集上的模型比。我拆 arimax.py 时看到它把搜索结果写成了 CSV 存下来，再选 AICc 最小的组合，这个习惯值得保留，后续换数据时能直接翻历史记录。

3.4 残差诊断：Ljung-Box 白噪声检验

选好阶数后必须做残差诊断，否则模型的良好拟合可能只是把噪声也学进去了。最核心的检验是 Ljung-Box，原假设是残差序列不存在自相关。

from statsmodels.stats.diagnostic import acorr_ljungbox resid = fit.resid.dropna() lb_test = acorr_ljungbox(resid, lags=[10], return_df=True) print(lb_test)

lags=[10]表示检验前 10 阶滞后的联合自相关，return_df=True让结果以 DataFrame 形式输出，方便直接读列名。输出的 p 值如果大于 0.05，说明残差近似白噪声，模型已经把自相关信息提取干净了。如果 p 值小于 0.05，先检查两件事：一是差分是否真的消除了趋势，二是外生变量里是否漏掉了某个重要的周期变量，比如星期几。这两个问题用 aumentar p、q 是修不好的。

4. 多变量预测的评估：时间切分、误差指标与差分还原

4.1 按时间顺序切分训练集和测试集

时序预测的评估切分和普通机器学习完全不同。不能用train_test_split的默认随机抽样，那会把未来数据混进训练集，造成评估结果虚高。正确做法是固定最后 N 期作为测试集，N 最好覆盖一个完整业务周期，比如有周规律的数据至少留 7 天，有月规律的数据至少留 30 天。

cut = int(len(df_model) * 0.8) train = df_model.iloc[:cut] test = df_model.iloc[cut:] print(f"train: {train.shape[0]} rows, test: {test.shape[0]} rows")

这里用 80% 训练、20% 测试的比例。样本量本身就不大时，我一般会固定最后 12 到 24 期做测试，而不是按比例切。原因很简单：比例切分在短序列上可能让训练集只有三四十行，SARIMAX 在这种样本量下估计出来的外生变量系数方差很大，说服力不够。

4.2 RMSE、MAE、MAPE 怎么配合着看

模型预测完必须计算误差指标。这套项目里最常用的三个指标分别是 RMSE、MAE 和 MAPE。RMSE 对大的偏差敏感，能暴露某些极端点预测失败的情况；MAE 更稳健；MAPE 是相对误差，方便跟其他模型直接对比。

from sklearn.metrics import mean_absolute_error, mean_squared_error import numpy as np mae = mean_absolute_error(test["y_model"], yhat) rmse = np.sqrt(mean_squared_error(test["y_model"], yhat)) mape = np.mean(np.abs((test["y_model"] - yhat) / test["y_model"])) * 100 print(f"MAE = {mae:.3f}") print(f"RMSE = {rmse:.3f}") print(f"MAPE = {mape:.2f}%")

这里直接使用差分后的y_model计算指标，等同于在增量层面评估预测。RMSE 和 MAE 的单位与目标列一致，MAPE 是无量纲的。注意如果y_model里有接近 0 的值，MAPE 会被极小的分母放大到几百甚至几千，这时候要改用 MAE 做对比。

4.3 差分还原与 datacf.csv 对照

如果建模用的是差分后的序列，预测结果也要还原到原始量纲。还原逻辑不复杂：差分序列的预测值逐项累加，再加上差分前的最后一个真实值。

yhat_level = yhat.cumsum() + train["y"].iloc[-1] compare = pd.DataFrame({ "true": test["y"].iloc[: len(yhat_level)], "pred": yhat_level, }) compare.to_csv("predict_result.csv", encoding="utf-8-sig") print(compare.head())

如果你在前面建模时用了df["y_model"] = df["y"].diff(d)，那么预测值是步长上的增量。cumsum()把增量累加起来，train["y"].iloc[-1]是差分起点之前的最后一个真实值，两者相加才是原始尺度上的预测结果。这里最容易犯的错是拿差分后的yhat直接和原始test["y"]对齐，误差看起来会非常小，因为增量波动往往远小于绝对水平，评估结果失真。

datacf.csv 在这一步就有用了。它是一个备用的对照序列，读出来后把三列放一起：test["y"]、yhat_level、datacf的某一列。如果预测曲线和 datacf 的走势基本重合，说明预处理和模型口径是自洽的；如果方向相反，回查是否差分被做了两次，或者外生变量没有对齐。

4.4 滚动回测：单次切分之外的验证方式

单次切分只能说明模型在这一个时间窗口上表现好，换一段历史窗口未必稳定。滚动回测的做法是：第一次用前 80% 训练、预测下一个点，然后把该点真实值并入训练集，再预测再往后一个点，重复直到覆盖整个测试期。

history = df_model.iloc[:cut].copy() rolling_preds = [] for i in range(cut, len(df_model)): row = df_model.iloc[[i]] trial = SARIMAX( history["y_model"], exog=history[x_cols], order=(1, 1, 2), enforce_stationarity=False, enforce_invertibility=False, ).fit(disp=False) next_x = row[x_cols] fc = trial.get_forecast(steps=1, exog=next_x) rolling_preds.append(fc.predicted_mean.iloc[0]) history = pd.concat([history, row])

这个循环里每次只预测一步，然后把真实值加进history重新拟合。SARIMAX的重新拟合开销不小，所以只推荐在数据量几百行以内时使用。滚动回测得到的误差比单次切分更接近真实上线表现，因为它模拟了「每个预测日只能拿到过去数据」的场景。

5. 小样本数据下 ARIMAX 的调参边界与实际使用技巧

5.1 外生变量的未来已知性决定预测步长

get_forecast要求预测期内的外生变量全部已知。如果 x 本身需要另一个模型先预测出来，误差会层层传递，ARIMAX 的优势就被抵消了。所以选 x 时优先选可计划变量：营销活动排期、节假日、天气预报这类能提前拿到的数据，而不是价格、竞品销量这类需要再预测的变量。

5.2 样本量小于 100 时不要只看 AIC

AIC 是渐近指标，样本量不大时更倾向于选择较复杂的模型。statsmodels 的拟合结果里提供了aicc，这是小样本修正版，选阶时以它为准，甚至可以拿最后 5 到 10 期做 holdout 验证。把 AIC 排名前三个候选阶分别拟合，看它们在 holdout 上的 RMSE 是否与 AIC 排名一致。如果不一致，以 holdout 表现为准。

5.3 外生变量共线性与系数符号检查

多变量模型里外生变量之间如果高度相关，β 的估计会不稳定。在建完模型后，先算外生变量相关矩阵，再检查各个系数的符号是否符合业务直觉。促销变量的系数应该是正数，如果变成负数，说明它跟另一个外生变量有共线性，考虑删掉其中一个。单看拟合优度很容易忽略这个问题，但系数方向错了的模型上线后业务方是不会接受的。

print(df[x_cols].corr()) for i, c in enumerate(x_cols): print(c, "=", round(fit.params[i + 1], 4))

fit.params里第一个元素通常是常数项，后面按exog列的顺序对应系数。相关系数超过 0.8 的变量对只保留一个，或者对其中一个做差分后再进入模型。模型输出里的 z 值和 p 值也可以辅助判断，系数不显著的外生变量留着只会增加预测期的数据采集成本。把这几个检查点过一遍，再用滚动回测做最终验证，模型就能从「跑通」往「真实可用」再走一步。

本文还有配套的精品资源，点击获取