第9章 基准与统计点预测方法

前面的章节详细描述了开发时间序列预测所需的步骤，包括：如何生成有用的解释变量；如何训练模型；如何避免过拟合；以及如何评估模型的准确性。尚未研究的是模型本身。本章将是三章中的第一章，介绍广泛的模型及其一些特性。

本章和下一章将讨论点预测方法，然后在第11章中，将研究概率预测，它提供了处理高度不确定性数据的模型，这对于低压馈线和变电站通常是必需的（第2章）。

在点预测章节中，本章讨论传统的统计方法，而第10章将讨论有时被称为机器学习模型的内容。每种模型都有优缺点，其中一些已在5.3节中描述，但进一步的标准将在12.2节中描述。简而言之，统计模型通常更透明，更易于解释和理解。这不仅使它们有助于研究数据的核心特征，而且使它们成为良好的基准候选者。

本章介绍的大多数模型可以通过开源科学计算编程语言（如Python和R）以及流行的专有软件（如MATLAB）中的包轻松实现。然而，这些模型也可以从头开始轻松推导和训练（因为它们通常是线性函数，因此可以使用例如线性最小二乘法轻松训练，参见第8.2节），当你想要扩展模型或进行定制调整时，这可能更可取。

本章从一些简单模型开始，然后逐步引入更复杂的模型（涉及更多参数和计算成本），从指数平滑（第9.2节）、多元线性回归模型（第9.3节）、ARIMA和SARIMA模型（分别为第9.4节和第9.5节），最后到广义加性模型（第9.6节）。

在深入模型之前，有必要强调这些预测的语境：短期负荷预测（STLF）。一种常见的分类方式是根据预测时域。短期预测估计未来一天到一周（有时两周）的需求。相比之下，从1周到1年的预测称为中期负荷预测，而超过1年的称为长期负荷预测。注意，这些定义可能因情境略有不同，但通常在上述范围内。适用于STLF的模型可能不适用于中期和长期负荷预测，反之亦然。因此，这里介绍的模型是专门针对短期负荷预测选择的，它们通常高度依赖最近观测到的信息。

9.1 基准方法

本节将从考虑基本且常用的基准方法开始。如第8.1.1节所述，为任何设计良好的预测实验制定合适的基准是至关重要的。与本章通常做法一致，将考虑形如$L_{1},L_{2},\ldots$的时间序列，目标是为时间步$N+k$生成估计$\hat{L}_{N+k}$，其中$N\in\mathbb{N}$是预测原点，$k \in \mathbb {N}$是预测时域（有关这些术语的更多细节，请参见第5.2节）。

最简单的基准之一是持久性模型，可以描述为

持久性：$\hat{L}_{N+1}=L_{N}$

该模型假设下一个时间步的需求就是当前负荷。可以通过简单地重复该值来估计更远的未来时间步。

如果数据在滞后一阶处存在单一强自相关（见第6.2.2节），这可能是一个有效的模型。但对于大多数需求变化较大的应用，该方法精度很低。相反，能源需求通常具有很强的日、周和年季节性成分（第5.1节）。因此，可以对简单的持久性模型进行有效调整，产生更准确的预测模型，该模型假设当前步的行为与恰好一个季节性周期前的行为相同。这些称为季节性持久性模型，形式如下：

季节性持久性：$\hat{L}_{N+k} = L_{N+k-s_1}$

其中$\hat{L}_{N+k}$是k步超前预测，N是预测原点，$s_1$表示季节性周期（注意，假设最后一个季节性点被观测到，即出现在预测原点之前，因此$N+k-s_1 \le N$）。例如，考虑半小时数据的情况，通过将$s_1$设置为48、336或$52\times 336$，可以分别产生日、周和年季节性的季节性持久性模型。这些季节性持久性预测非常容易实现，且不需要任何训练数据。图9.1显示了日前持久性和季节性持久性模型的示例，其中季节性持久性模型使用日季节性（即今天与昨天相同）。

对于季节性数据，这些方法的一种扩展（通常是改进）是包含同一时期的多个历史观测值，并取一个简单的季节性移动平均。即：

季节性移动平均（SMA）：

$$\hat{L}_{N+k}=\frac{1}{p}\sum_{i=1}^{p}L_{N+k-is_{1}}.$$

与季节性持久性类似，通常使用一周周期（对于半小时数据为 $s_1=336$）。周简单平均通常比等效的季节性持久性模型表现更好，因为它平滑了围绕期望值的随机周际波动，从而更好地复现典型的周行为。图9.2展示了一个每日季节性数据的例子。第3天出现了异常大的需求，因此每日季节性持久性模型不如前一个图9.1示例中那么准确。而第1-3天的简单平均则降低了异常第3天的影响，从而为第4天提供了更好的估计。对于简单平均方法，需要比持久性模型稍多的训练数据，此外还需要一个验证期来选择最合适的超参数 p 的值（第8.1.3节）。然而，该模型计算非常快，实践中通常只需设置 $p=4$ 或5周即可优化模型，并能显著改进持久性模型。

9.2 指数平滑

尽管简单，第9.1节介绍的基准方法，尤其是季节性移动平均，可能出奇地准确。然而，它们的缺点之一是所使用的每个历史周被赋予相同的权重，而实际上较旧的数据与当前预测期的相关性较小。换句话说，较旧的数据对最终预测的贡献应小于较新数据。这对于负荷数据尤其相关，因为它强烈受季节性和趋势驱动。例如，几个月前（比如夏季）的数据与冬季期的相关性较小。

指数平滑方法取过去观测值的加权平均，但较旧观测值的权重衰减。为了说明这一点，考虑最简单的指数平滑形式，它生成一个平滑的1步前向输出 $\hat{L}_{N+1}$，该输出在每一步使用最新观测值 $L_N$ 按以下方式更新：

