Python现代时间序列预测——分析和可视化时间序列数据

在上一章中，我们学习了从哪里获取时间序列数据集，以及如何使用 pandas 操作时间序列数据，处理缺失值等等。现在我们已经拥有了处理好的时间序列数据，是时候理解数据集了，数据科学家称之为探索性数据分析（EDA）。这是一个过程，数据科学家通过查看汇总统计、特征分布、可视化等来分析数据，试图发现可以在建模中利用的数据模式。在本章中，我们将探讨分析时间序列数据集的几种方法，一些专为时间序列设计的特定技术，并回顾一些时间序列数据的可视化技术。

在本章中，我们将涵盖以下主题：

• 时间序列的组成部分
• 时间序列数据的可视化
• 时间序列的分解
• 异常值的检测与处理
• 技术要求

技术要求：

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

您需要运行 Chapter02 文件夹中的 02-Preprocessing_London_Smart_Meter_Dataset.ipynb 笔记本。

本章的代码可在 github.com/PacktPublis… 找到。

时间序列的组成部分

在开始分析和可视化时间序列之前，我们需要了解时间序列的结构。任何时间序列都可能包含以下一些或全部组成部分：

• 趋势（Trend）
• 季节性（Seasonal）
• 周期性（Cyclical）
• 不规则（Irregular）

这些组成部分可以以不同的方式混合，但两种非常常见的假设方式是加法模型（Y = 趋势 + 季节性 + 周期性 + 不规则）和乘法模型（Y = 趋势 × 季节性 × 周期性 × 不规则），其中 Y 是时间序列。

趋势组成部分

趋势 是时间序列均值的长期变化。它是时间序列在特定方向上的平滑且稳定的运动。当时间序列向上移动时，我们称之为上升趋势或增长趋势；当其向下移动时，我们称之为下降趋势或减少趋势。在撰写本文时，如果我们考虑特斯拉多年来的收入，如下图所示，我们可以看到它在过去几年中一直在持续增长：

看前面的图，我们可以说特斯拉的收入呈上升趋势。趋势不需要是线性的，也可以是非线性的。

季节性组成部分

当一个时间序列表现出规律的、重复的、上下波动的变化时，我们称之为季节性。例如，零售销售通常在假日期间，尤其是在西方国家的圣诞节期间大幅上升。类似地，热带地区的电力消耗在夏季达到峰值，而在较冷的国家则在冬季达到峰值。在所有这些例子中，您可以看到每年重复的特定上下波动模式。另一个例子是太阳黑子，如下图所示：

周期性组成部分

周期性组成部分经常与季节性混淆，但由于一个非常微妙的区别，它们是不同的。像季节性一样，周期性组成部分也会围绕趋势线展示类似的上下波动模式，但这种周期的时间长度不是固定的，并且在一个大致的时间框架内会有所变化。一个很好的例子是经济衰退，它通常以10年的周期发生。然而，这并不是像时钟一样精确的，偶尔可能会少于或多于10年发生一次。

不规则组成部分

这是在从时间序列中去除趋势、季节性和周期性之后剩下的部分。传统上，这一部分被认为是不可预测的，也被称为残差、误差项或噪声项。在常见的经典统计学模型中，任何“模型”的目标都是捕捉其他所有组成部分，直到唯一没有被捕捉到的部分就是不规则组成部分。

在现代机器学习中，我们并不认为这一部分完全不可预测。我们通过使用外生变量来尝试捕捉这一部分。例如，零售销售的不规则组成部分可能可以通过他们进行的不同促销活动来解释。当我们获得这些额外的信息时，原本“不可预测”的部分开始变得可预测。但无论你向模型中添加多少额外变量，总会有一些部分是无法减少的误差，这是被称为不可约误差的部分。这是时间序列中的一部分，无论模型多么强大，或是你向模型中添加多少额外信息，都永远无法解释这一部分。

现在我们已经了解了时间序列的不同组成部分，接下来让我们看看如何可视化它们。

可视化时间序列数据

在第2章《获取与处理时间序列数据》中，我们学习了如何准备数据模型，这是分析新数据集的第一步。如果准备数据模型就像是接近一个你喜欢的人并做第一次接触，那么探索性数据分析（EDA）就像是与那个人约会。在这个阶段，你已经拥有了数据集，正在试图了解它，弄清楚它的规律，了解它的喜好和不喜欢的东西，等等。

