ML-For-Beginners 时间序列实战：用 statsmodels 构建并评估 ARIMA/SARIMA 电力负荷预测模型

本篇指南基于 ML-For-Beginners 课程第 7 周「TimeSeries」的第 2 课（7-TimeSeries/2-ARIMA/README.md），讲解如何用 ARIMA（自回归积分滑动平均）类模型对非平稳时间序列建模。以 GEFCom2014 电力负荷数据集为例，你将完整掌握从数据加载、训练/测试集切分、MinMax 缩放，到 SARIMAX 参数配置、walk-forward（滚动前推）验证和 MAPE 精度评估的全流程，最终得到一个单步预测误差约 0.56%、多步预测误差约 1.15% 的负荷预测模型。

一、核心概念：平稳性与差分

ARIMA 模型特别适用于拟合表现出**非平稳性（non-stationarity）**的数据。在动手之前，先明确两个统计概念：

• 平稳性（Stationarity）：从统计角度看，平稳数据的分布在时间平移后不发生改变。非平稳数据则表现出由趋势引起的波动，必须先做变换才能分析。例如季节性（seasonality）会给数据引入周期性波动，可以通过「季节性差分（seasonal-differencing）」过程来消除。
• 差分（Differencing）：将非平稳数据转换为平稳数据的过程，核心是去除其不恒定的趋势。差分「消除了时间序列水平上的变化，去掉趋势和季节性，从而稳定了时间序列的均值」。

这两个概念直接对应后面模型参数中的 d（差分次数）与 D（季节性差分次数），理解它们是理解 ARIMA 参数体系的前提。

二、拆解 ARIMA：AR、I、MA 三部分各做什么

ARIMA 是 AutoRegressive Integrated Moving Average（自回归积分滑动平均）的缩写。把它拆开看，每一部分对应时间序列的一个建模侧面：

组成	全称	作用
AR	AutoRegressive（自回归）	向「过去」看：分析数据中的历史值（称为 lags，滞后值），并将「关注的演化变量回归于它自身的滞后值」来建模
I	Integrated（积分）	指数据经过差分步数处理以消除非平稳性的过程，与 ARMA 模型的区别就在这里
MA	Moving Average（滑动平均）	输出变量由观测当前及过去各滞后值来确定

一句话总结：ARIMA 的作用是让模型尽可能贴合时间序列这种特殊形式的数据，从而做出预测。

当数据有季节性：SARIMA 与 P、D、Q 参数

本课程的电力负荷数据具有明显的季节性（白天/夜间用电模式），因此实际使用的是季节性 ARIMA 模型（SARIMA），需要额外引入一组季节参数：

• p：自回归部分参数，引入过去值；
• d：积分部分参数，决定对时间序列施加多少差分；
• q：滑动平均部分参数；
• P、D、Q：与 p、d、q 关联相同，但对应模型的季节性成分。

三、数据集与环境准备

数据集：GEFCom2014 电力负荷

本课程使用的数据来自 GEFCom2014 预测竞赛，包含 2012 年至 2014 年共 3 年的小时级电力负荷（load）与温度（temp）数值。任务是根据历史负荷模式预测未来的电力负荷。

数据文件为 7-TimeSeries/2-ARIMA/working/data/energy.csv，包含 timestamp,load,temp 三列，从 2012-01-01 00:00:00 到 2014-12-31 23:00:00 连续的小时记录（约 26,304 条），前几行为：

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

该 CSV 由配套脚本 extract_data.py 从 GEFCom2014 竞赛原始压缩包中提取生成：脚本读取 Excel 数据，用 Date + (Hour - 1) 构造 timestamp，把 T 列重命名为 temp，并丢弃无负荷数据的时段，最终保留 timestamp、load、temp 三列。

环境依赖

working/common/environment.yaml 定义了课程配套的 conda 环境 dlts，与本课直接相关的关键依赖包括：statsmodels==0.9.0（ARIMA 模型实现）、pandas、numpy、scikit-learn==0.20.3（MinMaxScaler）、matplotlib，以及通过 pip 安装的 pyramid-arima（auto_arima() 自动选参）。创建环境：

conda env create -f 7-TimeSeries/2-ARIMA/working/common/environment.yaml

该环境文件锁定了较旧的 Python 3.6.6 版本，是课程当年的配置；在新环境中运行时需自行调整版本，但 statsmodels 的 SARIMAX API 与本文代码保持一致。

公共工具函数：load_data 与 mape

数据加载与评估函数封装在 7-TimeSeries/common/utils.py 中（本课 working 目录下有一份同名拷贝）：

• load_data(data_dir)：读取 energy.csv，把 timestamp 解析为日期并设为索引，再用 pd.date_range(min, max, freq='H') 按小时频率重索引（reindex）整个索引区间，这样如果数据缺失时段会显式体现出来（本数据集恰好没有缺失）；
• mape(predictions, actuals)：实现 MAPE 公式 (np.absolute(predictions - actuals) / actuals).mean()，即逐点绝对百分比误差的均值。

四、加载数据与划分训练/测试集

在 working/notebook.ipynb 中，先安装并加载依赖：

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

加载全部 2012 年 1 月至 2014 年 12 月的负荷数据并可视化：

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

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

按时间切分，而非随机切分

时间序列的训练/测试切分有一个铁律：测试集必须覆盖晚于训练集的时间段，否则模型会从「未来」泄漏信息。本课的切分方式：

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

即训练集取 2014-11-01 至 2014-12-29（按小时共 1416 行），测试集取 2014-12-30 至 2014-12-31（48 行）。由于该数据是日消费模式很强、且近期消耗最接近当前消耗的序列，使用相对较短的时间窗口训练就足够了。

