用 ARIMA/SARIMA 进行时间序列预测——ML-For-Beginners 电力负荷预测实战指南

本篇技术指南围绕 ML-For-Beginners 仓库中《7-TimeSeries/2-ARIMA》课程展开，以 GEFCom2014 电力负荷数据（小时级，2012–2014）为研究对象，系统讲解 ARIMA（自回归积分滑动平均） 与 SARIMA（季节性 ARIMA） 的建模原理、数据准备、walk-forward 滚动验证与 MAPE 精度评估全流程。读完本文，你将能用 statsmodels 从零构建一个可复现的负荷预测模型，理解 p、d、q 与 P、D、Q 六类超参数的含义，并掌握时间序列模型“防止未来信息泄漏”的切分与评估规范。

一、课程背景与本讲目标

在时间序列系列课程的前一课《7-TimeSeries/1-Introduction》中，你已经初步了解时间序列预测，并加载了一份展示某时间段内电力负荷波动的数据集。本讲继续沿用这份数据，聚焦一个经典的统计建模方法：ARIMA（AutoRegressive Integrated Moving Average）。ARIMA 尤其适合拟合具有 非平稳性（non-stationarity） 的数据——这恰好是电力负荷数据的基本特征：负荷水平会随时间、季节与昼夜节律起伏。

读者在本讲将依次完成以下四个动作：

1. 用 statsmodels 实现 ARIMA/SARIMA 模型；
2. 完成训练集/测试集的时间顺序切分与 (0,1) 区间缩放；
3. 使用 walk-forward（前向滚动）验证逐时间步重训模型并预测；
4. 用 MAPE（平均绝对百分比误差） 定量评估并可视化结果。

仓库中可直接运行的示例位于 7-TimeSeries/2-ARIMA/working/notebook.ipynb（配套完整解在 7-TimeSeries/2-ARIMA/solution/notebook.ipynb），数据文件为 7-TimeSeries/data/energy.csv。

二、动手前必须理解的两个概念

要驾驭 ARIMA，先要掌握两个统计基础概念，它们是本讲“为什么电力负荷适合 ARIMA”的理论支点。

🎓 Stationarity（平稳性）。从统计学视角看，平稳性指数据的分布在时间平移后不改变。反过来说，非平稳数据会因趋势（trend）而呈现持续波动，必须先经变换才能被常规模型分析。例如**季节性（seasonality）**会在数据中引入周期性波动，可以通过“季节性差分（seasonal-differencing）”来消除。

🎓 Differencing（差分）。差分是将非平稳数据变换为平稳数据的过程，通过去除非常数趋势来实现。原文档引用的表述非常精辟：“差分移除了时间序列水平的变化，消除了趋势与季节性，从而稳定了时间序列的均值。”在 ARIMA 记号中，进行差分的阶数正是由参数 d 控制。

三、拆解 ARIMA：AR、I、MA 与季节性扩展

ARIMA 由三部分拼装而成，理解每一部分能帮你回答“它为什么能预测时间序列”：

• AR —— AutoRegressive（自回归）。自回归模型“向后看”，分析并利用数据的历史取值做推断，这些历史取值称为 lags（滞后项）。以“月度铅笔销量”为例，每个月的销量在数据集中是一个“演化变量（evolving variable）”，模型将该变量回归到它自身的滞后（即先前）取值上。
• I —— Integrated（积分）。与相似的 ARMA 模型不同，ARIMA 中的 “I” 指其*integrated（积分）*属性：当对数据施加差分步骤以消除非平稳性时，就称数据被“积分”了。
• MA —— Moving Average（滑动平均）。模型的滑动平均部分表示：输出变量由对当前及过去滞后值的观测来决定。

一句话总结：ARIMA 的目标就是让模型尽可能紧密地拟合时间序列数据所呈现的特殊形态。

需要补充的是，当数据存在季节效应时（电力负荷正是如此），需要使用季节 ARIMA（SARIMA），即引入一组大写参数：P、D、Q，它们与 p、d、q 描述相同的关联，但针对的是模型的季节分量。SARIMAX() 正是 statsmodels 中同时支持两者（并可选外生变量）的接口。