$$\begin{aligned} \hat{L}_{N+1}= \alpha L_{N} +(1-\alpha ) \hat{L}_{N} = \hat{L}_{N} +\alpha (L_N-\hat{L}_{N}), \end{aligned}\tag{9.1}$$

其中 $\alpha\in(0,1)$ 是一个平滑常数，需要在验证期内优化 α（参见第8.1.3节）。下一个观测值 $L_{N+1}$ 的估计 $\hat{L}_{N+1}$ 是当前估计 ${\hat{L}}_{N}$ 和最新观测值 $L_{N}$ 的加权平均。类似地，前一个估计也是之前观测值 $L_{N-1}$ 和更早估计 $\hat{L}_{N-1}$ 的加权平均，依此类推。换句话说，公式(9.1)可以写成展开形式：

$$\begin{aligned} \hat{L}_{N+1}= \alpha ( L_{N} +(1-\alpha )L_{N-1}+(1-\alpha )^2 L_{N-2} + \ldots (1-\alpha )^{N-1}L_2 ) +(1-\alpha )^{N}L_1, \end{aligned}\tag{9.2}$$

一个几何级数。由于 $1-\alpha\in(0,1)$，因此较旧的观测值被赋予较小的权重 α，从而对最终估计贡献较小。在特殊情况下，当 $\alpha=1$ 时，预测就是最后一个观测值，等价于第9.1节给出的简单持久性模型。该方法是一种1步前向预测，因此如果需要多步预测，则将预测值反馈回模型以代替未观测到的值。最优参数可以通过最小化1步前向预测的平方误差和（在验证期内）找到，但除了 α 外，还必须产生初始估计。初始估计可以通过对先前值取简单平均生成。由于平滑常数的嵌套应用，平方误差和是一个非线性方程，因此必须使用数值方法进行优化，而不是直接求解。

为了说明指数平滑方法，考虑图9.3给出的基本示例。使用两个不同的 $\alpha$ 值应用两个指数模型，以生成1步前向预测。使用 $\alpha=0.7$ 的模型不太平滑 α，主要由最近的点驱动。使用 $\alpha=0.2$ 的模型最平滑 α，其加权平均中来自较旧历史值的贡献更多。在这种情况下，较少平滑（较高的 $\alpha$ 值）对预测更有用 α，因为数据具有下降趋势，而较旧的点与最近数据相关性更小。

在这种基本形式下，指数平滑相对有限，因为它忽略了作为重要需求组成部分的趋势或季节性。一个更先进的考虑季节性的指数平滑算法是

Holt-Winters-Taylor (HWT) 指数平滑方法并建模两个季节周期。该方法使用以下方程组估计时刻 t 的负荷 $\hat{L}_{N+1}$：

$$\hat{L}_{N+1}=l_{N}+d_{N+1-s_{1}}+w_{N+1-s_{2}}+\phi e_{N}$$

$$e_{N+1}=\hat{L}_{N+1}-(l_{N}+d_{N+1-s_{1}}+w_{N+1-s_{2}})$$

$$l_{N+1}=l_{N}+\lambda e_{N+1}$$

$$d_{N+1}=d_{N+1-s_{1}}+\delta e_{N+1}$$

$$\begin{aligned} \hat{L}_{N+1}= & {} l_{N} + d_{N+1-s_1}+ w_{N+1-s_2} + \phi e_{N} \nonumber \\ e_{N+1}= & {} \hat{L}_{N+1} - (l_{N} + d_{N+1-s_1}+w_{N+1-s_2}) \nonumber \\ l_{N+1}= & {} l_{N} + \lambda e_{N+1} \nonumber \\ d_{N+1}= & {} d_{N+1-s_1} + \delta e_{N+1} \nonumber \\ w_{N+1}= & {} w_{N+1-s_2} + \omega e_{N+1}, \end{aligned}\tag{9.3}$$

其中参数 $\phi,\lambda,\delta,\omega$ 必须根据历史数据训练。负荷被分解为三个核心部分：水平分量 $l_t$（对应一阶相关性）以及两个季节项 $d_t$ 和 $w_t$，在负荷预测中分别对应日内和周内季节性（当然，根据数据也可以使用不同周期）。日内季节周期 $s_1$ 和周内周期 $s_2$ 是覆盖一天或一周的时间步长数，对于小时数据分别为24和168。注意每个水平项和季节项都有自己的简单指数平滑方程，如式(9.1)所示。误差项 $\epsilon_{N+1}=L_{N+1}-\hat{L}_{N+1}$ 假设服从均值为零的正态分布。在每个时间步 $N+1$ 上，使用当前水平项和季节项的值 $l_{N},d_{N},w_{N}$ 以及一阶误差项 $e_{N}$ 做出预测 $\hat{L}_{N+1}$。给定新估计值后，其他项的值可以使用各自的平滑方程（如式(9.3)所述）进行更新。由于该算法的递归性质，较旧的值对更新的贡献较小，而贡献大小由各自参数 $\phi,\lambda,\delta,$ 和 $\omega.$ 的取值决定。

φ λ δ ω 模型参数可以通过在训练数据上优化一步前向预测的平方和误差（即式(8.5)）（第8.2节）来数值求解，与之前相同。但注意，在训练参数之前，必须有水平分量和季节分量的初始估计。有几种方法可以做到这一点，但一个简单的方法是取最旧观测值的平均值，以确保有初始数据来训练算法。双重季节指数平滑模型的示例将在第14.2节的案例研究中给出。

9.3 多元线性回归