注意：因为拟合 ARIMA 模型的函数（SARIMAX.fit()）在拟合过程中使用样本内（in-sample）验证，所以本课省略了单独的验证集。

五、数据准备：过滤与 0–1 缩放

训练前需要两步预处理：按集合过滤时间区间与列，以及缩放，确保数据落在 0–1 区间内。

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)

然后用 MinMaxScaler 缩放。关键细节：fit_transform 只在训练集上拟合，测试集只能 transform，避免数据泄漏：

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

# 之后再用训练集拟合出的 scaler 变换测试集
test['load'] = scaler.transform(test)

可以通过直方图对比缩放前后的分布形态（原数据集中在 2000–4000 区间，缩放后压到 0–1）：

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

六、实现 ARIMA：SARIMAX 建模三步法

使用 statsmodels 实现 ARIMA 遵循三个步骤：

1. 调用 SARIMAX() 并传入模型参数 p, d, q 与 P, D, Q，定义模型；
2. 调用 fit() 用训练数据拟合模型；
3. 调用 forecast() 指定预测步数（horizon，预测视野）做出预测。

先设定预测视野。本课试 3 小时：

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

手动选取一组参数来尝试：

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

seasonal_order 的第四元 24 表示季节性周期为 24（小时），与负荷数据的日周期吻合；d=1 表示对序列做一次差分以消除趋势，D=1 表示再做一次季节性差分。拟合后会打印完整的回归结果表（系数、标准误、AIC/BIC 等信息）。

🎓 课程建议：为 ARIMA 挑选最优参数是主观且耗时的过程，可以考虑使用 pyramid 库（pmdarima）的 auto_arima() 函数 做自动网格搜索（注意该库现名已改为 pmdarima）。

七、模型评估：Walk-Forward 滚动前推验证

评估时间序列模型的推荐做法是 walk-forward 验证：模拟实践中「每来一个新数据点就重新训练一次模型」的流程——从时间序列开头出发，在训练集上训练，对下一个时间步做预测并与真实值比较；随后把训练集扩展纳入该真实值，重复上述过程。这种方式能更稳健地估计模型在真实场景下的表现，代价是要训练大量模型——数据量小或模型简单时可接受，大规模时则需注意计算成本。

为提升训练效率，应保持训练窗口大小固定：每向训练集尾部追加一个新观测，就从头部移除一个观测（滑动窗口）。

步骤 1：按视野构造测试目标

为每个 HORIZON 步构造「t+1 … t+H」的未来目标列：

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

步骤 2：滚动窗口循环预测

用长度为测试集大小的循环执行滑动窗口预测，固定 720 小时（30 天）的训练窗口，每轮重新拟合 SARIMAX：

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

注意此处使用的阶数从演示时的 (4, 1, 0) 调整为 (2, 1, 0)。运行输出可实时看到训练过程：

2014-12-30 00:00:00
1 : predicted = [0.32 0.29 0.28] expected = [0.33, 0.29, 0.27]

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

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

步骤 3：把预测与真实值合并并反变换回原尺度

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

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

预测值与真实负荷逐小时贴合，误差都在几十 MW 量级。

八、精度检查：MAPE 计算与可视化

用**平均绝对百分比误差（MAPE，Mean Absolute Percentage Error）**检验模型精度。其定义为：

即对每个预测点计算 |actual_t - predicted_t| / actual_t，再对所有 n 个预测点取平均。MAPE 以百分比直观表达预测偏差——例如 MAPE 为 10 意味着预测平均偏离 10%。数值越低越好。

代码实现分两步。先看每个视野步的逐点误差均值：

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

再分别计算单步与多步 MAPE（mape() 即 utils.py 中 mean(|pred - actual| / actual) 的实现）：

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 %

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

Multi-step forecast MAPE:  1.1460048657704118 %

单步预测误差仅约 0.56%，多步（3 小时视野）约 1.15%。最后用图直观呈现：红色粗线为真实负荷，蓝色线为各视野步的预测（随步数增大逐渐变细、变透明，表达不确定性递增）：

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，你还可用 MAE、MSE、RMSE、sMAPE 等指标，研究并注释它们的适用场景与局限。
• 课后作业：用新数据构建一个 ARIMA 模型——换一个数据集（可参考杜克大学的时间序列数据集合集），在 notebook 中记录过程、可视化数据与模型、并用 MAPE 测试精度。评分标准要求呈现一个已构建、已测试、已解释的新 ARIMA 模型 notebook，含可视化与精度说明。
• 配套材料：本课的完整可运行 notebook 见 7-TimeSeries/2-ARIMA/working/notebook.ipynb，参考答案见 7-TimeSeries/2-ARIMA/solution/notebook.ipynb（另提供 R 与 Julia 语言版本的 solution）。
• 下一步：本模块第 3 课将改用支持向量回归（SVR）做时间序列预测，对比不同经典模型在同一负荷数据上的表现；也可以深入微软 forecasting 仓库了解其他时间序列模型类型。

关键要点回顾

环节	关键点
数据	GEFCom2014 小时级负荷，load_data() 按 freq='H' 重索引保证无缺口
切分	按时间先后切分（2014-11-01 / 2014-12-30），杜绝未来信息泄漏
预处理	MinMaxScaler 仅在训练集 fit，测试集只 transform
建模	SARIMAX(endog, order=(p,d,q), seasonal_order=(P,D,Q,m))，m=24 对应日周期
评估	Walk-forward 滑动窗口（720 小时）+ forecast(steps=HORIZON) + MAPE
结果	单步 MAPE ≈ 0.56%，多步 MAPE ≈ 1.15%