第11章 概率预测方法

前两章讨论的是点预测，它只对预测范围内的每个时间步生成一个单一估计值，即每个时间段 $t=N+1,N+2,\ldots,N+k$ 对应一个值 $L_{t}$（假设预测范围为从预测原点 N 开始的 k 步向前）。点估计在描述未来需求方面存在局限，特别是当原始数据具有较大不确定性时。通过估计预测范围内每个时间段的需求分布，可以更细致地描述需求可能的值。估计分布范围的预测通常称为概率预测。这正是本章的主题。

11.1 概率预测的不同形式

如第5.2节和图5.4所述，本书将探讨概率预测的三种核心形式：分位数预测、密度预测和集合预测（不要与第10.3.2节中的随机森林等集成机器学习模型混淆）。这些可以归纳为两个核心类别：单变量（分位数和密度）和多变量（集合预测）。为理解这些类型，考虑尝试估计时间步 $t=N+1,N+2,\ldots,N+k$ 的数据分布的场景。

对于单变量预测，目标是估计每个时间步的需求分布。换言之，估计总共k个单变量分布。不同时间步的分布相互独立，这意味着一个时间步的变量分布不受其他时间步需求信息的影响。另一种说法是没有建模相互依赖关系。第3.1节介绍的单变量高斯分布是单变量分布函数最著名的例子。单变量分布的估计有两种主要形式：完整的连续密度函数或离散值（分位数）。

定义了等概率的均匀间隔水平。这两个形式在图11.1中针对标准高斯分布进行了说明。密度估计通常更可取，因为它描述了整个分布，但需要知道/假设分布来自特定参数族（例如此处的高斯分布），或者需要训练相对昂贵的方法，例如将在第11.5节中描述的内核密度估计。图11.1中的分位数估计（在第3.2节中更详细地描述）显示了高斯密度的5%、10%、...、95%分位数，显然对分布的描述性较差，但计算成本低得多，且不依赖对数据特定分布的假设。由于需要在预测范围内的每个k个时间步估计一个单变量分布，这种计算成本的降低特别有利。

对于多变量预测，任务是在预测范围内为所有k个需求变量估计一个单一的多变量分布（关于多变量分布的更多信息见第3.3节）。多变量分布的优势在于它们考虑了整个预测时间范围内的相互依赖关系。为了说明这一点，考虑家庭需求的例子。这主要由住户的行为决定。如果一个人上班迟到，那么他们很可能下班也晚，因此早上的需求变化将对应于晚上的变化。换句话说，由于这两个活动之间的联系，早上和晚上的需求之间存在相互依赖关系。因此，可以对多变量预测进行采样，以产生包含这些相关性的需求概况场景。从而可以模拟并利用更复杂和现实的相互依赖行为来优化诸如存储控制（第15.1节）等应用。

与单变量情况类似，可以估计完整的多元密度，但通常更复杂且难以精确建模。这些方法通常计算成本更高，用于拟合的包/资源更少，并且很少有标准的参数多元分布函数可用于拟合数据。相反，通常估计来自分布的有限样本。这些集成是来自分布的实现或“代表”。为了说明这一点，考虑一个基本情况$k=2$，其中两个连续时间步$t=1,2$的估计由二元高斯分布联合描述，如图11.2左侧所示。从该分布中抽取30个随机样本，得到右侧的时间相关二元集成。注意在第3.3节介绍的分布术语中，多变量概率预测是联合分布，而单变量概率预测是完整联合分布的边缘分布。

11.2 估计未来分布

如前所述，概率预测的目标是估计需求的未来分布，无论是单个时间步（单变量）还是多个时间步（多变量）。要估计不确定性，需要精确建模变异性。有几种标准的实用方法，将在本节中概述，这些方法是后续许多技术的基础。

第一种方法尝试通过对观测数据进行训练来直接建模分布。与大多数点预测一样，这些模型使用历史数据来捕捉变异性，并且通常假设过去的分布与未来的需求相关。参数模型（第11.3节）、核密度估计（第11.5节）和分位数回归（第11.4节）都以这种方式建模分布。

目标是直接训练分布模型的参数或超参数（例如高斯模型），或使用将估计分布（例如分位数）的模型。这些方法的优点是，只要使用正确的模型，并且它们在足够的目标分布数据上进行训练，就可以准确捕捉不确定性。例如，如果我们正在建模下午2点的需求，并且知道历史上的下午2点数据都来自同一分布，那么这些数据可用于估计真实分布。不幸的是，通常不能确定哪些数据来自同一分布，因此需要基于数据分析做出某些假设。这种方法的另一个缺点是其准确性与可用数据量相关。少量数据可能导致潜在不准确的估计。

第二种类型的模型不直接建模变异性，而是通过调整输入变量和/或模型参数将变异性引入模型。例如，假设有一个仅依赖于温度的需求模型。那么可以通过将不同的温度值插入模型来建模需求的变异性。通常这些是通过对温度的单个估计进行微调并模拟需求对温度的敏感性而形成的。

这种方法用于数值天气预报以产生预测集成/场景。对大气最可能的状态施加小的偏差，并将数值天气预报模型重新应用于调整后的状态，以产生一系列天气场景。对最终集成接近程度的分析可以指示对未来天气状态的置信度，而广泛分散的集成可能意味着未来天气高度不确定。

或者可以对模型参数进行小幅调整。这可以解释模型中的错误设定，并生成其他可能的未来状态。因此，多次调整可以产生一系列输出，从而估计未来的分布。这两种调整方法的难点在于，必须对输入/参数施加正确的偏差，才能产生准确的分布估计。在输入情况下，可以通过从历史观测中随机抽样或从可抽样的分布中进行估计来辅助实现这一点。

该模型的另一个缺点是需求变化并非直接模拟，而是估计模型对输入或参数的敏感性。考虑上面的温度示例。需求可能随温度变化，但实际上，人们对固定温度下的需求变化更感兴趣（假设温度可以准确预测）。关键在于对温度添加调整，以便捕捉这种变化。同样，交叉验证是一种可用于确定输入/参数的适当调整的方法。

以下各节将主要关注第一种生成概率预测的方法，并直接在历史观测上训练模型。

11.2.1 符号说明

