Python-现代时间序列预测第二版-二-

Python 现代时间序列预测第二版（二）

原文：annas-archive.org/md5/22eab741fce9c15dfad894ecf37bdd51

译者：飞龙

协议：CC BY-NC-SA 4.0

第二部分

时间序列的机器学习

在这一部分，我们将探讨如何应用现代机器学习技术进行时间序列预测。本部分还涵盖了强大的预测组合方法以及令人兴奋的全球模型新范式。通过本部分的学习，你将能够使用现代机器学习技术为时间序列预测建立模型流水线。

本部分包括以下章节：

• 第五章，将时间序列预测视为回归问题

• 第六章，时间序列预测的特征工程

• 第七章，时间序列预测的目标转换

• 第八章，使用机器学习模型进行时间序列预测

• 第九章，集成与堆叠

• 第十章，全球预测模型

第五章：将时间序列预测作为回归问题

在本书的前一部分中，我们对时间序列有了基本的理解，并装备了分析和可视化时间序列的工具和技术，甚至生成了我们第一个基线预测。到目前为止，我们主要介绍了经典和统计技术。现在，让我们深入了解现代机器学习，并学习如何利用这个相对较新的领域来进行时间序列预测。机器学习是近年来迅速发展的一个领域，能够利用这些新技术进行时间序列预测，将在今天的世界中成为一项无价的技能。

在本章中，我们将讨论以下主要主题：

• 了解机器学习的基础

• 将时间序列预测作为回归问题

• 局部模型与全局模型

了解机器学习的基础

我们希望使用机器学习进行时间序列预测。但在开始之前，让我们花一些时间来了解什么是机器学习，并建立一个框架来展示它的功能（如果你已经对机器学习非常熟悉，可以跳到下一节“将时间序列预测作为回归问题”，或者继续跟我们一起复习这些概念）。1959 年，Arthur Samuel 将机器学习定义为“一种使计算机能够在没有明确编程的情况下学习的研究领域。”传统上，编程是一种我们知道一套规则/逻辑来执行某个动作，并且在给定的数据上执行该动作以获得我们想要的输出的范式。而机器学习则颠覆了这一点。

在机器学习中，我们从数据和输出开始，要求计算机告诉我们通过哪些规则可以从数据中获得期望的输出：

图 5.1：传统编程与机器学习

机器学习中有许多种问题设置，如监督式学习、无监督式学习、自监督学习等，但我们将专注于监督式学习，这是最常见的，也是本书内容最适用的。监督式学习指的是我们在程序、数据和输出的范式转换示例中已经提到的内容。我们使用一个包含输入和预期输出的配对样本的数据集，并要求模型学习它们之间的关系。

让我们从一个小的讨论开始，逐步构建整个示意图，它包含了监督式机器学习问题的主要关键组件：

图 5.2：监督式机器学习示意图，部分 1—理想函数

正如我们已经讨论过的，我们希望机器学习从数据中学习并得出一套规则/逻辑。在数学中，与逻辑/规则最接近的类比是函数，它接受一个输入（这里是数据）并提供一个输出。从数学上看，可以写作如下：

y = g(X)

其中，X 是特征集合，g 是理想目标函数（在图 5.2中用1表示），它将 X 输入（在示意图中用2表示）映射到目标（理想）输出，y（在示意图中用3表示）。理想目标函数在很大程度上是一个未知函数，类似于我们在第一章《引入时间序列》中看到的数据生成过程（DGP），它不在我们的控制之下。

图 5.3：监督式机器学习示意图，第二部分—学到的近似值

但我们希望计算机能够学习这个理想目标函数。这个理想目标函数的近似值用另一个函数 h 表示（在示意图中用4表示），它接受相同的特征集 X，并输出预测的目标，（在示意图中用5表示）。 是 h 函数的参数（或模型参数）：

图 5.4：监督式机器学习示意图，第三部分—将所有内容整合

现在，我们如何找到这个近似 h 函数及其参数，？通过示例数据集（在示意图中用6表示）。监督式机器学习问题的前提是我们能够收集一组包含特征 X 和相应目标 y 的示例，这些目标在文献中也称为标签。计算机就是从这组示例（数据集）中学习近似函数 h 以及最优模型参数，。在之前的示意图中，唯一真正未知的实体是理想目标函数 g。因此，我们可以使用训练数据集 D 来为数据集中的每个样本预测目标。我们已经知道所有示例的理想目标。我们需要一种方法来比较理想目标和预测目标，这就是损失函数（在示意图中用7表示）发挥作用的地方。损失函数告诉我们，使用近似函数 h 时，我们离真实结果有多远。

尽管 h 可以是任何函数，但它通常是从一个著名函数类别 H 中选择的。H 是一个有限的函数集，可以拟合数据。这个函数类别就是我们口语中所说的模型。例如，h 可以从所有线性函数或所有基于树的函数中选择，等等。从 H 中选择一个 h 是通过超参数（由模型设计者指定）和模型参数（从数据中学习）来组合完成的。

现在，剩下的就是运行不同的函数，找到最好的近似函数 h，它给我们带来最低的损失。这是一个优化过程，我们称之为训练。

让我们还看看一些关键概念，这些概念将在接下来的讨论中非常重要。

有监督的机器学习任务

机器学习可以用来解决各种任务，如回归、分类和推荐。但由于分类和回归是最常见的两类问题，我们将花一点时间回顾它们是什么。

分类任务和回归任务的区别非常简单。在机器学习框架图（图 5.2）中，我们讨论了 y，即目标。这个目标可以是一个实数值，或者是一个项的类别。例如，我们可能在预测下周的股价，或者我们只预测股价是上涨还是下跌。在第一种情况下，我们是在预测一个实数值，这叫做回归。在另一种情况下，我们预测的是两个类别中的一个（上涨 或 下跌），这叫做分类。

过拟合和欠拟合

机器学习系统面临的最大挑战是，我们训练出的模型必须在新的、未见过的数据集上表现良好。一个机器学习模型在这方面的能力被称为泛化能力。机器学习中的训练过程类似于数学优化，但有一个微妙的区别。数学优化的目标是到达提供的数据集中的全局最大值。而在机器学习中，目标是通过使用训练误差作为代理，达到较低的测试误差。一个机器学习模型在训练误差和测试误差上的表现与过拟合和欠拟合的概念密切相关。让我们通过一个例子来理解这些术语。

机器学习模型的学习过程与人类学习有很多相似之处。假设三名学生，A、B 和 C，正在为考试做准备。A 是个懒汉，前一天晚上去夜店了。B 决定死记硬背，把教科书从头到尾都背一遍。C 上课时认真听讲，理解了考试的内容。

正如预期的那样，A 没有通过考试，C 得了最高分，B 还行。

A因为没有学到足够的内容而挂科。这在机器学习模型中也会发生，当它们没有学到足够的模式时，这被称为欠拟合。其特点是训练误差和测试误差都很高。

B的成绩没有预期的高；毕竟，他们确实把整篇文章逐字逐句地记住了。但考试中的许多问题并不是直接来自教科书，B没有能够正确回答这些问题。换句话说，考试中的问题是新的和未见过的。而且，由于B记住了所有内容，但没有努力理解基础概念，B没有能够将所学知识推广到新的问题上。在机器学习中，这种情况被称为过拟合。其特点通常是训练误差和测试误差之间的差距很大。通常我们会看到非常低的训练误差和很高的测试误差。

第三个学生，C，学得正确并理解了基础概念，因此能够推广到新的和未见过的问题。这也是机器学习模型的理想状态。这种状态的特点是合理较低的测试误差，以及训练误差和测试误差之间的小差距。

我们刚刚看到了机器学习中的两个最大挑战。现在，让我们也来看一些可以用来应对这些挑战的方法。

模型的容量与欠拟合或过拟合之间有着密切的关系。一个模型的容量是指它足够灵活，能够拟合各种不同函数的能力。容量较低的模型可能很难拟合训练数据，导致欠拟合。容量较高的模型可能通过过度记忆训练数据而发生过拟合。为了更好地理解这个容量的概念，我们来看一个例子。当我们从线性回归转向多项式回归时，我们正在增加模型的容量。我们不仅仅拟合直线，而是让模型也能拟合曲线。

当机器学习模型的容量与当前学习问题相匹配时，通常表现良好。

图 5.5：欠拟合与过拟合

图 5.5 显示了一个非常常见的例子，用于说明过拟合和欠拟合。我们通过已知函数创建一些随机点，并尝试使用这些数据样本进行学习。我们可以看到，作为最简单模型之一的线性回归，通过在这些点之间画一条直线，未能充分拟合数据。多项式回归是线性回归，但加入了一些更高阶的特征。现在，你可以将从线性回归到多项式回归（增加高阶）的转变，视为增加模型的容量。因此，当我们使用 4 次时，可以看到所学函数很好地拟合了数据，并与我们的理想函数匹配。但是，如果我们继续增加模型的容量，达到 degree = 15 时，我们会看到所学的函数仍然通过训练样本，但它已经学到了一个完全不同的函数，导致了对训练数据的过拟合。找到能够学习出具有良好泛化能力的函数的最佳容量，是机器学习中的核心挑战之一。

容量是模型的一个方面，另一个方面是正则化。即使在相同容量下，模型也可以从所有函数的假设空间中选择多个函数。通过正则化，我们试图在假设空间中对某些函数给予偏好，而非其他函数。

尽管所有这些函数都是有效的，可以选择的函数，我们会通过某种方式推动优化过程，使其最终趋向我们偏好的某种函数。尽管正则化是一个广泛的术语，用于指代我们在学习过程中施加的任何约束，以减少所学函数的复杂性，但更常见的是将其以权重衰减的形式使用。我们以线性回归为例，线性回归是通过学习与每个特征相关联的权重，将一条直线拟合到输入特征上。

线性回归模型可以用数学公式表示如下：

在这里，N 是特征的数量，c 是截距，x[i] 是第 i 个特征，w[i] 是与第 i 个特征相关的权重。我们通过将这个问题视为优化问题，最小化 与 y（真实输出）之间的误差，从而估计出正确的权重 (L)。

现在，通过正则化，我们向 L 添加了一个额外的项，强制权重变得更小。通常，这是通过使用 L1 或 L2 正则化器来完成的。L1 正则化器是将权重的平方和添加到 L 上：

其中， 是正则化系数，决定了我们对权重的惩罚强度。L2 正则化器是将权重的平方和加到 L 上：

在这两种情况下，我们都在强制要求更小的权重优于更大的权重，因为这样可以避免函数过度依赖于机器学习模型中任何一个特征。正则化是一个独立的话题；如果你想了解更多，请前往进一步阅读部分查看一些关于正则化的资源。

另一个有效减少过拟合的方法是简单地使用更多的数据来训练模型。通过使用更大的数据集，模型过拟合的可能性会减少，因为大数据集能够捕获更多的多样性。

现在，我们如何调整参数，以在欠拟合和过拟合之间取得平衡呢？让我们在下一节中详细探讨。

超参数和验证集

几乎所有的机器学习模型都有一些超参数。超参数是模型的参数，它们不是从数据中学习到的，而是在训练开始之前就已经设置好的。例如，正则化的权重就是一个超参数。大多数超参数要么帮助我们控制模型的容量，要么对模型应用正则化。通过控制容量、正则化或两者，我们可以找到在欠拟合和过拟合模型之间的边界，得到一个恰到好处的模型。

但是，由于这些超参数必须在算法外部设置，我们如何估计最佳的超参数呢？虽然它不是核心的学习过程的一部分，但我们也可以从数据中学习超参数。不过，如果我们仅仅使用训练数据来学习超参数，那么它会选择最大的可能模型容量，这会导致过拟合。这就是我们需要验证集的原因，验证集是训练过程中无法访问的部分数据。借用之前的类比，验证集就像学生参加的模拟考试，用来检查他们是否已经学得足够好。但当数据集较小（不是成千上万的样本）时，单一验证集上的表现并不能保证公平的评估。在这种情况下，我们依赖于交叉验证。常见的做法是对原始数据集的不同子集重复进行训练和评估程序。常见的一种方法叫做k 折交叉验证，它将原始数据集分成k个相等的、互不重叠且随机的子集，每个子集在训练其他子集后都会被评估。如果你想了解更多关于交叉验证的技术，我们在进一步阅读部分提供了相关链接。稍后在本书中，我们也会从时间序列的角度讲解这个话题，它与标准的交叉验证方式略有不同。

建议阅读：

尽管本书已经触及了机器学习的皮毛，但还有很多内容，若要更好地理解本书的其余部分，建议更深入地学习机器学习。我们建议从*斯坦福大学的机器学习课程（Andrew Ng）*开始—www.coursera.org/learn/machine-learning。如果你时间紧迫，谷歌的机器学习速成课程也是一个不错的起点—developers.google.com/machine-learning/crash-course/ml-intro。

近些年，机器学习取得了很大进展，伴随着这些进展，能够从数据中学习复杂模式的强大模型也随之出现。当我们将这些模型与经典的时间序列预测模型进行比较时，我们可以看到这些新型模型有着巨大的潜力。但机器学习与时间序列预测之间仍然存在一些根本的差异。在下一节中，我们将了解如何克服这些差异，并使用机器学习进行时间序列预测。

将时间序列预测视为回归

正如我们在第一章《介绍时间序列》中所看到的，时间序列是按时间顺序采集的一组观察值。通常，时间序列预测是关于尝试预测这些观察值在未来将会是什么样。给定一段任意长度的历史观察序列，我们可以预测未来某个时间点的值。

我们看到回归，或者说机器学习用于预测连续变量，是在一组示例数据集上进行的，每个示例由输入特征和目标组成。我们可以看出，回归任务是基于一组输入来预测单一输出，这与预测任务本质上是不同的，预测任务是基于一组历史值来预测未来值。这种时间序列与机器学习回归模型之间的根本不兼容性，就是我们不能直接使用回归来进行时间序列预测的原因。

此外，时间序列预测从定义上来说是一个外推问题，而回归大多数情况下是一个插值问题。外推通常比插值更难通过数据驱动方法来解决。回归问题中的另一个关键假设是训练所使用的样本是独立同分布（iid）的。但时间序列破坏了这个假设，因为时间序列中的后续观察值显示出显著的依赖性。

然而，为了利用机器学习的各种技术，我们需要将时间序列预测转化为回归问题。幸运的是，有方法可以将时间序列转换为回归，并通过引入一些特征来为机器学习模型添加记忆，从而克服 IID 假设。让我们看看如何做到这一点。

时间延迟嵌入

我们在第四章《设定强基准预测》中讨论了 ARIMA 模型，并看到了它是一个自回归模型。我们可以使用相同的概念，将一个时间序列问题转换为回归问题。让我们通过以下图示来明确这一概念：

图 5.6：使用滑动窗口将时间序列转换为回归

假设我们有一个时间序列，包含L个时间步长，就像图示中所示。我们有T作为最新的观测值，T - 1、T - 2，依此类推，随着时间倒退，一直到T - L。在理想的世界中，每个观测值在进行预测时应该以所有先前的观测值为条件。但这是不切实际的，因为L可以非常长。我们通常会限制预测函数，只使用序列中最新的M个观测值，其中M < L。这些被称为有限记忆模型，或马尔科夫模型，而M被称为自回归的阶数、记忆大小或感受野。

因此，在时间延迟嵌入中，我们假设一个任意长度的窗口M < L，并通过将窗口在时间序列的长度上滑动，提取固定长度的子序列。

在图示中，我们采用了一个记忆大小为3的滑动窗口。所以，首先提取的子序列（如果从最新时间点开始，按时间倒序提取）是T – 3、T – 2、T - 1。而T是紧接在该子序列之后的观测值。这将成为数据集中的第一个例子（图示中表格的第1行）。

现在，我们将窗口向左滑动一个时间步长（即倒退到过去），并提取新的子序列，T – 4、T – 3、T - 2。对应的目标将变为T - 1。我们在回到时间序列的开始时重复这一过程，在每一步滑动窗口时，我们都会向数据集中添加一个新的例子。

最终，我们得到了一个对齐的数据集，特征的固定向量大小（即等于窗口大小）和一个单一目标，这就是典型的机器学习数据集的样子。

现在我们有一个包含三个特征的表格，我们也给这三个特征赋予了语义意义。如果我们查看图示中表格的最右列，我们可以看到该列中的时间步长总是比目标落后一个时间步长。我们称之为滞后 1。从右数第二列总是比目标滞后两个时间步长，我们称之为滞后 2。一般化地说，特征中观测值比目标滞后n个时间步长时，我们称之为滞后 n。

通过时间延迟嵌入将时间序列转换为回归模型，能够以一种标准回归框架能够利用的方式，编码时间序列的自回归结构。我们还可以考虑另一种使用回归进行时间序列预测的方法，那就是对时间进行回归。

时间嵌入

如果我们依赖自回归模型中的先前观察值，那么我们在时间嵌入模型中则依赖于时间的概念。核心思想是，我们忽略时间序列的自回归特性，并假设时间序列中的任何值仅仅依赖于时间。我们从与时间序列相关的时间戳中提取能够捕捉时间、时间的流逝、时间的周期性等特征，然后使用这些特征通过回归模型来预测目标值。实现这一点的方法有很多，从简单地对齐一个单调且均匀递增的数值列来捕捉时间的流逝，到使用复杂的傅里叶项来捕捉时间中的周期性成分。我们将在第六章，时间序列预测的特征工程中详细讨论这些技术。

在我们结束本章之前，让我们讨论一个在时间序列预测领域逐渐受到关注的关键概念。本书的大部分内容都采纳了这种新的预测范式。

全球预测模型——范式转变

传统上，每个时间序列都是孤立地处理的。正因为如此，传统的预测方法总是仅仅基于单一时间序列的历史来拟合预测函数。但近年来，由于在当今以数字为主的世界中，收集数据变得更加容易，许多公司开始收集来自相似来源或相关时间序列的大量数据。

例如，零售商如沃尔玛会收集跨千家商店的数百万种产品的销售数据。像 Uber 和 Lyft 这样的公司会收集城市中所有区域的乘车需求。在能源领域，能源消费数据会跨所有消费者进行收集。所有这些时间序列数据集都有共同的行为，因此被称为相关时间序列。

我们可以认为，所有相关时间序列中的时间序列都来自不同的 DGP（数据生成过程），因此可以将它们分别建模。我们称这些为局部预测模型。该方法的另一种替代方式是假设所有时间序列都来自同一个 DGP。我们不为每个时间序列单独拟合预测函数，而是为所有相关时间序列拟合一个单一的预测函数。这种方法在文献中被称为全局或跨学习。大多数现代深度学习模型以及机器学习方法都采用了全局模型的范式。我们将在接下来的章节中详细看到这些内容。

参考检查：

全球一词由David Salinas 等人在DeepAR论文（参考文献1）中提出，跨学习则由Slawek Smyl（参考文献2）提出。

我们之前看到，拥有更多数据将减少过拟合的可能性，因此可以降低泛化误差（训练误差与测试误差之间的差异）。这是局部方法的一个缺点。传统上，时间序列数据通常不长，且在很多情况下，收集更多数据既困难又耗时。在小数据上拟合机器学习模型（具有强大的表达能力）容易导致过拟合。这也是为什么传统上用于预测此类时间序列的时间序列模型会强制施加强先验的原因。但这些限制传统时间序列模型拟合的强先验，也可能导致一种形式的欠拟合，限制了准确性。

强大且具有表达能力的数据驱动模型，如机器学习模型，需要大量数据来生成能够对新数据进行泛化的模型。时间序列本质上是与时间相关的，有时收集更多数据意味着需要等待数月甚至数年，这是不理想的。因此，如果我们无法增加时间序列数据集的长度，我们可以增加时间序列数据集的宽度。如果我们将多个时间序列加入数据集，我们就增加了数据集的宽度，从而增加了模型训练时可用的数据量。

图 5.7 展示了通过视觉化增加时间序列数据集宽度的概念：

图 5.7：时间序列数据集的长度和宽度

这对机器学习模型有利，因为在拟合预测函数时，机器学习模型具有更高的灵活性，且能够利用更多数据，从而能够学习出比传统时间序列模型更复杂的预测函数，而传统时间序列模型通常是通过与相关时间序列共享的方式进行的，完全基于数据驱动。

局部方法的另一个缺点是可扩展性问题。以我们之前提到的沃尔玛为例，需要预测数百万个时间序列，并且无法对所有这些模型进行人工监督。如果从工程角度考虑，在生产系统中训练和维护数百万个模型对任何工程师来说都是一场噩梦。但在全球方法下，我们只需为所有这些时间序列训练一个模型，这大大减少了我们需要维护的模型数量，同时还能够生成所有所需的预测。

这种新的预测范式已经获得了广泛关注，并且在多个时间序列竞赛中持续证明能够改进局部方法，尤其是在相关时间序列数据集上。在 Kaggle 竞赛中，例如Rossman 商店销售（2015 年）、维基百科网页流量时间序列预测（2017 年）、Favorita 公司杂货销售预测（2018 年）以及M5 竞赛（2020 年），获胜的参赛作品都是全球模型——无论是机器学习、深度学习还是二者的结合。Intermarché预测竞赛（2021 年）的获胜作品也是全球模型。有关这些竞赛的链接可以在进一步阅读部分找到。

尽管我们有许多经验性发现表明，全球模型在相关时间序列预测中优于局部模型，但全球模型仍然是一个相对较新的研究领域。Montero-Manson 和 Hyndman (2020) 展示了一些非常有趣的结果，并表明任何局部方法都可以通过具备必要复杂度的全球模型进行逼近，他们提出的最有趣的发现是，即使面对不相关的时间序列，全球模型也能表现得更好。我们将在第十章中进一步讨论全球模型及其策略，即全球预测模型。

参考检查：

Montero-Manson 和 Hyndman (2020) 的研究论文在参考文献中被引用，参考文献编号为3。

总结

我们已经开始了超越基准预测方法的探索，并初步涉足机器学习领域。在简要回顾机器学习的基础知识后，我们了解了诸如过拟合、欠拟合、正则化等关键概念，之后我们看到如何将时间序列预测问题转化为机器学习中的回归问题。我们还对不同的嵌入方法，如时间延迟嵌入和时间嵌入，进行了概念性的理解，这些方法可以用于将时间序列问题转化为回归问题。最后，我们还了解了一种新的时间序列预测范式——全球模型，并在概念层面上将其与局部模型进行了对比。在接下来的几章中，我们将开始实践这些概念，学习特征工程技术和全球模型的策略。

参考文献

以下是我们在本章中使用的参考文献：

1. David Salinas, Valentin Flunkert, Jan Gasthaus, Tim Januschowski (2020). DeepAR: 基于自回归递归网络的概率预测. 国际预测学杂志. 36-3. 1181–1191: doi.org/10.1016/j.ijforecast.2019.07.001

2. Slawek Smyl (2020). 指数平滑与递归神经网络的混合方法用于时间序列预测. 国际预测学杂志. 36-1: 75–85 doi.org/10.1016/j.ijforecast.2019.03.017

3. Pablo Montero-Manso, Rob J Hyndman (2020), Principles and Algorithms for Forecasting Groups of Time Series: Locality and Globality. arXiv:2008.00444[cs.LG]： arxiv.org/abs/2008.00444

进一步阅读

你可以查看以下资源以便进一步阅读：

• Regularization for Sparsity from Google Machine Learning Crash Course： developers.google.com/machine-learning/crash-course/regularization-for-sparsity/l1-regularization

• L1 & L2 Regularization, Inside Bloomberg： www.youtube.com/watch?v=d6XDOS4btck

