ML-For-Beginners 时间序列预测实战：用 ARIMA/SARIMA 模型预测电力负荷数据

本篇技术指南聚焦 ML-For-Beginners 课程库 7-TimeSeries/2-ARIMA 一课的完整知识体系：从平稳性、差分等时间序列核心概念入手，逐步讲解如何用 statsmodels 的 SARIMAX 为小时级电力负荷数据（GEFCom 2014）搭建 ARIMA 模型，并借助滑动窗口式的 Walk-Forward 验证与 MAPE 指标量化预测精度。读完你将掌握一套可直接照搬到其他时序数据上的 ARIMA 建模、训练、验证与评估实战流程。

为什么用 ARIMA：从非平稳数据到可预测的时间序列

ARIMA（AutoRegressive Integrated Moving Average，自回归积分滑动平均）是一类经典的统计时间序列模型，特别适合处理存在趋势与季节波动、分布随时间漂移的非平稳数据。本课程使用 2012—2014 年逐小时记录的电力负荷数据作为练习对象，预测未来时段的用电负荷，这正是工业与商业中典型的预测场景。

在该课的前一讲（7-TimeSeries/1-Introduction）中，你已经了解了时间序列预测的基础概念并加载了负荷数据；这一讲则把重点放到「如何真正构造一个 ARIMA 模型」。

两个绕不开的概念：平稳性与差分

• 平稳性（Stationarity）：从统计意义上看，平稳数据指其分布不随时间平移而改变。非平稳数据往往由于趋势（trend）或季节性（seasonality）而产生波动，需要先做变换才能被分析。例如季节性会给数据引入周期性波动，可通过「季节性差分」来消除。
• 差分（Differencing）：差分是将非平稳数据转化为平稳数据的核心手段——通过移除不恒定的趋势来稳定时序的均值。正如文献（Shixiong 等）所述：“差分移除了时间序列的水平变化，消除了趋势与季节性，从而使时间序列的均值趋于稳定。”ARIMA 中的字母 I（Integrated，积分） 正是指应用差分步数来消除非平稳性的过程，这也是 ARIMA 区别于 ARMA 等模型的关键。

拆解 ARIMA：AR、I、MA 分别做了什么

组成	全称	含义	通俗理解
AR	Autoregressive（自回归）	用目标变量**过去的滞后值（lags）**回归自身	模型“向后看”，依据历史观测推断当前值。例如月销量的当前值可由前几月的销量来解释
I	Integrated（积分/差分）	对序列施加差分步数以消除非平稳性	对应差分次数，决定对原始数据做几阶差分
MA	Moving Average（滑动平均）	输出变量由当前及过去滞后期的误差项决定	模型通过历史预测误差的组合来修正当前预测

一句话总结：ARIMA 的目标就是构建一个模型，尽可能贴合给定时间序列数据特有的形态（趋势、季节与噪声），并据此外推未来。

实战环境与数据准备

动手前先确认运行环境。本课配套的 conda 环境定义见 7-TimeSeries/2-ARIMA/working/common/environment.yaml，其中锁定的关键依赖包括 statsmodels==0.9.0、scikit-learn==0.20.3、pandas==0.23.4、numpy==1.16.2，并通过 pip 额外安装 pyramid-arima==0.8.1（供可选的 auto_arima 使用）。环境创建命令：

conda env create -f environment.yaml
conda activate dlts
python -m ipykernel install --user --name dlts --display-name "Python (dlts)"

说明：环境锁定版本较旧，在较新的 Python/依赖环境下若遇兼容问题，可仅安装本课必需的核心库，例如 pip install statsmodels scikit-learn pandas matplotlib，再运行下方代码。

加载依赖库

进入本课工作目录 7-TimeSeries/2-ARIMA/working，其中的 notebook.ipynb 是本节练习的起始骨架——它的首个单元即安装 statsmodels。随后在 notebook 中加载全部依赖：

import os
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import datetime as dt
import math

from pandas.plotting import autocorrelation_plot
from statsmodels.tsa.statespace.sarimax import SARIMAX
from sklearn.preprocessing import MinMaxScaler
from common.utils import load_data, mape
from IPython.display import Image

%matplotlib inline
pd.options.display.float_format = '{:,.2f}'.format
np.set_printoptions(precision=2)
warnings.filterwarnings("ignore")  # 忽略告警信息

其中 SARIMAX 来自 statsmodels 状态空间模块，是后续建模的核心；MinMaxScaler 用于数据归一化；load_data 与 mape 来自本课公共工具模块 common/utils.py。

加载能源负荷数据

