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

于 2025-07-20 01:07:52 发布
阅读量750 收藏 8
点赞数 22
CC 4.0 BY-SA版权
分类专栏： 默认分类 文章标签： 默认分类
License CC BY-NC-SA 4.0 / 自豪地采用谷歌翻译
本文链接：https://blog.youkuaiyun.com/wizardforcel/article/details/149473906

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

译者：飞龙

协议：CC BY-NC-SA 4.0

第十一章：11 种时间序列的附加统计建模技术

加入我们的 Discord 书籍社区

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file0.png

packt.link/zmkOY

在第十章，使用统计方法构建单变量时间序列模型中，你学习了流行的预测技术，如指数平滑法、非季节性 ARIMA，和季节性 ARIMA。这些方法通常被称为经典的统计预测方法，具有快速、易实现和易解释的特点。

在本章中，你将深入学习基于上一章基础的额外统计方法。本章将介绍一些可以自动化时间序列预测和模型优化的库——Facebook（Meta）的Prophet库。此外，你将探索statsmodels的向量自回归（VAR）类，用于处理多变量时间序列，以及arch库，它支持用于金融数据波动性建模的GARCH模型。

本章的主要目标是让你熟悉自动化预测工具（如 Prophet），并介绍多变量时间序列建模的概念，使用 VAR 模型。你还将了解如何在金融时间序列中建模和预测波动性，这对于风险管理和财务决策至关重要。

在本章中，我们将涵盖以下食谱：

  • 使用 Facebook Prophet 进行时间序列数据预测

  • 使用 VAR 预测多变量时间序列数据

  • 评估向量自回归（VAR）模型

  • 使用 GARCH 预测金融时间序列数据中的波动性

技术要求

你可以从本书的 GitHub 仓库下载 Jupyter Notebooks 以便跟随学习，并下载本章所需的数据集：

在本章的食谱中，你将使用一些常见的库。你可以通过以下代码提前导入它们：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')
plt.rc("figure", figsize=(16, 5))

使用 Facebook Prophet 进行时间序列数据预测

Prophet 库是一个流行的开源项目，最初由 Facebook（现为 Meta）开发，基于 2017 年的一篇论文，提出了一种时间序列预测算法，标题为《大规模预测》。该项目因其简单性、能够创建高效的预测模型、以及处理复杂季节性、假期效应、缺失数据和异常值的能力而广受欢迎。Prophet 自动化了许多设计预测模型的过程，同时提供了丰富的内置可视化功能。其他功能包括构建增长模型（如饱和预测）、处理趋势和季节性的 不确定性，以及检测变化点。

在本教程中，您将使用 Milk Production 数据集进行性能基准测试。这是 第十章，使用统计方法构建单变量时间序列模型 中介绍的相同数据集。使用相同的数据集有助于理解和比较不同的方法。

Prophet 是一个加性回归模型，可以处理非线性趋势，特别是在存在强烈季节性效应时。该模型将时间序列数据分解为三个主要组成部分：趋势、季节性和假期，形式如下所示：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file244.png

其中：

  • 是趋势函数，

  • 代表周期性季节性函数，

  • 考虑到假期的影响，

  • 是残差误差项。

Prophet 使用贝叶斯推断自动调整和优化模型组件。其背后依赖于Stan，一个最先进的贝叶斯建模平台，当前的 Python 接口为 cmdstand 和 cmdstanpy（取代了早期的 PyStan）。最近的更新提高了与 Apple M1/M2 芯片的兼容性，并增强了模型自定义选项，例如调整假期处理和缩放输出。

准备工作

您可以在本书的 GitHub 仓库中找到 Jupyter 笔记本和必要的数据集。有关更多信息，请参阅本章的 技术要求 部分。我们为本教程使用 Prophet 1.0 版本。

在安装像 Prophet 这样的新库时，创建一个新的 Python 环境总是一个好主意。如果您需要快速回顾如何创建虚拟 Python 环境，可以查看 第一章，时间序列分析入门 中的 开发环境设置教程。该章节介绍了两种方法：使用 conda 和 venv。

例如，您可以像以下示例一样使用 conda 创建环境：

conda create -n prophet python=3.11 -y

上述代码将创建一个名为 prophet 的新虚拟环境。要使新的 prophet 环境在 Jupyter 中可见，可以运行以下代码：

python -m ipykernel install --user --name prophet --display-name "Prophet"
conda activate prophet

一旦新环境被激活，您就可以安装 Prophet 库。要使用 pip 安装，您可以运行以下命令：

pip install prophet

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

conda install -c conda-forge prophet

如何做到这一点…

Prophet 要求输入数据为一个包含两个特定列的 pandas DataFrame：一个名为ds的datetime列和一个名为y的目标变量列——这是你希望预测的变量。Prophet 不能直接处理Datetimeindex。因此，确保这些列明确命名非常重要。如果数据包含超过两列，Prophet 只会识别ds和y，并默认忽略其他列。然而，如果你想添加额外的预测变量（回归器），可以使用add_regressor方法。如果这些列缺失或命名不正确，Prophet 会抛出错误。

为确保数据格式正确，请按照以下步骤操作：

  1. 首先读取milk_productions.csv文件并重命名列为ds和y：
from prophet import Prophet
milk_file = Path('../../datasets/Ch11/milk_production.csv')
milk = pd.read_csv(milk_file, parse_dates=['month'])
milk.columns = ['ds', 'y']
milk.info()
>>
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 168 entries, 0 to 167
Data columns (total 2 columns):
 #   Column  Non-Null Count  Dtype        
---  ------  --------------  -----        
 0   ds      168 non-null    datetime64[ns]
 1   y       168 non-null    int64        
dtypes: datetime64ns, int64(1)
memory usage: 2.8 KB
  1. 将数据拆分为测试集和训练集。我们使用90/10的比例进行拆分，代码如下：
idx = round(len(milk) * 0.90)
train = milk[:idx]
test = milk[idx:]
print(f'Train: {train.shape}')
print(f'Test: {test.shape}')
>>
Train: (151, 2)
Test: (17, 2)
  1. 你可以通过以下一行代码创建一个 Prophet 类的实例，并使用fit方法在训练集上进行拟合。牛奶生产的时间序列是按月记录的，具有趋势和稳定的季节性波动（加性）。Prophet 的默认seasonality_mode是additive，所以保持原样即可：
from prophet import Prophet
model = Prophet().fit(train)
  1. 在你使用模型进行预测之前，需要进行一些设置。使用make_future_dataframe方法将train DataFrame 扩展到特定的预测期数，并按指定频率生成数据：
future = m_milk.make_future_dataframe(len(test), freq='MS')

这会将训练数据延长 17 个月（test集中的期数）。总之，你应该拥有与牛奶 DataFrame（训练和测试）中相同数量的期数。频率设置为月初，即freq='MS'。future对象只包含一个列ds，其类型为datetime64[ns]，用于填充预测值：

len(milk) == len(future)
>> True
print(future.tail())
>>
            ds
163 1975-08-01
164 1975-09-01
165 1975-10-01
166 1975-11-01
167 1975-12-01
  1. 使用predict方法对future DataFrame 进行预测。结果将是一个与forecast长度相同的 DataFrame，但现在会增加一些额外的列：
forecast = model.predict(future)
forecast.columns.tolist()
>>
['ds',
 'trend',
 'yhat_lower',
 'yhat_upper',
 'trend_lower',
 'trend_upper',
 'additive_terms',
 'additive_terms_lower',
 'additive_terms_upper',
 'yearly',
 'yearly_lower',
 'yearly_upper',
 'multiplicative_terms',
 'multiplicative_terms_lower',
 'multiplicative_terms_upper',
 'yhat']

请注意，Prophet 返回了许多细节，帮助你了解模型的表现。重点是ds和预测值yhat。yhat_lower和yhat_upper分别表示预测的不确定性区间（yhat）。

  1. model对象提供了两个绘图方法：plot和plot_components。首先使用plot方法来可视化 Prophet 的预测结果：
model.plot(forecast,
           ylabel='Milk Production in Pounds',
           include_legend=True);

这应该会生成一个 Prophet 预测图：图中的点代表训练数据点，点上的线代表历史数据的估计预测，线条延伸到训练点之外，反映未来的预测。

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file245.png

图 11.4 – 使用 Prophet 绘制预测图（历史与未来）

如果你只想展示超出训练集的预测期，可以使用以下代码：

predicted = model.predict(test)
model.plot(predicted,
           ylabel='Milk Production in Pounds',
          include_legend=True);

这里你只预测了测试数据集的长度。预测方法将仅捕获 ds 列（日期时间）。这将生成一个类似于图 11.4所示的图表，但它只会显示未来数据点（训练数据点之后）的预测线——即未来的预测。

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file246.png

图 11.5 – 使用 Prophet 绘制预测（历史和未来）

图 11.4 中的阴影区域表示不确定性区间。这由 forecast 数据框中的 yhat_lower 和 yhat_upper 列表示。

  1. 接下来，重要的图表涉及预测组件。使用 plot_components 方法绘制这些组件：
model.plot_components(forecast);

子图的数量将取决于预测中已识别的组件数量。例如，如果包括了假期，那么将会显示holiday组件。在我们的示例中，将有两个子图：trend和yearly：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file247.png

图 11.6 – 绘制显示趋势和季节性（年度）的组成部分

图 11.6 展示了训练数据的趋势和季节性。如果你查看图 11.6，你会看到一个持续上升的趋势，自 1972 年以来逐渐稳定（放缓）。此外，季节性模式显示出夏季时生产量的增加。

趋势图中的阴影区域表示估算趋势的不确定性区间。数据存储在 forecast 数据框的 trend_lower 和 trend_upper 列中。

  1. 最后，与样本外数据（测试数据）进行比较，看看模型的表现如何：
ax = test.plot(x='ds', y='y',
                    label='Actual',
                    style='-.',
                    figsize=(12,4))
predicted.plot(x='ds', y='yhat',
               label='Predicted',
               ax=ax,
               title='Milk Production Actual vs Forecast');

将以下图表与第十章的图 10.43进行比较，看看Prophet与使用 auto_arima 获得的SARIMA模型有何不同：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file248.png

图 11.7 – 将 Prophet 的预测与测试数据进行比较

注意，对于高度季节性的牛奶生产数据，模型表现非常出色。通常，Prophet 在处理强季节性时间序列数据时表现优异。

工作原理…

Prophet 精简了构建和优化时间序列模型的多个方面，但在开始时需要一些关键指令，以便 Prophet 正确调整模型。例如，在初始化模型时，你需要决定季节性效应是加法模型还是乘法模型。你还需要指定诸如数据频率（例如，对于按月数据的freq='MS'）等参数，并使用 make_future_dataframe 方法来扩展预测范围。

当你实例化模型 model = Prophet() 时，Prophet 会使用默认的参数值，例如 yearly_seasonality='auto'、weekly_seasonality='auto' 和 daily_seasonality='auto'。这使得 Prophet 可以自动根据数据来确定需要包含的季节性分量。在牛奶生产数据集中，只检测到年度季节性，如下所示：

model.seasonalities
>>
OrderedDict([('yearly',
              {'period': 365.25,
               'fourier_order': 10,
               'prior_scale': 10.0,
               'mode': 'additive',
               'condition_name': None})])

Prophet 还提供了不确定性区间来预测结果，这些区间受到三个因素的影响：

  • 观测噪声：指的是无法通过模型解释的观测数据中的随机变化。

  • 参数不确定性：指的是在估计模型参数时的的不确定性。例如，调整 mcmc_samples 参数（马尔可夫链蒙特卡罗 或 MCMC 采样）来获取季节性分量的不确定性。默认值为零 (0)。

  • 未来趋势不确定性：指的是基于历史数据对未来趋势变化的预期不确定性。例如，增加 changepoint_prior_scale 参数值可以增加预测的不确定性。默认值设置为 0.05。此外，区间宽度也可以通过 interval_width 参数进行调整（默认为 0.80 或 80%）。

默认情况下，uncertainty_samples 参数设置为 1000，这意味着 Prophet 会使用**哈密顿蒙特卡罗（HMC）**算法进行 1000 次模拟来估计不确定性。你可以调整此值以控制模拟次数，或者通过设置 uncertainty_samples=0 或 uncertainty_samples=False 来完全关闭不确定性估计。如果禁用不确定性样本，Prophet 会从预测结果中省略像 yhat_lower 和 yhat_upper 这样的不确定性区间。

model = Prophet(uncertainty_samples=False).fit(train)
forecast = model.predict(future)
forecast.columns.tolist()
>>
['ds', 'trend', 'additive_terms', 'yearly', 'multiplicative_terms', 'yhat']

Prophet 的优势在于它能够自动检测变化点，这些是趋势发生显著变化的时间点。默认情况下，Prophet 会在训练数据的前 80% 中识别 25 个潜在变化点。你可以通过调整 n_changepoints 参数来修改此行为，或者通过 changepoint_range 控制用于变化点检测的历史数据量，默认值为 0.8（即 80%）。

你可以通过模型的 changepoints 属性检查检测到的变化点。例如，以下代码显示前五个变化点：

model.changepoints.shape
>>
(25,)
model.changepoints.head()
>>
5    1962-06-01
10   1962-11-01
14   1963-03-01
19   1963-08-01
24   1964-01-01
Name: ds, dtype: datetime64[ns]

这些变化点也可以在图表上进行可视化。以下代码将变化点叠加到原始时间序列数据上：

ax = milk.set_index('ds').plot(figsize=(12,5))
milk.set_index('ds').loc[model.changepoints].plot(style='X', ax=ax)
plt.legend(['original data', 'changepoints']);

