我：“XGBoost时序好强！”导师：“那还要ARIMA干什么？”我：“一起用？”导师：“先搞懂XGBoost为什么要学ARIMA的残差！”

哈喽，大家好~

最近，有同学聊到了，ARIMA已经能预测时间序列了，为什么还要再加一个XGBoost？

很常见的一个问题，其实两者擅长的东西不一样：ARIMA负责抓住趋势和自相关，XGBoost负责学习ARIMA没解释掉的复杂波动。

今天我们就用一个完整案例，把这套融合方法讲清楚~

融合逻辑

假设我们要预测每天的商品销量~

销量通常包含几类信息：

• 整体趋势，比如随着时间不断增长；
• 周期性，比如周末销量更高；
• 自相关，比如昨天销量高，今天通常也不会突然变得很低；
• 复杂波动，比如促销、天气、节假日带来的非线性影响。

ARIMA比较擅长前三类中的“趋势”和“时间依赖关系”。它会根据历史值以及历史误差，预测下一个时间点。

但ARIMA本质上是线性模型，它很难准确表示：

• 某个变量达到阈值后产生的突变；
• 多个因素组合产生的非线性关系；
• 复杂的局部波动。

所以我们可以分两步：

第一步，让ARIMA先给出一个基础预测：

第二步，计算真实值和ARIMA预测值之间的残差：

然后让XGBoost学习这个残差：

最终预测结果为：

这里的意思很简单：ARIMA先完成“基础分”，XGBoost再负责“改卷子”。

一个小例子

假设某天真实销量是100。

ARIMA根据趋势和历史规律，预测销量为95，那么残差就是：

如果我们发现最近几天的残差分别是：

2、4、5、6

说明ARIMA最近一直低估销量，而且误差还有逐渐增大的趋势。

XGBoost看到这些历史残差后，可能预测下一次残差为4，于是最终结果就是：

所以，XGBoost不是直接替代ARIMA，而是在学习：ARIMA通常会在哪些情况下预测偏低，哪些情况下预测偏高。

一句话概括的话，就是ARIMA负责解释“规律”，XGBoost负责修正“规律之外的部分”。

案例数据和整体流程

下面我们构造一份日销量数据，数据中包含：

• 线性增长趋势；
• 周期波动；
• 非线性波动；
• 随机噪声；
• 一个局部冲击。

整个流程是：

1. 按时间顺序切分训练集和测试集；
2. 用训练集拟合ARIMA；
3. 得到ARIMA在训练集上的残差；
4. 用残差的历史滞后项和时间特征训练XGBoost；
5. 先用ARIMA预测测试集；
6. XGBoost递归预测残差；
7. 将两部分相加得到最终预测。

这里必须注意，时间序列不能随机打乱切分，否则会把未来信息泄露给过去。

完整代码实现

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import acf
from sklearn.metrics import mean_absolute_error, mean_squared_error
from xgboost import XGBRegressor

np.random.seed(42)

# 1. 时间序列数据
n = 720
t = np.arange(n)

trend = 0.04 * t
weekly = 8 * np.sin(2 * np.pi * t / 7)
long_cycle = 5 * np.sin(2 * np.pi * t / 45)

# 非线性波动：超过某个周期位置后，波动幅度发生变化
nonlinear = 0.015 * (t % 30) ** 2

noise = np.random.normal(0, 2.5, n)

# 局部冲击
shock = np.zeros(n)
shock[500:530] += 15

y = 100 + trend + weekly + long_cycle + nonlinear + noise + shock
dates = pd.date_range("2022-01-01", periods=n, freq="D")

data = pd.DataFrame({"y": y}, index=dates)

# 2. 按时间切分数据
train_size = 560
train = data.iloc[:train_size]
test = data.iloc[train_size:]

y_train = train["y"].values
y_test = test["y"].values

# 3. ARIMA建模
arima_model = ARIMA(y_train, order=(2, 1, 2))
arima_fit = arima_model.fit()

