ML-For-Beginners 时序预测实战：用 ARIMA/SARIMA 建模电力负荷数据并评估精度

本指南对应 ML-For-Beginners 机器学习入门课程 7-TimeSeries 章节的第二课（含中文内容映射的原 阿拉伯语译文 与 英文原文），以 GEFCom 2014 电力负荷小时数据为对象，完整演示如何在 Jupyter Notebook 中借助 statsmodels 的 SARIMAX 从零构建自回归积分滑动平均模型：先厘清平稳性、差分、AR/I/MA 与 p、d、q / P、D、Q 参数等核心概念，再按“切分训练测试集 → MinMax 归一化 → 拟合 → walk-forward 滚动验证 → MAPE 评估 → 可视化”的完整链路跑通预测管线。读完你不仅能跑通 working/notebook.ipynb，还能举一反三用同一套方法论在自己的时序数据集上建模并量化精度。

课程背景与前置准备

本课之前，课程已介绍了时序预测的基本概念并加载了展示电力负荷波动的数据集。本课则聚焦于一种非常经典的建模方法——ARIMA（AutoRegressive Integrated Moving Average，自回归积分滑动平均）。ARIMA 模型特别适合拟合呈现出**非平稳性（non-stationarity）**的数据，而电力负荷数据恰恰因为存在趋势与周期性而普遍非平稳，因此是天然的实验对象。

运行环境：一份来自仓库的 conda 依赖清单

仓库在 environment.yaml 中给出了本课的运行环境，可用下述命令创建：

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

从文件内容可以看到课程环境的依赖版本组合，其中与本课直接相关的核心包是：

• python==3.6.6、pandas==0.23.4、numpy==1.16.2、matplotlib==3.0.0
• scikit-learn==0.20.3（提供 MinMaxScaler）
• statsmodels==0.9.0（提供 SARIMAX）
• pip 依赖中的 pyramid-arima==0.8.1（提供自动选参的 auto_arima 工具）

注意：上述版本为课程编写时的锁定版本，在你的环境中安装时，应使用与当前 Python 版本兼容的新版本组合；课程代码的核心 API（SARIMAX、MinMaxScaler）在后续版本中保持稳定。本课全部实操代码集中在 working/notebook.ipynb（另有已跑通的 solution/notebook.ipynb 可供对照）。

数据来源：GEFCom 2014 电力负荷

课程使用的 energy.csv 来自 GEFCom 2014 全球能源预测竞赛的负荷预测扩展赛道数据。从文件表头可见其结构为每小时一条记录：

timestamp,load,temp
2012-01-01 00:00:00,2698.0,32.0
2012-01-01 01:00:00,2558.0,32.666666667
...

其中 load 为电力负荷（本课只取该列建模），temp 为温度，时间范围从 2012 年 1 月起逐小时记录。若需了解原始数据如何由 GEFCom2014 Excel 抽取为 energy.csv，可查看仓库中的 extract_data.py。

预备概念：平稳性与差分

要理解 ARIMA，必须先掌握两个统计学概念：

• 平稳性（Stationarity）：从统计视角看，平稳指数据的分布在时间平移后不发生变化。反过来，非平稳数据会因趋势而表现出波动，必须先做变换才能分析。例如季节性会给数据引入波动，可通过“季节性差分（seasonal-differencing）”消除。
• 差分（Differencing）：差分是将非平稳数据变换为平稳数据的过程——通过去掉非常数趋势使序列平稳。其效果可以概括为“差分移除时间序列水平的变化，消除趋势与季节性，从而使时间序列的均值稳定”。ARIMA 中的“积分（Integrated）”指的就是执行差分步骤以消除非平稳性的过程。

拆解 ARIMA：三个字母分别代表什么

ARIMA 之所以能贴合时序数据，是因为它的三个组成部分各自处理序列的一个侧面：

• AR — 自回归（AutoRegressive）：自回归模型“回看”过去，用数据中先前的观测值（这些滞后值称为 lags，如“滞后 1 期”即上一小时/上一个月）对当前值做回归。例如某商品月度销量的数据中，每月总销量被视为一个“演化变量（evolving variable）”，模型构建的思路就是让“关注的演化变量对自身的滞后值做回归”。
• I — 积分（Integrated）：区别于形态相近的 ARMA 模型，ARIMA 中的字母 I 指其“积分”侧面。当通过差分步骤消除了非平稳性时，数据即被称为被“积分”了，对应的差分阶数就是参数 d。
• MA — 滑动平均（Moving Average）：模型的这一侧面指输出变量由当前值和对滞后项的当前/过去观测共同决定，用于刻画随机扰动（噪声）的累积影响。

一句话总结： ARIMA 的目标是让模型尽可能贴合时间序列数据特有的形态（趋势、季节性、噪声），从而获得可靠的预测。