这应该会生成一个图表，展示原始时间序列及 25 个潜在变化点，显示 Prophet 识别的趋势变化时刻。

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file249.png

图 11.8 – Prophet 所识别的 25 个潜在变化点

这些潜在变化点是从训练数据的前 80% 中估算出来的。

在接下来的章节中，你将更详细地探索变化点检测。

还有更多内容…

要绘制捕捉趋势中影响较大的变化点，你可以使用add_changepoints_to_plot函数，如以下代码所示：

from prophet.plot import add_changepoints_to_plot
fig = model.plot(forecast, ylabel='Milk Production in Pounds')
add_changepoints_to_plot(fig.gca(), model, forecast);

这将生成类似于图 11.8的图表，但会有额外的变化点线和趋势线。图 11.9 中的 25 个变化点中有 10 个（垂直线是显著的变化点）。线性趋势线应与图 11.6中显示的趋势成分相同：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file250.png

图 11.9 – 显示十个显著的变化点和趋势线

注意趋势线在识别出的变化点处的变化。这就是 Prophet 如何检测趋势变化的方式。由于应用了分段回归来构建趋势模型，趋势线并不是完全直线。当考虑分段线性模型时，你可以把它看作是多个线性回归线段（在显著变化点之间），然后将这些线段连接起来。这使得模型能够灵活地捕捉到趋势中的非线性变化，并进行未来预测。

Prophet 还包含交叉验证功能，以更好地评估模型在预测未来数据时的表现。交叉验证用于确保模型能够很好地泛化到未见过的数据，而不是出现过拟合或欠拟合的情况。此外，交叉验证还可以帮助你确定预测结果在多长时间内仍然可靠。

Prophet 提供了cross_validation和performance_metrics函数，允许你将数据拆分成训练集和测试集，并跨多个时段进行验证。这种方法通过在不同时间点进行预测并与实际值进行比较，帮助评估模型的准确性。

这是在 Prophet 中实现交叉验证的方法。

from prophet.diagnostics import cross_validation, performance_metrics
df_cv = cross_validation(model, initial='730 days', period='180 days', horizon='365 days')
df_cv.head()
>>
         ds        yhat  yhat_lower  yhat_upper    y     cutoff
0 1964-03-01  689.889300  685.612439  694.504671  688 1964-02-19
1 1964-04-01  701.435214  697.157285  706.019257  705 1964-02-19
2 1964-05-01  776.047139  771.707065  780.994528  770 1964-02-19
3 1964-06-01  735.045494  730.374821  739.547374  736 1964-02-19
4 1964-07-01  671.333097  666.625404  675.994830  678 1964-02-19

在前面的代码中，我们指定了以下参数：

initial：初始训练期的大小。在我们的代码中，我们指定了 730 天，大约是 2 年，或者是 24 个月，用于我们的月度牛奶生产数据。

period：进行预测时切分点之间的间隔。在这个例子中，我们指定了 180 天（大约 6 个月）。这意味着在初始训练后，Prophet 会按 6 个月的增量向前推进。

horizon：预测的时间跨度（预测未来的时间长度）。在这个例子中，我们指定了 365 天或 1 年。这意味着 Prophet 将在每个切分点之后的 12 个月（1 年）内进行预测。

在执行交叉验证之后，你可以使用performance_metrics函数评估模型的表现。

df_p = performance_metrics(df_cv)
print(df_p.iloc[: , 0:-1].head())
>>
horizon         mse       rmse        mae      mape     mdape     smape
0 41 days  226.788248  15.059490  12.300991  0.016356  0.016894  0.016345
1 42 days  220.336066  14.843721  11.849186  0.015699  0.015678  0.015694
2 45 days  214.385008  14.641892  11.647620  0.015503  0.015678  0.015503
3 46 days  207.646253  14.409936  11.380352  0.015170  0.014446  0.015164
4 47 days  242.132208  15.560598  12.179413  0.015953  0.014446  0.015986

该函数计算多个指标，如平均绝对误差（MAE）、均方根误差（RMSE）等。

你可以使用plot_cross_validation_metric函数来可视化交叉验证预测的表现：

from prophet.plot import plot_cross_validation_metric
fig = plot_cross_validation_metric(df_cv, metric='rmse');

这应该会绘制不同预测时段的 RMSE。

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file251.png

图 11.10 – 显示模型在不同预测期的 RMSE 曲线

对图形的解释表明，该模型在短期预测期内表现较好，因为 RMSE 相对较低。随着预测期的增加，RMSE 似乎普遍增加。通常，我们期望这些模型在短期预测（1-3 个月）内表现更好，误差更小。

另请参见

到目前为止，你一直在处理单变量时间序列。接下来的章节将教你如何处理多元时间序列。

使用 VAR 预测多元时间序列数据

在本章节中，你将探索用于处理多元时间序列的向量自回归（VAR）模型。在第十章，《使用统计方法构建单变量时间序列模型》中，我们讨论了 AR、MA、ARIMA 和 SARIMA 作为单变量单向模型的示例。而 VAR 则是双向且多元的。

VAR 与 AR 模型的比较

你可以将阶数为 p 的 VAR，或称VAR§，视为单变量 AR§ 的一种推广，用于处理多个时间序列。多个时间序列被表示为一个向量，因此命名为向量自回归。滞后为一（1）的 VAR 可以表示为 VAR(1)，涉及两个或更多变量。

还有其他形式的多元时间序列模型，包括向量移动平均（VMA）、向量自回归移动平均（VARMA）和向量自回归积分移动平均（VARIMA），它们是对其他单变量模型的推广。在实践中，你会发现由于其简单性，VAR 使用得最多。VAR 模型在经济学中非常流行，但你也会在其他领域看到它的应用，如社会科学、自然科学和工程学。

多变量时间序列的前提是，当利用多个时间序列（或输入变量）而不是单一时间序列（单一变量）时，可以增强预测的能力。简而言之，当你有两个或更多的时间序列，并且它们彼此间有（或假设有）相互影响时，VAR 就可以被使用。这些通常被称为内生变量，其关系是双向的。如果变量或时间序列之间没有直接关系，或者我们不知道它们是否存在直接的影响，则称其为外生变量。

外生与内生变量

当你开始深入研究 VAR 模型时，你会遇到关于内生和外生变量的引用。从宏观角度来看，这两者是相互对立的，在statsmodels中，它们分别被称为endog和exog。

内生变量受到系统内其他变量的影响。换句话说，我们预期一个状态的变化会影响其他状态。有时，这些变量在机器学习文献中被称为依赖变量。你可以使用Granger 因果检验来确定多个内生变量之间是否存在这种关系。例如，在statsmodels中，你可以使用grangercausalitytests。

另一方面，外生变量位于系统外部，并且不直接影响其他变量。它们是外部影响因素。有时，这些变量在机器学习文献中被称为独立变量。

类似于 AR 模型，VAR 模型假设时间序列变量的平稳性。这意味着每个内生变量（时间序列）需要是平稳的。

为了说明 VAR 是如何工作的以及背后的数学公式，我们从一个简单的 VAR(1)模型开始，该模型包含两个（2）内生变量，表示为(https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file252.png)。回顾一下第十章，使用统计方法构建单变量时间序列模型，AR(1)模型的形式如下：https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file253.png

一般来说，AR§模型是自我过去值的线性模型，其中§参数告诉我们应该追溯多远。现在，假设你有两个AR(1)模型，分别对应两个不同的时间序列数据。其形式如下：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file254.jpg

然而，这两个模型是独立的，没有展示出任何相互关系或相互影响。如果我们创建这两个模型的线性组合（即自身过去值和另一个时间序列的过去值），我们将得到以下公式：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file255.jpg

图 11.11 – 带滞后一期的 VAR 模型公式，或 VAR(1)

上述方程看起来可能很复杂，但最终，像 AR 模型一样，它仍然只是过去滞后值的线性函数。换句话说，在 VAR(1) 模型中，每个变量都将有一个滞后(1)的线性函数。当拟合 VAR 模型时，正如你将在本食谱中看到的那样，使用最小二乘法(OLS)方法来估计每个方程的 VAR 模型。

准备开始

在本食谱中，你将使用 pandas_datareader 库从 FRED（联邦储备经济数据）下载数据。相关文件也可以从本书的 GitHub 仓库下载。

使用conda安装：

conda install -c anaconda pandas-datareader

使用pip安装：

pip install pandas-datareader

如何操作……

在本食谱中，你将使用 pandas_datareader 库中的 FredReader() 函数来提取两个不同的时间序列数据集。正如 FRED 网站上所提到的，第一个符号 FEDFUNDS 是联邦基金有效利率，它*“是存款机构之间互相过夜交易联邦基金（存放在联邦储备银行的余额）的利率”。简单来说，联邦基金有效利率影响借贷成本。它是由联邦公开市场委员会（FOMC）设定的目标利率，用于确定银行向其他机构收取的利率，特别是借出其储备余额中的多余现金的利率。* 第二个符号是 unrate，代表失业率，它是指总劳动人口中正在积极寻找工作或愿意工作的失业人员的比例。

引用

美国联邦储备系统理事会，联邦基金有效利率 [FEDFUNDS]，来自 FRED，圣路易斯联邦储备银行； https://fred.stlouisfed.org/series/FEDFUNDS，2024 年 10 月 6 日。

美国劳动统计局，失业率 [UNRATE]，来自 FRED，圣路易斯联邦储备银行； https://fred.stlouisfed.org/series/UNRATE，2024 年 10 月 6 日。

按照以下步骤进行：

  1. 开始时加载必要的库并提取数据。注意，FEDFUNDS 和 unrate 都是按月报告的：
import pandas_datareader.data as web
import pandas as pd
from statsmodels.tsa.api import VAR,adfuller, kpss
from statsmodels.tsa.stattools import grangercausalitytests
  1. 使用 FredReader 提取数据，FredReader 封装了 FRED API，并返回一个 pandas DataFrame。对于 FEDFUNDS 和 unrate 符号，你将提取大约 34 年的数据（417 个月）：
import pandas_datareader.data as web
start = "01-01-1990"
end = "01-09-2024"
economic_df = web.FredReader(
        symbols=["FEDFUNDS", "unrate"],
        start=start,
        end=end).read()
file = '../../datasets/Ch11/economic_df.pickle'
economic_df.to_pickle(file)

将数据框存储为 pickle 对象，如前面的代码最后一行所示。这样，你就不需要再次调用 API 来重新运行示例。你可以使用 economic_df = pd.read_pickle(file) 读取 economic_df.pickle 文件。

  1. 检查数据，确保没有空值：
economic_df.isna().sum()
>>
FEDFUNDS    0
unrate      0
dtype: int64
  1. 将数据框的频率更改为月初（MS），以反映数据存储的方式：
economic_df.index.freq = 'MS'
  1. 绘制数据集进行可视化检查和理解：
economic_df.plot(subplots=True); plt.show()

由于 subplots 设置为 True，这将为每一列生成两个子图：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file256.png

图 11.12 – 绘制联邦基金有效利率和失业率

FEDFUND和unrate之间存在某种反向关系——随着FEDFUNDS的增加，unrate减少。从 2020 年开始，由于 COVID-19 疫情，出现了有趣的异常行为。我们可以进一步检查这两个变量之间的相关性：

correlation_matrix = economic_df.corr()
correlation_matrix
>>
          FEDFUNDS    unrate
FEDFUNDS  1.000000 -0.435171
unrate   -0.435171  1.000000

FEDFUNDS和unrate之间的相关系数为-0.435，表明存在中等强度的负相关关系（反向相关）。这表明，随着联邦基金利率的上升，失业率趋向下降，反之亦然。

进一步，你可以执行互相关函数（CCF）分析FEDFUNDS和unrate之间不同滞后的相关性。CCF 有助于识别两个时间序列之间的滞后关系或时间依赖性。结果主要告诉我们一个系列是否领先或滞后于另一个系列。这并不是一个正式的因果性检验，你将在后续步骤中进行更深入的调查。

from statsmodels.graphics.tsaplots import plot_ccf
import numpy as np
lags = np.arange(-12, 13)
plot_ccf(economic_df['FEDFUNDS'], economic_df['unrate'], lags=lags)
plt.grid(True)

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file257.png

图 11.13 – 前后 12 个月的互相关函数图，便于更好的解释

图中显示了超出阴影置信区间的尖峰，这表明存在显著的负相关。这暗示着FEDFUNDS的变化可能与unrate的相反变化有关，且这些变化发生在同一时段内。

  1. VAR 模型中的一个重要假设是平稳性。两个变量（这两个内生时间序列）需要是平稳的。创建一个check_stationarity()函数，该函数返回扩展型迪基-富勒（adfuller）检验的平稳性结果：
from statsmodels.tsa.stattools import adfuller
def check_stationarity(df):
    adf_pv = adfuller(df)[1]
    result = 'Stationary' if adf_pv < 0.05 else "Non-Stationary"
    return result

使用check_stationarity函数评估每个内生变量（列）：

for i in economic_df:
    adf = check_stationarity(economic_df[i])
    print(f'{i} adf: {adf}')
>>
FEDFUNDS adf: Stationary
unrate adf: Stationary

总体而言，这两个时间序列被证明是平稳的。

  1. 绘制ACF和PACF图，以便对每个变量及其所属的过程（例如 AR 或 MA 过程）获得直观理解：
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
for col in economic_df.columns:
    fig, ax = plt.subplots(1,2, figsize=(18,4))
    plot_acf(economic_df[col], zero=False,
             lags=30, ax=ax[0], title=f'ACF - {col}')
    plot_pacf(economic_df[col], zero=False,
              lags=30, ax=ax[1], title=f'PACF - {col}');