课程数据为逐小时的电力负荷记录（列 load，另含温度列 temp），存放于 7-TimeSeries/data/energy.csv，时间跨度从 2012-01-01 到 2014 年末。工具函数 load_data() 的实现（见 utils.py）值得留意：它在读取 CSV 后会以 pd.date_range(..., freq='H') 按小时重新索引，将时间轴补全为最小到最大时间戳的完整连续序列——这样任何缺失时段都会暴露出来（本数据集恰好无缺失），为后续按时间对齐的差分与移位运算打好了基础：

energy = load_data('./data')[['load']]
energy.head(10)

绘制 2012 年 1 月至 2014 年 12 月的全部负荷曲线：

energy.plot(y='load', subplots=True, figsize=(15, 8), fontsize=12)
plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()

划分训练集与测试集：杜绝未来数据泄漏

时间序列建模与普通监督学习最大的不同在于：训练集必须在时间上先于测试集，否则模型“看到了未来”，评估结果将毫无意义。本课将 2014 年 9 月 1 日至 10 月 31 日划为训练区间，测试区间取 2014 年 11 月 1 日至 12 月 31 日。代码中用两个时间戳作为分界标记（注意：训练起点直接取 '2014-11-01 00:00:00'，而测试起点 test_start_dt 被刻意延后到 2014-12-30，其意图是把测试集限定在更靠后的近期小窗口内——因为负荷数据表现出强季节性，近期样本与当前形态最相似，小训练窗口足以捕捉规律）：

train_start_dt = '2014-11-01 00:00:00'
test_start_dt = '2014-12-30 00:00:00'

可视化两个集合的时间关系（蓝色为训练段、橙色为测试段），可直观看到训练窗口相对整个历史确实很小：

energy[(energy.index < test_start_dt) & (energy.index >= train_start_dt)][['load']].rename(columns={'load':'train'}) \
    .join(energy[test_start_dt:][['load']].rename(columns={'load':'test'}), how='outer') \
    .plot(y=['train', 'test'], figsize=(15, 8), fontsize=12)
plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()

补充说明：由于用于拟合 ARIMA 的函数在拟合过程中会执行样本内验证（in-sample validation），本课不再单独切出验证集。

过滤与缩放数据

按时间区间过滤出 train 与 test 两个 DataFrame（仅保留 load 列与日期索引）：