EDA通常使用可视化技术来揭示模式、发现异常、形成并检验假设，等等。花些时间去理解你的数据集会在你试图从模型中挤出每一分性能时帮助你很多。你可能会理解必须创建什么样的特征，应该应用什么样的建模技术，等等。

在本章中，我们将介绍几种适合时间序列数据集的可视化技术。

笔记本提示：

要跟随完整的可视化时间序列代码，请使用Chapter03文件夹中的01-Visualizing_Time_Series.ipynb笔记本。

折线图

这是最基本和常用的时间序列可视化方式。我们只需将时间放在x轴，将时间序列的值放在y轴上。如果我们绘制数据集中的一个家庭，看看它是什么样子的：

当你拥有一个变化较大的长时间序列时（就像我们有的那样），折线图可能会显得有些混乱。为了从趋势和变化的宏观视角来看时间序列的整体情况，一个选项是绘制时间序列的平滑版本。让我们看看时间序列的滚动月平均值是什么样子的：

现在我们可以更清楚地看到宏观模式。季节性非常明显——该系列在冬季达到峰值，在夏季达到低谷。如果仔细想想，这就有道理了。我们说的是伦敦，冬季的能源消耗会更高，因为低温和随之而来的供暖系统使用。

例如，对于热带地区的家庭，模式可能是相反的，峰值出现在夏季，当空调开始使用时。

折线图的另一个用途是将两个或更多的时间序列一起可视化，并调查它们之间的相关性。在我们的例子中，试着将温度与能源消耗一起绘制，看看我们关于温度影响能源消耗的假设是否成立：

在这里，我们可以看到能源消耗与温度之间的年度分辨率的明显负相关关系。冬季在宏观尺度上显示出较高的能源消耗。我们还可以看到与温度 loosely 相关的每日模式，但这可能是由于其他因素，如人们下班后回家等原因。

还有一些其他的可视化方法更适合突出时间序列中的季节性。让我们来看看。

季节性图（Seasonal Plots）

季节性图与折线图非常相似，但关键区别在于，x轴表示“季节”，y轴表示时间序列值，不同的季节性周期以不同的颜色或线型表示。例如，按月分辨率的年度季节性可以通过将月份放在x轴上，并用不同的颜色表示不同的年份来绘制。

让我们看看这对我们提到的家庭是怎样的。这里，我们绘制了多年来的平均每月能源消耗：

我们可以立即看到这种可视化的吸引力，因为它让我们轻松地可视化季节性模式。我们可以看到，夏季月份的能源消耗下降，而且我们还可以看到这种模式在多个年份中一致发生。在我们拥有数据的这两年中，我们可以看到2013年10月的行为略有偏离2012年。也许还有其他因素可以帮助我们解释这个差异——温度呢？

我们还可以将季节性图与另一个感兴趣的变量一起绘制，比如温度：

注意到十月了吗？在2013年10月，温度持续较高，持续了一个月，这就是为什么能源消耗模式与去年略有不同的原因。

我们还可以在其他分辨率下绘制这些图表，例如每小时的季节性。我们需要做的就是计算每小时和每月每天的平均消耗量，然后将它们绘制成图，x轴是小时，不同的日期用不同的颜色表示（图3.8（上））。但当季节性周期过多时，绘图会变得混乱。季节性图的替代方案是季节性箱线图。

季节性箱线图

与将不同的季节性周期用不同颜色或线型绘制不同，我们可以将它们表示为箱线图（图3.8（下））。这可以立即清除图中的杂乱。通过这种表示方式，我们还能获得额外的好处——它使我们能够理解不同季节性周期之间的变异性：

在这里，我们可以看到，在这种分辨率下的季节性图表过于混乱，难以看出模式和季节性周期之间的变异性。然而，季节性箱线图则更具信息量。箱线图中的水平线告诉我们中位数，箱体表示四分位数范围（IQR），而标记的点则是异常值。通过观察中位数，我们可以看到能源消耗的峰值出现在上午9点之后。但上午9点之后的变异性也更高。如果你为每一周单独绘制箱线图，你会看到周日的模式略有不同（更多的可视化内容在配套的笔记本中）。

然而，还有一种可视化方法可以让你从两个维度检查这些模式。

日历热力图

与为每一天的每一周绘制单独的箱线图或折线图不同，如果我们能够将这些信息压缩成一个单一的图表，那将会更有用。这就是日历热力图的作用。热力图可视化使用颜色渐变来表示数据值，不同的颜色表示不同的强度或频率。日历热力图在一个矩形块中使用彩色单元格来表示信息。在矩形的两侧，我们可以找到两个不同的时间粒度，如月份和年份。在每个交点处，单元格的颜色与该交点处时间序列的值相关。