在后续小节中，将介绍至少涵盖第11.1节中三种类型各一种的概率预测方法：分位数预测、密度预测和集合预测。对于后续章节，值得考虑以下符号和条件。

11.3 参数模型

1. 如前所述，设需求由时间序列 $L_{1},L_{2},...,$ 表示，其中 $L_t$ 是时间步 t 的需求。
2. 不失一般性，假设目标是预测时间戳 $t=N+1,N+2,\ldots,N+k$ 的 k 步超前需求。
3. 对于单变量概率预测：用 CDF $F_{1}(L_{N+1}|\mathbf{Z}),F_{2}(L_{N+2}|\mathbf{Z}),\dots,F_{k}(L_{N+k}|\mathbf{Z})$ 表示真实分布，每个预测时点对应一个函数，即 $F_{t}$ 是时间步 t 的需求的单变量分布。每个预测都依赖于先验信息 Z，Z 表示所有所需因变量的集合，如天气、历史需求等，这些变量决定了未来的需求。对于每个时间步 $t\in\left\{{1,2,\dots,k}\right\}$，相应的 CDF 预测记为 $\hat{F}_{t}(L_{N+t}|\mathbf{Z})$。为简洁起见，符号中可能不包含 Z。
4. 对于多变量概率预测，真实分布可以用单个 CDF $F_{t=1,\ldots,k}(\mathbf{L}|\mathbf{Z})$ 表示，它描述了多变量随机变量 $\mathbf{L}=(L_{N+1},L_{N+2},\ldots,L_{N+k})^{T}$ 的分布。先验信息 Z 包含所有因变量以及截至时间步 N 的历史负荷。为清晰起见，通常不包含 Z。
5. 集合预测的第 m 个集合通常记为 $\hat{\mathbf{L}}^{(m)}=$ $(\hat{L}_{N+1}^{(m)},\hat{L}_{N+2}^{(m)},\dots,\hat{L}_{N+k}^{(m)})^{T}$。

11.3 参数模型

参数分布模型是可取的，因为它们通常只需使用少量参数就能完整描述数据的分布。本节首先通过一个单变量分布的简单示例（第11.3.1节）讨论参数模型。个别的单变量参数模型通常过于僵化，无法精确拟合分布，但简单单变量分布的族可以“混合”起来估计更一般的形状，这将在第11.3.2节中介绍。

11.3.1 简单单变量分布

一些简单的单变量分布已在第3.1节中介绍。最常见的是高斯（或正态）分布，但也介绍了对数正态分布和伽马分布。除了这些函数外，还可以使用其他类似的函数对一系列分布进行建模。这类模型的优势在于，它们仅需训练少量参数即可完全估计分布。然而，受限于特定的函数形式，简单的参数模型无法估计更复杂的分布。例如，上述分布函数都是单峰的，这意味着它们描述的是具有单一众数（最大值）的分布。无法对多峰分布（具有多个不同最大值的分布）进行建模。单峰、双峰和三峰分布的示例如图11.3所示。

尽管单变量模型不太可能产生最精确的单变量概率预测，但它们可以作为基准模型，与本章后面描述的更复杂方法进行比较。此外，由于它们由相对较少的参数描述，因此可能比非参数模型更容易训练。训练参数模型需要估计描述所选分布族的每个单独参数。例如，高斯分布需要估计均值和标准差，而伽马分布则需要估计形状和尺度参数。对于高斯分布，均值和标准差可以通过最大似然估计（第8.2.1节）得到，这些值恰好分别是样本均值和样本标准差（第3.5节）。为了确保产生最佳可能的估计，需要仔细选择最合适的输入数据来训练参数（这与数据驱动的机器学习技术不同，后者会从所有数据中学习）。数据可以通过第6章概述的分析方法进行识别。例如，假设发现某些每小时数据具有强烈的日周期性，那么可能适合训练24个不同的模型，每个模型仅使用一天中特定小时的数据。

参数模型也存在于多变量模型中。特别地，存在高斯分布的多变量版本。如第11.1节所述，这些参数模型可用于生成多步预测的集合概率预测（通过估计一个 h 维参数多变量分布）。

11.3 参数模型

然而，能够准确捕捉各种概率预测行为的明确多变量分布要少得多。这使得它们与更通用的方法相比不那么适用，后者将在第11.6节和第11.7节中介绍，这些方法也能捕获预测范围内时间步长之间的相互依赖性。

11.3.2 混合模型

通过组合第11.3.1节中讨论的简单参数分布的混合，可以建模更通用的单变量分布。随机变量 $\mathbf{x}\in\mathbb{R}^{p}$ 的有限混合模型的概率密度函数的一般形式为

$$\begin{aligned} f({\textbf {x}} ) = \sum _{k=1}^K \pi _k g_k({\textbf {x}} , \theta _k), \end{aligned}\tag{11.131}$$