标准回归是估计单个或多个变量之间关系的统计过程。最简单且最常见的模型之一是多元线性回归，因为它易于解释、计算快速且非常通用。假设有 $n\geq 1$ 个输入变量 $X_{1,t},X_{2,t},\ldots,X_{n,t}$，这些变量与时刻 t 的负荷 $L_t$ 呈线性关系，换句话说，构建了如下预测模型：

$$\begin{aligned} \hat{L}_{N+1} = \sum _{k=1}^{n} \beta _k X_{k,N+1}. \end{aligned}\tag{9.4}$$

系数（或回归参数）$\phi _k$ 描述了每个变量在建模负荷 $L_t$ 时的解释能力（但这取决于每个变量具有相似的量级）。自变量之间假设不相关¹，一步前向预测误差 $\epsilon_{t}=L_{N+1}-\hat{L}_{N+1}$ 也假设不相关，且通常服从均值为零、方差恒定的高斯分布（见第3章式(3.5)）（方差恒定意味着误差是同方差的——见第11.6.1节）。这些假设简化了系数的训练和预测区间的建模。然而，和往常一样，通过绘制残差及其ACF图来检查这些假设是一个好主意（更多细节见第7.5节）。

如果只有一个解释变量，则模型简称为线性回归；如果有多个，则称为多元线性回归。多元线性回归通常可以写成更简洁的向量化形式：

$$\begin{aligned} \hat{L}_{N+1} = \boldsymbol{\beta }^T \textbf{X}_{N+1}, \end{aligned}\tag{9.5}$$

其中 $\mathbf{X}_{t}=(X_{1,t},X_{2,t},...,X_{n,t})^{T}$ 和 $\beta=(\beta_{1},\ldots,\beta_{n})^{T}$ 分别是自变量向量和回归参数向量。

尽管式(9.4)和(9.5)只显示了与因变量 $\hat{L}_{t}$ 在同一时间步 t 的自变量 $X_{k,t}$，但方程当然可以包含滞后时间点和自回归变量。

例如，考虑图9.4中的情况，观测值的最佳回归拟合曲线为 $Y=(X-2)^{2}+1.2=X^{2}-4X+5.4$（红色）。注意，尽管函数包含二次项 $X^{2}$，但它在系数上仍然是线性的，自变量为 $\mathbf{X}=(X^{2},X,1)^{T}$，对应的回归参数为 $\beta=(1,-4,5.4)^{T}$。因此，理解非线性关系仍然可以在线性回归框架内建模是很重要的。在需求预测的例子中，注意到第6.2.2节图6.7中需求与温度之间的非线性关系可以通过多项式（如果选择足够阶数）的线性回归来建模。

线性回归同样非常适合通过使用虚拟变量（参见第6.2.6节）来建模分类/离散变量的影响。这在负荷预测中尤为有用，因为负荷预测通常需要考虑星期几或一年中不同时间的影响。例如，一周中不同的日子往往有不同的需求特性，此时模型应该包含不同日子的影响。在多元线性回归中，这是通过引入虚拟变量 $D_{j}(k)$ 对 $j=1,\dots,7$ 来实现的（一周中的每一天各对应一个虚拟变量——星期一由 $j=1$ 表示，星期日由 $j=7\mathrm{etc.})$ 表示，这些变量指示星期几，定义由