• Cross-validation: evaluating estimator performance from scikit-learn： scikit-learn.org/stable/modules/cross_validation.html

• Rossmann Store Sales： www.kaggle.com/c/rossmann-store-sales

• Web Traffic Time Series Forecasting： www.kaggle.com/c/web-traffic-time-series-forecasting

• Corporación Favorita Grocery Sales Forecasting： www.kaggle.com/c/favorita-grocery-sales-forecasting

• M5 Forecasting—Accuracy： www.kaggle.com/c/m5-forecasting-accuracy

加入我们的 Discord 社区

加入我们社区的 Discord 空间，与作者和其他读者进行讨论：

packt.link/mts

第六章：时间序列预测中的特征工程

在上一章中，我们开始将机器学习（ML）作为解决时间序列预测问题的工具进行探讨。我们还讨论了一些技术，如时间延迟嵌入和时间嵌入，它们将时间序列预测问题视为从 ML 范式出发的经典回归问题。在本章中，我们将详细介绍这些技术，并通过实践操作，使用我们在本书中一直使用的数据集进行讲解。

在本章中，我们将讨论以下主题：

• 理解特征工程

• 避免数据泄露

• 设置预测视野

• 时间延迟嵌入

• 时间嵌入

技术要求

你需要按照书籍《前言》中的说明，设置Anaconda环境，以便创建一个包含所有代码所需的库和数据集的工作环境。在运行笔记本时，任何额外需要的库都会自动安装。

在使用本章的代码之前，你需要先运行以下笔记本：

• 02-Preprocessing_London_Smart_Meter_Dataset.ipynb 来自第二章

• 01-Setting_up_Experiment_Harness.ipynb 来自第四章

本章的代码可以在github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-2E/tree/main/notebooks/Chapter06找到。

理解特征工程

特征工程顾名思义，是从数据中提取特征的过程，主要通过领域知识来使学习过程更加顺畅高效。在典型的机器学习设置中，工程化良好的特征对于获得优秀的模型性能至关重要。特征工程是机器学习中一个高度主观的部分，每个具体问题都有不同的解决路径——这一路径是为该问题量身定制的。假设你有一个房价数据集，并且有一个特征，建造年份，它告诉你房屋的建造年份。那么，为了优化这些信息，我们可以根据建造年份特征创建另一个特征，房屋年龄。这可能会为模型提供更好的信息，这就是特征工程的应用。

当我们将时间序列问题转化为回归问题时，有一些标准的技术可以应用。这是过程中的关键步骤，因为机器学习模型如何理解时间，取决于我们如何设计特征来捕捉时间。我们在第四章《设定强基准预测》中讨论的基线方法是针对时间序列预测的特定用例设计的，因此问题的时间维度已经内建于这些模型中。例如，ARIMA 模型不需要任何特征工程来理解时间，因为它已经内建于模型中。然而，标准回归模型并没有明确的时间理解，因此我们需要创建良好的特征来嵌入问题的时间维度。

在上一章（第五章，《将时间序列预测作为回归问题》）中，我们讨论了在回归框架中编码时间的两种主要方法：时间延迟嵌入和时间嵌入。虽然我们在高层次上触及了这些概念，但现在是时候深入探讨并看到它们的实际应用了。

笔记本提醒：

要跟随完整代码，请使用01-Feature_Engineering.ipynb笔记本，它位于Chapter06文件夹中。

我们已经将正在使用的数据集拆分为训练集、验证集和测试集。然而，由于我们正在生成基于先前观察的特征，从操作角度来说，当我们将训练集、验证集和测试集结合时会更好。稍后会更清楚为什么，但现在，让我们先相信这一点，继续前进。现在，让我们将这两个数据集合并：

# Reading the missing value imputed and train test split data
train_df = pd.read_parquet(preprocessed / "selected_blocks_train_missing_imputed.parquet")
val_df = pd.read_parquet(preprocessed / "selected_blocks_val_missing_imputed.parquet")
test_df = pd.read_parquet(preprocessed / "selected_blocks_test_missing_imputed.parquet")
#Adding train, validation and test tags to distinguish them before combining
train_df['type'] = "train"
val_df['type'] = "val"
test_df['type'] = "test"
full_df = pd.concat([train_df, val_df, test_df]).sort_values(["LCLid", "timestamp"])
del train_df, test_df, val_df 

现在，我们有了一个full_df，它结合了训练、验证和测试数据集。你们中的一些人可能已经在脑海中敲响了警钟，合并训练集和测试集会带来什么问题呢？那就是数据泄露。让我们来检查一下。

避免数据泄露

数据泄露发生在模型训练时包含了在预测时无法获得的信息。通常，这会导致训练集中的性能很高，但在未见过的数据中表现非常差。数据泄露有两种类型：

• 目标泄露是指关于目标（我们试图预测的内容）的信息泄露到模型中的某些特征中，导致模型过度依赖这些特征，最终导致泛化能力差。这包括以任何方式使用目标的特征。

• 训练-测试污染是指训练集和测试集之间存在一些信息泄露。这可能是由于数据处理不当或拆分数据时的疏忽造成的。但也可能以更微妙的方式发生，例如在拆分训练集和测试集之前对数据集进行缩放。

当我们处理时间序列预测问题时，最大的且最常见的错误是目标泄漏。我们必须对每个特征进行深思熟虑，确保我们不会使用那些在预测时无法获得的数据。以下的图示可以帮助我们记住并内化这个概念：

图 6.1：可用信息与不可用信息，以避免数据泄漏

为了使这个概念在时间序列预测的背景下更加清晰和相关，我们来看一个例子。假设我们正在预测洗发水的销量，而我们使用护发素的销量作为一个特征。我们开发了模型，在训练数据上进行了训练，并在验证数据上进行了测试。模型表现得非常好。然而，当我们开始预测未来时，就会看到一个问题。我们也不知道护发素的未来销量是多少。虽然这个例子比较简单，但有时这种问题并不那么显而易见。这就是为什么我们在创建特征时需要格外小心，并且始终从*这个特征在预测时是否可用？*的角度来评估特征。

最佳实践：

除了对特征进行深思熟虑外，还有很多方法可以识别目标泄漏：

• 如果你构建的模型好得让人难以置信，那么你很可能存在泄漏问题

• 如果任何单一特征在模型的特征重要性中占比过大，那么这个特征可能存在泄漏问题

• 仔细检查与目标高度相关的特征

尽管我们在本书前面已经生成了预测，但我们从未明确讨论过预测视野。这是一个重要的概念，对于我们接下来的讨论至关重要。让我们花点时间来理解预测视野。

设置预测视野

预测视野是指我们在任何时间点希望预测的未来时间步数。例如，如果我们想要预测在过去的电力消耗数据集上接下来的 24 小时的情况，那么预测视野就是 48（因为数据是按半小时记录的）。在第五章，时间序列预测作为回归问题中，我们生成了基准模型，我们一次性预测了所有的测试数据。在这种情况下，预测视野等于测试数据的长度。

在此之前我们从未需要担心这一点，因为在经典的统计预测方法中，这个决策与建模是分开的。如果我们训练了一个模型，就可以用它来预测任何未来的点，而无需重新训练。但是在时间序列回归预测中，我们对预测范围有所限制，这与数据泄漏有关。现在这可能对你来说还不清楚，因此我们将在学习特征工程技术后再回顾这一点。目前，我们只关注单步预测。在我们正在使用的数据集上下文中，这意味着我们将回答的问题是：下一小时的能耗是多少？ 我们将在第四部分，预测的机制 中讨论多步预测和预测的其他机制。

现在我们已经设定了一些基本规则，让我们开始看看不同的特征工程技术。为了跟随 Jupyter notebook 的操作，前往 Chapter06 文件夹并使用 01-Feature_Engineering.ipynb 文件。

时间延迟嵌入

时间延迟嵌入的基本思想是通过最近的观察值将时间嵌入其中。在第五章，时间序列回归预测中，我们讨论了将时间序列的前几次观察作为滞后项（第 5.6 图，位于时间延迟嵌入小节下）。

然而，使用这个概念还有一些其他方法可以捕捉近期和季节性的信息。

• 滞后项

• 滚动窗口聚合

• 季节性滚动窗口聚合

• 指数加权移动平均

让我们来看看。

滞后项或回溯

假设我们有一个时间序列，时间步为 Y[L]。假设当前时间为 T，并且我们有一个历史长度为 L 的时间序列。所以我们的时间序列中，y[T] 是最新的观察值，接着是 y[T-1]、y[T-2] 等等，随着时间的推移向后。正如第五章，时间序列回归预测中所解释的，滞后项是包含时间序列中先前观察值的特征，如下图所示：

图 6.2：滞后特征

我们可以通过包括比当前时间点早的时间步（y[T-a]）来创建多个滞后项；我们将其称为滞后 a。在前面的图示中，我们展示了滞后 1、滞后 2 和 滞后 3。然而，我们可以根据需要添加任意数量的滞后项。现在让我们来学习如何在代码中实现这一点：

df["lag_1"]=df["column"].shift(1) 

还记得我们合并训练集和测试集时，我让你以诚信对待吗？现在是回报这份信任的时候了。如果我们考虑滞后操作（或任何自回归特征），它依赖于时间轴上连续的表示。如果我们考虑测试数据集，对于前几行（或最早的日期），滞后值将会缺失，因为它们属于训练数据集的一部分。因此，通过合并这两个数据集，我们创建了一个沿时间轴的连续表示，其中可以利用 pandas 中的标准函数，如 shift，轻松高效地创建这些特征。

就是这么简单，但我们需要针对每个 LCLid 单独执行滞后操作。我们在 src.feature_engineering.autoregressive_features 中包含了一个名为 add_lags 的有用方法，可以快速高效地为每个 LCLid 添加所有需要的滞后。让我们看看如何使用它。

我们将导入该方法，并使用其中的一些参数来配置我们希望的滞后操作：

from src.feature_engineering.autoregressive_features import add_lags
# Creating first 5 lags and then same 5 lags but from previous day and previous week to capture seasonality
lags = (
    (np.arange(5) + 1).tolist()
    + (np.arange(5) + 46).tolist()
    + (np.arange(5) + (48 * 7) - 2).tolist()
)
full_df, added_features = add_lags(
    full_df, lags=lags, column="energy_consumption", ts_id="LCLid", use_32_bit=True
) 

现在，让我们来看一下我们在前面代码片段中使用的参数：

• lags：该参数接收一个整数列表，表示我们需要创建为特征的所有滞后。

• column：要进行滞后操作的列名。在我们的例子中，这是 energy_consumption。

• ts_id：包含时间序列唯一 ID 的列名。如果为 None，则假定数据框仅包含单一的时间序列。在我们的例子中，LCLid 就是该列的名称。

• use_32_bit：该参数在功能上没有任何作用，但通过牺牲浮动点数的精度，使得数据框在内存中的大小变得更小。

该方法返回添加滞后后的数据框（DataFrame），以及包含新添加特征的列名的列表。

滚动窗口聚合

使用滞后时，我们将当前点与过去的单个点连接，而使用滚动窗口特征时，我们将当前点与过去窗口的汇总统计量连接。我们不是查看前几个时间步的观察值，而是查看过去三个时间步的观察值的平均值。请查看下面的图表，以更好地理解这一点：

图 6.3：滚动窗口聚合特征

我们可以使用不同的窗口来计算滚动统计量，每个窗口将捕捉历史的略有不同的方面。在前面的图表中，我们可以看到一个窗口为三和一个窗口为四的示例。当我们处于时间步* T *时，窗口大小为三的滚动窗口将具有 y[T] [– 3]，y[T] [– 2]，y[T] [– 1] 作为过去观察值的向量。一旦我们拥有这些数据，就可以应用任何聚合函数，如均值、标准差、最小值、最大值等。通过聚合函数得到标量值后，我们可以将其作为时间步 t 的特征。

我们不包括 y[T]在过去观测值的向量中，因为那会导致数据泄漏。

让我们看看如何通过 pandas 来实现这一操作：

# We shift by one to make sure there is no data leakage
df["rolling_3_mean"] = df["column"].shift(1).rolling(3).mean() 

与滞后类似，我们需要为每个LCLid列分别执行此操作。我们在src.feature_engineering.autoregressive_features中提供了一个有用的方法add_rolling_features，它可以快速有效地为每个LCLid添加所需的所有滚动特征。让我们看看如何使用它。

我们将导入这个方法，并使用它的一些参数来按照我们希望的方式配置滚动操作：

from src.feature_engineering.autoregressive_features import add_rolling_features
full_df, added_features = add_rolling_features(
    full_df,
    rolls=[3, 6, 12, 48],
    column="energy_consumption",
    agg_funcs=["mean", "std"],
    ts_id="LCLid",
    use_32_bit=True,
) 

现在，让我们来看一下在前面代码片段中使用的参数：

• rolls: 该参数接受一个整数列表，表示我们需要计算聚合统计量的所有窗口。

• column: 要进行滞后操作的列名。在我们的案例中，这列是energy_consumption。

• agg_funcs: 这是一个聚合函数列表，用于对rolls中声明的每个窗口进行处理。允许的聚合函数包括{mean, std, max, min}。

• n_shift: 这是在进行滚动操作之前需要移动的时间步数。此参数可以避免数据泄漏。虽然我们在这里移动了一个时间步，但也有需要移动多个时间步的情况。通常这用于多步预测，我们将在第四部分，预测的机制中讨论。

• ts_id: 包含时间序列唯一 ID 的列名。如果为None，则假定数据框架只有一个时间序列。在我们的案例中，LCLid就是该列的名称。

• use_32_bit: 该参数在功能上没有任何作用，但能使数据框架在内存中占用更小的空间，牺牲浮动点数的精度。

该方法返回添加了滚动特征的数据框架，并且返回一个包含新添加特征列名的列表。

季节性滚动窗口聚合

季节性滚动窗口聚合与滚动窗口聚合非常相似，但它们不同的是，季节性窗口不会采用过去的n个连续观测值，而是采用一个季节性窗口，在窗口中的每个项之间跳过一个恒定数量的时间步。下面的图表将使这一点更加清晰：

图 6.4：季节性滚动窗口聚合

这里的关键参数是季节性周期，通常称为 M。这是我们预计季节性模式会重复的时间步数。在时间步 T 时，大小为三的滚动窗口将包含 y[T] [– 3]、y[T] [– 2]、y[T] [– 1]，作为过去观察值的向量。但是，季节性滚动窗口将在窗口中的每个元素之间跳过 m 个时间步数。这意味着季节性滚动窗口中的观察值将是 y[T] [–] [M]、y[T] [– 2][M]、y[T] [– 3][M]。此外，像往常一样，一旦我们获得窗口向量，我们只需应用聚合函数以获得标量值，并将其作为特征包含在内。

我们 不包括 y[T] 作为季节性滚动窗口向量中的一个元素，以避免数据泄漏。

这是一个你不能使用 pandas 轻松高效完成的操作。需要一些高级的 NumPy 索引和 Python 循环来实现这个功能。我们将使用来自 github.com/jmoralez/window_ops/ 的实现，它利用 NumPy 和 Numba 来使操作更快速高效。

就像我们之前看到的特征一样，我们需要为每个 LCLid 单独执行此操作。我们在 src.feature_engineering.autoregressive_features 中提供了一个有用的方法 add_seasonal_rolling_features，它可以快速有效地为每个 LCLid 添加你需要的所有季节性滚动特征。让我们看看如何使用它。

我们将导入该方法并使用该方法的几个参数来配置我们想要的季节性滚动操作：

from src.feature_engineering.autoregressive_features import add_seasonal_rolling_features
full_df, added_features = add_seasonal_rolling_features(
    full_df,
    rolls=[3],
    seasonal_periods=[48, 48 * 7],
    column="energy_consumption",
    agg_funcs=["mean", "std"],
    ts_id="LCLid",
    use_32_bit=True,
) 

现在，让我们看看在之前的代码片段中使用的参数：

• seasonal_periods：这是一个季节性周期的列表，应在季节性滚动窗口中使用。在多季节性的情况下，我们可以包含所有季节性滚动特征。

• rolls：该参数接受一个整数列表，表示我们需要计算聚合统计量的所有窗口。

• column：要滞后的列名。在我们的例子中，这是 energy_consumption。

• agg_funcs：这是我们希望对在 rolls 中声明的每个窗口进行的聚合操作的列表。允许的聚合函数为 {mean, std, max, min}。

• n_shift：这是在进行滚动操作之前需要移动的季节性时间步数。该参数可以防止数据泄漏。

• ts_id：包含时间序列唯一 ID 的列名。如果为 None，则假设 DataFrame 仅包含单个时间序列。在我们的例子中，LCLid 就是该列名。

• Use_32_bit：该参数在功能上没有任何作用，但通过牺牲浮动点数的精度，使 DataFrame 在内存中变得更小。

和往常一样，该方法返回包含季节性滚动特征的 DataFrame 以及一个包含新添加特征列名的列表。

指数加权移动平均（EWMA）

在滚动窗口均值操作中，我们计算了窗口的平均值，这与移动平均是同义的。EWMA 是移动平均的稍微聪明一些的“亲戚”。移动平均考虑了一个滚动窗口，并且对窗口中的每个项在计算的平均值中赋予相同的权重，而 EWMA 则尝试对窗口进行加权平均，并且权重按指数速率衰减。有一个参数，，决定了权重衰减的速度。因为这个原因，我们可以将所有可用的历史数据视为一个窗口，让参数决定 EWMA 中包括多少近期的数据。这个过程可以简单地递归表示如下：

在这里，我们可以看到，的值越大，平均值越偏向于最近的值（查看图 6.6以获得关于权重如何分布的直观印象）。如果我们展开递归，每一项的权重将会是：

其中，k是相对于T的时间步长。如果我们绘制权重，就能看到它们呈指数衰减；决定了衰减的速度。另一种理解的方式是从跨度的角度来看。跨度是指衰减后的权重接近零的周期数（这不是严格的数学意义，而是直观的理解）。和跨度通过以下方程相关：

这一点将在下图中更清晰地展示，我们绘制了不同值下权重衰减的情况 ：

图 6.5：不同值下的指数权重衰减

在这里，我们可以看到，当我们达到跨度时，权重变得很小。

直观地，我们可以将 EWMA 视为整个时间序列历史的平均值，但通过和跨度等参数，我们可以使不同时期的历史对平均值的影响更加显著。如果我们定义了一个 60 期的跨度，我们可以认为最后的 60 个时间期主要驱动了平均值。因此，通过使用不同跨度或值的 EWMA，我们可以得到具有代表性的特征，捕捉不同历史时期的特征。

整个过程在下图中进行了展示：

图 6.6：EWMA 特征

现在，让我们看看如何在 pandas 中实现这一点：

df["ewma"]=df['column'].shift(1).ewm(alpha=0.5).mean() 

和我们之前讨论的其他特性一样，EWMA 也需要对每个LCLid单独进行处理。我们在src.feature_engineering.autoregressive_features中包含了一个有用的方法，叫做add_ewma，它可以快速高效地为每个LCLid添加你所需要的所有 EWMA 特征。让我们看看如何使用它。

我们将导入这个方法，并使用该方法的几个参数来配置我们想要的 EWMA：

from src.feature_engineering.autoregressive_features import add_ewma
full_df, added_features = add_ewma(
    full_df,
    spans=[48 * 60, 48 * 7, 48],
    column="energy_consumption",
    ts_id="LCLid",
    use_32_bit=True,
) 

现在，让我们看一下在前面的代码片段中使用的参数：

• alphas：这是我们需要计算 EWMA 特征的所有的列表。

• spans：作为替代，我们可以用它列出我们需要计算 EWMA 特征的所有跨度。如果使用此特性，alphas将被忽略。

• column：需要进行滞后处理的列名。在我们的案例中，这一列是energy_consumption。

• n_shift：这是我们在进行滚动操作之前需要移动的季节性时间步数。这个参数可以避免数据泄漏。

• ts_id：包含时间序列唯一 ID 的列名。如果为None，则假设数据框仅包含单一时间序列。在我们的案例中，LCLid是该列的名称。

• use_32_bit：这个参数在功能上没有任何作用，但可以使数据框在内存中占用更少的空间，代价是牺牲浮点数的精度。

一如既往，该方法返回包含 EWMA 特征的数据框，并返回一个包含新添加特征列名的列表。

这些是将时间延迟嵌入到机器学习模型中的几种标准方式，但你并不局限于这些。像往常一样，特征工程是一个没有固定规则的领域，我们可以尽情发挥创意，并将领域知识注入到模型中。除了我们已经看到的特征，我们还可以包括滞后差异作为自定义滞后特征来注入领域知识，等等。在大多数实际情况下，我们最终会使用不止一种方式将时间延迟嵌入到模型中。滞后特征在大多数情况下是最基本和最重要的，但我们确实会通过季节性滞后、滚动特征等编码更多信息。和机器学习中的一切一样，没有万能的解决方案。每个数据集都有其独特性，这使得特征工程对于每种情况都非常重要且不同。

现在，让我们看一下我们可以通过时间嵌入添加的另一类特征。

时间嵌入

在第五章，时间序列预测作为回归问题中，我们简要讨论了时间嵌入的过程，即我们尝试将时间嵌入到机器学习模型可以利用的特征中。如果我们稍微思考一下时间，我们会发现，在时间序列预测的背景下，时间有两个方面对我们来说至关重要——时间的流逝和时间的周期性。

有一些特征可以帮助我们在机器学习模型中捕捉这些方面：

• 日历特征

• 时间流逝

• 傅里叶项

让我们逐一看一下它们。

日历特征

我们可以提取的第一类特征是基于日历的特征。尽管时间序列的严格定义是按时间顺序获取的一组观测值，但我们通常会在这些观测值的时间戳旁边收集时间序列。我们可以利用这些时间戳并提取日历特征，例如月份、季度、年中的第几天、小时、分钟等。这些特征捕捉了时间的周期性，帮助机器学习模型有效地捕捉季节性。只有比时间序列频率更高的日历特征才有意义。例如，在一个每周频率的时间序列中，小时特征是没有意义的，但月份和周数特征则有意义。我们可以利用 pandas 中的内置日期时间功能来创建这些特征，并在模型中将它们视为分类特征。

时间流逝

这是另一个在机器学习模型中捕捉时间流逝的特征。随着时间的推移，该特征单调增加，给机器学习模型提供时间流逝的感知。创建这个特征有很多方法，但最简单且高效的方法之一是使用 NumPy 中日期的整数表示：

df['time_elapsed'] = df['timestamp'].values.astype(np.int64)/(10**9) 

我们在 src.feature_engineering.temporal_features 中包含了一个有用的方法 add_temporal_features，该方法会自动添加所有相关的时间特征。让我们看看如何使用它。

我们将导入该方法，并使用该方法的一些参数来配置和创建时间特征：

full_df, added_features = add_temporal_features(
    full_df,
    field_name="timestamp",
    frequency="30min",
    add_elapsed=True,
    drop=False,
    use_32_bit=True,
) 

现在，让我们看看在前面的代码片段中使用的参数：

• field_name：这是包含应当用于创建特征的日期时间的列名。

• frequency：我们应当提供时间序列的频率作为输入，以便方法自动提取相关特征。这些是标准的 pandas 频率字符串。

• add_elapsed：此标志用于开启或关闭时间流逝特征的创建。

• use_32_bit：该参数在功能上没有任何作用，但使得 DataFrame 在内存中占用更小的空间，牺牲了浮点数的精度。

就像我们讨论的前几种方法一样，这个方法也会返回一个新的 DataFrame，添加了时间特征，并返回一个包含新增特征列名的列表。

傅里叶项

