Python-时间序列分析秘籍第二版-四-

Python 时间序列分析秘籍第二版（四）

原文：annas-archive.org/md5/7277c6f80442eb633bdbaf16dcd96fad

译者：飞龙

协议：CC BY-NC-SA 4.0

第八章：8 使用统计方法进行异常值检测

加入我们的书籍社区，参与 Discord 讨论

packt.link/zmkOY

除了缺失数据，如第七章《缺失数据处理》所讨论的内容，您可能遇到的另一个常见数据问题是异常值。异常值可以是点异常值、集体异常值或上下文异常值。例如，点异常值发生在某个数据点偏离其余数据集时——有时被称为全局异常值。集体异常值是指一组观察值，它们与总体不同，并且不遵循预期的模式。最后，上下文异常值是指当某个观察值在特定条件或上下文下被认为是异常值时，例如偏离邻近数据点的情况。需要注意的是，对于上下文异常值，如果上下文发生变化，同一观察值可能不再被视为异常值。

在本章中，您将接触到一些实用的统计技术，包括参数方法和非参数方法。在第十四章《使用无监督机器学习进行异常值检测》中，您将深入了解基于机器学习和深度学习的更高级技术。

在文献中，您会发现另一个常用的术语，异常检测，它与异常值检测是同义的。识别异常或异常值的方式和技术是相似的；它们的区别在于上下文和识别后采取的措施。例如，金融交易中的异常交易可能被视为异常，并触发反欺诈调查，以防止其再次发生。在不同的上下文中，研究人员在检查保留与移除这些点的总体影响后，可能会选择简单地删除调查数据中的异常数据点。有时，您可能会决定保留这些异常值，如果它们是自然过程的一部分。换句话说，它们是合法的，并选择使用不受异常值影响的稳健统计方法。

另一个与异常值检测相关的概念是变点检测（CPD）。在变点检测中，目标是预测时间序列数据中的剧烈且有影响的波动（增加或减少）。变点检测包含特定技术，例如累积和（CUSUM）和贝叶斯在线变点检测（BOCPD）。变点检测在许多情况下至关重要。例如，当机器的内部温度达到某个点时，可能会发生故障，或者如果您想了解打折价格是否促进了销售增长时，这也很重要。异常值检测和变点检测之间的区别至关重要，因为有时您可能需要后者。当这两种学科相交时，根据上下文，突发的变化可能表明异常值（异常）的潜在存在。

本章中您将遇到的食谱如下：

• 重新采样时间序列数据

• 使用可视化方法检测异常值

• 使用 Tukey 方法检测异常值

• 使用 z-score 检测异常值

• 使用修改后的 z-score 检测异常值

• 使用其他非参数方法检测异常值

技术要求

你可以从 GitHub 仓库下载 Jupyter Notebooks 和所需的数据集：

• Jupyter Notebook：github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook-Second-Edition/blob/main/code/Ch8/Chapter%208.ipynb

• 数据集：github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook-Second-Edition/tree/main/datasets/Ch8

在本章中，你将使用Numenta 异常基准（NAB）提供的数据集，它提供了用于异常值检测的基准数据集。有关 NAB 的更多信息，请访问他们的 GitHub 仓库：github.com/numenta/NAB。

纽约出租车数据集记录了特定时间戳下纽约市出租车乘客的数量。该数据包含已知的异常值，用于评估我们异常值检测器的性能。数据集包含10,320条记录，时间范围是2014 年 7 月 1 日到2015 年 5 月 31 日。这些观测值是以 30 分钟的间隔记录的，对应于freq = '30T'。

为了准备异常值检测步骤，首先加载你将在整个章节中使用的库：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
plt.rcParams[“figure.figsize”] = [16, 3]

将nyc_taxi.csv数据加载到 pandas DataFrame 中，因为它将在本章中贯穿始终：

file = Path(“../../datasets/Ch8/nyc_taxi.csv”)
nyc_taxi = pd.read_csv(folder / file,
                     index_col=‘timestamp’,
                     parse_dates=True)
nyc_taxi.index.freq = ‘30T’

你可以存储包含异常值的已知日期，这些异常值也被称为真实标签：

nyc_dates =  [
        ”2014-11-01”,
        ”2014-11-27”,
        ”2014-12-25”,
        ”2015-01-01”,
        ”2015-01-27”
]

如果你研究这些日期以获得更多的见解，你会发现类似以下总结的信息：

• 2014 年 11 月 1 日，星期六，是在纽约马拉松赛之前，官方的马拉松赛事是在 2014 年 11 月 2 日，星期日。

• 2014 年 11 月 27 日，星期四，是感恩节。

• 2014 年 12 月 25 日，星期四，是圣诞节。

• 2015 年 1 月 1 日，星期四，是元旦。

• 2015 年 1 月 27 日，星期二，是北美暴风雪事件，在 2015 年 1 月 26 日至 1 月 27 日之间，所有车辆被命令停驶。

你可以绘制时间序列数据，以便对你将在异常值检测中使用的数据有一个直观的了解：

nyc_taxi.plot(title=“NYC Taxi”, alpha=0.6)

这应该生成一个频率为 30 分钟的时间序列：

图 8.1：纽约市出租车时间序列数据的图示

最后，创建一个plot_outliers函数，你将在整个过程中使用该函数：

def plot_outliers(outliers, data, method=‘KNN’,
                 halignment = ‘right’,
                 valignment = ‘bottom’,
                 labels=False):
    ax = data.plot(alpha=0.6)

    if labels:
        for I in outliers’'valu’'].items():
            plt.plot(i[0], i[1],’'r’')
            plt.text(i[0], i[1], ‘'{i[0].date()’',
                         horizontalalignment=halignment,
                         verticalalignment=valignment)
    else:
        data.loc[outliers.index].plot(ax=ax, style’'r’')

    plt.title(‘'NYC Taxi–- {method’')
    plt.xlabel’'dat’'); plt.ylabel’'# of passenger’')
    plt.legend(‘'nyc tax’'‘'outlier’'])
    plt.show()

随着我们推进异常值检测的步骤，目标是看看不同的技术如何捕捉异常值，并将其与真实标签进行比较。

理解异常值

异常值的存在需要特别处理和进一步调查，不能草率决定如何处理它们。首先，你需要检测并识别它们的存在，这一章正是讲解这个内容。领域知识在确定这些被识别为异常值的点、它们对分析的影响以及如何处理它们方面可以发挥重要作用。

异常值可能表示由于过程中的随机变化（称为噪声）、数据输入错误、传感器故障、实验问题或自然变异等原因导致的错误数据。如果异常值看起来像是人工合成的，通常是不希望看到的，例如错误数据。另一方面，如果异常值是过程中的自然一部分，你可能需要重新考虑是否删除它们，而选择保留这些数据点。在这种情况下，你可以依赖于非参数统计方法，这些方法不对底层分布做假设。

一般来说，异常值会在基于强假设的模型构建过程中引发副作用；例如，数据来自高斯（正态）分布。基于假设底层分布的统计方法和测试称为参数方法。

处理异常值没有固定的协议，异常值的影响程度会有所不同。例如，有时你可能需要分别在有异常值和没有异常值的情况下测试你的模型，以了解它们对分析的整体影响。换句话说，并非所有异常值都是相同的，也不应当一视同仁。然而，正如前面所述，拥有领域知识在处理这些异常值时至关重要。

现在，在使用数据集构建模型之前，你需要测试是否存在异常值，以便进一步调查它们的重要性。识别异常值通常是数据清洗和准备过程的一部分，之后才会深入进行分析。

处理异常值的常见方法是删除这些数据点，使它们不成为分析或模型开发的一部分。或者，你也可以希望使用第七章，处理缺失数据中提到的类似技术，如插补和插值，来替换异常值。其他方法，如平滑数据，可能有助于最小化异常值的影响。平滑方法，如指数平滑法，在第十章，使用统计方法构建单变量时间序列模型中有讨论。你也可以选择保留异常值，并使用对它们影响更为坚韧的算法。

有许多知名的异常值检测方法。该研究领域正在发展，范围从基本的统计技术到更先进的方法，利用神经网络和深度学习。在统计方法中，你可以使用不同的工具，例如可视化工具（箱型图、QQ 图、直方图和散点图）、z 分数、四分位距（IQR）和图基围栏，以及诸如 Grubbs 检验、Tietjen-Moore 检验或广义极值学生化偏差（ESD）检验等统计检验。这些方法是基本的、易于解释且有效的。

在你的第一个实例中，你将会在深入异常值检测之前，学习一个重要的时间序列转换技术——重采样。

重采样时间序列数据

在时间序列数据中，一个典型的转换操作是重采样。该过程意味着改变数据的频率或粒度水平。

通常，你对时间序列生成的频率控制有限。例如，数据可能在小的时间间隔内生成和存储，如毫秒、分钟或小时。在某些情况下，数据可能是按较大的时间间隔生成的，如按天、按周或按月。

对时间序列进行重采样的需求可能来源于你的分析性质以及数据需要达到的粒度级别。例如，你可以拥有按天记录的数据，但分析要求数据以周为单位，因此你需要进行重采样。这个过程称为下采样。进行下采样时，你需要提供某种聚合方式，如均值、总和、最小值或最大值等。另一方面，一些情况下需要将数据从每日重采样为每小时。这一过程称为上采样。在上采样时，可能会引入空值行，你可以使用插补或插值技术来填充这些空值。详见第七章，处理缺失数据，其中详细讨论了插补和插值方法。

在这个实例中，你将会使用pandas库来探索如何进行重采样。

如何操作……

在这个实例中，你将使用技术要求部分中早些时候创建的nyc_taxis DataFrame。该数据捕获了30 分钟间隔内的乘客数量。

1. 下采样数据至每日频率。当前，你有10,320条记录，当你将数据重采样至每日时，你需要对数据进行聚合。在这个例子中，你将使用.mean()函数。这将把样本数减少到 215 条记录。术语下采样表示样本数量减少，具体来说，我们正在减少数据的频率。

检查原始 DataFrame 的前五行：

nyc_taxi.head()
>>
            value
timestamp   
2014-07-01 00:00:00    10844
2014-07-01 00:30:00    8127
2014-07-01 01:00:00    6210
2014-07-01 01:30:00    4656
2014-07-01 02:00:00    3820

重采样通过DataFrame.resample()函数来完成。对于按天重采样，你将使用'D'作为日期偏移规则，然后使用.mean()作为聚合方法：

df_downsampled = nyc_taxi.resample('D').mean()
df_downsampled.head()
>>
            value
timestamp   
2014-07-01    15540.979167
2014-07-02    15284.166667
2014-07-03    14794.625000
2014-07-04    11511.770833
2014-07-05    11572.291667

请注意，现在DatetimeIndex的频率为每日频率，且乘客数量现在反映的是每日平均数。检查第一个DatetimeIndex以查看其频率：

df_downsampled.index[0]
>>
Timestamp('2014-07-01 00:00:00', freq='D')

您还可以直接使用.freq属性检查频率：

df_downsampled.index.freq
>>
<Day>

检查下采样后的记录数：

df_downsampled.shape
>>
(215, 1)

事实上，您现在有215条记录。

1. 再次重采样nyc_taxi数据，这次使用 3 天频率。您可以使用偏移字符串'3D'来做到这一点。这次，请使用.sum()方法：

df_downsampled = nyc_taxi.resample('3D').sum()
df_downsampled.head()
>>
            value
timestamp   
2014-07-01   15540.979167
2014-07-02   15284.166667
2014-07-03   14794.625000
2014-07-04   11511.770833
2014-07-05   11572.291667

检查DatetimeIndex的频率：

df_downsampled.index.freq
>>
<3 * Days>

如果您使用df_downsampled.shape，您会注意到记录的数量减少到了 72 条。

1. 现在，将频率改为三（3）个工作日。pandas的默认值是从星期一到星期五。在第六章，处理日期和时间的使用自定义工作日的技巧中，您已经学习了如何创建自定义工作日。目前，您将使用默认的工作日定义。如果您查看上一步骤的输出，2014-07-13是星期天。使用'3B'作为DateOffset会将其推到下一个星期一，即2014-07-14：

df_downsampled = nyc_taxi.resample('3B').sum()
df_downsampled.head()
>>
             value
timestamp   
2014-07-01  745967
2014-07-04  1996347
2014-07-09  3217427
2014-07-14  3747009
2014-07-17  2248113

有趣的输出显示，它跳过了从 2014-07-04 到 2014-07-09 的 5 天时间，然后又从 2014-07-09 跳到 2014-07-14。原因是工作日规则，它规定了我们有两（2）天作为周末。由于该函数是日历感知的，它知道在第一次 3 天的增量之后会有周末，所以它会加上 2 天的周末来跳过它们，从而形成了一个 5 天的跳跃。从 2014-07-04 开始，跳到 2014-07-09，然后从 2017-07-09 跳到 2017-07-14。类似的，从 2014-07-17 跳到 2014-07-22，依此类推。

1. 最后，让我们上采样数据，从 30 分钟间隔（频率）变为 15 分钟频率。这将会在每两个数据之间创建一个空条目（NaN）。您将使用偏移字符串'T'表示分钟，因为'M'用于按月聚合：

df_upsampled = nyc_taxi.resample('15T').mean()
df_upsampled.head()
>>
                     value
timestamp   
2014-07-01 00:00:00    10844.0
2014-07-01 00:15:00    NaN
2014-07-01 00:30:00    8127.0
2014-07-01 00:45:00    NaN
2014-07-01 01:00:00    6210.0

请注意，上采样会创建NaN行。与下采样不同，在上采样时，您需要提供如何填充NaN行的说明。您可能会想知道为什么我们这里使用了.mean()？简单的答案是，无论您是使用.sum()、.max()还是.min()，例如，结果都不会有所不同。在上采样时，您需要使用插补或插值技术来填充缺失的行。例如，您可以使用.fillna()指定插补值，或者使用.ffill()或.bfill()方法：

df_upsampled = nyc_taxi.resample('15T').ffill()
df_upsampled.head()
>>
                    value
timestamp   
2014-07-01 00:00:00    10844
2014-07-01 00:15:00    10844
2014-07-01 00:30:00    8127
2014-07-01 00:45:00    8127
2014-07-01 01:00:00    6210

前五条记录显示了使用前向填充的方法。如需更多关于使用.fillna()、.bfill()或.ffill()，或一般插补的相关信息，请参考第七章，处理缺失数据。

总体而言，pandas 中的重采样非常方便和直接。当您想要改变时间序列的频率时，这是一个很实用的工具。

它是如何工作的…

DataFrame.resample() 方法允许你在指定的时间框架内对行进行分组，例如按天、周、月、年或任何 DateTime 属性。.resample() 的工作原理是通过使用 DatetimeIndex 和提供的频率对数据进行分组，因此，此方法特定于时间序列 DataFrame。

.resample() 函数的工作方式与 .groupby() 函数非常相似；不同之处在于，.resample() 专门用于时间序列数据，并按 DatetimeIndex 进行分组。

还有更多...

在进行降采样时，你可以一次提供多个聚合操作，使用 .agg() 函数。

例如，使用'M'表示按月，你可以为 .agg() 函数提供一个你想要执行的聚合操作列表：

nyc_taxi.resample('M').agg(['mean', 'min',
                            'max', 'median', 'sum'])

这应该生成一个包含五列的 DataFrame，每个聚合方法对应一列。索引列和timestamp列将按月进行分组：

图 8.2：使用 .agg() 方法进行多重聚合

请注意，'M' 或月频率的默认行为是在月末（例如 2014-07-31）。你也可以改为使用 'MS' 来设置为月初。例如，这将生成 2014-07-01（每月的第一天）。

另请参见

要了解更多关于 pandas DataFrame.resample() 的信息，请访问官方文档：pandas.pydata.org/docs/reference/api/pandas.DataFrame.resample.html。

使用可视化方法检测异常值

使用统计技术检测异常值有两种常见方法：参数性方法和非参数性方法。参数性方法假设你知道数据的基础分布。例如，如果你的数据遵循正态分布。另一方面，非参数性方法则没有这样的假设。

使用直方图和箱型图是基本的非参数技术，可以提供数据分布和异常值存在的洞察。更具体地说，箱型图，也叫做 箱线图，提供了五个数值概述：最小值、第一四分位数（25 百分位数）、中位数（50 百分位数）、第三四分位数（75 百分位数）和 最大值。在箱线图中，须的扩展范围有不同的实现方式，例如，须可以延伸到 最小值 和 最大值。在大多数统计软件中，包括 Python 的 matplotlib 和 seaborn 库，须会延伸到所谓的 Tukey's Lower and Upper Fences。任何超出这些边界的数据点都被认为是潜在的异常值。在 使用 Tukey 方法检测异常值 这一配方中，你将深入了解实际的计算和实现。现在，让我们关注分析中的可视化部分。

在这个例子中，你将使用seaborn，这是另一个基于matplotlib的 Python 可视化库。

准备工作

你可以从 GitHub 仓库下载 Jupyter Notebook 和所需的数据集。请参考本章的技术要求部分。你将使用在技术要求部分中加载的 nyc_taxi 数据框。

你将使用 seaborn 版本 0.11.2，这是截至目前的最新版本。

要使用 pip 安装 seaborn，请使用以下命令：

pip install seaborn

要使用 conda 安装 seaborn，请使用以下命令：

conda install seaborn

如何操作...

在本例中，你将探索来自 seaborn 的不同图表，包括 histplot()、displot()、boxplot()、boxenplot() 和 violinplot()。

你会注意到，这些图表传达了相似的故事，但在视觉上，每个图表的呈现方式不同。最终，在调查数据时，你将对这些图表中的某些图表产生偏好，用于自己的分析：

1. 首先，导入 seaborn 库，开始探索这些图表如何帮助你检测异常值：

Import seaborn as sns
sns.__version__
>> '0.13.1'

1. 回想一下图 8.1，nyc_taxi 数据框包含每 30 分钟记录的乘客数量。请记住，每个分析或调查都是独特的，因此你的方法应该与解决的问题相匹配。这也意味着你需要考虑数据准备方法，例如，确定需要对数据应用哪些转换。

对于本例，目标是找出哪些天存在异常观察值，而不是哪一天的哪个时间段，因此你将重采样数据以按天进行分析。你将首先使用 mean 聚合来下采样数据。尽管这样的转换会平滑数据，但由于 mean 对异常值敏感，因此你不会失去太多与查找异常值相关的细节。换句话说，如果某个特定时间段内存在极端异常值（一天有 48 个时间段），mean 仍然会保留这一信息。

在此步骤中，你将使用 .resample() 方法将数据下采样到每日频率，并使用 ‘D’ 作为偏移字符串。这将把观察值的数量从 10,320 减少到 215（即 10320/48 = 215）：

tx = nyc_taxi.resample('D').mean()

1. 绘制新的 tx 数据框，并使用真实标签作为参考。你将调用你在技术要求部分中创建的 plot_outliers 函数：

known_outliers= tx.loc[nyc_dates]
plot_outliers(known_outliers, tx, 'Known Outliers')

这应该会生成一个带有 X 标记的时间序列图，标记出已知的异常值。

图 8.3：使用真实标签（异常值）对纽约出租车数据进行下采样后的绘图

1. 现在，让我们从第一个图开始，使用 histplot() 函数检查你的时间序列数据：

sns.histplot(tx)

这应该会生成以下结果：

图 8.4：显示极端每日平均乘客乘车次数的直方图

在图 8.4中，标记为1、2、3、4和5的观测值似乎代表了极端的乘客数值。回想一下，这些数字表示的是重采样后的平均每日乘客数量。你需要问自己的是，这些观测值是否是离群值。直方图的中心接近 15,000 每日平均乘客数量。这应该让你质疑接近 20,000 的极端值（标签 5）是否真的那么极端。类似地，标记为3和4的观测值（因为它们接近分布的尾部），它们是否真的是极端值？那么标记为1和2的观测值呢？它们的平均乘客数量分别为 3,000 和 8,000 每日乘客，似乎比其他值更极端，可能是实际的离群值。再一次，确定什么是离群值、什么不是离群值需要领域知识和进一步分析。没有具体的规则，你会发现在本章中，一些被普遍接受的规则实际上是任意的、主观的。你不应急于得出结论。

1. 你可以使用displot()绘制类似的图表，它具有一个kind参数。kind参数可以取三个值之一：hist表示直方图，kde表示核密度估计图，ecdf表示经验累积分布函数图。

你将使用displot(kind='hist')来绘制一个类似于图 8.4中的直方图：

sns.displot(tx, kind='hist', height=3, aspect=4)

1. 箱型图提供的信息比直方图更多，可以更好地帮助识别离群值。在箱型图中，超出须状线或边界的观测值被认为是潜在的离群值。须状线表示上界和下界的视觉边界，这是数学家约翰·图基于 1977 年提出的。

以下代码展示了如何使用seaborn创建箱型图：

sns.boxplot(tx['value'], orient='h')

以下图表显示了潜在的离群值：

图 8.5：显示可能的离群值超出边界（须状线）的箱型图

箱体的宽度（Q1到Q3）称为四分位数范围（IQR），其计算方式为 75th 和 25th 百分位数的差值（Q3 – Q1）。下边界的计算公式为Q1 - (1.5 x IQR)，上边界为Q3 + (1.5 x IQR)。任何小于下边界或大于上边界的观测值都被视为潜在的异常值。有关更多信息，请参阅使用 Tukey 方法检测异常值一文。boxplot()函数中的whis参数默认为1.5（1.5 倍 IQR），控制上下边界之间的宽度或距离。较大的值意味着较少的观测值会被视为异常值，而较小的值则会让非异常值看起来超出边界（更多的异常值）。

sns.boxplot(tx['value'], orient='h', whis=1.5)

1. seaborn中有两种额外的箱型图变种（boxenplot() 和 violinplot()）。它们提供与箱型图相似的见解，但呈现方式不同。boxen 图（在文献中称为字母值图）可以视为对常规箱型图的增强，旨在解决它们的一些不足之处，具体描述见论文Heike Hofmann, Hadley Wickham & Karen Kafadar (2017) Letter-Value Plots: Boxplots for Large Data, Journal of Computational and Graphical Statistics, 26:3, 469-477。更具体地说，boxen（字母值）图更适合用于处理较大数据集（用于显示数据分布的观测值较多，更适合区分较大数据集中的异常值）。seaborn实现的boxenplot基于该论文。

以下代码展示了如何使用seaborn创建 boxen（字母值）图：

sns.boxenplot(tx['value'], orient='h')

这将生成类似于图 8.5中箱型图的图形，但箱体会延伸超出四分位数（Q1、Q2和Q3）。25th 百分位数位于 14,205 每日平均乘客数的位置，75th 百分位数位于 16,209 每日平均乘客数的位置。

图 8.6：纽约市每日出租车乘客的 boxen（字母值）图

百分位值

你是否好奇我如何确定 25th、50th 和 75th 百分位的精确值？

你可以通过describe方法为 DataFrame 或序列获取这些值。例如，如果你运行tx.describe()，你应该能看到一张描述性统计表，其中包括数据集的计数、均值、标准差、最小值、最大值、25th、50th 和 75th 百分位数值。

在图 8.6中，你可以获得超越分位数的乘客分布的额外洞察。换句话说，它扩展了箱形图以显示额外的分布，从而为数据的尾部提供更多的洞察。理论上，这些框可以一直延伸，以容纳所有数据点，但为了显示异常值，需要有一个停止点，这个点称为深度。在seaborn中，这个参数被称为k_depth，它可以接受一个数值，或者你可以指定不同的方法，如tukey、proportion、trustworthy或full。例如，k_depth=1的数值将显示与图 8.5中的箱形图相似的框（一个框）。作为参考，图 8.6显示了使用Tukey方法确定的四个框，这是默认值（k_depth="tukey"）。使用k_depth=4将产生相同的图形。

这些方法在*Heike Hofmann、Hadley Wickham & Karen Kafadar（2017）*的参考论文中有所解释。要探索不同的方法，你可以尝试以下代码：