这应该生成FEDFUNDS和unrate的 ACF 和 PACF：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file258.png

图 11.11 – FEDFUNDS和unrate的 ACF 和 PACF 图

FEDFUNDS和unrate的 ACF 和 PACF 图表明我们正在处理一个自回归（AR）过程。ACF 图显示出逐渐缓慢衰减，而 PACF 图则显示出在滞后 1 之后的急剧截断。FEDFUNDS的 PACF 图在滞后 2 和滞后 3 时显示出略微显著的（超出阴影区域）滞后，而unrate的 PACF 图在滞后 24 时显示出一定的显著性。我们可以暂时忽略这些。

  1. 将数据分为训练集和测试集：
train = economic_df.loc[:'2022']
test = economic_df.loc['2023':]
print(f'Train: {len(train)}, Test: {len(test)}')
>>
Train: 396, Test: 21

  1. 你可以对数据进行缩放（标准化），尽管 VAR 模型并不要求标准化，且这两个时间序列的规模差异不大。在这一步中，你将进行缩放操作，仅用于演示目的，以展示如何进行此操作，并如何在后续步骤中反向转换结果以便进行解释。

    在实现 VAR 时，是否需要对变量进行标准化（缩放）常常存在争议。然而，VAR 算法本身并不要求变量必须进行标准化，因为它是与规模无关的。在实践中，通常保持变量的原始单位，以保留系数和残差的可解释性，这样可以用有意义的术语（例如百分比点）来解释。然而，如果变量在规模上差异显著，则可以应用标准化，以便更容易比较系数并更有效地检测异常值。

    如果你选择标准化（例如，使用Scikit-Learn中的StandardScalar），需要注意，结果将以标准差为单位。你可以使用inverse_transform方法将数据还原到其原始单位，以便在解释结果时使用。这对于需要在原始变量规模的上下文中解释结果时非常有用。

  2. 使用StandardScalar对数据进行缩放，通过fit方法拟合训练集。然后，使用transform方法对训练集和测试集应用缩放变换。变换将返回一个 NumPy ndarray，然后你可以将其作为 pandas DataFrame 返回。这样，你可以更方便地绘制图形并进一步检查 DataFrame：

from sklearn.preprocessing import StandardScaler
scale = StandardScaler()
scale.fit(train)
train_sc = pd.DataFrame(scale.transform(train),
                        index=train.index,
                        columns=train.columns)
test_sc = pd.DataFrame(scale.transform(test),
                       index=test.index,
                       columns=test.columns)
  1. 那么，如何确定 VAR 模型的最优滞后期数（p）呢？幸运的是，statsmodels中的 VAR 实现会自动选择最佳的 VAR 滞后期数。你只需要定义滞后的最大数量（阈值）；模型将确定最小化每个四种信息准则得分的最佳p值：AIC、BIC、FPE和HQIC。select_order方法会计算每个滞后期数的得分，而summary方法会显示每个滞后的得分。这些结果将有助于在训练（fit）模型时指定算法应使用哪些信息准则：
model = VAR(endog=train_sc)
res = model.select_order(maxlags=10)
res.summary()

这应该显示所有10个滞后的结果。最低的得分会用*标记：

VAR Order Selection (* highlights the minimums) 
==================================================
       AIC         BIC         FPE         HQIC  
--------------------------------------------------
0      -0.2862     -0.2657      0.7511     -0.2780
1       -7.345      -7.283   0.0006461      -7.320
2       -7.874      -7.772   0.0003805      -7.833
3       -7.960     -7.817*   0.0003491     -7.903*
4       -7.960      -7.775   0.0003492      -7.887
5       -7.951      -7.726   0.0003523      -7.862
6       -7.967      -7.701   0.0003467      -7.861
7      -7.974*      -7.667  0.0003443*      -7.852
8       -7.957      -7.609   0.0003502      -7.819
9       -7.947      -7.557   0.0003539      -7.792
10      -7.931      -7.501   0.0003593      -7.761
--------------------------------------------------

res对象是LagOrderResults类的一个实例。你可以使用selected_orders属性打印每个信息准则的选定滞后期数，该属性返回最优滞后期数的字典，包括AIC、BIC、FPE和HQIC：

print(res.selected_orders)
>>
{'aic': 7, 'bic': 3, 'hqic': 3, 'fpe': 7}

在滞后期 7 时，AIC和FPE得分最低。另一方面，BIC和HQ得分在滞后期 3 时最低，表明该模型可能更简洁，滞后期数较少。

通常，BIC和HQIC倾向于偏好更简洁的模型（即滞后期数较少的模型），因为它们对模型复杂度施加更高的惩罚。相反，AIC和FPE则较为宽松，可能建议使用更高的滞后期数，因为它们对复杂度的惩罚较少。

在这种情况下，滞后 3似乎是更安全的选择，旨在追求简洁性并避免过拟合。然而，选择滞后 7将会得到一个更复杂的模型，可能捕捉到数据中的更多动态，但也带有更高的过拟合风险。

  1. 要训练模型，必须使用BIC分数（或其他你选择的标准）。你可以通过更新ic参数来尝试不同的信息准则。maxlags参数是可选的——如果你留空，模型会根据所选信息准则自动确定最优滞后阶数，在此案例中为'bic'。
results = model.fit(ic='bic')

这将自动选择最小化 BIC 分数的滞后阶数（即滞后 3）。你可以通过设置ic='aic'来尝试其他标准，如 AIC。如果你更倾向于手动指定最大滞后阶数，可以使用maxlags参数。请记住，如果同时使用了maxlags和ic，则ic将优先考虑。

  1. 运行results.summary()将打印详细的总结，包括每个自回归过程的回归结果。

总结中包括了方程数量（对应变量的数量）、信息准则（AIC、BIC、HQIC 和 FPE）及其他细节，如对数似然。这些指标有助于评估模型的整体拟合度和复杂度。

 Summary of Regression Results  
==================================
Model:                         VAR
Method:                        OLS
Date:           Mon, 07, Oct, 2024
Time:                     14:49:51
--------------------------------------------------------------------
No. of Equations:         2.00000    BIC:                   -7.83732
Nobs:                     393.000    HQIC:                  -7.92278
Log likelihood:           466.565    FPE:                0.000342624
AIC:                     -7.97888    Det(Omega_mle):     0.000330737
--------------------------------------------------------------------

在总结的最后，有残差相关矩阵。该矩阵显示了每个方程中残差（误差）之间的相关性。理想情况下，你希望这些非对角线值尽可能接近零。较低的相关性表明模型已经捕捉到了大部分变量之间的关系，且残差中剩余的相关性较小。

Correlation matrix of residuals
            FEDFUNDS    unrate
FEDFUNDS    1.000000 -0.115904
unrate     -0.115904  1.000000

FEDFUNDS和unrate之间的残差相关性为-0.1159，相对较小，表明模型很好地捕捉到了这两个变量之间的关系。

  1. 使用k_ar属性存储 VAR 滞后阶数，以便在以后使用预测方法时能再次使用：
lag_order = results.k_ar
lag_order
>> 3

这显示选择的最优滞后阶数为 3。

最后，你可以使用forecast或forecast_interval方法进行预测。后者将返回预测值，以及上和下置信区间。这两种方法都需要过去的值和预测的步数。先前的值（past_y）将作为预测的初始值：

past_y = train_sc[-lag_order:].values
n = test_sc.shape[0]
forecast, lower, upper = results.forecast_interval(past_y, steps=n)
  1. 由于你对数据集应用了StandardScalar，预测值已被缩放。你需要应用inverse_transform将它们转换回原始数据的尺度。你还可以将预测和置信区间转换为 pandas DataFrame，方便使用：
forecast_df = pd.DataFrame(scale.inverse_transform(forecast),
                           index=test_sc.index,
                           columns=test_sc.columns)
lower_df = pd.DataFrame(scale.inverse_transform(lower),
                        index=test_sc.index,
                        columns=test_sc.columns)
upper_df = pd.DataFrame(scale.inverse_transform(upper),
                        index=test_sc.index,
                        columns=test_sc.columns)
  1. 绘制实际值与预测值unrate及其置信区间：
idx = test.index
plt.figure(figsize=(10, 6))
plt.plot(idx, test['unrate'], label='Actual unrate')
plt.plot(idx, forecast_df['unrate'], label='Forecasted unrate', linestyle='dashed')
plt.fill_between(idx, lower_df['unrate'], upper_df['unrate'], alpha=0.2, label='Confidence Interval')
plt.title('Actual vs Forecasted unrate with Confidence Intervals')
plt.legend()
plt.show()

这应该绘制实际的测试数据，unrate 的预测（中间虚线），以及上下置信区间：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file259.png

图 11.12 – 绘制带有置信区间的 unrate 预测

在 VAR 模型中，所有变量都同时进行预测，因为该模型假设系统中所有变量的过去值都会对每个变量的未来值产生影响。因此，即使你只关心预测 unrate 来自 FEDFUNDS，该模型也会基于 unrate 预测 FEDFUNDS 的未来值。你可以像绘制 unrate 的预测一样，绘制 FEDFUNDS 的预测，如下所示：

idx = test.index
plt.figure(figsize=(10, 6))
plt.plot(idx, test['FEDFUNDS'], label='Actual FEDFUNDS')
plt.plot(idx, forecast_df['FEDFUNDS'], label='Forecasted FEDFUNDS', linestyle='dashed')
plt.fill_between(idx, lower_df['FEDFUNDS'], upper_df['FEDFUNDS'], alpha=0.2, label='Confidence Interval')
plt.title('Actual vs Forecasted FEDFUNDS with Confidence Intervals')
plt.legend()
plt.show()

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file260.png

图 11.13 – 绘制带有置信区间的 FEDFUNDS 预测

请记住，训练数据是设置到 2022 年底的，包括 2020 年 COVID-19 导致的重大经济冲击。因此，由于数据的巨大且突然而来的变化，模型的预测可能无法如预期般表现。

你可能需要调整训练-测试划分来适应这一事实，并对模型进行不同的实验。当处理像 COVID-19 这样的重要异常时，必须评估事件是一次性异常，其影响将逐渐消退，还是它代表了一种持久的变化，需要在模型中进行特殊处理。

在这里，领域知识变得至关重要：理解经济背景将帮助你决定是否将此类事件建模，或将其视为不应影响未来预测的异常值。

它是如何工作的…

向量自回归模型（VAR）在计量经济学中非常有用。在评估 VAR 模型的性能时，通常会报告几个重要的分析结果，例如Granger 因果关系检验、残差（或误差）分析和冲击响应分析。这些对理解变量之间的相互作用至关重要。你将在下一个章节中探索这些内容，评估向量自回归（VAR）模型。

在图 11.11中，两个变量的 VAR(1) 模型产生了两个方程。每个方程将一个变量的当前值建模为其自身滞后值和另一个变量滞后值的函数。矩阵符号显示了四个系数和两个常数。每个方程有两个系数，分别用于两个变量的滞后值（例如，对于和对于）。常数和代表每个方程的截距。

在本配方中，选择的模型是针对两个变量的 VAR(3)模型。VAR(1)和 VAR(3)之间的主要区别在于由于附加滞后项，复杂性有所增加。在 VAR(3)模型中，你将需要求解十二个系数（每个变量三个滞后期）和两个常数。每个方程都是使用OLS估计的，如results.summary()输出所示。

图 11.11中表示的两个函数显示了 VAR(1)中每个变量不仅受到自身过去值的影响，还受到第二个（内生）变量的过去值的影响。这是 VAR 模型和 AR 模型（仅考虑变量自身滞后）的一个关键区别。VAR 模型捕捉了多个变量之间随时间变化的动态交互。

在下一个配方中，既然模型已经拟合完成，你将花时间绘制 VAR 特定的输出——例如，脉冲响应（IRs）和预测误差方差分解（FEVD）——以更好地理解这些相互作用以及变量之间的影响。

在本配方中，我们专注于使用模型进行预测；接下来，我们将专注于理解因果关系，使用 Granger 因果关系检验等工具，并通过额外的诊断评估模型的整体性能。

VARIMA（向量 ARIMA）模型扩展了 VAR，用于通过加入差分来处理非平稳数据。我们在这里没有使用它，因为两个时间序列都是平稳的，所以不需要差分。在处理多个需要差分以实现平稳的非平稳变量时，可以考虑使用 VARIMA。

还有更多……

由于我们关注的是比较预测结果，观察包含两个内生变量的 VAR(3)模型是否比单变量 AR(3)模型表现更好是很有趣的。将多元 VAR(3)模型与更简单的 AR(3)（或 ARIMA(3,0,0)）模型进行比较有助于评估第二个变量（FEDFUNDS）的加入是否改善了对unrate的预测。

你可以尝试为unrate时间序列拟合 AR(3)或 ARIMA(3,0,0)模型，并使用相同的滞后值以保持一致性。由于unrate序列是平稳的，因此无需进行差分。回想一下，从之前的活动中，lag_order等于 3。

from statsmodels.tsa.arima.model import ARIMA
model = ARIMA(train['unrate'],
              order=(lag_order,0,0)).fit()

你可以使用model.summary()查看 ARIMA 模型的摘要。拟合模型后，你可以通过绘制诊断图表来评估残差，以检查模型拟合是否存在问题：

fig = model.plot_diagnostics(figsize=(12,6));
fig.tight_layout()
plt.show()

这将生成四个诊断子图：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file261.png

图 11.14 – AR(3)或 ARIMA(3,0,0)模型诊断图输出