动手实战：在 Notebook 中构建 ARIMA 模型

打开 working/notebook.ipynb 并逐步运行即可复现完整流程。核心步骤依次为：加载库 → 读数据并绘图 → 划分训练/测试集 → 归一化 → 定义并拟合 SARIMAX → 滚动预测 → 评估。

第一步：加载全部依赖库

statsmodels 是 ARIMA 建模的核心库。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") # specify to ignore warning messages

其中 load_data 与 mape 两个函数来自课程自带的工具模块 common/utils.py。

第二步：加载并观察数据

load_data 函数在 utils.py 中实现。它读取 CSV 后将 timestamp 设为索引，并用 pd.date_range(..., freq='H') 以小时为频率重建索引——这保证时间轴上每个小时都有记录，从而暴露数据中缺失的时间段（该数据集恰好无缺失）：

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

本课仅取 load 一列。随后绘制 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()

图上可清楚看到明显的年度季节性与日周期波动——这正是后面要引入季节性差分（seasonal_order）的原因。

第三步：切分训练集与测试集

必须确保测试集在时间上晚于训练集，这样模型才不会从“未来”时段泄露信息。按原始代码的定义：

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

（说明：课程文案与代码在起止日期描述上略有出入，请以代码为准——训练窗口实际从 2014-11-01 起，测试集从 2014-12-30 起覆盖约两天逐小时数据。）由于数据反映的是逐小时电力消耗，存在很强的季节性模式，且近期的消耗与更早的历史最相似，因此用一段相对较短的近期时间窗口训练通常就够了。

用下面的代码可视化训练与测试两段数据的衔接：

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 的 fit() 在拟合期间会做样本内验证，因此这里不再单独划分验证集。

第四步：过滤与归一化数据

训练前需要对数据做两步处理：过滤——只保留所需时间段与所需列；归一化——用 MinMaxScaler 将数据压缩到区间 (0, 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)

输出：

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

训练集 1416 行对应约 59 天逐小时数据，测试集 48 行对应两天。随后在训练集上拟合缩放器并做归一化：

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

对比原始数据与归一化后数据的分布直方图：

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()

左图为原始数据，右图为缩放到 (0,1) 区间后的数据。可以看到归一化只改变数值刻度、不改变分布形状。

关键点是：测试集必须使用训练集拟合出的同一个 scaler 进行 transform，绝不能对测试集重新 fit（否则会引入未来信息）：

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

参数体系：p、d、q 与季节性 P、D、Q

在实现前，先吃透 SARIMAX 的参数含义。普通 ARIMA 有三个参数，分别对应时序的三大侧面（季节性、趋势、噪声）：

参数	全称关联	含义
p	AutoRegressive 阶数	对应模型的自回归侧面，纳入序列过去的滞后值
d	Integrated（差分）阶数	对应模型的积分侧面，决定对序列施加多少阶差分以消除非平稳
q	Moving Average 阶数	对应模型的滑动平均侧面，刻画过去预测误差（噪声）的影响

季节性情形下的扩展（SARIMA）：如果数据带有季节性侧面——本课数据正是如此——应改用季节性 ARIMA 模型（SARIMA），此时需要追加第二组参数：

• P、D、Q：与 p、d、q 描述相同的关联，但作用于模型的季节性分量；
• 外加一个季节性周期长度 s（seasonal length）。本课数据为逐小时记录、以日为周期，故 s = 24。

由此 seasonal_order = (P, D, Q, s)。这也是本课代码中 (1, 1, 0, 24) 的由来：对每日同刻做 1 阶季节性差分以剔除日周期。

实现 SARIMAX：定义、拟合、预测三步走

执行 ARIMA 的标准流程为：

1. 用 SARIMAX() 定义模型，传入 endog（内生变量/训练序列）、order=(p,d,q) 与 seasonal_order=(P,D,Q,s)；
2. 用 fit() 让模型在训练数据上完成参数估计；
3. 用 forecast(steps=...) 指定预测步数（即 horizon，预测视界）产出预测值。

首先设定预测视界，这里尝试预测未来 3 小时：

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

参数寻优提示：为 ARIMA 挑选最优参数带有主观性且耗时。可以借助 pmdarima/pyramid-arima 提供的 auto_arima() 自动搜索（仓库环境文件 environment.yaml 中已锁定 pyramid-arima==0.8.1），本课先用手工试探的方式找一组不错的参数。

先手工选一组参数拟合“首个模型”，用于查看模型概要：

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())

fit() 会打印一张结果概要表，其中包含系数估计、标准误与信息准则（如 AIC/BIC）等诊断信息，可用来对比不同参数组合的优劣。

模型评估：Walk-Forward 滚动验证