for k in ["tukey", "proportion", "trustworthy", "full"]:
    sns.boxenplot(tx['value'], orient='h', k_depth=k, orient='h')
    plt.title(k)
    plt.show()

这将产生四个图形；请注意，每种方法确定的框数不同。回顾一下，你也可以按数值指定k_depth：

图 8.7：seaborn 中 boxenplot 函数可用的不同 k_depth 方法

1. 现在，最后的变体是提琴图，你可以使用violinplot函数来显示：

sns.violinplot(tx['value'], orient='h')

这将产生一个介于箱形图和核密度估计（KDE）之间的混合图。核是一个估算概率密度函数的函数，较大的峰值（较宽的区域），例如，表示大多数点集中所在的区域。这意味着一个数据点出现在该区域的概率较高，而在较薄的区域，概率则较低。

图 8.8：纽约市每日出租车乘客的提琴图

请注意，图 8.8显示了整个数据集的分布。另一个观察点是峰值的数量；在这种情况下，我们有一个峰值，这使得它成为一个单峰分布。如果有多个峰值，我们称之为多峰分布，这应该引发对数据的进一步调查。KDE 图将提供与直方图类似的洞察，但其曲线更加平滑。

它是如何工作的...

在本食谱中，我们介绍了几种图表，它们有助于可视化数据的分布并显示异常值。通常，直方图很适合显示分布，但箱形图（及其变体）在异常值检测方面要好得多。我们还探讨了箱线（字母值）图，它更适用于较大的数据集，比普通的箱形图更为合适。

还有更多...

在 seaborn 中，histplot() 和 displot() 都支持通过参数 kde 添加核密度估计（KDE）图。以下代码片段应该会生成一个相似的图形：

sns.histplot(tx, kde=True)
sns.displot(tx, kind='hist', height=3, aspect=4, kde=True)

这应该会生成一个结合直方图和 KDE 图的图形，如下所示：

图 8.9：带有 KDE 图的直方图

此外，另一个有用的可视化方法用于发现异常值是 滞后图。滞后图本质上是一个散点图，但与绘制两个变量以观察相关性不同，举个例子，我们将同一变量与其滞后版本进行比较。这意味着，它是一个使用相同变量的散点图，但 y 轴表示当前时刻 (t) 的乘客数量，x 轴显示的是前一个时刻 (t-1) 的乘客数量，这被称为 滞后。滞后参数决定了回溯的周期数；例如，滞后 1 表示回溯一个周期，滞后 2 表示回溯两个周期。在我们的重采样数据（降采样至每日）中，滞后 1 代表前一天。

pandas 库提供了 lag_plot 函数，你可以像下面的示例所示使用它：

from pandas.plotting import lag_plot
lag_plot(tx)

这应该会生成以下散点图：

图 8.10：纽约市日均出租车乘客的滞后图

被圈出的数据点突出显示了可能的异常值。有些点看起来比其他点更加极端。此外，你还可以看到乘客数量与其滞后版本（前一天）的线性关系，表明存在自相关性。回顾基础统计学中的相关性，相关性显示了两个独立变量之间的关系，因此你可以把自相关看作是一个变量在某一时刻 (t) 和它在前一时刻 (t-1) 的版本之间的相关性。更多内容请参见 第九章，探索性数据分析与诊断 和 第十章，使用统计方法构建单变量时间序列模型。

在 图 8.10 中，x 轴和 y 轴的标签可能会有些混淆，y 轴 被标记为 y(t+1)。实际上，这表示了我们之前描述的同一个意思：x 轴表示的是先前的值（预测变量），而 y 轴表示的是它的未来值 t+1。为了更清晰地理解，你可以像下面的代码所示，手动使用 seaborn 重现 lag_plot 所生成的可视化效果：

y = tx[1:].values.reshape(-1)
x = tx[:-1].values.reshape(-1)
sns.scatterplot(x=x, y=y)

这应该会生成一个与 图 8.10 相似的图形。

注意代码中，y 值从 t+1 开始（我们跳过了索引 0 处的值），直到最后一个观测值，而 x 值从索引 0 开始，到索引 -1（我们跳过了最后一个观测值）。这使得 y 轴上的值领先一个周期。

在下一个食谱中，我们将进一步探讨 IQR 和 Tukey fences，这两个概念我们在讨论箱形图时简要提到过。

另见

你可以通过seaborn文档了解我们使用的图表以及不同的选项。要了解更多信息，请访问相关网址：

• 对于箱型图（boxplot），你可以访问 seaborn.pydata.org/generated/seaborn.boxplot.html。

• 对于箱型图（boxenplot），你可以访问 seaborn.pydata.org/generated/seaborn.boxenplot.html。

• 对于小提琴图（violinplot），你可以访问 seaborn.pydata.org/generated/seaborn.violinplot.html#seaborn.violinplot。

• 对于直方图（histplot），你可以访问 seaborn.pydata.org/generated/seaborn.histplot.html。

• 对于分布图（displot），你可以访问 seaborn.pydata.org/generated/seaborn.displot.html#seaborn.displot。

使用 Tukey 方法检测异常值

这个示例将扩展前一个示例，使用可视化检测异常值。在图 8.5中，箱型图显示了四分位数，须状线延伸至上下边界。这些边界或围栏是使用 Tukey 方法计算得出的。

让我们在图 8.5的基础上，扩展一些其他组件的信息：

图 8.11：每日平均出租车乘客数据的箱型图

可视化图表非常有助于你对所处理数据的整体分布和潜在异常值有一个高层次的了解。最终，你需要通过编程识别这些异常值，以便隔离这些数据点，进行进一步的调查和分析。这个示例将教你如何计算 IQR，并定义落在 Tukey 围栏外的点。

如何操作...

大多数统计方法允许你发现超出某个阈值的极端值。例如，这个阈值可能是均值、标准差、第 10 或第 90 百分位数，或者其他你想要比较的值。你将通过学习如何获取基本的描述性统计量，特别是分位数，来开始这个示例。

1. DataFrame 和 Series 都有describe方法，用于输出总结性描述性统计量。默认情况下，它显示四分位数：第一四分位数，即第 25 百分位数，第二四分位数（中位数），即第 50 百分位数，第三四分位数，即第 75 百分位数。你可以通过向percentiles参数提供一个值列表来定制百分位数。以下代码演示了如何获取额外百分位数的值：

percentiles = [0, 0.05, .10, .25, .5, .75, .90, .95, 1]
tx.describe(percentiles= percentiles)

这应该生成以下 DataFrame：

图 8.12：包含自定义百分位数的每日出租车乘客数据描述性统计

分位数与四分位数与百分位数

这些术语可能会让人混淆，但本质上，百分位数和四分位数都是分位数。有时你会看到人们更宽松地使用百分位数，并将其与分位数交替使用。

四分位数将数据分为四个部分（因此得名），分别标记为Q1（第 25 百分位）、Q2（第 50 百分位或中位数）和Q3（第 75 百分位）。百分位数则可以取 0 到 100 之间的任何范围（在 pandas 中为 0 到 1，在 NumPy 中为 0 到 100），但最常见的是将数据分为 100 个部分。这些部分称为分位数。

这些名称基本上表示应用于数据分布的分割类型（分位数的数量）；例如，四个分位数称为四分位数，两个分位数称为中位数，十个分位数称为十分位数，100 个分位数称为百分位数。

1. NumPy 库还提供了percentile函数，该函数返回指定百分位数的值。以下代码解释了如何使用该函数：

percentiles = [0, 5, 10, 25, 50, 75, 90, 95, 100]
np.percentile(tx, percentiles)
>>
array([ 4834.54166667, 11998.18125   , 13043.85416667, 14205.19791667,
       15299.9375    , 16209.42708333, 17279.3       , 18321.61666667,
       20553.5       ])

1. 在图 8.11中，注意到大多数极端值和潜在异常值位于低界限以下（低界限计算公式为Q1 – (1.5 x IQR)）或位于高界限以上（高界限计算公式为Q3 + (1.5 x IQR)）。IQR 是Q3与Q1的差值（IQR = Q3 – Q1），它决定了箱形图中箱体的宽度。这些上下界限被称为Tukey 界限，更具体地说，它们被称为内界限。外界限也有更低的Q1 - (3.0 x IQR)和更高的Q3 + (3.0 x IQR)界限。我们将重点关注内界限，并将任何超出这些界限的数据点视为潜在异常值。

你将创建一个名为iqr_outliers的函数，该函数计算 IQR、上界（内界限）、下界（内界限），然后过滤数据以返回异常值。这些异常值是任何低于下界或高于上界的数据点：

def iqr_outliers(data):
    q1, q3 = np.percentile(data, [25, 75])
    IQR = q3 - q1
    lower_fence = q1 - (1.5 * IQR)
    upper_fence = q3 + (1.5 * IQR)
    return data[(data.value > upper_fence) | (data.value < lower_fence)]

1. 通过传递tx数据框来测试该函数：

outliers = iqr_outliers(tx)
outliers
>>    
                  value
timestamp              
2014-11-01  20553.500000
2014-11-27  10899.666667
2014-12-25   7902.125000
2014-12-26  10397.958333
2015-01-26   7818.979167
2015-01-27   4834.541667 

这些日期（点）与图 8.5和图 8.11中根据 Tukey 界限识别的异常值相同。

1. 使用技术要求部分中定义的plot_outliers函数：

plot_outliers(outliers, tx, "Outliers using IQR with Tukey's Fences")

这应生成类似于图 8.3中的图表，只不过x标记是基于 Tukey 方法识别的异常值：

图 8.13：使用 Tukey 方法识别的每日平均出租车乘客数和异常值

比较图 8.13和图 8.3，你会看到这个简单的方法在识别已知的五个异常值中的四个方面做得非常好。此外，Tukey 方法还识别出了在2014-12-26和2015-01-26的两个额外的异常值。

它是如何工作的……

使用 IQR 和 Tukey’s fences 是一种简单的非参数统计方法。大多数箱线图实现使用 1.5x(IQR) 来定义上下界限。

还有更多……

使用 1.5x(IQR) 来定义异常值是很常见的；尽管有很多关于其理由的讨论，但这个选择仍然是任意的。你可以更改这个值进行更多实验。例如，在 seaborn 中，你可以通过更新 boxplot 函数中的 whis 参数来改变默认的 1.5 值。当数据符合高斯分布（正态分布）时，选择 1.5 是最有意义的，但这并不总是如此。一般来说，值越大，你捕获的异常值就越少，因为你扩展了边界（界限）。同样，值越小，更多的非异常值会被定义为异常值，因为你缩小了边界（界限）。

让我们更新 iqr_outliers 函数，接受一个 p 参数，这样你就可以尝试不同的值：

def iqr_outliers(data, p):
    q1, q3 = np.percentile(data, [25, 75])
    IQR = q3 - q1
    lower_fence = q1 - (p * IQR)
    upper_fence = q3 + (p * IQR)
    return data[(data.value > upper_fence) | (data.value < lower_fence)]

在不同值上运行函数：

for p in [1.3, 1.5, 2.0, 2.5,  3.0]:
    print(f'with p={p}')
    print(iqr_outliers(tx, p))
    print('-'*15)
>>
with p=1.3
                   value
timestamp              
2014-07-04  11511.770833
2014-07-05  11572.291667
2014-07-06  11464.270833
2014-09-01  11589.875000
2014-11-01  20553.500000
2014-11-08  18857.333333
2014-11-27  10899.666667
2014-12-25   7902.125000
2014-12-26  10397.958333
2015-01-26   7818.979167
2015-01-27   4834.541667
---------------
with p=1.5
                   value
timestamp              
2014-11-01  20553.500000
2014-11-27  10899.666667
2014-12-25   7902.125000
2014-12-26  10397.958333
2015-01-26   7818.979167
2015-01-27   4834.541667
---------------
with p=2.0
                   value
timestamp              
2014-11-01  20553.500000
2014-12-25   7902.125000
2015-01-26   7818.979167
2015-01-27   4834.541667
---------------
with p=2.5
                  value
timestamp             
2014-12-25  7902.125000
2015-01-26  7818.979167
2015-01-27  4834.541667
---------------
with p=3.0
                  value
timestamp             
2014-12-25  7902.125000
2015-01-26  7818.979167
2015-01-27  4834.541667
---------------

最佳值将取决于你的数据以及你需要多敏感的异常值检测。

另请参见

要了解更多关于 Tukey’s fences 用于异常值检测的内容，可以参考这个维基百科页面：en.wikipedia.org/wiki/Outlier#Tukey's_fences。

我们将在接下来的配方中探索另一种基于 z-score 的统计方法。

使用 z-score 检测异常值

z-score 是一种常见的数据标准化变换。当你需要比较不同的数据集时，这种方法非常常见。例如，比较来自两个不同数据集的两个数据点相对于它们的分布会更容易。之所以能做到这一点，是因为 z-score 将数据标准化，使其围绕零均值居中，并且单位表示的是偏离均值的标准差。例如，在我们的数据集中，单位是按日计的出租车乘客数量（以千人计）。一旦应用了 z-score 变换，你将不再处理乘客数量，而是单位表示的是标准差，这告诉我们一个观测值距离均值有多远。以下是 z-score 的公式：

其中 是一个数据点（观测值），mu () 是数据集的均值，sigma () 是数据集的标准差。

请记住，z-score 是一种无损变换，这意味着你不会丢失诸如数据分布（数据形状）或观测之间关系的信息。唯一改变的是度量单位，它们正在被缩放（标准化）。

一旦使用 z-score 转换数据，您就可以选择一个阈值。所以，任何高于或低于该阈值的数据点（以标准差为单位）都被视为异常值。例如，您的阈值可以设置为离均值+3和-3标准差。任何低于-3或高于+3标准差的点可以视为异常值。换句话说，点离均值越远，它作为异常值的可能性就越大。

z-score 有一个主要的缺点，因为它是基于假设的参数统计方法。它假设数据服从高斯（正态）分布。那么，假设数据不是正态分布，您将需要使用修改版的 z-score，这将在下一个部分中讨论，使用修改版 z-score 检测异常值。

操作方法...

您将首先创建 zscore 函数，该函数接受一个数据集和一个阈值，我们称其为 degree。该函数将返回标准化后的数据和识别出的异常值。这些异常值是指任何高于正阈值或低于负阈值的点。

1. 创建 zscore() 函数来标准化数据，并根据阈值过滤掉极端值。请记住，阈值是基于标准差的：

def zscore(df, degree=3):
    data = df.copy()
    data['zscore'] = (data - data.mean())/data.std()
    outliers = data[(data['zscore'] <= -degree) | (data['zscore'] >= degree)]

    return outliers['value'], data

1. 现在，使用 zscore 函数并存储返回的对象：

threshold = 2.5
outliers, transformed = zscore(tx, threshold)

1. 要查看 z-score 转换的效果，您可以绘制一个直方图。转换后的 DataFrame 包含两列数据，原始数据标记为 value，标准化数据标记为 zscore：

transformed.hist()

这应该会生成两个直方图，分别对应两列数据：

图 8.14：原始数据和标准化数据分布对比的直方图

注意数据的形状没有变化，这也是为什么 z-score 被称为 无损转换 的原因。两者唯一的区别是尺度（单位）。

1. 您使用阈值 2.5 运行了 zscore 函数，意味着任何与均值相距 2.5 标准差的数据点（无论正负方向）都将被视为异常值。例如，任何高于 +2.5 标准差或低于 -2.5 标准差的数据点都将被视为异常值。打印出捕获在 outliers 对象中的结果：

print(outliers)
>>
timestamp
2014-11-01    20553.500000
2014-12-25     7902.125000
2015-01-26     7818.979167
2015-01-27     4834.541667
Name: value, dtype: float64

这种简单的方法成功地捕捉到了五个已知异常值中的三个。

1. 使用在 技术要求 部分中定义的 plot_outliers 函数：

plot_outliers(outliers, tx, "Outliers using Z-score")

这应该会生成类似于 图 8.3 中的图形，只不过 x 标记是基于使用 z-score 方法识别的异常值：

图 8.15：每日平均出租车乘客数和使用 z-score 方法识别的异常值

你需要多尝试几次，确定最佳阈值。阈值越大，你捕捉到的异常值就越少；阈值越小，更多的非异常值会被标记为异常值。

1. 最后，让我们创建一个plot_zscore函数，该函数接受标准化数据，并使用阈值线绘制数据。这样，你可以直观地看到阈值如何隔离极端值：

def plot_zscore(data, d=3):
    n = len(data)
    plt.figure(figsize=(8,8))
    plt.plot(data,'k^')
    plt.plot([0,n],[d,d],'r--')
    plt.plot([0,n],[-d,-d],'r--')

使用阈值2.5运行该函数：

data = transformed['zscore'].values
plot_zscore(data, d=2.5)

这应该生成一个包含两条水平线的散点图：

图 8.16：基于阈值线的标准化数据和异常值的图示

四个被圈出的数据点代表了zscore函数返回的异常值。使用不同的阈值运行该函数，以更深入理解这个简单的技术。

它是如何工作的...

z 分数方法是一种非常简单且易于解释的方法。z 分数被解释为与均值的标准差单位距离，均值是分布的中心。由于我们从所有观测值中减去均值，本质上是在进行数据的均值中心化。我们还通过标准差进行除法，以标准化数据。

图 8.15几乎解释了这种方法的理解。一旦数据被标准化，使用标准差阈值就变得非常容易。如果数据没有标准化，可能很难基于日常乘客量来确定阈值。

还有更多...

z 分数是一种参数化方法，假设数据来自高斯（正态）分布。statsmodels库中有几种测试可以检测数据是否符合正态分布。这些测试中的一个是Kolmogorov-Smirnov检验。零假设是数据来自正态分布。该检验返回检验统计量和p值；如果p值小于0.05，你可以拒绝零假设（数据不服从正态分布）。否则，你将无法拒绝零假设（数据服从正态分布）。

你将使用来自statsmodels库的kstest_normal函数。为了让结果更易于解释，创建test_normal函数如下：

from statsmodels.stats.diagnostic import kstest_normal
def test_normal(df):
    t_test, p_value = kstest_normal(df)
    if p_value < 0.05:
        print("Reject null hypothesis. Data is not normal")
    else:
        print("Fail to reject null hypothesis. Data is normal")

使用test_normal函数运行检验：

test_normal(tx)
>>
Reject null hypothesis. Data is not normal

正如预期的那样，该数据集不符合正态分布。在第九章《探索性数据分析与诊断》中，你将学习更多的正态性检验方法，具体内容在应用幂次变换配方中。但是要小心，这些检验通常在存在异常值时会失败。如果你的数据未通过正态性检验，那么可以使用《使用可视化检测异常值》一节中讨论的一些绘图方法，检查任何可能导致检验失败的异常值。

另请参见

要了解更多关于 z-score 和标准化的内容，可以参考这篇维基百科页面：en.wikipedia.org/wiki/Standard_score。

在接下来的方法中，你将探索一种与 z-score 非常相似的方法，这种方法对异常值更为鲁棒，并且更适合非正态数据。

使用修改版 z-score 检测异常值

在使用 z-score 检测异常值这一方法中，你体验了该方法的简单性和直观性。但是它有一个主要缺点：假设你的数据是正态分布的。

但是，如果你的数据不是正态分布的怎么办？幸运的是，存在一种修改版的 z-score，适用于非正态数据。常规 z-score 和修改版 z-score 之间的主要区别在于，我们用中位数代替均值：

其中(tilde x)是数据集的中位数，MAD 是数据集的中位绝对偏差：

0.6745值是标准差单位，表示高斯分布中第 75 百分位数（Q3），用于作为归一化因子。换句话说，它用于近似标准差。这样，你从该方法中获得的单位是以标准差为度量的，类似于你如何解释常规 z-score。

你可以使用 SciPy 的百分位点函数（PPF），也称为累积分布函数（CDF）的反函数，来获得此值。只需为 PPF 函数提供一个百分位数，例如 75%，它将返回对应的下尾概率的分位数。

import scipy.stats as stats
stats.norm.ppf(0.75)
>>
0.6744897501960817

这是公式中使用的归一化因子。

最后，修改版 z-score 有时也被称为鲁棒 z-score。

如何实现...

总体而言，该方法的工作原理与使用常规 z-score 方法时的步骤完全相同。你将首先创建modified_zscore函数，该函数接收一个数据集和一个我们称之为degree的阈值，然后该函数将返回标准化后的数据以及识别出的异常值。这些异常值是指超出正阈值或低于负阈值的点。

1. 创建modified_zscore ()函数来标准化数据，并根据阈值过滤掉极端值。回想一下，阈值是基于标准差的：

def modified_zscore(df, degree=3):
    data = df.copy()
    s = stats.norm.ppf(0.75)
    numerator = s*(data - data.median())
    MAD = np.abs(data - data.median()).median()
    data['m_zscore'] = numerator/MAD
    outliers = data[(data['m_zscore'] > degree) | (data['m_zscore'] < -degree)]

    return outliers['value'], data 

1. 现在，使用modified_zscore函数并存储返回的对象：

threshold = 3
outliers, transformed = modified_zscore (tx, threshold)

1. 为了查看修改版 z-score 转换的效果，让我们绘制一个直方图。转换后的 DataFrame 包含两列数据，原始数据列标为value，标准化数据列标为zscore。

transformed.hist()

这应该会生成两个直方图，分别对应两列数据：

图 8.17：比较原始和修改版 z-score 标准化数据分布的直方图

比较图 8.16与图 8.13中的结果。两种方法，z-score 和修改后的 z-score 方法，都没有改变数据的形状。不同之处在于缩放因子。

1. 使用modified_zscore函数，设置阈值为3，这意味着任何数据点距离中位数三倍标准差的距离（无论方向如何）都会被视为异常值。例如，任何高于+3标准差或低于-3标准差的数据点都将被视为异常值。打印出outliers对象中捕获的结果：

print(outliers)
>>
timestamp
2014-11-01    20553.500000
2014-11-27    10899.666667
2014-12-25     7902.125000
2014-12-26    10397.958333
2015-01-26     7818.979167
2015-01-27     4834.541667
Name: value, dtype: float64

有趣的是，修改后的 z-score 在捕捉五个已知异常值中的四个时表现得更好。

1. 使用前面在技术要求部分定义的plot_outliers函数：

plot_outliers(outliers, tx, "Outliers using Modified  Z-score")

这应该会生成类似于图 8.3中的图表，只不过x标记是基于使用修改后的 z-score 方法识别的异常值：

图 8.18：使用修改后的 z-score 方法识别的每日平均出租车乘客及异常值

你需要调整阈值来确定最佳值。阈值越大，你捕获的异常值就越少；阈值越小，更多的非异常值会被标记为异常值。

1. 最后，让我们创建一个plot_m_zscore函数，接受标准化数据并绘制带有阈值线的数据。通过这种方式，你可以直观地看到阈值如何隔离极端值：

def plot_m_zscore(data, d=3):
    n = len(data)
    plt.figure(figsize=(8,8))
    plt.plot(data,'k^')
    plt.plot([0,n],[d,d],'r--')
    plt.plot([0,n],[-d,-d],'r--')

使用阈值为3来运行该函数：

data = transformed['m_zscore'].values
plot_m_zscore(data, d=3)

这应该会生成一个带有两条水平线的散点图：

图 8.19：基于阈值线的标准化数据和异常值的图

六个圈中的数据点代表了由modified_score函数返回的异常值。使用不同的阈值运行此函数，以便对这种简单的技术有更深的直觉理解。

注意在图 8.19中，我们有一个数据点正好位于阈值线处。你认为这是异常值吗？通常，异常值检测时，你仍然需要进行尽职调查来检查结果。

它是如何工作的...

修改后的 z-score（稳健 z-score）方法与 z-score 方法非常相似，因为它依赖于定义标准差阈值。使得这种方法对异常值更加稳健的是它使用了中位数而不是均值。我们还使用了中位数绝对偏差（MAD）来代替标准差。

还有更多...

在前一个配方中，使用 z-score 检测异常值，我们使用了statsmodels中的kstest_normal来测试正态性。

另一个有用的图形是专门用来测试正态性并且有时可以帮助检测异常值的分位数-分位数图（QQ 图）。