让我们看看一个日历热力图，展示了不同工作日的每小时平均能源消耗（请参考彩色图像文件：packt.link/gbp/9781835…）：

从右侧的颜色尺度中，我们可以知道，较浅的颜色表示较高的值。我们可以看到，周一到周六的峰值类似——即早上一次，晚上一次。然而，周日的模式略有不同，全天的能源消耗都较高。

到目前为止，我们回顾了很多可以揭示季节性的可视化方法。现在，让我们看看用于检查自相关的可视化方法。

自相关图

如果相关性表示两个变量之间线性关系的强度和方向，那么自相关是时间序列中连续时期之间的值的相关性。大多数时间序列在很大程度上依赖于前一个时期的值，这在许多我们将要看到的预测模型中也是一个关键组成部分。

例如，ARIMA（我们将在第4章《设置强基线预测》中简要介绍）就是基于自相关构建的。因此，直观地了解前期时间步的依赖性有助于更好地理解模型。

这时，自相关图就非常有用。在这些图中，x轴表示不同的滞后（t-1，t-2，t-3，依此类推），y轴表示t和不同滞后之间的相关性。除了自相关外，我们还可以查看部分自相关，它与自相关非常相似，但有一个关键区别：部分自相关去除了任何可能存在的间接相关性后，再呈现相关性。让我们通过一个例子来理解这一点。如果t是当前的时间步，假设t-1与t高度相关。那么，按照这个逻辑，t-2将与t-1高度相关，因此，由于这种相关性，t和t-2之间的自相关性会很高。然而，部分自相关会纠正这一点，并提取出仅能归因于t-2和t的相关性。

我们需要记住的一点是，自相关和部分自相关分析在时间序列平稳时效果最好（我们将在第6章《时间序列预测的特征工程》中详细讨论平稳性）。

最佳实践：

有很多方法可以使序列平稳，但一种快速且简单的方法是使用季节性分解，并仅选择残差。残差应该不含有趋势和季节性，而这些通常是导致时间序列非平稳的主要因素。但正如我们稍后将在本书中看到的那样，这并不是使序列在最真实的意义上平稳的万无一失的方法。

现在，让我们看看这些图表在将数据集中的家庭数据（平稳化后）应用时的效果：

在这里，我们可以看到，第一滞后（t-1）具有最大的影响力，而且在部分自相关图中，它的影响力迅速下降到接近零。这意味着一天的能源消耗与前一天的能源消耗高度相关。

如果你以前见过这样的图表，你可能会看到图表上方有一个包络线，显示置信区间，用于作为选择显著自相关性的指导。虽然这是一种好的经验法则，但这里没有包含置信区间，因为我不希望你把它当作一个规则。置信区间的相关性取决于一些假设（如正态性等），这些假设并不是总能满足，特别是在实际应用中。

有了这些，我们已经了解了时间序列的不同组成部分，并学会了如何可视化其中的一些。现在，让我们看看如何将时间序列分解为它的组成部分。

分解时间序列

季节性分解是将时间序列分解为其组成部分的过程——通常是趋势、季节性和残差。分解时间序列的基本方法如下：

1. 去趋势（Detrending） ：在这里，我们估计趋势成分（即时间序列的平滑变化），并将其从时间序列中去除，从而得到去趋势的时间序列。
2. 去季节性（Deseasonalizing） ：在这里，我们从去趋势后的时间序列中估计季节性成分。去除季节性成分后，剩下的就是残差。

让我们详细讨论这些过程。

去趋势（Detrending）

去趋势可以通过几种不同的方法来实现。两种常见的去趋势方法是使用移动平均和局部估计散点平滑回归（LOESS回归）。

移动平均（Moving Averages）

估计趋势的最简单方法之一是使用时间序列上的移动平均。它可以看作是一个沿时间序列滑动的窗口，在每一步中，窗口内所有值的平均值会被记录下来。这个移动平均是一个平滑后的时间序列，有助于我们估计时间序列中的缓慢变化，也就是趋势。缺点是这种方法非常容易受到噪声的影响。即使使用这种方法平滑了时间序列，提取出的趋势仍然会有噪声，无法做到完全平滑。理想情况下，噪声应该集中在残差部分，而不是趋势部分（参见图3.13中的趋势线）。

LOESS

