时间序列数据的分析与可视化

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

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

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

技术要求

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

你需要从Chapter02文件夹运行02-Preprocessing_London_Smart_Meter_Dataset.ipynb笔记本。

本章的代码可在以下链接找到：https://github.com/PacktPublishing/Modern-Time-Series-Forecasting-with-Python-2E/tree/main/notebooks/Chapter03。

时间序列的组成部分

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

趋势
季节性
循环
不规则成分

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

趋势成分

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

图3.1：特斯拉收入（百万美元）

从上图可以看出，特斯拉的收入呈增长趋势。趋势不必是线性的，也可以是非线性的。

季节性成分

当一个时间序列呈现出规律性、重复性的上下波动时，我们称之为季节性。例如，零售额通常在节假日期间飙升，尤其是在西方国家的圣诞节期间。类似地，热带地区夏季和寒冷国家冬季的用电量达到峰值。在这些例子中，你可以看到每年重复的特定上下波动模式。另一个例子是太阳黑子，如下图所示：

图3.2：1749年至2017年的太阳黑子数量

如你所见，太阳黑子每11年达到一次峰值。

周期性分量

周期性分量常与季节性混淆，但由于一个非常微妙的差异而有所区别。与季节性类似，周期性分量也围绕趋势线呈现类似的上下波动模式，但该周期运行的时间并非固定，而是在一个大致时间框架内略有变化。一个很好的例子是经济衰退，它大约每10年发生一次。然而，这并非如时钟般精准；有时可能少于或多于10年。

不规则分量

该分量是从时间序列中去除趋势、季节性和周期性后剩余的部分。传统上，该分量被认为是不可预测的，也被称为残差、误差项或噪声项。在常见的经典统计模型中，任何“模型”的目标都是捕捉所有其他分量，使得未被捕捉的部分仅剩下不规则分量。

在现代机器学习中，我们并不认为该分量完全不可预测。我们试图通过使用外生变量来捕捉该分量的部分。例如，零售额的不规则分量可能由其开展的促销活动来解释。当我们拥有这些额外信息时，“不可预测”的分量开始变得可预测。但无论您向模型添加多少额外变量，总会留下一些分量，即不可约误差。这是时间序列中无论模型多强或添加多少额外信息都无法解释的部分。

现在我们已经了解了时间序列的不同分量，接下来看看如何对其进行可视化。

可视化时间序列数据

在第2章《获取与处理时间序列数据》中，我们学习了如何准备数据模型，作为分析新数据集的第一步。如果将准备数据模型比作接近你喜欢的人并初次接触，那么EDA就像与那个人约会。此时，你拥有数据集，正在试图了解它，弄清楚它如何运作、喜好什么等等。

EDA通常使用可视化技术来发现模式、发现异常、形成和检验假设等。花一些时间理解你的数据集，在试图从模型中榨取每一分性能时，会大有帮助。你可能会理解必须创建哪些特征，应该应用哪种建模技术等。

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

笔记本提示：

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

折线图

这是用于理解时间序列的最基本且最常见的可视化。我们只需将时间绘制在 x 轴上，时间序列值绘制在 y 轴上。让我们看看如果绘制数据集中一个住户的情况会是什么样子：

图3.3：住户MAC000193的折线图

当你有像我们这样具有高变异性的长时间序列时，折线图可能会变得有点混乱。为了从趋势和运动角度宏观查看时间序列，一种选择是绘制时间序列的平滑版本。让我们看看时间序列的滚动月平均值是什么样的：

图3.4：住户MAC000193的滚动月平均能耗

现在我们可以更清楚地看到宏观模式。季节性很明显——序列在冬季达到峰值，在夏季达到低谷。如果你认真思考，这很有道理。我们说的是伦敦，由于温度较低和随后的供暖系统使用，冬季的能耗会更高。

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

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

图3.5：温度与能源消耗（底部为缩放图）

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

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

季节性图

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

让我们看看这对我们讨论的家庭来说是什么样的。在这里，我们绘制了多年间的平均月能源消耗：

图3.6：月度分辨率的季节性图

我们可以立即看到这种可视化的吸引力，因为它让我们很容易看到季节性模式。我们可以看到夏季月份的消耗量下降，并且这种现象在多年间一致出现。在我们有数据的两年中，可以看到2013年10月的行为与2012年略有偏差。也许还有其他因素可以帮助解释这种差异——温度呢？