$$D_{j}(k)=\left\{\begin{array}{ll}1,&\text{if time step k occurs on day j of the week}\\0,&\text{otherwise}\end{array}\right.$$

在线性回归模型中，我们经常使用六个虚拟变量作为输入以避免虚拟变量陷阱（见第6.2.6节），因为事实上我们可以用其他六个变量来建模某一天的影响（第七天的影响可以通过将其他六个变量设为零来建模，假设至少存在另一个项（如常数项）来确保其效应可以被建模）。

线性回归的另一个有用特性是我们能够包含交互项。这指的是我们对两个或多个变量对因变量的影响进行建模。例如，温度 $T_{k}$ 可能对需求有影响，但仅限于一天中的特定小时，比如下午2–3点。在这种情况下，我们可以包含一个温度变量项，但乘以一个表示时间段的虚拟变量，该变量在除 $2\mathrm{-}3\mathrm{pm}$ 小时之外的所有时间均为零。在线性回归中，交互项通常表示为两个项的乘积，例如 $T_{k}D_{j}(k)$ 或 $T_{k}*D_{j}(k)$。如果要建模两个以上变量的同时影响，情况类似。第14.2节的案例研究将给出线性回归模型中交互项的一个例子。

给定关于误差的假设，线性回归模型的系数通常通过最小化最小二乘估计（见第8.2节）来求得，因此训练起来非常容易且快速。回忆，由于误差被假定为具有恒定方差的高斯分布，模型的最小二乘估计也是最大似然估计，如第8.2节所示。这特别方便，因为对数似然（见式(8.8)）以及贝叶斯信息准则（BIC）和赤池信息准则（AIC）都易于计算。从第8.2.2节回顾，识别具有最小AIC或BIC值的模型是在训练数据上选择最佳模型的一种方法，它在准确性和模型复杂度之间进行权衡，有助于限制过拟合的可能性。

如第8.2.2节所述，线性模型可以很容易地适应正则化框架，如LASSO和岭回归。与AIC和BIC类似，这些技术通过在普通最小二乘回归中加入惩罚项来惩罚系数的数量和/或大小。特别是LASSO可以用作模型选择技术，因为它倾向于将不相关（或影响较小）的解释变量的系数设为零。最后，当然，与所有方法一样，模型也可以通过交叉验证来选择，找到在验证集上误差最小的模型。如果考虑的自变量很多，这可能效率很低。

给定最终训练好的模型，系数的简单线性结构提供了一种有用的方式来解释每个变量的影响（假设它们独立）。本质上，它们告诉你，在其他自变量固定的情况下，自变量每变化一个单位，因变量的期望值会变化多少。当存在交互项时，解释会变得稍微复杂，因为效应大小现在取决于其他变量的值。在这些情况下，为这些其他变量插入一系列合理值可能有助于显示效应的范围。

9.4 ARIMA和ARIMAX方法

自回归滑动平均（ARMA）技术是一种传统的线性时间序列模型，已广泛用于时间序列预测。时间序列的ARMA $(p,q)$ 模型是由以下公式描述的线性模型：

$$\begin{aligned} \hat{L}_N= \ C+ \sum _{i=1}^{p} {\psi _i}{L_{N-i}}+ \sum _{j=1}^{q} {\varphi _j}{\epsilon _{N-j}}, \end{aligned}\tag{9.6}$$

其中 $\epsilon_{t}$ 是误差项的时间序列，C是常数。ARMA模型依赖于时间序列 $L_{t}$ 是平稳的（见第5.1节），但情况未必总是如此。当序列不平稳时，可以对时间序列进行差分，直到最终序列平稳。如果应用了d次差分，可以写成：

$$\begin{aligned} L^{(d)}_{N}= L^{(d-1)}_{N} - L^{(d-1)}_{N-1}, \end{aligned}\tag{9.7}$$

其中差分迭代应用d次。当使用差分时，$\mathbf{ARMA(p,q)}$ 模型现在被称为ARIMA(p, d, q)模型（自回归积分滑动平均，积分部分表示差分），可以写成：

$$\begin{aligned} \hat{L}^{(d)}_N= \ C+ \sum _{i=1}^{p} {\psi _i}{L^{(d)}_{N-i}}+ \sum _{j=1}^{q} {\varphi _j}{\epsilon _{N-j}}+ \epsilon _N. \end{aligned}\tag{9.8}$$

写ARIMA模型的一种方便且简洁的方式是使用后移算子B（也称为滞后算子），其中B定义在时间序列的元素上：$BL_{t}=L_{t-1}$。根据定义，滞后算子因此可以写成 $B^{k}L_{t}=L_{t-k}$。于是 $\mathbf{ARMA(p,q)}$ 模型可以写成如下形式：

$$\begin{aligned} \left( 1-\sum _{i=1}^{p}\psi _{i}B^{i}\right) L_{N}=\left( 1+\sum _{j=1}^{q}\varphi _{j}B^{i}\right) \epsilon _{N} + C, \end{aligned}\tag{9.9}$$

并且ARIMA(p, d, q)模型可以写成：

$$\begin{aligned} \left( 1-\sum _{i=1}^{p}\psi _{i}B^{i}\right) (1-B)^d L_{N}=\left( 1+\sum _{j=1}^{q}\varphi _{j}B^{i}\right) \epsilon _{N} + C, \end{aligned}\tag{9.10}$$

其中 $(1-B)^{d}L_{N}$ 是d阶差分。

ARIMA模型非常通用，能够估计广泛的时间序列。它由三个主要组成部分构成：差分项$d,$、自回归项$\left(\operatorname{AR}(p)\right)$（阶数为$\sum_{i=1}^{p}\psi_{i}L_{t-i}^{(d)}$）、以及移动平均项$p_{:}$（阶数为$MA(q)$，使用历史白噪声误差项$\sum_{j=1}^{q}\varphi_{j}\epsilon_{t-j}$ψ）。与多元线性回归类似，误差项通常假设服从高斯分布（尽管也可以使用其他分布），均值为零且彼此不相关。

自回归模型AR(p)和移动平均模型MA(q)是ARMA模型的特例（实际上是$\mathbf{ARMA}(\mathbf{p},0)$和$\mathbf{ARMA}(0,\mathbf{q})$模型），在讨论完整的ARMA模型之前，值得更详细地考虑它们。自回归模型是ARMA过程，但$\varphi_{j}=0$对所有$j$成立，因此这些模型简单，其中时间序列的$p$ϕ过去值影响当前值。AR(p)可以写成

$$\begin{aligned} \hat{L}_N= \ C+ \sum _{i=1}^{p} {\psi _i}{L_{N-i}}. \end{aligned}\tag{9.11}$$

回顾第3.5节，偏自相关是时间序列与其滞后值之间的自相关度量，但去除了中间滞后的影响。这意味着PACF是识别AR过程阶数p的自然方式，因为对于滞后$k>p$，PACF应为零。在实践中，可以通过观察可用时间序列的样本PACF图，并识别滞后值何时实际上为零（即持续落在通常与PACF一起绘制的95%置信区间内）来检测阶数。此外，对于AR(p)过程，ACF也应指数衰减至零。简单AR(4)模型的示例如图9.5所示。注意，在PACF中，超过滞后4之后没有超出置信区间的相关性，符合预期。

相比之下，移动平均模型受误差的过去值影响，因此过去的大偏差可能影响当前时间序列值。纯MA(q)过程的一个有用性质是，自相关函数从滞后$q+1$开始应为零。因此，可以利用ACF图来识别MA时间序列及其阶数。需要注意的是，尽管ACF和PACF可用于识别AR和MA模型及其阶数，但在实践中使用的是这些函数的样本版本，应用于实际观测数据，因此结果可能偏离更清晰的理论解。换句话说，自相关可能超出置信区间，但这可能是虚假的，仅出于随机偶然性。

现在考虑ARIMA模型，必须找到最优阶数p、q和d，以训练最终模型的系数$\psi_{i},\varphi_{j}$。通常，这是通过比较不同p、q和d选择下的AIC值（见第8.2.2节）来完成的。比较所有可能值是不可行的，因此通常仅在阶数的良好近似附近进行搜索。为历史时间序列数据找到ARIMA模型最佳阶数的常用方法是Box-Jenkins方法，该方法利用ACF和PACF来识别自回归和移动平均阶数，如上所述。该过程通常包括以下步骤：

1. 检查时间序列是否平稳。如果不是，则进行差分直到最终序列平稳。平稳性可以通过多种方式检查。除了时间序列图外，非平稳时间序列的另一个指示是自相关函数随滞后缓慢衰减（见第3章）。然而，也有平稳性检验方法，如附录A所述。
2. 识别自回归（AR）项和移动平均（MA）项的阶数。这可以通过检查自相关函数（ACF）和偏自相关函数（PACF）图（见第3章和第6.2.4节）来估计。特别地，如果模型有阶数为p的AR分量，则PACF应从滞后$p+1$及以上开始实际上为零。类似地，对于阶数为q的MA模型，ACF应从滞后$q+1$及以上开始实际上为零。在实践中，可以通过查看相应的图并考虑它们是否超过95%置信区间（通常包含在图中，见第6.2.4节）来找到这些阶数。
3. 使用ACF和PACF作为正确阶数的近似，检查一组具有不同p、d、q值（在近似值附近）的ARIMA模型的AIC（或BIC）值。最终阶数是那些给出最小AIC（BIC）值的阶数。

需要强调的是，ACF和PACF通常不能给出关于正确阶数的明确答案，因此在实践中，它们用于近似正确阶数，然后在步骤3中使用AIC/BIC进行测试。

这里通过一个具体示例说明Box-Jenkins方法，该示例使用ARIMA(3,0,1)模型（等价于ARMA(3,1)），由$y_t = 0.14+0.609y_{t-1}-0.5y_{t-2}+0.214y_{t-3}+0.624e_{t-1}+e_t$给出。时间序列如图9.6所示，使用Matlab的simulate函数生成。$e_{t}$是误差序列，服从标准正态分布。在这种情况下，序列是平稳的，因此无需差分。为了检查自回归和移动平均阶数，考虑ACF和PACF图，如图9.7所示，同时给出了95%显著性水平的置信区间。ACF（上图）指示MA阶数，显示最大相关性在滞后1处，符合预期，但在滞后16和17处也存在显著相关性（明显超出置信区间）。注意，ACF没有随滞后逐渐减小，这支持了时间序列平稳的结论。PACF指示AR阶数，在此示例中显示在滞后1–4处有显著峰值，表明阶数略大于预期。此外，在较大滞后处也有较小的峰值超出置信区间。这一分析表明，ACF和PACF分析在给出精确阶数的完整答案方面存在局限性。事实上，这些图有局限性，因为预计会有5%的自相关由于随机偶然性而超出置信区间。这意味着必须谨慎解释ACF和PACF，并结合AIC使用。

表9.1文中ARIMA示例的不同AR(p)和MA(q)值的Akaike信息准则结果

利用相关性分析有助于定位正确阶次的大致范围。在本示例中，ACF和PACF提示阶次约为$q = 1$和$p = 4$，应当针对接近这些值的自回归和移动平均阶次的各种组合进行AIC检验。由于$ARIMA(p,0,q)$模型的参数个数为$p + q + 1$（其中一个是常数项），因此Akaike信息准则（AIC）具有特别简单的形式

$$\begin{aligned} AIC = 2(p+q+1) - 2\ln (L), \end{aligned}\tag{9.12}$$

其中L是ARIMA模型的似然函数。检查所有阶次组合的AIC，其中$p=1,2,3,4$、4和$q=1,2,3,4$。每个p和q组合的结果如表9.1所示，该表显示$p=3$和$q=1$时达到最小AIC值50.45，正确识别了ARIMA(3,0,1)模型。

应当注意，任何MA模型都可以通过具有足够多滞后期数（p值）的AR模型来估计。由于AR模型的系数计算速度远快于完整的ARIMA模型，因此更推荐将任何ARIMA模型替换为阶数足够大的AR模型（若非平稳则进行差分）。这还可以简化模型的分析与解释。然而，这可能需要相对较大的阶数，从而相比简单的ARMA模型引入更多参数，降低了简约性和可解释性。

ARIMA模型有许多有用的扩展。对于负荷预测而言，最重要的扩展之一是将模型扩展到包含其他解释变量。该模型被称为带解释变量的自回归积分移动平均（ARIMAX）模型。ARIMAX $(p,d,q)$模型包含额外外部变量，由式(9.13)描述

$$\begin{aligned} \begin{aligned} \hat{L}^{(d)}_N= C+{\sum _{i=0}^{h} {\mu _i}{X_{N-i}}}+ \sum _{i=1}^{p} {\psi _i}L^{(d)}_{N-i}+{ \sum _{i=1}^{q} {\varphi _i}{\epsilon _{N-i}}} \end{aligned} \end{aligned}\tag{9.13}$$

差分与(9.7)中相同。$\sum_{i=0}^{h}\mu_{i}X_{t-i}$是解释变量项。分析该模型时，可先考虑不含外生输入的ARIMA模型以确定方程阶次，然后拟合包含外生变量的完整模型。注意，在最终确定AR和MA阶次时，应对完整的ARIMAX方程应用AIC。

9.5 SARIMA与SARIMAX模型

ARIMA模型的一个重要扩展是包含季节性。季节性ARIMA（SARIMA）包含一组额外的超参数，记为P、D、Q，将模型扩展为在指定季节性水平上包含自回归、差分和移动平均项。这些模型写作$\mathrm{ARIMA}(p,d,q)(P,D,Q)_{S}$，其中S表示季节性。例如，对于具有日季节性的小时数据，SARIMA模型可写作$\mathrm{ARIMA}(p,d,q)(P,D,Q)_{24}$。对于季节性ARIMA模型，ACF和PACF的解释有所不同。考虑一个简单情形：$d=D=q=Q=0$但$p=2$和$P=1$。这意味着时间序列将在滞后1、2、24、25、26处有自回归项。注意，$p$和P项的组合意味着季节内滞后（1,2）应用于季节性滞后24上。$\mathrm{ARIMA}(2,0,0)(1,0,0)_{10}$的一个示例如图9.8所示，同时给出了偏自相关函数。注意在周期性间隔10、20以及11、12、21、22处的PACF有显著尖峰。

后移算子对于表示SARIMA模型特别有用。例如，对于小时季节性数据，$\mathrm{ARIMA}(p,d,q)(P,D,Q)_{24}$可以表示为（为清晰起见，此处未包含常数项）

$$\begin{aligned} \left( 1-\sum _{i=1}^{p}\psi _{i}B^{i}\right) \left( 1-\sum _{j=1}^{P}\zeta _{j}B^{24j}\right) (1-B)^d (1-B^{24})^D L_{N}= \nonumber \\ \left( 1+\sum _{i=1}^{q}\varphi _{i}B^{i}\right) \left( 1+\sum _{j=1}^{Q}\theta {j}B^{24j}\right) \epsilon _{N}, \end{aligned}\tag{9.14}$$

其中 $\psi_{i}$ 是非季节性AR分量的系数，$\zeta_{j}$ 是季节性AR分量的系数，$\varphi_{i}$ ζ是非季节性MA分量的系数，$\theta {j}$ ϕ是季节性MA分量的系数。注意 $(1-B^{24})$ θ) 代表季节性差分，即 $(1-B^{24})L_N = L_N - L_{N-24}$ 的季节性差分通常 $D=1$ 就够了。