其中 $g_{k}(\mathbf{x})$ 是通常来自同一族（例如高斯分布）的概率密度函数，具有各自的参数 $\theta_{k}\left(\mathrm{e.g}\right.$（例如高斯分布的均值和标准差）。$\pi_{k}$ 是满足 $\sum_{k=1}^{K}\pi_{k}=1$ 的权重，通常称为混合概率。混合模型常用于聚类，在这种情况下，每个概率密度函数定义了来自某个簇的点的分布，权重表示每个簇中观测值所占的比例。

对于连续变量，最常见的混合模型使用高斯分量，即

$$\begin{aligned} g_k(\textbf{x}) = \frac{1}{(2\pi )^{N/2} det(\boldsymbol{\Sigma }_k)^{1/2}}exp \left( (\textbf{x}- \boldsymbol{\mu }_k )^T\boldsymbol{\Sigma }_k^{-1}(\textbf{x}- \boldsymbol{\mu }_k) \right) , \end{aligned}\tag{11.132}$$

具有协方差 $\Sigma_{k}\in\mathbb{R}^{p\times p}$ 和均值向量 $\boldsymbol{\mu_{k}}\in\mathbb{R}^{p}$。这被称为高斯混合模型（GMM）。图11.4展示了一个包含三个簇的GMM $(p=1)$ 的简单示例，混合概率分别为0.5、0.25和0.25，均值分别为1、3、6，所有标准差均为1。很容易看出，通过增加更多的组/簇，可以估计更复杂的分布。

尽管GMM相比单个高斯模型需要拟合更多参数，但可以通过一种称为期望最大化算法（EM）的迭代过程相对高效地求解，该过程为最大似然函数找到最优估计2（参见第8.2.1节）。

考虑观测值 $\mathbf{x_{1}},\mathbf{x_{2}},\ldots,\mathbf{x_{N}}\in\mathbb{R}^{p}$。不深入EM算法的细节，该过程在期望步骤（E步）和最大化步骤（M步）之间迭代：E步使用当前参数估计计算对数似然函数的期望，M步更新参数以最大化当前期望对数似然函数。对于GMM，这转化为以下步骤（每次迭代计算）：

1. 计算每个观测值 $\mathbf{X}_{i}$ 属于每个组 $k=1,\ldots,K$ 的后验概率 $\tau_{ik}$，

$$\begin{aligned} \tau _{ik} = \frac{\pi _k g_k(\textbf{x}_i, \theta _k)}{\sum _{k=1}^K \pi _k g_k(\textbf{x}_i, \theta _k)}. \end{aligned}\tag{11.133}$$

这是E步。

$\pi_{k}^{new}=\sum_{i=1}^{N}\tau_{k,i}$ $k=$ $1,\ldots,K$。

1. 更新每个分量 $k=1,\ldots,K$ 的均值，

$$\begin{aligned} \boldsymbol{\mu }_k = \frac{\sum _{i=1}^N \tau _{k, i}{} {\textbf {x}}_i}{\sum _{i=1}^N \tau _{k, i}}. \end{aligned}\tag{11.134}$$

注意，这是一个加权平均，权重基于隶属概率。

1. 更新每个分量 $k=1,\ldots,K$ 的协方差矩阵，

$$\begin{aligned} \boldsymbol{\Sigma }_k = \frac{\sum _{i=1}^N \tau _{k, i}({\textbf {x}}_i- \boldsymbol{\mu }_k)({\textbf {x}}_i- \boldsymbol{\mu }_k)^T}{\sum _{i=1}^N \tau _{k, i}}. \end{aligned}\tag{11.135}$$

即样本协方差的加权版本。

组的数量 K 是一个必须选择的超参数。虽然可以在交叉验证中选择，但似然函数框架意味着也可以使用信息准则（如第8.2.2节所述）来找到最合适的簇数。对于不同大小的簇，计算BIC（或AIC）。绘制BIC与簇数的关系图，可以找到增加簇数导致BIC下降收益递减的点。该点是图的“肘点”（因此这种启发式方法被称为“肘部方法”），并指示了一个合适的簇数选择。这里的“合适”是一个相对主观的术语，因为可能有其他几个原因使得不同的数量更合适或更有用。

该方法的一个示例如图11.5所示。这里，最佳聚类数约为四，因为切线在该值附近相交。切线通常用于更容易地识别肘部，从而确定聚类数。

拟合高斯混合模型是估计分布 $if$ 的一种简单方法，前提是所有数据来自同一分布。对于时间序列，这意味着用于训练的数据构成一个平稳序列。不幸的是，通常情况并非如此。一天中的不同时段可能具有不同的分布，时间序列可能依赖于天气、一年中的时间或许多其他变量。因此，尽管EM算法允许相对快速地训练高斯混合模型，但可能没有足够的数据来精确训练多个混合模型。

11.4 分位数回归与估计

用于概率负荷预测的大多数模型是非参数的，并且因其在建模分布方面提供了更大的灵活性而受欢迎。生成单变量概率预测的最简单、最常用的方法之一是分位数回归，这是本节的主题。该方法的一个优点是它是标准最小二乘回归的简单改编。

考虑估计时间步长 $t=N+1,N+2,\ldots,N+k$ 的 q 分位数（分位数介绍见第3.2节）。常用的选择是十分位数（10-分位数）或半十分位数（20-分位数），这样分布分别被分成10个或20个等概率区域。

考虑历史时间序列 $L_{1},L_{2},\dots,L_{N}$ 。回顾第9.3节，标准线性回归的目标是通过最小化与观测值的最小二乘差来找到某个预测模型 $f_{t}(\mathbf{Z},\beta)$ 的参数 $\beta$ 。用数学术语可以写成

$$\begin{aligned} \hat{\boldsymbol{\beta }} = \arg \min _{\boldsymbol{\beta }\in \mathcal {B}}\left( \sum _{t=1}^N(L_t - f_t(\textbf{Z}, \boldsymbol{\beta }))^2 \right) . \end{aligned}\tag{11.136}$$

这里 $\boldsymbol{B}$ 表示参数可以取值的可行集，通常是多维实空间 $\mathbb{R}^{p}$ ，其中 $p$ 是所选预测模型的参数个数。一旦找到参数，模型就可以用于使用新输入生成预测值。

分位数回归的原理相同，只是现在对于每个分位数 $\tau\in\{1,2,\ldots,q\}$ ，需要找到一组参数 $\hat{\boldsymbol{\beta }}_\tau$ τ，使得根据分位数损失函数，模型与观测值之间的差异最小化，即

$$\begin{aligned} \hat{\boldsymbol{\beta }}_\tau =\arg \min _{\boldsymbol{\beta }\in \mathcal {B}} \left( \min \sum _{t=1}^N c_{\tau }(L_t, f_t(\textbf{Z}, \boldsymbol{\beta })) \right) , \end{aligned}\tag{11.137}$$

其中代价函数 $c_{\tau}(x,y)$ 定义为

$$c_{\tau}(x,y)=\left\{\begin{array}{ll}\tau(x-y)&x\geq y\\(1-\tau)(y-x)&x<y\end{array}\right.,$$

对所有分位数重复这一过程。回想一下，这与第7章中介绍的钉板损失分数相同。分位数回归的过程比最小二乘回归稍微复杂一些，因为代价函数不可微。然而，该问题可以很容易地重新表述为线性规划问题，并且可以非常快速地求解。分位数回归仅适用于模型 $f_{t}(\mathbf{Z},\beta)$ ，该模型是参数的 β线性组合。这仍然允许在可以建模的关系类型中具有很大的灵活性。

为了说明该过程，考虑一个非常简单的例子。从均值为 $\mu=2$ 且标准差为 $\sigma=3$ μ 的高斯分布（见第3.1节）生成400个点，以表示一个包含400个点的时间序列。这里，y轴值是随机点，时间顺序就是它们被采样的顺序。时间序列在图11.6中以黑色显示。现在考虑一个形式为 $f_{t}(\beta)=at+b$ 的简单线性模型（即参数 $\beta=(a,b)^{T})$ ），由于该特定模型中没有依赖关系，因此没有其他输入Z。对该线性模型应用分位数回归，针对每个十分位数（或10-分位数）按照式(11.137)应用于时间序列。这些在图11.6中以红色虚线显示。注意，理论上，随着样本数量增加，最终的分位数应该是平坦的水平线，并且应该描述均值为 σ 的高斯分布的分位数

$\mu=2$和标准差$\sigma=3$。然而，在这种情况下存在轻微的梯度，因为只有相对少量的数据，抽样中的偏斜可能对模型拟合产生较大影响。

回顾第7章，概率积分变换（PIT）可用于评估概率预测的校准。分位数回归线应将数据分割成等概率的观测，这意味着每相邻十分位之间预期有40（400/10）个观测。图11.7中的PIT显示了这一点。注意，由于样本数量相对较少，某些分位数实际有39或41个观测。当选择了恰当的模型时，训练集上的PIT应呈现均匀分布。模型的真正评估一如既往地取决于在未见测试集（而非训练集）上的表现。此外，对于概率预测，校准和锐度都是重要属性，因此应使用第7.2节介绍的恰当评分函数来评估预测，而不仅仅是PIT。

最后，需要注意的是，由于每个分位数是独立训练的，有时分位数可能会相互交叉，这显然是不一致的，例如5%分位数可能高于10%分位数。为防止这种情况，可以在每个时间步重新排序每个分位数，以确保当$p<q$时，$p$百分位数低于$q$百分位数。

11.5 核密度估计方法

如第3.4节所示，核密度估计（KDE）可以看作是直方图的平滑版本。然而，不是在不同桶中添加离散计数，而是通过在观测位置添加核函数来估计连续分布。考虑随机变量X的观测$X_{j}$，j=1,...,N，则概率密度函数的KDE定义为

$$\begin{aligned} \hat{F}(X) = \frac{1}{Nh} \sum _{j=1}^{N} K \left( \frac{X-X_j}{h} \right) , \end{aligned}\tag{11.138}$$

其中h是带宽，估计的平滑参数，$K()$是某个核函数。核函数的一个流行例子是高斯核，定义为

$$\begin{aligned} K(x) = \frac{1}{\sqrt{2\pi }} \exp \left( -\frac{1}{2}x\right) . \end{aligned}\tag{11.139}$$

所选的核函数通常不如带宽的正确训练重要。另外注意，核函数的选择与真实分布无关，即高斯核并不意味着数据呈高斯分布。带宽的重要性如图11.8所示。该图比较了200个观测的直方图（左）与相同观测的KDE（右，三种不同带宽）。选择带宽过小，KDE会过拟合观测；带宽过大，KDE会欠拟合并具有更高偏差（回忆第8.1.2节偏差-方差权衡）。虽然选择带宽有一些经验法则，但这些通常基于对底层分布的强假设（如高斯分布），因此对于负荷预测而言过于严格。相反，可以通过交叉验证选择带宽，最小化验证集上估计与观测之间的拟合（通常根据最小化概率评分函数如CRPS来定义，见第7章）。这可以通过超参数空间搜索（例如网格搜索，第8.2.3节）实现。尽管最简单形式的KDE只有一个参数，但计算过程可能很昂贵。

KDE可以很容易地适用于概率时间序列预测。再次考虑历史观测$L_{t}$，对于$t=1,\ldots$，N，来自某个随机变量，假设它们来自同一分布，目标是为预测范围内的每个时间步$t=N+1,\ldots,N+k$生成k步超前密度预测。在这种简化情况下，最基本的核密度估计可定义为

$$\begin{aligned} \hat{F}_i(L_{N+i}) = \frac{1}{Nh} \sum _{t=1}^{N} K \left( \frac{L-L_t}{h} \right) , \end{aligned}\tag{11.140}$$

对于预测范围$i=1,\ldots,k$中的每个时间步，以及带宽h。换句话说，假设每个时间步的分布相同。这显然是不现实的，原因有几个。首先，较旧的数据可能不如较新的信息相关。此外，方程(11.140)中的KDE估计也独立于任何其他输入，例如温度或一天/一周中的时间。为了纠正这些缺点，可以使用简单KDE估计的修改版本。

为了降低较旧点的影响，可以引入一个简单的衰减因子 $0<\lambda\leq 1$，其中 λ。这可以降低较旧数据对整体分布函数的贡献。一种可能的实现方式是

$$\begin{aligned} \hat{F}_i(L_{N+i})=\sum _{t=1}^{N} w_t K \left( \frac{L-L_t}{h} \right) , \end{aligned}\tag{11.141}$$

其中指数衰减权重为 $w_{t}=\frac{\lambda^{N-t}}{h\sum_{l=1}^{N}\lambda^{N-l}}$ N −t λ。在这种情况下，衰减因子 λ 和带宽都需要进行优化（同样，通常通过交叉验证）。当然，也可以使用其他权重，如下文的条件KDE形式所示。如果更多历史数据可能相关，那么较慢的衰减（例如线性衰减）可能更值得关注。唯一的限制是权重之和应为1，以确保最终函数仍然是一个定义良好的概率分布。

另一个简单的修改是仅训练特定的历史点。例如，负荷数据通常具有整数周期 s 的简单季节性（如每日或每周）。在这种情况下，可以使用以下公式估计密度

$$\begin{aligned} \hat{F}_i(L_{N+i}) = \frac{1}{N_s h} \sum _{t \in \mathcal {I}_i} K \left( \frac{L-L_t}{h} \right) , \end{aligned}\tag{11.142}$$

其中 $N_{s}$ 是集合 $\mathcal{T}_{i}=\{j|N+i-j=sk$ 的元素个数，对于某个整数 $k\}.\\mathcal{T}_{i}$ 是周期 s 的常数倍索引。因此，对于具有每日季节性的小时数据，要预测 $i=2$ 时间段（例如，如果预测原点在午夜，则为凌晨2点），用于构建KDE估计的历史数据将仅使用每个历史日期的凌晨2点数据。回顾一下，周期性和季节性可以通过第5章介绍的方法找到。

另一个对简单KDE的流行更新是条件核密度估计。它估计变量 $L_{i}$ 在依赖于某些自变量（例如 T, S）的条件下的分布。现在，我们可以利用独立-依赖观测对作为 $(T_{t},S_{t},L_{t})$，并定义 $\hat{F}_{i}(L_{i}|T,S)$ 的条件分布为

$$\begin{aligned} \hat{F}_i(L_{N+i}| T, S)=\sum _{t=1}^{N} \frac{ K((T_t-T)/h_T) K((S_t-S)/h_S)}{\sum _{l=1}^{n} K((T_l-T)/h_T) K((S_l-S)/h_S)} K \left( \frac{L-L_t}{h} \right) , \end{aligned}\tag{11.143}$$

其中 $h_{T},h_{S}$ 分别是表示 T 和 S 的分布的带宽，现在必须与依赖序列的带宽 h 一起确定。与之前一样，这仅仅是类似于式(11.141)的加权和，但权重是基于自变量的核函数。通常，自变量的常用选择是天气变量，以及一周中的时间段。这也可以扩展或简化，以分别考虑更少或更多的变量。然而，每个带宽都会显著增加模型训练的计算成本，在超过两个条件独立变量的情况下可能不切实际。

为了加速优化，通常将变量归一化（例如到[0,1]）以缩小搜索空间（见6.1.3节）。训练完成后可以对预测进行重新缩放。如3.4节所述，有不同核函数的选择，可以测试不同的核函数，尽管通常选择对预测精度的影响很小[1]。

这里介绍的不同修改显然可以组合起来创建其他模型。例如，式(11.143)中所示的条件核密度形式可以扩展，以包含式(11.141)中的衰减因子，或者像式(11.142)那样对输入施加限制。与许多KDE方法一样，其缺点是每次修改通常会增加训练复杂度和计算成本。

11.6 集成方法

本节介绍集合预测，即来自同一预测原点的一组点预测，以相同的预测水平（长度为 h 个时间步）估计每个时间步。这些点预测是从一个 h 维多元分布中抽取的等概率样本，该分布表示预测水平上的联合分布（有关联合分布的更多信息，请参阅第3.3节）。换句话说，每个集成代表一条等可能的负荷轨迹。本节描述的方法无需生成完整的联合分布即可产生这些集成。

11.6.1 残差自助法集成（同方差）

本节描述了一种生成集合预测的方法，在假设时间序列具有固定方差的情况下，这些集合预测可视为完整多元概率预测的实现。首先考虑一个单步向前点预测模型，例如，这可以是第9.2节中的指数平滑模型或第9.4节中的ARIMA模型。将其记为 $f(L_{t}|\mathbf{Z},\beta)$，它可能使用历史数据以及任何解释性输入（全部由变量集 $\mathbf{Z})$ 描述），以生成真实值 $L_{t+1},\text{ at }\t+1$ 的估计值 $\hat{L}_{t+1|t}$。$\beta$ 是模型的参数。下一个时间步的预测可以通过迭代应用模型，并将前一个时间步的预测作为伪历史输入包含在模型中来获得。因此，对于下一个时间步