你可以使用 SciPy 或 statsmodels 绘制 QQ 图。两者都会生成相同的图。以下代码展示了你如何使用其中任何一个绘图。

这展示了如何使用 SciPy 绘图：

import scipy
import matplotlib.pyplot as plt
res = scipy.stats.probplot(tx.values.reshape(-1), plot=plt)

这展示了如何使用statsmodels绘图：

from statsmodels.graphics.gofplots import qqplot
qqplot(tx.values.reshape(-1), line='s')
plt.show()

无论是 SciPy 还是statsmodels，都会生成以下图形：

图 8.20：QQ 图，比较出租车乘客数据与假设的正态分布

实线代表正态分布数据的参考线。如果你比较的数据是正态分布的，所有数据点将会落在这条直线上。在图 8.19中，我们可以看到分布几乎是正态的（虽然不完美），并且在分布的尾部存在一些问题。这与我们在图 8.16和图 8.13中看到的情况一致，显示大多数异常值位于底部尾部（低于-2标准差）。

另见

若想了解更多关于 MAD 的信息，你可以参考维基百科页面：en.wikipedia.org/wiki/Median_absolute_deviation。

第九章：9 探索性数据分析与诊断

加入我们在 Discord 上的书籍社区

packt.link/zmkOY

到目前为止，我们已经涵盖了从不同来源提取数据的技术。这些内容在第二章，从文件读取时间序列数据和第三章，从数据库读取时间序列数据中已经讨论过了。第六章，在 Python 中处理日期和时间，以及第七章，处理缺失数据，介绍了几种有助于准备、清理和调整数据的技术。

你将继续探索额外的技术，以便更好地理解数据背后的时间序列过程。在建模数据或进行进一步分析之前，一个重要步骤是检查手头的数据。更具体地说，你需要检查一些特定的时间序列特征，如平稳性、趋势和季节性的影响、以及自相关等。这些描述你所处理的时间序列过程的特征需要与该过程背后的领域知识相结合。

本章将建立在你从前几章学到的知识基础上，帮助你为从第十章开始创建和评估预测模型做准备，构建单变量时间序列模型（使用统计方法）。

在本章中，你将学习如何可视化时间序列数据，如何将时间序列分解为其组成部分（趋势、季节性和残差随机过程），如何检验模型可能依赖的不同假设（如平稳性、正态性和同方差性），以及如何探索数据转换技术以满足其中的一些假设。

本章中你将遇到的食谱如下：

• 使用 pandas 绘制时间序列数据

• 使用 hvPlot 绘制交互式时间序列数据

• 分解时间序列数据

• 检测时间序列的平稳性

• 应用幂次变换

• 测试时间序列数据的自相关性

技术要求

你可以从 GitHub 仓库下载所需的 Jupyter 笔记本和数据集，进行跟随学习：

• Jupyter 笔记本： github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook./blob/main/code/Ch9/Chapter%209.ipynb

• 数据集： github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook./tree/main/datasets/Ch9

从本章起，我们将广泛使用 pandas 2.2.0（2024 年 1 月 20 日发布）。这适用于本章中的所有食谱。

我们将使用四个额外的库：

• hvplot 和 PyViz

• seaborn

• matplotlib

如果你使用的是pip，你可以通过终端安装这些包，命令如下：

pip install hvplot seaborn matplotlib jupyterlab

如果你使用的是conda，你可以通过以下命令安装这些包：

conda install jupyterlab matplotlib seaborn
conda install -c pyviz hvplot

HvPlot 库将用于在 JupyterLab 中构建交互式可视化。如果你使用的是最新版本的 JupyterLab（jupyterlab >= 3.0），那么所有所需的扩展都会自动安装，因为它们已经捆绑在 pyviz_comms 包中。如果你使用的是较旧版本的 JupyterLab（jupyterlab < 3.0），那么你需要手动安装 jupyterlab_pyviz 扩展，具体步骤如下：

jupyter labextension install @pyviz/jupyterlab_pyviz

在本章中，你将使用三个数据集（Closing Price Stock Data、CO2 和 Air Passengers）。CO2 和 Air Passengers 数据集由 statsmodels 库提供。Air Passengers 数据集包含了 1949 年至 1960 年的每月航空乘客人数。CO2 数据集包含了 1958 年至 2001 年间，位于毛纳罗亚山的每周大气二氧化碳浓度。Closing Price Stock Data 数据集包括 2019 年 11 月到 2021 年 11 月的微软、苹果和 IBM 的股票价格。

为了开始，你需要加载数据集，并将其存储为 pandas DataFrame，并加载在整个过程中需要的任何库或方法：

import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
from statsmodels.datasets import co2, get_rdataset
file = Path(‘../../datasets/Ch9/closing_price.csv’)
closing_price = pd.read_csv(file,
                            index_col=‘Date’,
                            parse_dates=True)
co2_df = co2.load_pandas().data
co2_df = co2_df.ffill()
air_passengers = get_rdataset(“AirPassengers”)
airp_df = air_passengers.data
airp_df.index = pd.date_range(‘1949’, ‘1961’, freq=‘M’)
airp_df.drop(columns=[‘time’], inplace=True)

现在，你应该已经有了三个数据框：airp_df、closing_price 和 co2_df。

使用 pandas 绘制时间序列数据

可视化是数据分析中的一个关键方面，尤其在处理时间序列数据时显得尤为重要。在前面的章节和示例中，你已经遇到过多次绘制数据的实例，这些实例对于突出特定点或得出关于时间序列的结论至关重要。可视化我们的时间序列数据使我们能够一眼识别出模式、趋势、异常值以及其他关键信息。此外，数据可视化有助于跨不同小组之间的沟通，并通过提供一个共同的平台促进各利益相关方（如商业专家和数据科学家）之间的建设性对话。

在时间序列分析以及机器学习领域，我们在探索性数据分析（EDA）中优先可视化数据，以全面理解我们正在处理的数据。在评估模型时，我们也依赖于可视化来比较它们的表现，并识别需要改进的地方。可视化在模型可解释性方面发挥着关键作用，使利益相关者能够理解模型如何进行预测。此外，在部署模型后，我们依赖可视化进行持续监控，寻找任何性能下降的迹象，如模型漂移。

pandas 库提供了内置的绘图功能，用于可视化存储在 DataFrame 或 Series 数据结构中的数据。在后台，这些可视化由 Matplotlib 库支持，后者也是默认选项。

pandas 库提供了许多方便的绘图方法。简单调用 DataFrame.plot() 或 Series.plot() 默认将生成折线图。你可以通过两种方式更改图表的类型：

• 使用plot方法中的kind参数，如.plot(kind=)，通过替换<charttype>为图表类型来指定绘图类型。例如，.plot(kind=“hist”)将绘制一个直方图，而.plot(kind=“bar”)将绘制柱状图。

• 或者，你可以扩展plot方法。这可以通过链式调用特定的绘图函数来实现，例如.hist()或.scatter()，比如使用.plot.hist()或.plot.line()。

本配方将使用标准的 pandas .plot() 方法，并支持 Matplotlib 后端。

准备工作

你可以从 GitHub 仓库下载所需的 Jupyter 笔记本和数据集。请参考本章节的技术要求部分。

你将使用 Microsoft、Apple 和 IBM 的Closing Price Stock数据集，数据文件为closing_price.csv。

如何操作…

在这个配方中，你将探索如何绘制时间序列数据、改变主题、生成子图以及自定义输出的可视化效果：

1. 在 pandas 中绘图可以通过简单地在 DataFrame 或 Series 名称后添加.plot()来实现：

closing_price.plot();

这将生成一张线性图，它是kind参数的默认选项，类似于.plot(kind=“line”)：

图 9.1：使用 pandas 绘制的多线时间序列图

你可以通过添加标题、更新轴标签以及自定义 x 和 y 轴刻度等方式进一步自定义图表。例如，添加标题和更新 y 轴标签，可以使用title和ylabel参数，如下所示：

start_date = ‘2019’
end_date = ‘2021’
closing_price.plot(
          title=f’Closing Prices from {start_date} – {end_date}’,
          ylabel= ‘Closing Price’);

1. 如果你想看看价格之间的波动（上下浮动），一种简单的方法是标准化数据。为此，只需将每个股票的股价除以第一天的价格（第一行）。这将使所有股票具有相同的起点：

closing_price_n = closing_price.div(closing_price.iloc[0])
closing_price_n.plot(
            title=f’Closing Prices from {start_date} – {end_date}’,
            ylabel= ‘Normalized Closing Price’);

这将生成如下图表：

图 9.2：使用简单的标准化技术，使得价格波动更加直观易比对

从输出中可以观察到，现在这些线条的起点（原点）已经设置为1。图表显示了时间序列图中价格之间的偏差：

closing_price_n.head()

请注意，输出表格的第一行已设置为1.0：

图 9.3：标准化时间序列输出，所有序列的起点为 1

1. 此外，Matplotlib 允许你更改图表的样式。为此，你可以使用style.use()函数。你可以从现有模板中指定一个样式名称，或者使用自定义样式。例如，以下代码展示了如何从default样式切换到ggplot样式：

plt.style.use('ggplot')
closing_price_n.plot(
          title=f'Closing Prices from {start_date} - {end_date}',
          ylabel= 'Normalized Closing Price');

前面的代码应该生成相同数据内容的图表，但样式有所不同。

图 9.4：使用来自 Matplotlib 的 ggplot 风格

ggplot 风格的灵感来自于 R 中的 ggplot2 包。

您可以探索其他吸引人的样式：fivethirtyeight，灵感来自 fivethirtyeight.com，dark_background，dark-background，和 tableau-colorblind10。

要查看可用的样式表的完整列表，您可以参考 Matplotlib 文档：matplotlib.org/stable/gallery/style_sheets/style_sheets_reference.html

如果您想恢复到原始主题，可以指定 plt.style.use("default")。

1. 您可以将您的图表保存为 jpeg、png、svg 或其他文件类型。例如，您可以使用 .savefig() 方法将文件保存为 plot_1.jpg 文件，并指定更高的 dpi 分辨率以保证打印质量。默认的 dpi 值为 100：

plot = closing_price_n.plot(
            title=f'Closing Prices from {start_date} - {end_date}',
            ylabel= 'Norm. Price')
plot.get_figure().savefig('plot_1.jpg', dpi=300)

图表应作为 plot_1.jpg 图像文件保存在您的本地目录中。

它是如何工作的……

pandas 和 Matplotlib 库之间的协作非常好，双方有着集成并在 pandas 中增加更多绘图功能的雄心。

在 pandas 中，您可以通过向 kind 参数提供一个值来使用多种绘图样式。例如，您可以指定以下内容：

• line 用于常用于展示时间序列的折线图

• bar 或 barh（水平）用于条形图

• hist 用于直方图

• box 用于箱形图

• kde 或 density 用于核密度估计图

• area 用于面积图

• pie 用于饼图

• scatter 用于散点图

• hexbin 用于六边形图

还有更多……

如前一节所示，我们将时间序列中的所有三列绘制在一个图中（同一图中的三条线图）。如果您希望每个符号（列）单独绘制呢？

这可以通过简单地将 subplots 参数设置为 True 来完成：

closing_price.plot(
           subplots=True,
           title=f'Closing Prices from {start_date} - {end_date}');

上述代码将在 DataFrame 中为每一列生成一个子图。使用 closing_price DataFrame，这将生成三个子图。

图 9.5：使用 pandas 子图功能

另见

若要了解更多有关 pandas 绘图和可视化的功能，请访问官方文档：pandas.pydata.org/pandas-docs/stable/user_guide/visualization.html

使用 hvPlot 绘制时间序列数据并进行交互式可视化

交互式可视化比静态图像更高效地分析数据。简单的交互，例如放大缩小或切片操作，可以揭示出更多有价值的见解，供进一步调查。

在本食谱中，我们将探索 hvPlot 库来创建交互式可视化。HvPlot 提供了一个高级 API 用于数据可视化，并且与多种数据源无缝集成，包括 pandas、Xarray、Dask、Polars、NetworkX、Streamlit 和 GeoPandas。利用 hvPlot 和 pandas 渲染交互式可视化需要最少的工作量，只需对原始代码做少量修改，就能创建动态可视化。我们将使用 'closing_price.csv' 数据集来探索该库在本食谱中的应用。

准备工作

你可以从 GitHub 仓库下载所需的 Jupyter 笔记本和数据集。请参阅本章的 技术要求 部分。

如何操作…

1. 首先导入所需的库。注意，hvPlot 有一个 pandas 扩展，使得使用更为方便。这将允许你使用与前面示例相同的语法：

import pandas as pd
import hvplot.pandas
import hvplot as hv
closing_price_n = closing_price.div(closing_price.iloc[0])

使用 pandas 绘图时，你会使用 .plot() 方法，例如，closing_price_n.plot()。类似地，hvPlot 允许你通过将 .plot() 替换为 .hvplot() 来渲染交互式图表。如果你有内容密集的图表，这非常有用。你可以缩放到图表的特定部分，然后使用平移功能移动到图表的不同部分：

closing_price_n.hvplot(
    title=f'Closing Prices from {start_date} - {end_date}')

通过将 .plot 替换为 .hvplot，你可以获得一个带有悬停效果的交互式可视化：

图 9.6：hvPlot 交互式可视化

同样的结果可以通过简单地切换 pandas 绘图后端来实现。默认后端是 matplotlib。要切换到 hvPlot，只需更新 backend='hvplot'：

closing_price_n.plot(
    backend='hvplot',
    title=f'Closing Prices from {start_date} - {end_date}'
)

这应该会生成与图 9.6相同的图表。

注意右侧的小部件栏，其中包含一组交互模式，包括平移、框选缩放、滚轮缩放、保存、重置和悬停。

图 9.7：带有六种交互模式的小部件栏

1. 你可以将每个时间序列按符号（列）分割成独立的图表。例如，将数据分为三列，每列对应一个符号（或股票代码）：MSFT、AAPL 和 IBM。子图可以通过指定 subplots=True 来完成：

closing_price.hvplot(width=300, subplots=True)

这应该会生成每列一个子图：

图 9.8：hvPlot 子图示例

你可以使用 .cols() 方法来更精确地控制布局。该方法允许你控制每行显示的图表数量。例如，.cols(1) 表示每行一个图表，而 .cols(2) 表示每行两个图表：

closing_price.hvplot(width=300, subplots=True).cols(2)

这应该会生成一个图表，第一行有两个子图，第二行有第三个子图，如下所示：

图 9.9：使用 .col(2) 每行两个列的示例 hvPlot

请记住，.cols() 方法仅在 subplots 参数设置为 True 时有效，否则会导致错误。

工作原理…

由于 pandas 被广泛使用，许多库现在都支持将 pandas DataFrame 和 Series 作为输入。此外，Matplotlib 和 hvPlot 之间的集成简化了与 pandas 一起使用的绘图引擎的更换过程。

HvPlot 提供了几种方便的选项来绘制你的 DataFrame：你可以轻松切换后端，使用 DataFrame.hvplot() 扩展 pandas 功能，或利用 hvPlot 的原生 API 进行更高级的可视化。

还有更多内容…

hvPlot 允许你使用两个算术运算符，+ 和 *，来配置图表的布局。

加号 (+) 允许你将两个图表并排显示，而乘号 (*) 则使你能够合并图表（将一个图表与另一个合并）。在以下示例中，我们将两个图表相加，使它们在同一行上并排显示：

(closing_price_n['AAPL'].hvplot(width=400) +
 closing_price_n['MSFT'].hvplot(width=400))

这应该生成如下图所示的结果：

图 9.10：使用加法运算符并排显示两个图表

请注意，这两个图表将共享同一个小部件条。如果你对一个图表进行筛选或缩放，另一个图表将应用相同的操作。

现在，让我们看看如何通过乘法将两个图表合并成一个：

(closing_price_n['AAPL'].hvplot(width=500, height=300) *
 closing_price_n['MSFT'].hvplot()).opts(legend_position='top_left')

上述代码应该生成一个结合了 AAPL 和 MSFT 的图表：

图 9.11：使用乘法运算符将两个图表合并为一个

最后，为了创建子组（类似于 "group by" 操作），其中每个组由不同的颜色表示，你可以使用下方演示的 by 参数：

closing_price['AAPL'].hvplot.line(by=['index.year'])

这段代码生成了一个如预期的折线图，并按年份分段（分组显示），如图 9.12 所示：

图 9.12：按子组（按年份）绘制的折线图。

鉴于我们的数据涵盖了三年，你将在图表中看到三种不同的颜色，每种颜色对应不同的年份，如图例所示。

另请参阅

有关 hvPlot 的更多信息，请访问其官方网站：hvplot.holoviz.org/。

时间序列数据的分解

在进行时间序列分析时，一个关键目标通常是预测，即构建一个能够做出未来预测的模型。在开始建模过程之前，至关重要的一步是提取时间序列的各个组成部分进行分析。这个步骤对整个建模过程中的决策至关重要。

一个时间序列通常由三个主要组成部分构成：趋势、季节性和残差随机过程。对于需要时间序列平稳的统计模型，可能需要估计并随之去除时间序列中的趋势和季节性成分。时间序列分解的技术和库通常提供趋势、季节性和残差随机过程的可视化表示和识别。

趋势成分反映了时间序列的长期方向，可能是向上、向下或水平的。例如，销售数据的时间序列可能显示出向上的趋势，表明销售随着时间的推移在增加。季节性是指在特定间隔内重复出现的模式，例如每年圣诞节前后的销售增长，这是一个随着假日季节的临近每年都会出现的模式。残差随机过程表示在去除趋势和季节性后，时间序列中剩余的部分，包含了数据中无法解释的变动性。

时间序列的分解是将其分离为三个组成部分的过程，并将趋势和季节性成分作为各自的模型来估计。分解后的组成部分可以根据它们之间的交互性质进行加法或乘法建模。

当你使用加法模型时，可以通过将所有三个组成部分相加来重建原始时间序列：

当季节性变化不随时间变化时，加法分解模型是合理的。另一方面，如果时间序列可以通过将这三部分相乘来重建，则使用乘法模型：

当季节性变化随时间波动时，使用乘法模型是合适的。

此外，你可以将这些成分分为可预测与不可预测的成分。可预测成分是稳定的、重复的模式，可以被捕捉和建模。季节性和趋势就是其中的例子。另一方面，每个时间序列都有一个不可预测的成分，它表现出不规则性，通常称为噪声，但在分解的上下文中被称为残差。

在本教程中，你将探索使用seasonal_decompose、季节-趋势分解（LOESS）（STL）和hp_filter方法来分解时间序列，这些方法在statsmodels库中可用。

准备工作

你可以从 GitHub 仓库下载所需的 Jupyter 笔记本和数据集。请参考本章节中的技术要求部分。

如何操作…

你将探索在 statsmodels 库中提供的两种方法：seasonal_decompose 和 STL。

使用 seasonal_decompose

seasonal_decompose函数依赖于移动平均法来分解时间序列。你将使用技术要求部分中的 CO2 和航空乘客数据集。

1. 导入所需的库并设置rcParams，使可视化图表足够大。通常，statsmodels 生成的图表较小。你可以通过调整rcParams中的figure.figsize来修复这个问题，使本食谱中的所有图表都应用相同的大小：

from statsmodels.tsa.seasonal import seasonal_decompose, STL
plt.rcParams["figure.figsize"] = [10, 5]

这将使所有图表的大小一致：宽度为 10 英寸，高度为 3 英寸（W x H）。

你可以应用诸如灰度主题之类的样式

plt.style.use('grayscale')

1. 你可以使用seasonal_decompose()函数分解这两个数据集。但在此之前，你应该绘制你的时间序列，以了解季节性是否表现出乘法或加法特征：

co2_df.plot(title=co2.TITLE);

这应该会显示一张显示 1960 年到 2000 年每周二氧化碳水平的线图，单位为百万分之一（ppm）。使用.plot()方法时，默认的图表类型是线图，参数为kind="line"。关于 pandas 绘图功能的更多信息，请参考使用 pandas 绘制时间序列数据的食谱。

图 9.13：显示上升趋势和恒定季节性变化的 CO2 数据集

co2_df数据展示了一个长期的线性趋势（向上），并且有一个以恒定速率重复的季节性模式（季节性变化）。这表明 CO2 数据集是一个加法模型。

同样，你可以探索airp_df数据框，观察航空乘客数据集中的季节性是否表现为乘法或加法行为：

airp_df['value'].plot(title=air_passengers['title']);

这应该生成一张显示 1949 年至 1960 年每月乘客数量的线图：

图 9.14：显示趋势和逐渐增加的季节性变化的航空乘客数据集

airp_df数据展示了一个长期的线性趋势和季节性变化（向上）。然而，季节性波动似乎也在增加，这表明是一个乘法模型。

1. 对这两个数据集使用seasonal_decompose。对于 CO2 数据，使用加法模型，而对于航空乘客数据，使用乘法模型：

co2_decomposed = seasonal_decompose(co2_df,model='additive')
air_decomposed = seasonal_decompose(airp_df,model='multiplicative')

co2_decomposed和air_decomposed都可以访问几种方法，包括.trend、.seasonal和.resid。你可以通过使用.plot()方法绘制所有三个组成部分：

air_decomposed.plot();

以下是结果图：

图 9.15：航空乘客的乘法分解为趋势、季节性和残差

让我们将结果图分解成四个部分：

1. 这是我们正在分解的原始观察数据。

2. 趋势组件显示出上升的方向。趋势指示是否存在正向（增加或上升）、负向（减少或下降）或恒定（没有趋势或水平）长期运动。

3. 季节性组件显示季节性效应，表现为高低交替的重复模式。

4. 最后，残差（有时称为噪声）组件显示去除趋势和季节性后的数据中的随机变化。在这种情况下，使用了乘法模型。

同样，你可以绘制 CO2 数据集的分解图：

co2_decomposed.plot();

这应该会生成以下图形：

图 9.16：CO2 加性分解为趋势、季节性和残差

1. 在重建时间序列时，例如，在乘法模型中，你将会将三个组件相乘。为了演示这一概念，使用air_decomposed，这是DecomposeResult类的一个实例。该类提供了seasonal、trend和resid属性以及.plot()方法。

在下面的代码中，你可以将这些组件相乘来重建时间序列：

(air_decomposed.trend *
 air_decomposed.seasonal *
 air_decomposed.resid).plot() ;

它输出如下图：

图 9.17：重建航空乘客时间序列数据集

使用 STL

另一种在statsmodels中的分解选项是STL。STL 代表使用 LOESS 的季节性-趋势分解，这是一种更先进的分解技术。在statsmodels中，STL类比seasonal_decompose函数需要更多的参数。你将使用的另外两个参数是seasonal和robust。seasonal参数用于季节性平滑器，并且只能接受大于或等于 7 的奇数整数值。同样，STL函数有一个趋势平滑器（trend参数）。

第二个参数是robust，它接受一个布尔值（True 或 False）。设置robust=True有助于消除异常值对季节性和趋势组件的影响。

1. 你将使用STL来分解co2_df数据框：

co2_stl = STL(
    co2_df,
    seasonal=13,
    robust=True).fit()
co2_stl.plot(); plt.show()

这应该会生成与seasonal_decompose函数类似的子图，显示趋势、季节性和残差：

图 9.18：使用 STL 分解 CO2 数据集

比较图 9.16中的输出与图 9.18中的输出。

请注意，当你使用STL时，你提供了seasonal=13，因为数据具有年度季节性效应。季节性参数仅接受大于或等于 7 的奇数整数。

1. 季节性分解和 STL 都生成DecomposeResult类的实例，可以直接访问残差。你可以比较seasonal_decompose和 STL 的残差，如下所示：

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 4))
co2_stl.resid.plot(ax=ax1, title='Residual Plot from Seasonal Decomposition')
co2_decomposed.resid.plot(ax=ax2, title='Residual Plot from STL');

这应该会生成如下图，包含两个子图

图 9.19：比较季节分解与 STL 的残差图