关于ARIMA和SARIMA模型的更多细节，请参阅文献[2]以及附录D中列出的其他文献。

9.6 广义加性模型

第9.3节中指定的线性模型存在各种局限性。两个最强且最常见的假设是误差服从高斯分布，并且模型是各个输入变量的简单线性组合。

广义线性模型（GLM）是对简单多元线性模型的扩展，包含一个连接函数，可以允许更多样化的关系类型。使用第9.3节中的符号，如果对于 $n\geq 1$ 个输入变量 $X_{1,t},X_{2,t},\dotsc,X_{n,t}$ ，时间t的因变量 $L_{t}$ 服从GLM，则有

$$\begin{aligned} g(\mathbb {E}(\hat{L}_{N+1})) = \sum _{k=1}^{n} \beta _k X_{k,N+1}, \end{aligned}\tag{9.15}$$

对于某个（可能非线性的）连接函数 $g(.)$ ，并且响应变量 g 来自指数族中的概率分布（例如高斯分布、二项分布或伽马分布——参见第3.1节）。换句话说，对于GLM，通过 $g)$ 对因变量期望值的变换是一个线性 g 模型。注意，与线性模型一样，所有线性系数 $\beta_{k}\mathbf{\dot{s}}$ 必须被估计，但此外，还必须选择连接函数和误差的概率分布模型。当连接函数是恒等函数 $(g(x)=x)$ ，并且因变量被假定为高斯分布时，式(9.15)退化为第9.3节中介绍的简单多元线性回归模型。所需的连接函数和分布的选择取决于所考虑的问题。例如，如果因变量是非负的，则对数连接函数可能是有效的。