$$\begin{aligned} \hat{L}_{t+2|t} = f(\hat{L}_{t+1|t}| L_t, \textbf{Z}, \beta ) = f(f(L_t| \textbf{Z}, \beta )| L_t, \textbf{Z}, \beta ). \end{aligned}\tag{11.144}$$

该过程显然可以重复进行以生成k步超前预测。现在回想一下，存在一个误差过程描述了观测值与1步超前预测之间的差异，该差异由残差$\boldsymbol{\epsilon}_{t+1}=\boldsymbol{L}_{t+1}-\boldsymbol{\hat{L}}_{t+1|t}$描述。通过将残差序列所描述的小偏差纳入预测，可以创建不同的轨迹，这些轨迹代表了不同但同样可能的结果。

为了更详细地描述该算法，考虑从一个单步超前模型（例如$\hat{L}_{t+1}=f(L_{t}|\mathbf{Z},\beta)$）生成的k步超前预测。假设预测原点位于时间步$t=N$。因此，目标是生成覆盖时间段$N+1,\ldots,N+k$的集合预测。残差序列$\boldsymbol{\epsilon}_{t}=\boldsymbol{L}_{t+1}-\boldsymbol{\hat{L}}_{t+1}$（在整个训练集上计算）被假定具有固定方差（该序列称为同方差）且彼此不相关。利用该残差序列，使用单步超前预测模型$f_{:}$生成k步超前预测的新集合的过程相对简单。对于每个集合b，过程如下：