标准化残差图显示在 2020 年左右出现显著峰值，这可能是由于 COVID-19 的经济影响。根据 Q-Q 图，还可以观察到一些偏离正态分布的迹象，表明该模型可能未完全捕捉数据尾部行为。总体而言，AR 模型根据自相关图捕捉到了必要的信息。

现在，使用 AR(3)模型预测未来步骤，并将其与我们已经拟合并应用了inverse_transform的 VAR(3)模型进行比较：

# Forecast from AR(3) model
ar_forecast = pd.Series(model.forecast(n), index=test.index)
# Plot the forecast for AR(3) against the actual data
idx = test.index
plt.figure(figsize=(10, 4))
plt.plot(idx, test['unrate'], label='Actual unrate')
plt.plot(idx, ar_forecast, label='Forecasted unrate', linestyle='dashed')
plt.title('AR(3) model - Actual vs Forecasted unrate')
plt.legend()
plt.show()

这应该会生成以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file262.png

图 11.15 – AR(3)预测与实际结果对比（测试集）

同样的图表也可以为 VAR(3)模型绘制，该模型考虑了FEDFUNDS对unrate的影响：

# plotting VAR(3) same code as before but without confidence intervals
idx = test.index
plt.figure(figsize=(10, 4))
plt.plot(idx, test['unrate'], label='Actual unrate')
plt.plot(idx, forecast_df['unrate'], label='Forecasted unrate', linestyle='dashed')
plt.title('VAR(3) model - Actual vs Forecasted unrate')
plt.legend()
plt.show()

这应该会生成以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file263.png

图 11.16 – VAR(3)预测与实际结果对比（测试集）

最后，我们将计算均方根误差（RMSE）得分，以比较两个模型（VAR 和 AR）在测试数据上的表现：

from statsmodels.tools.eval_measures import rmse
rmse_var = rmse(test['FEDFUNDS'], forecast_df['unrate'])
print('VAR(3) RMSE = ', rmse_var)
rmse_ar = rmse(test['FEDFUNDS'], ar_forecast)
print('AR(3) RMSE = ', rmse_ar)
>>
VAR(3) RMSE =  0.9729655920788434
AR(3) RMSE =  0.7693416723850359

AR(3)模型相较于 VAR(3)模型具有更低的 RMSE，表明仅考虑unrate过去值的简单 AR(3)模型在此案例中表现更佳。包含unrate和FEDFUNDS的 VAR(3)模型并未显著提高unrate的预测精度。

根据结果，可能更倾向于选择 AR(3)模型，因为其复杂度较低且预测精度较高。然而，当变量之间存在更强的相互作用时，VAR 模型可以提供额外的洞察力和更精确的预测。

另见…

要了解有关 statsmodels 中 VAR 类的更多信息，请访问官方文档：www.statsmodels.org/dev/generated/statsmodels.tsa.vector_ar.var_model.VAR.html。

评估向量自回归（VAR）模型

在拟合 VAR 模型之后，下一步是评估模型如何捕捉不同内生变量（多个时间序列）之间的相互作用和动态关系。了解这些关系有助于您评估因果关系，变量如何相互影响，以及一个变量的冲击如何在系统中传播。

在本配方中，您将继续上一个配方《多元时间序列数据预测》中使用VAR的内容，并探索各种诊断工具以加深您对 VAR 模型的理解。具体来说，测试格兰杰因果关系，分析残差自相关函数（ACF）图，使用脉冲响应函数（IRF）以及进行预测误差方差分解（FEVD）。

这些评估步骤将帮助你理解系统中变量之间的因果关系和相互依赖关系，确保你的模型正确捕捉到潜在的动态。

如何操作…

以下步骤接着上一个步骤进行。如果你没有执行过这些步骤，可以运行附带的Jupyter Notebook中的代码来继续进行。你将专注于诊断之前创建的 VAR(3)模型：

  1. 首先，进行格兰杰因果关系检验，以确定一个时间序列是否可以用来预测另一个。在本例中，你需要找出FEDFUNDS是否可以用来预测unrate。

根据之前的步骤，你已经使用BIC准则选择了3 个滞后期。然而，你可以根据数据的特性或具体的研究问题，调整滞后期的顺序，并在需要时测试更高的滞后期。

  1. 格兰杰因果关系检验在statsmodels中通过grangercausalitytests()函数实现，该函数在每个过去的滞后期上执行四项检验。你可以通过maxlag参数控制滞后期的数量。格兰杰因果关系检验用于确定一个变量的过去值是否影响另一个变量。

格兰杰因果关系检验中的原假设是第二个变量（在本例中为FEDFUNDS）不对第一个变量（在本例中为unrate）产生格兰杰因果关系。换句话说，假设在影响或效应方面没有统计学意义。要测试FEDFUNDS是否影响unrate，你需要在应用检验之前交换列的顺序：

granger = grangercausalitytests(
            x=economic_df[['unrate', 'FEDFUNDS']],
            maxlag=3)

根据之前选择的 BIC 准则，测试的最大滞后期设为 3。输出将显示每个滞后期的检验统计量、p 值和自由度。重点关注p 值，以决定是否拒绝或接受原假设。

Granger Causality
number of lags (no zero) 1
ssr based F test:         F=0.5680  , p=0.4515  , df_denom=413, df_num=1
ssr based chi2 test:   chi2=0.5721  , p=0.4494  , df=1
likelihood ratio test: chi2=0.5717  , p=0.4496  , df=1
parameter F test:         F=0.5680  , p=0.4515  , df_denom=413, df_num=1
Granger Causality
number of lags (no zero) 2
ssr based F test:         F=21.9344 , p=0.0000  , df_denom=410, df_num=2
ssr based chi2 test:   chi2=44.4039 , p=0.0000  , df=2
likelihood ratio test: chi2=42.1852 , p=0.0000  , df=2
parameter F test:         F=21.9344 , p=0.0000  , df_denom=410, df_num=2
Granger Causality
number of lags (no zero) 3
ssr based F test:         F=21.6694 , p=0.0000  , df_denom=407, df_num=3
ssr based chi2 test:   chi2=66.1262 , p=0.0000  , df=3
likelihood ratio test: chi2=61.3477 , p=0.0000  , df=3
parameter F test:         F=21.6694 , p=0.0000  , df_denom=407, df_num=3

由于滞后期 2 和 3 的 p 值都小于 0.05，这表明FEDFUNDS对unrate有格兰杰因果关系。换句话说，考虑到 2 个月或 3 个月的滞后期时，FEDFUNDS的过去值对unrate具有显著的预测能力。

  1. 接下来，你将探讨残差图。之前步骤中的results对象results = model.fit(ic='bic')是VARResultsWrapper类型，与VARResults类相同，拥有相同的方法和属性。首先使用plot_acorr方法绘制残差的ACF图：
results.plot_acorr(resid=True)
plt.show();

这将生成四个图（2x2 图——每个变量各两个）来展示残差的自相关和交叉相关。回忆图 11.11，对于两个变量，你将有两个函数，这将转化为 2x2 的残差子图。如果有三个变量，则会有一个 3x3 的子图：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file264.png

图 11.17 – 残差自相关和交叉相关图

不幸的是，图 11.17中的图表没有正确的标签。第一行的图表对应数据框架中的第一个变量（FEDFUNDS），第二行对应第二个变量（unrate）。

你需要观察残差中没有显著的自相关性，这在图 11.17中是成立的。这将确认 VAR(3)模型已经有效捕捉了变量之间的动态关系（FEDFUNDS 和 unrate），并且残差中没有模型未能捕捉和考虑的剩余模式。

  1. 如果你希望单独查看每个变量的残差，你可以通过提取resid属性来实现。这将返回一个每个变量残差的 DataFrame。你可以使用标准的plot_acf函数来创建自相关函数（ACF）图：
for col in results.resid.columns:
    fig, ax = plt.subplots(1,1, figsize=(10,2))
    plot_acf(results.resid[col], zero=False,
             lags=10, ax=ax, title=f'ACF - {col}')

这将生成两个 ACF 图——一个用于FEDFUNDS，另一个用于unrate。这将确认相同的发现——即残差中没有显著的自相关性。

  1. 接下来，你将进行**脉冲响应函数（IRF）**分析。IRF 帮助你可视化并理解一个变量的冲击如何随时间影响另一个变量（或其自身），这对于评估 VAR 模型中变量之间的动态相互作用至关重要。

你将使用irf方法分析系统对冲击的响应：

irf_output = results.irf()
irf_output.plot()
plt.show()

irf_output对象是IRAnalysis类型，可以访问plot()函数：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file265.png

图 11.18 – 显示一个变量的单位变化对另一个变量的冲击响应

脉冲响应函数（IRF）分析计算动态脉冲响应和近似标准误差，并显示一个变量的冲击（或脉冲）对其他变量随时间（滞后）变化的影响。在VAR 模型中，所有变量相互影响，IRF 追踪一个变量变化对其他变量随时间的影响。

例如，FEDFUNDS → unrate（位于图 11.18的左下角）显示了

unrate如何对FEDFUNDS增加一个单位在 10 个滞后期中的反应。立刻可以观察到unrate有明显的急剧负响应，这表明联邦基金利率的提高会降低失业率。该效应持续存在，但随着时间的推移逐渐减弱，这与货币政策如何影响经济中的就业一致。从滞后期 2 到 滞后期 3，我们可以看到有一个稳定期，响应水平有所平稳，这就是滞后效应开始显现的地方；初始下降发生得很快，但总的调整需要更长时间，因为unrate仍然低于起始点。在滞后期 3 之后，我们看到unrate略微上升，这种逐步调整表明unrate开始缓慢恢复，但在若干期内仍未完全回到冲击前的水平。

另一方面，unrate → FEDFUNDS（图 11.18 的右上角）显示了一个相对较小且略微正向的响应。这表明失业率（unrate）的提高会导致FEDFUNDS略微上升，但该效应会随着时间的推移减弱。

如果你只想查看一个变量对另一个变量的响应（而不是所有子图），你可以指定 impulse 和 response 变量：

fig = irf_output.plot(impulse='FEDFUNDS', response='unrate', figsize=(5, 7))
fig.tight_layout()
plt.show()
  1. 绘制累积响应效应：
irf.plot_cum_effects()
plt.show()

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file266.png

图 11.19 – 累积脉冲响应函数图

累积响应图展示了冲击效果如何随着时间的推移逐渐累积。例如，FEDFUNDS → unrate 图（图 11.19 的左下角），你可以看到 FEDFUNDS 增加一个单位对unrate的累积效应会导致unrate持续下降。累积响应有助于你评估这些冲击的长期影响。

  1. 接下来，你将进入预测误差方差分解（FEVD），以量化系统中每个变量的预测误差方差有多少是由自身冲击和其他变量的冲击造成的。换句话说，你需要了解每个变量冲击的贡献。

你可以使用 fevd() 方法计算 FEVD：

fv = results.fevd()

你可以使用 fv.summay() 和 fv.plot() 方法探索 FEVD 结果。它们提供相似的信息：

fv.plot()
plt.show()

这将生成两个 FEVD 图（图 11.20），每个变量一个。

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file267.png

图 11.20 – FEDFUNDS 和 unrate 变量的 FEVD 图

x 轴表示从 0 到 9 的时期（滞后期），y 轴表示每个冲击对预测误差方差的贡献百分比（0 到 100%）。从视觉上看，每个冲击的贡献在每个时期的总条形中占据一定部分。

FEDFUNDS 图（顶部）显示，FEDFUNDS的整个预测误差方差由其自身的冲击解释。这意味着unrate的冲击对FEDFUNDS的预测误差方差没有影响。

另一方面，失业率图（底部）开始时显示，unrate的预测误差方差主要由其自身的冲击解释。然而，从滞后 4 期开始，FEDFUNDS的贡献开始增大。到滞后 9 期时，大约 30%的unrate方差可以归因于FEDFUNDS的冲击。这表明，在较长的时间范围内，FEDFUNDS成为unrate预测误差方差的更重要驱动因素。

它是如何工作的……

来自 statsmodels 的 VAR 实现提供了多个工具，帮助你理解不同时间序列（或内生变量）之间的动态关系。这包括格兰杰因果检验、脉冲响应函数（IRFs）和预测误差方差分解（FEVD）等方法。这些工具为你提供了不同的视角，帮助你理解系统中变量如何随时间互动，不论是通过因果关系还是一个变量的冲击如何影响其他变量。

除了理解变量之间的关系外，还提供了诊断工具（如残差分析和自相关图），帮助你评估模型拟合度，并判断是否需要进行调整（如模型阶数、差分或增加变量）。

还有更多……

在之前的使用 VAR 预测多元时间序列数据示例中，你手动创建了预测图。然而，results对象（即VARResults类的实例）提供了一种便捷的方式，使用plot_forecast方法可以快速绘制你的预测图。此函数会自动生成预测值及其置信区间。

n = len(test)
results.plot_forecast(steps=n, plot_stderr=True);

在这里，n是未来的步数。这个操作应该会生成两个子图，每个变量都有一个预测图，如下图所示：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file268.png

图 11.21 – FEDFUNDS 和失业率的预测图

另见…

要了解更多关于VARResults类及其所有可用方法和属性的信息，请访问官方文档：www.statsmodels.org/dev/generated/statsmodels.tsa.vector_ar.var_model.VARResults.html。

使用 GARCH 模型预测金融时间序列数据的波动性

在处理金融时间序列数据时，一个常见的任务是衡量波动性，它代表了未来回报的不确定性。一般来说，波动性衡量的是在特定时期内回报的概率分布范围，通常通过方差或标准差（即方差的平方根）来计算。它被用作量化风险或不确定性的代理指标。换句话说，它衡量了金融资产回报围绕预期值的分散程度。更高的波动性意味着更高的风险。这有助于投资者理解他们能期望的回报水平，以及回报与预期值之间的差异频率。