本文不研究一般的GLM。相反，重点是一种非常具体且强大的GLM形式，称为广义加性模型（GAM），它在负荷预测中非常成功。4 GAM的一般形式为

$$\begin{aligned} g(\mathbb {E}(\hat{L}_{N+1})) = \sum _{k=1}^{n} f_k(X_{k,N+1}), \end{aligned}\tag{9.16}$$

对于某些（可能非线性的）光滑函数 $f_{k}$ 。

GAM相比GLM有几个优势，首先，函数 $f_{k}$ 允许对更多样化的、可能非线性的关系进行建模，而GLM只有 $f_{k}(X_{k},N+1)=\beta X_{k,N+1}$ 形式。此外，这些函数 β 通常是非参数建模的，而GLM通常假设参数变换和分布（GAM也可以利用常见的参数形式，例如对数函数或每个 $f_{k})$ 的多项式）。注意，GAM仍然使用连接函数 $g$ ，它可用于将因变量变换为更适于训练的 g 形式。

对加性模型中的每个函数 $(f_{k})$ 采用非参数方法，使算法能够从观测数据中学习每个输入变量 $X_{k,N+1}$ 之间的关系。最常用的方法是使用基函数对每个函数进行建模（参见第6.2.5节）。因此，每个函数 $f_k$ 被建模为

$$\begin{aligned} f_k(X_{k, N+1}) = \sum _{i=1}^m \alpha _{k, i} \phi _{k, i}(X_{k, N+1}), \end{aligned}\tag{9.17}$$

对于基函数 $\phi_{k,i}(X)$ 。注意，这种形式将GAM (9.16) 转换为 φGLM，因为加性函数的和现在是基中线性函数的和。

对于GAM，通常选择样条作为这些基函数。样条是一种分段连续函数，由其他更简单的多项式函数组成。样条最简单的例子之一是分段线性组合。线性和三次样条的示例见图9.9。注意，由于样条是连续的，一个多项式的末端必须与下一个多项式的开始连接。节点指定了多项式相互连接的位置。三次版本对节点间的观测值（红点）进行回归，以确定每个三次多项式中的另外两个系数（其中两个系数已经由插值约束得到）。