1. 随机有放回地抽取一个残差（这称为自助样本），$\hat{e}_{1}^{(b)}$ $\{\epsilon_{1},\epsilon_{2},\ldots,\epsilon_{N}\}$

2. 将该残差加到当前的1步超前预测值${\tilde{L}}_{N+1|N}$上，以产生一个新值$\hat{L}_{N+1|N}^{(b)}=\tilde{L}_{N+1|N}+\hat{e}_{1}^{(b)}$。

$\hat{L}_{N+1|N}^{(b)}$ $\tilde{L}_{N+2|N+1}^{(b)}=f(\hat{L}_{N+1|N}^{(b)}|L_{N},\mathbf{Z},\beta)$

1. 使用来自残差序列的另一个自助样本来更新该值，得到$\hat{L}_{N+2|N+1}^{\dot{(}b)}=\tilde{L}_{N+2|N+1}^{(b)}+\hat{e}_{2}^{(b)}$
2. 继续该过程直到达到第k步。

$\hat{L}_{N+1|N}^{(b)},\hat{L}_{N+2|N+1}^{(b)},\ldots\hat{L}_{N+k|N+k-1}^{(b)}$ ble.

这也称为残差自助预测。该过程可以重复以生成任意数量的集合。生成的集合越多，捕获的负荷多样性就越多。生成更多样本会增加计算成本，但由于每个集合彼此独立，因此可以并行生成。请注意，该方法强烈假设未来的1步超前误差将与过去的1步超前误差相似。图11.9a展示了具有同方差性的简单周期序列的示例。