我们之前讨论的大多数模型（如 ARIMA、SARIMA 和 Prophet）都专注于基于过去的值预测观察到的变量。然而，这些模型假设方差是恒定的（同方差性），并没有考虑方差随时间的变化（异方差性）。

在本教程中，你将处理另一种类型的预测：预测和建模方差随时间变化的情况。这就是所谓的波动性。一般而言，当存在不确定性时，波动性是衡量风险的一个重要指标，也是处理金融数据时的一个重要概念。

为了实现这一点，你将接触到一类与自回归条件异方差（ARCH）相关的新算法。ARCH 模型捕捉了方差随时间的变化，它将变化建模为来自过去时间点的平方误差项（创新）的函数。ARCH 的扩展是GARCH模型，代表广义自回归条件异方差。它通过增加一个移动平均成分来考虑过去的方差，从而提供了一个更灵活、更持久的波动性结构。

GARCH 在计量经济学中很受欢迎，金融机构广泛使用它来评估投资和市场风险。预测市场动荡和波动性与预测未来价格同样重要，许多交易策略——例如均值回归——都将波动性作为一个关键因素。

在继续之前，让我们分解一下一个典型 ARCH 模型的组成部分：

  • 自回归 - 这是我们在第十章，使用统计方法构建单变量时间序列模型中探讨过的一个概念，意味着变量的当前值受到其过去值的影响。在 ARCH 模型中，这意味着当前的波动性（方差）受过去的平方误差项（创新）的影响。

  • 异方差性 - 意味着模型在不同的时间点可能具有不同的幅度或波动性（方差随时间变化）。

  • 条件性 - 由于方差（波动性）不是固定的，它依赖于过去的信息。换句话说，在某一时点的方差条件性地依赖于该序列中过去的值，意味着随着新信息的到来，波动性会被更新。

在本例中，您将创建一个GARCH 模型，其阶数为（p, q），其中p表示滞后方差的数量（GARCH 项），q表示滞后平方残差的数量（ARCH 项）。

使用arch Python 库，您将通过arch_model函数实现这一点。该函数中的参数可以根据 GARCH 模型的假设分为三个主要组件：

  • dist ：控制创新项（残差）的分布假设，默认为 'normal'。

  • mean ：控制条件均值的模型，默认为'Constant'，即假设均值为常数。

  • vol ：控制波动率模型（条件方差），默认为'GARCH'。

在大多数文献中，q通常表示ARCH 阶数（滞后平方残差或创新项的数量），而p表示GARCH 阶数（滞后方差的数量）。然而，在arch Python 库中，p和q的角色是颠倒的。

在arch包中，p表示平方残差（ARCH组件）的滞后阶数，而q表示滞后方差（GARCH组件）的滞后阶数。

这个区别在使用arch库指定和解释 GARCH 模型时非常重要，因为它与学术文献中使用的常规符号有所不同。

例如，在arch_model中指定p=2和q=1，实际上会根据常规符号拟合一个**GARCH(1, 2)**模型（包含 1 个滞后方差和 2 个滞后残差）。

准备工作

在本例中，您将使用arch库，它包含了多个波动率模型以及金融计量经济学工具。该库的输出与 statsmodels 库的输出类似。撰写本文时，最新版本为7.1.0。

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

pip install arch

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

conda install -c conda-forge arch-py

操作步骤…

当使用arch库中的arch_model函数构建GARCH模型时，有几个关键参数：默认的分布是normal（dist='normal'），均值模型是常数（mean='Constant'），波动率模型是 GARCH（默认vol='GARCH'）。根据上下文，您还可以指定其他均值模型，如自回归模型（'AR'）。同样，您也可以选择其他波动率模型，如'GARCH'、'ARCH'、'EGARCH'、'FIGARCH'和'APARCH'。您还可以选择不同的分布，如'gaussian'、't'、'studentst'、'skewstudent'。

本例中，您将使用微软的日收盘价数据集。请按照以下步骤进行：

  1. 首先加载本例所需的库：
from arch import arch_model
import pandas as pd
  1. 加载MSFT.csv数据集：
msft = pd.read_csv('../../datasets/Ch11/MSFT.csv',
                   index_col='date',
                    usecols=['date', 'close'],
                   parse_dates=True)
  1. 你需要将每日股价转换为每日股价收益。这可以通过计算当前时间的价格与前一天的价格之间的比例来完成。在 pandas 中可以使用 DataFrame.pct_change() 函数，之后乘以 100 将收益表示为百分比。

pct_change() 函数接受一个 periods 参数，默认值为 1。如果你想计算 30 天的收益率，那么你需要将该值改为 pct_change(periods=30)：

msft['returns'] = 100 * msft.pct_change()
msft.dropna(inplace=True, how='any')

dropna 步骤是必要的，因为计算 returns 会对第一行产生 NaN。

  1. 你可以同时绘制每日股价和每日收益：
title = 'Microsoft Daily Closing Price and Daily Returns'
msft.plot(subplots=True,
          title=title);

这将产生以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file269.png

图 11.22 – 微软每日收盘价和日收益

  1. 将数据拆分为训练集和测试集。对于短期波动率预测，你将检查模型如何预测接下来几天的波动率。在拆分数据时，你将留出最后 5 天作为测试数据。
train = msft.returns[:-5] 
test = msft.returns[-5:]  
print(f'Train: {train.shape}')
print(f'Test: {test.shape}')
print(f'Train: {train.shape}')
print(f'Test: {test.shape}')
>>
Train: (1253,)
Test: (5,)
  1. 拟合一个GARCH(p, q) 模型。从简单的 GARCH(1,1) 模型开始，使用所有默认选项——即 mean='Constant'，分布为 dist='normal'，波动率为 vol='GARCH'，p=1，q=1。GARCH(1,1) 是最常用的模型，也是波动率预测的一个很好的起点，因其简洁性和有效性在捕捉波动率集群方面表现出色。
model = arch_model(train,
                   p=1, q=1,
                   mean='Constant',
                   vol='GARCH',
                   dist='normal')
results = model.fit(update_freq=5)
>>
Iteration:      3,   Func. Count:     22,   Neg. LLF: 2419.4197011866704
Iteration:      6,   Func. Count:     41,   Neg. LLF: 2409.599553434422
Iteration:      9,   Func. Count:     55,   Neg. LLF: 2409.592672855605
Optimization terminated successfully    (Exit mode 0)
            Current function value: 2409.59267285418
            Iterations: 9
            Function evaluations: 55
            Gradient evaluations: 9

拟合过程在九次（9）迭代中完成。使用 update_freq=3 会影响拟合过程中打印进度的频率。默认值是 1，意味着每次迭代都会打印输出。通过设置为三（3），我们每三次迭代打印一次输出。要打印总结，使用 summary() 方法：

Print(results.summary())

这将产生以下输出：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file270.png

图 11.22 – GARCH(1,1) 模型的总结

GARCH(1,1) 模型的omega、alpha 和 beta 参数（符号）是通过最大似然方法估计的。系数的 p-value 表明它们是统计学上 显著的。

你可以通过从 results 对象中调用相应的属性来访问总结表中看到的几个组件——例如，results.pvalues，results.tvalues，results.std_err，或 results.params。

print(results.params)
>>
mu          0.144878
omega       0.055152
alpha[1]    0.095416
beta[1]     0.891086
Name: params, dtype: float64
  1. 接下来，你需要评估模型的表现。首先，使用 plot 方法绘制标准化残差和条件波动率：
results.plot();

这应该产生以下图表：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file271.png

图 11.23 – GARCH(1,1) 模型的诊断

如果模型拟合良好，那么标准化残差图应该呈现出白噪声的样子，表明残差中没有残余模式。该图似乎表现出随机性，没有明显的模式，表明拟合良好。条件波动率图显示了模型估计的时间变化波动率，反映了波动率随时间变化的情况。你可以看到反映市场状况的高波动和低波动时期。

绘制标准化残差的直方图。你可以使用std_resid属性来获取该数据：

results.std_resid.hist(bins=20)
plt.title('Standardized Residuals')

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file272.png

图 11.24 – GARCH 模型标准化残差的直方图

直方图表明标准化残差大致符合正态分布，尽管你可以使用额外的统计检验，如Jarque-Bera检验，进行更正式的正态性评估。

  1. 你还可以通过使用Ljung-Box 检验来测试自相关性，该检验检验无自相关性的原假设。在前 10 个滞后期中，使用acorr_ljungbox。p 值大于 0.05 表示无法拒绝原假设，表明该滞后期没有显著的自相关性。
from statsmodels.stats.diagnostic import acorr_ljungbox
acorr_ljungbox(results.std_resid,
               lags=10,
               return_df=True)['lb_pvalue']
>>
1     0.050203
2     0.038038
3     0.077375
4     0.136003
5     0.195838
6     0.157237
7     0.201474
8     0.248204
9     0.153473
10    0.210838
Name: lb_pvalue, dtype: float64

来自 Ljung-Box 检验的 p 值表明，在大多数滞后期，标准化残差没有显著的自相关性。由于大多数 p 值大于 0.05，我们未能拒绝无自相关性的原假设，表明模型拟合良好。

  1. 要进行预测，请使用forecast()方法。默认情况下，它会生成一步（1 步）预测。要获取n步预测，你需要更新horizon参数：
msft_forecast = results.forecast(horizon=test.shape[0])
  1. 要访问预测的未来方差（或波动率），请使用msft_forecast对象中的variance属性：
forecast = msft_forecast.variance
print(forecast)
>>
                 h.1       h.2       h.3       h.4       h.5
date                                                       
2024-08-27  1.623692  1.656928  1.689714  1.722059  1.753967

你还可以评估预测的均值。回想一下，当你拟合模型时，你指定了均值为Constant。如果你检查均值，这一点将得到进一步验证：

print(msft_forecast.mean)
>>
                 h.1       h.2       h.3       h.4       h.5
date                                                       
2024-08-27  0.144878  0.144878  0.144878  0.144878  0.144878

这将在所有时间段中输出一个常数值0.144878。

工作原理……

GARCH(p,q)可以写成如下形式：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file273.jpgOmega、alpha 和 beta（https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file274.png）是此处的参数。p阶通常被称为 GARCH 阶，而q被称为 ARCH 阶。

GARCH 模型的参数表示如下：

  • Omega：常数或基准方差

  • Alpha：滞后平方残差的系数（ARCH 项），衡量近期冲击对波动率的影响。

  • Beta：滞后方差的系数（GARCH 项），表示波动率随时间的持续性。

在arch Python 包中，p和q的角色被交换，p表示 ARCH 成分（滞后平方残差），q表示 GARCH 成分（滞后方差）。

你实现的 GARCH(1,1)模型可以写成如下形式：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file275.png

时间序列中的创新与误差

在时间序列文献中，你会遇到创新这个术语，指的是无法通过过去的信息预测的意外和不可预测的新信息、冲击或误差。简而言之，你可以将创新视为每个时间步的预测误差或意外情况。虽然在机器学习中，我们通常称这些为预测误差，但在像 ARCH/GARCH 这样的时间序列模型中，我们使用创新这个术语来描述影响模型的新、出乎意料的信息。

还有更多…

之前，在实现 GARCH 模型时，你将均值设置为'Constant'。现在，让我们探讨将均值更改为'Zero'的影响，这实际上会将均值模型从方程中移除。

让我们从设置mean='Zero'开始：

model = arch_model(train,
                   p=1, q=1,
                   mean='Zero',
                   vol='GARCH',
                   dist='normal')
results = model.fit(disp=False)

这将生成一个表格格式的 GARCH 总结：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file276.png

图 11.25 – 带零均值的 GARCH(1, 1)

请注意，在图 11.25 中，没有像图 11.21 中那样的均值模型。使用零均值在你想将波动性建模与均值分开时很常见。当你主要关注建模波动性，而不需要用于基础回报的均值模型时，这种方法非常有帮助。

另见

要了解更多关于arch库的信息，请访问官方文档：arch.readthedocs.io/en/latest/univariate/introduction.html。

第十二章：14 无监督机器学习中的异常值检测

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

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file0.png

packt.link/zmkOY

在第八章，使用统计方法的异常值检测中，你探讨了使用参数和非参数统计技术来识别潜在的异常值。这些方法简单、可解释，而且相当有效。

异常值检测并非简单直观，主要是由于异常值的定义存在模糊性，这种定义取决于你的数据或你试图解决的问题。例如，尽管常见，但第八章中使用的某些阈值，使用统计方法的异常值检测，仍然是任意的，并非你必须遵循的规则。因此，拥有领域知识或能够访问主题专家（SMEs）对于正确判断异常值至关重要。

在本章中，你将会接触到基于机器学习的几种异常值检测方法。大多数机器学习异常值检测技术被视为无监督异常值检测方法，例如孤立森林（iForest）、无监督K-近邻算法（KNN）、局部异常因子（LOF）以及基于 Copula 的异常值检测（COPOD）等。

通常情况下，异常值（或异常情况）被认为是一种罕见现象（在本章稍后你会看到，这被称为污染率）。换句话说，你会假设在一个大型数据集中，只有一小部分数据是异常值。例如，1%的数据可能是潜在的异常值。然而，这种复杂性需要设计出能够发现数据模式的方法。无监督的异常值检测技术擅长于发现罕见现象中的模式。

在调查异常值之后，你将拥有一个历史标签数据集，允许你使用半监督的异常值检测技术。本章重点讲解无监督异常值检测。