之前，我们提取了一些日历特征，如月份、年份等，并讨论了将它们作为分类变量用于机器学习模型。另一种表示相同信息的方式，是使用 傅里叶项，并以连续尺度表示。我们在第三章《分析与可视化时间序列数据》中讨论了傅里叶级数。为了重申，傅里叶级数的正弦-余弦形式如下：

在这里，S[N]是信号S的N项近似。理论上，当N为无限时，得到的近似值等于原始信号。P是周期的最大长度，a[n]和b[n]分别是余弦项和正弦项的系数，a[0]是截距。

我们可以将这些余弦和正弦函数作为特征来表示季节性循环。如果我们对月份进行编码，我们知道它的范围从 1 到 12，然后会重复。因此，在这种情况下，P将是 12，x将是 1, 2, …12。因此，对于每个x，我们可以计算出余弦和正弦项，并将它们作为特征添加到机器学习模型中。从直观上看，我们可以认为模型会根据数据推断出系数，从而帮助模型更容易地预测时间序列。

下图展示了按序数尺度表示的月份与作为傅里叶级数表示之间的差异：

图 6.7：作为序数步进函数的月份（上图）与傅里叶项（下图）

上图展示了单一傅里叶项；我们可以添加多个傅里叶项来帮助捕捉复杂的季节性。

我们不能简单地说季节性的连续表示优于类别表示，因为这取决于你使用的模型类型和数据集。这是我们需要通过经验来发现的。

为了简化添加傅里叶特征的过程，我们提供了一些易于使用的方法，这些方法位于src.feature_engineering.temporal_features中，文件名为bulk_add_fourier_features，它会自动为我们需要的所有日历特征添加傅里叶特征。让我们来看看如何使用它。

我们将导入该方法，并使用其中的一些参数来配置和创建基于傅里叶级数的特征：

full_df, added_features = bulk_add_fourier_features(
    full_df,
    ["timestamp_Month", "timestamp_Hour", "timestamp_Minute"],
    max_values=[12, 24, 60],
    n_fourier_terms=5,
    use_32_bit=True,
) 

现在，让我们来看一下我们在前面的代码片段中使用的参数：

• columns_to_encode：这是我们需要使用傅里叶项进行编码的日历特征列表。

• max_values：这是一个季节性循环的最大值列表，对于日历特征，按columns_to_encode中给出的顺序排列。例如，对于month要作为列进行编码时，我们给出12作为对应的max_value。如果没有给出，max_value将会被推断。只有在你的数据至少包含一个完整的季节性循环时，才建议使用此方法。

• n_fourier_terms：这是要添加的傅里叶项数量。这与前面提到的傅里叶级数公式中的n是同义词。

• use_32_bit：这个参数在功能上没有任何作用，但它会使 DataFrame 在内存中变得更小，从而牺牲浮点数的精度。

就像我们之前讨论过的方法一样，它也会返回一个新的 DataFrame，添加了傅里叶特征，同时返回一个包含新增特征列名的列表。

在执行Chapter06中的01-Feature_Engineering.ipynb笔记本后，我们将生成以下特征工程文件并保存到磁盘：

• selected_blocks_train_missing_imputed_feature_engg.parquet

• selected_blocks_val_missing_imputed_feature_engg.parquet

• selected_blocks_test_missing_imputed_feature_engg.parquet

在本节中，我们探讨了一些流行且有效的时间序列特征生成方法。但实际上还有很多其他方法，具体选择哪些方法取决于你的问题和领域，很多方法都将适用。

附加信息：

特征工程的领域非常广泛，已有一些开源库使得探索这一领域变得更加容易。其中一些库包括github.com/Nixtla/tsfeatures、tsfresh.readthedocs.io/en/latest/以及github.com/DynamicsAndNeuralSystems/catch22。Ben D. Fulcher 撰写的预印本论文《基于特征的时间序列分析》(arxiv.org/abs/1709.08055)也提供了对这一领域的良好总结。

一个名为 functime 的新库（github.com/functime-org/functime）也提供了快速的特征工程例程，它是用 Polars 编写的，值得一试。书中讨论的许多特征工程方法，通过使用 functime 和 Polars 可以使处理速度更快。

总结

在上一章简要概述了时间序列预测的机器学习范式后，本章我们从实践角度深入探讨，了解如何准备带有所需特征的数据集，以开始使用这些模型。我们回顾了几种时间序列特定的特征工程技术，如滞后、滚动窗口和季节性特征。本章中学习的所有技术都是我们可以快速迭代实验的工具，帮助我们找到最适合我们数据集的方案。然而，我们仅讨论了特征工程，它只影响标准回归方程的一部分（y = mX + c）。另一个部分，即我们预测的目标（y），同样重要。下一章我们将探讨一些概念，比如平稳性以及一些影响目标的转换。

加入我们在 Discord 上的社区

加入我们社区的 Discord 空间，与作者和其他读者进行讨论：

packt.link/mts

第七章：时间序列预测的目标转换

在上一章中，我们探讨了如何通过利用特征工程技术进行时间嵌入和时间延迟嵌入。但那只是回归方程的一方面——特征。通常，我们会发现方程的另一侧——目标，并没有按照我们预期的方式表现出来。换句话说，目标并没有具备一些理想的特性，这些特性能够让预测变得更加容易。这个领域的主要罪魁祸首之一是平稳性——或者更具体地说，缺乏平稳性。这会给我们在开发机器学习（ML）/统计模型时所做的假设带来问题。在本章中，我们将讨论一些处理目标相关问题的技术。

本章将覆盖以下主题：

• 处理时间序列中的非平稳性

• 检测和修正单位根

• 检测和修正趋势

• 检测和修正季节性

• 检测和修正异方差性

• AutoML 方法用于目标转换

技术要求

你需要按照本书前言中的说明设置Anaconda环境，以便获得一个包含所有所需库和数据集的工作环境，供本书代码使用。任何额外的库将在运行笔记本时自动安装。

在使用本章中的代码之前，你需要运行以下笔记本：

• Chapter02文件夹中的02-Preprocessing_London_Smart_Meter_Dataset.ipynb

• Chapter04文件夹中的01-Setting_up_Experiment_Harness.ipynb

• Chapter06文件夹中的01-Feature_Engineering.ipynb

本章的相关代码可以在github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-/tree/main/notebooks/Chapter07找到。

检测时间序列中的非平稳性

平稳性是大多数计量经济学模型中普遍的假设，是一个严谨且数学化的概念。但在不涉及太多数学的情况下，我们可以直观地理解平稳性为时间序列所采样的分布的统计特性随着时间保持不变的状态。这在时间序列回归中也很相关，因为我们在跨时间估计一个单一的预测函数。如果时间序列的行为随时间发生变化，那么我们估计的单一函数可能并不总是适用。例如，如果我们把附近公园每日的游客数量看作一个时间序列，我们知道这些模式在疫情前后会有很大不同。在机器学习领域，这种现象被称为概念漂移。

直观上，我们可以理解，平稳的序列比非平稳的序列更容易预测。但这里有个关键点：在现实世界中，几乎所有的时间序列都不满足平稳假设——更具体地说，严格平稳假设。严格平稳是指所有统计性质，如均值、方差、偏度等，都不随时间变化。很多时候，这一严格平稳假设被放宽，采用弱平稳，即只要求时间序列的均值和方差不随时间变化。

我们可以问自己四个主要问题来检查我们的时间序列是否平稳：

• 均值是否随时间变化？换句话说，时间序列中是否有趋势？

• 方差是否随时间变化？换句话说，时间序列是否异方差？

• 时间序列的均值是否表现出周期性的变化？换句话说，时间序列中是否存在季节性？

• 时间序列是否具有单位根？

其中前三个问题可以通过简单的视觉检查来确定。单位根较难理解，我们将稍后深入探讨单位根。现在，让我们看几个时间序列，检查我们是否能通过视觉检查来判断它们是否平稳（你可以记录你的答案）：

图 7.1：测试你对平稳性的理解

现在，检查你的回答，看看你猜对了多少个。如果你答对了六个中的至少四个，你对平稳性的直觉理解就非常好：

• 时间序列 1是平稳的，因为它是一个白噪声过程，按照定义，它的均值为零，方差恒定。它通过了我们对前三个问题的检查。

• 时间序列 2是非平稳的，因为它有一个明显的下降趋势。这意味着序列在开始时的均值与结束时不同。因此，它没有通过我们的第一个检查问题。

• 时间序列 3乍一看可能是平稳的，因为它基本上围绕 0 波动，但随着时间的推移，波动幅度逐渐增大。这意味着它的方差在增大——换句话说，它是异方差的。因此，尽管它回答了我们的第一个问题，但未能通过第二个检查，即方差是否恒定。所以，它是非平稳的。

• 现在，我们要讨论的难题是——时间序列 4。乍一看，我们可能认为它是平稳的，因为尽管它一开始有趋势，但它也发生了反转，使得均值几乎保持不变。而且，方差的变化也不那么明显。但这实际上是一个含有单位根的时间序列（我们将在稍后的单位根部分详细讨论），通常，含单位根的时间序列很难通过视觉检查来判断。

• 时间序列 5 回答了前两个问题——常数均值和常数方差——但它有一个非常明显的季节性模式，因此是非平稳的。

• 时间序列 6 是另一个白噪声过程，仅用于迷惑你。这也是平稳的。

当我们拥有成百上千，甚至是数百万个时间序列时，我们在实际操作中无法通过视觉检查来判断它们是否平稳。因此，现在让我们来看看使用统计测试检测这些关键特性的一些方法，以及如何尝试修正它们。

虽然我们讨论的是修正或使时间序列平稳，但在机器学习范式下，并不总是必须这样做，因为一些问题可以通过使用合适的特征在模型中得到处理。是否将一个序列转换为平稳序列是一个我们需要在实验技术后做出的决定。因为正如你将看到的，虽然将序列平稳化有其优点，但使用这些技术也有其缺点，我们将在详细讨论每种转换时看到这些问题。

笔记本提醒：

若要跟随完整的代码，请使用 Chapter06 文件夹中的 02-Dealing_with_Non-Stationarity.ipynb 笔记本。

检测和修正单位根

让我们先谈谈单位根，因为这是最常见的平稳性测试。时间序列分析起源于计量经济学和统计学，而单位根是直接源于这些领域的一个概念。

单位根

单位根相当复杂，但为了培养一些直觉，我们可以简化来看。让我们考虑一个一阶自回归模型（AR(1) 模型）：

，其中 是白噪声， 是 AR 系数。

如果我们考虑方程中不同的 值，我们可以提出三种情境（图 7.2）：

1. ：当 大于 1 时，时间序列中的每个后续值都会乘以大于 1 的数字，这意味着它将具有强烈且迅速增加/减少的趋势，因此是非平稳的。

2. ：当 小于 1 时，时间序列中的每个后续值都会乘以小于 1 的数字，这意味着从长远来看，序列的均值趋向于零，并围绕零振荡。因此，它是平稳的。

3. ：当 等于 1 时，情况变得更加复杂。当 为 AR(1) 模型时，这被称为具有单位根，方程变为 。这在计量经济学中被称为随机漫步，是金融和经济领域中一种非常流行的时间序列类型。从数学上讲，我们可以证明这样的序列将具有常数均值，但方差却是非恒定的：

图 7.2：具有不同 参数的自回归时间序列。顶部：情境 1，phi <1；中间：情境 2，phi = 1；底部：情境 3，phi>1

虽然我们讨论了AR(1)过程中的单位根问题，但我们可以将相同的直觉扩展到多个滞后项或AR(p)模型。在那里，计算和检验单位根更加复杂，但仍然是可能的。

既然我们已经知道什么是单位根，如何通过统计方法进行检验呢？这时 Dickey-Fuller 检验就派上用场了。

扩展的 Dickey-Fuller（ADF）检验

这个检验中的零假设是AR(1)模型中的 等于 1，从而使序列非平稳。备择假设是AR(1)模型中的 小于 1。ADF 检验将 Dickey-Fuller 检验扩展到AR(p)模型，因为大多数时间序列不仅仅由一个滞后项定义。这个检验是检查单位根的标准且最流行的统计方法。该检验的核心是对时间序列的滞后进行回归，并计算残差方差的统计量。

让我们看看如何在 Python 中使用statsmodels来实现这一点：

from statsmodels.tsa.stattools import adfuller
result = adfuller(y) 

adfuller的result是一个元组，包含了检验统计量、p 值和不同置信度下的临界值。在这里，我们最关注的是 p 值，这是一个简单实用的方法来检查零假设是否被拒绝。如果p < 0.05，则有 95%的概率表明序列没有单位根；从单位根的角度来看，序列是平稳的。

为了让这个过程更加简便，我们在src.transforms.stationary_utils中包含了一个名为check_unit_root的方法，它会为你进行推断（通过将返回的概率与置信度进行比较并拒绝或接受零假设），并返回一个包含布尔属性stationary的namedtuple，以及来自statsmodels的完整结果在results中：

from src.transforms.stationary_utils import check_unit_root
# We pass the time series along with the confidence with which we need the results
check_unit_root(y, confidence=0.05) 

现在我们已经学会了如何检查一个序列是否有单位根，那么如何使其平稳呢？我们来看看几个可以帮助我们实现这一目标的变换方法。

差分变换

差分变换是使时间序列平稳，或至少消除单位根的非常流行的变换方法。这个概念很简单：我们将时间序列从观察领域转化为观察变化领域。差分变换通过将后续观察值相互减去来实现：

z[t] = y[t] – y[t][-1]

差分有助于稳定时间序列的均值，并因此减少或消除趋势和季节性。让我们看看差分是如何使序列变得平稳的。

让我们考虑问题中的时间序列 ，其中 和 是系数， 是白噪声。从这个方程中，我们可以看到时间 t 是方程的一部分，这使得 y[t] 成为一个具有趋势的时间序列。因此，差分后的时间序列 z 将如下所示：

我们需要在这个新方程中寻找的是没有提到 t。这意味着，创建趋势的 t 依赖性已经被去除，现在时间序列在任何时间点都有恒定的均值和方差。

差分并不能去除所有种类的非平稳性，但适用于大多数时间序列。然而，这种方法也有一些缺点。其中一个缺点是我们在建模过程中会丢失时间序列的尺度。很多时候，时间序列的尺度包含一些对预测有用的信息。例如，在供应链中，高销售量的 SKU 表现出不同于低销售量 SKU 的模式，而当我们进行差分时，这种区分信息会丢失。

另一个缺点更多是从操作角度来看。当我们使用差分进行预测时，我们还需要在从模型获取差分后的输出后对其进行逆变换。这是我们必须管理的一个额外复杂性。解决方法之一是将最新的观察值保存在内存中，并不断将差分添加到其中以逆变换。另一种方法是为每个需要逆变换的 t 保留 y[t][-1]，并将差分不断添加到 y[t][-1] 中。

我们已经使用 datetime 索引作为关键字来对齐并获取src.transforms.target_transformations.py中* y *[t][-1]的观察值，该文件位于本书的 GitHub 仓库中。让我们看看如何使用它：

from src.transforms.target_transformations import AdditiveDifferencingTransformer
diff_transformer = AdditiveDifferencingTransformer()
# [1:] because differencing reduces the length of the time series by one
y_diff = diff_transformer.fit_transform(y, freq="1D")[1:] 

y_diff将包含变换后的序列。要返回原始时间序列，我们可以调用inverse_transform并使用diff_transformer。关联的笔记本有示例和图表，展示差分如何改变时间序列。

在这里，我们将差分看作是从时间序列中减去后续值的过程。但我们也可以使用其他运算符进行差分，例如除法（y[t]/ y[t][-1]），这在src.transforms.target_transformations.py文件中实现为MultiplicativeDifferencingTransformer。我们还可以尝试这些变换，以检查它们是否对你的数据集最有效。

尽管差分解决了大多数平稳性问题，但并不能保证解决所有类型的趋势（非线性或分段趋势）、季节性等问题。有时候，我们可能不想对序列进行差分，但仍然需要处理趋势和季节性。那么，让我们看看如何检测并去除时间序列中的趋势。

检测并纠正趋势

在第五章，《作为回归的时间序列预测》中，我们讨论了预测是一个困难的问题，因为它本质上是一个外推问题。趋势是导致预测成为外推问题的主要因素之一。如果我们有一个上升趋势的时间序列，任何试图预测它的模型都需要在训练过程中所看到的值范围之外进行外推。ARIMA 通过自回归来处理这个问题，而指数平滑通过显式建模趋势来处理它。但标准回归可能并不天然适合外推。然而，通过合适的特征，如滞后，它也可以开始进行外推。

但是，如果我们能自信地估计并提取时间序列中的趋势，我们可以通过去趋势化时间序列来简化回归问题。

但在继续之前，值得了解两种主要的趋势类型。

确定性和随机性趋势

让我们再次使用之前看到的简单*AR(1)模型来培养对这个问题的直觉。之前我们看到在AR(1)*模型中有时，时间序列中会出现趋势。但我们也可以通过将时间作为有序变量包含在定义时间序列的方程中来考虑趋势时间序列。例如，让我们考虑两个时间序列：

时间序列 1：

时间序列 2：

这些可以通过以下图表看到：

图 7.3：顶部：随机趋势；底部：确定性趋势

我们之前看到过这两个方程；时间序列 1是AR(1)模型，而时间序列 2是我们选择用来说明差分的时间序列方程。我们已经知道，对于，时间序列 1和时间序列 2都有趋势。但这两种趋势之间有差异。

在时间序列 2中，趋势是恒定的，可以完美地建模。在这种情况下，简单的线性拟合就能完美地解释趋势。但在时间序列 1中，趋势并不是可以通过简单的线性拟合来解释的。它本质上依赖于时间序列的前一个值，因此是随机的。因此，时间序列 2具有确定性趋势，而时间序列 1具有随机性趋势。

我们可以使用本章中早些时候看到的相同 ADF 检验来检查时间序列是否具有确定性或随机性趋势。在不深入讨论统计检验的数学原理的情况下，我们知道它通过将AR(p)模型拟合到时间序列上来检验单位根。我们可以通过statsmodels实现中的regression参数指定该检验的几种变体。该参数接受以下值：

• c：这意味着我们在AR(p)模型中包括了常数截距项。实际上，这意味着即使序列不以零为中心，我们也会认为时间序列是平稳的。这是statsmodels中的默认设置。

• n：这意味着我们在AR(p)模型中甚至不包括常数截距项。

• ct：如果我们提供此选项，则AR(p)模型将包含常数截距项和线性、确定性趋势分量。这意味着即使时间序列中存在确定性趋势，它也会被忽略，并且该序列将被测试为平稳序列。

• ctt：这是当我们包括常数截距项，即线性和二次趋势时。

因此，如果我们使用regression="c"运行 ADF 检验，它将显示为非平稳序列。而如果我们使用regression="ct"运行 ADF 检验，它将显示为平稳序列。这意味着，当我们从时间序列中移除确定性趋势后，它变得平稳。这项测试可以帮助我们判断时间序列中观察到的趋势是确定性趋势还是随机趋势。在进一步阅读部分，我们提供了一个由Fabian Kostadinov撰写的博客链接，他通过对一些时间序列进行实验，清楚地说明了不同 ADF 检验变种的区别。

我们已经在src.transforms.stationary_utils中实现了这个测试，命名为check_deterministic_trend，它会为您进行推断，并返回一个包含布尔属性deterministic_trend的namedtuple。如果您想进一步调查，namedtuple还包含我们在adf_res和adf_ct_res中进行的两个adfuller检验的原始结果。让我们看看如何使用这个测试：

check_deterministic_trend(y, confidence=0.05) 

这将告诉我们趋势是平稳的还是确定性的。接下来，让我们看看几种方法，如何识别并对时间序列中的趋势（无论是否为确定性趋势）进行统计检验。

Kendall 的 Tau

Kendall 的 Tau 是相关性的一种衡量方式，但它是在数据的等级上进行的。Kendall 的 Tau 是一个非参数检验，因此不对数据做任何假设。相关系数 Tau 的值介于-1 和 1 之间，其中 0 表示无关系，1 或-1 表示完全关系。我们不会深入探讨 Kendall 的 Tau 是如何计算的，也不会讨论其显著性检验，因为这超出了本书的范围。进一步阅读部分包含了一个很好的链接来解释这一点。

在本节中，我们将看到如何使用 Kendall 的 Tau 来衡量我们时间序列中的趋势。如前所述，Kendall 的 Tau 计算两个变量之间的等级相关性。如果我们选择其中一个变量作为时间序列，并将另一个变量设置为时间的序数表示，则得到的 Kendall 的 Tau 将表示时间序列中的趋势。另一个好处是，Kendall 的 Tau 值越高，我们预计趋势越强。

scipy 提供了一个 Kendall’s Tau 的实现，我们可以如下使用：

import scipy.stats as stats
tau, p_value = stats.kendalltau(y, np.arange(len(y))) 

我们可以将返回的 p 值与所需的置信度（通常为 0.05）进行比较，并得出结论：如果 p_value < confidence，我们可以认为趋势在统计上是显著的。tau 的符号告诉我们这是一个上升趋势还是下降趋势。

我们在 src.transforms.stationary_utils 中实现了 Kendall’s Tau，命名为 check_trend，它会帮助你检查是否存在趋势。我们只需要提供以下参数：

• y：需要检查的时间序列

• confidence：对结果 p 值进行检查的置信度水平

还有一些参数，但它们是用于 Mann-Kendall（M-K）检验的，接下来会做详细解释。

让我们看看如何使用这个检验：

check_trend(y, confidence=0.05) 

此方法还检查已识别的趋势是确定性还是随机性，并计算趋势的方向。结果会作为 namedtuple 返回，包含以下参数：

• trend：这是一个布尔标志，表示是否存在趋势。

• direction：这将是 increasing 或 decreasing。

• slope：这是估算趋势线的斜率。对于 Kendall’s Tau，它将是 Tau 值。

• p：这是统计检验的 p 值。

• deterministic：这是一个布尔标志，表示确定性趋势。

现在，让我们来看看 Mann-Kendall 检验。

Mann-Kendall 检验（M-K 检验）

Mann-Kendall 检验用于检查是否存在单调上升或下降的趋势。由于 M-K 检验是非参数检验，就像 Kendall’s Tau 一样，不需要假设数据的正态性或线性。检验通过分析时间序列中连续点之间的符号来进行。检验的核心思想是，在存在趋势的情况下，如果将符号值求和，它会不断增加或减少。

尽管是非参数检验，但原始检验中有几个假设：

• 时间序列中没有自相关

• 时间序列中没有季节性

多年来，已经对原始检验进行了许多修改，以解决这些问题，很多这种修改，包括原始检验，已在 github.com/mmhs013/pyMannKendall 上实现。它们也作为 pymannkendall 在 pypi 上提供。

预白化 是一种常用的技术，用于去除时间序列中的自相关。简而言之，基本思路如下：

1. 使用 AR(1) 模型识别

2. 

M. Bayazit 和 B. Önöz（2007）建议，如果样本量大于 50 且趋势足够强（slope>0.01），在进行 M-K 检验时不需要使用预白化。对于季节性数据，pymannkendall 也实现了 M-K 检验的季节性变体。

参考检查：

M. Bayazit 和 B. Önöz 的研究论文在 参考文献 部分被引用，引用号为 1。

我们之前讨论过的相同方法check_trend，也实现了 M-K 检验，可以通过设置mann_kendall=True来启用。然而，我们需要记住的一点是，M-K 检验比 Kendall’s Tau 要慢得多，特别是对于长时间序列。M-K 检验还有一些特定的参数：