LOESS算法，也叫局部加权多项式回归，是由Bill Cleveland在70年代到90年代开发的。它是一种非参数方法，用于将平滑曲线拟合到噪声信号上。我们使用一个在时间序列间移动的序列变量作为自变量，时间序列信号作为因变量。

对于序列变量中的每一个值，算法使用离该点最近的一部分数据点，并仅使用这些点通过加权回归来估计平滑的趋势。加权回归中的权重是与当前点最近的数据点。离该点最近的点会被赋予最高的权重，随着距离的增加，权重逐渐减小。这为我们提供了一个非常有效的工具，用于建模时间序列中的平滑变化（趋势）（参见图3.14中的趋势线）。

去季节性（Deseasonalizing）

季节性成分也可以通过几种不同的方法来估计。最常见的两种方法是使用周期调整平均值或傅里叶级数。

周期调整平均值（Period Adjusted Averages）

这是一种相当简单的技术，其中我们通过计算预期周期中每个周期的季节性指数，方法是取所有周期中该周期所有值的平均值。为了使其更清晰，我们来看一个月度时间序列，假设该时间序列有年度季节性。因此，起伏模式将在12个月内完成一个完整的周期，或者说季节性周期是12。换句话说，时间序列中的每12个数据点有相似的季节性成分。所以，我们取所有1月值的平均数作为1月的周期调整平均值。以同样的方式，我们计算12个月的周期平均值。在这个过程中，我们最终会得到12个周期平均值，并且我们还可以计算一个平均周期平均值。现在，我们可以通过以下方式将这些周期平均值转换为一个指数：要么从每个周期平均值中减去所有周期平均值的平均值（加法模型），要么从每个周期平均值中除以所有周期平均值的平均值（乘法模型）。

傅里叶级数（Fourier Series）

在18世纪末，数学家和物理学家Joseph Fourier在研究热流时意识到了一些深刻的事情——任何周期函数都可以分解为一组简单的正弦和余弦波。让我们稍作停顿，思考一下。任何周期函数，无论它的形状如何，是否有曲线，或者它如何围绕轴线大幅振荡，都可以分解成一系列正弦和余弦波。

附加信息：

对于更有数学背景的人来说，原始理论提出将任何周期函数分解为指数的积分。使用欧拉公式，我们可以将它们视为正弦和余弦波的总和。如果你想深入研究并探索相关概念（如傅里叶变换），进一步阅读部分包含了一些资源。

正是利用这一特性，我们提取时间序列中的季节性成分，因为季节性是一个周期性函数，任何周期性函数都可以用正弦和余弦波的组合来近似。傅里叶级数的正弦-余弦形式如下：

在这里，SN 是信号 S 的 N 项近似值。理论上，当 N 无穷大时，得到的近似值等于原始信号。P 是周期的最大长度。

我们可以使用这个傅里叶级数，或者傅里叶级数中的几个项，来建模我们的季节性。在我们的应用中，P 是我们尝试建模的周期的最大长度。例如，对于每年季节性的月度数据，周期的最大长度（P）是 12。x 将是一个从 1 到 P 增长的序列变量。在这个例子中，x 将是 1, 2, 3, … 12。现在，利用这些项，我们要做的就是找到 an 和 bn，方法是对信号进行回归。

我们已经看到，通过正确的傅里叶项组合，我们可以复制任何信号。但问题是，我们应该这么做吗？我们想从数据中学习的是一个通用的季节性剖面，这个剖面不仅能很好地适应已有数据，还能适应未知数据。因此，我们使用 N 作为一个超参数，从数据中提取我们想要的复杂信号。

这是一个很好的时机，复习一下你的三角学并记住正弦和余弦波的形状。第一个傅里叶项（n=1）是你熟悉的正弦和余弦波，它们在最大周期长度（P）内完成一个完整的周期。随着 n 的增加，我们得到的正弦和余弦波将在最大周期长度（P）内完成多个周期。这可以在以下图中看到：

正弦波和余弦波是互补的，如下图所示：

现在，让我们看看如何在实际中使用这一方法。

实现

笔记本提示：

要跟随完整的时间序列分解代码，请使用Chapter03文件夹中的02-Decomposing_Time_Series.ipynb笔记本。

我们将在以下小节中介绍四种实现方法。

使用 statsmodels 的 seasonal_decompose

statsmodels.tsa.seasonal 中有一个名为 seasonal_decompose 的函数。该函数使用移动平均法来处理趋势成分，并使用周期调整的平均值来处理季节性成分。它支持加法和乘法两种分解模式。然而，它不支持缺失值。让我们看看如何使用它：