更精确地说，考虑一维情况，目标是近似定义在区间 $[a,b]\subset\mathbb{R}$ 上的函数 $f:[a,b]\longrightarrow\mathbb{R}$。对于在 $a=z_{1}<z_{2}<\cdot\cdot\cdot<z_{m-1}<z_{m}=b$ 处的 m 个节点，在每个子区间 $[z_{i},z_{i+1}]$ 上通过多项式 $s_{i}(z)$ 将样条拟合到某些数据上。此外，$s_{i}(z_{i+1})=s_{i+1}(z_{i+1})$ 因为样条在节点处应该是连续的。

可以对样条施加其他约束，以使其更易于训练或满足其他标准。GAM最常见的要求之一是确保样条具有特定的光滑度。可以看出，图9.9中的三次插值在节点之间是光滑的，但在节点本身处不光滑。在跨节点插值时，约束三次样条光滑意味着所有系数都可以唯一确定。另一种表示样条光滑的方式是，导数（到足够阶数）在节点处连续。图9.10展示了一个跨节点光滑的三次样条示例。

注意，预测的目标是对数据进行回归，因此没有必要（也不希望）严格地插值通过观测值。然而，原理仍然是相同的，最终的样条应该在包括节点在内的整个区间上连续。这是通过将关系的基函数版本对观测值进行回归来实现的（参见上文式(9.17)）。

某些基函数，如B样条，具有非常理想的性质，例如在节点处提供光滑性。此外，尽管基函数的数量和类型应足够灵活以适应数据，但如果没有额外的约束或正则化（第8.2.4节），大量的节点和高的多项式次数会增加对噪声过拟合的机会。此外，这将意味着多项式会变得非常“弯曲”。为了防止这种情况，一种方法是包含一个额外的项，通常用于惩罚最终解中光滑性的缺乏。回想一下，这很像LASSO（第8.2.4节）方法和其他用于防止过拟合的正则化技术。

一个简单的例子是链接函数为恒等函数的情况。由于GAM在基函数上是线性的，因此可以考虑对N个观测到的因变量值 $\mathbf{L}=(L_{1},\ldots,L_{N})^{T}$ 进行最小二乘拟合（见第8.2.4节）。换句话说，目标是极小化

$$\begin{aligned} \sum _{l=1}^N\left( L_l - \sum _{k=1}^n \sum _{i=1}^m \alpha _{k, i} \phi _{k, i}(X_{k, l}) \right) ^2 \end{aligned}\tag{9.18}$$

通过训练参数 $\alpha_{k,i}$ 对于 $k=1, \ldots , n$、n 和 $i=1,\ldots,m$。对于大量的基函数，该模型很可能会过拟合数据。为了防止这一点，可以应用惩罚，即

$$\begin{aligned} \left( \sum _{l=1}^N\left( L_l - \sum _{k=1}^n \sum _{i=1}^m \alpha _{k, i} \phi _{k, i}(X_{k, l}) \right) ^2\right) + K(f_1, \ldots , f_n). \end{aligned}\tag{9.19}$$

函数 $K(f_1, \ldots , f_n)$ 是基于各个函数 $f_k$ 的惩罚。为了惩罚对光滑性的偏离，通常考虑如下惩罚：

$$\begin{aligned} K(f_1, \ldots , f_n)= \sum _{k=1}^n\lambda _{k}\int f_{k}^{\prime \prime }(x_k)^{2}dx_k, \end{aligned}\tag{9.20}$$

其中每个变量的惩罚大小由光滑参数 $\lambda _k$ 控制。极小化每个函数的二阶导数之和会减少函数的弯曲程度，即根据lambda的值鼓励更多的光滑性。通常，这些光滑参数像往常一样通过交叉验证（见第8.1.3节）或优化信息准则（第8.2.2节）来确定。

(a) 工作日项

(b) 小时与温度的交互作用。

由于(9.17)中的基函数表示，可以证明惩罚项具有特别方便的二次型形式

$$\begin{aligned} \int f_{k}^{\prime \prime }(x_k)^{2}dx_k = \boldsymbol{\alpha }_k ^T \textbf{S}_k \boldsymbol{\alpha }_k, \end{aligned}\tag{9.21}$$

其中$\boldsymbol{\alpha }_k = (\alpha _{k,1}, \ldots , \alpha _{k, m})^T$和$\textbf{S}_k \in \mathbb {R}^{m \times m}$是由基函数在输入值$X_{k,l}$处的导数αα构成的矩阵。

与多元线性回归模型类似，GLM和GAM可用于建模两个或多个特征的交互作用，例如

$$\begin{aligned} g(\mathbb {E}(\hat{L}_{N+1})) = f_1(X_{1,N+1})+ f_2(X_{2,N+1})+ f_3(X_{1,N+1}, X_{3,N+1}) \end{aligned}\tag{9.22}$$

在这种情况下，前两个函数$f_1(X_{1,N+1}), f_2(X_{2,N+1})$分别建模单个变量，但第三个函数建模$X_{1,N+1},X_{3,N+1}$交互作用的效果。在这些情况下，可以使用多维版本的样条函数。

GAM模型的加性特性使得模型可解释，因为即使使用复杂的非线性函数，也可以分析和可视化单个特征及交互作用的贡献。图9.11a和b展示了各术语对最终预测贡献的示例性可视化。图9.11a显示了工作日$(W_{k})$对需求的贡献，表明对于该特定模型，周末负荷较低，周四最高。图9.11b显示了小时$\left(\boldsymbol{H_{k}}\right)$和室外温度$(T_{k}^{out})$交互作用的综合效果，例如，夜间和低温时影响最小，中午高温时影响最大。绘制完整输入变量的较小子集（通常为一个或两个）的图称为部分依赖图，使我们能够检查并更好地解释不同分量的整体效果。