• seasonal_period：默认值为None。但如果存在季节性，我们可以在此提供seasonal_period，并且会检索到 M-K 检验的季节性变体。

• prewhiten：这是一个布尔标志，用于在应用 M-K 检验之前对时间序列进行预白化处理。默认值为None。在这种情况下，我们会根据之前讨论的条件（N>50）来决定是否进行预白化处理。如果我们在此明确传递True或False，将会尊重此设置。

让我们看看如何使用这个检验：

check_trend(y, confidence=0.05, mann_kendall=True) 

结果将作为一个namedtuple返回，包含以下参数：

• trend：这是一个布尔标志，表示是否存在趋势。

• direction：这将是increasing（增加）或decreasing（减少）。

• slope：这是估算趋势线的斜率。对于 M-K 检验，它将是使用 Theil-Sen 估算器估算的斜率。

• p：这是统计检验的 p 值。

• deterministic：这是一个布尔标志，表示确定性趋势。

让我们来看一个示例，在其中我们对一个时间序列应用了这两个检验（完整代码请参见02-Dealing_with_Non-Stationarity.ipynb）：

# y_unit_root is the a synthetic unit root timeseries
y_unit_root.plot()
plt.show() 

图 7.4：M-K 检验

kendall_tau_res = check_trend(y_unit_root, confidence=0.05)
mann_kendall_res = check_trend(y_unit_root, confidence=0.05, mann_kendall=True)
print(f"Kendalls Tau: Trend: {kendall_tau_res.trend} | Direction: {kendall_tau_res.direction} | Deterministic: {kendall_tau_res.deterministic}")
print(f"Mann-Kendalls: Trend: {mann_kendall_res.trend} | Direction: {mann_kendall_res.direction} | Deterministic: {mann_kendall_res.deterministic}")
## Output 

>> Kendalls Tau: Trend: True | Direction: decreasing | Deterministic: False
>> Mann-Kendalls: Trend: True | Direction: decreasing | Deterministic: False 

如果你能生成一些时间序列，甚至选择一些你遇到的时间序列，使用这些函数来看它是如何工作的，以及结果如何帮助你，那将对你有很大帮助。相关的笔记本中有一些示例可以帮助你入门。你可以观察不同类型趋势下方向和斜率的不同。

现在我们知道如何检测趋势，让我们来看一下去趋势化。

去趋势转换

如果趋势是确定性的，去除趋势将为建模过程增添一些价值。在第三章，分析与可视化时间序列数据中，我们讨论了去趋势化，因为它是我们进行分解时的一个重要部分。但像移动平均或 LOESS 回归这类技术有一个缺点——它们无法进行外推。然而，如果我们考虑的是确定性的线性（甚至是多项式）趋势，那么可以通过线性回归轻松估算。这里的一个额外优点是，识别出的趋势可以轻松进行外推。

这个过程很简单：我们将时间序列对时间的序数表示进行回归，并提取参数。一旦我们获得这些参数，就可以利用日期将趋势外推到未来的任何一点。Python 中的核心逻辑如下：

# y is the time series we are detrending
x = np.arange(len(y))
# degree is the degree of trend we are estimating. Linear, or polynomial
# Fitting a regression on y using a linearly increasing x
linear_params = np.polyfit(x=x, y=y, deg=degree)
# Extract trend using fitted parameters
trend = get_trend(y)
# Now this extracted trend can be removed from y
detrended = y - trend 

我们已经将这个去趋势器作为一个变换器实现，并且放在src.transforms.target_transformations.py中，名为DetrendingTransformer。你可以在 GitHub 仓库中查看我们是如何实现的。现在，让我们看看如何使用它：

from src.transforms.target_transformations import DetrendingTransformer
detrending_transformer = DetrendingTransformer(degree=1)
y_detrended = detrending_transformer.fit_transform(y, freq="1D") 

y_detrended将包含去趋势后的时间序列。要恢复原始时间序列，我们可以使用detrending_transformer调用inverse_transform。相关的笔记本中有示例和图表，展示去趋势如何改变时间序列。

最佳实践：

我们必须小心趋势假设，特别是当我们进行长期预测时。即使是线性趋势假设，也可能导致不现实的预测，因为现实世界中的趋势不会永远以相同的方式延续。通常建议通过某种因子来减弱趋势，，以便在趋势外推时保持保守。这种减弱可以像一样简单。

使时间序列非平稳的另一个关键因素是季节性。让我们来看一下如何识别季节性并将其去除。

检测和修正季节性

绝大多数现实世界的时间序列具有季节性，比如零售销售、能源消耗等。通常，季节性的存在与否是领域知识的一部分。但是，当我们处理时间序列数据集时，领域知识会有所稀释。大多数时间序列可能表现出季节性，但这并不意味着数据集中的每个时间序列都有季节性。例如，在一个零售数据集中，可能有季节性商品，也可能有非季节性商品。因此，当处理时间序列数据集时，能够判断某个特定时间序列是否具有季节性是有价值的。

检测季节性

除了肉眼观察，还有两种常见的检测季节性的方法：自相关和快速傅里叶变换。两者都能自动识别季节性周期。为了讨论的方便，我们将介绍自相关方法，并探讨如何使用该方法来确定季节性。

自相关，如第三章《时间序列数据的分析与可视化》中所解释的，是时间序列与其滞后值之间的相关性。通常，我们期望在立即的滞后期（lag 1，lag 2等）中相关性较高，随着时间推移，相关性逐渐减弱。但对于具有季节性的时间序列，我们也会看到在季节性周期中出现相关性峰值。

让我们通过一个例子来理解这一点。考虑一个合成时间序列，它只是白噪声与一个具有 25 周期季节性信号的正弦波组合（与我们之前在图 7.1中看到的季节性时间序列相同）：

#WhiteNoise + Seasonal
y_random = pd.Series(np.random.randn(length), index=index)
t = np.arange(len(y_random))
y_seasonal = (y_random+1.9*np.cos((2*np.pi*t)/(length/4))) 

如果我们绘制这个时间序列的自相关函数（ACF），它将如下所示（计算和绘制此图的代码可以在02-Dealing_with_Non-Stationarity.ipynb笔记本中找到）：

图 7.5：具有 25 周期季节性信号的合成时间序列自相关图

我们可以看到，除了前几个滞后期外，自相关在接近季节周期时增加，并在季节性达到峰值时达到最高。我们可以利用自相关函数（ACF）这一特性来检测季节性。darts 是一个时间序列预测库，它实现了这种检测季节性的技术。但由于它是为 darts 的时间序列数据结构设计的，因此我们已将相同的逻辑改编为适用于常规的 pandas 序列，代码位于 src.transforms.stationary_utils.py，命名为 check_seasonality。该实现可以执行两种季节性检查。它可以接受一个 seasonality_period 作为输入，并验证数据中是否存在与该 seasonality_period 对应的季节性。如果我们没有提前提供 seasonality_period，它将返回一个统计学上显著的最短 seasonality_period。

这个过程从高层次上来说，执行以下操作：

1. 它计算 ACF。

2. 它会在 ACF 中找到所有的相对最大值。相对最大值是指函数从增加变为减少的转折点。

3. 它检查提供的 seasonal_period 是否为相对最大值。如果不是，我们就得出结论，表示该 seasonal_period 没有相关的季节性。

4. 现在，我们假设 ACF 服从正态分布，并计算指定置信度下的上限。上限由以下公式给出：

其中 r[h] 是滞后 h 时的估计自相关，SE 是标准误差，而 是基于所需置信度的正态分布分位数，。SE 使用 Bartlett 公式近似（有关数学原理，请参考 进一步阅读 部分）。

1. 每一个 seasonality_period 的候选值都会与这个上限进行比较，超过该上限的值被认为是统计上显著的。

除了时间序列本身，这个函数只有三个参数：

• max_lag：指定应该包含在 ACF 中的最大滞后期，并在随后的季节性搜索中使用。这个值应该至少比预期的季节性周期多一个。

• seasonal_period：这是我们从领域知识中给出的季节性周期的直觉，函数会验证这一假设。

• confidence：这是标准的统计置信度水平，默认值为 0.05。

让我们看看如何在我们之前在图 7.4中看到的相同数据上使用这个函数（季节周期为 25）。这将返回一个包含seasonal（表示季节性的布尔标志）和seasonal_periods（具有显著季节性的季节周期）的namedtuple：

# Running the function without specifying seasonal period to identify the seasonality
seasonality_res = check_seasonality(y_seasonal, max_lag=60, confidence=0.05)
print(f"Seasonality identified for: {seasonality_res.seasonal_periods}")
## Output 

>> Seasonality identified for: 25 

This function can also be used to verify if your assumption about the seasonality is right.
# Running the function specifying seasonal period to verify
seasonality_res = check_seasonality(y_seasonal, max_lag=30, seasonal_period=25, confidence=0.05)print(f"Seasonality Test for 25th lag: {seasonality_res.seasonal}")
## Output 

>> Seasonality Test for 25th lag: True 

既然我们已经知道如何识别和检测季节性，接下来我们来谈谈去季节化。

去季节化转换

在第三章，时间序列数据的分析与可视化中，我们回顾了季节性分解的技术。我们可以在这里使用相同的技术，但只需稍作调整。之前，我们并不关心将季节性投射到未来。但在进行预测时，去季节化处理是至关重要的，必须能够将季节性投射到未来。幸运的是，季节性周期的前向投射是非常简单的。这是因为我们正在处理一个固定的季节性配置，它将在季节性周期中不断重复。例如，如果我们为一年的 12 个月份（按月频率的数据）识别了季节性配置，那么为这 12 个月提取的季节性将会在每 12 个月的周期中重复出现。

使用这一特性，我们在src.transforms.target_transformations.py中实现了一个转换器，命名为DeseasonalizingTransformer。我们需要注意一些参数和属性：

• seasonality_extraction：此转换器支持两种提取季节性的方法——"period_averages"，即通过季节性平均来估计季节性配置，以及"fourier_terms"，即通过对傅里叶项进行回归来提取季节性。

• seasonality_period：根据我们用于提取季节性的方法，这个参数可以是整数或字符串。如果是"period_averages"，此参数表示季节性周期重复的周期数。如果是"fourier_terms"，则表示从日期时间索引中提取的季节性。可以使用pandas datetime的属性，如week_of_day、month等来指定最显著的季节性。类似于我们之前看到的FourierDecomposition，我们也可以省略此参数，并在fit/transform方法中提供自定义季节性。

• n_fourier_terms：此参数指定在回归中要包含的傅里叶项的数量。增加此参数会使拟合的季节性变得更加复杂。

• 该实现中没有去趋势处理，因为我们之前已经看到过DetrendingTransformer。该实现假设在使用fit函数之前，任何趋势都已经被去除。

让我们来看一下如何使用它：

from src.transforms.target_transformations import DeseasonalizingTransformer
deseasonalizing_transformer = DeseasonalizingTransformer(seasonality_extraction="period_averages",seasonal_period=25)
y_deseasonalized = deseasonalizing_transformer.fit_transform(y, freq="1D") 

y_deseasonalized将包含去季节化的时间序列。为了恢复到原始时间序列，我们可以使用inverse_transform函数。通常，这可以在做出预测后用于将季节性加回。

最佳实践：

季节性建模可以单独进行，如这里所讨论的，或者使用我们在本章前面讨论的季节性特征。尽管最终的评估需要通过实证方法找到哪个方法更有效，但我们可以遵循一些经验法则/指导方针来决定优先级。

当我们有足够的数据时，让模型将季节性作为主要预测问题的一部分来学习似乎效果更好。但在数据不丰富的情况下，单独提取季节性然后再输入到机器学习模型中也能取得良好的效果。

当数据集具有不同的季节性（即不同时间序列有不同的季节性周期）时，应该根据情况进行处理。可以分别去季节化每个时间序列，或将全局机器学习模型拆分为多个本地模型，每个本地模型都有其独特的季节性模式。

我们之前讨论的最后一个方面是异方差性。让我们也快速看一下这个内容。

检测并修正异方差性

尽管这个名字听起来有些吓人，异方差性其实是一个足够简单的概念。它源自古希腊，hetero 意为 不同，skedasis 意为 离散。正如它的名字所示，当一个变量的变异性在另一个变量中不同时，我们就定义为异方差性。在时间序列的背景下，当时间序列的变异性或离散度随着时间变化时，我们就说时间序列是异方差的。例如，让我们想象一个家庭在若干年的支出情况。在这些年里，这个家庭经历了从贫困到中产阶级，再到上层中产阶级的变化。当家庭贫困时，支出较少，仅限于必需品，因此支出的变异性也较小。但当他们接近上层中产阶级时，家庭有能力购买奢侈品，这在时间序列中造成了波动，从而导致了更高的变异性。如果我们回顾图 7.1，可以看到异方差时间序列的表现。

但是，除了视觉检查外，如果我们能够进行自动化的统计测试来验证异方差性，那就更好了。

检测异方差性

检测异方差性有很多方法，但我们将使用其中最流行的技术之一，即 1980 年由 Halbert White 提出的White 检验。White 检验通过辅助回归任务来检查方差是否恒定。我们先使用一些协变量进行初步回归，并计算该回归的残差。然后，我们用这些残差作为目标，协变量、协变量的平方以及协变量的交叉乘积作为特征，再进行一次回归。最终的统计量是通过这个辅助回归的R²值来估计的。要了解更详细的测试过程，请参考进一步阅读部分；对于严格的数学过程，相关研究论文已在参考文献部分中引用。

参考检查：

要了解 White 检验的严格数学过程，请查看参考文献部分中标注为2的研究论文。

在时间序列的背景下，我们通过使用一个确定性趋势模型来调整这种公式。初步回归是通过将时间作为序数变量进行的，残差用于进行 White 检验。White 检验在statsmodels的het_white中有实现，我们将使用它来执行这个检验。het_white检验返回两个统计量和 p 值——拉格朗日乘数（Lagrangian Multiplier）和 F 统计量。拉格朗日乘数检验残差的方差与回归模型中的自变量之间是否存在关系。F 统计量比较原始模型与允许误差方差变化的模型的拟合优度。任何一个检验的 p 值小于置信度都表示存在异方差性。为了更保守起见，我们也可以使用两个检验，只有当两个 p 值都小于置信度时，才标记为异方差性。

我们已经将这一切封装到src.transforms.stationary_utils中的一个有用函数check_heteroscedasticity里，该函数只有一个附加参数——confidence。让我们来看一下该方法的核心实现（Python 代码）：

import statsmodels.api as sm
# Fitting a linear trend regression
x = np.arange(len(y))
x = sm.add_constant(x)
model = sm.OLS(y,x)
results = model.fit()
# Using the het_white test on residuals
lm_stat, lm_p_value, f_stat, f_p_value = het_white(results.resid, x)
# Checking if both p values are less than confidence
if lm_p_value<confidence and f_p_value < confidence:
        hetero = True
    else:
        hetero = False 

现在，让我们看看如何使用这个函数：

from src.transforms.stationary_utils import check_heteroscedastisticity
check_heteroscedastisticity(y, confidence=0.05) 

它返回一个namedtuple，包含以下参数：

• Heteroscedastic: 一个布尔值标志，指示是否存在异方差性

• lm_statistic: 拉格朗日乘数（LM）统计量

• lm_p_value: 与 LM 统计量相关的 p 值

最佳实践：

我们正在进行的异方差性检验仅考虑回归中的趋势，因此在季节性存在的情况下，可能效果不佳。建议在应用该函数之前先去季节性化数据。

检测异方差性是较为简单的部分。有一些变换方法试图去除异方差性，但每种方法都有其优缺点。我们来看看几种这样的变换。

对数变换

如其名称所示，对数变换是指对时间序列应用对数变换。对数变换有两个主要特性——方差稳定性和减少偏态——从而使数据分布更接近正态分布。在这两点中，我们更关注第一个特性，因为它有助于对抗异方差性。

对数变换通常被认为能减少数据的方差，从而消除数据中的异方差性。直观地看，我们可以将对数变换视为将直方图右侧的极端值拉回，同时拉伸直方图左侧的极低值。

但是已经证明，对数变换并不总是能稳定方差。除此之外，对数变换在机器学习中还带来了另一个挑战。现在，损失的优化发生在对数尺度上。由于对数变换在值范围的较低端压缩得比高端更多，因此学到的模型可能对较低范围的错误不如对较高范围的错误敏感。另一个主要缺点是，对数变换只能应用于严格正的数据。如果你的数据中有零或负值，你将需要通过加上某个常数 M 来偏移整个分布，然后应用变换。这也会在数据中引入一些扰动，可能会产生不利影响。

关键是我们在应用对数变换时需要小心。我们在 src.transforms.target_transformations.py 中实现了一个名为 LogTransformer 的变换器，它只有一个参数 add_one，该参数在变换之前加一，在逆变换后减一。Python 中的关键逻辑就像在变换中应用 np.log1p 或 np.log 函数，分别用 np.expm1 或 np.exp 进行逆变换：

# Transform
np.log1p(y) if self.add_one else np.log(y)
# Inverse Transform
np.expm1(y) if self.add_one else np.exp(y)y_log = log_transformer.fit_transform(y) 

我们所做的就是将它封装成一个漂亮且易于使用的变换器。让我们看看如何使用它：

from src.transforms.target_transformations import LogTransformer
log_transformer = LogTransformer(add_one=True)
y_log = log_transformer.fit_transform(y) 

y_log 是对数变换后的时间序列。我们可以调用 inverse_transform 来恢复原始时间序列。

Box-Cox 变换

尽管对数变换有效且常见，但它是非常 强 的。但对数变换并不是我们能使用的唯一单调变换。还有许多其他的变换，如 y²，， 等等，它们统称为幂变换家族。这个家族中非常著名且广泛使用的一类变换就是 Box-Cox 变换：

和，

参考检查：

Box 和 Cox 的原始研究论文在 参考文献 部分被引用，作为参考文献 3。

直观地，我们可以看到 Box-Cox 变换是一个广义的对数变换。对数变换只是 Box-Cox 变换的特例（当 时）。在不同的 值下，它近似其他变换，如当 时为 y²，当 时为 ，当 时为 ，依此类推。当 时，没有重大变换。

我们之前提到的对数变换的许多缺点在这里也适用，但这些效果的程度有所不同，我们有一个参数，，可以帮助我们决定这些效果的合适程度。像对数变换一样，Box-Cox 变换也仅适用于严格正值的数据。同样，必须对数据分布加上一个常数来进行偏移。该参数的另一面是，它增加了一个需要调节的超参数。

有一些自动化方法可以找到任何数据分布的最佳。其中之一是通过最小化数据分布的对数似然，假设数据服从正态分布。因此，基本上，我们将做的事情是找到最佳的，使数据分布尽可能接近正态分布。这种优化已经在流行的实现中实现，例如scipy中的scipy.special模块中的boxcox函数。

另一种找到最佳的方法是使用 Guerrero 方法，这种方法通常适用于时间序列。在此方法中，我们不是试图将数据分布调整为正态分布，而是尝试最小化时间序列中不同子序列之间的变异性，这些子序列在时间序列中是同质的。子序列的定义有些主观，但通常我们可以安全地认为子序列是季节性长度。因此，我们将要做的是最小化时间序列在不同季节性周期之间的变异性。

参考文献检查：

提出 Guerrero 方法的研究论文已在参考文献部分的参考文献4中引用。

这两种优化方法的工作方式有显著差异，我们在使用它们时需要小心。如果我们的主要关注点是去除时间序列的异方差行为，Guerrero 方法是我们可以使用的方法。

我们在src.transforms.target_transformations.py中提供了一个名为BoxCoxTransformer的转换器。我们需要注意一些参数和属性：

• box_cox_lambda：这是 Box-Cox 变换中使用的参数。如果设置为None，实现将自动找到最佳的。

• optimization：可以是guerrero（默认设置）或loglikelihood。这决定了参数的估计方式。

• seasonal_period：这是使用 Guerrero 方法寻找最佳参数的输入。严格来说，这是子序列的长度，通常取为季节性周期。

• bounds：这是另一个参数，用于通过 Guerrero 方法控制优化。这是一个包含下界和上界的元组，用于搜索最佳的参数。

• add_one：这是一个标志，在应用对数变换之前将一加到序列中，以避免对数为零。

Transformer 中实现的核心逻辑如下：

## Fit Process
# Add one if needed
y = self._add_one(y)
# Find optimum box cox lamda if optimization is Guerrero
self.boxcox_lambda = self._optimize_lambda(y)
## Transform Process
boxcox(y.values, lmbda=self.boxcox_lambda)
## Inverse Transform
self._subtract_one(inv_boxcox(y.values, self.boxcox_lambda)) 

现在，让我们看看如何使用它：

from src.transforms.target_transformations import BoxCoxTransformer
boxcox_transformer = BoxCoxTransformer()
y_boxcox = boxcox _transformer.fit_transform(y) 

y_boxcox 将包含经过 Box-Cox 变换的时间序列。若要恢复到原始时间序列，可以使用 inverse_transform 函数。

Box-Cox 变换和对数变换都可以用于修正异方差性。但是，如前所述，对数变换是一种强烈的变换，而 Box-Cox 变换则为我们提供了另一种手段，允许我们调整和优化变换，以适应我们的数据。我们可以将 Box-Cox 看作是一种灵活的对数变换，可以根据需要调整，以实现适合我们数据的正确变换。请查看笔记本，在其中你可以查看并尝试这些不同的变换，感受它们对数据的影响。

当我们面对大规模的预测问题时，我们将需要分析成百上千甚至数百万个时间序列，才能进行预测。在这种情况下，AutoML 方法显得尤为重要，才能保持实用性。

AutoML 方法用于目标变换

到目前为止，我们已经讨论了许多使序列 更 稳定的方法（这里我们使用稳定的非数学意义），例如去趋势、去季节性、差分和单调变换。我们还查看了统计测试，以检查时间序列中是否存在趋势、季节性等。因此，下一步自然是将这些方法整合起来，以自动化的方式执行这些变换，并在可能的情况下选择合适的默认值。这正是我们所做的，并在 src.transforms.target_transformations 中实现了一个 AutoStationaryTransformer。

以下流程图以自动化方式解释了这一逻辑：

图 7.6：AutoStationaryTransformer 的流程图

我们在这个实现中排除了差分，原因有两个：

• 在预测的背景下，差分带来了相当大的技术债务。如果进行差分，你本质上会使得进行多步预测变得更加困难。虽然可以做到，但它更加困难且灵活性较差。

• 差分可以看作是我们所做的事情的另一种方式。因为差分去除了线性趋势，而季节性差分也去除了季节性。因此，对于自回归时间序列来说，差分可以做很多事情，值得作为独立的变换来使用。

现在，让我们看看可以使用哪些参数来调整 AutoStationaryTransformer：

• confidence：这是统计测试的置信水平，默认值为 0.05。

• seasonal_period：这是季节性周期重复的周期数。如果设置为 None，则 seasonal_period 将从数据中推断出来，默认值为 None。

• seasonality_max_lags: 仅在未提供seasonality_period时使用。这设置了我们搜索季节性的最大滞后。默认为None。

• trend_check_params: 这些是用于趋势统计检验的参数。check_trend默认为{"mann_kendall": False}。

• detrender_params: 这些是传递给DetrendingTransformer的参数。默认为{"degree":1}。

• deseasonalizer_params: 传递给DeseasonalizingTransformer的参数。seasonality_extraction被固定为period_averages。

• box_cox_params: 这些是传递给BoxCoxTransformer的参数。它们默认为{"optimization": "guerrero"}。

让我们将这个AutoStationaryTransformer应用于一个合成时间序列，看看它的效果如何（完整代码在相关笔记本中）：