时序模型不应使用普通随机切分的交叉验证（会破坏时间顺序）。业界公认的标准做法是 walk-forward validation（前向滚动验证）：实践中，每当新数据到来，时序模型都会被重新训练一次，从而使每个时间步都给出最新、最好的预测。

过程如下：从序列起点开始，先在训练集上训练，接着预测下一个时间步；将该步预测与已知真实值比对；然后把该真实值并入训练集、重复上述过程。

效率提示：为了更高效地训练，应保持训练窗口大小固定——每次向训练集末尾加入一个新观测时，同时从集合开头移除一个最旧的观测（即“滑窗”）。

这种评估方式能更稳健地反映模型在真实部署中的表现；代价是每一步都要新建并拟合一个模型，计算成本随数据量与模型复杂度上升，小数据或简单模型下完全可以接受。

构造多步预测的标签矩阵

由于 HORIZON = 3，需要把测试集中的每个时刻同时变成它未来 1、2、3 小时真实值的“标签列”，实现方式是向前移动（shift）得到 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)

得到的 DataFrame 中每一行都是一条评估样本：

		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

数据按视界点水平位移后，load 列作为“当前已知观测”，load+1/load+2 就是 1 小时、2 小时后要验证的目标真值。

滚动预测主循环

核心评估循环：固定 30 天（720 小时）训练窗口，每步用窗口内历史重新拟合并预测未来 3 小时，随后将真实观测追加到窗口末尾、弹出最旧观测，实现滑动更新：

%%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]

可以看到每步的 3 个预测值与 3 个期望值（当前时刻及其后两小时的真值）非常接近。

整理预测结果并与真实负荷对比

把滚动预测组织成便于评估的长表（melt 后的 t+1/t+2/t+3 行），并用之前拟合好的 scaler.inverse_transform 把归一化后的预测值与真实值还原回原始负荷量纲（MW 级），方便人读与业务解释：

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

逐小时预测值与真实负荷已经相当贴合（误差在几十 MW 以内，相对负荷量级约 3000 仅为 1% 上下）。

精度度量：MAPE 平均绝对百分比误差

课程采用 MAPE（Mean Absolute Percentage Error，平均绝对百分比误差） 作为精度指标。它把预测精度表示为一个比率：先求每个时点 实际值 - 预测值 的绝对值、除以实际值得到该点的绝对百分比误差（APE），再对所有拟合点求和后除以点数 n（即取平均）：

APE = |prediction - actual| / actual
MAPE = mean(APE)

MAPE 为 10 意味着平均偏差 10%，因此数值越低越好。该指标在 common/utils.py 中与课程代码完全一致地实现为：

def mape(predictions, actuals):
    """Mean absolute percentage error"""
    return ((predictions - actuals).abs() / actuals).mean()

在 Notebook 中，先对 HORIZON > 1 的多步情形按视界分组查看各步平均绝对百分比误差：

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 %

也就是说，该模型单步预测平均只偏 0.56%，多步平均偏 1.15%，精度相当不错。可以预期越往前的步数预测越准，因为误差会随视界累积。

可视化评估结果

数值之外，用绘图能更直观地确认预测与实际的吻合程度。代码分别处理单步与多步两种绘图路径——多步情形下用红色粗线画实际负荷，用蓝色、透明度与线宽随步数衰减的曲线画 t+1、t+2、t+3 的预测，直观体现“越远期越不确定”：

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()

从图中可以直观看到预测曲线与红色真实曲线高度重合——一个精度良好的时序模型就这样诞生了。

挑战与作业

课程挑战：本课只涉及 MAPE 一种精度度量。请继续调研时序模型精度测试的其它方法（例如 MASE、RMSE、MAD、MSD 等）并逐一说明其适用场景，可参考《Forecasting: Principles and Practice》中关于预测精度评估的章节。

课后作业：为巩固所学，官方作业要求你用一份全新数据（课程提供的 Duke 大学时间序列数据集列表是一个不错的来源）构建一个全新的 ARIMA 模型，并在 Notebook 中完成数据可视化、建模、预测，用 MAPE 检验精度。完整的评分细则（Exemplary / Adequate / Needs Improvement 三档标准）参见 assignment.md。

延伸阅读

本课只触及 ARIMA 时序预测的基础。若想加深理解，可以：

• 阅读本课程时序预测章节的前置课 1-Introduction，理解时序数据的加载与绘图，以及后续课程 3-SVR 用支持向量回归做预测的对照思路；
• 结合本课配套的 sketchnotes/ml-timeseries.png 图形化笔记复习整体概念；
• 更进一步，探索更多时间序列模型类型（如 Prophet、向量自回归 VAR、深度学习时序模型等），比较它们在同一份 energy.csv 数据上的表现差异。