我们还可以绘制另一个感兴趣变量（如温度）的季节性图：

图3.7：月度分辨率的季节性图（能耗与温度的关系）

注意到十月了吗？2013年10月，气温持续偏暖多了一个月，因此能耗模式与去年略有不同。

我们也可以在其他分辨率下绘制这类图，比如每小时季节性。我们只需计算每个月每小时的平均能耗，并以小时为x轴，月份的不同日期用不同颜色绘制（图3.8（上））。但当要绘制的季节性周期过多时，会增加视觉混乱。季节性图的替代方案是季节性箱线图。

季节性箱线图

与其用不同颜色或线型绘制不同的季节性周期，不如将其表示为箱线图（图3.8（下））。这能立即消除图中的杂乱。这种表示的额外好处是让我们了解各季节性周期之间的变异性：

图3.8：小时分辨率的季节性图（上）和季节性箱线图（下）

这里可以看出，该分辨率的季节性图过于杂乱，难以看清模式和季节性周期间的变化。而季节性箱线图信息量更大。箱中的水平线表示中位数，箱体表示四分位距（IQR），标记的点为异常值。通过观察中位数，我们可以看到高峰能耗从上午9点开始。但变异性也从上午9点起更高。例如，若为每周单独绘制箱线图，您会发现星期日的模式略有不同（更多可视化内容见配套笔记本）。

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

日历热力图

比起为每天每星期分别绘制箱线图或折线图，如果能将信息压缩到一张图中会更有用。这时候就需要用到 日历热力图。热力图可视化使用颜色梯度来表示数据值，不同颜色代表不同的强度或频率。日历热力图在矩形块中使用彩色单元格来表示信息。在矩形的两侧，可以找到两个不同的时间粒度，例如月份和年份。在每个交叉点，单元格的颜色根据该交叉点的时间序列值而定。

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

图3.9：能源消耗的日历热力图

从右侧的颜色刻度可知，浅色代表较高值。我们可以看到周一至周六有相似的峰值——即一次在早晨，一次在傍晚。然而，周日的模式略有不同，全天的能耗较高。

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

自相关图

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

诸如 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章 时间序列预测的特征工程 中详细讨论平稳性）。

最佳实践：

有很多方法可以使序列平稳，但一种快速粗糙的方法是使用季节分解并只取残差。它应当没有趋势和季节性，这两者是时间序列非平稳的主要驱动因素。但正如我们将在本书后面看到的，这并不是一种真正意义上使序列平稳的万无一失的方法。

现在，让我们看看数据集中的家庭在（使其平稳后）这些图是什么样的：

图3.10：自相关和偏自相关图

在这里，我们可以看到第一个滞后（t-1）影响最大，并且其影响在偏自相关图中迅速下降到接近零。这意味着某天的能耗与前一天的能耗高度相关。

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

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

分解时间序列

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

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

让我们详细讨论它们。

去趋势

去趋势可以通过几种不同的方式完成。两种流行的方法是使用移动平均和局部估计散点图平滑（LOESS）回归。

移动平均

估计趋势的最简单方法之一是使用沿时间序列的移动平均。它可以看作是一个窗口，沿着时间序列逐步移动，在每个步骤中，记录窗口中所有值的平均值。这个移动平均是一个平滑后的时间序列，帮助我们估计时间序列中的缓慢变化，即趋势。缺点是这种方法噪声很大。即使使用这种方法平滑时间序列后，提取的趋势也不会平滑；它会有噪声。理想情况下，噪声应该存在于残差中，而不是趋势中（参见图3.13中显示的趋势线）。

LOESS

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

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

去季节化

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

周期调整平均值

这是一种相当简单的技术，我们通过取所有周期中每个周期的平均值来计算预期周期中每个时期的季节指数。为了说明清楚，让我们看一个月度时间序列，其中我们预期有年度季节性。因此，上下波动模式将在12个月内完成一个完整周期，即季节周期为12。换句话说，时间序列中每12个点具有相似的季节成分。所以，我们取所有1月份值的平均值作为1月份的周期调整平均值。同样地，我们计算所有12个月的周期平均值。最后，我们得到12个周期平均值，我们还可以计算一个平均周期平均值。现在，我们可以通过从每个周期平均值中减去所有周期平均值的平均值（加法）或除以所有周期平均值的平均值（乘法）来将这些周期平均值转换为指数。

