时间序列：季节性多因子时间序列预测模型SARIMAX

写在前面：
首先感谢兄弟们的订阅，让我有创作的动力，在创作过程我会尽最大能力，保证作品的质量，如果有问题，可以私信我，让我们携手共进，共创辉煌。

路虽远，行则将至；事虽难，做则必成。只要有愚公移山的志气、滴水穿石的毅力，脚踏实地，埋头苦干，积跬步以至千里，就一定能够把宏伟目标变为美好现实。

Hello，大家好。前面我们完整讲解了 ARIMA（基础时序） 、ARIMAX（带外部因子时序）和SARIMA（带季节性），今天来看看SARIMAX。

在进入主题之前，给大家分享一下今日语录（摘抄自网络）：
你活着的时候，没有几个人在意你，你死了以后，也没几个人记得你。所以，请你务必记住，除了生病以外，你所感受到的痛苦，都是你的情绪带给你的，而非真实存在的。你现在要在乎的，无非就三件事：你的钱、你的健康、你内心的平静。

以下内容均结合参考资料与个人实操理解原创撰写，不存在滥用原创问题。

什么是SARIMAX

SARIMAX = SARIMA + 外生变量= 季节性时间序列 + 外部影响因子

• ARIMA：只看历史数据，无周期、无外部变量
• ARIMAX：历史数据 + 外部变量，无周期
• SARIMA：历史数据 + 周期规律，无外部变量
• SARIMAX：历史数据 + 周期规律 + 外部变量
  它能同时解决数据的趋势、季节性波动、外部因素干扰三大问题，是工业界最常用的时间序列预测模型之一。

SARIMAX表达式：
$$SARIMAX(p,d,q)(P,D,Q,s)+X\beta$$

• $(p,d,q)$：非季节性部分（和 ARIMA 一样）
• $(P,D,Q,s)$：季节性部分（和 SARIMA 一样）
• $X$：外生变量矩阵（和 ARIMAX 一样）
• $\beta$：外生变量的影响系数（模型自动学习）

适用场景

真实业务中，很多数据同时受两类因素影响：
周期性因素：如月度销量的淡旺季、日客流的周度规律
外部因素：如促销活动、天气温度、节假日政策
普通ARIMA无法捕捉周期和外部因素，ARIMAX无法捕捉周期，SARIMA无法捕捉外部因素，而SARIMAX同时支持两者。在以下业务场景适用：

• 月度 / 季度销量预测（淡旺季 + 促销活动）
• 电网负荷预测（季节变化 + 温度影响）
• 景区客流预测（周度周期 + 节假日）
• 电商平台流量预测（日 / 周周期 + 营销活动）

代码实战

下面的代码包含数据构建、自动定阶、训练、预测、评估、可视化。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.metrics import mean_squared_error, mean_absolute_error
import joblib
import warnings

warnings.filterwarnings('ignore')

# ===================== 1. 构造业务数据（带季节性 + 外生变量） =====================
np.random.seed(42)
n = 120  # 10年月度数据
dates = pd.date_range("2015-01-01", periods=n, freq="M")

# 外生变量：促销力度、平均气温（可替换为真实业务数据）
promote = np.linspace(0, 10, n) + np.random.normal(0, 0.5, n)
temp = 20 + 10 * np.sin(np.linspace(0, 4 * np.pi, n)) + np.random.normal(0, 1, n)

# 目标序列：销量 = 趋势 + 季节性 + 促销影响 + 温度影响 + 噪声
trend = np.linspace(10, 50, n)
season = 10 * np.sin(np.linspace(0, 10 * np.pi, n))
y = trend + season + 2 * promote + 0.5 * temp + np.random.normal(0, 2, n)

df = pd.DataFrame({
    "date": dates,
    "sales": y,
    "promote": promote,
    "temperature": temp
})
df.set_index("date", inplace=True)

df.to_csv("data_sarimax.csv")