你会注意到，残差图看起来不同，这表明两种方法使用不同的机制捕获了相似的信息。

它是如何工作的……

你使用了两种不同的时间序列分解方法。这两种方法都将时间序列分解为趋势、季节性和残差成分。

STL 类使用 LOESS 季节性平滑器，代表的是 局部加权散点图平滑。与 seasonal_decompose 相比，STL 在测量非线性关系方面更加稳健。另一方面，STL 假设成分是加性组合的，因此你不需要像 seasonal_decompose 那样指定模型。两种方法都可以从时间序列中提取季节性，以更好地观察数据的整体趋势。

seasonal_decompose 函数执行以下 简化 逻辑：

• 对时间序列数据进行平滑处理，以观察趋势。这是通过应用卷积滤波器来估计趋势。

• 一旦趋势被估计出来，它就会从时间序列中被移除（去趋势）。

• 剩下的去趋势时间序列会根据每个周期或季节组进行平均。季节性平均会将时间序列中对应每个季节的所有数据值取平均（例如，我们取每年的每个一月并对其进行平均）。

• 一旦季节性成分被估计出来，它就会被移除，剩下的就是残差。

STL 函数执行以下 简化 逻辑：

• 类似于 seasonal_decompose，时间序列通过平滑处理来估计趋势，但在 STL 中，这一过程是通过 LOESS 平滑来实现的。

• 一旦趋势被估计出来，它就会从时间序列中被移除（去趋势）。

• 对于季节性成分，STL 对去趋势后的数据应用了 Loess 平滑，但每个季节单独处理。

• 一旦季节性成分被估计出来，它就会被移除，剩下的就是残差。

还有更多……

霍德里克-普雷斯科特滤波器 是一种平滑滤波器，用于将短期波动（周期变化）与长期趋势分离。它在 statsmodels 库中实现为 hp_filter。

回顾一下，STL 和 seasonal_decompose 返回了三个成分（趋势、季节性和残差）。另一方面，hp_filter 只返回两个成分：周期成分和趋势成分。

首先从 statsmodels 库中导入 hpfilter 函数：

from statsmodels.tsa.filters.hp_filter import hpfilter

要将你的时间序列分解为趋势和周期成分，只需将你的时间序列 DataFrame 提供给 hpfilter 函数，如下所示：

co2_cyclic, co2_trend = hpfilter(co2_df)

hpfilter 函数返回两个 pandas Series：第一个 Series 是周期成分，第二个 Series 是趋势成分。将 co2_cyclic 和 co2_trend 并排绘制，以便更好地理解霍德里克-普雷斯科特滤波器从数据中提取了哪些信息：

fig, ax = plt.subplots(1,2, figsize=(16, 4))
co2_cyclic.plot(ax=ax[0], title='CO2 Cyclic Component')
co2_trend.plot(ax=ax[1], title='CO2 Trend Component');

这应该会在同一行中产生两个子图（并排显示），如所示：

图 9.20：使用赫德里克-普雷斯科特滤波器的周期性和趋势成分

请注意，hpfilter 得到的两个成分是加性的。换句话说，要重构原始时间序列，你需要将 co2_cyclic 和 co2_trend 相加。

(co2_cyclic + co2_trend).plot();

图 9.21：从 hpfilter 函数给出的趋势和周期成分重构 CO2 数据集

你可以将图 9.21 中从趋势和周期性成分重构的 CO2 图与图 9.13 中的原始 CO2 图进行比较。

另见

• 要了解更多关于hpfilter()的信息，请访问官方文档页面：www.statsmodels.org/dev/generated/statsmodels.tsa.filters.hp_filter.hpfilter.html#statsmodels.tsa.filters.hp_filter.hpfilter。

• 要了解更多关于seasonal_decompose()的信息，请访问官方文档页面：www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.seasonal_decompose.html。

• 要了解更多关于STL()的信息，请访问官方文档页面：www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.STL.html#statsmodels.tsa.seasonal.STL。

检测时间序列的平稳性

一些时间序列预测技术假设时间序列过程是平稳的。因此，确定你正在处理的时间序列（无论是观察到的时间序列还是你所拥有的实现）是否来源于平稳或非平稳过程是至关重要的。

一个平稳的时间序列表明特定的统计性质随时间不会发生变化，保持稳定，这使得建模和预测过程更加容易。相反，非平稳过程由于其动态特性和随时间的变化（例如，存在趋势或季节性）更难建模。

定义平稳性的方式有多种；有些方法严格，可能在现实数据中无法观察到，称为强平稳性。与之相反，其他定义则在标准上更为宽松，能够在现实数据中观察到（或通过变换得到），称为弱平稳性。

在本教程中，出于实际考虑，平稳性一词指的是“弱”平稳性，即定义为具有恒定均值mu（），恒定方差sigma squared（）和一致的协方差（或自相关）在相同距离的周期（lags）之间的时间序列。均值和方差作为常数有助于简化建模，因为你不需要将它们作为时间的函数来求解。

一般来说，具有趋势或季节性的时间序列可以认为是非平稳的。通常，通过图形直观识别趋势或季节性有助于判断时间序列是否平稳。在这种情况下，一个简单的折线图就足够了。但在本教程中，你将探索统计检验，帮助你从数值上识别时间序列是平稳还是非平稳。你将探索平稳性检验及使时间序列平稳的技术。

你将使用 statsmodels 库探索两种统计检验：扩展迪基-福勒（ADF）检验和克维亚特科夫-菲利普斯-施密特-辛（KPSS）检验。ADF 和 KPSS 都用于检验单变量时间序列过程中的单位根。需要注意的是，单位根只是时间序列非平稳的原因之一，但通常单位根的存在表明时间序列是非平稳的。

ADF 和 KPSS 都基于线性回归，并且是统计假设检验的一种。例如，ADF 的零假设表示时间序列中存在单位根，因此是非平稳的。另一方面，KPSS 的零假设相反，假设时间序列是平稳的。因此，你需要根据检验结果来判断是否可以拒绝零假设或未能拒绝零假设。通常，可以依赖返回的 p 值来决定是否拒绝或未能拒绝零假设。

请记住，ADF 和 KPSS 检验的解释不同，因为它们的零假设相反。如果 p 值小于显著性水平（通常为 0.05），则可以拒绝零假设，说明时间序列没有单位根，可能是平稳的。如果 KPSS 检验中的 p 值小于显著性水平，则表示可以拒绝零假设，表明该序列是非平稳的。

准备工作

你可以从 GitHub 仓库下载所需的 Jupyter notebook 和数据集。请参阅本章节的技术要求部分。

在本教程中，你将使用 CO2 数据集，该数据集之前已作为 pandas DataFrame 在本章节的技术要求部分加载。

如何操作…

除了通过时间序列图形的视觉解释来判断平稳性之外，一种更为具体的方法是使用单位根检验，例如 ADF KPSS 检验。

在图 9.13中，你可以看到一个上升的趋势和一个重复的季节性模式（年度）。然而，当存在趋势或季节性（在这种情况下，两者都有）时，这会使时间序列非平稳。并不是总能通过视觉轻松识别平稳性或其缺乏，因此，你将依赖统计检验。

你将使用 statsmodels 库中的 adfuller 和 kpss 检验，并在理解它们有相反的原假设的前提下解释它们的结果：

1. 首先从 statsmodels 导入 adfuller 和 kpss 函数：

from statsmodels.tsa.stattools import adfuller, kpss

为了简化对检验结果的解释，创建一个函数以用户友好的方式输出结果。我们将该函数称为 print_results：

def print_results(output, test='adf'):   
    pval = output[1]
    test_score = output[0]
    lags = output[2]   
    decision = 'Non-Stationary'
    if test == 'adf':
        critical = output[4]
        if pval < 0.05:
            decision = 'Stationary'
    elif test=='kpss':
        critical = output[3]
        if pval >= 0.05:
            decision = 'Stationary'           
    output_dict = {
    'Test Statistic': test_score,
    'p-value': pval,
    'Numbers of lags': lags,
    'decision': decision
    }
    for key, value in critical.items():
        output_dict["Critical Value (%s)" % key] = value

    return pd.Series(output_dict, name=test)

该函数接受 adfuller 和 kpss 函数的输出，并返回一个添加了标签的字典。

1. 运行 kpss 和 adfuller 检验。对于这两个函数，使用默认的参数值：

adf_output = adfuller(co2_df)
kpss_output = kpss(co2_df)

1. 将两个输出传递给 print_results 函数，并将它们合并为 pandas DataFrame，便于比较：

pd.concat([
print_results(adf_output, 'adf'),
print_results(kpss_output, 'kpss')
], axis=1)

这应该生成以下 DataFrame：

图 9.22：来自 ADF 和 KPSS 单位根检验的结果输出

对于 ADF，p 值为 0.96，大于 0.05，因此你不能拒绝原假设，因此时间序列是非平稳的。对于 KPSS，p 值为 0.01，小于 0.05，因此你拒绝原假设，因此时间序列是非平稳的。

接下来，你将探索六种使时间序列平稳的方法，如变换和差分。所涵盖的方法有一阶差分、二阶差分、减去移动平均、对数变换、分解和霍德里克-普雷斯科特滤波器。

本质上，通过去除趋势（去趋势化）和季节性效应可以实现平稳性。对于每个变换，你将运行平稳性检验，并比较不同方法的结果。为了简化解释和比较，你将创建两个函数：

• check_stationarity 接受一个 DataFrame，执行 KPSS 和 ADF 检验，并返回结果。

• plot_comparison 接受一个方法列表并比较它们的图表。该函数接受一个 plot_type 参数，因此你可以探索折线图和直方图。该函数调用 check_stationarity 函数，以捕捉子图标题的结果。

创建 check_stationarity 函数，这是对之前使用的 print_results 函数的简化重写：

 def check_stationarity(df):
    kps = kpss(df)
    adf = adfuller(df)

    kpss_pv, adf_pv = kps[1], adf[1]
    kpssh, adfh = 'Stationary', 'Non-stationary'

    if adf_pv < 0.05:
        # Reject ADF Null Hypothesis
        adfh = 'Stationary'
    if kpss_pv < 0.05:
        # Reject KPSS Null Hypothesis
        kpssh = 'Non-stationary'
    return (kpssh, adfh) 

创建 plot_comparison 函数：

def plot_comparison(methods, plot_type='line'):
    n = len(methods) // 2
    fig, ax = plt.subplots(n,2, sharex=True, figsize=(20,10))
    for i, method in enumerate(methods):
        method.dropna(inplace=True)
        name = [n for n in globals() if globals()[n] is method]
        v, r = i // 2, i % 2
        kpss_s, adf_s = check_stationarity(method)
        method.plot(kind=plot_type,
                    ax=ax[v,r],
                    legend=False,
                    title=f'{name[0]} --> KPSS: {kpss_s}, ADF {adf_s}')
        ax[v,r].title.set_size(20)
        method.rolling(52).mean().plot(ax=ax[v,r], legend=False) 

让我们实现一些使时间序列平稳或提取平稳成分的方法。然后，将这些方法合并成一个 Python 列表：

1. 一阶差分：也称为去趋势化，它通过将时间 t 的观测值减去时间 t-1 的前一个观测值来计算（

   

   在 pandas 中，可以使用.diff()函数来实现，这个函数的默认参数是period=1。注意，差分后的数据将比原始数据少一个数据点（行），因此需要使用.dropna()方法：

first_order_diff = co2_df.diff().dropna()

1. 二阶差分：如果存在季节性或一阶差分不足时，这非常有用。这本质上是进行两次差分——第一次差分去除趋势，第二次差分去除季节性趋势：

differencing_twice = co2_df.diff(52).diff().dropna()

1. 从时间序列中减去移动平均（滚动窗口），使用DataFrame.rolling(window=52).mean()，因为这是每周数据：

rolling_mean = co2_df.rolling(window=52).mean()
rolling_mean = co2_df - rolling

1. 对数变换，使用np.log()，是一种常见的技术，用于稳定时间序列的方差，有时足以使时间序列平稳。简单来说，它所做的就是用每个观测值的对数值替代原始值：

log_transform = np.log(co2_df)

1. 使用时间序列分解方法去除趋势和季节性成分，例如seasonal_decompose。从图 9.13来看，似乎该过程是加法性的。这是seasonal_decompose的默认参数，因此这里无需进行任何更改：

decomp = seasonal_decompose(co2_df)
seasonal_decomp = decomp.resid

让我们也添加 STL 作为另一种分解方法

co2_stl = STL(co2_df, seasonal=13,robust=True).fit()
stl_decomp = co2_stl.resid

1. 使用霍德里克-普雷斯科特滤波器去除趋势成分，例如使用hp_filter：

hp_cyclic, hp_trend = hpfilter(co2_df)

现在，让我们将这些方法组合成一个 Python 列表，然后将该列表传递给plot_comparison函数：

methods = [first_ord_diff, second_ord_diff,
           diseasonalize, rolling_mean,
           log_transform,  seasonal_decomp,
           stl_decomp, hp_cyclic]
plot_comparison(methods)

这应显示如图所示的 4 x 2 子图：

图 9.23：绘制不同方法使 CO2 时间序列平稳

通常，你不想对时间序列进行过度差分，因为一些研究表明，基于过度差分数据的模型准确性较低。例如，first_order_diff已经使时间序列平稳，因此无需进一步进行差分。换句话说，differencing_twice是不必要的。此外，注意到log_transform仍然是非平稳的。

注意中间线代表时间序列的平均值（移动平均）。对于平稳时间序列，均值应保持恒定，且更像是一条水平直线。

它是如何工作的…

平稳性是时间序列预测中的一个重要概念，尤其在处理金融或经济数据时尤为相关。我们之前将平稳性（弱平稳）定义为均值恒定、方差恒定、协方差一致。

如果时间序列是平稳的，则均值被认为是稳定且恒定的。换句话说，存在一种平衡，尽管值可能会偏离均值（高于或低于），但最终它总会回到均值。一些交易策略依赖于这一核心假设，正式称为均值回归策略。

statsmodels库提供了多个平稳性检验方法，例如adfuller和kpss函数。两者都被认为是单位根检验，用于确定是否需要差分或其他转换策略来使时间序列平稳。

记住，ADF 和 KPSS 检验基于不同的原假设。例如，adfuller和kpss的原假设是相反的。因此，你用来拒绝（或无法拒绝）原假设的 p 值在这两者之间的解读会有所不同。

在图 9.22中，测试返回了额外的信息，具体包括以下内容：

• 检验统计量值为 ADF 的 0.046 和 KPSS 的 8.18，均高于 1%的临界值阈值。这表明时间序列是非平稳的，确认你不能拒绝原假设。ADF 的临界值来自 Dickey-Fuller 表格。幸运的是，你不必参考 Dickey-Fuller 表格，因为所有提供 ADF 检验的统计软件/库都会在内部使用该表格。KPSS 也是如此。

• p 值结果与检验统计量相关。通常，当 p 值小于 0.05（5%）时，你可以拒绝原假设。再次强调，在使用 ADF、KPSS 或其他平稳性检验时，确保理解原假设，以便准确解读结果。

• 滞后期数表示检验中自回归过程使用的滞后期数（ADF 和 KPSS）。在这两个检验中，使用了 27 个滞后期。由于我们的 CO2 数据是按周记录的，一个滞后期表示 1 周。因此，27 个滞后期代表我们数据中的 27 周。

• 使用的观测值数量是数据点的数量，排除了滞后期数。

• 最大化的信息准则基于autolag参数。默认值是autolag="aic"，对应赤池信息准则（AIC）。其他可接受的autolag参数值有bic（对应贝叶斯信息准则（BIC））和t-stat。

你还探索了一些时间序列去趋势（移除趋势）的技术，以使其平稳。例如，你使用了一级差分、分解和对数变换来去除趋势的影响。去趋势稳定了时间序列的均值，有时这就是使其平稳所需的全部。当你决定去趋势时，本质上是移除一个干扰因素，这样你就能专注于那些不那么明显的潜在模式。因此，你可以构建一个模型来捕捉这些模式，而不会被长期趋势（向上或向下的波动）所掩盖。一个例子是一级差分方法。

然而，在存在季节性模式的情况下，您还需要去除季节性效应，这可以通过季节性差分来完成。这是在去趋势的第一阶差分之外进行的；因此，它可以称为二阶差分、双重差分，或者称为“差分两次”，因为您首先使用差分去除趋势效应，然后再次去除季节性。这假设季节性差分不足以使时间序列平稳。您的目标是使用最少的差分，避免过度差分。您很少需要超过两次差分。

还有更多内容…

在本食谱的介绍部分，我们提到过，ADF 和 KPSS 都使用线性回归。更具体地说，普通最小二乘法（OLS）回归用于计算模型的系数。要查看 ADF 的 OLS 结果，您需要使用store参数并将其设置为True：

adf_result = adfuller(first_order_diff, store=True)

上面的代码将返回一个包含测试结果的元组。回归摘要将作为最后一项附加到元组中。元组中应该有四项：第一项，adf_result[0]，包含t 统计量；第二项，adf_result[1]，包含p 值；第三项，adf_result[2]，包含 1%、5%和 10%区间的临界值。最后一项，adf_result[3]，包含ResultStore对象。您也可以使用adf_result[-1]来访问最后一项，如以下代码所示：

adf_result[-1].resols.summary()

ResultStore对象让您可以访问.resols，其中包含.summary()方法。这将生成如下输出：

图 9.24：ADF OLS 回归摘要及前五个滞后项及其系数

另见

要了解更多关于平稳性和去趋势的内容，请访问官方 statsmodels 页面：www.statsmodels.org/dev/examples/notebooks/generated/stationarity_detrending_adf_kpss.html。

应用幂次变换

时间序列数据可能很复杂，数据中嵌入着您需要理解和分析的重要信息，以便确定构建模型的最佳方法。例如，您已经探索了时间序列分解，理解了趋势和季节性的影响，并测试了平稳性。在前面的食谱中，检测时间序列的平稳性，您研究了将数据从非平稳转变为平稳的技术。这包括去趋势的概念，旨在稳定时间序列的均值。

根据你所进行的模型和分析，可能需要对观察数据集或模型的残差进行额外假设检验。例如，检验同方差性（也拼作 homoscedasticity）和正态性。同方差性意味着方差在时间上是稳定的。更具体地说，它是残差的方差。当方差不是常数，而是随时间变化时，我们称之为异方差性（也拼作 heteroscedasticity）。另一个需要检验的假设是正态性；该特定观察值是否来自正态（高斯）分布？有时，你可能还需要检查残差的正态性，这通常是模型诊断阶段的一部分。因此，了解特定模型或技术所做的假设非常重要，这样你才能确定使用哪些测试以及针对哪个数据集进行测试。如果不进行这些检查，你可能会得到一个有缺陷的模型或一个过于乐观或悲观的结果。

此外，在本教程中，你将学习Box-Cox 变换，它可以帮助你转换数据，以满足正态性和同方差性的要求。Box-Cox 变换的形式如下：

图 9.25：Box-Cox 变换

Box-Cox 变换依赖于一个参数，lambda（），并涵盖了对数变换和幂变换。如果为 0，则得到自然对数变换；否则，则为幂变换。该方法是尝试不同的值，然后测试正态性和同方差性。例如，SciPy 库有一个 boxcox 函数，你可以通过 lmbda 参数指定不同的值（有趣的是，在实现中这样拼写，因为 lambda 是 Python 的保留关键字）。如果将 lmbda 参数设置为 None，该函数将为你找到最佳的 lambda（）值。

准备工作

你可以从 GitHub 仓库下载所需的 Jupyter notebook 和数据集。请参考本章的技术要求部分。

在本教程中，你将使用空气乘客数据集，这个数据集已在本章技术要求部分作为 pandas DataFrame 加载。

你将使用 SciPy 和 statsmodels。

对于 pip 安装，请使用以下命令：

> pip install scipy

对于 conda 安装，请使用以下命令：

> conda install -c anaconda scipy

除了在技术要求部分中提到的准备工作外，你还需要导入在本教程中将使用到的常用库：

import numpy as np
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
import statsmodels.api as sm

为了让图表更大、更易读，可以使用以下命令设置固定的尺寸（20, 8）——宽度为 20 英寸，高度为 8 英寸：

plt.rcParams["figure.figsize"] = (20,8)

如何操作…

在本节中，你将扩展从上一节检测时间序列平稳性中学到的内容，并测试两个额外的假设：正态性和同方差性。

通常，平稳性是你需要关注的最关键假设，但熟悉额外的诊断技术将对你大有帮助。

有时，你可以通过图形来判断正态性和同方差性，例如，通过直方图或Q-Q 图。本节旨在教你如何在 Python 中编程执行这些诊断检验。此外，你还将了解White 检验和Breusch-Pagan Lagrange统计检验，用于同方差性。

对于正态性诊断，你将探索Shapiro-Wilk、D'Agostino-Pearson和Kolmogorov-Smirnov统计检验。总体而言，Shapiro-Wilk 倾向于表现最好，且能够处理更广泛的情况。

正态性检验

statsmodels 库和 SciPy 库有重叠的实现。例如，Kolmogorov-Smirnov 检验在 SciPy 中实现为ktest，而在 statsmodels 中实现为ktest_normal。在 SciPy 中，D'Agostino-Pearson 检验实现为normaltest，Shapiro-Wilk 检验实现为shapiro：

1. 首先导入 SciPy 和 statsmodels 库提供的正态性检验：

from scipy.stats import shapiro, kstest, normaltest
from statsmodels.stats.diagnostic import kstest_normal

1. 正态性诊断是一种基于原假设的统计检验，你需要确定是否可以接受或拒绝该假设。方便的是，以下你将实现的检验有相同的原假设。原假设声明数据是正态分布的；例如，当 p 值小于 0.05 时，你将拒绝原假设，表示时间序列不服从正态分布。让我们创建一个简单的函数is_normal()，它将根据 p 值返回Normal或Not Normal：

def is_normal(test, p_level=0.05):
    stat, pval = test
    return 'Normal' if pval > 0.05 else 'Not Normal'

运行每个检验以检查结果：

normal_args = (np.mean(co2_df),np.std(co2_df))
print(is_normal(shapiro(co2_df)))
print(is_normal(normaltest(co2_df)))
print(is_normal(normal_ad(co2_df)))
print(is_normal(kstest_normal(co2_df)))
print(is_normal(kstest(co2_df,
                   cdf='norm',
                   args=(np.mean(co2_df), np.std(co2_df)))))
>>
Not Normal
Not Normal
Not Normal
Not Normal
Not Normal

检验结果显示数据不来自正态分布。你不需要运行那么多检验。例如，shapiro检验是一个非常常见且受欢迎的检验，你可以依赖它。通常，像任何统计检验一样，你需要阅读关于该检验实现的文档，以了解检验的具体内容。更具体地，你需要了解检验背后的原假设，以便决定是否可以拒绝原假设或未能拒绝原假设。

1. 有时，你可能需要在模型评估和诊断过程中测试正态性。例如，你可以评估残差（定义为实际值和预测值之间的差异）是否符合正态分布。在第十章，使用统计方法构建单变量时间序列模型中，你将探索使用自回归和移动平均模型构建预测模型。目前，你将运行一个简单的自回归（AR(1)）模型，演示如何将正态性检验应用于模型的残差：

from statsmodels.tsa.api import AutoReg
model = AutoReg(co2_df.dropna(), lags=1).fit()

你可以对残差运行shapiro检验。要访问残差，你可以使用.resid属性，如model.resid。这是你在第十章中构建的许多模型中的常见方法，使用统计方法构建单变量时间序列模型：

print(is_normal(shapiro(model.resid)))
>>
'Not Normal'

输出结果表明残差不是正态分布的。残差不正态分布这一事实本身不足以确定模型的有效性或潜在改进，但结合其他测试，这应该有助于你评估模型的优劣。这是你将在下一章进一步探讨的话题。

测试同方差性

你将测试模型残差的方差稳定性。这将是与之前正态性测试中使用的相同的 AR(1)模型：

1. 让我们首先导入本章所需的方法：

from statsmodels.stats.api import (het_breuschpagan,
                                   het_white)