四、环境与运行入口

本讲所有练习代码都组织在两个目录中：带空格的 working（练习）与 solution（参考答案）。打开 7-TimeSeries/2-ARIMA/working/notebook.ipynb 即可开始。仓库同时提供了 conda 环境文件 7-TimeSeries/2-ARIMA/solution/common/environment.yaml（Python 3.6.6、statsmodels、scikit-learn、pandas 等固定版本），可用 conda env create -f environment.yaml 复现环境。

notebook 的第一段只需执行一行，加载本讲的核心依赖库 statsmodels：

pip install statsmodels

随后导入绘图与建模所需的库。注意其中两个自研工具 load_data 与 mape 来自仓库共享模块 7-TimeSeries/common/utils.py：

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") # specify to ignore warning messages

关于 load_data 的底层行为值得说明：从 7-TimeSeries/common/utils.py 的源码可见，它会读取 energy.csv 并以 parse_dates=['timestamp'] 解析时间列，然后以小时为频率（freq='H'）对时间戳做全区间重索引（reindex），从而可以显式暴露数据中缺失的时间点（该数据集中没有缺失）。同样的模块里还提供了 mape(predictions, actuals)，其实现为 (np.absolute(predictions - actuals) / actuals).mean()，即后文评估阶段直接复用的 MAPE 公式。

五、加载并可视化电力负荷数据

从 /data/energy.csv 加载数据到 Pandas DataFrame 并预览前 10 行：

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

由于数据文件实际位于仓库的 7-TimeSeries/data/energy.csv，若不在 notebook 目录下运行，可相应调整路径。预览可见每个小时一条记录，load 表示该小时电网负荷值（例如 2012-01-01 00:00 为 2,698，随后逐时波动）。

接着绘制 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()

现在，开始建模吧！

六、构造训练集与测试集

数据加载完成后，下一步是切分。基本规则与其他机器学习任务一致：在训练集上训练模型，训练结束后用测试集评估精度。但时间序列有一项特殊铁律：测试集必须覆盖比训练集更晚的时间段，否则模型会“偷看”到未来信息，评估结果将失去意义。

原文档与代码使用如下两个时间界标（注意：文档正文提到的“9 月至 10 月”与代码略有出入，实际以代码为准，训练数据从 2014-11-01 开始）：

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），因此这里不再单独划分验证集。

6.1 按时间段与列过滤

将原始数据集按时间与所需列 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)

输出确认了切分的规模（训练 1416 个点 ≈ 59 天 × 24 小时，测试 48 个点 = 2 天）：

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

6.2 缩放到 (0, 1) 区间

统计模型的数值稳定性高度依赖特征尺度，因此这里用 sklearn 的 MinMaxScaler 把负荷压到 [0,1]：先对训练集 fit_transform（拟合 min/max 并变换），得到缩放器后再用它变换测试集——绝不能用测试集的数据去拟合缩放器，否则同样属于信息泄漏。

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

对比缩放前后的分布（各画 100 个桶的直方图）：

energy[(energy.index >= train_start_dt) & (energy.index < test_start_dt)][['load']].rename(columns={'load':'original load'}).plot.hist(bins=100, fontsize=12)
train.rename(columns={'load':'scaled load'}).plot.hist(bins=100, fontsize=12)
plt.show()

原始负荷取值大体落在 2000–3500 之间，且分布形态（偏态）在缩放后得到保留：

原始数据

缩放后的数据

缩放器校准完毕后，用它变换测试数据：

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

观察表头可以看到，测试数据的 load 已落入 0.2–0.4 区间（例如 2014-12-30 00:00 为 0.33），与训练集同尺度。

七、实现 ARIMA：定义模型、拟合、预测

现在调用前面安装好的 statsmodels。实现流程分三步：

1. 定义模型：调用 SARIMAX()，传入 order=(p,d,q) 与 seasonal_order=(P,D,Q,m)；
2. 拟合：调用 fit() 让模型适配训练数据；
3. 预测：调用 forecast() 并指定要预测的步数（即 horizon，预测视野）。