在本章中，你将接触到PyOD库，它被描述为*“一个全面且可扩展的 Python 工具包，用于检测多变量数据中的异常对象。”*该库提供了一个广泛的实现集合，涵盖了异常值检测领域的流行算法和新兴算法，你可以在此阅读更多内容：github.com/yzhao062/pyod。

你将使用相同的纽约出租车数据集，这样可以更方便地比较本章不同机器学习方法和第八章，使用统计方法的异常值检测中的统计方法的结果。

本章中你将遇到的内容如下*：

  • 使用KNN检测异常值

  • 使用LOF检测异常值

  • 使用iForest检测异常值

  • 使用一类支持向量机（OCSVM）检测异常值

  • 使用COPOD检测异常值

  • 使用PyCaret检测异常值

技术要求

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

你可以使用 pip 或 Conda 安装 PyOD。对于 pip 安装，运行以下命令：

pip install pyod

对于 Conda 安装，运行以下命令：

conda install -c conda-forge pyod

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

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')
plt.rcParams["figure.figsize"] = [16, 3]

将 nyc_taxi.csv 数据加载到 pandas DataFrame 中，因为它将在整个章节中使用：

file = Path("../../datasets/Ch14/nyc_taxi.csv")
nyc_taxi_2 = pd.read_csv(file,
                     index_col='timestamp',
                     parse_dates=True)
nyc_taxi_2.index.freq = '30T'

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

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

创建你将在整个方法中使用的 plot_outliers 函数：

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

    if labels:
        for i in outliers['value'].items():
            plt.plot(i[0], i[1], 'v', markersize=8, markerfacecolor='none', markeredgecolor='k')
            plt.text(i[0], i[1]-(i[1]*0.04), f'{i[0].strftime("%m/%d")}',
                         horizontalalignment=halignment,
                         verticalalignment=valignment)
    else:
        data.loc[outliers.index].plot(ax=ax, style='rX', markersize=9)       
    plt.title(f'NYC Taxi - {method}')
    plt.xlabel('date'); plt.ylabel('# of passengers')
    plt.legend(['nyc taxi','outliers'])
    plt.show()

在进行异常值检测时，目标是观察不同技术如何捕捉异常值，并将其与真实标签进行比较，如下所示：

tx = nyc_taxi.resample('D').mean()
known_outliers = tx.loc[nyc_dates]
plot_outliers(known_outliers, tx, 'Known Outliers')

上述代码应生成一个时间序列图，其中标记 X 表示已知的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file278.jpg

图 14.1：对下采样后的纽约出租车数据进行绘图，并附上真实标签（异常值）

PYOD 的训练与预测方法

像 scikit-learn 一样，PyOD 提供了熟悉的方法来训练模型并进行预测，方法包括：model.fit()、model.predict() 和 model.fit_predict()。

在这些方法中，我们将过程分为两步，首先使用 .fit() 拟合模型（训练），然后使用 .predict() 进行预测。

除了 predict 方法，PyOD 还提供了两个附加方法：predict_proba 和 predict_confidence。

在第一个方法中，你将探索 PyOD 如何在幕后工作，并介绍基本概念，例如 contamination 的概念，以及如何使用 threshold_ 和 decision_scores_ 来生成二元标签（异常 或 正常）。这些概念将在接下来的方法中详细讨论。

使用 KNN 检测异常值

KNN 算法通常用于监督学习环境，其中先前的结果或输出（标签）是已知的。

它可以用来解决分类或回归问题。这个思想很简单；例如，你可以根据最近邻来分类一个新的数据点 Y。例如，如果 k=5，算法会通过与点 Y 的距离找到五个最近的邻居，并根据多数邻居的类别来确定 Y 的类别。如果五个邻居中有三个是蓝色，两个是红色，那么 Y 将被分类为蓝色。KNN 中的 K 是一个参数，你可以修改它以找到最佳值。

在离群点检测的情况下，算法的使用方式有所不同。由于我们事先不知道离群点（标签），KNN 以无监督的方式进行学习。在这种情况下，算法会为每个数据点找到最近的K个邻居，并计算它们的平均距离。与数据集其余部分的距离最远的点将被视为离群点，更具体地说，它们被认为是全局离群点。在这种情况下，距离成为确定哪些点是离群点的评分依据，因此 KNN 是一种基于邻近的算法。

一般来说，基于邻近的算法依赖于离群点与其最近邻之间的距离或接近度。在 KNN 算法中，最近邻的数量，k，是你需要确定的一个参数。PyOD 支持 KNN 算法的其他变种，例如，平均 KNN（AvgKNN），它使用与 KNN 的平均距离进行评分；以及中位数 KNN（MedKNN），它使用中位数距离进行评分。

如何实现…

在这个示例中，你将继续使用在技术要求部分创建的tx数据框，通过 PyOD 的KNN类来检测离群点：

  1. 首先加载KNN类：
from pyod.models.knn import KNN
  1. 你需要熟悉一些参数来控制算法的行为。第一个参数是contamination，一个数字（浮动）值，表示数据集中离群点的比例。这是 PyOD 中所有不同类别（算法）通用的参数。例如，contamination值为0.1表示你期望数据中 10%是离群点。默认值为contamination=0.1。contamination的值可以在0到0.5（即 50%）之间变化。你需要通过实验来调整contamination值，因为该值会影响用于确定潜在离群点的评分阈值，以及返回多少潜在离群点。你将在本章的*它是如何工作的…*部分深入了解这一点。

例如，如果你怀疑数据中的离群点比例为 3%，你可以将其作为contamination值。你可以尝试不同的contamination值，检查结果，并确定如何调整contamination水平。我们已经知道在 215 个观测值中有 5 个已知的离群点（约 2.3%），在这个示例中，你将使用 0.03（或 3%）。

第二个参数，特定于 KNN 的是method，其默认值为method='largest'。在本教程中，你将把它更改为mean（所有k邻居距离的平均值）。第三个参数，亦为 KNN 特有的是metric，它告诉算法如何计算距离。默认值是minkowski距离，但它可以接受来自 scikit-learn 或 SciPy 库的任何距离度量。最后，你需要提供邻居的数量，默认值为n_neighbors=5。理想情况下，你将希望使用不同的 KNN 模型并比较不同k值的结果，以确定最佳邻居数量。

  1. 使用更新后的参数实例化 KNN，然后训练（拟合）模型：
knn = KNN(contamination=0.03,
          method='mean',
          n_neighbors=5)
knn.fit(tx)
>>
KNN(algorithm='auto', contamination=0.05, leaf_size=30, method='mean',
  metric='minkowski', metric_params=None, n_jobs=1, n_neighbors=5, p=2,
  radius=1.0)
  1. predict方法将为每个数据点生成二进制标签，1或0。值为1表示是异常值。将结果存储在 pandas Series 中：
predicted = pd.Series(knn.predict(tx),
                      index=tx.index)
print('Number of outliers = ', predicted.sum())
>>
Number of outliers =  6
  1. 过滤predicted Series，仅显示异常值：
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>> 
Timestamp  value
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

总体来看，结果很有前景；已识别出五个已知日期中的四个。此外，算法还识别了圣诞节后的那一天，以及 2015 年 1 月 26 日，那一天由于北美暴风雪，所有车辆都被命令驶离街道。

  1. 使用在技术要求部分创建的plot_outliers函数来可视化输出，以获得更好的洞察：
plot_outliers(outliers, tx, 'KNN')

上述代码应生成类似于图 14.1的图表，不同之处在于x标记是基于使用 KNN 算法识别的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file279.jpg

图 14.2：使用 KNN 算法识别的潜在异常值标记

要打印标签（日期）和标记，只需再次调用plot_outliers函数，但这次要设置labels=True：

plot_outliers(outliers, tx, 'KNN', labels=True)

上述代码应生成一个类似于图 14.2的图表，并附加文本标签。

工作原理…

KNN 算法的无监督方法计算一个观测值与其他邻近观测值的距离。PyOD 中 KNN 的默认距离是闵可夫斯基距离（p-范数距离）。你可以更改为不同的距离度量，例如使用euclidean或l2表示欧几里得距离，或使用manhattan或l1表示曼哈顿距离。可以使用metric参数来实现这一点，metric参数可以接受字符串值，例如metric='l2'或metric='euclidean'，也可以是来自 scikit-learn 或 SciPy 的可调用函数。这是一个需要实验的参数，因为它影响距离的计算方式，而异常值分数正是基于此进行计算的。

传统上，当人们听到 KNN 时，他们立即认为它仅仅是一个监督学习算法。对于无监督 KNN，有三种常见的算法：ball tree、KD tree 和暴力搜索。PyOD 库支持这三种算法，分别为 ball_tree、kd_tree 和 brute。默认值设置为 algorithm="auto"。

PyOD 使用特定于每个算法的内部评分，对训练集中的每个观测值进行评分。decision_scores_ 属性将显示每个观测值的这些评分。较高的评分表示该观测值更有可能是异常值：

knn_scores = knn.decision_scores_

你可以将其转换为 DataFrame：

knn_scores_df = (pd.DataFrame(scores,
             index=tx.index,
             columns=['score']))
knn_scores_df

由于所有数据点都被评分，PyOD 会确定一个阈值来限制返回的异常值数量。阈值的大小取决于你之前提供的 污染 值（你怀疑的异常值比例）。污染值越高，阈值越低，因此返回的异常值越多。污染值越低，阈值会提高。

你可以使用模型在拟合训练数据后，从 threshold_ 属性中获取阈值。以下是基于 3% 污染率的 KNN 阈值：

knn.threshold_
>> 225.0179166666657

这是用来过滤显著异常值的值。以下是如何重现这一点的示例：

knn_scores_df[knn_scores_df['score'] >= knn.threshold_].sort_values('score', ascending=False)

输出结果如下：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file280.jpg

图 14.3：显示来自 PyOD 的决策评分

请注意，最后一个观测值 2014-09-27 稍微高于阈值，但在你使用 predict 方法时并没有返回。如果使用污染阈值，你可以得到一个更好的截止点：

n = int(len(tx)*0.03)
knn_scores_df.nlargest(n, 'score')

另一个有用的方法是 predict_proba，它返回每个观测值是正常值和异常值的概率。PyOD 提供了两种方法来确定这些概率：linear 或 unify。这两种方法都会在计算概率之前对异常值评分进行缩放。例如，在 linear 方法的实现中，使用了 scikit-learn 的 MinMaxScaler 来缩放评分，然后再计算概率。unify 方法则使用 z-score（标准化）和 SciPy 库中的高斯误差函数（erf）（scipy.special.erf）。

你可以比较这两种方法。首先，使用 linear 方法计算预测概率，你可以使用以下代码：

knn_proba = knn.predict_proba(tx, method='linear')
knn_proba_df = (pd.DataFrame(np.round(knn_proba * 100, 3),
            index=tx.index,
            columns=['Proba_Normal', 'Proba_Anomaly']))
knn_proba_df.nlargest(n, 'Proba_Anomaly')

对于 unify 方法，你只需将 method='unify' 更新即可。

要保存任何 PyOD 模型，你可以使用 joblib Python 库：

from joblib import dump, load
# save the knn model
dump(knn, 'knn_outliers.joblib')
# load the knn model
knn = load('knn_outliers.joblib')

还有更多内容…

在之前的步骤中，当实例化 KNN 类时，你将计算异常值 评分 的 method 值更改为 mean：

knn = KNN(contamination=0.03,
          method='mean',
          n_neighbors=5)

让我们为 KNN 算法创建一个函数，通过更新 method 参数为 mean、median 或 largest，以训练模型并检查这些方法对决策评分的影响：

  • largest 使用到 k 邻居的最大距离作为异常值评分。

  • mean 使用与 k 邻居的距离的平均值作为异常值分数。

  • median 使用与 k 邻居的距离的中位数作为异常值分数。

创建 knn_anomaly 函数，包含以下参数：data、method、contamination 和 k：

def knn_anomaly(df, method='mean', contamination=0.05, k=5):
    knn = KNN(contamination=contamination,
              method=method,
              n_neighbors=5)
    knn.fit(df)   
    decision_score = pd.DataFrame(knn.decision_scores_,
                          index=df.index, columns=['score'])
    n = int(len(df)*contamination)
    outliers = decision_score.nlargest(n, 'score')
    return outliers, knn.threshold_

你可以通过不同的方法、污染度和 k 值来运行该函数进行实验。

探索不同方法如何生成不同的阈值，这会影响异常值的检测：

for method in ['mean', 'median', 'largest']:
    o, t = knn_anomaly(tx, method=method)
    print(f'Method= {method}, Threshold= {t}')
    print(o)

前述代码应该会打印出每种方法的前 10 个异常值（污染度为 5%）：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file281.jpg

图 14.4：使用不同 KNN 距离度量比较决策分数

注意，前六个（表示 3% 污染度）对于三种方法是相同的。顺序可能会有所不同，各方法的决策分数也不同。请注意，各方法之间的差异在前六个之后更为明显，如 图 14.4 所示。

另见

查看以下资源：

使用 LOF 检测异常值

在前一个实例中，使用 KNN 检测异常值，KNN 算法通过观测点之间的距离来为检测异常值计算决策分数。与 KNN 距离较远的数据点可以被认为是异常值。总体来说，该算法在捕捉全局异常值方面表现不错，但对于那些远离周围点的数据点，可能无法很好地识别局部异常值。

这时，LOF（局部异常因子）就能解决这个问题。LOF 不使用邻近点之间的距离，而是通过密度作为基础来为数据点评分并检测异常值。LOF 被认为是一个基于密度的算法。LOF 的理念是，异常值会离其他数据点较远且更加孤立，因此会出现在低密度区域。