1. 你将对模型的残差进行同方差性测试。正如前面所述，理解这些统计测试背后的假设非常重要。原假设表示数据是同方差的，适用于这两个测试。例如，如果 p 值小于 0.05，你将拒绝原假设，这意味着时间序列是异方差的。

让我们创建一个小函数，命名为het_test(model, test)，该函数接受一个模型和一个测试函数，并根据 p 值返回Heteroskedastic或Homoskedastic，以确定是否接受或拒绝原假设：

def het_test(model, test=het_breuschpagan):
    lm, lm_pvalue, fvalue, f_pvalue = (
        het_breuschpagan(model.resid,
                         sm.add_constant(
                             model.fittedvalues)
                        ))

    return "Heteroskedastic" if f_pvalue < 0.05 else "Homoskedastic"

1. 从 Breusch-Pagan 拉格朗日乘数检验开始诊断残差。在 statsmodels 中，你将使用het_breuschpagan函数，该函数接受resid（模型的残差）和exog_het（提供与残差异方差相关的原始数据，即解释变量）：

het_test(model, test=het_breuschpagan)
>> 'Homoskedastic'

这个结果表明残差是同方差的，具有恒定方差（稳定性）。

1. 一个非常相似的测试是怀特（White）的拉格朗日乘数检验。在 statsmodels 中，你将使用het_white函数，它有两个参数，与你使用het_breuschpagan时相同：

het_test(model, test=het_white)
>> 'Homoskedastic'

两个测试都表明自回归模型的残差具有恒定方差（同方差）。这两个测试都估计了辅助回归，并使用了残差的平方以及所有解释变量。

请记住，正态性和同方差性是你在诊断模型时可能需要对残差进行的一些测试。另一个重要的测试是自相关性测试，相关内容将在接下来的章节中讨论，时间序列数据中的自相关性测试。

应用 Box-Cox 变换

Box-Cox 变换可以是一个有用的工具，熟悉它是很有帮助的。Box-Cox 将一个非正态分布的数据集转换为正态分布的数据集。同时，它稳定了方差，使数据呈同方差性。为了更好地理解 Box-Cox 变换的效果，你将使用包含趋势和季节性的航空乘客数据集：

1. 从 SciPy 库中导入 boxcox 函数：

from scipy.stats import boxcox

1. 回顾本配方引言部分和 图 9.22，有一个 lambda 参数用于确定应用哪种变换（对数变换或幂变换）。使用默认的 lmbda 参数值为 None 的 boxcox 函数，只需提供数据集以满足所需的 x 参数：

xt, lmbda = boxcox(airp_df['passengers'])
xts = pd.Series(xt, index=airp_df.index)

通过不提供 lmbda 的值并将其保持为 None，该函数将找到最佳的 lambda () 值。根据本配方的引言，你会记得在 boxcox 实现中 lambda 被拼写为 lmbda。该函数返回两个值，xt 用于保存变换后的数据，lmda 用于保存找到的最佳 lambda 值。

直方图可以直观地展示变换的影响：

fig, ax = plt.subplots(1, 2, figsize=(16,5))
airp_df.hist(ax=ax[0])
ax[0].set_title('Original Time Series')
xts.hist(ax=ax[1])
ax[1].set_title('Box-Cox Transformed');

这应该生成以下两个图：

图 9.26：Box-Cox 变换及其对分布的影响

第二个直方图显示数据已被变换，总体分布发生了变化。将数据集作为时间序列图进行分析会很有意思。

绘制两个数据集，以比较变换前后的效果：

fig, ax = plt.subplots(1, 2, figsize=(16,5))
airp_df.plot(ax=ax[0])
ax[0].set_title('Original Time Series')
xts.plot(ax=ax[1])
ax[1].set_title('Box-Cox Transformed');

这应该生成以下两个图：

图 9.27：Box-Cox 变换及其对时间序列数据的总体影响

注意观察变换后的数据集，季节性效应看起来比之前更稳定。

1. 最后，构建两个简单的自回归模型，比较变换前后对残差的影响：

fig, ax = plt.subplots(1, 2, figsize=(16,5))
model_airp.resid.plot(ax=ax[0])
ax[0].set_title('Residuals Plot - Regular Time Series')
model_bx.resid.plot(ax=ax[1])
ax[1].set_title('Residual Plot - Box-Cox Transformed');

这应该生成以下两个图：

图 9.28：Box-Cox 变换及其对残差的影响

它是如何工作的…

Box-Cox 使我们能够将数据转换为正态分布且具有同方差性，它是一个包含对数变换和平方根变换等变换的幂变换家族的一部分。Box-Cox 是一种强大的变换方法，因为它支持根变换和对数变换，并且通过调整 lambda 值可以实现其他变换。

需要指出的是，boxcox 函数要求数据为正数。

还有更多…

AutoReg 模型有两个有用的方法：diagnostic_summary() 和 plot_diagnostics()。它们可以节省你编写额外代码的时间，测试模型残差的正态性、同方差性和自相关性。

以下代码展示了如何获取model_bx的诊断摘要：

print(model_bx.diagnostic_summary())

这应显示 Ljung-Box 自相关检验结果和模型残差的同方差性检验。

图 9.29：自相关的 diagnostic_summary

要获得视觉摘要，可以使用以下代码：

model_bx.plot_diagnostics(); plt.show()

.plot_diagnostics()函数将显示四个图表，您可以检查模型的残差。主要地，图表将显示残差是否从 Q-Q 图和直方图中呈现正态分布。此外，自相关函数图（ACF）将允许您检查自相关。您将在第十章的“绘制 ACF 和 PACF”一节中更详细地研究 ACF 图。

图 9.30：来自 plot_diagnostics()方法的输出

参见

要了解更多关于boxcox函数的信息，请访问官方 SciPy 文档：docs.scipy.org/doc/scipy/reference/generated/scipy.stats.boxcox.html。

测试时间序列数据中的自相关

自相关类似于统计学中的相关性（可以理解为高中学过的皮尔逊相关），用于衡量两个变量之间线性关系的强度，不同之处在于我们衡量的是滞后时间序列值之间的线性关系。换句话说，我们是在比较一个变量与其滞后的版本之间的关系。

在本节中，您将执行Ljung-Box 检验，检查是否存在直到指定滞后的自相关，以及它们是否显著偏离 0。Ljung-Box 检验的原假设表示前期滞后与当前期无关。换句话说，您正在检验自相关的不存在。

使用 statsmodels 中的acorr_ljungbox运行检验时，您需要提供一个滞后值。该检验将在所有滞后值（直到指定的最大滞后）上进行。

自相关测试是另一种有助于模型诊断的测试。如前面“应用幂变换”一节所讨论的那样，模型的残差需要进行假设检验。例如，当对残差进行自相关测试时，期望残差之间没有自相关。这确保了模型已经捕获了所有必要的信息。残差中的自相关可能表示模型未能捕捉到关键信息，需进行评估。

准备工作

您可以从 GitHub 仓库下载所需的 Jupyter 笔记本和数据集。请参考本章的技术要求部分。

你将使用来自 statsmodels 库的acorr_ljungbox。

如何操作…

你将使用存储在co2_df DataFrame 中的 CO2 数据集：

1. 从 statsmodels 库加载acorr_ljungbox：

from statsmodels.stats.diagnostic import acorr_ljungbox

1. 由于数据不是平稳的（请回顾检测时间序列平稳性食谱），这次你将进行对数变换（对数差分）：

co2_diff= np.log(co2_df).diff().dropna()

1. 运行 Ljung-Box 检验。首先设置lags=10：

acorr_ljungbox(co2_diff, lags=10, return_df=True)

这将打印出前 10 个滞后的结果。

图 9.31：自相关检验的前 10 个滞后

这表明，对于所有滞后期直到滞后 10，自检验统计量都显著（p 值 < 0.05），因此你可以拒绝原假设。拒绝原假设意味着你拒绝没有自相关的假设。

它是如何工作的…

acorr_ljungbox是一个累积自相关直到指定滞后的函数。因此，它有助于确定结构是否值得建模。

还有更多…

我们将对model_bx模型的残差使用 Ljung-Box 检验，该模型是在应用幂变换食谱中创建的：

acorr_ljungbox(model_bx.resid, return_df=True, lags=10)

这将打印出前 10 个滞后的结果：

图 9.32：自相关检验的前 10 个滞后，针对残差

从前面的例子中，p 值小于 0.05，因此你拒绝原假设，并且存在自相关。

另见

要了解更多关于acorr_ljungbox函数的信息，请访问官方文档：www.statsmodels.org/dev/generated/statsmodels.stats.diagnostic.acorr_ljungbox.html。

第十章：10 使用统计方法构建单变量时间序列模型

加入我们的书籍社区，加入 Discord

packt.link/zmkOY

在第九章，探索性数据分析与诊断中，你已经接触了帮助你理解时间序列过程的几个概念。这些方法包括时间序列数据分解、检测时间序列的平稳性、应用幂次变换和测试时间序列数据的自相关性。这些技术将在本章讨论的统计建模方法中派上用场。

在处理时间序列数据时，可以根据你所处理的时间序列是单变量还是多变量、季节性还是非季节性、平稳还是非平稳、线性还是非线性，采用不同的方法和模型。如果你列出需要考虑和检查的假设——例如平稳性和自相关性——就会明显看出，为什么时间序列数据被认为是复杂和具有挑战性的。因此，为了建模这样一个复杂的系统，你的目标是获得一个足够好的近似，捕捉到关键的兴趣因素。这些因素会根据行业领域和研究目标的不同而有所变化，例如预测、分析过程或检测异常。

一些流行的统计建模方法包括指数平滑、自回归积分滑动平均（ARIMA）、季节性 ARIMA（SARIMA）、向量自回归（VAR）以及这些模型的其他变种，如 ARIMAX、SARIMAX、VARX 和 VARMA。许多实践者，如经济学家和数据科学家，仍然使用这些统计“经典”模型。此外，这些模型可以在流行的软件包中找到，如 EViews、MATLAB、Orange、KNIME 和 Alteryx，以及 Python 和 R 的库中。

在本章中，你将学习如何在 Python 中构建这些统计模型。换句话说，我只会简要介绍理论和数学，因为重点在于实现。如果你对这些模型的数学和理论感兴趣，我会在适当的地方提供参考文献，供你深入研究。

在本章中，我们将涵盖以下内容：

• 绘制 ACF 和 PACF

• 使用指数平滑法预测单变量时间序列数据

• 使用非季节性 ARIMA 预测单变量时间序列数据

• 使用季节性 ARIMA 预测单变量时间序列数据

• 使用 Auto_Arima 预测单变量时间序列数据

在深入这些方法之前，请特别注意即将到来的技术要求部分，在这一部分你将进行前期准备。这将消除任何干扰和重复编码，以便你可以专注于方法的核心目标以及每个实现背后的概念。

技术要求

你可以从本书的 GitHub 仓库下载 Jupyter Notebooks 和必要的数据集：

• Jupyter Notebook: github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook./blob/main/code/Ch10/Chapter%2010.ipynb

• 数据集: github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook./tree/main/datasets/Ch10

在开始学习本章的配方之前，请运行以下代码以加载将在全章中引用的数据集和函数：

1. 首先导入本章所有配方中将共享的基本库：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import warnings
from statsmodels.tsa.api import (kpss, adfuller,
                              seasonal_decompose, STL)
from statsmodels.tools.eval_measures import rmspe, rmse
from sklearn.metrics import mean_absolute_percentage_error as mape
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from itertools import product
from pathlib import Path
warnings.filterwarnings('ignore')
plt.rcParams["figure.figsize"] = [12, 5]
plt.style.use('grayscale')
warnings.filterwarnings('ignore')

1. 本章中您将使用两个数据集：Life Expectancy from Birth 和 Monthly Milk Production。将这两个以 CSV 格式存储（life_expectancy_birth.csv 和 milk_production.csv）的数据集导入到 pandas DataFrame 中。每个数据集来自不同的时间序列过程，因此它们将包含不同的趋势或季节性。一旦导入数据集，您将得到两个名为 life 和 milk 的数据框：

life_file = Path('../../datasets/Ch10/life_expectancy_birth.csv')
milk_file = Path('../../datasets/Ch10/milk_production.csv')
life = pd.read_csv(life_file,
                   index_col='year',
                   parse_dates=True,
                  skipfooter=1)
milk = pd.read_csv(milk_file,
                   index_col='month',
                   parse_dates=True)

直观地检查数据，观察时间序列是否包含任何趋势或季节性。您可以随时返回本节中展示的图表作为参考：

fig, ax = plt.subplots(2, 1, figsize=(16, 12))
life.plot(title='Annual Life Expectancy',
                        legend=False, ax=ax[0])
milk.plot(title='Monthly Milk Production',
                        legend=False, ax=ax[1]);­

这将显示两个时间序列图：

图 10.1：年度预期寿命和每月牛奶生产的时间序列图

上图展示了life expectancy 数据框的时间序列图，显示了一个积极（向上）的趋势，且没有季节性。该预期寿命数据包含了 1960 年至 2018 年（59 年）的每年出生时的预期寿命记录。原始数据集包含了各国的记录，但在本章中，您将使用全球的记录。

monthly milk production 数据框的时间序列图显示了一个积极（向上）的趋势，并且呈现出周期性（每年夏季）。该牛奶生产数据从 1962 年 1 月到 1975 年 12 月每月记录（共 168 个月）。季节性幅度和随时间变化的波动似乎保持稳定，表明其为加法模型。季节性分解明确了加法模型的水平、趋势和季节性，也能反映这一点。

如需深入了解季节性分解，请参考第九章中的季节性分解时间序列数据配方，位于探索性数据分析与诊断部分。

1. 你需要将数据分割成 test 和 train 数据集。你将在训练数据集上训练模型（拟合），并使用测试数据集来评估模型并比较预测结果。基于用于训练的数据所创建的预测称为样本内预测，而对未见数据（如测试集）进行的预测称为样本外预测。在评估不同模型时，你将使用样本外或测试集。

创建一个通用函数 split_data，根据测试拆分因子分割数据。这样，你也可以尝试不同的拆分方式。我们将在本章中引用这个函数：

 def split_data(data, test_split):
    l = len(data)
    t_idx = round(l*(1-test_split))
    train, test = data[ : t_idx], data[t_idx : ]
    print(f'train: {len(train)} , test: {len(test)}')
    return train, test

1. 调用 split_data 函数，将两个数据框分割为 test 和 train 数据集（开始时使用 15% 测试集和 85% 训练集）。你可以随时尝试不同的拆分因子：

test_split = 0.15
milk_train, milk_test = split_data(milk, test_split)
life_train, life_test = split_data(life, test_split)
>>
train: 143 , test: 25
train: 50 , test: 9

1. 你将经常检查平稳性，因为这是你将构建的许多模型的一个基本假设。例如，在第九章，探索性数据分析与诊断，在检测时间序列平稳性的食谱中，我们讨论了检验平稳性的重要性以及使用**扩展的迪基-富勒（Augmented Dickey-Fuller）**检验。创建一个可以在本章中引用的函数，用于执行检验并解释结果：

def check_stationarity(df):
    results = adfuller(df)[1:3]
    s = 'Non-Stationary'
    if results[0] < 0.05:
        s = 'Stationary'
    print(f"'{s}\t p-value:{results[0]} \t lags:{results[1]}")
    return (s, results[0])

1. 在某些食谱中，你将运行模型的多个变体，寻找最优配置，这一做法通常被称为超参数调优。例如，你可以训练一个 ARIMA 模型并使用不同的参数值，从而生成多个 ARIMA 模型的变体（多个模型）。

get_top_models_df 函数将比较不同的模型——例如，多个 ARIMA 模型——以选择最佳模型及其相关的参数集。get_top_models_df 函数将接受一个字典，包含生成的模型、相关的参数和每个模型的评分。它返回一个数据框，详细列出表现最好的模型，便于比较。该函数允许你指定返回的最佳模型数量和用于选择模型的 criterion，例如均方根百分比误差（RMSPE）、均方根误差（RMSE）、平均绝对百分比误差（MAPE）、赤池信息量准则（AIC）、修正赤池信息量准则（AICc）或贝叶斯信息量准则（BIC）。这些指标在确定最适合你的数据分析需求的模型时至关重要。

例如，默认情况下你可能会选择根据 AIC 分数评估模型，但如果 RMSPE 或 RMSE 更适合你的具体情况，你可以轻松切换到这些指标。这种灵活性确保你可以根据数据的分析需求和复杂性量身定制模型选择过程。

def get_top_models_df(scores, criterion='AIC', top_n=5):
    sorted_scores = sorted(scores.items(),
                           key=lambda item: item[1][criterion])

    top_models = sorted_scores[:top_n]
    data = [v for k, v in top_models]
    df = pd.DataFrame(data)

    df['model_id'] = [k for k, v in top_models]
    df.set_index('model_id', inplace=True)
    return df

1. 创建plot_forecast函数，该函数接受你已经训练的模型对象、起始位置以及训练和测试数据集，生成一个将预测值与实际值进行比较的图表。当你深入本章的配方时，情况将变得更加清晰：

def plot_forecast(model, start, train, test):
    forecast = pd.DataFrame(model.forecast(test.shape[0]),
                            index=test.index)
    ax = train.loc[start:].plot(style='--')
    test.plot(ax=ax)
    forecast.plot(ax=ax, style = '-.')
    ax.legend(['orig_train', 'orig_test', 'forecast'])
    plt.show()

1. 最后，创建一个combinator工具函数，该函数接受一组参数值并返回这些选择的笛卡尔积。在进行超参数调优时，你将使用这个函数进行网格搜索。在网格搜索中，你会指定一组参数值的组合，分别训练多个模型，然后使用get_top_models_df函数评估最佳模型。例如，假设你的列表中包含三个不同参数的三个可能值。在这种情况下，combinator函数将返回一个包含 3x3 或九种可能组合的列表。当你深入本章的配方时，情况将变得更加清晰：

def combinator(items):
    combo = [i for i in product(*items)]
    return combo

我们可以像图 10.2中所示的那样表示整体流程，图中展示了你如何利用刚刚创建的函数。

图 10.2：本章利用已准备好的辅助函数创建的整体过程

现在，让我们深入了解这些配方。

在第一个配方中，你将接触到ACF和PACF图，它们用于评估模型拟合、检查平稳性，并确定本章中将使用的一些模型（如 ARIMA 模型）的阶数（参数）。

绘制 ACF 和 PACF

在构建任何统计预测模型之前，如自回归（AR）、移动平均（MA）、自回归移动平均（ARMA）、自回归积分移动平均（ARIMA）或季节性自回归积分移动平均（SARIMA），你需要确定最适合你数据的时间序列模型类型。此外，你还需要确定一些必需参数的值，这些参数称为阶数。更具体地说，这些包括自回归（AR）或移动平均（MA）成分的滞后阶数。这个过程将在本章的*“使用 ARIMA 预测单变量时间序列数据”*部分进一步探讨。例如，自回归移动平均（ARMA）模型表示为ARMA(p, q)，其中'p'表示自回归阶数或 AR(p)成分，'q'表示移动平均阶数或 MA(q)成分。因此，ARMA 模型结合了 AR(p)和 MA(q)模型。

这些模型的核心思想基于这样一个假设：可以通过过去的数值来估计当前特定变量的值！。例如，在一个自回归模型（AR）中，假设当前值！，在时间！可以通过其过去的值（）估算，直到p，其中p决定了我们需要回溯多少个时间步。如果！，这意味着我们必须使用前两个周期！来预测！。根据时间序列数据的粒度，p=2可能代表 2 小时、2 天、2 个月、2 季度或 2 年。

要构建 ARMA(p,q)模型，你需要提供p和q阶数（即滞后）。这些被视为超参数，因为它们是由你提供的，用于影响模型。

参数和超参数这两个术语有时被交替使用。然而，它们有不同的解释，你需要理解它们之间的区别。

参数与超参数

在训练 ARMA 或 ARIMA 模型时，结果将生成一组被称为系数的参数——例如，AR 滞后 1 的系数值或 sigma——这些系数是算法在模型训练过程中估算出来的，并用于进行预测。它们被称为模型的参数。

另一方面，（p，d，q）参数是 ARIMA(p，q，d)中的 AR、差分和 MA 的阶数。这些被称为超参数。它们在训练过程中提供，并影响模型生成的参数（例如系数）。这些超参数可以通过网格搜索等方法进行调优，以找到生成最佳模型的最佳参数组合。

现在，你可能会问，如何找到 AR 和 MA 模型的显著滞后值？

这就是自相关函数（ACF）和偏自相关函数（PACF）以及它们的图形发挥作用的地方。可以绘制 ACF 和 PACF 图，以帮助你识别时间序列过程是 AR、MA 还是 ARMA 过程（如果两者都存在），并识别显著的滞后值（对于p和q）。ACF 和 PACF 图都被称为相关图，因为这些图表示相关性统计。

ARMA 和 ARIMA 的区别，在于平稳性假设。ARIMA 中的d参数用于差分阶数。ARMA 模型假设过程是平稳的，而 ARIMA 模型则不假设，因为它处理差分问题。ARIMA 模型是一个更为广泛的模型，因为通过将差分因子d=0，它可以满足 ARMA 模型的需求。因此，ARIMA(1,0,1)就是ARMA(1,1)。

AR阶数与MA阶数

你将使用 PACF 图来估计 AR 阶数，并使用 ACF 图来估计 MA 阶数。ACF 和 PACF 图的纵轴（y 轴）显示从 -1 到 1 的值，横轴（x 轴）表示滞后的大小。显著的滞后是任何超出阴影置信区间的滞后，正如你在图中看到的。

statsmodels 库提供了两个函数：acf_plot 和 pacf_plot。零滞后的相关性（对于 ACF 和 PACF）始终为 1（因为它表示第一个观测值与自身的自相关）。因此，这两个函数提供了 zero 参数，该参数接受一个布尔值。因此，为了在可视化中排除零滞后，可以传递 zero=False。

在 第九章，探索性数据分析与诊断 中，测试时间序列数据中的自相关 配方中，你使用了 Ljung-Box 测试来评估残差的自相关。在本例中，你还将学习如何使用 ACF 图来直观地检查 残差自相关。

如何操作…

在本例中，你将探索 statsmodels 库中的 acf_plot 和 pacf_plot。我们开始吧：

1. 你将在本例中使用寿命预期数据。如 图 10.1 所示，由于存在长期趋势，数据不是平稳的。在这种情况下，你需要对时间序列进行差分（去趋势），使其平稳，然后才能应用 ACF 和 PACF 图。

从差分开始，然后创建不包含零滞后的图：

life_diff = life.diff().dropna()
fig, ax = plt.subplots(2,1, figsize=(12,8))
plot_acf(life_diff, zero=False, ax=ax[0])
plot_pacf(life_diff, zero=False, ax=ax[1]);

这将生成以下两个图：

图 10.3：差分后寿命预期数据的 ACF 和 PACF 图

如果你想查看更多滞后的计算 PACF 和 ACF，可以更新 lags 参数，如下所示：

fig, ax = plt.subplots(2,1, figsize=(12,8))
plot_acf(life_diff, lags=25, zero=False, ax=ax[0])
plot_pacf(life_diff, lags=25,  zero=False, ax=ax[1]);

这将生成以下两个图：

图 10.4：前 25 个滞后的 ACF 和 PACF 图

ACF 图在滞后（阶数）1 处显示了一个显著的峰值。当滞后（垂直线）超出阴影区域时，表示显著性。阴影区域代表置信区间，默认设置为 95%。在 ACF 图中，只有第一个滞后显著，低于下置信区间，并且在此之后 迅速消失。其余的滞后均不显著。这表明是一个一阶移动平均 MA(1)。

PACF 图呈现 渐变 衰减并伴随震荡。通常，如果 PACF 显示渐变衰减，表示使用了移动平均模型。例如，如果你使用 ARMA 或 ARIMA 模型，一旦数据差分使其平稳后，它将表现为 ARMA(0, 1) 或 ARIMA(0, 1, 1)，表示一阶差分，d=1。在 ARMA 和 ARIMA 模型中，AR 阶数为 p=0，MA 阶数为 q=1。