# 不支持缺失值，因此使用插补后的时间序列
res = seasonal_decompose(ts, period=7*48, model="additive", extrapolate_trend="freq")

需要记住的几个关键参数如下：

• period：是你预计模式重复的季节性周期。
• model：接受“additive”或“multiplicative”作为参数，用于确定分解的类型。
• filt：传入一个数组，用作移动平均中的权重（具体来说是卷积）。它还可以用于定义我们需要计算移动平均的窗口。我们可以增加它，以在某种程度上平滑趋势成分。
• extrapolate_trend：是一个可以用来将趋势成分延伸到两侧的参数，从而避免应用移动平均滤波器时产生的缺失值。
• two_sided：是一个允许我们定义如何计算移动平均的参数。如果为True（默认值），则移动平均使用过去和未来的值计算，因为移动平均的窗口是居中的。如果为False，它只使用过去的值来计算移动平均。

让我们看看我们在数据集中的一个时间序列分解得如何。我们使用 period=7*48 来捕捉工作日每小时的模式，并使用 filt=np.repeat(1/(30*48), 30*48) 来对30天进行移动平均，权重均匀分配：

我们无法看到季节性模式，因为它在图表的整体尺度下太小。相关的笔记本中有放大的图表，帮助你理解季节性模式。即使使用较大的窗口（例如20天）进行平滑，趋势中仍然有一些噪声。我们或许可以通过增加窗口来减少噪声，但正如我们将看到的，还有一个更好的替代方案。

使用 LOESS 进行季节性和趋势分解（STL）

正如我们之前看到的，LOESS 更适合用于趋势估计。STL 是一种实现，它使用 LOESS 进行趋势估计，并使用周期平均值来估计季节性。尽管 statsmodels 也有一个实现，但我们重新实现了它，以提高性能和灵活性。

这个实现可以在本书的 GitHub 仓库中找到，路径为 src.decomposition.seasonal.py。它期望一个带有 datetime 索引的 pandas DataFrame 或 series 作为输入。让我们看看如何使用它：

stl = STL(seasonality_period=7*48, model = "additive")
res_new = stl.fit(ts_df.energy_consumption)

这里需要注意的关键参数如下：

• seasonality_period：是你预计模式重复的季节性周期。
• model：接受“additive”或“multiplicative”作为参数，用于确定分解的类型。
• lo_frac：表示用于拟合 LOESS 回归的数据比例。
• lo_delta：是一个比例距离，用于指定我们在哪些情况下使用线性插值而不是加权回归。使用非零的 lo_delta 会显著减少计算时间。

让我们看看这种分解的效果。这里，我们使用 seasonality_period=7*48 来捕捉工作日每小时的模式：

让我们也看一下仅针对一个月的分解，以更清晰地看到提取的季节性模式：

趋势现在已经足够平滑，季节性也已经被提取出来。这里，我们可以清楚地看到每小时的峰值和谷值，以及周末的较高峰值。然而，由于我们依赖于平均值来推导季节性，它也容易受到异常值的影响。时间序列中一些非常高或非常低的值会扭曲基于周期平均值提取的季节性模式。这个方法的另一个缺点是，当数据的分辨率与期望的季节性周期之间的差异较大时，提取的季节性的“好度”会受到影响。例如，当在日常或子日数据上提取年度季节性时，这会使提取的季节性非常嘈杂。如果你拥有的季节性周期少于两个周期，这个技术也不起作用——例如，如果我们想提取年度季节性，但我们只有不到2年的数据。

傅里叶分解

我们可以在 src.decomposition.seasonal.py 中找到用于通过傅里叶项分解时间序列的 Python 实现。它使用 LOESS 进行趋势检测，使用傅里叶项提取季节性。我们可以通过两种方式使用它。首先，我们可以将 seasonality_period 指定为 pandas datetime 属性之一（如小时和星期几）：

stl = FourierDecomposition(seasonality_period="hour", model = "additive", n_fourier_terms=5)
res_new = stl.fit(pd.Series(ts.squeeze(), index=ts_df.index))

另外，我们可以创建一个自定义的季节性数组，它的长度与时间序列相同，并且具有季节性的序数表示。如果是日数据的年度季节性，数组的最小值为1，最大值为365，每天增加一次：