train = energy.copy()[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']]
test = energy.copy()[energy.index >= test_start_dt][['load']]

print('Training data shape: ', train.shape)
print('Test data shape: ', test.shape)

输出形状为：

Training data shape:  (1416, 1)
Test data shape:  (48, 1)

训练集 1416 行即 59 天 × 24 小时，测试集 48 行即 2 天 × 24 小时。接着用 MinMaxScaler 将训练数据缩放到 (0, 1) 区间，再用同一 scaler 变换测试数据（注意：fit_transform 只应在训练集上调用，测试集仅 transform，避免统计信息泄漏）：

scaler = MinMaxScaler()
train['load'] = scaler.fit_transform(train)
train.head(10)

test['load'] = scaler.transform(test)
test.head()

对比原始与缩放后训练数据的直方图分布（两者形状一致、仅横轴刻度不同）：

实现 ARIMA / SARIMA 模型

建模分三步走：定义模型 → fit 拟合 → forecast 预测。

理解 p、d、q 与季节性参数 P、D、Q

ARIMA 的三个核心超参数分别刻画时序的趋势、季节与噪声形态：

• p：对应模型的自回归部分，涉及对过去值的依赖阶数；
• d：对应模型的积分/差分部分，表示对时间序列施加的差分次数（还记得上文提到的差分吗？）；
• q：对应模型的滑动平均部分，表示对过去误差项依赖的阶数。

当数据存在季节性（本例正是如此）时，需要升级为季节性 ARIMA，即 SARIMA，并额外指定一组与 p、d、q 同义、但作用于季节成分的参数：

• P：季节性自回归阶数；
• D：季节性差分次数；
• Q：季节性滑动平均阶数；
• 外加一个季节周期长度 s（本课为 24，即一天的 24 个小时）。

设定预测视界并手动调参

先用一个「预测视界（horizon）」变量声明要向前预测几步。下面尝试预测未来 3 小时：

# 指定向前预测的步数
HORIZON = 3
print('Forecasting horizon:', HORIZON, 'hours')

参数选择是本课点名的难点：它带有一定主观性且耗时。课程提示，如果想自动化地搜索最优参数，可借助 pyramid（pmdarima 前身）提供的 auto_arima() 函数。这里先手动尝试一组参数来寻找较优模型：

order = (4, 1, 0)
seasonal_order = (1, 1, 0, 24)

model = SARIMAX(endog=train, order=order, seasonal_order=seasonal_order)
results = model.fit()

print(results.summary())

上述代码中：endog 传入训练序列；order=(4, 1, 0) 表示 4 阶自回归、1 次差分、无滑动平均项；seasonal_order=(1, 1, 0, 24) 表示季节自回归 1 阶、季节差分 1 次、季节周期为 24 小时。fit() 完成后打印的 summary() 会输出包含系数估计、显著性（P 值）、AIC/BIC 等信息在内的模型诊断表，可据此判断模型质量。

用 Walk-Forward 验证模拟真实预测节奏

单个模型一次拟合显然不足以评估泛化能力。时间序列模型在真实生产中往往每来一条新数据就重训一次，这正是 Walk-Forward（向前滚动）验证所模拟的过程：从序列起点开始训练 → 预测下一步 → 用已知真实值评估该步预测 → 把真实值并入训练窗口 → 重复，直至覆盖整个测试期。

提示：为兼顾训练效率，应当保持训练窗口长度固定——每向窗口尾部追加一个新观测，就从窗口头部移除一个最老的观测。

Walk-Forward 验证被本课称为时序模型评估的“金标准”，适合规模较小的数据；缺点是需反复重训大量模型，在大规模场景下有可观的计算开销。

第一步，把测试序列按 HORIZON 平移出多个“未来列”，构造监督学习所需的输入-目标对（load+1、load+2 分别代表 t+1、t+2 时刻的目标值）：

test_shifted = test.copy()

for t in range(1, HORIZON+1):
    test_shifted['load+'+str(t)] = test_shifted['load'].shift(-t, freq='H')

test_shifted = test_shifted.dropna(how='any')
test_shifted.head(5)

平移后的数据形态如下——每一行都携带了接下来 3 个小时的“答案”：

		load	load+1	load+2
2014-12-30	00:00:00	0.33	0.29	0.27
2014-12-30	01:00:00	0.29	0.27	0.27
2014-12-30	02:00:00	0.27	0.27	0.30
2014-12-30	03:00:00	0.27	0.30	0.41
2014-12-30	04:00:00	0.30	0.41	0.57

第二步，在滑动窗口循环中逐条滚动训练并预测（这里改用 order=(2, 1, 0)；training_window = 720 表示保留最近 30 天即 720 小时作为固定训练窗口）：

%%time
training_window = 720  # 用 30 天（720 小时）做训练

train_ts = train['load']
test_ts = test_shifted

history = [x for x in train_ts]
history = history[(-training_window):]

predictions = list()

order = (2, 1, 0)
seasonal_order = (1, 1, 0, 24)

for t in range(test_ts.shape[0]):
    model = SARIMAX(endog=history, order=order, seasonal_order=seasonal_order)
    model_fit = model.fit()
    yhat = model_fit.forecast(steps=HORIZON)
    predictions.append(yhat)
    obs = list(test_ts.iloc[t])
    # 移动训练窗口：末尾加入真实观测，头部弹出最旧观测
    history.append(obs[0])
    history.pop(0)
    print(test_ts.index[t])
    print(t+1, ': predicted =', yhat, 'expected =', obs)

可以看到训练逐步进行：每一时刻的 predicted 都是长度为 3 的数组（对应 t+1、t+2、t+3），而 expected 是同一时刻的 3 个真实目标值：

2014-12-30 00:00:00
1 : predicted = [0.32 0.29 0.28] expected = [0.32945389435989236, 0.2900626678603402, 0.2739480752014323]

2014-12-30 01:00:00
2 : predicted = [0.3  0.29 0.3 ] expected = [0.2900626678603402, 0.2739480752014323, 0.26812891674127126]

2014-12-30 02:00:00
3 : predicted = [0.27 0.28 0.32] expected = [0.2739480752014323, 0.26812891674127126, 0.3025962399283795]

汇总预测结果并还原到原始量纲

由于模型在 (0, 1) 区间上训练，评估前必须把预测与真实值通过 scaler.inverse_transform 还原为真实负荷量纲（兆瓦级），同时用 pd.melt 把 t+1/t+2/t+3 各列展开为长表，方便后续按 horizon 分组统计：

eval_df = pd.DataFrame(predictions, columns=['t+'+str(t) for t in range(1, HORIZON+1)])
eval_df['timestamp'] = test.index[0:len(test.index)-HORIZON+1]
eval_df = pd.melt(eval_df, id_vars='timestamp', value_name='prediction', var_name='h')
eval_df['actual'] = np.array(np.transpose(test_ts)).ravel()
eval_df[['prediction', 'actual']] = scaler.inverse_transform(eval_df[['prediction', 'actual']])
eval_df.head()

还原后的结果（预测值与实际负荷非常接近，差距通常在几十个单位以内）：

		timestamp	h	prediction	actual
0	2014-12-30	00:00:00	t+1	3,008.74	3,023.00
1	2014-12-30	01:00:00	t+1	2,955.53	2,935.00
2	2014-12-30	02:00:00	t+1	2,900.17	2,899.00
3	2014-12-30	03:00:00	t+1	2,917.69	2,886.00
4	2014-12-30	04:00:00	t+1	2,946.99	2,963.00

用 MAPE 量化模型精度

肉眼观察曲线吻合度之外，还需要一个客观指标。本课采用 MAPE（平均绝对百分比误差）——它把每一时刻的真实值与预测值之差的绝对值除以真实值，再对所有预测点取平均：

对应代码实现如下（本课公共工具中的 mape() 函数即此公式的封装，见 utils.py）：

if(HORIZON > 1):
    eval_df['APE'] = (eval_df['prediction'] - eval_df['actual']).abs() / eval_df['actual']
    print(eval_df.groupby('h')['APE'].mean())

单步（t+1）预测的 MAPE：

print('One step forecast MAPE: ', (mape(eval_df[eval_df['h'] == 't+1']['prediction'], eval_df[eval_df['h'] == 't+1']['actual']))*100, '%')

One step forecast MAPE:  0.5570581332313952 %

多步预测（t+1、t+2、t+3 全体）的 MAPE：

print('Multi-step forecast MAPE: ', mape(eval_df['prediction'], eval_df['actual'])*100, '%')

Multi-step forecast MAPE:  1.1460048657704118 %

解读口径：MAPE 越低越好——例如 MAPE 为 10 表示预测平均偏离真实值约 10%。上面单步约 0.56%、多步约 1.15% 的成绩说明该模型已相当精准，且符合“预测视界越长、误差累积越大”的直觉。

可视化实际负荷与预测曲线

将实际值与各 horizon 的预测绘在同一张图中：红色粗线为真实负荷，蓝色半透明曲线族为 t+1/t+2/t+3 的预测，透明度与线宽随 horizon 增加而递减，可清楚看出预测的衰减：

if(HORIZON == 1):
    ## 绘制单步预测
    eval_df.plot(x='timestamp', y=['actual', 'prediction'], style=['r', 'b'], figsize=(15, 8))

else:
    ## 绘制多步预测
    plot_df = eval_df[(eval_df.h=='t+1')][['timestamp', 'actual']]
    for t in range(1, HORIZON+1):
        plot_df['t+'+str(t)] = eval_df[(eval_df.h=='t+'+str(t))]['prediction'].values

    fig = plt.figure(figsize=(15, 8))
    ax = plt.plot(plot_df['timestamp'], plot_df['actual'], color='red', linewidth=4.0)
    ax = fig.add_subplot(111)
    for t in range(1, HORIZON+1):
        x = plot_df['timestamp'][(t-1):]
        y = plot_df['t+'+str(t)][0:len(x)]
        ax.plot(x, y, color='blue', linewidth=4*math.pow(.9,t), alpha=math.pow(0.8,t))

    ax.legend(loc='best')

plt.xlabel('timestamp', fontsize=12)
plt.ylabel('load', fontsize=12)
plt.show()

红蓝曲线几乎贴合，说明这个 ARIMA 模型在 2014-12-30 至 12-31 的测试窗口上表现出良好的预测精度——至此你已完整走通一遍时序建模闭环。

进阶练习与任务

• 挑战：本课只用 MAPE 衡量精度，时序预测领域还有其他评估方法（如 MAD、MSD、MASE 等）。请自行查阅相关准确性度量资料并描述它们各自的适用场景与优缺点。
• 课后作业：参考本课作业说明 assignment.md——用一份全新的时间序列数据（课程推荐 Duke 大学的公开时间序列数据集）重新构建一个 ARIMA 模型，在 notebook 中完成数据可视化、建模、MAPE 精度评估并配以解释性注释。
• 延伸阅读：学有余力可继续探索本课程库时间序列单元的其他建模路线，例如用支持向量回归做时序预测的 3-SVR 一课；7-TimeSeries/common/utils.py 中的 TimeSeriesTensor 与 create_evaluation_df 还展示了面向 RNN/深度学习时序模型的张量化数据预处理方案，可作为对照了解经典统计模型之外的另一条技术路线。

小结

本课的核心方法论可以沉淀为一条可复用的 ARIMA 实战链路：理解平稳性并决定差分阶数（d）→ 依据序列特征指定自回归（p）与滑动平均（q）阶数 → 存在周期性时扩展为 SARIMA 并配置季节参数（P、D、Q、s）→ 按时间先后划分训练/测试集并做 MinMax 缩放 → 用 Walk-Forward 滚动重训做预测 → 以 MAPE 等指标评估并在原始量纲下可视化验证。这套流程的全部代码骨架见 working/notebook.ipynb，完整可运行实现见同目录 solution/notebook.ipynb，数据文件位于 7-TimeSeries/data/energy.csv，非常适合作为你处理电力、销量、流量等带强周期性的小时级/日级数据的起点。