plt.rcParams["font.family"] = ["SimHei"]
plt.rcParams["axes.unicode_minus"] = False
plt.figure(figsize=(12,6))
plt.plot(df["sales"], label="数据", color="#1f77b4")

# ===================== 2. 划分训练集、测试集 =====================
train = df.iloc[:-12]  # 前108个月训练
test = df.iloc[-12:]  # 最后12个月测试

# 分离内生变量（目标）和外生变量（影响因子）
train_endog = train["sales"]
train_exog = train[["promote", "temperature"]]
test_exog = test[["promote", "temperature"]]


# ===================== 3. 自动定阶（SARIMAX 最优参数） =====================
def best_sarimax_order(endog, exog, s=12):
    best_aic = float("inf")
    best_order = (1, 1, 1)
    best_seasonal_order = (1, 1, 1, s)

    # 遍历常见参数组合（可根据需求扩大范围）
    for p in range(3):
        for d in range(2):
            for q in range(3):
                for P in range(2):
                    for D in range(2):
                        for Q in range(2):
                            try:
                                model = SARIMAX(
                                    endog=endog,
                                    exog=exog,
                                    order=(p, d, q),
                                    seasonal_order=(P, D, Q, s),
                                    enforce_stationarity=False,
                                    enforce_invertibility=False
                                )
                                res = model.fit(disp=False)
                                if res.aic < best_aic:
                                    best_aic = res.aic
                                    best_order = (p, d, q)
                                    best_seasonal_order = (P, D, Q, s)
                            except:
                                continue
    return best_order, best_seasonal_order


# 月度数据，周期s=12
order, seasonal_order = best_sarimax_order(train_endog, train_exog, s=12)
print("最优非季节阶数 order：", order)
print("最优季节阶数 seasonal_order：", seasonal_order)

# ===================== 4. 训练 SARIMAX 模型 =====================
model = SARIMAX(
    endog=train_endog,
    exog=train_exog,
    order=order,
    seasonal_order=seasonal_order,
    enforce_stationarity=False,
    enforce_invertibility=False
)
model_fit = model.fit(disp=False)
print(model_fit.summary())

# ===================== 5. 预测未来12期 =====================
forecast = model_fit.get_forecast(steps=len(test), exog=test_exog)
pred_mean = forecast.predicted_mean
true_y = test["sales"]

# 修复点：置信区间正确取值方式
conf_int = forecast.conf_int()
lower = conf_int.iloc[:, 0]
upper = conf_int.iloc[:, 1]

# ===================== 6. 模型评估 =====================
rmse = np.sqrt(mean_squared_error(true_y, pred_mean))
mae = mean_absolute_error(true_y, pred_mean)
print(f"\n【模型评估指标】")
print(f"RMSE: {rmse:.2f}")
print(f"MAE: {mae:.2f}")

# ===================== 7. 可视化预测效果 =====================
plt.figure(figsize=(12, 6))
plt.plot(train.index, train["sales"], label="训练数据", color="#1f77b4")
plt.plot(test.index, true_y, label="真实销量", color="#2ca02c")
plt.plot(test.index, pred_mean, label="SARIMAX预测销量", color="#d62728", linewidth=2)
# 使用拆分后的上下界
plt.fill_between(test.index, lower, upper, color="pink", alpha=0.3, label="95%置信区间")
plt.title("SARIMAX 销量预测结果（季节性 + 外生变量）")
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

# ===================== 8. 模型保存与加载（工程部署） =====================
# joblib.dump(model_fit, "sarimax_best_model.pkl")
# print("✅ SARIMAX模型已保存成功")

# 加载模型（下次直接使用）
# loaded_model = joblib.load("sarimax_best_model.pkl")
# future_exog = pd.DataFrame(...)  # 准备未来外生变量
# pred = loaded_model.get_forecast(steps=12, exog=future_exog)

原始数据：

预测效果：

SARIMAX预测未来时，必须提供未来时间段的外生变量数据（比如未来12个月的促销计划、气温预测），这是它和SARIMA的核心区别。

说明：文章首发在微信公众号