傅里叶级数

18世纪末，数学家兼物理学家约瑟夫·傅里叶在研究热流时，意识到一个深刻的事实——任何周期函数都可以分解为一系列简单的正弦波和余弦波。让我们仔细思考一下。任何周期函数，无论其形状、曲线、缺失情况，或围绕轴线的振荡多么剧烈，都可以分解为一系列正弦波和余弦波。

附加信息：

对于数学爱好者和希望深入研究的读者，原始理论提出将任何周期函数分解为指数函数的积分。利用欧拉恒等式，$e^{iy} = \cos(y) + i \cdot \sin(y)$，我们可以将其视为正弦波和余弦波的求和。进一步阅读部分包含一些资源，如果你希望深入了解傅里叶变换等相关概念。

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

$$S_{N(x)} = \frac{a_0}{2} + \sum_{n=1}^{N} \left( a_n \cdot \cos\left(\frac{2\pi}{P} \cdot n \cdot x\right) + b_n \cdot \sin\left(\frac{2\pi}{P} \cdot n \cdot x\right) \right)$$

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

我们可以使用这个傅里叶级数，或其中的几项，来建模我们的季节性。在我们的应用中，P 是我们试图建模的周期的最大长度。例如，对于月度数据的年季节性，周期的最大长度（P）为12。x 是一个从 1 到 P 的序数变量。在这个例子中，x 取值为 1，2，3，… 12。现在，有了这些项，剩下的就是找到a_n 和 b_n，我们可以通过对信号进行回归来实现。

我们已经看到，通过合适的傅里叶项组合，我们可以复制任何信号。但问题是，我们应该这样做吗？我们希望从数据中学习的是一个广义的季节性轮廓，它在未见过的数据上也能表现良好。因此，我们将 N 作为一个超参数，从数据中提取所需复杂程度的信号。

现在是复习三角学并回想正弦波和余弦波形貌的好时机。第一个傅里叶项（n=1）是古老的正弦波和余弦波，它们在最大周期长度（P）内完成一个完整周期。随着 n 的增加，我们得到在最大周期长度（

图3.11: 余弦傅里叶项 (n=1, 2, 3)

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

图3.12: 正弦和余弦傅里叶项（n = 1）

现在，让我们看看如何在实际中应用它。

实现

笔记本提示：

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

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

来自 statsmodel 的 seasonal_decompose

statsmodels.tsa.seasonal 有一个名为 seasonal_decompose 的函数。该实现使用滑动平均来处理趋势分量，并使用周期调整平均值来处理季节分量。它支持加法和乘法分解模式。但是，它不能容忍缺失值。让我们看看如何使用它：

#Does not support missing values, so using imputed ts instead
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天进行移动平均：

图3.13：使用 statsmodels 进行季节性分解

我们看不到季节性模式，因为在图的整体尺度中它太小了。随附的笔记本中有放大后的图，以帮助您理解季节性模式。即使使用较大的平滑窗口（例如20天），趋势中仍然存在一些噪声。我们可以通过增大窗口来进一步减少噪声，但有更好的替代方法，我们将马上看到。

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

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

此实现可以在本书的 GitHub 仓库的 src.decomposition.seasonal.py 下找到。它需要一个带有日期时间索引的 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 来捕获工作日-小时模式：

图3.14：STL分解

让我们再看一个月的分解，以便更清晰地观察提取的季节性模式：

图3.15：STL分解（一个月放大）

趋势现在足够平滑，季节性也被捕获了。在这里，我们可以清楚地看到每小时的峰值和谷值，以及周末更高的峰值。然而，由于我们依赖平均值来推导季节性，它也会受到异常值的很大影响。时间序列中几个非常高或非常低的值会使从周期平均值推导出的季节性特征发生偏移。这种技术的另一个缺点是，当数据分辨率与预期季节性周期之间的差异较大时，提取的季节性的“质量”会下降。例如，在日或低于日的数据上提取年季节性，会使提取的季节性非常嘈杂。如果预期季节性的周期少于两个，该技术也将失效——例如，我们想提取年季节性，但数据不足2年。

傅里叶分解

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

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，并在一年中每天增加1：