# 创建自定义的季节性项
ts_df["dayofweek"] = ts_df.index.dayofweek
ts_df["hour"] = ts_df.index.hour
# 创建一个排序后的唯一组合数据框
map_df = ts_df[["dayofweek","hour"]].drop_duplicates().sort_values(["dayofweek", "hour"])
# 分配一个序数变量来捕捉顺序
map_df["map"] = np.arange(1, len(map_df)+1)
# 将序数映射回原始数据框并获取季节性数组
seasonality = ts_df.merge(map_df, on=["dayofweek","hour"], how='left', validate="many_to_one")['map']
stl = FourierDecomposition(model = "additive", n_fourier_terms=50)
res_new = stl.fit(pd.Series(ts, index=ts_df.index), seasonality=seasonality)

此过程涉及的关键参数如下：

• seasonality_period：从 datetime 索引中提取的季节性。可以使用 pandas datetime 属性（如 week_of_day 和 month）来指定最突出季节性。如果将其设置为 None，则在调用 fit 时需要提供季节性数组。
• model：接受 "additive" 或 "multiplicative" 作为参数，用于确定分解的类型。
• n_fourier_terms：确定用于提取季节性的傅里叶项的数量。增加这个参数的值，可以提取更复杂的季节性。
• lo_frac：用于拟合 LOESS 回归的数据比例。
• lo_delta：是一个比例距离，用于指定在哪些情况下使用线性插值而不是加权回归。使用非零的 lo_delta 会显著减少计算时间。

让我们看看使用 FourierDecomposition 进行分解后的放大图：

趋势将与 STL 的结果相同，因为我们在这里也使用了 LOESS。季节性剖面可能略有不同，并且对异常值具有更强的鲁棒性，因为我们通过傅里叶项在信号上进行正则化回归。另一个优点是，我们已经解耦了数据的分辨率和预期的季节性。现在，在子日数据上提取年度季节性不再像使用周期平均值时那样具有挑战性。

到目前为止，我们只看到了提取每个序列一个季节性的技术；大多数情况下，我们提取的是主要季节性。那么，当我们有多个季节性模式时，应该怎么做呢？

使用 LOESS 进行多季节性分解（MSTL）

具有高频数据（如日数据或小时数据）的时间序列容易表现出多个季节性模式。例如，可能存在小时季节性模式、周季节性模式和年季节性模式。但如果我们只提取主导模式并将其余部分留给残差，我们就没有做到真正的分解。Kasun Bandara 等人提出了 STL 分解的一个扩展，称为 MSTL，用于多季节性分解，并且在 R 生态系统中有相应的实现。在 Python 中，类似的实现可以在 src.decomposition.seasonal.py 中找到。除了 MSTL，该实现还使用傅里叶项提取多个季节性。

参考文献检查：

Kasun Bandara 等人的研究论文已在参考文献部分作为参考文献1引用。

让我们看看如何使用这个方法的例子：

stl = MultiSeasonalDecomposition(seasonal_model="fourier", seasonality_periods=["day_of_year", "day_of_week", "hour"], model = "additive", n_fourier_terms=10)
res_new = stl.fit(pd.Series(ts, index=ts_df.index))

这里涉及的关键参数如下：

• seasonality_periods：预期季节性的列表。对于 STL，它是一个季节性周期的列表，而对于 FourierDecomposition，它是一个字符串列表，表示 pandas datetime 属性。
• seasonality_model：接受 "fourier" 或 "averages" 作为参数，用于确定季节性分解的类型。
• model：接受 "additive" 或 "multiplicative" 作为参数，用于确定分解的类型。
• n_fourier_terms：确定用于提取季节性的傅里叶项的数量。增加该参数值，提取的季节性将更加复杂。
• lo_frac：用于拟合 LOESS 回归的数据比例。
• lo_delta：是一个比例距离，用于指定在哪些情况下使用线性插值而不是加权回归。使用非零的 lo_delta 会显著减少计算时间。

让我们看看使用傅里叶分解时分解结果的样子：

在这里，我们可以看到 day_of_week 的季节性已经被提取出来。为了查看 day_of_week 和 hour 的季节性成分，我们需要稍微放大一些：

在这里，我们可以观察到小时季节性已被很好地提取出来，并且它也成功地隔离了 day_of_week 的季节性成分，该成分在周末达到峰值。day_of_week 季节性成分的离散步骤性质是因为数据的频率是半小时一次，对于 48 个数据点，day_of_week 会保持不变。

我们总结了我们所覆盖的四种技术，并整理在下面的表格中：

实现方法	趋势	季节性	是否支持多季节性？	是否支持缺失值？
移动平均	LOESS	周期调整平均值	否	否
STL	LOESS	周期调整平均值	否	是
傅里叶分解	LOESS	傅里叶项	否	否
多季节性分解	LOESS	周期调整平均值 / 傅里叶项	是	否