如果不是从实际残差中抽样，而是从假设或拟合的分布中抽样，那么该方法可称为蒙特卡洛预测。例如，通常假设残差服从均值为零的高斯分布，因此不必从残差集合中抽样，而是可以从基于残差训练的高斯分布中抽样。图11.10展示了一个由蒙特卡洛模拟生成的集合预测示例，针对一个简单的ARIMA(4,1,1)模型，包含100个集合。注意，误差随着预测视界长度的增加而变宽（具有更大的变异性）。这是由于误差从一步到下一步的累积。这直观上是合理的，因为预测越远，不确定性应该越大。

值得注意的是，每个时间步的预测可用于估计单变量估计。这可以通过在每个时间步对集成点集合拟合分位数或密度估计来实现。

11.6.2 残差自助法集成（异方差性）

第11.6.1节所述自助法的一个优点是，仅需训练原始点预测模型，便能以极低的计算成本生成多元预测。此外，如果模型包含自回归特征（如ARIMA和指数平滑），则集成还能保留时间序列的相互依赖性。该方法的缺点是对残差序列的强同方差假设。实际上，高需求时期可能具有更大的变异性。方差随时间变化的时间序列称为具有异方差性。

图11.9b展示了一个具有异方差的简单周期序列示例，其中需求的最大变化与周期的最大幅度一致。当时间序列存在异方差时，可以使用所谓的GARCH类模型来纳入变异性，这些模型可扩展第11.6.1节所述的自助法。本节概述了这些方法，但由于它们相对复杂，故省略细节。感兴趣的读者可参考附录D中的进一步阅读材料。

考虑一个负荷时间序列模型，包含某个均值过程（例如第9章和第10章中给出的任何点预测），并像之前一样令$\epsilon_{t}$为1步超前误差/残差项。现在假设标准差$\sigma_{t}$也随时间t变化。将残差写成以下形式是有帮助的：

$$\begin{aligned} \epsilon _{t} = \sigma _{t} Z_t, \end{aligned}\tag{11.145}$$

其中$\sigma_{t}>0$是时变条件标准差，而$(Z_{t})_{t_{\in}\mathbb{Z}}$是一个σ随机变量，它是平稳的、时间独立的，并满足条件$\mathbb{E}(Z_{t})=0$且$\mathbb{V}ar(Z_{t})=1$。将数据按标准差拆分有助于将实际变化简化为由标准差表示的变化幅度以及由于缩放而变得平稳的随机分量。一旦找到式(11.145)的分量，就可以通过对$Z_{t}$的分布进行随机抽样，然后用预测时间步的建模标准差$\sigma_{N+k}$重新缩放，从而轻松应用改编的自助法。

为此，首先必须为标准差选择一个模型。在经济模型中，标准选择是ARCH和GARCH模型。ARCH和GARCH模型本质上是第9.4节中引入的点预测AR和ARMA模型的方差对应物。GARCH(p, q)模型形式如下：

$$\begin{aligned} \sigma ^2_t = \alpha _0 +\sum _{i=1}^q \alpha _i\epsilon ^2_{t-i}+\sum _{j=1}^p \beta _j \sigma ^2_{t-j}, \end{aligned}\tag{11.146}$$

其中$q$是滞后残差项（称为ARCH项），$p$是滞后方差项（称为GARCH项）。ARCH模型只有残差项的$q$阶滞后，没有方差项。这些模型通常与相应的AR或ARIMA模型（用于均值过程）结合，并使用普通最小二乘法或最大似然法求解。

ARCH和GARCH是方差的具体形式，通常适用于金融时间序列应用。实际上，标准差可以更一般地建模，这些方法将被称为GARCH类方法。在负荷预测中，需求的变化通常在需求较高的时间段更大（当然，对于每个新时间序列，在选择模型之前应分析方差模式）。因此，通常适合选择与所选点预测模型相似的标准差模型。

对于ARCH和GARCH问题，程序相对成熟，因此有许多自动选择模型阶数和系数的软件包。在更一般的情况下，可遵循改编自文献[2]的简单程序：

1. 由于式(11.145)给出的残差分量未知，因此考虑绝对残差或平方残差$\left|\epsilon_{t}\right|$或$\epsilon_{t}^{2}$。前者将是重点，因为它能更好地捕捉能源需求的厚尾性。然而，相同的程序适用于两种形式。
2. 由于$\mathbb{E}(|\epsilon_{t}|)=\sigma_{t}\mathbb{E}(|Z_{t}|))$，对$\lvert\epsilon_{t}\rvert$拟合模型等价于对标准差的缩放版本$C\sigma_{t}$拟合模型，其中$C>0$为某个常数，因为$Z_{t}$是平稳变量（因此具有恒定期望）。通常使用普通最小二乘法（第8.2节）进行拟合。
3. 现在必须通过考虑标准化残差$\epsilon_{t}/|\epsilon_{t}|=\epsilon_{t}/C\sigma_{t}=Z_{t}/C$（即$Z_{t}$的缩放版本）来估计常数C。由于$Z_{t}$的方差等于1，因此可以通过计算标准化残差的样本标准差来估计缩放因子C，这告诉我们$C=1/\sqrt{\alpha}$。

对于这种方法，注意没有对潜在分布（仅对其方差）做出假设。第11.6.1节中自助法预测的更新版本通过从Z的分布中抽取样本（从标准化残差$\epsilon_{t}/\sigma_{t})$形成的经验分布$Z$中抽取，然后由当前时刻的$\sigma_{t}$模型缩放）来更新一步提前预测。注意假定误差$\epsilon_{t}$和标准化残差在时间上不相关，以便它们可以在自助法中随机抽样。通常应检查残差（这里指标准化残差），以确保它们满足关于相关性、固定方差和零均值的假设。如果这些假设不成立，则可以应用进一步的更新（如第7.5节所述）。上述GARCH型模型的一个具体示例将在第14.2节的LV案例研究中给出。