1. 现在，让我们看看如何将 PACF 和 ACF 用于包含强烈趋势和季节性的更复杂的数据集。在图 10.1中，月度牛奶生产图显示了年度季节性效应和一个正向上升的趋势，表示这是一个非平稳时间序列。它更适合使用 SARIMA 模型。在 SARIMA 模型中，你有两个部分：非季节性部分和季节性部分。例如，除了之前看到的表示非季节性部分的 AR 和 MA 过程（分别由小写p和q表示），你还会有季节性部分的 AR 和 MA 阶数，分别由大写P和Q表示。这个模型可以表示为SARIMA(p,d,q)(P,D,Q,S)。你将在使用季节性 ARIMA 预测单变量时间序列数据的食谱中了解更多关于 SARIMA 模型的内容。

为了使此类时间序列平稳，你需要从季节性差分开始，以去除季节性效应。由于观测是按月进行的，因此季节性效应是按年观察的（每 12 个月或每个周期）：

milk_diff_12 = milk.diff(12).dropna()

1. 使用你在本章早些时候创建的check_stationarity函数，执行增广的迪基-富勒（Augmented Dickey-Fuller）检验以检查平稳性：

check_stationarity(milk_diff_12)
>> 'Non-Stationary     p-value:0.16079880527711382      lags:12

1. 差分后的时间序列仍然不是平稳的，因此你仍然需要进行第二次差分。这一次，你必须执行一阶差分（去趋势）。当时间序列数据包含季节性和趋势时，你可能需要进行两次差分才能使其平稳。将结果存储在milk_diff_12_1变量中，并再次运行check_stationarity：

milk_diff_12_1 = milk.diff(12).diff(1).dropna()
check_stationarity(milk_diff_12_1)
>> 'Stationary     p-value:1.865423431878876e-05     
lags:11

太棒了——现在，你有了一个平稳过程。

1. 绘制milk_diff_12_1中平稳时间序列的 ACF 和 PACF 图：

fig, ax = plt.subplots(1,2)
plot_acf(milk_diff_12_1, zero=False, ax=ax[0], lags=36)
plot_pacf(milk_diff_12_1, zero=False, ax=ax[1], lags=36);

这应该会生成以下 ACF 和 PACF 图：

图 10.5: 经差分后的月度牛奶生产的 PACF 和 ACF 图

对于季节性阶数P和Q，你应该诊断在滞后s、2s、3s等位置的峰值或行为，其中s是一个季节中的周期数。例如，在牛奶生产数据中，s=12（因为一个季节有 12 个月）。然后，我们观察 12（s）、24（2s）、36（3s）等位置的显著性。

从自相关函数（ACF）图开始，滞后 1 处有一个显著的峰值，表示 MA 过程的非季节性阶数为q=1。滞后 12 处的峰值表示 MA 过程的季节性阶数为Q=1。注意，在滞后 1 之后有一个截断，然后在滞后 12 处有一个峰值，之后又是截断（没有其他显著的滞后）。这些现象表明该模型是一个移动平均模型：非季节性部分是 MA(1)，季节性部分是 MA(1)。偏自相关函数（PACF）图也证实了这一点；在滞后 12、24 和 36 处的指数衰减表明这是一个 MA 模型。因此，SARIMA 模型应该是ARIMA (0, 1,1)(0, 1, 1, 12)。

尽管使用 ACF 和 PACF 图对于识别 ARIMA 阶数 p 和 q 很有用，但不应单独使用。在本章中，您将探索不同的技术，帮助您通过 AIC 和 BIC 等模型选择技术来确定阶数。

工作原理…

ACF 和 PACF 图可以帮助您理解过去观测值之间线性关系的强度以及在不同滞后期的显著性。

ACF 和 PACF 图显示出显著的自相关或部分自相关，超出了 置信区间。阴影部分表示置信区间，它由 pacf_plot 和 acf_plot 函数中的 alpha 参数控制。statsmodels 中 alpha 的默认值是 0.05（95% 置信区间）。显著性可以是任意方向；如果自相关强烈为正，它会接近 1（上方），如果强烈为负，则会接近 -1（下方）。

下表展示了一个示例指南，用于从 PACF 和 ACF 图中识别平稳的 AR 和 MA 阶数：

表 10.1：使用 ACF 和 PACF 图识别 AR、MA 和 ARMA 模型

还有更多…

在本食谱中，您使用了 ACF 和 PACF 图来估计应该为季节性和非季节性 ARIMA 模型使用哪些阶数（滞后）。

让我们看看 ACF 图如何用于诊断模型的 残差。检查模型的残差是评估模型的重要组成部分。这里的假设很简单：如果模型正确地捕捉到了所有必要的信息，那么残差中不应该包含任何在任意滞后期有相关性的点（即没有自相关）。因此，您会期望残差的 ACF 图显示出接近零的自相关。

让我们构建之前在本食谱中识别的季节性 ARIMA 模型 SARIMA(0,1,1)(0,1,1,12)，然后使用 ACF 来诊断残差。如果模型捕捉到了时间序列中嵌入的所有信息，您会期望残差 没有自相关：

from statsmodels.tsa.statespace.sarimax import SARIMAX
model = SARIMAX(milk, order=(0,1,1),
                seasonal_order=(0,1,1, 12)).fit(disp=False)
plot_acf(model.resid, zero=False, lags=36);

这将生成以下自相关图：

图 10.6：SARIMA 残差的自相关图

总体来说，SARIMA(0,1,1)(0,1,1,12) 很好地捕捉到了必要的信息，但仍然有提升的空间。存在一个显著的滞后期（在滞后=12 时，超出了置信阈值），表明残差中存在某些自相关。

您可以进一步调整模型并尝试不同的季节性和非季节性阶数。在本章及后续食谱中，您将探索一种网格搜索方法来选择最佳的超参数，以找到最优模型。

如果您想进一步诊断模型的残差，您可以使用 plot_diagnostics 方法进行：

model.plot_diagnostics(figsize=(12,7), lags=36);

这将生成以下图形：

图 10.7：使用 plot_diagnostics 方法进行的残差分析

请注意，生成的诊断图是基于标准化残差的，这是一种常见的技术，因为它使得在不同模型之间比较残差更加容易，因为它们是归一化的，并以标准差的形式表示。

你可以通过访问model.standardized_forecasts_error来复制相同的图表，如下所示：

plot_acf(model.standardized_forecasts_error.ravel(), lags=36,
         title='Standardized Residuals ACF Plot');

和

pd.DataFrame(model.standardized_forecasts_error.ravel(),
             index=milk.index).plot(title='Standardized Residuals Plot',
                                    legend=False);

生成的两个图应类似于以下图形，展示标准化残差的自相关和时间序列模式：

图 10.8：SARIMA 模型标准化残差的 ACF 图

图 10.9：SARIMA 模型标准化残差图

图 10.6与图 10.8之间的差异是由于尺度归一化。标准化可以减少异常值的影响，因为残差被缩小。这就是为什么在图 10.8 中，尽管 Lag 12 的自相关在两张图中都可见，但它位于置信区间的边界，而在图 10.6中则明显不同。同时，标准化可能会放大一些较小的自相关，这些自相关在原始的 ACF 图中可能最初是不可见的。在整个食谱中，你将依赖于plot_diagnostics方法进行残差诊断。

另见

• 若要了解更多关于 ACF 图的内容，请访问官方文档 www.statsmodels.org/dev/generated/statsmodels.graphics.tsaplots.plot_acf.html.

• 若要了解更多关于 PACF 图的内容，请访问官方文档 www.statsmodels.org/dev/generated/statsmodels.graphics.tsaplots.plot_pacf.html.

这样，你就知道在构建 ARIMA 模型及其变体（例如，ARMA 或 SARIMA 模型）时，如何使用 ACF 和 PACF 图。在接下来的食谱中，你将学习本章的第一个时间序列预测技术。

使用指数平滑法预测单变量时间序列数据

在本食谱中，你将探索使用statsmodels库的指数平滑技术，它提供了与 R 语言forecast包中的流行实现（如ets()和HoltWinters()）类似的功能。在 statsmodels 中，指数平滑有三种不同的实现（类），具体取决于你所处理数据的性质：

• SimpleExpSmoothing：当时间序列过程缺乏季节性和趋势时，使用简单指数平滑法。这也被称为单一指数平滑法。

• Holt：Holt 指数平滑是简单指数平滑的增强版本，用于处理只包含趋势（但没有季节性）的时间序列过程。它被称为双指数平滑。

• ExponentialSmoothing：Holt-Winters 指数平滑是 Holt 指数平滑的增强版本，用于处理同时具有季节性和趋势的时间序列过程。它被称为三重指数平滑。

你可以像下面这样导入这些类：

from statsmodels.tsa.api import (ExponentialSmoothing,
                                SimpleExpSmoothing,
                                Holt)

statsmodels 的实现遵循*《预测：原理与实践》（作者：Hyndman, Rob J., and George Athanasopoulos）*中的定义，你可以在这里参考：otexts.com/fpp3/expsmooth.html。

如何实现…

在这个示例中，你将对本章介绍的两个数据集进行指数平滑处理（技术要求）。由于Holt类和SimpleExpSmoothing类是ExponentialSmoothing类的简化版本，因此你可以使用后者来进行简单操作。与使用这三者不同，你可以使用ExponentialSmoothing类运行这三种不同类型，因为ExponentialSmoothing是更通用的实现。这种方法允许你使用单一更多功能的实现来管理不同类型的时间序列，无论它们是否表现出趋势、季节性或两者都有。让我们开始吧：

1. 导入ExponentialSmoothing类：

from statsmodels.tsa.api import ExponentialSmoothing

1. 你将从寿命数据集开始，并使用ExponentialSmoothing类。

ExponentialSmoothing 需要多个参数（称为超参数），可以分为两种类型：在构建模型时指定的参数和在拟合模型时指定的参数。

1. 模型构建：

   • trend：选择从‘multiplicative’（别名‘mul’）、‘additive’（别名‘add’）或 None 中进行选择。

   • seasonal：选择从‘multiplicative’（别名‘mul’）、‘additive’（别名‘add’）或 None 中进行选择。

   • seasonal_periods：代表季节周期的整数；例如，对于月度数据使用 12，季度数据使用 4。

   • damped_trend：布尔值（True 或 False），指定趋势是否应该被阻尼。

   • use_boxcox：布尔值（True 或 False），确定是否应用 Box-Cox 变换。

2. 模型拟合：

   • smoothing_level：浮点数，指定作为alpha的平滑水平的平滑因子。

     

     ), 其有效值在 0 和 1 之间 (

     

     ).

   • smoothing_trend：浮点数，指定作为beta的趋势的平滑因子 (

     

     ), 其有效值在 0 和 1 之间 (

     

     ).

   • smoothing_seasonal：浮点数，指定作为gamma的季节性趋势的平滑因子。

     

     ), 其有效值在 0 和 1 之间 (

     

     ).

在后续的*如何运作...*部分，你将探索**霍尔特-温特斯（Holt-Winters）**的水平、趋势和季节性公式，以及这些参数是如何应用的。

1. 创建一个包含不同超参数值组合的列表。这样，在下一步中，你可以在每次运行中评估不同的超参数值组合。从本质上讲，你将训练不同的模型，并在每次迭代中记录其得分。一旦每个组合都被评估，你将使用get_top_models_df函数（来自技术要求部分）来确定表现最好的模型及其最佳超参数值，通过这个详尽的网格搜索过程。这个过程可能需要耗费一些时间，但幸运的是，有一种混合技术可以缩短搜索时间。

你可以使用ExponentialSmoothing类来找到alpha、beta和gamma的最佳值（）。这种方法无需在网格中指定这些值（尽管如果你更愿意控制过程，仍然可以指定）。这种简化意味着你只需提供剩余超参数的值，如trend和seasonal。你可以通过使用seasonal_decompose()函数绘制它们的分解，初步判断这些组件是乘法型还是加法型。如果仍不确定，详尽的网格搜索仍然是一个可行的替代方案。

对于life数据框，只有trend，因此你只需要探索两个参数的不同值；即trend和damped：

trend = ['add', 'mul']
damped = [True, False]
life_ex_comb = combinator([trend, damped])
life_ex_comb
[('add', True), ('add', False), ('mul', True), ('mul', False)]

在这里，我们有两个参数，每个参数有两个不同的值，这为我们提供了 2x2 或四个总的组合来评估。

1. 遍历组合列表，并在每次迭代中训练（拟合）一个不同的模型。将评估指标捕捉到字典中，以便稍后比较结果。你将捕捉到的示例得分包括 RMSE、RMSPE、MAPE、AIC 和 BIC 等。请记住，大多数自动化工具和软件会在后台使用 AIC 和 BIC 分数来确定最佳模型：

train = life_train.values.ravel()
y = life_test.values.ravel()
score = {}
for i, (t, dp) in enumerate(life_ex_comb):
    exp = ExponentialSmoothing(train,
                     trend=t,
                     damped_trend=dp,
                     seasonal=None)
    model = exp.fit(use_brute=True, optimized=True)   
    y_hat = model.forecast(len(y))
    score[i] = {'trend':t,
                'damped':dp,
                'AIC':model.aic,
                'BIC':model.bic,
                'AICc':model.aicc,
                'RMSPE': rmspe(y, y_hat),
                'RMSE' : rmse(y, y_hat),
                'MAPE' : mape(y, y_hat),
                'model': model}

在前面的函数中，你使用life_train来训练不同的模型，使用life_test来评估错误度量，如 RMSPE、RMSE 和 MAPE。

要使用get_top_models_df函数获取前几个模型，只需传递得分字典。目前，保持默认标准设置为c=AIC以保持一致性：

model_eval = get_top_models_df(score, 'AIC', top_n=5)

1. get_top_models_df函数将返回一个 DataFrame，显示排名前五的模型（默认为 5），根据所选标准进行排名，例如在此案例中是 AIC 分数。DataFrame 不仅包含所有附加得分，还将模型实例本身存储在名为'model'的列中。

要查看排名和各种得分，你可以执行以下代码：

model_eval.iloc[:, 0:-1]

上述代码排除了最后一列，该列包含模型实例，因此显示的 DataFrame 包括 AIC、BIC、RMSE 等每个评估指标的列。

图 10.10：基于 AIC 分数排名的寿命预期数据的指数平滑模型

通常，对于基于信息准则（如 AIC、BIC 和 AICc）进行的模型选择，较低的值更好，表示模型拟合和复杂度之间的更优平衡。在我们的案例中，我们选择使用 AIC。如果您检查图 10.10 中的 DataFrame，会发现一些可以观察到的现象：

• 如果优先考虑信息准则（AIC、BIC、AICc），模型 1（趋势：加性，阻尼：False）会被视为最佳模型，因为它在所有三个信息准则中得分最低。该模型可能在模型复杂度和拟合度之间提供了最佳的折衷。

• 如果优先考虑误差指标（RMSPE、RMSE、MAPE），这些指标用于衡量预测精度，模型 0（趋势：加性，阻尼：True）会被认为是更优的，因为它的预测误差较小。

选择“获胜”模型将取决于您的具体目标和模型将被使用的上下文。如果您需要两者之间的平衡，您可能需要考虑其他因素或进一步验证，以决定选择模型 1 还是模型 0。

我们将继续使用 AIC 作为我们的选择标准。

1. 存储在 DataFrame 中的模型是 HoltWintersResultsWrapper 类的实例。您可以直接从 DataFrame 访问顶部模型，这样可以利用与该模型相关的其他方法和属性，如 summary、predict 和 forecast。要提取并与第一行的获胜模型进行交互，请使用以下代码：

top_model = model_eval.iloc[0,-1]

您可以像下面这样访问 summary() 方法：

top_model.summary()

上述代码将生成一个总结输出，提供一个表格布局，详细展示模型——例如，使用的参数值和计算出的系数：

图 10.11：寿命预期数据的指数平滑总结

总结将显示关键的信息，例如通过拟合过程自动推导出的alpha（平滑水平）和beta（平滑趋势）的最佳值。

1. 您可以使用 forecast 方法预测未来的值，然后将结果与测试集（模型未见过的数据）进行比较。我们在本章的技术要求部分介绍的 plot_forecast() 函数将用于生成并绘制预测结果，同时显示测试数据。要执行此可视化，将存储在 top_model 中的模型对象与 training 和 test 数据集一起传递给 plot_forecast()：

plot_forecast(life_best_model, '2000', life_train, life_test)

plot_forecast函数中的start参数将数据从该点开始切片，以便更容易比较结果。可以将其视为聚焦于时间线的特定片段。例如，不是显示 1960 到 2018 年（59 个月）的数据，而是只请求从 2000 年开始的这一段数据。

这将生成一个图，其中 x 轴从 2000 年开始。应有三条线：一条线表示训练数据，另一条线表示测试数据，还有一条线表示预测值（预测值）：

图 10.12：将指数平滑预测与生命预期数据集的实际数据进行对比

简单指数平滑的预测结果是延伸自训练数据的上升趋势的直线。

1. 重复之前的过程，但使用milk数据框。请记住，这里最重要的区别是添加了季节性参数。这意味着你将添加两个额外的超参数来进行评估——即seasonal和seasonal_periods。

为不同选项构建笛卡尔积。对于seasonal_periods，你可以探索三个周期——4、6 和 12 个月。这应该会给你提供 24 个需要评估的模型：

trend , damped= ['add', 'mul'], [True, False]
seasonal, periods = ['add' , 'mul'], [4, 6, 12]
milk_exp_comb = combinator([trend, damped, seasonal, periods])

循环遍历组合列表，训练多个模型并捕获它们的得分：

train = milk_train.values.ravel()
y = milk_test.values.ravel()
milk_model_scores = {}
for i, (t, dp, s, sp) in enumerate(milk_exp_comb):
    exp = ExponentialSmoothing(train,
                        trend=t,
                        damped_trend=dp,
                        seasonal=s,
                        seasonal_periods=sp)
    model = exp.fit(use_brute=True, optimized=True)   
    y_hat = model.forecast(len(y))
    score[i] = {'trend':t,
                'damped':dp,
                'AIC':model.aic,
                'BIC':model.bic,
                'AICc': model.aicc,
                'RMSPE': rmspe(y, y_hat),
                'RMSE' : rmse(y, y_hat),
                'MAPE' : mape(y, y_hat),
                'model': model}

1. 训练完成后，运行get_top_models_df函数，根据 AIC 得分识别出最佳模型：

model_eval = get_top_models_df(score, 'AIC', top_n=5)
model_eval.iloc[:, 0:-1]

这将显示以下数据框：

图 10.13：基于 AIC 得分排名的牛奶生产数据前 5 个指数平滑模型

要从结果中确定获胜模型，通常会查看各种度量标准，如 AIC、BIC、AICc（信息准则）以及误差度量，如 RMSPE、RMSE 和 MAPE。AIC、BIC 和 AICc 的较低值表示模型在拟合优度和复杂性之间有更好的平衡。RMSPE、RMSE 和 MAPE 的较低值表示更好的预测精度。

如果你检查图 10.13 中的数据框，你会发现有一些观察结果：

如果优先考虑信息准则（AIC、BIC、AICc），模型 8（趋势：加性，阻尼：False）似乎是最佳模型，因为它在所有信息准则中具有最低的值。这表明它在很好地拟合数据和保持简单性之间提供了一个有利的平衡。

如果优先考虑误差度量（RMSPE、RMSE、MAPE），模型 2（趋势：加性，阻尼：True）在预测精度方面表现更好。该模型具有最低的误差率，表明它在列出的模型中最准确地预测了未来值。

选择“获胜”模型将取决于您的具体目标以及模型将用于的上下文。如果您需要在两种方法之间找到平衡，您可能需要考虑其他因素或进一步验证，以决定选择模型 8 还是模型 2。

我们将继续使用 AIC 作为选择标准。

您可以使用以下命令显示最佳模型的汇总信息：

top_model = model_eval.iloc[0,-1]
top_model.summary()

这应该生成一个汇总最佳模型的表格布局——例如，构建模型时使用的参数值和计算出的系数：

图 10.14：月度牛奶生产数据的指数平滑总结

请注意，趋势、季节性和季节周期的超参数值的最佳组合。最佳的 季节周期 为 12 个月或滞后期。汇总结果表将显示所有这些滞后的系数，并且这将是一个很长的列表。前面的截图仅显示了顶部部分。

此外，汇总还将显示关键的优化信息，如 alpha（平滑水平）、beta（平滑趋势）和 gamma（平滑季节性）的最佳值。

请记住，最佳模型是根据 AIC 分数选定的。因此，您应探索已捕捉到的不同指标，例如使用 get_top_models_df(score, 'MAPE', top_n=5)。

1. 将您使用最佳模型的预测与测试数据进行比较：

plot_forecast(top_model, '1969', milk_train, milk_test);

这应该会生成一个从 1969 年开始的图表，展示三条线，分别表示训练数据、测试数据和预测（预测值）：

图 10.15：绘制指数平滑预测与实际月度牛奶生产数据的对比

总体而言，模型有效地捕捉了趋势和季节性，且与测试集中的实际值高度吻合。

它是如何工作的……

平滑时间序列数据有多种方法，包括简单移动平均法、简单指数平滑法、霍尔特指数平滑法和霍尔特-温特指数平滑法。

移动平均模型将过去的值视为相同，而指数平滑模型则更加注重（加权）最近的观察数据。在指数平滑中，较早观察数据的影响逐渐减少（加权衰减），因此得名“指数”。这一方法基于这样一个逻辑假设：较新的事件通常比较旧的事件更为重要。例如，在日常时间序列中，昨天或前天的事件通常比两个月前的事件更为相关。

简单指数平滑（单一）的公式，适用于没有趋势或季节性的时间序列过程，如下所示：

这里，ExponentialSmoothing 类的目标是找到平滑参数 alpha 的最佳值（）。在这个公式中， 代表当前时刻的期望（平滑）水平， 和 分别是当前时刻和之前时刻的平滑水平值， 是当前时刻的观测值（）。alpha（）参数至关重要，它作为水平平滑参数，决定了模型是否更信任过去（）还是当前（）。因此，当 趋近于零时，第一项（）趋近于零，更多的权重会被放在过去；而当 趋近于一时， 项趋近于零，更多的权重则会放在当前。选择 的因素有很多，包括系统中的随机性程度。系数 的输出值决定了模型如何加权当前和过去的观测值来预测未来事件（）。

这个解释与类似公式中所呈现的主题一致；虽然我们不会深入探讨每个细节，但整体概念保持一致。

Holt 的指数平滑（双重）公式包含了趋势（）及其平滑参数 beta（）。因此，一旦加入趋势，模型将输出两个系数的值——即 alpha 和 beta（）： Holt-Winters 指数平滑（三重）公式同时包含趋势（）和季节性（）。以下公式展示了 乘法 季节性作为示例： 当使用 ExponentialSmoothing 寻找最佳参数值时，它是通过最小化误差率（误差平方和 或 SSE）来实现的。因此，在每次循环中，你传入新参数值（例如，damped 为 True 或 False），模型通过最小化 SSE 来求解最优的 系数值。这个过程可以写成如下公式：

在一些教科书中，你会看到不同的字母用于表示水平、趋势和季节性，但公式的整体结构是相同的。

通常，指数平滑是一种快速且有效的技术，用于平滑时间序列以改善分析，处理异常值、数据插补和预测（预测）。

还有更多…

一个名为 Darts 的令人兴奋的库提供了一个 ExponentialSmoothing 类，它是基于 statsmodels 的 ExponentialSmoothing 类的封装。

要使用 pip 安装 Darts，请运行以下命令：

pip install darts

要使用 conda 安装，请运行以下命令：

conda install -c conda-forge -c pytorch u8darts-all

加载 ExponentialSmoothing 和 TimeSeries 类：

from darts.models import ExponentialSmoothing
from darts import TimeSeries