# 训练集拟合值与残差
arima_train_fit = np.asarray(arima_fit.fittedvalues)
train_residual = y_train - arima_train_fit

# 预测测试集的基础结果
arima_test_pred = np.asarray(
    arima_fit.get_forecast(steps=len(y_test)).predicted_mean
)

# 4. 构造XGBoost训练特征
def time_features(index):
    """
    使用时间位置构造周期特征。
    这里没有使用未来销量，只使用已知的时间信息。
    """
    return np.column_stack([
        np.sin(2 * np.pi * index / 7),
        np.cos(2 * np.pi * index / 7),
        np.sin(2 * np.pi * index / 30),
        np.cos(2 * np.pi * index / 30)
    ])

lag = 7
X_train, y_res_train = [], []

for i in range(lag, len(train_residual)):
    residual_lags = train_residual[i-lag:i][::-1]
    calendar = time_features(np.array([i]))[0]
    X_train.append(np.r_[residual_lags, calendar])
    y_res_train.append(train_residual[i])

X_train = np.array(X_train)
y_res_train = np.array(y_res_train)

# 5. 训练XGBoost残差模型
xgb_model = XGBRegressor(
    n_estimators=400,
    max_depth=4,
    learning_rate=0.03,
    subsample=0.85,
    colsample_bytree=0.85,
    objective="reg:squarederror",
    random_state=42
)

xgb_model.fit(X_train, y_res_train)

# 6. 递归预测测试集残差
residual_history = list(train_residual)
xgb_residual_pred = []

for i in range(len(y_test)):
    current_time = train_size + i

    residual_lags = np.array(residual_history[-lag:][::-1])
    calendar = time_features(np.array([current_time]))[0]

    x_input = np.r_[residual_lags, calendar].reshape(1, -1)
    residual_pred = xgb_model.predict(x_input)[0]

    xgb_residual_pred.append(residual_pred)

    # 将预测残差加入历史，供下一步递归预测
    residual_history.append(residual_pred)

xgb_residual_pred = np.array(xgb_residual_pred)

# 最终融合预测
hybrid_pred = arima_test_pred + xgb_residual_pred

# 7. 评估结果
def rmse(y_true, y_pred):
    return np.sqrt(mean_squared_error(y_true, y_pred))

print("ARIMA MAE:", round(mean_absolute_error(y_test, arima_test_pred), 3))
print("ARIMA RMSE:", round(rmse(y_test, arima_test_pred), 3))

print("Hybrid MAE:", round(mean_absolute_error(y_test, hybrid_pred), 3))
print("Hybrid RMSE:", round(rmse(y_test, hybrid_pred), 3))

# 8. 图一：预测结果与滚动误差
test_dates = test.index

arima_error = y_test - arima_test_pred
hybrid_error = y_test - hybrid_pred

rolling_window = 14
arima_roll_mae = pd.Series(np.abs(arima_error)).rolling(
    rolling_window
).mean()
hybrid_roll_mae = pd.Series(np.abs(hybrid_error)).rolling(
    rolling_window
).mean()

fig, axes = plt.subplots(
    2, 1, figsize=(14, 9), sharex=True,
    gridspec_kw={"height_ratios": [2, 1]}
)

axes[0].plot(train.index, y_train, label="训练集", color="gray")
axes[0].plot(test_dates, y_test, label="真实值", color="black", linewidth=2)
axes[0].plot(
    test_dates, arima_test_pred,
    label="ARIMA预测", linestyle="--", color="#1f77b4"
)
axes[0].plot(
    test_dates, hybrid_pred,
    label="ARIMA + XGBoost",
    color="#d62728", linewidth=2
)
axes[0].axvline(
    test_dates[0], color="green",
    linestyle=":", label="训练/测试分界线"
)
axes[0].set_title("ARIMA与融合模型的预测结果")
axes[0].set_ylabel("销量")
axes[0].legend()
axes[0].grid(alpha=0.3)