表 3.1：不同的季节性分解技术

MSTL 也已经在 statsmodels 中实现，相关的笔记本中有使用它的代码。statsmodels 和本书代码中实现的关键区别在于，本书中的实现还提供了基于傅里叶级数的分解选项。

现在，让我们开始理解和分析一个时间序列数据集。

检测和处理异常值

如其名所示，异常值 是指与其余观测值之间距离异常的观测值。如果我们将数据生成过程（DGP）看作是一个生成时间序列的随机过程，那么异常值就是那些从 DGP 中生成的概率最小的点。这可能是由多种原因引起的，包括故障的测量设备、错误的数据输入、黑天鹅事件等。能够检测和处理这些异常值可能有助于预测模型更好地理解数据。

异常值/异常检测本身就是时间序列中的一个专门领域，但在本书中，我们将局限于一些更简单的技术来识别和处理异常值。这是因为我们的主要目标不是检测异常值，而是清理数据，使我们的预测模型表现得更好。如果你想深入学习异常检测，可以查阅进一步阅读部分，了解一些入门资源。

现在，让我们看看一些识别异常值的技术。

笔记本提示：

要跟随完整的检测异常值代码，请使用 Chapter03 文件夹中的 03-Outlier_Detection.ipynb 笔记本。

标准差（Standard Deviation）

这是一个几乎每个从事数据工作的人都听说过的经验法则——如果 $μ$ 是时间序列的均值，$σ$ 是标准差，那么任何落在$μ±3σ$之外的观测值就是异常值。其背后的理论深深根植于统计学。如果我们假设时间序列的值服从正态分布（正态分布是一种具有非常理想性质的对称分布），根据概率理论，我们可以推导出正态分布下68%的区域位于均值的一个标准差内，约95%的区域位于两个标准差内，约99%的区域位于三个标准差内。因此，当我们使用三倍标准差作为界限时，我们的意思是，如果某个观测值属于该概率分布的概率低于1%，那么它就是一个异常值。

稍微偏向一些实际问题时，这个三倍标准差的截断值并不是不可动摇的。我们需要尝试不同的倍数值，并通过主观评估结果来确定合适的倍数。倍数越大，异常值越少。

对于高度季节性的数据，直接将规则应用于原始时间序列的方法并不好。对于这种情况，我们必须先使用我们之前讨论的任意方法去季节化数据，然后将异常值检测应用于残差。如果不这样做，我们可能会将季节性峰值错误地标记为异常值，这显然不是我们想要的。

这里的另一个关键假设是正态分布。然而，实际上，我们遇到的许多时间序列可能不是正态分布，因此这个规则会迅速失去其理论保证。

四分位距（IQR）

另一种非常相似的技术是使用四分位距（IQR）代替标准差来定义超出边界的异常值。分位数将所有数据按顺序排列，然后将其分成相等的部分，使得每部分有相同数量的元素。四分位数做的是相同的事，但特别是将数据分成四等份。IQR 是第 3 四分位数（或第 75 百分位数或 0.75 分位数）与第 1 四分位数（或第 25 百分位数或 0.25 分位数）之间的差值。上下界定义如下：

• 上界 =$Q3+n×IQR$
• 下界 =$Q1−n×IQR$

其中，IQR =$Q3−Q1$，n 是 IQR 的倍数，决定了可接受区域的宽度。

对于我们观察到有大量异常值和剧烈波动的数据集，这比标准差方法更具鲁棒性。这是因为标准差和均值会受到数据集中的个别点的影响。如果在之前的方法中使用的是标准差，那么在这里是 1.5 倍的 IQR。这也与正态分布假设相关，而 1.5 倍的 IQR 大约等于 ~3标准差（精确值）。在这里应用规则之前去季节化的建议也适用。它适用于我们将要看到的所有技术。

隔离森林（Isolation Forest）

隔离森林是一种基于决策树的无监督异常检测算法。典型的异常检测算法是通过建模正常数据点并将不符合正常模式的点标记为异常点。然而，隔离森林采用了不同的方式，直接建模异常点。它通过随机分割特征空间来创建决策树森林。该技术假设异常点位于外部边缘，更容易进入树的叶节点。因此，你可以在短分支中找到异常点，而正常点由于彼此靠近，需要更长的分支。任何点的“异常分数”是通过计算在到达该点之前需要遍历的树的深度来确定的。scikit-learn 在 sklearn.ensemble.IsolationForest 中有该算法的实现。除了决策树的标准参数外，关键参数是 contamination。默认值为 auto，但可以设置为 0 和 0.5 之间的任何值。该参数指定你期望数据集中有多少百分比是异常的。