六个超参数的作用归纳如下：

参数	含义
p	与模型自回归部分关联的参数，纳入过去的取值（滞后阶数）
d	与模型积分（差分）部分关联的参数，决定对序列施加差分的次数
q	与模型滑动平均部分关联的参数
P	对应 p 的季节版本，描述季节自回归成分
D	对应 d 的季节版本，描述季节差分的次数
Q	对应 q 的季节版本，描述季节滑动平均成分
m	季节周期的长度（对小时级数据取 24，即“天”这一季节周期）

7.1 设定预测视野（horizon）

# Specify the number of steps to forecast ahead
HORIZON = 3
print('Forecasting horizon:', HORIZON, 'hours')

7.2 手动选择一组参数

为 ARIMA 挑选最优参数本身略带主观且耗时。实践中可借助 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())

运行后会打印一张 SARIMAX Results 汇总表。以仓库 solution notebook 的输出为例：模型记号为 SARIMAX(4, 1, 0)x(1, 1, 0, 24)，1416 个观测，对数似然 3477.239，AIC −6942.477；系数表中 ar.L1 约 0.84、ar.L2 约 −0.52、ar.S.L24 约 −0.23（季节滞后 24 的系数显著），下方还附带了 Ljung-Box、Jarque-Bera 等残差诊断统计量。这些信息告诉你每一阶滞后项对当前负荷的贡献方向与显著性。

第一个模型建好了，接下来需要评估它的好坏。

八、Walk-forward 验证：时间序列评估的黄金标准

对时间序列模型，实践中每得到一个新观测都应重新训练，以便在每个时刻给出最佳预测。walk-forward 验证（前向滚动验证）正是模拟这一机制：从序列起点开始，用当前训练数据训练模型 → 预测下一个时间步 → 用已知真实值评估预测 → 把该真实值并入训练集 → 重复上述过程。

提示：为了更高效训练，应保持训练窗口长度固定——每往训练集末尾加入一个新观测，就从集合开头移除一个旧观测。

该过程对模型“在真实环境中的表现”给出更稳健的估计；代价是需要反复建模型，计算开销大。数据量小或模型简单时可接受，规模大时会成为瓶颈。尽管有成本，walk-forward 仍是时间序列模型评估的黄金标准，推荐在你的项目中使用。

8.1 为每个 horizon 步构造“未来真实值”列

用 shift 把测试序列按小时频率向前平移，构造出每个时刻对应的 t+1、t+2 真实值（列名 load+1、load+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)

得到的表结构如下（每行代表一个预测起点，load 是当前值，load+1/load+2 是未来 1、2 小时真实值）：

		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

可见数据按 horizon 目标发生了“水平错位”。

8.2 滚动窗口循环预测

保持训练窗口固定为 720（约 30 天），每轮用 SARIMAX 重新拟合并前向预测 HORIZON 步，然后滑动窗口：

%%time
training_window = 720 # dedicate 30 days (720 hours) for training

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])
    # move the training window
    history.append(obs[0])
    history.pop(0)
    print(test_ts.index[t])
    print(t+1, ': predicted =', yhat, 'expected =', obs)

注意此处把演示模型的阶数降为 order=(2,1,0)，以降低每轮重训的成本。你可以实时观察逐点滚动训练的过程：

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]

从第 1–3 行可见，缩放空间里的预测值（如 [0.32, 0.29, 0.28]）与真实值（0.329…, 0.290…, 0.273…）已经相当接近。作为参照，仓库 solution notebook 完整跑完 46 个测试点耗时约 2 分 36 秒（Wall time），这印证了 walk-forward 反复重训的计算代价。

8.3 将预测与真实负荷对齐比较

把所有预测整理成 DataFrame，并用 pd.melt 把“t+1、t+2、t+3”多列转成长表（每行一个“预测起点 × 水平”组合），最后用 scaler.inverse_transform 把缩放后的预测与真实值还原回原始负荷量纲（兆瓦级数字），便于人类阅读：

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