axes[1].plot(
    test_dates, arima_roll_mae,
    label="ARIMA 14日滚动MAE",
    color="#1f77b4"
)
axes[1].plot(
    test_dates, hybrid_roll_mae,
    label="融合模型14日滚动MAE",
    color="#d62728"
)
axes[1].fill_between(
    test_dates, 0, np.abs(hybrid_error),
    color="#d62728", alpha=0.15
)
axes[1].set_title("预测误差的时间变化")
axes[1].set_ylabel("误差")
axes[1].legend()
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

# 9. 图二：残差自相关与残差修正效果
hybrid_test_error = y_test - hybrid_pred
arima_test_error = y_test - arima_test_pred

acf_arima = acf(train_residual, nlags=30, fft=True)
acf_hybrid = acf(hybrid_test_error, nlags=30, fft=True)

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

axes[0].stem(
    range(len(acf_arima)), acf_arima,
    linefmt="C0-", markerfmt="C0o",
    basefmt="gray", label="ARIMA训练残差"
)
axes[0].stem(
    range(len(acf_hybrid)), acf_hybrid,
    linefmt="C1--", markerfmt="C1s",
    basefmt="gray", label="融合模型测试误差"
)
axes[0].axhline(0, color="black", linewidth=0.8)
axes[0].set_title("误差自相关变化")
axes[0].set_xlabel("滞后阶数")
axes[0].set_ylabel("自相关系数")
axes[0].legend()
axes[0].grid(alpha=0.3)

axes[1].scatter(
    arima_test_error,
    xgb_residual_pred,
    c=np.arange(len(y_test)),
    cmap="viridis",
    alpha=0.75
)
axes[1].axhline(0, color="black", linewidth=0.8)
axes[1].axvline(0, color="black", linewidth=0.8)
axes[1].set_title("XGBoost对ARIMA误差的修正")
axes[1].set_xlabel("ARIMA实际误差")
axes[1].set_ylabel("XGBoost预测残差")
axes[1].grid(alpha=0.3)

plt.tight_layout()
plt.show()

第一，XGBoost的训练目标不是原始销量，而是：

这样做的好处是，XGBoost不用重新学习完整的趋势，只需要专注于ARIMA没有处理好的部分。

第二，代码中的residual_lags表示过去7天的预测误差。比如当前要预测第100天，就把第93到99天的残差作为输入。

第三，测试阶段采用递归预测。也就是说，第一个测试点预测出的残差，会被放入历史记录中，参与下一个测试点的预测。这样没有使用未来真实值，评估结果更加接近真实部署场景。

第四，时间特征使用了正弦和余弦转换：

它们可以表示7天周期。直接把星期几写成1到7，会让模型误以为周日和周一相距很远，但实际上它们在周期上是相邻的。

注意事项

ARIMA的阶数 不要盲目调得很复杂。可以先通过ACF、PACF和AIC选择一个合理范围，再用时间序列交叉验证确认效果。

XGBoost的lag也不是越大越好。如果数据周期是7天，可以优先尝试7、14、21这样的滞后长度。如果数据量比较小，滞后特征过多反而容易过拟合。

还要特别注意一个问题：如果有天气、价格、活动等外部变量，必须确认预测时这些变量是否提前已知。已知的变量可以加入ARIMA的外生变量或XGBoost特征，未来无法获得的变量不能直接使用。

另外，融合模型不一定永远比ARIMA好。最可靠的方式不是看训练集误差，而是使用滚动窗口验证，对ARIMA、XGBoost和融合模型进行公平比较。

最后

ARIMA+XGBoost的核心逻辑并不复杂：ARIMA先学习趋势和线性时间依赖，XGBoost再学习ARIMA残差中的非线性规律，最后把两部分结果加起来。

在实际项目中，你还可以继续尝试两条方向：一是加入节假日、价格、天气等外部特征；二是使用滚动预测和多个时间窗口，比较不同融合方式在长期预测中的稳定性。

总之就是说，这不是简单地把两个模型拼在一起，而是让它们各自负责自己擅长的部分。