#Making a custom seasonality term
ts_df["dayofweek"] = ts_df.index.dayofweek
ts_df["hour"] = ts_df.index.hour
#Creating a sorted unique combination df
map_df = ts_df[["dayofweek","hour"]].drop_duplicates().sort_values(["dayofweek", "hour"])
# Assigning an ordinal variable to capture the order
map_df["map"] = np.arange(1, len(map_df)+1)
# mapping the ordinal mapping back to the original df and getting the seasonality array
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 index 中提取的季节性。pandas 日期时间属性，如 week_of_day 和 month，可用于指定最显著的季节性。如果设置为 None，则需要在调用 fit 时提供季节性数组。
model接受additive或multiplicative作为参数，以确定分解的类型。
n_fourier_terms 决定用于提取季节性的傅里叶项数量。增加此参数越多，从数据中提取的季节性越复杂。
lo_frac是用于拟合LOESS回归的数据比例。
lo_delta是使用线性插值代替加权回归的分数距离。使用非零的lo_delta可以显著减少计算时间。

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

图3.16：使用傅里叶项进行分解（放大显示一个月）

趋势将与 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的日期时间属性。
seasonality_model接受fourier或averages作为参数，以确定季节性分解的类型。
model接受additive或multiplicative作为参数，以确定分解的类型。
n_fourier_terms决定了用于提取季节性的傅里叶项数量。随着该参数的增加，从数据中提取的季节性将更加复杂。
lo_frac是用于拟合LOESS回归的数据比例。
lo_delta是使用线性插值代替加权回归的分数距离。使用非零的lo_delta可以显著减少计算时间。

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

图3.17：使用傅里叶项的多重季节性分解

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

图3.18：使用傅里叶项的多重季节性分解（放大到一个月）

在这里，我们可以观察到hour的季节性已被很好地提取出来，并且它还隔离了day_of_week季节性成分，该成分在周末达到峰值。离散阶梯性质的day_of_week季节性成分是因为数据的频率是半小时一次，并且对于48个数据点，day_of_week将是相同的。

我们将所介绍的四种技术总结在下面的表格中：

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

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

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

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

检测与处理异常值

异常值（outlier），顾名思义，是与其他观测值相距异常远的观测值。如果我们把数据生成过程（DGP）视为生成时间序列的随机过程，那么异常值就是从DGP中生成概率最小的点。这可能由多种原因造成，包括测量设备故障、数据录入错误以及黑天鹅事件等。能够检测并处理这些异常值可能帮助您的预测模型更好地理解数据。

异常值/异常检测本身是时间序列中的一个专门领域，但在本书中，我们将局限于更简单的识别和处理异常值的技术。这是因为我们的主要目标不是检测异常值，而是清理数据以使预测模型表现更好。如果您想了解更多关于异常检测的内容，请参阅'延伸阅读'部分，那里有一些入门资源。

现在，让我们来看几种识别异常值的技术。

笔记本提示：

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

标准差

这是一个几乎所有处理过数据一段时间的人都听说过的经验法则——如果 $\mu$ 是时间序列的均值，$\sigma$ 是标准差，那么任何落在 $\mu \pm 3\sigma$ 之外的值就是异常值。其理论基础深深植根于统计学。如果我们假设时间序列的值服从正态分布（一种具有非常优良性质的对称分布），利用概率论，我们可以推导出正态分布曲线下68%的面积位于均值两侧一个标准差之内，约95%的面积位于两个标准差之内，约99%的面积位于三个标准差之内。因此，当我们使用经验法则将边界设为三个标准差时，我们所说的是，如果某个观测值属于该概率分布的概率小于1%，那么它就是异常值。

转向更实际的问题，这个三标准差的阈值绝非神圣不可侵犯。我们需要尝试不同的倍数，并通过主观评估所得结果来确定合适的倍数。倍数越高，异常值就越少。

对于高度季节性的数据，将经验法则直接应用于原始时间序列的简单方法效果不佳。在这种情况下，我们必须使用前面讨论过的任何一种技术对数据进行季节调整，然后将异常值检测应用于残差。如果不这样做，我们可能会将季节性峰值标记为异常值，这不是我们想要的。

这里的另一个关键假设是正态分布。然而，在现实中，我们遇到的许多时间序列可能并不服从正态分布，因此该经验法则会很快失去其理论保证。

IQR