Darts 期望数据是 TimeSeries 类的一个实例，因此在使用它来训练模型之前，你需要先将 pandas DataFrame 转换为 TimeSeries。TimeSeries 类提供了 from_dataframe 方法，你将在其中使用：

model = ExponentialSmoothing(seasonal_periods=12)
ts = TimeSeries.from_dataframe(milk.reset_index(),
                                    time_col='month', value_cols='production', freq='MS')

在创建 TimeSeries 对象时，你必须指定哪一列是日期，哪一列包含观测值（数据）。你可以使用 .fit() 方法训练模型。一旦训练完成，你可以使用 .predict() 方法进行预测。要绘制结果，可以使用 .plot() 方法：

train, test = split_data(ts, 0.15)
model.fit(train)
forecast = model.predict(len(test), num_samples=100)
train.plot()
forecast.plot(label='forecast', low_quantile=0.05, high_quantile=0.95)

图 10.16：使用 Darts 对每月牛奶生产数据进行 ExponentialSmoothing 预测

darts 库自动化了评估过程，以找到最佳配置（超参数）。Darts 的 ExponentialSmoothing 类是 statsmodels 的 ExponentialSmoothing 类的封装，这意味着你可以访问熟悉的方法和属性，例如 summary() 方法：

model.model.summary()

这应该会生成熟悉的 statsmodels 表格总结和优化后的参数值。作为挑战，请比较 Dart 的总结与 图 10.14 中的结果。尽管你会发现你达到了类似的结果，但使用 Darts 时付出的努力更少。它自动选择了 图 10.14 中确定的最佳超参数。

Darts 库还包含另一个有用的类，名为 StatsForecastAutoETS，其功能来源于 StatsForecast 库中的 AutoETS 实现。与传统的 ExponentialSmoothing 类相比，AutoETS 通常因其更快的性能而受到赞扬。

要探索 StatsForecastAutoETS 的功能，可以参考以下代码片段：

from darts.models import StatsForecastAutoETS
modelets = StatsForecastAutoETS(season_length=12)
modelets.fit(train)
etsforecast = modelets.predict(len(test))
train.plot()
etsforecast.plot(label='AutoETS');

图 10.17：使用 Darts 对每月牛奶生产数据进行 AutoETS 预测

你可以使用以下代码比较两个预测方法，ExponentialSmoothing 和 StatsForecastAutoETS：

forecast.plot(label='ExponentialSmoothing')
etsforecast.plot(label='StatsForecastAutoETS');

图 10.18：比较 AutoETS 和 ExponentialSmoothing

上面的线表示 ExponentialSmoothing 预测结果，而下面的线表示 StatsForecastAutoETS 预测结果。

指数平滑与 ETS

ETS 和指数平滑密切相关，因为它们都使用过去数据点的平滑方法来预测时间序列数据。然而，它们在方法上有所不同。指数平滑通过最小化平方误差和来估计参数，而 ETS 则通过最大化似然估计。此外，指数平滑提供点预测（预测值），ETS 也提供相同的点预测，但附带预测区间。

参见

要了解更多有关 ExponentialSmoothing 类的信息，您可以访问 statsmodels 的官方文档：www.statsmodels.org/dev/generated/statsmodels.tsa.holtwinters.ExponentialSmoothing.html。

您是否注意到在指数平滑中无需进行平稳性检验？指数平滑仅适用于非平稳时间序列（例如具有趋势或季节性的时间序列）。

在下一部分，在构建 ARIMA 模型时，您将进行平稳性检验，以确定差分阶数，并利用本章前面讨论过的 ACF 和 PACF 图。

使用 ARIMA 进行单变量时间序列数据预测

在本配方中，您将探索 ARIMA，使用 statsmodels 包进行实现。ARIMA 代表自回归集成滑动平均（Autoregressive Integrated Moving Average），它结合了三种主要成分：自回归或 AR(p) 模型、滑动平均或 MA(q) 模型和一个集成过程或 I(d)，它对数据应用差分。

ARIMA 模型通过 p、d 和 q 参数来表征，非季节性时间序列的 ARIMA 模型用符号 ARIMA(p, d, q) 来描述。p 和 q 参数表示阶数或滞后；例如，AR 的阶数为 p，MA 的阶数为 q。它们被称为滞后，因为它们表示我们需要考虑的“过去”时期的数量。您可能还会遇到另一种关于 p 和 q 的称呼，即多项式的阶数。

ARIMA 模型通过差分（一种时间序列转换技术）来处理非平稳时间序列数据，从而使非平稳时间序列变为平稳。差分的阶数或集成阶数 d 是构建模型时需要选择的参数之一。有关平稳性的复习，请参阅 第九章，探索性数据分析与诊断中的检测时间序列平稳性配方。

虽然 ARIMA 模型通过利用集成因子 'd' 设计用于处理趋势，但它们传统上假设数据集中没有季节性。然而，如果季节性是一个因素，那么季节性 ARIMA（SARIMA）模型就是合适的替代方案，因为它扩展了 ARIMA，包含季节性差分。

准备工作

首先，从 statsmodels 库加载本配方所需的类和函数：

from statsmodels.tsa.arima.model import ARIMA
from statsmodels.stats.diagnostic import acorr_ljungbox

如何操作……

不同的时间序列模型适用于各种类型的数据。因此，选择一个与数据集的特征以及你所解决的具体问题相符合的模型至关重要。在本例中，你将使用life数据框，它呈现出趋势但没有季节性。

你将结合视觉检查（使用 ACF 和 PACF 图）和统计检验，以便做出关于 AR 和 MA 模型成分（即p和q阶数）的明智决策。这些方法在第九章，探索性数据分析与诊断中已有介绍，包括自相关性检验、时间序列数据分解和时间序列平稳性检测的相关内容。让我们开始吧：

1. 首先通过分解数据集，将其分为三个主要成分：趋势、季节性和残差（通常被认为是噪声）。你可以使用seasonal_decompose函数来实现这一点。

decomposed = seasonal_decompose(life)
decomposed.plot();

你可以看到以下图表：

图 10.19：生命预期数据的分解

观察到分解结果显示数据集中有一个正向（上升）趋势。这表明数据随时间的一致增长。然而，数据中没有明显的季节性效应，这与我们对生命数据集的预期一致。

1. 你需要首先对数据进行去趋势处理。进行一次差分，然后使用本章前面创建的check_stationarity函数测试数据的平稳性：

check_stationarity(life)
>>
Non-Stationary   p-value:0.6420882853800064      lags:2
life_df1 = life.diff().dropna()
check_stationarity(life_df1)
>>
Stationary     p-value:1.5562189676003248e-14      lags:1

现在，数据是平稳的。p 值显著，可以拒绝原假设。请注意，diff()的默认periods值为1。通常，diff(periods=n)表示当前时期t的观测值与其滞后版本t-n之间的差异。对于diff(1)或diff()，滞后版本是t-1（例如，前一个月的观测值）。

你可以使用plot方法绘制差分后的时间序列数据：

life_df1.plot();

这将生成以下图表：

图 10.20：生命预期数据的一阶差分（去趋势）

接下来，你需要确定 ARIMA（p, d, q）模型的p和q阶数。

1. ACF 和 PACF 图将帮助你估计 AR 和 MA 模型的合适p和q值。对平稳化后的life_df1数据使用plot_acf和plot_pacf：

fig, ax = plt.subplots(1,2)
plot_acf(life_df1, ax=ax[0])
plot_pacf(life_df1, ax=ax[1]);

它会生成以下图表：

图 10.21：差分后的生命预期数据的 ACF 和 PACF 图

在前面的例子中，零滞后值被包含在图表中，以帮助你将其与过去的滞后期进行可视化比较。滞后期为 0 时，ACF 和 PACF 的值总是为 1；它们有时会被从图表中省略，因为它们并不提供有意义的信息。因此，更重要的是关注滞后期 1 及之后的滞后期，以确定它们的显著性。

ACF 图有助于识别 MA(q)组件的重要滞后。ACF 图在滞后 1 后显示截断，表示 MA(1)模型。相反，PACF 图有助于确定 AR(p)组件的重要滞后。你可以观察到滞后 1 后逐渐衰减并有振荡，表明在滞后 1 处是 MA 模型或 MA(1)。这表示没有 AR 过程，因此p阶数为零或 AR(0)。有关更多细节，请参考表 10.1。

MA(1)过程也叫做一阶移动平均过程，意味着当前值（在时间t时刻）受紧接其前一个值（在时间t-1时刻）的影响。

现在，你可以使用 p=0，q=1，d=1 的配置构建 ARIMA(p, d, q)模型，从而得到ARIMA(0,1,1)。通常，p 和 q 的最佳滞后值（阶数）并不一目了然，因此你需要通过评估多个 ARIMA 模型并使用不同的 p、d、q 参数来进行实验。这可以通过网格搜索等方法实现，类似于使用指数平滑法预测单变量时间序列数据中描述的方法。

1. 在训练集life_train上训练 ARIMA 模型，并查看模型的摘要。重要的是不要使用之前差分过的life_df1，因为 ARIMA 内部会根据d参数的值进行差分。在这个例子中，由于一阶差分（d=1）足以去趋势并使数据平稳，你将设置模型初始化中的d=1：

model = ARIMA(life_train, order=(0,1,1))
results = model.fit()
results.summary()

你将看到如下摘要：

图 10.22：ARIMA(0,1,1)模型在预期寿命数据上的摘要

注意，在模型摘要中提供了 AIC 和 BIC 评分。虽然这些评分很有用，但它们在用于比较多个模型时最有意义，因为它们有助于评估模型拟合度，同时惩罚过度复杂的模型。

在这个例子中，ARIMA 模型主要是一个 MA 过程，差分因子（d=1），摘要结果仅提供了 MA(1)组件的系数值。有关更多信息，请参阅*如何工作…*部分。

1. 你需要验证模型的残差，以确定 ARIMA(0, 1, 1)模型是否充分捕捉了时间序列中的信息。你会假设模型预测的残差是随机的（噪声），并且不遵循任何模式。更具体地说，你期望残差中没有自相关。你可以从acorr_ljungbox检验开始，然后查看残差的自相关函数（ACF）图。如果模型效果良好，你应该不期望看到任何自相关：

(acorr_ljungbox(results.resid,
                lags=25,
                return_df=True) < 0.05)['lb_pvalue'].sum()
>> 0

结果显示0，这是前 25 个滞后的结果的汇总，表示没有自相关。

也可以尝试查看 ACF 图：

plot_acf(results.resid, zero=False);

这应该生成一个 ACF 图。在这里，你应该期望图表不会显示出显著的滞后。换句话说，所有垂直线应该更接近零，或者对于所有滞后都应该为零：

图 10.23：显示残差没有自相关的 ACF 图

该图表确认了没有自相关的迹象（视觉上）。

1. 你还可以检查残差的分布。例如，你应该期望残差符合正态分布，均值为零。你可以使用 QQ 图和核密度估计（KDE）图来观察分布并评估正态性。你可以使用plot_diagnostics方法来实现：

results.plot_diagnostics();plt.show()

上面的代码将生成以下图表：

图 10.24：ARIMA(0,1,1)模型的视觉诊断

这些图表显示了与正态分布的轻微偏差。例如，一个完美的正态分布数据集将具有完美的钟形 KDE 图，并且所有点将在 QQ 图中完美对齐在直线上。

到目前为止，结果和诊断表明模型表现良好，尽管可能还有改进的空间。记住，构建 ARIMA 模型通常是一个迭代过程，涉及多轮测试和调整，以达到最佳结果。

1. ARIMA 建模过程中的最后一步是预测未来的数值，并将这些预测与测试数据集进行比较，测试数据集代表的是未见过的或超出样本的数据。使用plot_forecast()函数，这个函数是在本章技术要求部分中创建的：

plot_forecast(results, '1998', life_train, life_test)

执行此函数将生成一个图表，其中 x 轴从 1998 年开始。图表将显示三条线：实际数据分为两条线，一条是训练数据（orig_train），另一条是测试数据（orig_test），还有一条是forecast（预测值）。

这种视觉表现有助于评估 ARIMA 模型是否成功捕捉了数据集的潜在模式，以及它对未来数值的预测有多准确。

图 10.25：ARIMA(0,1,1)预测与实际预期寿命数据的对比

虚线（预测）似乎与预期的趋势不符，这与图 10.6 中的指数平滑模型结果不同，后者表现更好。为了解决这个问题，你可以运行多个具有不同(p, d, q)值的 ARIMA 模型，并比较它们的 RMSE、MAPE、AIC 或 BIC 分数，从而选择最佳拟合模型。你将在*后续内容...*部分探索这一选项。

它是如何工作的…

自回归模型或 AR(p) 是一个线性模型，利用前几个时间步的观测值作为回归方程的输入，来确定下一个步骤的预测值。因此，自回归中的auto部分表示自身，可以描述为一个变量对其过去版本的回归。典型的线性回归模型将具有以下方程：

这里，是预测变量，是截距，是特征或独立变量，是每个独立变量的系数。在回归分析中，你的目标是求解这些系数，包括截距（可以把它们当作权重），因为它们稍后会用于做预测。误差项，，表示残差或噪声（模型中未解释的部分）。

将其与自回归方程进行比较，你将看到相似之处：

这是一个阶数为p的 AR 模型，表示为AR(p)。自回归模型与回归模型的主要区别在于预测变量是，它是当前时刻的值，，而这些变量是的滞后（前）版本。在这个示例中，你使用了 ARIMA(0,1,1)，即 AR(0)，表示没有使用自回归模型。

与使用过去值的自回归模型不同，移动平均模型或 MA(q) 使用过去的误差（来自过去的估计）来进行预测：

将 AR(p) 和 MA(q) 模型结合会生成一个 ARMA(p,q) 模型（自回归滑动平均模型）。AR 和 ARMA 过程都假设时间序列是平稳的。然而，假设时间序列由于趋势的存在而不是平稳的，在这种情况下，除非进行一些转换，例如差分，否则不能在非平稳数据上使用 AR 或 ARMA 模型。这正是life数据的情况。

差分就是将当前值减去其前一个（滞后的）值。例如，差分阶数为一（lag=1）可以表示为。在 pandas 中，你使用了diff方法，默认设置为periods=1。

ARIMA 模型通过添加一个集成（差分）因子来使时间序列平稳，从而改进了 ARMA 模型。

你利用了 ACF 图和 PACF 图来估计 AR 和 MA 模型的阶数。自相关函数衡量当前观测值与其滞后版本之间的相关性。ACF 图的目的是确定过去的观测值在预测中的可靠性。

另一方面，偏自相关函数（PACF）类似于自相关，但去除了干预观测值之间的关系。

ACF 与 PACF 通过示例对比

如果在滞后 1、2、3 和 4 处存在强相关性，这意味着滞后 1 的相关性度量会受到与滞后 2 的相关性的影响，滞后 2 又会受到与滞后 3 的相关性的影响，依此类推。

在滞后 1 处的自相关函数（ACF）度量将包括这些先前滞后的影响，如果它们是相关的。相比之下，在滞后 1 处的偏自相关函数（PACF）会去除这些影响，以测量当前观察值与滞后 1 之间的纯关系。

ARIMA 之所以流行的原因之一是它能够推广到其他更简单的模型，如下所示：

• ARIMA(1, 0, 0) 是一个一阶自回归或 AR(1)模型

• ARIMA(1, 1, 0) 是一个差分的一级自回归模型

• ARIMA(0, 0, 1) 是一个一阶移动平均或 MA(1)模型

• ARIMA(1, 0, 1) 是一个 ARMA(1,1)模型

• ARIMA(0, 1, 1) 是一个简单的指数平滑模型

还有更多……

有时候，很难确定时间序列是 MA 过程还是 AR 过程，或者确定p或q的最优阶数（滞后值）。你可以参考下面这个例子，使用一种朴素的网格搜索方法，通过尝试不同的p，d和q组合来训练其他 ARIMA 模型，然后再选择最优模型。

在这里，你将利用你在技术要求部分创建的combinator()函数。你将训练多个 ARIMA 模型，并使用get_top_models_df()来找到最佳模型。作为起点，尝试对三个超参数（p，d，q）进行(0,1,2)组合。你将测试 3x3x3 或 27 个 ARIMA 模型：

pv, dv, qv = [list(range(3))]*3
vals = combinator([pv, dv, qv ])
score = {}
for i, (p, d, q) in enumerate(vals):
    m = ARIMA(life_train, order=(p,d,q))
    res = m.fit()
    y = life_train.values.ravel()
    y_hat = res.forecast(steps=len(y))
    score[i] = {'order': (p,d,q),
                'AIC':res.aic,
                'RMSPE': rmspe(y, y_hat),
                'BIC': res.bic,
                'AICc':res.aicc,
                'RMSE' : rmse(y, y_hat),
                'MAPE' : mape(y, y_hat),
                'model': res}
get_top_models_df(score, 'AIC')

这应该生成一个按 AIC 排序的 DataFrame。以下表格展示了前五个模型：

图 10.26：根据 AIC 得分排序的 27 个 ARIMA 模型的结果

你可以使用以下方法选择最佳模型：

best_m = get_top_models_df(score, 'AIC').iloc[0,-1]

如果你运行best_m.summary()来查看模型摘要，你会注意到它是一个ARIMA(0,2, 2)模型。这进一步确认了我们之前的观察，这实际上是一个移动平均（MA）过程，但我们错过了阶数。

赤池信息量准则（AIC）是一种衡量模型最大似然估计和模型简洁性之间平衡的指标。过于复杂的模型有时可能会过拟合，意味着它们看起来像是学习了，但一旦面对未见过的数据时，它们的表现会很差。AIC 得分随着参数数量的增加而受到惩罚，因为它们增加了模型的复杂性：

在这里，2k 被视为惩罚项。

贝叶斯信息准则（BIC）与 AIC 非常相似，但在模型的复杂性上有更高的惩罚项。通常，BIC 的惩罚项更大，因此它倾向于鼓励具有较少参数的模型，相较于 AIC。因此，如果你将排序或评估标准从 AIC 更改为 BIC，你可能会看到不同的结果。BIC 更偏向选择简单的模型：

在这里， 是最大似然估计，k 是估计的参数数量，n 是数据点的数量。

若要使用最佳模型绘制预测图表，可以运行以下命令：

plot_forecast(best_m, '1998', life_train, life_test);

这将产生以下图表：

图 10.26：ARIMA(0,2,2) 预测与实际的预期寿命数据对比

比较 图 10.25 中 ARIMA(0, 1, 1) 模型的输出和 图 10.26 中 ARIMA(0,2,2) 模型的输出。它们与图 10.12 中使用指数平滑法的结果相比如何？

在继续下一个示例之前，我们来对牛奶生产数据应用 ARIMA 模型，该数据与预期寿命数据不同，具有趋势和季节性：

pv, dv, qv = [list(range(3))]*3
vals = combinator([pv, dv, qv])
score = {}
for i, (p, d, q) in enumerate(vals):
    m = ARIMA(milk_train, order=(p,d,q))
    res = m.fit()
    y = milk_test.values.ravel()
    y_hat = res.forecast(steps=len(y))
    score[i] = {'order': (p,d,q),
                'AIC':res.aic,
                'BIC': res.bic,
                'AICc':res.aicc,
                'RMSPE': rmspe(y, y_hat),
                'RMSE' : rmse(y, y_hat),
                'MAPE' : mape(y, y_hat),
                'model': res}
model = get_top_models_df(score, 'AIC').iloc[0,-1]
plot_forecast(model, '1971', milk_train, milk_test);

运行代码并检查 model.summary() 输出的顶级模型后，你会发现它识别为 ARIMA(2,2,2)。结果的预测图可能整体上表现较差，如图 10.27 所示——预测与实际牛奶生产数据对比。

图 10.27：ARIMA(2,2,2) 预测与实际的牛奶生产数据对比

这个结果是预期中的，因为标准的 ARIMA 模型并不设计用来处理季节性。下一个示例将介绍 SARIMA（季节性 ARIMA），它更适合建模具有季节模式和趋势的时间序列数据。

另见

若要了解更多关于 ARIMA 类的信息，你可以访问 statsmodels 的官方文档 www.statsmodels.org/dev/generated/statsmodels.tsa.arima.model.ARIMA.html。

那么，具有趋势和季节性的 milk 数据怎么样呢？下一个示例将探索如何使用 SARIMA 模型处理此类数据。

FORECAST 与 PREDICT 方法的区别

在 plot_forecast 函数中，我们使用了 forecast 方法。在 statsmodels 中，SARIMA 模型家族（如 ARMA 和 ARIMA）有两种方法可以进行预测：predict 和 forecast。

predict 方法允许你同时进行样本内（历史）和样本外（未来）的预测，因此该方法需要 start 和 end 参数。另一方面，forecast 方法仅接受 steps 参数，表示样本外预测的步数，从样本或训练集的末尾开始。

使用季节性 ARIMA 进行单变量时间序列数据预测

在本示例中，你将接触到一种增强型 ARIMA 模型，用于处理季节性，称为季节性自回归积分滑动平均模型或 SARIMA。与 ARIMA(p, d, q) 类似，SARIMA 模型也需要(p, d, q)来表示非季节性阶数。此外，SARIMA 模型还需要季节性成分的阶数，表示为(P, D, Q, s)。将两者结合，模型可以写成 SARIMA(p, d, q)(P, D, Q, s)。字母含义保持不变，字母的大小写表示不同的成分。例如，小写字母代表非季节性阶数，而大写字母代表季节性阶数。新参数s表示每个周期的步数——例如，对于月度数据，s=12，对于季度数据，s=4。

在 statsmodels 中，你将使用 SARIMAX 类来构建 SARIMA 模型。

在本示例中，你将使用milk数据，该数据包含趋势性和季节性成分。此数据已在技术要求部分准备好。

如何实现…

按照以下步骤操作：

1. 首先导入必要的库：

from statsmodels.tsa.statespace.sarimax import SARIMAX

1. 从图 10.1中，我们确定了季节性和趋势性都存在。我们还可以看到季节性效应是加性的。一个季节的周期或期数为 12，因为数据是按月收集的。这可以通过 ACF 图进行确认：

plot_acf(milk, lags=40, zero=False);

这应该会生成一个milk数据集的 ACF 图，并在特定滞后处显示明显的周期性波动：

图 10.28: ACF 图显示滞后 1、12 和 24 时有显著的峰值

注意，每 12 个月（滞后）会有一个重复的模式。如果这种模式不容易发现，可以尝试在差分数据后查看 ACF 图——例如，先对数据进行去趋势（一次差分），然后再绘制 ACF 图：

plot_acf(milk.diff(1).dropna(), lags=40, zero=False);

这应该会生成一个差分数据的 ACF 图，使季节性峰值更加明显：

图 10.29 – 差分后的 ACF 图显示在滞后 1、12、24 和 36 时有显著的峰值

你也可以提取季节性成分并用其生成 ACF 图，如下代码所示：

decomposed = seasonal_decompose(milk, period=12, model='multiplicative')
milk_s = decomposed.seasonal
plot_acf(milk_s, zero=False, lags=40);

ACF 图将显示经过分解后的季节性成分的自相关，并会讲述与图 10.28和图 10.29类似的故事。

图 10.30: 分解后的季节性成分的 ACF 图显示在滞后 1、12、24 和 36 时有显著的峰值

通常，处理月度数据时可以假设一个 12 个月的周期。例如，对于非季节性 ARIMA 部分，可以从d=1开始去趋势化，对于季节性 ARIMA 部分，可以将 D=1 作为起始值，前提是 s=12。

1. 假设你不确定d（非季节性差分）和D（季节性差分）的值，在这种情况下，你可以在差分后使用check_stationarity函数来判断季节性差分是否足够。通常，如果时间序列既有趋势又有季节性，你可能需要进行两次差分。首先进行季节性差分，然后进行一阶差分以去除趋势。

通过使用diff(12)进行季节性差分（去季节化），并测试这是否足以使时间序列平稳。如果不够，那么你需要继续进行一阶差分diff()：