有许多不同的方法和参数可供选择，以及许多GAM编程包，例如R中的gam或mgcv，以及Python中的pygam，这些包适用于一系列样条、平滑参数选择方法和链接函数。通常这些包有自己的默认设置，但在许多情况下可以调整这些设置以确保更准确的拟合和更好的性能。特别是如果已知误差不是高斯分布，或者某个自变量与因变量仅呈线性关系，则可以在实现时进行指定。还应检查其他参数或数据假设，但如果不确定，可以通过交叉验证方法检查多个值。由于大多数包都使用了正则化，因此最好指定比实际需要更多由样条决定的自由度。通常，残差检查（第7.5节）可用于评估最终模型并识别错误假设或改进领域。

注意，可能对基/样条函数施加额外的约束，以更好地对需求数据中的特征进行建模。特别是，由于许多因变量（例如小时或星期几）经常具有周期性，可以选择基函数来包含这些特征，例如周期性B样条，这些可在上述某些包中使用。

以上是对GAM的基本介绍，对于这个非常复杂的领域，更详细的描述超出了本书的范围。附录D.2中提供了一些进一步阅读资料。

9.7 问题

对于需要使用实际需求数据的问题，请尝试使用附录D.4中列出的一些数据。最好选择至少包含一年小时或半小时数据的数据。在所有使用该数据的情况下，将其按3:1:1的比例划分为训练集、验证集和测试集（第8.1.3节）。

1. 选择一个需求时间序列。分析其季节性（见第6.2节）。为测试集生成一些简单的基准预测，包括持久性预测和季节性持久性预测（针对你发现的每个季节分别生成一个）。计算RMSE误差。哪个更低？这与你观察到的季节性趋势相比如何？将这些结果与时间序列的ACF和PACF图进行比较。
2. 延续上一节的实验，使用识别出的季节性生成季节性移动平均。使用验证集（第8.1.3节）确定在平均中应包含的季节项p的最优值。如果存在多个季节性，哪一个的整体误差最小？在测试集上，最优季节性平均预测的RMSE误差与前一个问题中的持久性预测相比如何？
3. 生成一个简单的1步前向指数平滑预测（第9.2节），用于负荷预测时间序列（最好是具有双重季节性模式的数据，通常是日和周）。手动选择不同的平滑参数值。绘制RMSE误差与平滑参数的关系图。执行网格搜索以找到最优平滑参数α（第8.2.3节）。最优预测与简单的持久性预测相比如何？现在考虑Holt-Winters-Taylor预测，并对四个参数 ， ， ， 执行网格搜索。

φ λ δ ω4. 研究线性模型的LASSO拟合。设置一个包含少量正弦项的模型的系数，例如$\sum_{k=0}^{N}\alpha_{k}$ sin kx，其中N约为5，以及$x\in[0$ , 4 ]。从这些数据中采样α π20个点（并添加少量高斯噪声）。现在使用最小二乘回归拟合形如$\sum_{k=0}^{50}\gamma_{k}$ sin kx的多元线性方程，以找到系数γ。然后在20个新的x ∈ [0, 4 ]值上绘制训练好的模型。γ π拟合是否良好？现在尝试使用不同正则化参数λ值最小化LASSO函数（见第8.2.4节）。随着参数λ的变化，拟合如何变化？有多少个系数为零（或非常小）？使用内置函数进行LASSO拟合，例如Python中的sklearn7或R中的glmnet8。

1. 证明对于GAM的基表示，二阶罚项（9.20）的形式为$\alpha_{k}^{T}\mathbf{S}_{k}\alpha_{k}$。
2. 尝试生成一个拟合需求曲线的线性模型。考虑使用哪些特征，如果一天中的时间很重要，考虑使用虚拟变量。如果天气数据可用，检查其与需求的关系（见第14章）。在第14.2节的案例研究中，将生成一个线性模型用于模拟低压需求。在阅读到该部分后回到这个问题，看看有哪些相似之处。你做了什么不同的事情？你想在模型中改变什么？分别使用Python和R中的标准包如sklearn9和lm10拟合线性模型。现在使用相同的特征实现一个GAM。同样，这些预测可以使用标准Python和R包如pygam1和mgcv12进行训练。这些包通常具有与线性模型相似的语法。现在比较预测和误差。对于GAM，查看部分依赖图。每个所选变量的关系是什么？

参考文献

1. D. Ruppert, S.M. David, Statistics and Data Analysis for Financial Engineering: With R Examples. Springer Texts in Statistics (2015)

1. R.J. Hyndman, G. Athanasopoulos，《预测：原理与实践》，第二版（OTexts，澳大利亚墨尔本，2018）。https://Otexts.com/fpp2。访问于2020年7月。

1. P. Gaillard, Y. Goude, R. Nedellec, Additive models and robust aggregation for gefcom2014 probabilistic electric load and electricity price forecasting. Int. J. Forecast. 32(3), 1038–1050 (2016)

1. M. Voss, J.F. Heinekamp, S. Krutzsch, F. Sick, S. Albayrak, K. Strunz, Generalized additive modeling of building inertia thermal energy storage for integration into smart grid control. IEEE Access 9, 71699–71711 (2021)

除非在材料的致谢行中另有说明，本章中的图像或其他第三方材料均包含在本章的知识共享许可中。如果材料未包含在本章的知识共享许可中，且您的预期使用不被法定法规允许或超出允许使用范围，您将需要直接获得版权持有人的许可。