另一种非常相似的技术是使用IQR代替标准差来定义将观测值标记为异常值的上下界。分位数将所有数据按顺序排列，然后将其分成相等的部分，因此每个部分具有相同数量的项。四分位数做同样的事情，但具体将数据分成四个相等的部分。IQR是第三个四分位数（即第75百分位数或0.75分位数）与第一个四分位数（即第25百分位数或0.25分位数）之间的差。上下界定义如下：

上界 = Q3 + n x IQR
下界 = Q1 - n x IQR

这里，IQR = Q3--Q2，n 是决定可接受区域宽度的IQR倍数。

对于观测到异常值频繁出现且变化剧烈的数据集，这种方法比标准差稍微稳健一些。这是因为标准差和均值受数据集中个别点的影响很大。如果在前一种方法中使用$2\sigma$作为经验法则，那么这里使用的是1.5倍IQR。这也与相同的正态分布假设相对应，1.5倍IQR大约相当于$3\sigma$（确切地说是$2.7\sigma$）。在应用该规则之前先行去除季节性的要点同样适用于这里。它适用于我们将要看到的所有技术。

Isolation Forest

Isolation Forest 是一种基于决策树的无监督异常检测算法。典型的异常检测算法会对正常点进行建模，并将任何不符合正常的点视为异常值。然而，Isolation Forest则选择了不同的路径，直接对异常值进行建模。它通过随机划分特征空间来创建决策树森林。该技术基于这样一个假设：异常点位于外围区域，更容易落入树的叶节点。因此，你可以在短分支中找到异常值，而彼此更靠近的正常点则需要更长的分支。任何点的“异常得分”由到达该特定点之前需要遍历的树的深度决定。scikit-learn在sklearn.ensemble.IsolationForest下提供了该算法的实现。除了决策树的标准参数外，这里的关键参数是contamination。默认设置为auto，但可以设置为0到0.5之间的任何值。该参数指定了数据集中你期望的异常值百分比。

但我们必须记住一点，IsolationForest完全不考虑时间，只是标记出那些偏离常态的值。

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

这种基于统计的技术比基本的$\pm \sigma$技术更复杂，但仍使用相同的正态性假设。它基于另一种称为Grubbs检验的统计检验，该检验用于在正态分布数据集中查找单个异常值。ESD通过每步识别并移除一个异常值来迭代地使用Grubbs检验。它还根据剩余点的数量调整临界值。要更详细地理解该检验，请参阅进一步阅读部分，我们提供了关于ESD和S-ESD的一些资源。2017年，来自Twitter Research的Hochenbaum等人提出使用带季节分解的广义ESD作为时间序列异常检测的方法。

我们改编了该算法的现有实现用于我们的案例，该实现可在本书的GitHub仓库中找到。虽然所有其他方法都通过调整几个参数让用户自行确定合适的异常值水平，但S-ESD仅接受期望异常值数量的上界，然后独立地识别异常值。例如，我们设置上界为800，算法在我们处理的数据中识别出大约400个异常值。

参考文献检查：

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

让我们看看如何使用我们回顾过的所有技术来检测异常值：

图3.19：使用不同技术检测到的异常值

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

处理异常值

我们必须回答的第一个问题是，是否应该纠正我们已识别的异常值。自动识别异常值的统计测试应该经过另一次人工验证。如果我们盲目地“处理”异常值，可能会砍掉一个有助于预测时间序列的有价值模式。如果你只预测少量时间序列，那么查看异常值并通过分析其原因将其与现实挂钩仍然是有意义的。

但是当你有数千个时间序列时，人类无法检查所有异常值，因此我们必须求助于自动化技术。常见的做法是用启发式值（如最大值、最小值和第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]. https://arxiv.org/abs/2107.13462.
Hochenbaum, J., Vallis, O., & Kejariwal, A. (2017). Automatic Anomaly Detection in the Cloud Via Statistical Learning. ArXiv, abs/1704.07706. https://arxiv.org/abs/1704.07706.

留下评论！

感谢您从Packt Publishing购买本书——希望您喜欢！您的反馈非常宝贵，有助于我们改进和成长。阅读完成后，请花一点时间在亚马逊上留下评论；只需一分钟，但对像您这样的读者意义重大。

扫描二维码或访问链接以免费获取您选择的一本电子书。

https://packt.link/NzOWQ