但我们必须记住的是，隔离森林完全不考虑时间，仅仅突出那些偏离正常模式的值。

极端学生化偏差（ESD）和季节性极端学生化偏差（S-ESD）

这种基于统计学的方法比基本方法更复杂，但仍然使用正态性假设。它基于另一个统计测试，称为 Grubbs 测试，用于在正态分布的数据集中找出一个单独的异常值。ESD 通过迭代使用 Grubbs 测试，每次识别并去除一个异常值。它还根据剩余点数调整临界值。要更详细地了解该测试，请参阅进一步阅读部分，在其中我们提供了有关 ESD 和 S-ESD 的一些资源。2017 年，Twitter 研究的 Hochenbaum 等人提出使用去季节化的广义 ESD 方法作为时间序列异常值检测方法。

我们已经为我们的用例调整了现有算法的实现，并且它可以在本书的 GitHub 仓库中找到。与其他方法不同，S-ESD 只需输入期望的异常值上限，然后独立识别异常值。例如，我们将上限设置为 800，算法在我们使用的数据中识别出了大约 400 个异常值。

参考文献检查：

Hochenbaum 等人的研究论文已在参考文献部分作为参考文献 2 引用。

让我们看看我们使用所有已讨论的技术检测到的异常值：

我们已经学会了如何检测异常值，接下来让我们讨论如何处理这些异常值并清理数据集。

处理异常值

我们必须首先回答的一个问题是，是否应该修正我们识别出来的异常值。自动识别异常值的统计测试应该经过另一轮人工验证。如果我们盲目地“处理”异常值，我们可能会切断一个有价值的模式，而这个模式可能对我们预测时间序列有帮助。如果你只预测少量时间序列，那么查看异常值并通过研究这些异常值的原因将其与现实对接是有意义的。

但当你有成千上万的时间序列时，一个人无法检查所有异常值，因此我们必须依赖自动化技术。一种常见做法是用一个启发式方法来替换异常值，比如使用最大值、最小值和第75百分位数。更好的方法是将异常值视为缺失数据，并使用我们之前讨论的任何方法来填补这些异常值。

需要记住的一点是，异常值修正并不是预测中的必要步骤，尤其是在使用机器学习或深度学习等现代方法时。是否进行异常值修正是我们需要通过实验来决定的。

恭喜你！

这是一个相当繁忙的章节，涉及了很多概念和代码，因此恭喜你完成了它。如果需要，随时可以回去复习几个主题。

总结

在本章中，我们学习了时间序列的关键组成部分，并熟悉了如趋势和季节性等术语。我们还回顾了几种适用于时间序列的可视化技术，这些技术在进行探索性数据分析（EDA）时会非常有用。接着，我们了解了几种将时间序列分解为其组成部分的技术，并学习了数据中异常值的检测方法。最后，我们学习了如何处理已识别的异常值。现在，你已经为开始预测时间序列做好了准备，我们将在下一章开始这一部分内容。

参考文献

以下是本章的参考文献：

• Kasun Bandara, Rob J Hyndman, Christoph Bergmeir. (2021). MSTL: A Seasonal-Trend Decomposition Algorithm for Time Series with Multiple Seasonal Patterns. arXiv:2107.13462 [stat.AP]. arxiv.org/abs/2107.13…
• Hochenbaum, J., Vallis, O., & Kejariwal, A. (2017). Automatic Anomaly Detection in the Cloud Via Statistical Learning. ArXiv, abs/1704.07706. arxiv.org/abs/1704.07…

进一步阅读

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

• 傅里叶级数：www.setzeus.com/public-blog…
• Khan Academy的傅里叶级数：www.youtube.com/watch?v=UKH…
• 傅里叶变换：betterexplained.com/articles/an…
• Ane Blázquez-García, Angel Conde, Usue Mori, 和 Jose A. Lozano. (2021). A Review on Outlier/Anomaly Detection in Time Series Data. arXiv:2002.04236. arxiv.org/abs/2002.04…
• Braei, M., & Wagner, S. (2020). Anomaly Detection in Univariate Time-series: A Survey on the State-of-the-Art. ArXiv, abs/2004.00433. arxiv.org/abs/2004.00…
• 广义ESD测试（异常值检测） ：www.itl.nist.gov/div898/hand…