观察小时级预测与真实负荷：预测 3,008.74 对真实 3,023.00、预测 2,955.53 对真实 2,935.00……误差基本在几十 MW 以内。这有多准？下面用指标说话。

九、用 MAPE 检验模型精度

MAPE（Mean Absolute Percentage Error，平均绝对百分比误差） 以比率形式展示预测精度，公式为：

$$MAPE = \frac{1}{n}\sum_{t=1}^{n}\left|\frac{actual_{t} - predicted_{t}}{actual_{t}}\right|$$

即真实值与预测值之差除以真实值；该差值的绝对值对每个被预测时刻求和，再除以拟合点数 n。仓库中的 mape() 实现（见 7-TimeSeries/common/utils.py）与此公式完全对应。指标越低越好——若某预测的 MAPE 为 10，意味着平均偏离真实值约 10%。

9.1 按水平分组计算 APE

对多步预测（HORIZON>1），先逐行算 APE，再按 h（t+1/t+2/t+3）分组求均值，以查看“越远越难预测”的衰减规律：

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

9.2 单步预测的 MAPE

只取 h=='t+1' 的行，乘以 100 转成百分比：

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 %

9.3 多步预测的整体 MAPE

对整张 eval_df（包含 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 约 0.56%、整体约 1.15%，对电力负荷这种数量级在数千 MW 的序列而言是相当漂亮的精度——即使整体 1.15% 也只是平均 1% 出头的相对误差。

9.4 可视化预测 vs 真实

数值指标之外，图形更直观。绘制实际负荷曲线（红色）与各 horizon 的预测曲线（蓝色，透明度/线宽随 t 增大而递减，体现置信度递减）：

if(HORIZON == 1):
    ## Plotting single step forecast
    eval_df.plot(x='timestamp', y=['actual', 'prediction'], style=['r', 'b'], figsize=(15, 8))

else:
    ## Plotting multi step forecast
    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()

红色实线是真实负荷，蓝色由深到浅依次是 t+1、t+2、t+3 的预测。曲线几乎重叠，只在预测线颜色变浅处出现可察觉的微小偏离——这是一张精度相当好的模型图。

十、深入阅读与自主练习

🚀 挑战：探索更多精度指标

本讲只触及 MAPE。时间序列精度度量远不止这一种，例如 MAD（平均绝对偏差）、MSD（均方偏差）等都能从不同侧面刻画误差。可以调研这些指标各自的适用场景（是否对量纲敏感、是否会被零值分母放大等），并为你的模型补充计算与注释。

复习与自学

ARIMA 只是时间序列预测的入门方法。后续课程《7-TimeSeries/3-SVR》会引入机器学习回归器做同类任务。官方 forecasting 仓库中还有更多模型家族，包括指数平滑、动态回归等，值得对比它们在趋势、季节性与噪声建模上的差异。

作业

本讲配有一份独立作业《构造一个新的 ARIMA 模型》，见 7-TimeSeries/2-ARIMA/assignment.md。建议在动手作业时实践以下要点：

1. 保证测试时段严格晚于训练时段，杜绝未来信息泄漏；
2. 只对训练集 fit 缩放器，再用同一缩放器 transform 测试集；
3. 用 walk-forward 而非一次性“train → predict”来评估；
4. 先打印 results.summary() 检查系数显著性，再谈精度；
5. 用单步与多步两组 MAPE 分别报告，避免被平均掩盖远步误差。

结语

通过本讲，你完成了一条完整的小时级电力负荷 ARIMA 预测流水线：从理解平稳性与差分入手，识别出 p、d、q / P、D、Q 六参数并借助 SARIMAX 落地，随后用 (0,1) 缩放与时间有序切分保护评估的公平性，以 walk-forward 滚动重训逼近线上部署的真实节奏，最终用 MAPE（单步 0.56%、多步 1.15%）与预测曲线双重验收模型。这套“概念—数据—建模—验证—度量”的框架，同样适用于客流、销量、能源、交通等任何带季节周期的时间序列预测任务。