11.7 用于多变量预测的Copula模型

本节介绍Copula，这是另一种生成多变量概率预测的流行方法，广泛应用于量化金融，但现在也广泛用于能源预测。这里仅描述基础知识。建议的进一步阅读材料见附录D。

Copula本质上是一个函数C，从N维单位立方体映射到一维单位立方体，即$C:[0,1]^{N}\longrightarrow[0,1]$，它描述了变量$U_{1},\dots,U_{N}$上的累积分布函数，其中每个变量$U_{i}$在[0,1]上均匀分布，即$C(u_{1},u_{2},\ldots,u_{N})=P(U_{1}\leq u_{1},U_{2}\leq u_{2},\ldots,U_{N}\leq u_{N})$。换句话说，Copula对具有均匀边际分布的变量$U_{1},\dots,U_{N}$的联合分布进行建模（关于联合分布和边际分布的更多细节见第3.3节）。Copula关注变量之间的相关性/相互依赖结构。这些方法的一个优点是存在多种不同的Copula函数可用于建模相互依赖关系。

Copula的强大之处在于，根据Sklar定理，任何多变量分布都可以仅使用边际分布和一个Copula来建模。考虑一个具有联合分布的随机变量$\mathbf{X}=(X_{1},X_{2},\ldots,X_{N})^{T}\in\mathbb{R}^{N}$

$F_{X_{1},\ldots,X_{N}}(x_{1},\ldots,x_{N})$及其边际CDF记为$F_{1}(x_{1}),\ldots,F_{N}(x_{N})$，那么由Sklar定理，存在一个Copula C使得

$$\begin{aligned} F_{X_1, \ldots , X_N}(x_1, \ldots , x_N) = C(F_1(x_1), \ldots , F_n(x_N)). \end{aligned}\tag{11.147}$$

注意，对于任何具有CDF F的随机变量X，当其通过CDF变换时是均匀分布的，即$F_{X}(X)$均匀分布。回忆这就是第7章描述的概率积分变换（PIT）。

式(11.147)的过程是可逆的，这意味着一旦在观测数据上训练了Copula模型，就可以很容易地生成具有相关依赖结构的多变量样本。首先从Copula分布生成样本$u_{1},u_{2},\dotsc,u_{N}$，然后使用样本点每个对应分量的边际$F_{i}^{-1}(u_{i})$的逆CDF将每个变量变换回原始随机变量空间。

假设$\mathbf{X}=(X_{1},X_{2},\ldots,X_{N})^{T}\in\mathbb{R}^{N}$具有多变量高斯分布，则相应的Copula完全由相关矩阵定义，因为这解释了整个依赖结构。这也意味着每个相关矩阵定义了一个特定的高斯Copula。当然，仅仅因为一个多变量随机变量具有高斯Copula并不意味着它遵循多变量高斯分布，因为边际分布不一定是高斯的。如果相关矩阵是单位矩阵，则称该Copula为独立Copula，定义为

$$\begin{aligned} C_0(u_1, u_2, \ldots , u_N) = u_1\ldots u_N, \end{aligned}\tag{11.148}$$

其中每个分量独立于其他分量。一般情况下，高斯Copula没有简单的解析形式，但可以表示为

$$\begin{aligned} C^{Gauss}_R(u_1, u_2, \ldots , u_N) = \Psi _R(\Psi ^{-1}(u_1), \ldots , \Psi ^{-1}(u_N)), \end{aligned}\tag{11.149}$$

其中Φ是标准高斯（均值为零、标准差为1）的单变量CDF，$\Psi_{R}$是均值为零、相关矩阵为R的多变量高斯分布。

另一类流行的copula是阿基米德copula，它们具有显式公式，并且仅使用一个参数即可表示多元分布。它们的一般形式为

$$\begin{aligned} C(u_{1},\dots ,u_{N};\theta )=\psi ^{-1}\left( \psi (u_{1};\theta )+\cdots +\psi (u_{N};\theta );\theta \right) \end{aligned}\tag{11.150}$$

其中$\psi:[0,1]\times\Theta[0,\infty)$是一个连续、严格递减的凸函数ψ，使得$\psi(1;\theta)=0$。例如，一种流行的阿基米德copula是Gumbel ψθ Copula，其定义为

$$\begin{aligned} C_{Gum}(u_1, u_2, \ldots , u_N| \theta ) = \exp \left[ -\left( (-\log (u_1))^{\theta }+\cdots +(-\log (u_N))^{\theta }\right) ^{1/\theta }\right] . \end{aligned}\tag{11.151}$$

注意，当θ=1时，这退化为独立copula。来自双变量Gumbel copula的不同θ值的样本示例如图11.11所示。注意，Gumbel copula永远无法表示负相关性。

显然，不同的copula适用于不同的依赖结构，并非所有copula都适用于所有类型的数据。例如，Gumbel copula不应与具有负相关性的数据一起使用。如何选择和拟合copula将在本节后面简要讨论。

Pearson相关系数的值依赖于边缘分布和copula，这意味着随机变量在使用边缘CDF变换后会有不同的Pearson值。对于copula的相关结构，更便捷的度量是秩相关系数，它仅依赖于数据的秩（参见第3.5节）。由于数据的秩在应用单调递增函数时保持不变，这意味着应用边缘CDF后秩相关系数不会改变。第3.5节给出了一个秩相关系数的例子：Spearman秩相关系数。

为了更好地理解copula的工作原理以及如何利用它们生成新样本，考虑一个简单的双变量分布，随机变量(X1, X2)的依赖结构由具有协方差的高斯copula描述

$$R=\left[\begin{array}{cc}1&0.8\\0.8&1\end{array}\right]$$

即它具有线性Pearson相关系数$\rho=0.8$。还假设第一个变量X1的边缘分布由Gamma分布描述，其中

$$\begin{aligned} Gamma(X, \alpha , \beta ) = \frac{1}{\beta ^\alpha \Gamma (\alpha )}\int _{0}^{X} t^{\alpha -1} e^{-t/\beta } dt \end{aligned}\tag{11.152}$$