用一个例子来说明这个概念会更清楚：假设一个人站在小而忙碌的星巴克队伍中，每个人几乎都靠得很近；那么，我们可以说这个人处于高密度区域，更具体地说是高局部密度。如果这个人决定在停车场里等着，直到队伍稍微缓解，他就被孤立了，处于低密度区域，因此被认为是异常值。从排队的人角度来看，他们可能不知道车里的人，但车里的人能看到所有排队的人。所以，从他们的角度来看，车里的人被认为是不可接近的。我们称之为逆向可达性（从邻居的角度看你有多远，而不仅仅是你自己的视角）。

和 KNN 一样，你仍然需要定义用于最近邻居数量的k参数。最近邻居是基于观察值之间的距离进行识别的（想想 KNN），然后对每个邻近点计算局部可达密度（LRD或简称局部密度）。这个局部密度是用于比较第k个邻近观察值的评分，那些局部密度低于第k个邻居的点被视为异常值（它们离邻居的范围更远）。

如何操作……

在这个实例中，你将继续使用在技术要求部分创建的tx DataFrame，使用 PyOD 中的LOF类来检测异常值：

  1. 首先加载LOF类：
from pyod.models.lof import LOF
  1. 你应该熟悉几个控制算法行为的参数。第一个参数是contamination，它是一个数值（浮动型），表示数据集中异常值的比例。例如，0.1表示你预计 10%的数据是异常值。默认值是contamination=0.1。在本例中，你将使用0.03（3%）。

第二个参数是邻居的数量，默认为n_neighbors=5，类似于 KNN 算法。理想情况下，你会使用不同的k（n_neighbors）值运行不同的模型，并比较结果以确定最佳的邻居数。最后，metric参数指定用来计算距离的度量。可以使用 scikit-learn 或 SciPy 库中的任何距离度量（例如，欧几里得距离或曼哈顿距离）。默认值是闵可夫斯基距离，metric='minkowski'。由于闵可夫斯基距离是欧几里得距离（https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file282.png）和曼哈顿距离（https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file283.png）的推广，你会看到一个p参数。默认情况下，p=2表示欧几里得距离，而p=1表示曼哈顿距离。

  1. 通过更新n_neighbors=5和contamination=0.03来实例化 LOF，同时保持其他参数为默认值。然后，训练（拟合）模型：
lof = LOF(contamination=0.03, n_neighbors=5)
lof.fit(tx)
>>
LOF(algorithm='auto', contamination=0.03, leaf_size=30, metric='minkowski',
  metric_params=None, n_jobs=1, n_neighbors=5, novelty=True, p=2)
  1. predict方法将为每个数据点输出1或0。值为1表示异常值。将结果存储在一个 pandas Series 中：
predicted = pd.Series(lof.predict(tx),
                      index=tx.index)
print('Number of outliers = ', predicted.sum())
>>
Number of outliers = 6
  1. 过滤预测的序列，仅显示异常值：
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>>       
Timestamp    value
2014-10-31  17473.354167
2014-11-01  20553.500000
2014-12-25  7902.125000
2014-12-26  10397.958333
2015-01-26  7818.979167
2015-01-27  4834.541667

有趣的是，它捕获了五个已知日期中的三个，但成功识别了感恩节后的那一天和圣诞节后的那一天为异常值。此外，10 月 31 日是星期五，那天是万圣节夜晚。

使用在技术要求部分创建的plot_outliers函数来可视化输出，以获得更好的洞察：

plot_outliers(outliers, tx, 'LOF')

前面的代码应生成一个类似于图 14.1的图表，只不过x标记是基于使用 LOF 算法识别的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file284.jpg

图 14.5：使用 LOF 算法识别的潜在异常值标记

若要打印带有标记的标签（日期），只需再次调用plot_outliers函数，但这次传入labels=True：

plot_outliers(outliers, tx, 'LOF', labels=True)

前面的代码应生成一个类似于图 14.5的图表，并附加文本标签。

它是如何工作的…

LOF是一种基于密度的算法，它假设异常点比邻居更加孤立，且具有较低的局部密度得分。

LOF 类似于 KNN，因为我们在计算局部密度之前，测量邻居之间的距离。局部密度是决策得分的基础，你可以通过decision_scores_属性查看这些得分：

timestamp  score
2014-11-01  14.254309
2015-01-27  5.270860
2015-01-26  3.988552
2014-12-25  3.952827
2014-12-26  2.295987
2014-10-31  2.158571

这些得分与图 14.3中的 KNN 得分非常不同。

若要更深入了解decision_得分、threshold_或predict_proba，请查看本章的第一篇教程——使用 KNN 检测异常值。

还有更多内容…

类似于 LOF，算法的另一个扩展是基于聚类的局部异常因子（CBLOF）。CBLOF 在概念上与 LOF 相似，因为它在计算得分以确定异常值时依赖于聚类大小和距离。因此，除了 LOF 中的邻居数（n_neighbors），我们现在有了一个新的参数，即聚类数（n_clusters）。

PyOD 中的默认聚类估计器clustering_estimator是 K 均值聚类算法。

你将使用 PyOD 中的 CBLOF 类，并保持大部分参数为默认值。更改n_clusters=8和contamination=0.03参数：

from pyod.models.cblof import CBLOF
cblof = CBLOF(n_clusters=4, contamination=0.03)
cblof.fit(tx)
predicted = pd.Series(lof.predict(tx),
                      index=tx.index)
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
plot_outliers(outliers, tx, 'CBLOF')

前面的代码应生成一个类似于图 14.1的图表，只不过x标记是基于使用 CBLOF 算法识别的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file285.jpg

图 14.6：使用 CBLOF 算法识别的潜在异常值标记

将图 14.6与图 14.5（LOF）进行比较，注意它们之间的相似性。

另见

要了解更多关于 LOF 和 CBLOF 算法的信息，你可以访问 PyOD 文档：

使用 iForest 检测异常值

iForest 与另一个流行的算法随机森林有相似之处。随机森林是一种基于树的监督学习算法。在监督学习中，你有现有的标签（分类）或值（回归）来表示目标变量。这就是算法学习的方式（它是监督学习）。

森林这个名称来源于算法工作原理的底层机制。例如，在分类中，算法随机采样数据来构建多个弱分类器（较小的决策树），这些分类器共同做出预测。最终，你得到一个由较小树（模型）组成的森林。这种技术优于单个可能会过拟合数据的复杂分类器。集成学习是多个弱学习者协作产生最优解的概念。

iForest，作为一种集成学习方法，是随机森林的无监督学习方法。iForest 算法通过随机划分（拆分）数据集成多个分区来隔离异常值。这个过程是递归进行的，直到所有数据点都属于一个分区。隔离一个异常值所需的分区数量通常比隔离常规数据点所需的分区数量要少。这个思路是，异常数据点距离其他点较远，因此更容易被分离（隔离）。

相比之下，正常的数据点可能会更接近较大的数据集，因此需要更多的分区（拆分）来隔离该点。因此，称之为隔离森林，因为它通过隔离来识别异常值。一旦所有的点都被隔离，算法会生成一个异常值评分。你可以把这些分区看作是创建了一个决策树路径。到达某个点的路径越短，异常的可能性越大。

如何实现…

在本例中，你将继续使用nyc_taxi数据框，利用 PyOD 库中的IForest类来检测异常值：

  1. 开始加载IForest类：
from pyod.models.iforest import IForest
  1. 有一些参数是你应该熟悉的，以控制算法的行为。第一个参数是contamination。默认值为contamination=0.1，但在本例中，你将使用0.03（3%）。

第二个参数是 n_estimators，默认为 n_estimators=100，即生成的随机树的数量。根据数据的复杂性，您可能希望将该值增大到更高的范围，如 500 或更高。从默认的小值开始，以了解基准模型的工作原理——最后，random_state 默认为 None。由于 iForest 算法会随机生成数据的划分，因此设置一个值有助于确保工作结果的可复现性。这样，当您重新运行代码时，能够获得一致的结果。当然，这个值可以是任何整数。

  1. 实例化 IForest 并更新 contamination 和 random_state 参数。然后，将该类的新实例（iforest）拟合到重采样数据上，以训练模型：
iforest = IForest(contamination=0.03,
                 n_estimators=100,
                 random_state=0)
iforest.fit(nyc_daily)
>>
IForest(behaviour='old', bootstrap=False, contamination=0.03,
    max_features=1.0, max_samples='auto', n_estimators=100, n_jobs=1,
    random_state=0, verbose=0)
  1. 使用 predict 方法来识别异常值。该方法会为每个数据点输出 1 或 0。例如，值为 1 表示是一个异常值。

让我们将结果存储在一个 pandas Series 中：

predicted = pd.Series(iforest.predict(tx),
                      index=tx.index)
print('Number of outliers = ', predicted.sum())
>>
Number of outliers =  7

有趣的是，与之前的算法 使用 KNN 检测异常值 不同，iForest 检测到了 7 个异常值，而 KNN 算法检测到了 6 个异常值。

Filter the predicted Series to only show the outlier values:
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>>     
timestamp  value
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

总体而言，iForest 捕获了已知的五个异常值中的四个。还有一些其他有趣的日期被识别出来，这些日期应触发调查，以确定这些数据点是否为异常值。例如，2014 年 11 月 8 日被算法检测为潜在异常值，而该日期并未被考虑在数据中。

  1. 使用 技术要求 部分中创建的 plot_outliers 函数来可视化输出，以便更好地理解：
plot_outliers(outliers, tx, 'IForest')

上述代码应生成类似于 图 14.1 中的图表，唯一不同的是 x 标记基于使用 iForest 算法识别的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file286.jpg

图 14.7：使用 iForest 算法标识潜在异常值的标记

要打印带有标记的标签（日期），只需再次调用 plot_outliers 函数，但这次将 labels=True：

plot_outliers(outliers, tx, 'IForest', labels=True)

上述代码应生成类似于 图 14.7 的图表，并且增加了文本标签。

它是如何工作的…

由于 iForest 是一种集成方法，您将创建多个模型（决策树学习器）。n_estimators 的默认值是 100。增加基础估计器的数量可能会提高模型的性能，但在某个程度上可能会影响计算性能。因此，您可以将估计器数量视为训练好的模型。例如，对于 100 个估计器，您实际上是创建了 100 个决策树模型。

还有一个值得一提的参数是bootstrap参数。默认值为False，它是一个布尔值。由于 iForest 会随机抽样数据，你有两个选择：带替换的随机抽样（称为自助抽样）或不带替换的随机抽样。默认行为是没有替换的抽样。

还有更多…

PyOD 中的 iForest 算法（IForest类）是 scikit-learn 中IsolationForest类的封装。这对上一食谱中使用的 KNN 也是如此，使用 KNN 检测异常值。

让我们进一步探索，使用 scikit-learn 实现 iForest 算法。你将使用fit_predict()方法作为一步训练和预测，这个方法也可以在 PyOD 的各种算法实现中找到：

from sklearn.ensemble import IsolationForest
sk_iforest = IsolationForest(contamination=0.03)
sk_prediction = pd.Series(sk_iforest.fit_predict(tx),
                      index=tx.index)
sk_outliers = sk_prediction[sk_prediction == -1]
sk_outliers = tx.loc[sk_outliers.index]
sk_outliers
>> 
timestamp   value
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

结果是一样的。但请注意，与 PyOD 不同，识别出的异常值在标记时为-1，而在 PyOD 中，异常值被标记为1。

参见

PyOD 的 iForest 实现实际上是 scikit-learn 中IsolationForest类的封装：

使用单类支持向量机（OCSVM）检测异常值

支持向量机（SVM） 是一种流行的有监督机器学习算法，主要用于分类，但也可以用于回归。SVM 的流行源于使用核函数（有时称为核技巧），例如线性、多项式、基于半径的函数（RBF）和 sigmoid 函数。

除了分类和回归，SVM 还可以以无监督的方式用于异常值检测，类似于 KNN。KNN 通常被认为是一种有监督的机器学习技术，但在异常值检测中它是以无监督的方式使用的，正如在使用 KNN 进行异常值检测食谱中所见。

如何做…

在这个食谱中，你将继续使用在技术要求部分创建的tx数据框，利用 PyOD 中的ocsvm类检测异常值：

  1. 首先加载OCSVM类：
from pyod.models.ocsvm import OCSVM
  1. 有一些参数你应该了解，以控制算法的行为。第一个参数是contamination。默认值为contamination=0.1，在这个食谱中，你将使用0.03（即 3%）。

第二个参数是 kernel，其值设为 rbf，你将保持其不变。

通过更新 contamination=0.03 来实例化 OCSVM，同时保持其余参数为默认值。然后，训练（拟合）模型：

ocsvm = OCSVM(contamination=0.03, kernel='rbf')
ocsvm.fit(tx)
>>
OCSVM(cache_size=200, coef0=0.0, contamination=0.03, degree=3, gamma='auto',
   kernel='rbf', max_iter=-1, nu=0.5, shrinking=True, tol=0.001,
   verbose=False)
  1. predict 方法将为每个数据点输出 1 或 0。值为 1 表示异常值。将结果存储在一个 pandas Series 中：
predicted = pd.Series(ocsvm.predict(tx),
                      index=tx.index)
print('Number of outliers = ', predicted.sum())
>>
Number of outliers =  5
  1. 筛选预测的 Series，仅显示异常值：
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>> 
timestamp  value
2014-08-09  15499.708333
2014-11-18  15499.437500
2014-11-27  10899.666667
2014-12-24  12502.000000
2015-01-05  12502.750000

有趣的是，它捕捉到了五个已知日期中的一个。

  1. 使用技术要求部分中创建的 plot_outliers 函数可视化输出，以便获得更好的洞察：