milk_dif_12 = milk.diff(12).dropna()
milk_dif_12_1 = milk.diff(12).diff(1).dropna()
sets = [milk, milk_dif_12, milk_dif_12_1]
desc = ['Original', 'Deseasonalize (Difference Once)', 'Differencing Twice']
fig, ax = plt.subplots(2,2, figsize=(20,10))
index, l = milk.index, milk.shape[0]
for i, (d_set, d_desc) in enumerate(zip(sets, desc)):
    v, r = i // 2, i % 2
    outcome, pval = check_stationarity(d_set)
    d_set.plot(ax= ax[v,r], title=f'{d_desc}: {outcome}', legend=False)
    pd.Series(d_set.mean().values.tolist()*l, index=index).plot(ax=ax[v,r])
    ax[v,r].title.set_size(20)
ax[1,1].set_visible(False)
plt.show()

这应该会生成 2x2 子图（每行两个图），其中额外的子图是隐藏的：

图 10.31：原始数据、季节差分数据和两次差分数据的平稳性比较

1. 现在，你需要估计非季节性（p, q）和季节性（P, Q）分量的 AR 和 MA 阶数。为此，你必须使用平稳数据上的 ACF 和 PACF 图，这些图可以在milk_dif_12_1数据框中找到：

fig, ax = plt.subplots(1,2)
plot_acf(milk_dif_12_1, zero=False, lags=36, ax=ax[0], title=f'ACF - {d_desc}')
plot_pacf(milk_dif_12_1, zero=False, lags=36, ax=ax[1], title=f'PACF - {d_desc}')
plt.show()

这应该会在同一行中生成 ACF 和 PACF 图：

图 10.32：乳制品数据变为平稳后的 ACF 和 PACF 图

从自相关函数（ACF）图开始，可以看到在滞后 1 处有一个显著的峰值，表示 MA 过程的非季节性阶数。滞后 12 处的峰值表示 MA 过程的季节性阶数。注意，在滞后 1 之后有一个截断，接着是滞后 12 的峰值，然后又是另一个截断（之后没有其他显著的滞后）。这些都是移动平均模型的指示——更具体地说，是q=1和Q=1的阶数。

PACF 图也确认了这一点；在滞后 12、24 和 36 处的指数衰减表明是 MA 模型。在这里，季节性 ARIMA 将是 ARIMA(0, 1,1)(0, 1, 1, 12)。

1. 根据最初提取的 AR 和 MA 阶数信息构建 SARIMA 模型。以下代码将在训练数据集上拟合一个 SARIMA(0, 1, 1)(0, 1, 1, 12)模型。请注意，结果可能与绘制 ACF 和 PACF中的结果有所不同，因为在该配方中数据没有被拆分，但在这里已经拆分：

sarima_model = SARIMAX(milk_train,
                       order=(0,1,1),
                       seasonal_order=(0,1,1,12))
model = sarima_model.fit(disp=0)

1. 现在，使用plot_diagnostics方法，它将在拟合模型后可用：

model.plot_diagnostics(figsize=(15,7));

这将提供四个图——一个标准化残差图，一个 QQ 图，一个 ACF 残差图和一个带有核密度图的直方图：

图 10.33：SARIMA(0,1,1)(0,1,1,12)诊断图

残差的 ACF 图（相关图）未显示自相关（忽略滞后 0 处的尖峰，因为它总是 1）。然而，直方图和 QQ 图显示残差并不符合完美的正态分布。与随机残差（无自相关）相比，这些假设并不关键。总体来说，结果非常有前景。

你可以使用summary方法获取摘要：

model.summary()

这应该以表格格式打印出关于模型的附加信息，包括季节性和非季节性成分的系数。回想一下，诊断图是基于标准化残差的。你可以绘制残差的 ACF 图而不进行标准化：

plot_acf(model.resid[1:])

这应该生成如下的 ACF 图，其中显示了阈值之外的一些滞后，表明存在自相关，并且有改进的潜力。

图 10.34：SARIMA(0,1,1)(0,1,1,12)的残差 ACF 图

1. 使用plot_forecast函数绘制 SARIMA 模型的预测，并与测试集进行比较：

plot_forecast(model, '1971', milk_train, milk_test);

这应该生成一个图表，x 轴从 1971 年开始：

图 10.35：使用 SARIMA(0,1,1)(0,1,1,12)的牛奶生产预测与实际生产对比

总体而言，SARIMA 模型在捕捉季节性和趋势效应方面表现得相当不错。你可以通过评估使用其他度量标准（例如 RMSE、MAPE 或 AIC 等）来迭代并测试(p, q)和(P, Q)的不同值。更多内容请参见*还有更多……*部分。

它是如何工作的……

SARIMA 模型通过加入季节性来扩展 ARIMA 模型，因此需要额外的一组季节性参数。例如，一个没有季节性成分的非季节性顺序为(1, 1, 1)的 SARIMA 模型会被指定为 SARIMA(1,1,1)(0,0,0,0)，这本质上简化为 ARIMA(1, 1, 1)模型。为了处理季节性，你需要将季节性顺序设置为非零值。

在 statsmodels 中，SARIMAX类将 AR、MA、ARMA、ARIMA 和 SARIMA 模型进行了一般化，使你能够拟合适合自己时间序列数据的模型，无论其是否具有季节性成分。同样，正如在使用指数平滑法预测单变量时间序列数据一节中所讨论的，ExponentialSmoothing类作为SimpleExpSmoothing和Holt模型的广义实现。

还有更多……

类似于使用 ARIMA 预测单变量时间序列数据一节中的方法，你可以执行一个朴素的网格搜索，以评估不同的非季节性(p, d, q)和季节性(P, D, Q, s)参数组合，从而确定最佳的 SARIMA 模型。

这可以通过利用combinator()函数来实现，遍历所有可能的参数组合，在每次迭代中拟合一个 SARIMA 模型。要确定前 N 个模型，可以使用get_top_models_df函数。

例如，考虑测试所有组合，其中非季节性（p, d, q）参数分别取值（0,1,2），季节性（P, D, Q）参数分别取值（0,1），同时将s（季节周期）保持为 12。这种设置将测试总共（3x3x3x2x2x2）= 216 个 SARIMA 模型。虽然这种暴力（朴素）方法在计算上可能非常密集，但它仍然是有效的方法。像Auto_ARIMA这样的自动时间序列库通常也采用类似的穷举网格搜索来优化模型参数：

P_ns, D_ns, Q_ns = [list(range(3))]*3
P_s, D_s, Q_s = [list(range(2))]*3
vals = combinator([P_ns, D_ns, Q_ns, P_s, D_s, Q_s])
score = {}
for i, (p, d, q, P, D, Q) in enumerate(vals):
    if i%15 == 0:
        print(f'Running model #{i} using SARIMA({p},{d},{q})({P},{D},{Q},12)')
    m = SARIMAX(milk_train,
                order=(p,d,q),
                seasonal_order=(P, D, Q, 12),
                enforce_stationarity=False)
    res = m.fit(disp=0)
    y = milk_test.values.ravel()
    y_hat = res.forecast(steps=len(y))
    score[i] = {'non-seasonal order': (p,d,q),
                'seasonal order': (P, D, Q),
                'AIC':res.aic,
                'AICc': res.aicc,
                'BIC': res.bic,
                'RMSPE': rmspe(y, y_hat),
                'RMSE' : rmse(y, y_hat),
                'MAPE' : mape(y, y_hat),
                'model': res}

执行前面的代码后，它应该每 15 次迭代打印一次状态输出，如下所示：

Running model #0 using SARIMA(0,0,0)(0,0,0,12)
Running model #15 using SARIMA(0,0,1)(1,1,1,12)
Running model #30 using SARIMA(0,1,0)(1,1,0,12)
Running model #45 using SARIMA(0,1,2)(1,0,1,12)
Running model #60 using SARIMA(0,2,1)(1,0,0,12)
Running model #75 using SARIMA(1,0,0)(0,1,1,12)
Running model #90 using SARIMA(1,0,2)(0,1,0,12)
Running model #105 using SARIMA(1,1,1)(0,0,1,12)
Running model #120 using SARIMA(1,2,0)(0,0,0,12)
Running model #135 using SARIMA(1,2,1)(1,1,1,12)
Running model #150 using SARIMA(2,0,0)(1,1,0,12)
Running model #165 using SARIMA(2,0,2)(1,0,1,12)
Running model #180 using SARIMA(2,1,1)(1,0,0,12)
Running model #195 using SARIMA(2,2,0)(0,1,1,12)
Running model #210 using SARIMA(2,2,2)(0,1,0,12)

请注意enforce_stationarity=False参数，以避免在进行初始网格搜索时可能出现的LinAlgError。

要识别按 AIC 排序的前 5 个模型，可以运行get_top_models_df函数：

get_top_models_df(score, 'AIC')

图 10.36：按 AIC 排名的前 5 个 SARIMA 模型，针对牛奶生产数据

注意到前两个模型的 AIC 得分相似。通常，当两个模型的 AIC 得分相似时，较简单的模型更受偏好。例如，SARIMA(0,2,2)(0,1,1)模型比 SARIMA(2,2,2)(0,1,1)模型更简单。还有其他指标，如 AICc 和 BIC，在这种情况下，它们会更倾向于选择第二个模型（model_id=67），而不是第一个模型（model_id=211）。类似地，如果考虑 RMSPE 或 MAPE，你可能会选择不同的模型。

奥卡姆剃刀

奥卡姆剃刀原理表明，当多个模型产生相似的质量，换句话说，能同样良好地拟合数据时，应该偏好较简单的模型。这一原理在评估多个模型时尤其有用，比如在评估多个 SARIMA 模型时。如果一些模型脱颖而出成为强有力的候选者，通常来说，简单的模型更受青睐，假设最初这些模型是等可能的候选者。

要根据 BIC 得分对模型进行排序，请重新运行该函数：

get_top_models_df(score, 'BIC')

请查看下图中的结果（图 10.36），它展示了按 BIC 排名的前 5 个 SARIMA 模型，针对牛奶生产数据。

图 10.37：按 BIC 排名的前 5 个 SARIMA 模型，针对牛奶生产数据

比较图 10.36 和图 10.35 中的结果，你会注意到图 10.36 中的第三个模型（model_id = 35），即之前推导出来的 SARIMA(0,1,1)(0,1,1)模型。

现在，让我们根据 BIC 得分选择最佳模型：

best_model = get_top_models_df(score, BIC).iloc[0,-1]

最后，你可以使用plot_forecast函数将模型的预测结果与实际数据一起可视化：

plot_forecast(best_model, '1962', milk_train, milk_test);

这应该生成一个图表，x 轴从 1962 年开始，如图 10.37 所示：

图 10.38：使用 SARIMA(0,2,2)(0,1,1,12)的牛奶生产预测与实际生产对比

参见

要了解更多关于 SARIMAX 类的信息，可以访问 statsmodels 的官方文档：www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html。

使用auto_arima进行单变量时间序列预测

在这个示例中，你需要安装pmdarima，这是一个 Python 库，其中包括auto_arima——一个旨在自动化优化和拟合 ARIMA 模型的工具。Python 中的auto_arima实现灵感来源于 R 中forecast包中的流行auto.arima。

正如你在前面的示例中看到的，确定 AR 和 MA 组件的正确阶数可能会很具挑战性。虽然像检查 ACF 和 PACF 图这样的技术很有帮助，但找到最优模型通常需要训练多个模型——这一过程称为超参数调优，可能非常耗时费力。这正是auto_arima的亮点，它简化了这个过程。

与简单粗暴的手动网格搜索逐一尝试每种参数组合的方法不同，auto_arima使用了一种更高效的方法来寻找最优参数。auto_arima函数采用了一种逐步 算法，其速度更快、效率更高，优于完整的网格搜索或随机搜索：

• 默认情况下，使用stepwise=True时，auto_arima会进行逐步搜索，逐步优化模型参数。

• 如果你设置stepwise=False，它将执行全面的“蛮力”网格搜索，遍历所有参数组合。

• 使用random=True时，它会执行随机搜索。

逐步 算法由 Rob Hyndman 和 Yeasmin Khandakar 于 2008 年在论文Automatic Time Series Forecasting: The forecast Package for R中提出，该论文发表在《统计软件杂志》27 卷 3 期（2008 年）（doi.org/10.18637/jss.v027.i03）。简而言之，逐步是一种优化技术，它更高效地利用了网格搜索。这是通过使用单位根检验和最小化信息准则（例如，赤池信息量准则（AIC）和最大似然估计（MLE）来实现的）。

此外，auto_arima还可以处理季节性和非季节性的 ARIMA 模型。对于季节性模型，设置seasonal=True以启用季节性参数（P、D、Q）的优化。

准备工作

在继续此示例之前，你需要安装pmdarima。

要使用pip安装，可以使用以下命令：

pip install pmdarima

要使用conda安装，可以使用以下命令：

conda install -c conda-forge pmdarima

你将使用在技术要求部分准备的“Milk Production”数据集。你将使用milk_train进行训练，使用milk_test进行评估。请回忆一下，数据包含了趋势和季节性，因此你将训练一个SARIMA模型。

如何做…

pmdarima库封装了statsmodels库，因此你会遇到熟悉的方法和属性。你将遵循类似的流程，首先加载数据，将数据拆分为训练集和测试集，训练模型，然后评估结果。这些步骤已在技术要求部分完成。

1. 首先导入pmdarima库

import pmdarima as pm

1. 使用pmdarima中的auto_arima函数来找到 SARIMA 模型的最佳配置。对 Milk Production 数据集的先验知识是从auto_arima获得最佳结果的关键。你知道数据存在季节性模式，因此需要为两个参数提供值：seasonal=True和m=12，其中m代表季节中的周期数。如果没有设置这些参数（seasonal和m），搜索将仅限于非季节性阶数(p, d, q)。

test参数指定用于检测平稳性并确定差分阶数(d)的单位根检验类型。默认的检验是kpss。你将把参数改为使用adf（以保持与之前配方一致）。类似地，seasonal_test用于确定季节性差分的阶数(D)。默认的seasonal_test是OCSB，你将保持不变：

auto_model = pm.auto_arima(milk_train,
                           seasonal=True,
                           m=12,
                           test='adf',
                           stepwise=True)
auto_model.summary()

摘要提供了所选 SARIMA 模型的详细配置，包括信息准则得分，如 AIC 和 BIC：

图 10.39：使用 auto_arima 选择的最佳 SARIMA 模型摘要

有趣的是，所选模型 SARIMA(0,1,1)(0,1,1,12)与在使用季节性 ARIMA 预测单变量时间序列数据配方中推导出的模型一致，在该配方中，你通过 ACF 和 PACF 图估算了非季节性阶数(p, q)和季节性阶数(P, Q)。

1. 为了监控在逐步搜索过程中评估的每个模型配置的性能，在auto.arima函数中启用trace=True参数：

auto_model = pm.auto_arima(milk_train,
                           seasonal=True,
                           m=12,
                           test='adf',
                           stepwise=True,
                           trace=True)

该设置将打印由逐步算法测试的每个 SARIMA 模型的 AIC 结果，如图 10.39 所示：

图 10.40：auto_arima 基于 AIC 评估不同的 SARIMA 模型

最佳模型是基于 AIC 选择的，AIC 由information_criterion参数决定。默认情况下，设置为aic，但可以更改为其他支持的准则之一：bic、hqic或oob。

在图 10.39中，两个突出的模型具有相似的 AIC 分数，但非季节性（p，q）阶数大相径庭。优选模型（标记为数字 1）缺少非季节性自回归 AR(p)成分，而是依赖于移动平均 MA(q)过程。相反，第二个突出模型（标记为数字 2）仅包含非季节性成分的 AR(p)过程。这表明，尽管auto_arima显著有助于模型选择，但仍需小心判断和分析，以有效解释和评估结果。

若要探索信息准则的选择如何影响模型选择，请将information_criterion更改为bic并重新运行代码：

auto_model = pm.auto_arima(milk_train,
                           seasonal=True,
                           m=12,
                           test='adf',
                           information_criterion='bic',
                           stepwise=True,
                           trace=True)

图 10.41：auto_arima 根据 BIC 评估不同的 SARIMA 模型

如图 10.40 所示，这将基于 BIC 从每次迭代中生成输出。值得注意的是，最终选择的模型与基于 AIC 从图 10.39 选择的模型相同。然而，请注意，第二个模型（标记为数字 2），尽管在图 10.39 中是一个接近的竞争者，但在 BIC 标准下已不再具有竞争力。

1. 使用plot_diagnostics方法评估模型的整体表现。这是你之前在 statsmodels 中使用过的相同方法。

auto_model.plot_diagnostics(figsize=(15,7));

这应该会生成选定的 SARIMA 模型的残差分析诊断图，如图 10.41 所示：

图 10.42：基于选定的 SARIMA 模型的残差分析诊断图

若要访问模型摘要，请使用summary方法。这将生成图 10.42，显示由auto_arima选择的 SARIMA 模型摘要：

图 10.43：基于 auto_arima 选择模型的 SARIMA 模型摘要

1. 若要预测未来的周期，请使用predict方法。你需要指定要预测的周期数：

n = milk_test.shape[0]
index = milk_test.index
ax = milk_test.plot(style='--', alpha=0.6, figsize=(12,4))
pd.Series(auto_model.predict(n_periods=n),
          index=index).plot(style='-', ax=ax)
plt.legend(['test', 'forecast']);

生成的图形，如图 10.43 所示，比较了预测值与测试数据。

图 10.44：将 auto_arima 的预测与实际测试数据进行对比

你可以通过将return_conf_int参数从False更新为True来获得预测的置信区间。这将允许你使用 matplotlib 的fill_between函数绘制上下置信区间。默认的置信区间设置为 95%（alpha=0.05）。

以下代码使用predict方法，返回置信区间，并将预测值与测试集进行对比：

n = milk_test.shape[0]
forecast, conf_interval = auto_model.predict(n_periods=n,
                                             return_conf_int=True,
                                            alpha=0.05)
lower_ci, upper_ci  = zip(*conf_interval)
index = milk_test.index
ax = milk_test.plot(style='--', alpha=0.6, figsize=(12,4))
pd.Series(forecast, index=index).plot(style='-', ax=ax)
plt.fill_between(index, lower_ci, upper_ci, alpha=0.2)
plt.legend(['test', 'forecast']);

这应该生成一个图，带有阴影区域，表示实际值落在此范围内的可能性。理想情况下，你会更倾向于选择较窄的置信区间范围。

阴影区域基于图 10.44 中所示的置信区间的上下界：

图 10.45：将来自 auto_arima 的预测与实际测试数据及置信区间进行对比的图示

请注意，预测线位于阴影区域的中间位置，表示上下界的均值。

以下代码检查预测值是否与置信区间的均值一致：

sum(forecast) == sum(conf_interval.mean(axis=1))
>> True

工作原理…

pmdarima 库中的 auto_arima 函数作为 statsmodels SARIMAX 类的封装器，旨在自动化识别最佳模型和参数的过程。此函数提供了三种主要方法来控制训练过程的优化，这些方法由 stepwise 和 random 参数决定：

• 朴素暴力网格搜索：通过设置 stepwise=False 和 random=False，对所有参数组合进行全面的网格搜索。该方法对每种组合进行详尽评估，虽然耗时，但非常彻底。

• 随机网格搜索：通过设置 stepwise=False 和 random=True 来激活参数空间中的随机搜索。这种方法随机选择组合进行测试，尤其在大参数空间中，虽然是随机的，但仍然有效，且速度较快。

• 逐步搜索算法：默认启用 stepwise=True，此方法使用启发式算法逐步探索参数空间，基于前一步的 AIC 或 BIC 分数，逐个添加或减少每个参数。与全面网格搜索相比，它通常更快、更高效，因为它智能地缩小了模型范围，聚焦于最可能提供最佳拟合的模型。

还有更多…

pmdarima 库提供了大量有用的函数，帮助你做出明智的决策，使你能够更好地理解你所处理的数据。

例如，ndiffs 函数执行平稳性检验，以确定差分阶数 d，使时间序列变得平稳。检验包括增广迪基-富勒（adf）检验、克维亚茨科夫-菲利普斯-施密特-辛（kpss）检验和菲利普斯-佩龙（pp）检验。

同样，nsdiffs 函数有助于估算所需的季节性差分阶数（D）。该实现包括两个检验——奥斯本-崔-史密斯-比尔钦霍尔（ocsb）检验和卡诺瓦-汉森（ch）检验：

from pmdarima.arima.utils import ndiffs, nsdiffs
n_adf = ndiffs(milk, test='adf')
# KPSS test (the default in auto_arima):
n_kpss = ndiffs(milk, test='kpss')
n_pp = ndiffs(milk, test='pp') 
n_ch = nsdiffs(milk, test='ocsb', m=10, max_D=12,)
n_ocsb = nsdiffs(milk, test='ch' , m=10, max_D=12,)

auto_arima函数通过设置各种参数的最小和最大约束，允许对模型评估过程进行详细控制。例如，你可以指定非季节性自回归阶数p或季节性移动平均Q的限制。以下代码示例演示了如何设置一些参数和约束条件：

model = pm.auto_arima(milk_train,
                      seasonal=True,
                      with_intercept=True,
                      d=1,
                      max_d=2,
                      start_p=0, max_p=2,
                      start_q=0, max_q=2,
                      m=12,
                      D=1,
                      max_D=2,
                      start_P=0, max_P=2,
                      start_Q=0, max_Q=2,
                      information_criterion='aic',
                      stepwise=False,
                      out_of_sample_siz=25,
                      test = 'kpss',
                      score='mape',
                      trace=True) 

如果你运行前面的代码，auto_arima将根据你提供的约束条件，为每个参数值的组合创建不同的模型。由于stepwise被设置为False，这就变成了一种暴力的网格搜索，其中每种不同变量组合的排列都会逐一测试。因此，这通常是一个比较慢的过程，但通过提供这些约束条件，你可以提高搜索性能。

通过启用trace=True，可以显示测试的每个模型配置的 AIC 得分。完成后，应该会打印出最佳模型。

这里采取的方法，设置stepwise=False，应该与你在使用季节性 ARIMA 进行单变量时间序列预测教程中采取的方法类似，该教程位于*更多内容...*部分。

使用 Darts 的 AutoArima

在使用指数平滑进行单变量时间序列预测的教程中，你在*更多内容...*部分介绍了Darts库。

Darts 库提供了AutoArima类，它是pmdarima的auto_arima的一个轻量级封装。以下代码演示了如何利用 Darts 执行相同的功能：

model = AutoARIMA(seasonal=True,
                           m=12,
                 stepwise=True)
ts = TimeSeries.from_dataframe(milk_train.reset_index(),
                                    time_col='month', value_cols='production', freq='MS')
darts_arima = model.fit(ts)
darts_forecast = model.predict(len(milk_test))
ts.plot(label='Training')
darts_forecast.plot(label='Forecast', linestyle='--');

这段代码生成的图表展示了预测结果，如图 10.45 所示：

图 10.46：使用 Darts 库的 AutoARIMA 进行预测

使用 Darts 的 StatsForecastAutoARIMA

Darts 还提供了一个Statsforecasts的包装器，封装了其 Auto_ARIMA，提供了比 AutoArima 可能更快的实现。以下代码演示了如何使用StatsForecastAutoARIMA来执行相同的功能：

from darts.models import StatsForecastAutoARIMA
model = StatsForecastAutoARIMA(season_length=12)
model.fit(ts)
pred = model.predict(len(milk_test))
ts.plot(label='Training')
darts_forecast.plot(label='AutoArima', linestyle='--');
pred.plot(label='StatsforecstsAutoArima');

这段代码生成的图表展示了 AutoARIMA 和 StatsForecastAutoARIMA 的预测对比，如图 10.47 所示：

图 10.47：使用 Darts 库的 StatsForecastAutoARIMA 进行预测

另见

若要了解有关auto_arima实现的更多信息，请访问官方文档：alkaline-ml.com/pmdarima/modules/generated/pmdarima.arima.auto_arima.html。

在下一个教程中，你将学习一种新的算法，该算法提供了一个更简单的 API 来进行模型调优和优化。换句话说，你需要关注的参数大大减少。