其中$\alpha=2$是形状参数（决定分布的形状），$\beta=1$是尺度参数（决定分布的离散程度）。$\Gamma(.)$是所谓的伽马函数。第二个变量X2由标准高斯分布描述（均值为零，标准差为1）。该分布的1000个观测值示例如图11.12所示。图的水平和垂直方向分别是每个变量的边缘直方图，分别显示Gamma分布和高斯分布。

应用每个边缘的CDF到各自变量后的分布如图11.13所示。边缘分布现在如预期般由均匀分布描述。注意在本例中，X1与X2之间的样本Pearson相关系数为$\rho=0.767$，但$U1$与$U2$之间的为$\rho=0.785$，并且该变换不保持ρ ρ（尽管在本例中它们接近）。相比之下，X1与X2之间以及U1与$U2$之间的Spearman相关系数均为$r=0.787$，正如预期。在实践中，目标是对这些变换后的数据拟合一个copula。在此情况下，假设已知copula是高斯copula，因此目标是找到Pearson相关系数$\rho$。实际上，在[3]双变量情况下，Spearman系数r与线性相关的关系为

$$\begin{aligned} \rho = 2\sin \left( r \frac{\pi }{6} \right) , \end{aligned}\tag{11.153}$$

因此，斯皮尔曼相关的不变值可用于估计皮尔逊相关，并得到 $\rho=0.801$（注意，由于样本中的数值偏差，它并不恰好是用于生成数据的 ρ 值 0.8，样本越大，样本值越接近原始参数）。

给定连接函数，可以生成新样本，然后使用边缘分布的逆累积分布函数将其转换回原始空间（保持相同的线性相关）。图11.14展示了使用连接函数生成并通过逆累积分布函数转换后的1000个新点示例。注意最终分布如何成功地与图11.12中的原始分布相似。

类似的过程可用于生成集合需求预测。在此情况下，考虑一个需求时间序列 $L_1, L_2, \ldots ,$，目标是生成一个日前多变量预测，例如预测原点 $t=N$。为简化起见，假设数据是每小时采样的，因此考虑的是24步向前预测。这里的目标是生成该日的多变量分布 $F_{N+1,\ldots,N+24}(L_{N+1},\ldots,L_{N+24})$，从而建模一天中不同时刻之间的相互依赖关系。假设一天中不同时刻的边缘分布的 CDF 已知，即 $F_{1}(L_{N+1}),\dots,F_{24}(L_{N+24})$ 已知。这些可以通过例如本章前面描述的单变量概率模型来估计。在这种情况下，可以使用连接函数，通过对经过边缘分布转换的日负荷曲线进行训练，来建模日内依赖结构。

存在一系列不同的连接函数模型，上文仅提及其中一部分。如何根据数据训练和选择连接函数仍然是问题。

理论上，对于连接函数和边缘分布的参数模型，可以使用最大似然估计（见章节8.2）来拟合连接函数，但对于高维问题，由于需要训练大量参数，这可能很复杂。相反，模型可以使用一种称为伪最大似然的两步过程来估计，其中首先估计边缘分布，然后最大化最大似然的简化形式，如下所示：

$$\begin{aligned} \sum _{k=1}^n \log \left[ c \{ \hat{F}_1(X_{1, k}), \ldots , \hat{F}_N(X_{N, k}) | \boldsymbol{\theta }_C \}) \right] , \end{aligned}\tag{11.154}$$

被最大化。其中 c 是与连接函数 CDF C 对应的连接函数密度，$\boldsymbol{\theta}_{C}$ 表示连接函数模型的参数。上述最大化仍然可能很复杂，特别是对于高维问题，具体取决于所考虑的连接函数模型。对于上面给出的简单高斯连接函数示例，已经提出了一种方法。可以使用每对变量的斯皮尔曼相关系数来估计相关矩阵，这既可以用作最终的相关矩阵，也可以用作式(11.154)中伪最大似然优化的初始猜测。

选择正确的连接函数取决于许多因素，详细研究超出了本书的范围。建议进一步阅读附录 D。总之，选择取决于所建模的相关类型以及分布尾部/极值内的依赖性。一种选择合适连接函数模型的可能方法是基于在验证集上的比较，如章节 8.1.3 所述。

11.8 习题

对于需要真实需求数据的问题，请尝试使用附录D.4中列出的一些数据。最好选择至少一年小时或半小时分辨率的数据。在所有使用这些数据的情况下，按照60%、20%、20%的比例划分训练集、验证集和测试集（见第8.1.3节）。

1. 从一个五维高斯分布中抽取20个样本。通过操作协方差矩阵中变量之间的相关性，确保某些变量比其他变量更相关。你可以将所有变量的方差固定为一，以简化模型。现在将高斯的每个维度视为一个五点时间序列中的不同时间步。绘制每个样本以创建类似图11.2的样本图。你在高度相关的变量之间看到了什么？如果改变不同变量的方差，集合图会如何变化？
2. 生成一个分位数回归。采用你在章节9.7中生成的线性预测模型。现在，使用内置包（如 R 中的 quantreg 包）对训练数据拟合分位数为 10、20、30、...、90 的分位数回归。应用于测试集，并统计落在每组分位数之间的值的数量。绘制概率积分变换图。它是什么形状？模型是否存在偏差？是欠分散还是过分散？对分位数进行何种调整有助于产生均匀的 PIT？
3. 对于具有日或周周期性的需求数据，为每个时间步（来自季节周期中的每个周期）拟合一个核密度估计。例如，如果数据是半小时频率且具有日季节性，则为一天中的每半小时训练48个模型。通过带宽的网格搜索来拟合模型。使用最终模型应用于测试集。生成估计的分位数，从而计算与前一个问题相同百分位点的PIT。PIT是均匀的、过度分散的还是欠分散的？

参考文献

1. J. Jeon, J.W. Taylor，使用条件核密度估计进行风电功率密度预测. J. Amer. Stat. Ass. 107(497), 66–79 (2012)
2. S. Haben, G. Giasemidis, F. Ziel, S. Arora, 短期负荷预测及低压等级温度的影响。Int. J. Forecast. 35, 1469–1484 (2019)
3. T.C. Headrick，关于皮尔逊积矩相关与斯皮尔曼等级相关系数之间关系的注记. Open J. Stat. 6, 1025–1027 (2016)

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