plot_outliers(outliers, tx, 'OCSVM')

上述代码应产生一个与图 14.1类似的图表，唯一不同的是 x 标记是基于 OCSVM 算法识别的异常值：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file287.jpg

图 14.8：使用 OCSVM 标记每个异常点的折线图

当查看图 14.8中的图表时，不清楚为什么 OCSVM 识别出这些日期为异常值。RBF 核函数可以捕捉非线性关系，因此它应该是一个强健的核函数。

这个不准确的原因在于 SVM 对数据缩放敏感。为了获得更好的结果，您需要首先对数据进行标准化（缩放）。

  1. 让我们解决这个问题，先对数据进行标准化，然后再次运行算法：
from pyod.utils.utility import standardizer
scaled = standardizer(tx)
predicted = pd.Series(ocsvm.fit_predict(scaled),
                      index=tx.index)
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>>
timestamp  value
2014-07-06  11464.270833
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

有趣的是，现在模型识别出了五个已知异常日期中的四个。

  1. 在新的结果集上使用 plot_outliers 函数：
plot_outliers(outliers, tx, 'OCSVM Scaled'))

上述代码应生成一个更合理的图表，如下图所示：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file288.jpg

图 14.9：使用标准化函数缩放数据后的 OCSVM

比较图 14.9和图 14.8中的结果，看看缩放如何在 OCSVM 算法识别异常值方面产生了显著差异。

工作原理…

PyOD 对 OCSVM 的实现是对 scikit-learn 中 OneClassSVM 实现的封装。

与 SVM 相似，OneClassSVM 对异常值以及数据的缩放非常敏感。为了获得合理的结果，在训练模型之前标准化（缩放）数据非常重要。

还有更多…

让我们探讨不同核函数在相同数据集上的表现。在下面的代码中，您将测试四种核函数：'linear'、'poly'、'rbf' 和 'sigmoid'。

回顾一下，当使用 SVM 时，您需要缩放数据。您将使用之前创建的缩放数据集：

for kernel in ['linear', 'poly', 'rbf', 'sigmoid']:
    ocsvm = OCSVM(contamination=0.03, kernel=kernel)
    predict = pd.Series(ocsvm.fit_predict(scaled),
                      index=tx.index, name=kernel)
    outliers = predict[predict == 1]
    outliers = tx.loc[outliers.index]
    plot_outliers(outliers, tx, kernel, labels=True)

上述代码应生成每个核函数的图表，以便您可以直观地检查并比较它们之间的差异：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file289.jpg

图 14.10：比较不同核函数与 OCSVM 的效果

有趣的是，每种核方法捕获的异常值略有不同。您可以重新运行之前的代码，通过传递 labels=True 参数来打印每个标记（日期）的标签。

参见

要了解有关 OCSVM 实现的更多信息，请访问官方文档：pyod.readthedocs.io/en/latest/pyod.models.html#module-pyod.models.ocsvm。

使用 COPOD 检测异常值

COPOD 是一个令人兴奋的算法，基于 2020 年 9 月发布的一篇论文，你可以在这里阅读：arxiv.org/abs/2009.09463。

PyOD 库提供了许多基于最新研究论文的算法，这些算法可以分为线性模型、基于邻近的模型、概率模型、集成模型和神经网络。

COPOD 属于概率模型，并被标记为 无参数 算法。它唯一的参数是 contamination 因子，默认为 0.1。COPOD 算法受统计方法启发，使其成为一个快速且高度可解释的模型。该算法基于 copula，一种通常用于建模相互独立的随机变量之间的依赖关系的函数，这些变量不一定服从正态分布。在时间序列预测中，copula 已被应用于单变量和多变量预测，这在金融风险建模中非常流行。copula 这个术语源自 copula 函数，它将单变量的边际分布连接（耦合）在一起，形成一个统一的多变量分布函数。

如何操作…

在本食谱中，您将继续使用 tx DataFrame，通过 PyOD 库中的 COPOD 类来检测异常值：

  1. 首先加载 COPOD 类：
from pyod.models.copod import COPOD
  1. 您需要考虑的唯一参数是 contamination。通常，将该参数（用于所有异常值检测实现）视为一个阈值，用于控制模型的敏感性并最小化假阳性。由于这是一个由您控制的参数，理想情况下，您希望运行多个模型，实验出适合您用例的理想阈值。

如需了解更多关于 decision_scores_、threshold_ 或 predict_proba 的信息，请查看本章的第一个食谱，使用 KNN 检测异常值。

  1. 实例化COPOD并将contamination更新为0.03。然后，在重新采样的数据上进行拟合，以训练模型：
copod = COPOD(contamination=0.03)
copod.fit(tx)
>>
COPOD(contamination=0.03, n_jobs=1)
  1. 使用 predict 方法识别异常值。该方法将为每个数据点输出 1 或 0。例如，1 表示异常值。

将结果存储在 pandas Series 中：

predicted = pd.Series(copod.predict(tx),
                      index=tx.index)
print('Number of outliers = ', predicted.sum())
>>
Number of outliers =  7

异常值的数量与使用 iForest 获得的数量相匹配。

  1. 仅筛选预测的 Series，以显示异常值：
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>>           
timestamp  value
2014-07-04  11511.770833
2014-07-06  11464.270833
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

与你迄今为止探索的其他算法相比，你会注意到使用 COPOD 捕获了一些有趣的异常点，这些异常点之前没有被识别出来。例如，COPOD 识别了 7 月 4 日——美国的国庆节（独立日）。那天恰好是周末（星期五是休息日）。COPOD 模型在 7 月 4 日和 7 月 6 日的整个周末期间捕获了异常。恰巧 7 月 6 日是由于纽约的一场棒球赛而成为一个有趣的日子。

  1. 使用技术要求部分中创建的plot_outliers函数来可视化输出，以便获得更好的洞察：
plot_outliers(outliers, tx, 'COPOD')

上述代码应该生成一个类似于图 14.1的图表，唯一的区别是x标记基于使用 COPOD 算法识别的异常点：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file290.jpg

图 14.11：使用 COPOD 算法识别的潜在异常点的标记

要打印带有标记的标签（日期），只需再次调用plot_outliers函数，但这次需要将labels=True：

plot_outliers(outliers, tx, 'COPOD', labels=True)

上述代码应该生成一个类似于图 14.11的图表，并加上文本标签。

它是如何工作的…

COPOD 是一个先进的算法，但它仍然基于概率建模和在数据中找到统计学上显著的极端值。使用 COPOD 的几项测试已经证明它在基准数据集上的卓越表现。使用 COPOD 的一个吸引力是它不需要调参（除了污染因子）。因此，作为用户，你无需担心超参数调优。

还有更多…

另一个简单且流行的概率算法是中位绝对偏差（MAD）。我们在第八章，使用统计方法进行异常值检测中探讨了 MAD，具体是在使用修改过的 z-score 进行异常值检测的食谱中，你是从零开始构建该算法的。

这是 PyOD 提供的一个类似实现，只有一个参数：阈值。如果你还记得第八章，使用统计方法进行异常值检测，阈值是基于标准差的。

以下代码展示了如何使用 PyOD 实现 MAD。你将使用threshold=3来复现你在第八章，使用统计方法进行异常值检测中的操作：

from pyod.models.mad import MAD
mad = MAD(threshold=3)
predicted = pd.Series(mad.fit_predict(tx),
                      index=tx.index)
outliers = predicted[predicted == 1]
outliers = tx.loc[outliers.index]
outliers
>> 
timestamp  value
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

这应该与你在第八章，使用统计方法进行异常值检测中，修改过的 z-score 实现的结果一致。

另见

若要了解更多关于 COPOD 及其在 PyOD 中的实现，请访问官方文档：pyod.readthedocs.io/en/latest/pyod.models.html?highlight=copod#pyod.models.copod.COPOD。

如果你有兴趣阅读COPOD: 基于 Copula 的异常值检测（2020 年 9 月发布）的研究论文，请访问 arXiv.org 页面：arxiv.org/abs/2009.09463。

使用 PyCaret 检测异常值

在本食谱中，你将探索PyCaret用于异常检测。PyCaret (pycaret.org) 被定位为“一个开源的、低代码的 Python 机器学习库，自动化机器学习工作流”。PyCaret 是 PyOD 的封装，你之前在食谱中使用过它进行异常检测。PyCaret 的作用是简化整个过程，用最少的代码进行快速原型设计和测试。

你将使用 PyCaret 来检查多种异常检测算法，类似于你在之前的食谱中使用过的算法，并查看 PyCaret 如何简化这个过程。

准备开始

探索 PyCaret 的推荐方式是为 PyCaret 创建一个新的虚拟 Python 环境，这样它可以安装所有所需的依赖项，而不会与当前环境发生任何冲突或问题。如果你需要快速回顾如何创建虚拟 Python 环境，请参考开发环境设置食谱，见第一章，时间序列分析入门。本章介绍了两种方法：使用conda和venv。

以下说明将展示使用conda的过程。你可以为环境命名任何你喜欢的名字；在以下示例中，我们将命名为pycaret：

>> conda create -n pycaret python=3.8 -y
>> conda activate pycaret
>> pip install "pycaret[full]"

为了使新的pycaret环境在 Jupyter 中可见，你可以运行以下代码：

python -m ipykernel install --user --name pycaret --display-name "PyCaret"

本食谱有一个单独的 Jupyter notebook，你可以从 GitHub 仓库下载：

github.com/PacktPublishing/Time-Series-Analysis-with-Python-Cookbook./blob/main/code/Ch14/Chapter%2014-pycaret.ipynb

如何操作…

在本食谱中，你将不会接触到任何新概念。重点是展示如何使用 PyCaret 作为实验的起点，并快速评估不同的模型。你将加载 PyCaret 并运行不同的异常检测算法：

  1. 从pycaret.anomaly模块加载所有可用的函数：
from pycaret.anomaly import *
setup = setup(tx, session_id = 1, normalize=True)

上述代码应该生成一个表格摘要，如图 14.12 所示

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file291.png

图 14.12 – PyCaret 摘要输出

  1. 要打印可用的异常检测算法列表，可以运行models()：
models()

这应该会显示一个 pandas DataFrame，如下所示：

https://github.com/OpenDocCN/freelearn-ds-pt4-zh/raw/master/docs/ts-anal-py-cb-2e/img/file292.png

图 14.14: PyCaret 可用的异常检测算法

请注意，这些算法都来自 PyOD 库。如前所述，PyCaret 是 PyOD 和其他库（如 scikit-learn）的封装。

  1. 让我们将前八个算法的名称存储在列表中，以便稍后使用：
list_of_models = models().index.tolist()[0:8]
list_of_models
>>
['abod', 'cluster', 'cof', 'iforest', 'histogram', 'knn', 'lof', 'svm']
  1. 遍历算法列表并将输出存储在字典中，以便稍后在分析中引用。要在 PyCaret 中创建模型，你只需使用create_model()函数。这类似于 scikit-learn 和 PyOD 中的fit()函数，用于训练模型。一旦模型创建完成，你可以使用该模型通过predict_model()函数预测（识别）离群值。PyCaret 将生成一个包含三列的 DataFrame：原始的value列，一个新的Anomaly列，存储结果为0或1，其中1表示离群值，另一个新的Anomaly_Score列，存储使用的得分（得分越高，表示越有可能是离群值）。

你只需要改变污染参数，以便与之前使用 PyOD 的配方匹配。在 PyCaret 中，污染参数被称为fraction，为了保持一致性，你需要将其设置为0.03或者 3%，即fraction=0.03：

results = {}
for model in list_of_models:
    cols = ['value', 'Anomaly_Score']
    outlier_model = create_model(model, fraction=0.03)
    print(outlier_model)
    outliers = predict_model(outlier_model, data=tx)
    outliers = outliers[outliers['Anomaly'] == 1][cols]
    outliers.sort_values('Anomaly_Score', ascending=False, inplace=True)
    results[model] = {'data': outliers, 'model': outlier_model}

results字典包含每个模型的输出（一个 DataFrame）。

  1. 要打印每个模型的离群值，你可以简单地遍历字典：
for model in results:
    print(f'Model: {model}')
    print(results[model]['data'], '\n')

这将打印出每个模型的结果。以下是列表中的前两个模型作为示例：

Model: abod
                   value  Anomaly_Score
timestamp                             
2014-11-01  20553.500000      -0.002301
2015-01-27   4834.541667      -0.007914
2014-12-26  10397.958333      -3.417724
2015-01-26   7818.979167    -116.341395
2014-12-25   7902.125000    -117.582752
2014-11-27  10899.666667    -122.169590
2014-10-31  17473.354167   -2239.318906
Model: cluster
                   value  Anomaly_Score
timestamp                             
2015-01-27   4834.541667       3.657992
2015-01-26   7818.979167       2.113955
2014-12-25   7902.125000       2.070939
2014-11-01  20553.500000       0.998279
2014-12-26  10397.958333       0.779688
2014-11-27  10899.666667       0.520122
2014-11-28  12850.854167       0.382981

它是如何工作的…

PyCaret 是一个出色的自动化机器学习库，最近他们在时间序列分析、预测和离群值（异常值）检测方面不断扩展其功能。PyCaret 是 PyOD 的封装库，你在本章之前的配方中也使用了 PyOD。图 14.14展示了 PyCaret 支持的 PyOD 算法数量，这是 PyOD 更广泛算法列表的一个子集：pyod.readthedocs.io/en/latest/index.html#implemented-algorithms。

另请参见

要了解更多关于 PyCaret 离群值检测的信息，请访问官方文档：pycaret.gitbook.io/docs/get-started/quickstart#anomaly-detection。

评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值