from src.transforms.target_transformations import AutoStationaryTransformer
auto_stationary = AutoStationaryTransformer(seasonal_period=25)
y_stat = auto_stationary.fit_transform(y_final) 

图 7.7：AutoStationaryTransformer—之前和之后

我们可以看到AutoStationaryTransformer已经对时间序列进行了去季节化和去趋势化处理。在这个特定的例子中，AutoStationaryTransformer应用了去趋势化、去季节化和 Box-Cox 转换。

现在，让我们将这种自动转换应用于我们一直在处理的数据集：

train_df = pd.read_parquet(preprocessed/"selected_blocks_train_missing_imputed_feature_engg.parquet")
transformer_pipelines = {}
for _id in tqdm(train_df["LCLid"].unique()):
    #Initialize the AutoStationaryTransformer with a seasonality period of 48*7
    auto_stationary = AutoStationaryTransformer(seasonal_period=48*7)
    #Creating the timeseries with datetime index
    y = train_df.loc[train_df["LCLid"]==_id, ["energy_consumption","timestamp"]].set_index("timestamp")
    #Fitting and transforming the train
    y_stat = auto_stationary.fit_transform(y, freq="30min")
    # Setting the transformed series back to the dataframe
    train_df.loc[train_df["LCLid"]==_id, "energy_consumption"] = y_stat.values
    #Saving the pipeline
    transformer_pipelines[_id] = auto_stationary 

执行此操作的代码分为两个笔记本，分别为02-Dealing_with_Non-Stationarity.ipynb和02a-Dealing_with_Non-Stationarity-Train+Val.ipynb，位于Chapter06文件夹中。前者对训练数据进行自动稳态转换，而后者对训练和验证数据合并进行转换。这是为了模拟我们如何对验证数据进行预测（仅使用训练数据进行训练），以及对测试数据进行预测（在训练中使用训练和验证数据）。

这个过程稍微耗时。我建议您运行笔记本，吃午餐或小吃，然后回来。一旦完成，02-Dealing_with_Non-Stationarity.ipynb笔记本将保存一些文件：

• selected_blocks_train_auto_stat_target.parquet: 一个 DataFrame，其索引为LCLid和timestamp，而转换后的目标

• auto_transformer_pipelines_train.pkl: 一个 Python 字典，包含每个LCLid的AutoStationaryTransformer，以便我们将来可以反转转换

02a-Dealing_with_Non-Stationarity-Train+Val.ipynb笔记本还保存了用于训练和验证数据集的相应文件。

我们正在处理的数据集几乎没有趋势，并且在整个过程中非常稳定。这些转换的影响将在具有强趋势和异方差性的时间序列中更为明显。

最佳实践：

这种在建模之前进行的显式去趋势和去季节性化也可以视为一种增强方法。这应该被视为另一种将所有这些因素一起建模的替代方法。在某些情况下，让模型通过数据驱动的方式从头到尾学习，可能比通过显式去趋势和去季节性化注入这些强烈的归纳偏见效果更好，反之亦然。交叉验证的测试得分应该始终拥有最终的发言权。

恭喜你顺利完成了这一章，里面充满了新概念、一些统计学内容和数学内容。从应用机器学习模型于时间序列的角度来看，本章的概念将对提升你的模型水平非常有帮助。

小结

在上一章深入实践后，我们停留在那儿，继续回顾了像平稳性这样的概念，以及如何处理这些非平稳的时间序列。我们学习了可以显式处理非平稳时间序列的技术，例如差分、去趋势、去季节性等。为了将这些内容结合起来，我们看到了自动转换目标的方式，学习了如何使用提供的实现，并将其应用于我们的数据集。现在我们掌握了将时间序列有效转换为机器学习数据集所需的技能，在下一章中，我们将开始使用我们创建的特征对数据集应用一些机器学习模型。

参考文献

以下是本章的参考文献：

1. Bayazit, M. 和 Önöz, B. (2007)，在趋势分析中是否需要预处理？，《水文学科学期刊》，52:4，611–624。 doi.org/10.1623/hysj.52.4.611.

2. White, H. (1980)，异方差一致的协方差矩阵估计量及异方差性的直接检验。《计量经济学》 第 48 卷，第 4 期（1980 年 5 月），817–838 页（22 页）。 doi.org/10.2307/1912934.

3. Box, G. E. P. 和 Cox, D. R. (1964)，变换分析。《皇家统计学会学报》B 系列，26，211–252。 www.ime.usp.br/~abe/lista/pdfQWaCMboK68.pdf.

4. Guerrero, Victor M. (1993), 通过幂变换支持的时间序列分析。《预测学期刊》，第 12 卷，第 1 期，37–48。 onlinelibrary.wiley.com/doi/10.1002/for.3980120104.

进一步阅读

要了解本章所涉及的更多主题，请参考以下资源：

• 时间序列分析中的平稳性，作者：Shay Palachy：towardsdatascience.com/stationarity-in-time-series-analysis-90c94f27322

• 比较 R 中的 ADF 检验函数，作者 Fabian Kostadinov（相同的概念也可以在 Python 中实现）：fabian-kostadinov.github.io/2015/01/27/comparing-adf-test-functions-in-r/

• 肯达尔的 tau 值：www.statisticshowto.com/kendalls-tau/

• Mann-Kendall 趋势检验：www.statisticshowto.com/wp-content/uploads/2016/08/Mann-Kendall-Analysis-1.pdf

• Theil-Sen 估计器：en.wikipedia.org/wiki/Theil%E2%80%93Sen_estimator

• 使用自相关图进行统计推断—维基百科：en.wikipedia.org/wiki/Correlogram#Statistical_inference_with_correlograms

• 白检验用于异方差性检测：itfeature.com/hetero/white-test-of-heteroscedasticity/

加入我们的 Discord 社区

加入我们社区的 Discord 空间，与作者和其他读者进行讨论：

packt.link/mts

第八章：使用机器学习模型进行时间序列预测

在上一章中，我们开始将机器学习作为解决时间序列预测问题的工具。我们讨论了几种技术，如时间延迟嵌入和时间嵌入，这两者将时间序列预测问题作为机器学习范式中的经典回归问题。在本章中，我们将详细探讨这些技术，并通过本书中一直使用的伦敦智能电表数据集，实践这些技术。

本章将涵盖以下主题：

• 使用机器学习模型进行训练和预测

• 生成单步预测基准

• 标准化代码用于训练和评估机器学习模型

• 为多个家庭进行训练和预测

技术要求

您需要按照本书前言中的说明设置Anaconda环境，以便为本书中的代码提供一个包含所有必需库和数据集的工作环境。在运行笔记本时，任何额外的库都会被安装。

在使用本章代码之前，您需要运行以下笔记本：

• 02-Preprocessing_London_Smart_Meter_Dataset.ipynb位于Chapter02中

• 01-Setting_up_Experiment_Harness.ipynb位于Chapter04中

• 01-Feature_Engineering.ipynb位于Chapter06中

• 02-Dealing_with_Non-Stationarity.ipynb位于Chapter07中

• 02a-Dealing_with_Non-Stationarity-Train+Val.ipynb位于Chapter07中

本章的代码可以在github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-/tree/main/notebooks/Chapter08找到。

使用机器学习模型进行训练和预测

在第五章，时间序列预测作为回归中，我们讨论了监督式机器学习的示意图（图 5.2）。在示意图中，我们提到监督学习问题的目的是找到一个函数，，其中 是预测值，X 是作为输入的特征集， 是模型参数，h 是理想函数的近似。在本节中，我们将更详细地讨论 h，并看看如何使用不同的机器学习模型来估计它。

h 是任何近似理想函数的函数，但可以视为一个函数族中所有可能函数的元素。更正式地，我们可以这样说：

这里，H 是我们也称之为模型的一个函数族。例如，线性回归是一种模型或函数族。对于每个系数值，线性回归模型都会给出一个不同的函数，H 就是线性回归模型可以生成的所有可能函数的集合。

有许多种函数或模型可供选择。为了更全面地理解这一领域，我们需要参考其他机器学习资源。进一步阅读部分包含了一些可能帮助您开始学习的资源。至于本书的范围，我们将其狭义地定义为应用机器学习模型进行预测，而不是机器学习的全部内容。尽管我们可以使用任何回归模型，但我们将仅回顾几个流行且有用的时间序列预测模型，并观察它们的实际应用。我们鼓励您自己探索其他算法，熟悉它们。然而，在查看不同模型之前，我们需要再次生成几个基准。

生成单步预测基准

我们在第四章，设置强基准预测中回顾并生成了几个基准模型。但有一个小问题——预测范围。在第六章，时间序列预测的特征工程中，我们讨论了机器学习模型一次只能预测一个目标，而且我们坚持使用单步预测。我们之前生成的基准并非单步预测，而是多步预测。生成单步预测的基准算法，如 ARIMA 或 ETS，需要我们根据历史数据进行拟合，预测一步，然后再次拟合，增加一天数据。以这种迭代方式对测试或验证期进行预测需要我们迭代约 1,440 次（每天 48 个数据点，30 天），并且需要对我们选择的数据集中的所有家庭（在我们的例子中是 150 个）进行此操作。这将需要相当长的时间来计算。

我们选择了简单方法和季节性简单方法（第四章，设置强基准预测），这两种方法可以作为原生的 pandas 方法实现，用作生成单步预测的基准方法。

简单预测在单步预测中表现异常好，可以被认为是一个强有力的基准。在Chapter08文件夹中，有一个名为00-Single_Step_Backtesting_Baselines.ipynb的笔记本，它生成这些基准并将其保存到磁盘。让我们现在运行这个笔记本。该笔记本会为验证集和测试集生成基准，并将预测、指标和聚合指标保存到磁盘。测试期的聚合指标如下：

图 8.1：单步基准的聚合指标

为了简化这些模型的训练和评估，我们在整个过程中使用了统一的结构。让我们快速回顾一下这个结构，以便您能够更好地跟随笔记本的内容。

标准化代码以训练和评估机器学习模型

训练机器学习模型时有两个主要因素 —— 数据 和 模型 本身。因此，为了标准化流程，我们定义了三个配置类（FeatureConfig、MissingValueConfig 和 ModelConfig），以及一个封装类（MLForecast），用于 scikit-learn 风格的估算器（.fit - .predict），以使过程更加顺畅。让我们逐一了解它们。

笔记本提示：

要跟随代码进行操作，请使用 Chapter08 文件夹中的 01-Forecasting_with_ML.ipynb 笔记本和 src 文件夹中的代码。

FeatureConfig

FeatureConfig 是一个 Python dataclass，定义了一些在处理数据时必需的关键属性和函数。例如，连续型、分类型和布尔型列需要分别进行预处理，才能输入机器学习模型。让我们看看 FeatureConfig 包含了什么：

• date：一个必填列，设置 DataFrame 中 date 列的名称。

• target：一个必填列，设置 DataFrame 中 target 列的名称。

• original_target：如果 target 包含已转换的目标（如对数转换、差分等），则 original_target 指定没有经过转换的目标列的名称。这对于计算依赖于训练历史的指标，如 MASE，非常重要。如果未提供此参数，则默认认为 target 和 original_target 是相同的。

• continuous_features：一个连续特征列表。

• categorical_features：一个分类特征列表。

• boolean_features：一个布尔特征列表。布尔特征是分类特征，但只有两个独特值。

• index_cols：在预处理过程中设置为 DataFrame 索引的列列表。通常，我们会将日期时间以及在某些情况下时间序列的唯一 ID 作为索引。

• exogenous_features：一个外生特征列表。DataFrame 中的特征可能来自特征工程过程，如滞后或滚动特征，也可能来自外部数据源，例如我们数据集中的温度数据。这是一个可选字段，允许我们将外生特征与其余特征区分开来。此列表中的项应是 continuous_features、categorical_features 或 boolean_features 的子集。

除了对输入进行一些验证外，类中还有一个名为 get_X_y 的有用方法，包含以下参数：

• df：一个包含所有必要列的 DataFrame，包括目标列（如果可用）。

• categorical：一个布尔标志，用于指示是否包含分类特征。

• exogenous：一个布尔标志，用于指示是否包含外生特征。

该函数返回一个元组 (features, target, original_target)。

我们只需像初始化任何其他类一样，初始化该类，并将特征名称分配给类的参数。包含所有特征的完整代码可在随附的笔记本中找到。

在设置FeatureConfig数据类后，我们可以将任何定义了特征的 DataFrame 传递给get_X_y函数，获取特征、目标和原始目标：

train_features, train_target, train_original_target = feat_config.get_X_y(
    sample_train_df, categorical=False, exogenous=False
) 

如你所见，我们在这里没有使用类别特征或外生特征，因为我想专注于核心算法，并展示它们如何成为我们之前看到的其他经典时间序列模型的直接替代。我们将在第十五章，全球深度学习预测模型策略中讨论如何处理类别特征。

MissingValueConfig

另一个关键设置是如何处理缺失值。我们在第三章，分析与可视化时间序列数据中看到了一些填充时间序列上下文中的缺失值的方法，并且我们已经填充了缺失值并准备好了数据集。但是，在将时间序列转换为回归问题所需的特征工程中，还会产生一些缺失值。例如，在创建滞后特征时，数据集中的最早日期将没有足够的数据来创建滞后，因此会留空。

最佳实践：

尽管用零或均值填充是大多数数据科学家社区的默认方法或常用方法，但我们应该始终尽力尽可能智能地填充缺失值。在滞后特征方面，填充零可能会扭曲特征。与其填充零，不如使用向后填充（使用列中的最早值向后填充），这可能会更合适。

一些机器学习模型自然处理空值或NaN特征，而对于其他机器学习模型，我们需要在训练之前处理这些缺失值。如果我们能定义一个config，为一些预期会有NaN信息的列设置如何填充这些缺失值，那将非常有帮助。MissingValueConfig是一个 Python dataclass，正是用来做这件事的。让我们看看它包含了什么：

• bfill_columns：需要使用向后填充策略来填充缺失值的列名列表。

• ffill_columns：需要使用向前填充策略来填充缺失值的列名列表。如果某列名同时出现在bfill_columns和ffill_columns中，则该列首先使用向后填充，剩余的缺失值则使用向前填充策略。

• zero_fill_columns：需要用零填充的列名列表。

填充缺失值的顺序是bfill_columns，然后是ffill_columns，最后是zero_fill_columns。作为默认策略，数据类使用列均值填充缺失值，因此即使你没有为某列定义任何策略，缺失值也会使用列均值填充。有一个叫做impute_missing_values的方法，它接受 DataFrame 并根据指定的策略填充空单元格。

ModelConfig

ModelConfig是一个 Python 的dataclass，它包含了关于建模过程的一些细节，比如是否对数据进行标准化、是否填充缺失值等等。让我们详细看一下它包含的内容：

• model: 这是一个必需的参数，可以是任何 scikit-learn 风格的估算器。

• name: 模型的字符串名称或标识符。如果未使用该参数，它将恢复为作为model传递的类名称。

• normalize: 一个布尔标志，用于设置是否对输入应用StandardScaler。

• fill_missing: 一个布尔标志，用于设置是否在训练之前填充空值。一些模型能够自然处理NaN，而其他模型则不能。

• encode_categorical: 一个布尔标志，用于设置是否在拟合过程中将类别列进行编码。如果为False，则期望单独进行类别编码，并将其作为连续特征的一部分包含。

• categorical_encoder: 如果encode_categorical为True，则categorical_encoder是我们可以使用的 scikit-learn 风格的编码器。

让我们看看如何定义ModelConfig数据类：

model_config = ModelConfig(
    model=LinearRegression(),
    name="Linear Regression",
    normalize=True,
    fill_missing=True,
) 

它只有一个方法，clone，该方法会将估算器和配置克隆到一个新的实例中。

MLForecast

最后但同样重要的是，我们有一个围绕 scikit-learn 风格模型的包装类。它使用我们之前讨论的不同配置来封装训练和预测函数。让我们看看初始化模型时有哪些可用的参数：

• model_config: 我们在ModelConfig部分讨论过的ModelConfig类的实例。

• feature_config: 我们之前讨论过的FeatureConfig类的实例。

• missing_config: 我们之前讨论过的MissingValueConfig类的实例。

• target_transformer: 来自src.transforms的目标转换器实例。它应支持fit、transform和inverse_transform。它还应返回带有日期时间索引的pd.Series，以确保无错误运行。如果我们单独执行了目标转换，那么在预测时也会使用它来执行inverse_transform。

MLForecast有一些功能，可以帮助我们管理模型生命周期，一旦初始化。让我们看一下。

拟合函数

fit函数的目的与 scikit-learn 的fit函数相似，但它额外处理了标准化、类别编码和目标转换，使用了三个配置中的信息。该函数的参数如下：

• X: 这是一个包含要在模型中使用的特征（作为列）的 pandas DataFrame。

• y: 这是目标，可以是 pandas DataFrame、pandas Series 或 numpy 数组。

• is_transformed: 这是一个布尔参数，告诉我们目标是否已经被转换。如果为True，即使我们已经用target_transformer初始化了对象，fit方法也不会转换目标。

• fit_kwargs：这是一个 Python 字典，包含需要传递给估算器fit函数的关键字参数。

predict 函数

predict函数处理推断。它是对 scikit-learn 估算器的predict函数的封装，但与fit一样，它还做了一些其他事情，比如标准化、分类编码和反转目标转换。这个函数只有一个参数：

• X：一个包含特征的 pandas DataFrame，作为模型的列。DataFrame 的索引将传递到预测中。

feature_importance 函数

feature_importance函数从模型中提取特征重要性（如果有的话）。对于线性模型，它提取系数；而对于基于树的模型，它提取内建的特征重要性并返回一个排序后的 DataFrame。

用于评估模型的辅助函数

虽然我们之前看到的其他函数处理了核心训练和预测，我们还希望评估模型、绘制结果等等。我们在笔记本或代码库中也定义了这些函数。下面的函数是用来在笔记本中评估模型的：

def evaluate_model(
    model_config,
    feature_config,
    missing_config,
    train_features,
    train_target,
    test_features,
    test_target,
):
    ml_model = MLForecast(
        model_config=model_config,
        feature_config=feat_config,
        missing_config=missing_value_config,
    )
    ml_model.fit(train_features, train_target)
    y_pred = ml_model.predict(test_features)
    feat_df = ml_model.feature_importance()
    metrics = calculate_metrics(test_target, y_pred, model_config.name, train_target)
    return y_pred, metrics, feat_df 

这为我们提供了评估所有不同模型的标准方法，并且能够在大规模上自动化这个过程。我们还有一个用于计算指标的函数calculate_metrics，它定义在src/forecasting/ml_forecasting.py中。

本书中提供的标准实现并非一刀切的方法，而是最适合本书的流程和数据集的方法。请不要将其视为一个强健的库，而是一个很好的起点和指南，帮助你开发自己的代码。

现在我们有了基准线和应用不同模型的标准方法，让我们回到讨论不同模型的内容。接下来的讨论中，我们暂时不考虑时间因素，因为我们已将时间序列预测问题转化为回归问题，并将时间作为问题的特征（滞后和滚动特征）进行处理。

线性回归

线性回归是一类函数，具有以下形式：

在这里，k是模型中的特征数量，且是模型的参数。每个特征都有一个，还有一个，我们称之为截距，它是从数据中估算出来的。实质上，输出是特征向量X[i]的线性组合。顾名思义，这是一个线性函数。

模型参数可以通过数据 D(X[i,] y[i])来估算，使用优化方法和损失函数，最常用的估算方法是普通最小二乘法（OLS）。在这里，我们找到模型参数，它最小化残差平方和（均方误差（MSE））：

这里的损失函数非常直观。我们本质上是在最小化训练样本与我们预测点之间的距离。平方项作为一种技术，能够避免正负误差相互抵消。除了损失函数的直观性之外，另一个广泛选择它的原因是最小二乘法有解析解，因此我们不需要使用像梯度下降这样计算密集型的优化技术。

线性回归深深植根于统计学，在正确的假设下，它可以成为一个强大的工具。线性回归通常与五个假设相关，如下所示：

• 自变量和因变量之间的关系是线性的。

• 误差服从正态分布。

• 误差的方差在所有自变量的值范围内是恒定的。

• 误差中没有自相关性。

• 自变量之间几乎没有相关性（多重共线性）。

但是，除非你担心使用线性回归来得出预测区间（预测结果可能出现在某个区间内的概率），否则我们可以在一定程度上忽略除第一个假设外的所有假设。

线性假设（第一个假设）很重要，因为如果变量之间不是线性相关的，会导致欠拟合，从而表现不佳。我们可以通过将输入投影到更高维空间中，在一定程度上解决这个问题。从理论上讲，我们可以将一个非线性问题投影到更高维空间，在那里问题变成线性。例如，考虑一个非线性函数，。如果我们在和的输入空间中运行线性回归，我们知道得到的模型将会严重欠拟合。但如果我们通过使用多项式变换将输入空间从和投影到，和，那么y的函数将成为完美的线性拟合。

多重共线性假设（最后一个假设）与线性函数的拟合部分相关，因为当我们有高度相关的自变量时，估计的系数会变得非常不稳定且难以解释。拟合函数仍然能够很好地工作，但由于存在多重共线性，即使是输入中的微小变化也会导致系数的幅度和符号发生变化。如果你正在使用纯线性回归，检查多重共线性是一个最佳实践。这通常是时间序列中的一个问题，因为我们提取的特征，如滞后和滚动特征，可能会相互关联。因此，在使用和解释时间序列数据上的线性回归时，我们需要格外小心。

现在，让我们看看如何使用线性回归并评估来自验证数据集的一个样本家庭的拟合效果：

from sklearn.linear_model import LinearRegression
model_config = ModelConfig(
    model=LinearRegression(),
    name="Linear Regression",
    # LinearRegression is sensitive to normalized data
    normalize=True,
    # LinearRegression cannot handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

单步预测看起来不错，且已经优于天真预测（MAE = 0.173）：

图 8.2：线性回归预测

模型的系数，（可以通过已训练的 scikit-learn 模型的coef_属性访问），展示了每个特征对输出的影响程度。因此，提取并绘制这些系数为我们提供了对模型的第一层可视化。让我们来看看模型的系数：

图 8.3：线性回归的特征重要性（前 15 名）

如果我们查看特征重要性图表中的Y轴，可以看到它是以十亿为单位的，因为某些特征的系数以十亿的数量级存在。我们还可以看到这些特征是基于傅里叶级数的特征，并且彼此相关。尽管我们有很多系数在十亿级别，但我们可以发现它们分布在零的两侧，因此它们在函数中实际上会相互抵消。这就是我们之前讨论的多重共线性问题。我们可以通过去除多重共线特征并进行某种特征选择（前向选择或后向消除）来使线性模型变得更好。

但是，既然如此，不如我们来看看可以对线性模型进行的一些修改，这些修改能使其在面对多重共线性和特征选择时更为稳健。

正则化线性回归

我们在第五章《时间序列预测与回归》中简要讨论了正则化，并提到正则化在一般意义上是我们在学习过程中对学习函数复杂度进行约束的任何手段。线性模型变得更加复杂的一种方式是系数的幅度很大。例如，在线性拟合中，我们的系数为 200 亿。这个特征的任何微小变化都会导致预测结果出现巨大波动。直观地看，如果我们有一个较大的系数，函数就会变得更加灵活和复杂。我们可以通过应用正则化（如权重衰减）来解决这个问题。权重衰减是指在损失函数中加入一个惩罚项，用于惩罚系数的幅度。损失函数，即残差平方和，现在变为如下形式：

这里，W 是权重衰减， 是正则化的强度。

W 通常是权重矩阵的范数。在线性代数中，矩阵的范数是衡量其元素大小的一种方式。矩阵有许多种范数，但用于正则化的两种最常见的范数是 L1 范数和 L2 范数。当我们使用 L1 范数来正则化线性回归时，我们称之为 lasso 回归，而当我们使用 L2 范数时，我们称之为 ridge 回归。当我们应用权重衰减正则化时，我们迫使系数变小，这意味着它也充当了内部特征选择的作用，因为那些没有提供太多价值的特征会得到非常小或零的系数（取决于正则化类型），这意味着它们在最终的函数中贡献甚微甚至没有贡献。

L1 范数定义为矩阵各元素绝对值的和。对于权重衰减正则化，L1 范数可以表示为如下形式：

L2 范数定义为矩阵各个值的平方和。对于权重衰减正则化，L2 范数可以表示为如下形式：

通过将这个项添加到线性回归的损失函数中，我们迫使系数变小，因为优化器在减少 RSS 的同时，还会激励它减少 W。

我们思考正则化的另一个方式是从线性代数和几何学的角度来看。

接下来的部分讨论正则化的几何直觉。虽然这将使你对正则化的理解更加深刻，但它并不是理解本书其余部分的必备条件。所以，如果你时间紧迫或者以后有时间再回来复习，你可以跳过接下来的部分，只阅读 关键点 提示。

正则化——几何视角

如果我们从稍微不同的角度来看 L1 和 L2 范数，我们会发现它们是距离的度量。

让B表示线性回归中所有系数的向量，。向量是一个数字数组，但从几何学的角度来看，它也是从原点到n维坐标空间中某一点的箭头。现在，L2 范数仅仅是该向量 B 所定义的空间中从原点到该点的欧几里得距离。L1 范数是从原点到该点的曼哈顿距离或出租车距离，该点也由向量 B 定义。让我们通过图示来看看这一点：

图 8.4：欧几里得距离与曼哈顿距离

欧几里得距离是从原点到该点的直线距离。但如果我们只能平行于两个坐标轴移动，那么我们首先需要沿一个坐标轴走距离 ，然后再沿另一个坐标轴走距离 。这就是曼哈顿距离。

假设我们在一个城市（例如曼哈顿），那里建筑物排布成方形街区，直路相交成直角，我们想从 A 点到 B 点。欧几里得距离是从 A 点到 B 点的直接距离，在现实中，只有在我们能在建筑物顶上进行跑酷的情况下才能实现。而曼哈顿距离则是出租车沿着直角道路从 A 点到 B 点的实际行驶距离。

为了进一步发展对 L1 和 L2 范数的几何直觉，我们进行一个思想实验。如果我们在二维空间中移动点 ，同时保持欧几里得距离或 L2 范数不变，我们将得到一个以原点为中心的圆。在三维空间中，这将变成一个球体，在 n 维空间中则是一个超球体。如果我们保持 L1 范数不变而追踪相同的路径，我们将得到一个以原点为中心的菱形。在三维空间中，这将变成一个立方体，在 n 维空间中是一个超立方体。

现在，当我们优化权重时，除了减少损失函数的主要目标外，我们还鼓励系数保持在一个定义的距离（范数）内，远离原点。从几何角度来看，这意味着我们要求优化找到一个向量 ，该向量在最小化损失函数的同时，也保持在由范数定义的几何形状（圆形或正方形）内。我们可以在以下图表中看到这一点：

图 8.5：使用 L1 范数（套索回归）与 L2 范数（岭回归）进行正则化

图中的同心圆是损失函数的等高线，最内圈为最低点。随着我们向外移动，损失增加。因此，正则化回归不会选择 ，而是选择一个与范数几何形状相交的 。

这种几何解释也使得理解岭回归和套索回归之间的另一个关键区别变得更容易。由于 L1 范数的作用，套索回归会产生一个稀疏解。之前，我们提到过权重衰减正则化会进行隐式特征选择。但根据你是应用 L1 还是 L2 范数，隐式特征选择的方式是不同的。

关键点：

对于 L2 范数，较不相关特征的系数会被推近零，但不会完全为零。该特征仍然会在最终的函数中发挥作用，但其影响会微乎其微。另一方面，L1 范数则将这些特征的系数完全推向零，从而得到一个稀疏解。因此，L1 正则化促进稀疏性和特征选择，而 L2 正则化通过将系数缩小到接近零来减少模型复杂度，但不一定会完全消除任何系数。

使用正则化的几何解释可以更好地理解这一点。在优化中，通常在极值点或角落处找到有趣的点。圆形没有角落，因此创建了 L2 范数；极小值可以位于圆的任何边缘。但对于菱形，我们有四个角，极小值会出现在这些角落。因此，使用 L2 范数时，解可以接近零，但不一定是零。然而，使用 L1 范数时，解会出现在角落处，在那里系数可以被压缩到绝对零。

现在，让我们看看如何使用岭回归并评估我们验证数据集中一个家庭的拟合情况：

from sklearn.linear_model import RidgeCV
model_config = ModelConfig(
    model=RidgeCV(),
    name="Ridge Regression",
    # RidgeCV is sensitive to normalized data
    normalize=True,
    # RidgeCV does not handle missing values
    fill_missing=True
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

让我们来看看RidgeCV的单步预测。它看起来与线性回归非常相似。即使是 MAE，对于这个家庭来说也完全相同：

图 8.6：岭回归预测

但用 L2 正则化模型来看系数是很有意思的。让我们来看看模型的系数：

图 8.7：岭回归特征重要性（前 15 名）

现在，Y轴看起来合理且较小。多重共线性特征的系数已经缩小到一个更合理的水平。像滞后特征这样的特征，理应具有较大影响力，已经占据了前几位。正如你可能记得的，在线性回归中（图 8.3），这些特征被傅里叶特征的巨大系数所压制。我们这里只绘制了前 15 个特征，但如果查看完整的列表，你会看到很多特征的系数接近零。

现在，让我们尝试对这个家庭样本进行套索回归：

from sklearn.linear_model import LassoCV
model_config = ModelConfig(
    model=LassoCV(),
    name="Lasso Regression",
    # LassoCV is sensitive to normalized data
    normalize=True,
    # LassoCV does not handle missing values
    fill_missing=True
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

让我们来看看LassoCV的单步预测。像岭回归一样，与线性回归几乎没有视觉差异：

图 8.8：套索回归预测

让我们来看看模型的系数：

图 8.9：套索回归特征重要性（前 15 名）

系数与岭回归非常相似，但如果查看完整的系数列表（在笔记本中），你会看到很多特征的系数为零。

即使在 MAE、MSE 等指标相同的情况下，岭回归或套索回归仍然比线性回归更受偏好，因为正则化回归带来了额外的稳定性和鲁棒性，尤其是在预测中，多重共线性几乎总是存在。但我们需要记住，所有的线性回归模型仍然只是捕捉线性关系。如果数据集具有非线性关系，那么线性回归的拟合效果不会太好，有时甚至会非常糟糕。

现在，让我们换个方向，看看另一类模型——决策树。

决策树

决策树是另一类函数，其表达能力远超过线性函数。决策树将特征空间划分为不同的子空间，并对每个子空间拟合一个非常简单的模型（例如平均值）。让我们通过一个例子来理解这种划分是如何工作的。我们考虑一个回归问题，目标是预测Y，并且只有一个特征X，如下图所示：

图 8.10：由决策树划分的特征空间

我们可以立即看到，拟合一个线性函数会导致欠拟合。但决策树的做法是将特征空间（在这里就是X）划分为不同的区域，其中目标Y是相似的，然后对每个区域拟合一个简单的函数，如平均值（因为这是一个回归问题）。在这种情况下，决策树将特征空间划分为 A、B 和 C 三个区域。现在，对于任何落入区域 A 的X，预测函数将返回区域 A 中所有点的平均值。

这些划分是通过使用数据创建决策树来形成的。直观地说，决策树通过创建一组 if-else 条件来划分特征空间，并尝试找到最佳的划分方式，以最大化目标变量在每个划分中的同质性。理解决策树作用的一种有帮助的方法是将数据点想象成珠子，沿着树流动，根据其特征选择路径，最终停留在一个最终位置。在我们讨论如何从数据中创建决策树之前，让我们先看看它的组成部分，并理解相关的术语：

图 8.11：决策树的结构

决策树中有两种类型的节点——决策节点 和 叶节点。决策节点是我们之前提到的 if-else 语句。这个节点将有一个基于数据点流向左支路或右支路的条件。位于最顶部的决策节点有一个特殊的名称——根节点。最后，根据条件划分数据点并将其导向右支路或左支路的过程被称为划分。叶节点是没有其他分支的节点。它们是“珠子流经树”的类比中的最终停靠点。这些就是我们在本节中讨论的划分。

正式地，我们可以定义由具有 M 个划分的决策树生成的函数 P[1]、P[2]、…、P[M]，如下所示：

在这里，x 是输入，c[m] 是该区域的常数响应，P[m]，而 I 是一个函数，当 时为 1；否则，它为 0。

对于回归树，我们通常采用平方损失作为损失函数。在这种情况下，c[m] 通常设置为所有 y 的平均值，其中对应的 x 落入 P[m] 划分中。

现在我们知道了决策树是如何工作的，剩下的就是理解如何决定在什么特征上进行划分，以及如何划分特征。

多年来，已经提出了许多算法来根据数据创建决策树，比如 ID3、C4.5、CART 等。在这些算法中，使用 分类与回归树（CART）是最流行的方法之一，并且它也支持回归。因此，本书中我们将坚持使用 CART。分类树用于目标变量为类别型（例如，预测类标签）时。回归树用于目标变量为连续型（例如，预测数值）时。

全局最优的二分划分集合，能够最小化平方和，通常是无法处理的。因此，我们采用贪心算法来创建决策树。贪心优化是一种启发式方法，它通过逐步构建解决方案，每一步选择局部最优解。因此，我们不会全局寻找最佳的特征划分，而是通过逐个节点创建决策树，在每个阶段选择最优的特征划分。对于回归树，我们选择一个划分特征 f 和划分点 s，使得它创建两个划分 P[1] 和 P[2]，其最小化如下：

在这里，c[1] 和 c[2] 是所有 y 的平均值，其中对应的 x 落在 P[1] 和 P[2] 之间。

因此，通过使用这个标准，我们可以不断地将区域进一步划分。每划分一次，我们就增加一次树的深度。在某个时刻，我们会开始过拟合数据集。但是如果我们没有做足够的划分，可能也会导致欠拟合数据。一个策略是当我们达到预定的深度时停止进一步的划分。在 DecisionTreeRegressor 的 scikit-learn 实现中，这对应于 max_depth 参数。这个超参数需要通过验证数据集进行估算。还有其他策略可以停止划分，比如设置划分所需的最小样本数（min_samples_split），或者进行划分所需的最小成本下降（min_impurity_decrease）。有关 DecisionTreeRegressor 中参数的完整列表，请参阅文档 scikit-learn.org/stable/modules/generated/sklearn.tree.DecisionTreeRegressor.html。

现在，让我们看看如何使用决策树并在验证数据集中的一个家庭样本上评估拟合情况：

from sklearn.tree import DecisionTreeRegressor
model_config = ModelConfig(
    model=DecisionTreeRegressor(max_depth=4, random_state=42),
    name="Decision Tree",
    # Decision Tree is not affected by normalization
    normalize=False,
    # Decision Tree in scikit-learn does not handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

让我们来看一下 DecisionTreeRegressor 的单步预测。它的表现不如我们迄今为止运行的线性回归或正则化线性回归模型：

图 8.12：决策树预测

对于线性模型，一些系数帮助我们理解每个特征在预测函数中的重要性。而在决策树中，我们没有系数，但特征的重要性依然是通过损失函数的平均下降进行估算的，这一变化归因于树构建过程中的每个特征。在 scikit-learn 模型中，我们可以通过训练模型的 feature_importance_ 属性来访问这个特征重要性。让我们来看看这个特征重要性：

图 8.13：决策树的特征重要性（前 15）

最佳实践：

尽管默认的特征重要性是快速且简便的方式来检查不同特征的使用情况，但在将其用于其他目的（如特征选择或业务决策）之前，应谨慎处理。这种评估特征重要性的方法会给一些连续特征和高基数类别特征误导性地赋予过高的值。建议使用置换重要性（sklearn.inspection.permutation_importance）来进行更好的特征重要性评估，这是一种简单但更好的评估方法。进一步阅读部分包含了一些关于模型可解释性的资源，这些资源可以作为理解模型影响因素的良好起点。

在这里，我们可以看到重要的特征，如滞后特征和季节性滚动特征，排在了最上面。

我们在第五章，时间序列预测作为回归中讨论了过拟合和欠拟合。这些也在机器学习术语中被称为高偏差（欠拟合）和高方差（过拟合）（进一步阅读部分包含了链接，如果你想了解更多关于偏差、方差及它们之间的权衡问题，可以参考这些资源）。

决策树是一种非常容易发生过拟合或高方差的算法，因为与线性函数不同，决策树在具有足够表现力的情况下，能够通过划分特征空间来记住训练数据集的特征。另一个关键缺点是决策树无法外推。假设我们有一个特征f，它线性地增加我们的目标变量y。我们拥有的训练数据中，f[max]是f的最大值，y[max]是y的最大值。由于决策树划分了特征空间，并为每个区域分配了一个常数值，即使我们提供了f > f[max]，它依然会仅仅给出一个预测值为的结果。

现在，让我们来看一个使用决策树模型，但在集成方法中使用的模型，并且不容易发生过拟合。

随机森林

随机森林是一种集成学习方法，在训练过程中构建多个决策树，并将它们的结果结合在一起，以提高准确性和鲁棒性。它在分类和回归任务中都表现出色，通过袋装法和特征随机性减少了过拟合，并提高了预测性能。

集成学习是一种使用多个模型或专家，并通过某种方式将它们结合来解决手头问题的过程。它借用了集体智慧的方法，这一方法认为一群人的决策通常比其中任何一个人做出的决策更为准确。在机器学习中，这些单独的模型被称为基本学习器。单一模型可能因为过拟合数据集而表现不佳，但当我们将多个这样的模型结合起来时，它们可以形成一个强大的学习器。

袋装法（Bagging）是一种集成学习方法，在这种方法中，我们使用自助采样（从总体中反复有放回地抽取样本）来绘制数据集的不同子集，在每个子集上训练弱学习器，然后通过平均或投票的方式将它们结合在一起（分别用于回归和分类）。袋装法最适合于高方差、低偏差的弱学习器，而决策树就是一个非常适合与袋装法配合使用的成功候选者。从理论上讲，袋装法保持弱学习器的偏差水平不变，但减少了方差，从而提高了模型的表现。但如果弱学习器之间存在相关性，袋装法的效果将会受到限制。

2001 年，Leo Brieman 提出了随机森林，它在标准的袋装法基础上作出了重要修改，通过构建大量去相关的决策树来改进。为了确保在自助采样数据集上生成的所有树不相互关联，他提出略微修改树的构建过程。

参考文献检查：

随机森林的原始研究论文在 参考文献 部分被引用为参考文献 1。

在随机森林算法中，我们决定构建多少棵树。我们称之为 M 棵树。现在，对于每棵树，重复以下步骤：

1. 从训练数据集中抽取一个自助采样样本。

2. 从所有特征中随机选择 f 个特征。

3. 仅使用 f 个特征选择最佳分裂，并将节点分裂成两个子节点。

4. 重复 步骤 2 和 步骤 3，直到达到任何定义的停止标准。

这组 M 棵树就是随机森林。与常规树的关键区别在于每次分裂时特征的随机抽样，这增加了随机性并减少了不同树输出之间的相关性。在预测时，我们使用每棵 M 树来获得一个预测结果。对于回归问题，我们取它们的平均值，而对于分类问题，我们取多数投票。我们从随机森林中学习到的回归预测函数如下：

这里，Tt 是随机森林中 t^(th) 棵树的输出。

我们需要控制决策树复杂度的所有超参数在这里同样适用（来自 scikit-learn 的 RandomForestRegressor）。除了这些之外，我们还有两个其他重要的参数——集成中要构建的树的数量（n_estimators）和每次分裂时随机选择的特征数量（max_features）。

现在，让我们看看如何使用随机森林并在验证数据集中的一个样本家庭上评估拟合效果：

from sklearn.ensemble import RandomForestRegressor
model_config = ModelConfig(
    model=RandomForestRegressor(random_state=42, max_depth=4),
    name="Random Forest",
    # RandomForest is not affected by normalization
    normalize=False,
    # RandomForest in scikit-learn does not handle missing values
    fill_missing=True,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

让我们看看 RandomForestRegressor 的这一步预测。它比决策树更好，但不如线性模型。然而，我们应该记住，我们还没有调优模型，可能通过设置正确的超参数能够获得更好的结果。

现在，让我们看看使用随机森林生成的预测：

图 8.14：随机森林预测

就像决策树中的特征重要性一样，随机森林也有一个非常类似的机制来估计特征重要性。由于我们在随机森林中有很多棵树，我们会累积所有树中的分裂标准的减少，并得出一个单一的特征重要性。我们可以通过使用训练模型的 feature_importance_ 属性来访问这一点。让我们来看一下特征重要性：

图 8.15：决策树的特征重要性（前 15 名）

在这里，我们可以看到特征重要性与决策树非常相似。关于这种特征重要性的警告同样适用于这里。这只是一种粗略的方式来看模型在内部使用了哪些特征。

通常，随机森林在许多数据集上表现良好，几乎不需要调整，因此随机森林在机器学习中非常流行。由于随机森林不容易过拟合，这也增加了它的吸引力。但由于随机森林使用决策树作为弱学习器，决策树无法外推的问题也会传递给随机森林。

scikit-learn中实现的随机森林在树的数量和数据量较大时可能会变得较慢。XGBoost库中的XGBRFRegressor提供了一个随机森林的替代实现，特别是在较大的数据集上，由于XGBoost优化的算法和并行化能力，它可能会更快。此外，XGBRFRegressor使用的超参数与scikit-learn的随机森林类似，使得在调整模型时切换实现方式相对简单。在大多数情况下，这是一个直接替代，几乎能得到相同的结果。细微的差异来源于一些小的实现细节。我们在笔记本中也使用了这个变体。由于显而易见的运行时考虑，未来我们更倾向于使用这种变体。它还原生支持缺失值处理，避免了额外的预处理步骤。有关实现的更多细节以及如何使用它，可以参见xgboost.readthedocs.io/en/latest/tutorials/rf.html。

现在，让我们看一下一种最后的函数族，它是最强大的学习方法之一，并且在各种数据集上得到了极好的验证——梯度提升。

梯度提升决策树

提升法，像袋装法一样，是另一种集成方法，使用一些弱学习器来生成强大的模型组合。袋装法和提升法的关键区别在于弱学习器的组合方式。与袋装法在自助采样数据集上并行构建不同的模型不同，提升法按顺序使用弱学习器，每次弱学习器都会应用于反复修改过的数据版本。

为了理解加法函数的公式，让我们考虑这个函数：

我们可以将这个函数拆分成f1= 25，f2= x²，f3= cos(x)，并将F(x)重新写成如下形式：

F(x) = f1 + f2 + f3

这是我们在提升法中学习的加法集成函数。尽管从理论上讲，我们可以使用任何弱学习器，但决策树是最常见的选择。因此，我们使用决策树来探讨梯度提升是如何工作的。

早些时候，当我们讨论决策树时，我们看到一个具有M个分区的决策树，P[1]，P[2]，…，P[M]，其形式如下：

在这里，x 是输入，C[m] 是该区域的常数响应，P[m]，I 是一个函数，如果 ，则为 1；否则为 0。一个增强决策树模型是这些树的总和：

由于为所有树在集成中找到最佳分区 P 和常数值 c 是一个非常困难的优化问题，我们通常采用一种次优的阶段性解决方案，在构建集成时逐步优化每一步。在梯度提升中，我们使用损失函数的梯度来指导优化过程，因此得名“梯度提升”。

让我们在训练中使用的损失函数为 。由于我们观察的是一种逐步加法的功能形式，我们可以将 替换为 ，其中 是前 k-1 棵树的预测总和，Tk 是第 k 阶段树的预测。我们来看一下用于训练数据 D（有 N 个样本）的梯度提升学习过程：

1. 通过最小化损失函数来初始化模型，使用常数值：

• b[0] 是在第 0 次迭代中最小化损失函数的模型预测。在这个迭代中，我们还没有任何弱学习者，这个优化与任何特征无关。

• 对于平方误差损失，这相当于所有训练样本的平均值，而对于绝对误差损失，它是中位数。

1. 现在我们有了初步解，可以开始树构建过程。对于 k=1 到 M，我们必须执行以下操作：

   1. 对所有训练样本计算 ：

      • r[k] 是损失函数关于 F(x) 相对于上一轮迭代的导数。它也叫做伪残差。

      • 对于平方误差损失，这只是残差（）。

   2. 建立一个常规回归树，针对 r[k] 值，使用 M[k] 个分区或叶节点，P[mk]。

   3. 计算

      • 是当前阶段叶节点或分区值的缩放因子。

      • 是当前阶段决策树学习到的函数。

   4. 更新

      • 是收缩参数或学习率。

这种“增强”前一个弱模型误差的过程赋予了算法它的名字——梯度提升，其中梯度指的是前一个弱模型的残差。

提升（Boosting）通常是一种高方差算法。这意味着过拟合训练数据集的风险相当高，因此需要采取足够的措施来确保不会发生过拟合。在梯度提升树中，已经实现了多种正则化和容量约束方法。与往常一样，决策树为了减少容量以适应数据的所有关键参数在这里依然有效，因为弱学习器是决策树。除此之外，还有两个其他的关键参数——树的数量，M（在 scikit-learn 中为 n_estimators），以及学习率，（在 scikit-learn 中为 learning_rate）。

当我们在加法形式中应用学习率时，本质上是在缩小每个弱学习器的影响力，从而减少任何一个弱学习器对整体函数的影响。这最初被称为收缩（shrinkage），但现在，在所有流行的梯度提升树实现中，它被称为学习率。树的数量和学习率是高度互相依赖的。对于相同的问题，如果我们减少学习率，就需要更多的树。经验表明，较低的学习率可以改善泛化误差。因此，一种非常有效且便捷的方法是将学习率设置为非常低的值（<0.1），将树的数量设置为非常高的值（>5,000），并通过早停法训练梯度提升树。早停法是指在训练模型时，我们使用验证数据集来监控模型的外部样本表现。当外部样本误差不再减少时，我们停止向集成中添加更多的树。

许多实现采用的另一个关键技术是子采样（subsampling）。子采样可以在行和列上进行。行子采样类似于自助采样（bootstrapping），即每个候选模型在数据集的子样本上训练。列子采样类似于随机森林中的随机特征选择。两种技术都为集成引入了正则化效果，有助于减少泛化误差。像 XGBoost 和 LightGBM 等一些梯度提升树的实现直接在目标函数中实现了 L1 和 L2 正则化。

有许多回归梯度提升树的实现。以下是一些流行的实现：

• scikit-learn 中的 GradientBoostingRegressor 和 HistGradientBoostingRegressor

• 由 T Chen 开发的 XGBoost

• 来自 Microsoft 的 LightGBM

• 来自 Yandex 的 CatBoost

这些实现中的每一种都在标准的梯度提升算法上做出了从细微到非常根本的改动。我们在进一步阅读部分包含了一些资源，您可以通过这些资源了解这些差异，并熟悉它们支持的不同参数。

在本次练习中，我们将使用微软研究院的 LightGBM，因为它是最快且表现最佳的实现之一。LightGBM 和 CatBoost 还原生支持类别特征并处理缺失值。

参考检查：

XGBoost、LightGBM 和 CatBoost 的原始研究论文在参考文献部分分别被引用为2、3和4。

现在，让我们看看如何使用 LightGBM 并在验证数据集中的一个家庭样本上评估拟合效果：

from lightgbm import LGBMRegressor
model_config = ModelConfig(
    model=LGBMRegressor(random_state=42),
    name="LightGBM",
    # LightGBM is not affected by normalization
    normalize=False,
    # LightGBM handles missing values
    fill_missing=False,
)
y_pred, metrics, feat_df = evaluate_model(
    model_config,
    feat_config,
    missing_value_config,
    train_features,
    train_target,
    test_features,
    test_target,
) 

让我们看看从 LGBMRegressor 获得的单步预测。它已经明显优于我们迄今为止尝试的所有其他模型：

图 8.16：LightGBM 预测

与决策树中的特征重要性类似，梯度提升法（Gradient Boosting）也有一个非常相似的机制来估计特征重要性。集成模型的特征重要性由所有树中每个特征的划分标准减少量的平均值给出。这可以通过训练模型的feature_importance_属性在 scikit-learn API 中访问。让我们来看一下特征重要性：

图 8.17：LightGBM 特征重要性（前 15 名）

有多种方式可以从模型中获取特征重要性，每种实现的计算方式略有不同。这些方式是由参数控制的。最常见的提取方式（使用 LightGBM 术语）是 split 和 gain。如果我们选择 split，特征重要性就是特征用于划分树中节点的次数。另一方面，gain 是划分标准的总减少量，可以归因于任何特征。图 8.17 显示了 split，这是 LightGBM 中的默认值。我们可以看到，特征重要性的顺序与决策树或随机森林非常相似，几乎相同的特征占据了前三名。

梯度提升决策树（GBDTs）通常在表格数据和时间序列上表现非常好，回归任务也不例外。这个非常强大的模型通常是近年来几乎所有 Kaggle 时间序列预测竞赛获胜作品的一部分。虽然它是最好的机器学习模型家族之一，但它仍然有一些缺点：

• GBDT 是高方差算法，因此容易过拟合。这就是为什么在大多数成功实现的 GBDT 中以不同的方式应用了各种正则化技术。

• GBDT 通常需要更长的训练时间（尽管许多现代实现已经加速了这一过程），并且不像随机森林那样容易进行并行化。在随机森林中，我们可以并行训练所有的树，因为它们彼此独立。但在 GBDT 中，算法的顺序性限制了并行化。所有成功的实现都有一些巧妙的方式来在创建决策树时启用并行化。LightGBM 具有多种并行化策略，例如特征并行、数据并行和投票并行。有关这些策略的详细信息可以在lightgbm.readthedocs.io/en/latest/Features.html#optimization-in-distributed-learning找到，值得理解。该库的文档还包含一个有用的指南，帮助选择这些并行化策略之间的优先级，表格中有详细说明：

图 8.18：LightGBM 中的并行化策略

• 外推是 GBDT 的一个问题，就像所有基于树的模型一样。GBDT 在外推上有一些非常微弱的潜力，但并没有解决这个问题的方法。因此，如果你的时间序列具有一些强趋势，基于树的方法很可能无法捕捉到这个趋势。解决方法是使用去趋势的数据来训练模型，或者转向另一种模型类别。一种简单的去趋势方法是使用 AutoStationaryTransformer，我们在第六章《时间序列预测的特征工程》中讨论过。

总结一下，让我们看看这些机器学习模型所采用的度量标准和运行时间。如果你在本章中运行了笔记本，那么你会在其中找到以下总结表格：

图 8.18：示例家庭的度量标准和运行时间总结

一开始，我们就能看到所有尝试过的机器学习模型在所有度量标准上都优于基线，除了预测偏差。三种线性回归模型在 MAE、MASE 和 MSE 上表现相当，且正则化模型的运行时间略有增加。决策树的表现不尽如人意，但这通常是可以预期的。决策树需要稍微调整一下，以减少过拟合。随机森林（包括 scikit-learn 和XGBoost 实现）提高了决策树的表现，这是我们所期望的。这里需要注意的一个关键点是，XGBoost 实现的随机森林比 scikit-learn 实现的快了将近六倍。最后，LightGBM 在所有度量标准上表现最佳，且运行时间更短。

现在，这只是所有选定家庭中的一个。为了查看这些模型的表现，我们需要在所有选定的家庭上进行评估。

为多个家庭进行训练和预测

我们选择了一些模型（LassoCV、XGBRFRegressor 和 LGBMRegressor），这些模型在指标和运行时间方面表现更好，来处理我们验证数据集中的所有家庭。这个过程非常简单：循环遍历所有独特的组合，内部循环遍历不同的模型进行运行，然后进行训练、预测和评估。代码可以在Chapter08中的01-Forecasting_with_ML.ipynb笔记本里找到，位于对所有消费者进行机器学习预测标题下。你可以运行代码然后休息一下，因为这个过程大约需要不到一小时。笔记本还会计算指标，并包含一个总结表格，当你回来时它已经准备好。

现在我们来看一下总结：

图 8.19：验证数据集上所有家庭的聚合指标

在这里，我们可以看到，即使在聚合级别，我们使用的不同模型也按预期表现。笔记本还将验证集的预测结果保存到磁盘。

笔记本提醒：

我们还需要运行另一个名为01a-Forecasting_with_ML_for_Test_Dataset.ipynb的笔记本，位于Chapter08中。这个笔记本遵循相同的流程，生成预测，并计算测试数据集上的指标。

测试数据集的聚合指标如下（来自笔记本）：

图 8.20：测试数据集上所有家庭的聚合指标

在第六章，时间序列预测的特征工程中，我们在所有家庭上使用了AutoStationaryTransformer（而不是转换器模型，我们将在第十四章中学习）并保存了转化后的数据集。

使用 AutoStationaryTransformer

这个过程与我们在本章早些时候做的非常相似，唯一的区别是有一些小的变化。我们读取了转化后的目标数据，并将其与常规数据集连接，原始目标命名为energy_consumption，而转化后的目标命名为energy_consumption_auto_stat：

#Reading the missing value imputed and train test split data
train_df = pd.read_parquet(preprocessed/"block_0-7_train_missing_imputed_feature_engg.parquet")
auto_stat_target = pd.read_parquet(preprocessed/"block_0-7_train_auto_stat_target.parquet")
transformer_pipelines = joblib.load(preprocessed/"auto_transformer_pipelines_train.pkl")
#Reading in validation as test
test_df = pd.read_parquet(preprocessed/"block_0-7_val_missing_imputed_feature_engg.parquet")
# Joining the transformed target
train_df = train_df.set_index(['LCLid','timestamp']).join(auto_stat_target).reset_index() 

在定义FeatureConfig时，我们使用energy_consumption_auto_stat作为target，energy_consumption作为original_target。

笔记本提醒：

02-Forecasting_with_ML_and_Target_Transformation.ipynb和02a-Forecasting_with_ML_and_Target_Transformation_for_Test_Dataset.ipynb笔记本使用这些转化后的目标生成验证和测试数据集的预测。

让我们看看这些笔记本在转化数据上生成的总结指标：

图 8.21：验证数据集中所有家庭的目标转化后的汇总指标

目标转化后的模型表现不如原始模型。这可能是因为数据集没有强烈的趋势。

恭喜你顺利完成了这一章，它包含了大量的理论和实践内容。我们希望这能增强你对机器学习的理解，并提升你将这些现代技术应用于时间序列数据的能力。

总结

这一章非常实用，具有操作性，我们开发了一些标准代码来训练和评估多个机器学习模型。接着，我们回顾了几个关键的机器学习模型，如岭回归、套索回归、决策树、随机森林和梯度提升树，并探讨了它们背后的工作原理。为了巩固和强化所学知识，我们将所学的机器学习模型应用于伦敦智能电表数据集，并观察它们的表现。本章为接下来的章节做好了准备，届时我们将利用标准化代码和这些模型深入探索使用机器学习进行预测。

在下一章，我们将开始将不同的预测结果合并为一个单一的预测，并探讨组合优化和堆叠等概念，以实现最先进的结果。

参考文献

本章提供了以下参考文献：

1. Breiman, L. 随机森林，机器学习 45，5–32（2001）：doi.org/10.1023/A:1010933404324。

2. Chen, Tianqi 和 Guestrin, Carlos. (2016). XGBoost: A Scalable Tree Boosting System. 2016 年 ACM SIGKDD 国际知识发现与数据挖掘会议论文集 (KDD ‘16)。计算机协会，纽约，NY，美国，785–794: doi.org/10.1145/2939672.2939785。

3. Ke, Guolin 等（2017），LightGBM: 高效的梯度提升决策树。神经信息处理系统进展，3149-3157 页：dl.acm.org/doi/pdf/10.5555/3294996.3295074。

4. Prokhorenkova, Liudmila 等（2018），CatBoost: 带有类别特征的无偏增强。2018 年神经信息处理系统国际会议论文集（NIPS’18）：dl.acm.org/doi/abs/10.5555/3327757.3327770。

进一步阅读

若想深入了解本章所涉及的主题，请查看以下资源：

• L1 和 L2 正则化的区别，作者：Terrence Parr：explained.ai/regularization/L1vsL2.html

• L1 范数与 L2 范数，作者：Aleksey Bilogur：www.kaggle.com/residentmario/l1-norms-versus-l2-norms

• 可解释性 – 打开黑箱，作者：Manu Joseph：deep-and-shallow.com/2019/11/13/interpretability-cracking-open-the-black-box-part-ii/

• 梯度提升算法 – 第三部分：XGBoost，作者：Manu Joseph：deep-and-shallow.com/2020/02/12/the-gradient-boosters-iii-xgboost/

• 梯度提升算法 – 第四部分：LightGBM，作者：Manu Joseph：deep-and-shallow.com/2020/02/21/the-gradient-boosters-iii-lightgbm/

• 梯度提升算法 – 第五部分：CatBoost，作者：Manu Joseph：deep-and-shallow.com/2020/02/29/the-gradient-boosters-v-catboost/

• 梯度提升算法 – 第二部分：正则化贪婪森林，作者：Manu Joseph：deep-and-shallow.com/2020/02/09/the-gradient-boosters-ii-regularized-greedy-forest/

• LightGBM 分布式学习指南：lightgbm.readthedocs.io/en/latest/Parallel-Learning-Guide.html

加入我们社区的 Discord

加入我们社区的 Discord 空间，与作者和其他读者进行讨论：

packt.link/mts

第九章：集成和堆叠

在上一章中，我们查看了几种机器学习算法，并使用它们对伦敦智能电表数据集进行了预测。现在，我们对数据集中的所有家庭生成了多个预测，如何通过选择或结合这些不同的预测来得出一个单一的预测呢？最终，我们只能拥有一个预测，用于规划你正在预测的任务。这就是我们在本章中要做的——我们将学习如何利用组合和数学优化来得出一个单一的预测。

本章将涵盖以下主题：

• 结合预测的策略

• 堆叠或混合

技术要求

你需要按照本书前言中的说明设置 Anaconda 环境，以获得一个包含所有代码所需库和数据集的工作环境。任何额外的库将在运行笔记本时安装。

你需要在使用本章代码之前运行以下笔记本：

• 02-处理伦敦智能电表数据集.ipynb 在 Chapter02 中

• 01-设置实验框架.ipynb 在 Chapter04 中

• 02-使用 darts 的基准预测.ipynb 在 Chapter04 中

• 01-特征工程.ipynb 在 Chapter06 中

• 02-处理非平稳性.ipynb 在 Chapter07 中

• 02a-处理非平稳性-训练+验证.ipynb 在 Chapter07 中

• 00-单步回测基准.ipynb 在 Chapter08 中

• 01-使用机器学习进行预测.ipynb 在 Chapter08 中

• 01a-使用机器学习进行测试数据集预测.ipynb 在 Chapter08 中

• 02-使用目标转换进行预测.ipynb 在 Chapter08 中

• 02a-使用目标转换进行预测（测试）.ipynb 在 Chapter08 中

本章的代码可以在 github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-/tree/main/notebooks/Chapter09 找到。

结合预测

我们已经使用许多技术生成了预测——有些是单变量的，有些是机器学习的，等等。但最终，我们需要一个单一的预测，这意味着选择一个预测或结合多种预测。最简单的选择是选择在验证数据集中表现最好的算法，在我们的案例中是 LightGBM。我们可以将这种选择看作是另一个函数，它接受我们生成的预测作为输入并将它们合并成一个最终的预测。从数学角度来看，可以表示如下：

Y = F(Y[1], Y[2], …, Y[N])

在这里，F 是一个结合 N 个预测的函数。我们可以使用 F 函数来选择验证数据集中表现最好的模型。然而，这个函数可以非常复杂，选择一个合适的 F 函数，同时平衡偏差和方差是必须的。

笔记本提示：

要跟随代码进行操作，请在Chapter09文件夹中使用01-Forecast_Combinations.ipynb笔记本。

我们将从加载所有预测（包括验证和测试预测）以及到目前为止生成的所有预测对应的指标开始，并将它们合并到pred_val_df和pred_test_df中。接下来，我们必须使用pd.pivot重塑 DataFrame，以便获得我们想要的形状。到目前为止，我们一直在跟踪多个指标。但为了实现这个目标，我们需要选择一个指标。在这个练习中，我们选择 MAE 作为指标。验证指标可以合并并重塑成metrics_combined_df：

图 9.1：重塑后的预测 DataFrame

现在，让我们看看结合预测的几种不同策略。

最佳拟合

这种选择最佳预测的策略迄今为止是最流行的，其方法非常简单：根据验证指标为每个时间序列选择最佳预测。这种策略已被许多自动化预测软件工具所采用，称其为“最佳拟合”预测。该算法非常简单：

1. 使用验证数据集找到每个时间序列的最佳预测。

2. 对于每个时间序列，选择相同模型的测试数据集中的预测。

我们可以轻松做到这一点：

# Finding the lowest metric for each LCLid
best_alg = metrics_combined_df.idxmin(axis=1)
#Initialize two columns in the dataframe
pred_wide_test["best_fit"] = np.nan
pred_wide_test["best_fit_alg"] = ""
#For each LCL id
for lcl_id in tqdm(pred_wide_test.index.get_level_values(0).unique()):
    # pick the best algorithm
    alg = best_alg[lcl_id]
    # and store the forecast in the best_fit column
    pred_wide_test.loc[lcl_id, "best_fit"] = pred_wide_test.loc[lcl_id, alg].values
    # also store which model was chosen for traceability
    pred_wide_test.loc[lcl_id, "best_fit_alg"] = alg 

这将创建一个名为best_fit的新列，其中包含根据我们讨论的策略选择的预测。现在，我们可以评估这个新的预测，并获得测试数据集的指标。下表显示了最佳单一模型（LightGBM）和新策略—best_fit：

图 9.2：最佳拟合策略的汇总指标

在这里，我们可以看到最佳拟合策略整体表现优于单一的最佳模型。然而，这种策略的一个缺点是它的基本假设——在验证期表现最佳的模型也将在测试期表现最佳。它没有考虑其他预测模型等的对冲策略。考虑到时间序列的动态特性，这并不总是最佳策略。这种方法的另一个缺点是最终预测的不稳定性。

当我们在实际环境中使用这种规则时，每周重新训练并重新运行最佳拟合时，任何时间序列的预测可能会在不同的预测模型之间来回跳动，产生截然不同的预测结果。因此，最终的预测表现出很大的周间不稳定性，这会妨碍我们使用这些预测进行的后续操作。我们可以考虑一些没有这种不稳定性的其他技术。

集中趋势度量

另一种显著的策略是使用均值或中位数来合并预测。这是一个函数，F，它不依赖于验证指标。这既是这种方法的吸引力，也是一种困扰。由于我们根本没有使用验证指标，所以不可能过度拟合验证数据。但另一方面，由于没有任何验证指标的信息，我们可能会包括一些非常差的模型，这会拉低整体预测效果。然而，经验表明，这种简单的平均或中位数合并方法已被证明是一种非常强大的预测组合方法，且很难被超越。让我们看看如何实现这一方法：

# ensemble_forecasts is a list of column names(forecast) we want to combine
pred_wide_test["average_ensemble"] = pred_wide_test[ensemble_forecasts].mean(axis=1)
pred_wide_test["median_ensemble"] = pred_wide_test[ensemble_forecasts].median(axis=1) 

上述代码将创建两个新列，分别称为average_ensemble和median_ensemble，用于存储合并后的预测。现在，我们可以评估这个新的预测，并获取测试数据集的指标。下表显示了最佳单个模型（LightGBM）和新的策略：

图 9.3：均值和中位数策略的汇总指标

在这里，我们可以看到，无论是均值策略还是中位数策略，都没有比最好的单个模型整体表现得更好。这可能是因为我们包含了像 Theta 和 FFT 这样的模型，它们的表现远不如其他机器学习方法。但由于我们没有使用验证数据集中的任何信息，所以我们并不知道这一点。我们可以做一个例外，假设我们使用验证指标来选择哪些模型包含在平均值或中位数中。但我们必须小心，因为现在我们越来越接近于假设在验证期有效的模型也会在测试期有效。

这里有几种手动技术可以使用，例如修剪（丢弃表现最差的模型）和筛选（只选择表现最好的几个模型）。虽然这些方法有效，但有些主观性，尤其是在我们需要从成百上千个模型中选择时，它们变得难以使用。

如果我们考虑这个问题，本质上是一个组合优化问题，我们需要选择最佳的模型组合来优化我们的指标。如果我们考虑通过取平均值来合并不同的预测，从数学角度看，可以表示为：

在这里，L 是我们尝试最小化的损失或指标。在我们的例子中，我们选择的是 MAE。 是每个基础预测的二进制权重。最后， 是N个基础预测集合，Y 是时间序列的实际观测值。

但与纯粹的优化不同，纯粹优化中没有偏差和方差的概念，我们需要一个能够泛化的最优解。因此，选择训练数据中的全局最小值并不可取，因为那样可能会进一步过拟合训练数据集，增加最终模型的方差。在这种最小化中，我们通常使用样本外预测，在这种情况下可以是验证期间的预测。

最直接的解决方案是找到* w *，使其在验证数据上最小化此函数。但这种方法有两个问题：

• 随着基准预测数量N的增加，可能的候选（基准预测的不同组合）呈指数增长。这很快就变得计算上难以处理。

• 在验证期间选择全局最小值可能不是最佳策略，因为可能会导致验证期间的过拟合。

现在，让我们来看看一些基于启发式的解决方案来解决这个组合优化问题。

启发式问题解决是一种策略，利用经验法则或捷径快速找到解决方案，即使这些解决方案可能不是最优的。当精确解法在计算上昂贵或耗时时，启发式方法可以非常有用。然而，在某些情况下，它们可能导致次优甚至错误的解决方案。

启发式方法通常与其他问题解决方法（如元启发式方法）结合使用，以提高搜索过程的效率和效果。元启发式方法是高层次的、与问题无关的策略，用于解决优化问题。它们提供了一个框架，用于开发启发式算法，能够高效地探索复杂的搜索空间并找到近似最优解。与传统优化方法不同，元启发式方法通常从自然现象或生物过程中汲取灵感。

元启发式方法的常见例子包括遗传算法（灵感来自自然选择）、模拟退火（灵感来自冶金学）、粒子群优化（灵感来自鸟群聚集）和蚁群优化（灵感来自蚂蚁觅食）。这些方法采用概率或随机方法来平衡探索与开发，使其能够避免陷入局部最优解，发现潜在的更好解决方案。

简单的爬山法

在讨论决策树以及梯度提升树时，我们简要介绍了贪婪算法。贪婪优化是一种启发式方法，通过逐步构建解决方案，在每一步选择一个局部最优解。在这两种机器学习模型中，我们采用了贪婪的、逐步的方式来解决计算上不可行的优化问题。为了选择最佳子集，给我们提供最佳的预测组合，我们可以采用一种简单的贪婪算法，称为爬坡算法。如果我们将目标函数的曲面看作一座山，为了找到最大值，我们需要爬上这座山。顾名思义，爬坡算法逐步上升，在每一步中，它选择最优路径，从而增加目标函数的值。下面的示意图可以帮助更清晰地理解。

图 9.4：一维目标的爬坡算法示意图

我们可以从 图 9.4 中看到，目标函数（我们需要优化的函数）有多个峰值（山峰），而在爬坡算法中，我们“爬”上山，逐步到达峰顶。我们还需要记住，根据我们开始爬坡的位置不同，可能会达到目标中的不同点。在 图 9.4 中，如果我们从 A 点开始爬坡，我们到达局部最优解，并错过全局最优解。现在，让我们看看该算法是如何以更严格的方式运作的。

在这里，C 是一组候选解（基础预测），O 是我们希望最小化的目标。简单爬坡算法如下：

1. 初始化起始解，C[best]，作为在 O 中给出最小值的候选解，即 O[best]，并从 C 中移除 C[best]。

2. 当 C 的长度大于 0 时，执行以下操作：

   1. 通过将 C[best] 中的基础预测与 C 中的每个元素进行平均，评估 C 中的所有成员，并选择最佳成员 (C[stage best])，将其添加到 C[best] 中，以最小化目标函数 O（即 O[stage best]）。

   2. 如果 O[stage best] > O[best]，则执行以下操作：

      1. C[best] = C[best] U C[stage best]。

      2. O[best] = O[stage best]。

      3. 从 C 中移除 C[stage best]。

   3. 否则，退出。

在运行结束时，我们得到 C[best]，这是通过贪婪优化得到的最佳预测组合。我们在 src.forecasting.ensembling.py 中实现了这个功能，位于 greedy_optimization 函数下。该函数的参数如下：

• objective：这是一个可调用函数，接受一个字符串列表作为候选解，并返回一个 float 类型的目标值。

• candidates：这是一个候选解列表，将被包括在优化中。

• verbose：一个标志，指定是否打印进度。

该函数返回一个元组，其中包含作为字符串列表的最佳解和通过优化得到的最佳评分。

让我们看看如何在我们的示例中使用这个算法：

1. 导入所有必需的库/函数：

   # Used to partially construct a function call
   from functools import partial
   # calculate_performance is a custom method we defined to calculate the MAE provided a list of candidates and prediction dataframe
   from src.forecasting.ensembling import calculate_performance, greedy_optimization 
   

2. 定义目标函数并运行贪心优化：

   # We partially construct the function call by passing the necessary parameters
   objective = partial(
       calculate_performance, pred_wide=pred_wide_val, target="energy_consumption"
   )
   # ensemble forecasts is the list of candidates
   solution, best_score = greedy_optimization(objective, ensemble_forecasts) 
   

3. 一旦我们得到了最佳解，我们可以在测试数据框中创建组合预测：

   pred_wide_test["greedy_ensemble"] = pred_wide_test[solution].mean(axis=1) 
   

一旦我们运行此代码，我们将在预测数据框（DataFrame）中得到名称为greedy_ensemble的组合预测。最优解中的候选模型包括 LightGBM、Lasso 回归和 LightGBM_auto_stat。接下来，让我们评估结果并查看汇总的度量指标：

图 9.5：基于简单爬山的集成方法汇总指标

如我们所见，简单的爬山算法的表现优于我们迄今为止看到的任何单一模型或集成技术。在这种情况下，贪心算法似乎运作良好。现在，让我们了解爬山算法的几个局限性，如下所示：

• 运行时考虑：由于简单的爬山算法需要在每一步评估所有候选者，这可能会导致运行时瓶颈。如果候选者数量很大，这种方法可能会花费更多时间才能完成。

• 短视性：爬山优化是短视的。在优化过程中，它每一步总是选择最佳的选项。有时，通过在某一步选择一个稍差的解决方案，我们可能会得到一个更好的整体解决方案。

• 只前进：爬山算法是一个只前进的算法。一旦一个候选者被纳入解决方案，我们就不能回头将其移除。

贪心算法并不总能为我们找到最优解，尤其是在需要组合多个模型时。因此，让我们来看看一种小的变种——爬山算法，它试图克服贪心算法的一些局限性。

随机爬山

简单的爬山算法和随机爬山算法的关键区别在于候选者的评估。在简单的爬山中，我们会评估所有可能的选项并从中挑选最佳的一个。然而，在随机爬山中，我们会随机挑选一个候选者，如果它比当前解更好，就将其添加到解决方案中。换句话说，在爬山算法中，我们总是逐步向上移动，但在随机爬山中，我们会神奇地瞬移到目标函数的不同点，检查自己是否比之前更高。这种加入随机性的做法有助于优化算法避免陷入局部最大值/最小值，但也引入了相当大的不确定性，可能无法达到任何最优解。接下来，我们来看看这个算法。

在这里，C 是候选集（基础预测），O 是我们希望最小化的目标，N 是我们希望运行优化的最大迭代次数。随机爬山算法如下：

1. 初始化起始解，C[best]，作为候选者。可以通过随机挑选一个候选者或选择表现最好的模型来完成。

2. 将 C[best] 的目标函数值 O 设置为 O[best]，并从 C 中移除 C[best]。

3. 对 N 次迭代重复以下步骤：

   1. 从 C 中随机抽取一个样本，将其加入 C[best]，并存储为 C[stage]。

   2. 在目标函数 O 上评估 C[stage]，并将其存储为 O[stage]。

   3. 如果 O[stage] > O[best]，则执行以下操作：

      1. C[best] = C[best] U C[stage]

      2. O[best] = O[stage]。

      3. 从 C 中移除 C[best]。

在运行结束时，我们得到的是 C[best]，它是通过随机爬山算法获得的最佳预测组合。我们已经在 src.forecasting.ensembling.py 中的 stochastic_hillclimbing 函数下实现了这一方法。该函数的参数如下：

• objective：这是一个可调用的函数，它接收一个包含候选字符串的列表并返回一个 float 类型的目标值。

• candidates：这是一个候选列表，将被包含在优化过程中。

• n_iterations：执行爬山算法的迭代次数。如果未给定该值，则使用启发式方法（候选数量的两倍）来设置该值。

• init：决定用于初始解的策略，可以是 random 或 best。

• verbose：一个标志，用来指定是否打印进度。

• random_state：一个种子，用于获得可重复的结果。

该函数返回一个元组，包含作为字符串列表的最佳解和通过优化获得的最佳得分。

这可以与 greedy_optimization 以非常相似的方式使用。我们这里只展示不同的部分，完整代码可在笔记本中查看：

from src.forecasting.ensembling import stochastic_hillclimbing
# ensemble forecasts is the list of candidates
solution, best_score = stochastic_hillclimbing(
    objective, ensemble_forecasts, n_iterations=10, init="best", random_state=9
) 

一旦我们运行这段代码，我们将在预测 DataFrame 中得到一个名为 stochastic_hillclimb__ensemble 的组合预测。作为最佳解的一部分的候选模型包括 LightGBM、Lasso 回归 _auto_stat、LightGBM_auto_stat 和 Lasso 回归。现在，让我们评估结果并查看聚合指标：

图 9.6：基于随机爬山的集成的聚合指标

随机爬山算法的效果不比贪婪算法更好，但却优于均值、媒体和最佳拟合集成。我们之前讨论了简单爬山法的三个缺点——运行时考虑、短视性和仅向前搜索。随机爬山解决了运行时考虑的问题，因为我们并没有评估所有的组合并选择最佳，而是通过随机评估组合并一旦找到一个表现更好的解就将其加入集成。它在一定程度上解决了短视性问题，因为算法中的随机性可能会导致每个阶段选择一个次优解，但它仍然只选择比当前解更好的解。

现在，让我们来看一下另一个改进版的爬山算法，它也解决了这个问题。

模拟退火

模拟退火是对爬山算法的一种改进，灵感来源于一种物理现象——退火固体。退火是将固体加热到预定温度（通常高于其再结晶温度，但低于其熔点），保持一段时间，然后冷却（可以慢慢冷却，也可以通过在水中淬火来快速冷却）。

这样做是为了确保原子达到新的全局最小能量状态，这会在某些金属中引入期望的性质，比如铁。

1952 年，Metropolis 提出了模拟退火作为一种优化技术。退火类比也适用于优化背景。当我们说加热系统时，实际上是指我们鼓励随机扰动。因此，当我们以高温开始优化时，算法会探索空间并得出问题的初始结构。随着温度的降低，结构会被细化，从而得到最终解决方案。这种技术有助于避免陷入任何局部最优解。局部最优解是目标函数表面上的极值，它比附近的其他值更好，但可能不是最优解。进一步阅读部分包含了简洁解释局部最优解和全局最优解的资源。

现在，让我们来看看这个算法。

这里，C是候选集（基本预测），O是我们要最小化的目标，N是我们希望运行优化的最大迭代次数，T[max]是最大温度，是温度衰减。模拟退火算法如下：

1. 初始化一个起始解，C[best]，作为候选解。这可以通过随机选择一个候选解或选择表现最好的模型来完成。

2. 为C[best]设置目标函数的值，O，作为O[best]，并从C中移除C[best]。

3. 将当前温度设置为t = T[max]。

4. 对N次迭代重复执行此操作：

   1. 从C中随机抽取一个样本，添加到C[best]，并将其存储为C[stage]。

   2. 在目标函数上评估C[stage]，O，并将其存储为O[stage]。

   3. 如果 O[stage] > O[best]，则执行以下操作：

      1. C[best].= C[best] U C[stage]

      2. O[best] = O[stage]

      3. 从C中移除C[best]。

   4. 否则，执行以下操作：

      1. 计算接受概率，。

      2. 从 0 到 1 之间随机抽取一个样本，记为 p。

      3. 如果 p < s，则执行以下操作：

      4. C[best] = C[best] U C[stage]

      5. O[best] = O[stage]。

      6. 从C中移除C[best]。

   5. t = t - （对于线性衰减）和 t = t/（对于几何衰减）。

   6. 当C为空时退出。

在运行结束时，我们得到了 C[best]，这是通过模拟退火获得的最佳预测组合。我们在 src.forecasting.ensembling.py 文件中的 simulated_annealing 函数下提供了该实现。将温度设置为合适的值对于算法的良好运行至关重要，而且通常是最难设置的超参数。更直观地，我们可以将温度视为开始时接受较差解的概率。在实现中，我们还使得可以输入接受较差解的起始和结束概率。

1989 年，D.S. Johnson 等人提出了一种从给定的概率范围估算温度范围的过程。该过程已经在 initialize_temperature_range 中实现。

总结一下，算法从随机解开始，并评估每个解的好坏。然后，它不断尝试新解，有时接受较差的解以避免陷入局部最优解，但随着时间的推移，它接受较差解的可能性会降低，因为它在“冷却”（就像金属冷却并硬化一样）。它重复这一过程，直到选项用尽或温度变得太低，最终保留找到的最佳解。

参考检查：

D.S. Johnson 的研究论文，标题为 模拟退火优化：实验评估；第一部分，图划分，在 参考文献 部分被引用为参考文献 1。

simulated_annealing 函数的参数如下：

• objective：这是一个可调用的函数，接受一个字符串列表作为候选项，并返回一个 float 类型的目标值。

• candidates：这是一个候选列表，用于包含在优化中。

• n_iterations：模拟退火运行的迭代次数。这个参数是必需的。

• p_range：起始和结束概率的元组。这是在模拟退火中接受较差解的概率。温度范围（t_range）将在优化过程中从 p_range 推断得出。

• t_range：如果我们想直接设置温度范围为一个元组（起始，结束），可以使用这个参数。如果设置了该值，p_range 会被忽略。

• init：这个参数决定了用于初始解的策略。可以是 random 或 best。

• temperature_decay：指定温度衰减的方式。可以是 linear 或 geometric。

• verbose：一个标志，指定是否打印进度。

• random_state：用于获取可重复结果的种子。

该函数返回一个元组，包含作为字符串列表的最佳解决方案和通过优化获得的最佳得分。

这可以像其他组合预测的方式一样使用。我们将在这里展示不同之处。完整代码可在笔记本中查看：

from src.forecasting.ensembling import simulated_annealing
# ensemble forecasts is the list of candidates
solution, best_score = simulated_annealing(
    objective,
    ensemble_forecasts,
    p_range=(0.5, 0.0001),
    n_iterations=50,
    init="best",
    temperature_decay="geometric",
    random_state=42,
) 

一旦我们运行这段代码，我们将在预测 DataFrame 中得到一个名为simulated_annealing_ensemble的组合预测。作为最优解的一部分的候选模型包括 LightGBM、Lasso 回归 _auto_stat、LightGBM_auto_stat 和 XGB 随机森林。让我们评估一下结果并查看汇总指标：

图 9.7：基于模拟退火的集成的汇总指标

模拟退火似乎比随机爬山表现得更好。我们之前讨论过简单爬山算法的三个缺点——运行时考虑、目光短浅以及仅向前搜索。模拟退火解决了运行时问题，因为我们不是评估所有组合并选择最优的，而是随机评估组合，并在发现更好的解时立即将其添加到集成中。它也解决了目光短浅的问题，因为通过使用温度，我们在优化的初期也会接受略差的解。然而，它仍然是一个仅向前搜索的过程。

到目前为止，我们已经看过组合优化问题，因为我们说过。但如果我们可以放宽这个约束，使得 。 （实数），那么组合优化问题可以放宽为一个一般的数学优化问题。让我们看看如何做到这一点。

最优加权集成

之前，我们将我们试图解决的优化问题定义为如下：

这里，L是我们试图最小化的损失或指标。在我们的例子中，我们选择了 MAE。 是 N 个基本预测集合，而 Y 是时间序列的真实观测值。我们不再定义 ，而是让 成为每个基本预测的连续权重。通过这个新放宽的约束，组合变成了不同基本预测之间的加权平均。现在，我们正在看不同预测的软混合，而不是基于硬选择的组合优化（这也是我们一直使用的方法）。

这是一个可以使用 scipy 中现成算法解决的优化问题。让我们看看如何使用 scipy.optimize 来解决这个问题。

首先，我们需要定义一个损失函数，该函数接受一组作为列表的权重，并返回我们需要优化的指标：

def loss_function(weights):
        # Calculating the weighted average
        fc = np.sum(pred_wide[candidates].values * np.array(weights), axis=1)
        # Using any metric function to calculate the metric
        return metric_fn(pred_wide[target].values, fc) 

现在，我们所需要做的就是用必要的参数调用scipy.optimize。让我们学习如何做这件事：

from scipy import optimize
opt_weights = optimize.minimize(
        loss_function,
        # set x0 as initial values, which is a uniform distribution over all the candidates
        x0=[1 / len(candidates)] * len(candidates),
        # Set the constraint so that the weights sum to one
        constraints=({"type": "eq", "fun": lambda w: 1 - sum(w)}),
        # Choose the optimization technique. Should be gradient-free and bounded.
        method="SLSQP",
        # Set the lower and upper bound as a tuple for each element in the candidate list.
        # We set the maximum values between 1 and 0
        bounds=[(0.0, 1.0)] * len(candidates),
        # Set the tolerance for termination
        options={"ftol": 1e-10},
    )["x"] 

优化通常很快，我们会得到作为浮动点数的权重列表。我们将其包装在src.forecasting.ensembling.py中的一个名为find_optimal_combination的函数中。该函数的参数如下：

• candidates：这是待纳入优化的候选项列表。它们的返回顺序与返回的权重顺序相同。

• pred_wide：这是我们需要学习权重的预测数据框。

• target：这是目标列的名称。

• metric_fn：这是任何具有metric(actuals, pred)签名的可调用对象。

该函数返回最优权重，作为一个浮动点数列表。让我们看看通过验证预测学习到的最优权重是什么：

图 9.8：通过优化学习得到的最优权重

在这里，我们可以看到优化自动学会忽略FFT、Theta、XGB Random Forest和XGB Random Forest_auto_stat，因为它们对集成模型贡献不大。它还学会了为每个预测分配一些非零的权重。这些权重已经与我们之前讨论的技术选择相似。现在，我们可以使用这些权重来计算加权平均值，并将其称为optimal_combination_ensemble。

聚合结果应该如下所示：

图 9.9：基于最优组合的集成的聚合指标

在这里，我们可以看到，这种软性混合预测的表现远远好于所有基于硬性选择的集成方法，在所有三个指标上都有明显的提升。

在我们讨论的所有技术中，我们使用的是 MAE 作为目标函数。但我们也可以使用任何指标、指标组合，甚至是带正则化的指标作为目标函数。当我们讨论随机森林时，我们提到去相关树对提高性能至关重要。一个非常相似的原则也适用于选择集成方法。拥有去相关的基础预测为集成模型增值。因此，我们可以使用任何多样性度量来对我们的指标进行正则化。例如，我们可以使用相关性作为度量，并创建一个正则化指标用于这些技术。Chapter09文件夹中的01-Forecast_Combinations.ipynb笔记本包含一个附加部分，展示了如何做到这一点。

我们一开始讨论的是如何通过数学公式来组合预测：

Y = F(Y[1], Y[2], …, Y[N])

在这里，F是将N个预测值组合起来的函数。

我们在寻找将该函数作为优化问题的解决方案时，使用了均值或中位数等方法来组合这些指标。但我们也看到了另一种从数据中学习这个函数F的方式，不是吗？让我们看看怎么做。

堆叠与混合

本章一开始我们讨论了机器学习算法，这些算法从一组输入和输出中学习一个函数。在使用这些机器学习算法的过程中，我们学习了预测时间序列的函数，现将其称为基预测。

为什么不使用相同的机器学习范式来学习我们也想学习的这个新函数 F 呢？

这正是堆叠中所做的（通常称为堆叠泛化），我们在一些基学习器的预测结果上训练另一个学习算法来结合这些预测。这种二级模型通常被称为堆叠模型或元模型。通常，这种元模型的表现与基学习器相当，甚至更好。这与混合（blending）非常相似，唯一的区别在于我们分割数据的方式。

尽管这一思想最早由沃尔珀特（Wolpert）于 1992 年提出，但莱奥·布雷曼（Leo Breiman）在 1996 年的论文《堆叠回归》（Stacked Regressions）中正式化了这一概念，成为如今的应用方式。并且在 2007 年，M. J. Van der Laan 等人建立了这一技术的理论基础，并提供了证明，表明这种元模型的表现至少与基学习器一样好，甚至更好。

参考检查：

莱奥·布雷曼（1996 年）和马尔科·J·范德兰（2007 年）等人的研究论文在参考文献部分被标记为2和3。

这是在机器学习竞赛中非常流行的一种技术，比如 Kaggle，被认为是机器学习从业者之间的一种秘密技巧。我们还讨论了其他一些技术，比如袋装法（bagging）和提升法（boosting），它们通过将基学习器组合成更复杂的模型来改进效果。但这些技术要求基学习器是一个弱学习器。而堆叠（stacking）则不同，因为堆叠尝试将一组多样化的强学习器进行组合。

堆叠的直觉是，不同的模型或函数族学习输出函数的方式略有不同，捕捉了问题的不同特征。例如，一个模型可能很好地捕捉了季节性，而另一个模型则可能更好地捕捉了与外生变量的某种交互。堆叠模型将能够将这些基模型结合成一个模型，其中一个模型关注季节性，另一个模型关注交互。这是通过让元模型学习基模型的预测结果来实现的。但为了防止数据泄漏并避免过拟合，元模型应在样本外的预测数据上进行训练。如今有两种小变体的技术——堆叠和混合。

堆叠是指元模型在整个训练数据集上进行训练，但使用的是样本外预测。堆叠过程包括以下步骤：

1. 将训练数据集分割成k部分。

2. 在k-1部分上迭代训练基本模型，在k^(th)部分上进行预测，并保存预测结果。完成此步骤后，我们有了来自所有基本模型的训练数据集的样本外预测。

3. 在这些预测上训练一个元模型。

混合与此类似，但在生成样本外预测的方式上略有不同。混合涉及以下步骤：

1. 将训练数据集分为两部分——训练和保留。

2. 在训练数据集上训练基本模型并在保留数据集上进行预测。

3. 在具有基本模型预测结果的验证数据集上训练元模型。

直觉上，我们可以看到堆叠可以工作得更好，因为它使用一个更大的数据集（通常是所有训练数据）作为样本外预测，所以元模型可能更加泛化。但是有一个警告：我们假设整个训练数据是独立同分布（iid）。这通常是一个很难在时间序列中满足的假设，因为数据生成过程可以随时改变（逐渐或急剧）。如果我们知道数据分布随时间发生了显著变化，那么混合保留期（通常是数据集的最近部分）更好，因为元模型只学习最新的数据，从而尊重数据分布的时间变化。

我们可以包含作为基本模型的模型数量没有限制，但通常会达到一个平台，额外的模型不会对堆叠集成产生太大的影响。我们还可以添加多个堆叠级别。例如，假设有四个基本学习器：B[1,] B[2,] B[3] 和 B[4]。我们还训练了两个元模型 M[1] 和 M[2]，在基本模型上。现在，我们可以在 M[1] 和 M[2] 的输出上训练第二级元模型 M，并将其用作最终预测。我们可以使用pystacknet Python 库（github.com/h2oai/pystacknet），这是一个名为stacknet的旧库的 Python 实现，以便轻松创建多级（或单级）堆叠集成的过程。

另一个要牢记的关键点是我们通常用作元模型的模型类型。假设大部分学习已经由基本模型完成，这些基本模型是用于预测的多维数据的模式。因此，元模型通常是简单模型，例如线性回归、决策树，甚至比基本模型低得多的随机森林。另一种思考这个问题的方式是从偏差和方差的角度来看。堆叠可能会过拟合训练或留出集，并通过包含具有更大灵活性或表达能力的模型族，我们正在促使这种过拟合发生。进一步阅读部分包含了一些链接，从通用机器学习的角度解释了不同的堆叠技术。

现在，让我们快速看看如何在我们的数据集中使用这个：

from sklearn.linear_model import LinearRegression
stacking_model = LinearRegression()
# ensemble_forecasts is the list of candidates
stacking_model.fit(
    pred_wide_val[ensemble_forecasts], pred_wide_val["energy_consumption"]
)
pred_wide_test["linear_reg_blending"] = stacking_model.predict(
    pred_wide_test[ensemble_forecasts]
) 

这将为线性回归保存混合预测为linear_reg_blending。我们可以使用相同的代码，但交换模型以尝试其他模型。

最佳实践：

当存在许多基本模型并且我们想要进行隐式基本模型选择时，我们可以选择其中一个正则化线性模型，例如岭回归或套索回归。在他最初的论文中，Breiman 提出了堆叠回归，建议使用具有正系数且没有截距的线性回归作为元模型。他认为这样可以理论上保证堆叠模型至少与任何最佳个体模型一样好。但在实践中，我们可以在实验中放宽这些假设。没有截距的非负回归与我们之前讨论过的最佳加权集成非常接近。最后，如果我们正在评估多个堆叠模型以选择哪个效果好，我们应该要么使用单独的验证数据集（而不是训练-验证-测试分割，我们可以使用训练-验证-验证元-测试分割），要么使用交叉验证估计。如果我们只是选择在测试数据集上表现最好的堆叠模型，那么我们就是在测试数据集上过拟合了。

现在，让我们看看混合模型在我们的测试数据上的表现：

图 9.10：混合模型的聚合指标

在这里，我们可以看到简单的线性回归学习了一个比我们任何平均集成方法都要好得多的元模型。而 Huber 回归（这是一种直接优化 MAE 的方法）在 MAE 基准测试上表现得更好。然而，请记住这并非普遍适用，必须针对遇到的每个问题进行评估。选择要优化的指标和要用于组合的模型会产生很大的差异。通常，简单的平均集成是组合模型的一个非常可观的基准。

Huber 回归是线性回归的另一种版本（如岭回归和套索回归），其损失函数是平方损失（用于常规线性回归）和绝对损失（用于 L1 方法）的组合。对于小残差，它表现得像平方损失，而对于大残差，它表现得像绝对损失。这使得它对异常值不太敏感。Scikit-Learn 提供了 HuberRegressor (scikit-learn.org/stable/modules/generated/sklearn.linear_model.HuberRegressor.html)，用于实现这一方法。

附加阅读：

还有其他更具创新性的方式来结合基础预测。这是一个活跃的研究领域。进一步阅读部分包含了两个非常相似的想法的链接。基于特征的预测模型平均法（FFORMA）从时间序列中提取一组统计特征，并用它来训练一个机器学习模型，预测基础预测应如何加权结合。另一种技术（用于快速且可扩展的时间序列超参数调整的自监督学习），来自 Facebook（Meta）研究，训练一个分类器，预测给定一组从时间序列中提取的统计特征时，哪个基础学习器表现最好。

摘要

延续上一章中的实用课程，我们又完成了一个动手实践的课程。在本章中，我们从上一章的不同机器学习模型中生成了预测结果。我们学会了如何将这些不同的预测结果结合成一个比任何单一模型表现都更好的预测。接着，我们探索了组合优化和堆叠/混合等概念，以实现最先进的结果。

在下一章中，我们将开始讨论全球预测模型，并探索策略、特征工程等，以便实现这种建模。

参考文献

本章提供了以下参考文献：

1. David S. Johnson，Cecilia R. Aragon，Lyle A. McGeoch，和 Catherine Schevon（1989），模拟退火优化：实验评估；第一部分，图形划分。运筹学，1989 年，第 37 卷，第 6 期，865-892：dx.doi.org/10.1287/opre.37.6.865

2. L. Breiman（1996），堆叠回归。机器学习 24，49-64：doi.org/10.1007/BF00117832

3. Mark J. van der Laan；Eric C. Polley；和 Alan E. Hubbard（2007），超级学习者。加利福尼亚大学伯克利分校生物统计学系工作论文系列。工作论文 222：biostats.bepress.com/ucbbiostat/paper222

进一步阅读

若想进一步了解本章所涉及的主题，请查阅以下资源：

• Kaggler 实践中的模型堆叠指南，由 Ha Nguyen 编写：datasciblog.github.io/2016/12/27/a-kagglers-guide-to-model-stacking-in-practice/

• Kai Ming Ting 和 Ian H. Witten (1997)，堆叠泛化：何时有效？：www.ijcai.org/Proceedings/97-2/Papers/011.pdf

• Pablo Montero-Manso, George Athanasopoulos, Rob J. Hyndman, Thiyanga S. Talagala (2020)，FFORMA: 基于特征的预测模型平均。《国际预测学杂志》，第 36 卷，第 1 期：robjhyndman.com/papers/fforma.pdf

• Peiyi Zhang 等人 (2021)，自监督学习用于快速且可扩展的时间序列超参数调优：www.ijcai.org/Proceedings/97-2/Papers/011.pdf

• 局部最优与全局最优：www.mathworks.com/help/optim/ug/local-vs-global-optima.html

加入我们的 Discord 社区

加入我们社区的 Discord 空间，与作者及其他读者进行讨论：

packt.link/mts