机器学习算法交易教程第二版（四）

原文：zh.annas-archive.org/md5/a3861c820dde5bc35d4f0200e43cd519

译者：飞龙

协议：CC BY-NC-SA 4.0

第九章：用于波动率预测和统计套利的时间序列模型

在第七章，线性模型 - 从风险因素到资产收益预测中，我们介绍了用于推断和预测的线性模型，从静态模型开始，考虑具有即时影响的横截面输入的同时关系。我们介绍了普通最小二乘法（OLS）学习算法，并且发现它对于正确指定的模型产生无偏系数，且残差与输入变量不相关。假设残差具有恒定方差，可以保证 OLS 在无偏估计器中产生最小的均方预测误差。

我们还遇到了既有横截面又有时间序列维度的面板数据，当我们学习如何使用 Fama-Macbeth 回归来估计随时间和跨资产的风险因素的价值时。然而，随着时间的推移，收益之间的关系通常相当低，因此这个过程可能会大部分忽略时间维度。

此外，我们介绍了正则化的岭回归和套索回归模型，它们产生了偏差系数估计，但可以减少均方预测误差。这些预测模型采用了更加动态的视角，将历史收益与其他输入结合起来，以预测未来的回报。

在本章中，我们将建立动态线性模型，以明确表示时间，并包括在特定间隔或滞后观察到的变量。时间序列数据的一个关键特征是它们的顺序：与横截面数据中的个体观测随机样本不同，我们的数据是一个无法重复的随机过程的单个实现。

我们的目标是识别时间序列中的系统模式，以帮助我们预测时间序列在未来的行为。更具体地说，我们将专注于从输出的历史序列中提取信号的模型，并且可以选择其他同时或滞后的输入变量来预测输出的未来值。例如，我们可以尝试使用过去的回报，结合基准或宏观经济变量的历史回报，来预测股票的未来回报。在转向第四部分之前，我们将专注于线性时间序列模型，然后转向非线性模型，如循环或卷积神经网络。

时间序列模型在交易中具有内在的时间维度，因此非常受欢迎。主要应用包括资产收益和波动率的预测，以及资产价格序列的共同波动的识别。随着越来越多连接设备收集具有潜在信号内容的定期测量，时间序列数据可能会变得更加普遍。

我们首先介绍可以用来诊断时间序列特征和提取捕捉潜在模式的特征的工具。然后，我们将介绍如何诊断和实现时间序列平稳性。接下来，我们将介绍单变量和多变量时间序列模型，并将其应用于预测宏观数据和波动性模式。最后，我们将介绍协整概念以及如何将其应用于开发配对交易策略。

具体来说，我们将涵盖以下主题：

• 如何使用时间序列分析来准备和指导建模过程

• 估计和诊断单变量自回归和移动平均模型

• 构建自回归条件异方差（ARCH）模型以预测波动性

• 如何构建多元向量自回归模型

• 使用协整来开发配对交易策略

您可以在 GitHub 仓库的相应目录中找到本章的代码示例和其他资源链接。这些笔记本包含了图片的彩色版本。关于本章主题的投资角度的全面介绍，请参阅 Tsay（2005 年）和 Fabozzi、Focardi 和 Kolm（2010 年）。

诊断和特征提取工具

时间序列是一系列以离散间隔分隔的值，通常是均匀间隔的（除了缺失值）。时间序列通常被建模为一个随机过程，由一组随机变量，，组成，每个时间点有一个变量，。单变量时间序列由每个时间点上的单个值y组成，而多变量时间序列由可以用向量表示的多个观察值组成。

两个不同时间点t[i]、t[j]之间的周期数，，称为滞后，每个时间序列有T-1 个不同的滞后。正如交叉模型中不同变量在给定时间点的关系对于横截面模型至关重要一样，分隔给定滞后的数据点之间的关系对于分析和利用时间序列中的模式至关重要。

对于横截面模型，我们用标签y和x来区分输入和输出变量，或目标和预测变量。在时间序列的情况下，结果y的一些或所有滞后值扮演横截面情况下输入或x值的角色。

如果一个时间序列是独立同分布（IID）的随机变量序列，且具有有限均值和方差，则称其为白噪声。特别地，如果这个序列的随机变量服从均值为零、方差为常数的正态分布，则称其为高斯白噪声。

如果时间序列可以写成过去扰动的加权和，，也称为创新，并且在此假定为代表白噪声的扰动序列的平均值，那么时间序列是线性的：

时间序列分析的一个关键目标是了解由系数驱动的动态行为，。时间序列分析提供了针对这种类型数据量身定制的方法，目的是提取有用的模式，进而帮助我们建立预测模型。

我们将介绍实现这一目的的最重要工具，包括分解为关键的系统元素、自相关分析以及滚动窗口统计，如移动平均值。

对于本章中的大多数示例，我们将使用联邦储备提供的数据，您可以使用我们在 第二章，市场和基本数据-来源和技术 中介绍的 pandas-datareader 访问。本节的代码示例在笔记本 tsa_and_stationarity 中可用。

如何分解时间序列模式

时间序列数据通常包含多种模式，可以将其分解为几个组件。特别是，时间序列通常将趋势、季节性和周期性等系统性组件与非系统性噪声结合在一起。这些组件可以建模为线性组合（例如，当波动不依赖于系列水平时）或非线性乘法形式。

根据模型假设，它们也可以自动拆分。Statsmodels 包括一种简单的方法，使用移动平均值将时间序列拆分为单独的趋势、季节性和残差组件。我们可以将其应用于包含强趋势组件和季节性组件的月度工业制造数据，如下所示：

import statsmodels.tsa.api as tsa
industrial_production = web.DataReader('IPGMFN', 'fred', '1988', '2017-12').squeeze()
components = tsa.seasonal_decompose(industrial_production, model='additive')
ts = (industrial_production.to_frame('Original')
      .assign(Trend=components.trend)
      .assign(Seasonality=components.seasonal)
      .assign(Residual=components.resid))
ts.plot(subplots=True, figsize=(14, 8)); 

图 9.1 显示了显示加法组件的结果图表。残差组件将成为后续建模工作的重点，假设趋势和季节性组件更具确定性并且更易于简单外推：

图 9.1：时间序列分解为趋势、季节性和残差

还有更复杂的基于模型的方法，例如 Hyndman 和 Athanasopoulos（2018）中的 第六章，机器学习过程。

滚动窗口统计和移动平均线

鉴于时间序列数据的顺序排列，自然而然地可以计算给定长度期间的熟悉描述性统计。目标是检测系列是否稳定或随时间变化，并获得捕获系统方面的平滑表示，同时过滤噪声。

滚动窗口统计信息服务于此过程：它们生成一个新的时间序列，其中每个数据点代表原始数据的一定期间内的摘要统计量。移动平均是最熟悉的例子。原始数据点可以以相等的权重进入计算，也可以，例如，强调更近期的数据点。指数移动平均递归地计算权重，这些权重对于更早的数据点衰减。新数据点通常是所有先前数据点的摘要，但也可以从周围的窗口计算得出。

Pandas 库包括滚动或扩展窗口，并允许使用各种权重分布。在第二步中，您可以对每个窗口捕获的数据集应用计算。这些计算包括用于单个系列的内置函数，如均值或总和，以及用于多个系列的相关性或协方差，以及用户定义的函数。

我们在第四章、金融特征工程 - 如何研究阿尔法因子和第七章、线性模型 - 从风险因素到收益预测中使用了这个功能来设计特征。例如，下一节中的移动平均和指数平滑示例也将应用这些工具。

早期的预测模型包括带有指数权重的移动平均模型，称为指数平滑模型。我们将再次遇到移动平均作为线性时间序列的关键构建模块。依赖指数平滑方法的预测使用过去观察值的加权平均值，其中权重随着观察值变老而指数衰减。因此，较近期的观察值会获得较高的相关权重。这些方法在时间序列没有非常复杂或突然的模式时很受欢迎。

如何测量自相关

自相关（也称为串行相关）将相关性的概念应用到时间序列的背景下：就像相关系数衡量两个变量之间线性关系的强度一样，自相关系数，，衡量了相隔给定滞后时间k的时间序列值之间的线性关系的程度：

因此，我们可以计算时间序列的T-1 个滞后期中的每个自相关系数。自相关函数（ACF）将相关系数作为滞后期的函数计算。

超过 1 个滞后（即，两个观察值相隔超过一个时间步长）的自相关反映了这些观察值之间的直接相关性以及介于其间数据点的间接影响。偏自相关消除了这种影响，仅测量给定滞后距离T处数据点之间的线性依赖关系。消除意味着使用线性回归的残差，其中结果x[t]和滞后值x[t][-1]、x[t][-2]，…，x[T][-1]作为特征（也称为AR(T-1)模型，我们将在下一节关于单变量时间序列模型中讨论此模型）。偏自相关函数（PACF）提供了一旦较短滞后的相关性效应被消除后所得到的所有相关性，如前述所述。

还有一些算法根据样本自相关来估计偏自相关，这些算法基于 PACF 和 ACF 之间的精确理论关系。

自相关图简单地是顺序滞后k=0,1,...,n的 ACF 或 PACF 的图。它允许我们一眼看出不同滞后的相关结构（有关示例，请参见图 9.3）。自相关图的主要用途是在去除确定性趋势或季节性后检测任何自相关性。ACF 和 PACF 都是线性时间序列模型设计的关键诊断工具，我们将在下一节关于时间序列变换的示例中回顾 ACF 和 PACF 图。

如何诊断和实现平稳性

平稳时间序列的统计特性，如均值、方差或自相关，与周期无关—即，它们随时间不变。因此，平稳性意味着时间序列没有趋势或季节效应。此外，它要求在不同滚动窗口计算的描述性统计量，如均值或标准差，是恒定的或随时间变化不显著。平稳时间序列会回归到其均值，并且偏差具有恒定的振幅，而短期波动在统计意义上始终是相同的。

更正式地说，严格平稳性要求时间序列观测的任何子集的联合分布在所有时刻上都与时间无关。因此，除了均值和方差外，高阶矩如偏度和峰度也需要保持恒定，而不受不同观测值之间滞后的影响。在大多数应用中，例如本章中大多数可用于建模资产收益率的时间序列模型中，我们将平稳性限制在一阶和二阶矩上，使得时间序列具有恒定的均值、方差和自相关性。然而，在构建建模波动率时，我们放弃了这一假设，并明确假设方差会以可预测的方式随时间变化。

请注意，我们特别允许输出值在不同滞后期之间存在依赖关系，就像我们希望线性回归的输入数据与结果相关一样。平稳性意味着这些关系是稳定的。平稳性是经典统计模型的关键假设。下面的两个小节介绍了可以帮助使时间序列平稳的转换，以及如何处理由单位根引起的随机趋势的特殊情况。

转换时间序列以实现平稳性

为了满足许多时间序列模型的平稳性假设，我们需要对原始序列进行转换，通常需要几个步骤。常见的转换包括（自然）对数，将指数增长模式转换为线性趋势并稳定方差。通货紧缩意味着将一个时间序列除以另一个导致趋势行为的序列，例如，将一个名义序列除以一个价格指数以将其转换为实际度量。

如果一个序列是趋势稳定的，那么它将恢复到一个稳定的长期线性趋势。通常可以通过使用线性回归拟合趋势线并使用残差来使其平稳。这意味着在线性回归模型中包括时间指数作为独立变量，可能结合对数化或通货紧缩。

在许多情况下，去趋势化不足以使序列平稳。相反，我们需要将原始数据转换为一系列周期间和/或季节间差异。换句话说，我们使用相邻数据点或季节滞后的值相减的结果。请注意，当这种差分应用于对数变换序列时，结果代表的是金融背景下的瞬时增长率或收益率。

如果单变量序列经过d次差分后变得平稳，就称为d阶整合，或者简称为一阶整合如果d=1。这种行为是由单位根引起的，接下来我们将解释。

处理而不是如何处理

单位根对确定使时间序列平稳的转换方法构成特殊问题。在讨论诊断测试和解决方案之前，我们将首先解释单位根的概念。

关于单位根和随机游走

时间序列经常被建模为以下自回归形式的随机过程，以使当前值成为过去值的加权和，再加上一个随机扰动：

在下一节中，我们将更详细地探讨这些模型作为 ARIMA 模型中 AR 的构建块的单变量时间序列模型。这样的过程具有以下形式的特征方程：

如果这个多项式的（最多）p个根之一等于 1，则过程被称为具有单位根。它将是非平稳的，但不一定会有趋势。如果特征方程的剩余根的绝对值小于 1，则该过程的一阶差分将是平稳的，并且过程是一阶整合的或 I(1)。如果额外的根的绝对值大于 1，则积分阶数更高，并且需要额外的差分。

在实践中，利率或资产价格的时间序列通常不是平稳的，因为没有价格水平可供系列回归到。非平稳序列的最突出例子是随机漫步。给定具有初始价格 p[0]（例如，股票的 IPO 价格）和白噪声扰动 的价格时间序列 p[t]，随机漫步满足以下自回归关系：

重复替换表明当前值 p[t] 是所有先前扰动或创新，，以及初始价格 p[0] 的总和。如果方程包括一个常数项，那么随机漫步被认为有漂移。

随机漫步因此是以下形式的自回归随机过程：

它的特征方程为 ，具有单位根，既是非平稳的，也是一阶整合的。一方面，鉴于 的 IID 特性，时间序列的方差等于 ，这是不是二阶平稳的，这意味着原则上，该序列可以随时间取任何值。另一方面，进行第一阶差分，得到 ，留下了 平稳的序列 ，鉴于对 的统计假设。

具有单位根的非平稳序列的定义特征是长期记忆：由于当前值是过去扰动的总和，大的创新比均值回归、平稳序列持续时间更长。

如何诊断单位根

统计单位根检验是一种确定是否需要（额外）差分的常见方法。这些是用于确定是否需要差分的平稳性的统计假设检验。

增广的迪基-富勒检验（ADF 检验）评估一个时间序列样本是否具有单位根的零假设，对立假设是平稳性。它将差分时间序列回归到时间趋势、第一滞后和所有滞后差分上，并从滞后时间序列值的系数值计算出一个检验统计量。statsmodels使其易于实现（参见笔记本tsa_and_stationarity）。

形式上，对于一个时间序列的 ADF 检验，，进行线性回归，其中 是一个常数， 是时间趋势上的系数，p 指的是模型中使用的滞后数：

约束 意味着一个随机游走，而仅有 意味着一个带漂移的随机游走。滞后阶数通常使用赤池信息准则（AIC）和贝叶斯信息准则（BIC）信息准则来决定，这些准则在 第七章，线性模型 - 从风险因素到回报预测 中介绍。

ADF 检验统计量使用样本系数 ，在单位根非平稳的零假设下，该系数为零，在其他情况下为负。其意图在于证明，对于一个积分序列，滞后序列值不应提供有用信息来预测超过滞后差异的一阶差异。

如何去除单位根并处理得到的序列

除了使用相邻数据点之间的差异来消除恒定变化模式外，我们还可以应用季节性差分来消除季节性变化模式。这涉及到在代表季节模式长度的滞后距离处取值的差异。对于月度数据，这通常涉及到滞后 12，而对于季度数据，它涉及到滞后 4，以消除季节性和线性趋势。

确定正确的转换，特别是适当的差分数量和滞后数，并不总是一目了然的。一些启发式方法已经被提出，可以总结如下：

• 滞后一期自相关接近零或为负，或者自相关普遍较小且无规律：无需进行更高阶的差分处理

• 正自相关延伸至 10+ 滞后：该序列可能需要更高阶的差分

• 滞后一期自相关 < -0.5：该序列可能存在过度差分

• 轻微的过度或不足的差分可以通过 AR 或 MA 项进行校正（请参阅下一节关于单变量时间序列模型）

一些作者建议使用分数差分作为使积分序列变得平稳的更灵活的方法，并且可能能够保留比简单或季节性差分在离散间隔更多的信息或信号。例如，参见 Marcos Lopez de Prado（2018）的 第五章，投资组合优化和绩效评估。

实践中的时间序列变换

图 9.2 中的图表显示了纳斯达克股票指数和工业生产的时间序列，经过对数变换和随后的应用第一次和季节性差分（在滞后 12 处），分别展示了它们的原始形式以及变换后的版本，跨越 2017 年的 30 年时间。

图表还显示了 ADF p 值，这使我们能够在两种情况下拒绝单位根非平稳性的假设：

图 9.2：时间序列转换和单位根检验结果

我们可以进一步分析使用 Q-Q 图和基于 ACF 和 PACF 的自相关图的转换系列的相关时间序列特征的特征。

对于图 9.3中的 NASDAQ 图表，我们可以看到虽然没有趋势，但方差不是恒定的，而是在 1980 年代末、2001 年和 2008 年的市场动荡期间显示出集中的尖峰。Q-Q 图突出了分布的厚尾，极端值比正态分布所暗示的更频繁。

ACF 和 PACF 显示类似的模式，几个滞后的自相关显得显著：

图 9.3：转换后的 NASDAQ 综合指数的描述统计

对于工业制造生产的月度时间序列，我们可以看到在 2008 年危机后有一个大的负异常值，以及 Q-Q 图中相应的偏度（见图 9.4）。自相关远高于 NASDAQ 收益，并且平稳下降。PACF 在滞后 1 和 13 处显示出明显的正自相关模式，并在滞后 3 和 4 处显示出显著的负系数：

图 9.4：转换后的工业生产数据的描述统计

单变量时间序列模型

多元线性回归模型将感兴趣的变量表示为输入的线性组合，加上一个随机扰动。相比之下，单变量时间序列模型将时间序列的当前值与系列滞后值的线性组合、当前噪声以及可能的过去噪声项相关联。

尽管指数平滑模型基于对数据中趋势和季节性的描述，但ARIMA 模型旨在描述数据中的自相关性。ARIMA(p, d, q) 模型需要平稳性，并利用两个构建块：

• 自回归（AR）项包含时间序列的p滞后值

• 包含q滞后扰动的移动平均（MA）项

I 代表 integrated，因为模型可以通过对系列进行d次差分来考虑单位根非平稳性。自回归一词强调了 ARIMA 模型意味着时间序列对其自身值的回归。

我们将介绍 ARIMA 的构建模块，AR 和 MA 模型，并解释如何将它们组合成可以像 ARIMA 模型那样考虑序列积分，或者包含外生变量像AR(I)MAX模型那样的自回归移动平均(ARMA)模型。此外，我们将说明如何包含季节性 AR 和 MA 项以扩展工具箱，使其还包括SARMAX模型。

如何构建自回归模型

p阶 AR 模型旨在捕捉不同滞后时期时间序列值之间的线性相关性，可以写成如下形式：

这与对y[t]的滞后值进行多元线性回归非常相似。该模型具有以下特征方程：

解这个关于x的p次多项式的倒数是特征根，如果所有根的绝对值都小于 1，则 AR(p)过程是稳定的，否则不稳定。对于稳定的序列，多步预测将收敛于序列的均值。

我们可以使用熟悉的最小二乘法估计模型参数，使用p+1，...，T的观测数据确保每个滞后项和结果都有数据。

如何确定滞后数量

在实践中，挑战在于确定滞后项的适当阶p。我们在如何测量自相关性部分讨论过的串行相关性时间序列分析工具在做出这一决定方面起着关键作用。

更具体地说，对自相关图的目视检查通常提供了有用的线索：

• ACF估计不同滞后时期观测值之间的自相关性，其结果既来自直接线性相关性，也来自间接线性相关性。因此，如果p阶 AR 模型是正确的模型，自相关函数将显示出直至滞后k的显著串行相关性，并且由于线性关系间接效应所造成的惯性，将延伸到随后的滞后直至最终随着效应的减弱而消失。

• PACF只测量了给定滞后间隔的观测值之间的直接线性关系，因此不会反映超过k的滞后的相关性。

如何诊断模型拟合度

如果模型正确捕捉了滞后间的线性依赖关系，则残差应类似于白噪声，而自相关函数应突出显示出显著的自相关系数的缺失。

除了残差图外，Ljung-Box Q 统计量还允许我们测试残差序列是否符合白噪声的假设。零假设是所有m个序列相关系数都为零，而备择假设是某些系数不为零。检验统计量是根据不同滞后k的样本自相关系数计算的，并且遵循X²分布：

正如我们将看到的，statsmodels 提供了关于不同滞后的系数显著性的信息，不显著的系数应该被移除。如果 Q 统计拒绝无自相关的零假设，你应该考虑额外的 AR 项。

如何构建移动平均模型

MA(q) 模型使用 q 个过去扰动，而不是时间序列的滞后值，作为回归模型的一部分，如下所示：

由于我们没有观察到白噪声扰动值，，MA(q) 不像我们迄今为止看到的那样是一个回归模型。MA(q) 模型不是使用最小二乘法估计的，而是使用最大似然估计（MLE），或者在系列开始时初始化或估计扰动，然后递归和迭代地计算剩余值。

MA(q) 模型得名于将 y[t] 的每个值表示为过去 q 个创新的加权移动平均值。换句话说，当前的估计值表示相对于模型过去错误的修正。MA(q) 模型中使用移动平均值的方式与指数平滑或季节时间序列成分的估计不同，因为 MA(q) 模型旨在预测未来值，而不是去噪声或估计过去值的趋势周期。

MA(q) 过程总是平稳的，因为它们是白噪声变量的加权和，这些变量本身就是平稳的。

如何确定滞后数

由 MA(q) 过程生成的时间序列由前 q 个模型预测的残差驱动。因此，MA(q) 过程的 ACF 将显示出直到滞后 q 的值的显著系数，然后急剧下降，因为这是模型假定系列值被生成的方式。

注意这与我们刚刚描述的 AR 情况不同，那里 PACF 会显示类似的模式。

AR 和 MA 模型之间的关系

无论何时，都可以使用重复替换的方法将 AR(p) 模型表示为一个 过程，例如“如何处理由单位根引起的随机趋势”部分中的随机游走示例。

当 MA(q) 过程的系数满足一定的大小约束时，它也变得可逆，并且可以表示为一个 过程（详见 Tsay，2005 年）。

如何构建 ARIMA 模型及其扩展

自回归积分滑动平均—ARIMA(p, d, q)—模型结合了 AR(p) 和 MA(q) 过程，以利用这些基本模型的互补性并简化模型开发。它们使用更紧凑的形式并减少参数数量，进而减少过度拟合的风险。

这些模型还通过使用时间序列值的d^(th)差分来消除单位根的非平稳性。一个 ARIMA(p, 1, q)模型相当于使用序列的第一阶差分的 ARMA(p, q)模型。用y'表示非季节性差分d次后的原始序列，ARIMA(p, d, q)模型简单地表示为：

ARIMA 模型也使用 MLE 进行估计。根据实现方式，更高阶的模型可能一般包含较低阶的模型。

例如，截至版本 0.11，statsmodels 包括所有低阶p和q项，并且不允许删除低于最高值的滞后的系数。在这种情况下，更高阶的模型将始终拟合得更好。使用太多项不要使模型过度拟合数据。在撰写时，最新版本为 0.11，增加了一个带有更灵活配置选项的实验性新 ARIMA 模型。

如何对差分序列建模

使用数据设计单变量时间序列模型时也有一些指导方针：

• 没有差分的模型假设原始序列是平稳的，包括均值回归。它通常包括一个常数项，以允许非零均值。

• 一个一阶差分的模型假设原始序列具有恒定的趋势，因此应包括一个常数项。

• 一个两阶差分的模型假设原始序列具有时间变化的趋势，并且不应包括一个常数。

如何确定 AR 和 MA 项的数量

由于 AR(p)和 MA(q)项相互作用，ACF 和 PACF 提供的信息不再可靠，只能用作起点。

传统上，AIC 和 BIC 信息准则被用于在选择模型设计时依赖样本内拟合。或者，我们可以依赖样本外测试来交叉验证多个参数选择。

在考虑 AR 和 MA 模型时，以下摘要提供了如何选择模型阶数的一些指导：

• PACF 截断的滞后阶数是表示 AR 项的数量。如果差分序列的 PACF 显示出明显的截断和/或滞后 1 的自相关为正，则添加一个或多个 AR 项。

• ACF 截断的滞后阶数是表示 MA 项的数量。如果差分序列的 ACF 显示出明显的截断和/或滞后 1 的自相关为负，则考虑将 MA 项添加到模型中。

• AR 和 MA 项可能会抵消彼此的影响，因此，如果您的模型同时包含两者，请始终尝试将 AR 和 MA 项的数量减少 1，以避免过度拟合，特别是如果更复杂的模型需要超过 10 次迭代才能收敛。

• 如果 AR 系数之和接近于 1，并且建议模型的 AR 部分存在单位根，则消除一个 AR 项并再次进行差分。

• 如果 MA 系数总和接近于一，并且暗示模型的 MA 部分存在单位根，则删除一个 MA 项并将差分阶数减少一个。

• 不稳定的长期预测表明模型的 AR 或 MA 部分可能存在单位根。

添加特征 – ARMAX

具有外生输入的自回归滑动平均模型（ARMAX）模型在 ARMA 时间序列模型的右侧添加输入变量或协变量（假设系列是平稳的，所以我们可以跳过差分）：

这类似于线性回归模型，但解释起来相当困难。这是因为对y[t]的影响并不是x[t]增加一个单位所产生的影响，如线性回归中那样。相反，方程右侧的y[t]的滞后值的存在意味着只有在给定响应变量的滞后值的情况下，系数才能被解释，这几乎是不直观的。

添加季节差分 – SARIMAX

对于具有季节效应的时间序列，我们可以包括捕捉季节性周期性的 AR 和 MA 项。例如，当使用月度数据且季节效应长度为 1 年时，季节性 AR 和 MA 项将反映这个特定的滞后长度。

然后 ARIMAX(p, d, q) 模型变成了 SARIMAX(p, d, q) × (P, D, Q) 模型，这个写法略显复杂，但 statsmodels 文档（GitHub 上的链接）详细提供了这些信息。

我们现在将使用宏观数据构建一个季节性 ARMA 模型，以说明其实现过程。

如何预测宏观基本面

我们将为 1988-2017 年期间的工业生产时间序列构建一个月度数据的 SARIMAX 模型。如在分析工具的第一节中所示，数据已经进行了对数变换，并且我们使用了季节性（滞后 12 个月）差分。我们对一系列普通和常规的 AR 和 MA 参数进行估计，使用了 10 年的训练数据的滚动窗口，并评估了 1 步预测的均方根误差（RMSE），如下所示的简化代码（详细信息请参见笔记本arima_models）：

for p1 in range(4):                # AR order
    for q1 in range(4):            # MA order
        for p2 in range(3):        # seasonal AR order
            for q2 in range(3):    # seasonal MA order
                y_pred = []
                for i, T in enumerate(range(train_size, len(data))):
                    train_set = data.iloc[T - train_size:T]
                    model = tsa.SARIMAX(endog=train_set, # model specification
                                        order=(p1, 0, q1),
                                        seasonal_order=(p2, 0, q2, 12)).fit()
                    preds.iloc[i, 1] = model.forecast(steps=1)[0]
                mse = mean_squared_error(preds.y_true, preds.y_pred)
                results[(p1, q1, p2, q2)] = [np.sqrt(mse),
                    preds.y_true.sub(preds.y_pred).std(),
                    np.mean(aic)] 

我们还收集了 AIC 和 BIC 准则，显示出 0.94 的非常高的秩相关系数，其中 BIC 更倾向于比 AIC 稍少参数的模型。根据 RMSE 最佳的五个模型是：

 RMSE         AIC         BIC
p1 q1 p2 q2                                  
2  3  1  0   0.009323 -772.247023 -752.734581
3  2  1  0   0.009467 -768.844028 -749.331586
2  2  1  0   0.009540 -770.904835 -754.179884
   3  0  0   0.009773 -760.248885 -743.523935
   2  0  0   0.009986 -758.775827 -744.838368 

我们重新估计 SARIMA(2, 0 ,3) × (1, 0, 0) 模型，如下所示：

best_model = tsa.SARIMAX(endog=industrial_production_log_diff, order=(2, 0, 3),
                         seasonal_order=(1, 0, 0, 12)).fit()
print(best_model.summary()) 

我们得到以下总结：

图 9.5：SARMAX 模型结果

系数是显著的，并且 Q 统计量拒绝了进一步自相关的假设。相关图表同样表明我们已成功消除了序列的自相关：

图 9.6：SARIMAX 模型诊断

如何使用时间序列模型预测波动率

在金融领域中，单变量时间序列模型的一个特别重要的应用是对波动性进行预测。这是因为波动性通常随时间而变化，并且会出现波动性集聚的情况。方差的变化给使用经典的假设平稳的 ARIMA 模型进行时间序列预测带来了挑战。为了解决这一挑战，我们现在将对波动性进行建模，以便我们可以预测方差的变化。

异方差性 是指变量方差的变化的技术术语。ARCH 模型将误差项的方差表示为前期误差的函数。更具体地说，它假设误差方差遵循一个 AR(p) 模型。

广义自回归条件异方差 (GARCH) 模型扩展了 ARCH 的范围，以允许 ARMA 模型。时间序列预测通常将 ARIMA 模型用于期望均值，并将 ARCH/GARCH 模型用于时间序列的预期方差。2003 年诺贝尔经济学奖授予了罗伯特·恩格尔和克莱夫·格兰杰，因为他们开发了这一类模型。前者还在纽约大学斯特恩商学院经营着波动率实验室（vlab.stern.nyu.edu），该实验室有许多关于我们将讨论的模型的在线示例和工具。

ARCH 模型

ARCH(p) 模型简单地是应用于时间序列模型残差方差的 AR(p) 模型，这使得在时期 t 的方差条件于滞后观测的方差。

更具体地说，误差项，，是原始时间序列上的线性模型（如 ARIMA）的残差，被分为一个时间相关的标准差，，和一个扰动，z[t]，如下所示：

ARCH(p) 模型可以使用 OLS 进行估计。恩格尔提出了一种方法来使用拉格朗日乘数检验来确定适当的 ARCH 阶数，该检验对应于线性回归中所有系数为零的假设的 F 检验（见 第七章，线性模型 - 从风险因素到收益预测）。

ARCH 模型的一个关键优点是，它产生具有正超额峰度的波动性估计 —— 即，相对于正态分布的 fat tails —— 这与有关收益的实证观察相一致。缺点包括假设正面和负面波动性冲击具有相同的效应，而资产价格往往会有不同的反应。它也不能解释波动性的变化，并且可能会过度预测波动性，因为它们对收益序列的大规模、孤立的冲击反应迟缓。

对于一个正确规范的 ARCH 模型，标准化残差（除以该时期标准差的模型估计值）应该类似于白噪声，并且可以进行 Ljung-Box Q 检验。

推广 ARCH —— GARCH 模型

ARCH 模型相对简单，但通常需要许多参数来捕捉资产回报序列的波动率模式。GARCH 模型适用于对数回报序列，r[t]，具有扰动项，如图所示：，如果：

GARCH(p, q) 模型假设误差项方差的 ARMA(p, q) 模型，如图所示：。

与 ARCH 模型类似，GARCH(1,1) 过程的尾部分布比正态分布更重。该模型遇到与 ARCH 模型相同的弱点。例如，它对正面和负面冲击的反应相同。

要为 ARCH 和 GARCH 模型配置滞后阶数，请使用训练好的时间序列的平均预测的平方残差。残差被零中心化，以便它们的平方也是方差。然后，检查平方残差的 ACF 和 PACF 图，以识别时间序列方差的自相关模式。

如何构建一个预测波动率的模型

资产回报序列的波动率模型开发包括四个步骤：

1. 基于 ACF 和 PACF 显示的序列依赖性，为金融时间序列构建 ARMA 时间序列模型

2. 再次依赖于序列的 ACF 和 PACF 来测试 ARCH/GARCH 效应的模型残差

3. 如果序列相关效应显著，则指定波动率模型，并联合估计平均值和波动率方程

4. 仔细检查拟合的模型，并在必要时进行调整

当将波动率预测应用于回报序列时，序列依赖性可能有限，因此可以使用恒定均值而不是 ARMA 模型。

arch 库（参见 GitHub 上的文档链接）提供了几种选项来估计波动率预测模型。您可以将预期平均值建模为恒定值，如 如何构建自回归模型 部分中讨论的 AR(p) 模型，或者作为更近期的异质自回归过程 (HAR)，它使用每日（1 天）、每周（5 天）和每月（22 天）滞后来捕捉短期、中期和长期投资者的交易频率。

平均模型可以与几种条件异方差模型一起定义和估计，除了 ARCH 和 GARCH 之外，还包括允许在正面和负面回报之间存在非对称效应的指数 GARCH (EGARCH) 模型，以及补充 HAR 平均模型的异质 ARCH (HARCH) 模型。

我们将使用 2000 年至 2020 年的每日纳斯达克回报来演示 GARCH 模型的使用方法（详见笔记本 arch_garch_models）：

nasdaq = web.DataReader('NASDAQCOM', 'fred', '2000', '2020').squeeze()
nasdaq_returns = np.log(nasdaq).diff().dropna().mul(100) # rescale to facilitate optimization 

重新缩放后的每日收益序列仅表现出有限的自相关性，但是相对于平均值的平方偏差确实具有显著的记忆，这反映在缓慢衰减的 ACF 和 PACF 中，对于前两个而言高，在第六个滞后之后才被截断：

plot_correlogram(nasdaq_returns.sub(nasdaq_returns.mean()).pow(2), lags=120,                  title='NASDAQ Daily Volatility') 

函数plot_correlogram生成如下输出：

图 9.7：每日纳斯达克综合波动性

因此，我们可以估计一个 GARCH 模型来捕捉过去波动性的线性关系。我们将使用滚动 10 年窗口来估计一个 GARCH(p, q)模型，其中p和q的范围为 1-4，以生成一步外样本预测。

然后，我们比较预测波动性的 RMSE 与实际收益偏离其均值的平方差，以确定最具预测性的模型。我们使用修剪数据来限制极端收益值的影响，这些值在波动率的正偏态中反映得非常高：

trainsize = 10 * 252  # 10 years
data = nasdaq_returns.clip(lower=nasdaq_returns.quantile(.05),
                           upper=nasdaq_returns.quantile(.95))
T = len(nasdaq_returns)
results = {}
for p in range(1, 5):
    for q in range(1, 5):
        print(f'{p} | {q}')
        result = []
        for s, t in enumerate(range(trainsize, T-1)):
            train_set = data.iloc[s: t]
            test_set = data.iloc[t+1]  # 1-step ahead forecast
            model = arch_model(y=train_set, p=p, q=q).fit(disp='off')
            forecast = model.forecast(horizon=1)
            mu = forecast.mean.iloc[-1, 0]
            var = forecast.variance.iloc[-1, 0]
            result.append([(test_set-mu)**2, var])
        df = pd.DataFrame(result, columns=['y_true', 'y_pred'])
        results[(p, q)] = np.sqrt(mean_squared_error(df.y_true, df.y_pred)) 

GARCH(2, 2)模型实现了最低的 RMSE（与 GARCH(4, 2)相同的值，但参数更少），所以我们继续估计这个模型以检查摘要：

am = ConstantMean(nasdaq_returns.clip(lower=nasdaq_returns.quantile(.05),
                                      upper=nasdaq_returns.quantile(.95)))
am.volatility = GARCH(2, 0, 2)
am.distribution = Normal()
best_model = am.fit(update_freq=5)
print(best_model.summary()) 

输出显示了最大化的对数似然以及 AIC 和 BIC 标准，这些标准通常在选择基于样本内性能的模型时被最小化（见第七章，线性模型 - 从风险因素到收益预测）。它还显示了均值模型的结果，这种情况下，只是一个常数估计，以及常数 omega 的 GARCH 参数，AR 参数，，和 MA 参数，，所有这些都具有统计学显著性：

图 9.8：GARCH 模型结果

现在让我们探索多时间序列模型和协整概念，这将使一种新的交易策略成为可能。

多元时间序列模型

多元时间序列模型旨在同时捕获多个时间序列的动态，并利用这些序列之间的依赖关系进行更可靠的预测。对这一主题最全面的介绍是 Lütkepohl（2005）。

方程组

单变量时间序列模型，就像我们刚刚讨论的 ARMA 方法一样，仅限于目标变量及其滞后值或滞后扰动以及外生序列之间的统计关系，在 ARMAX 的情况下。相比之下，多元时间序列模型还允许其他时间序列的滞后值影响目标。这种影响适用于所有系列，导致复杂的交互作用，如下图所示：

图 9.9：单变量和多变量时间序列模型中的交互作用

除了可能更好的预测之外，多元时间序列还用于获得对交叉系列依赖性的见解。例如，在经济学中，多元时间序列用于理解一个变量的政策变化，例如利率，如何在不同的视角影响其他变量。

多变量模型产生的冲击响应函数达到了这一目的，并允许我们模拟一个变量如何对其他变量的突然变化做出反应。格兰杰因果性概念分析了一个变量是否对另一个变量的预测有用（以最小二乘意义上）。此外，多元时间序列模型允许对预测误差方差进行分解，以分析其他系列如何做出贡献。

向量自回归（VAR）模型

我们将看到向量自回归 VAR(p) 模型如何通过创建一个包含所有k系列p滞后值的k方程组来扩展 AR(p) 模型。在最简单的情况下，k=2 的 VAR(1) 模型采用以下形式：

这个模型可以用矩阵形式更简洁地表达：

对输出滞后值的系数提供有关系列本身动态的信息，而交叉变量的系数则提供一些关于系列之间交互作用的见解。此符号扩展到k时间序列和p阶，如下所示：

VAR(p)模型还需要平稳性，以便从单变量时间序列建模的初步步骤延伸。首先，探索系列并确定必要的转换。然后，应用增广迪基-福勒检验来验证每个系列是否满足平稳性标准，否则应用进一步的转换。它可以根据初始信息进行 OLS 估计，也可以根据 MLE 进行估计，这是对正态分布误差而言的等价物，但对其他情况不是。

如果一些或全部的k系列是单位根非平稳的，它们可能是协整的（见下一节）。将单位根概念扩展到多个时间序列意味着两个或更多系列的线性组合是平稳的，因此是均值回归的。

VAR 模型无法处理这种情况，需要进行差分处理；而是使用向量误差修正模型（VECM，Johansen 和 Juselius 1990）。我们将进一步探讨协整性，因为如果存在并被认为是持续存在的，它可以用于配对交易策略。

滞后阶数的确定也从每个系列的 ACF 和 PACF 中获得线索，但受到相同滞后阶数适用于所有系列的限制。在模型估计之后，残差诊断还要求结果类似于白噪声，并且模型选择可以使用样本内信息准则，或者如果目标是使用模型进行预测，则可以使用样本外预测性能来交叉验证备选模型设计

如在单变量情况中提到的那样，对原始时间序列的预测要求我们在训练模型之前撤销应用于使序列稳定的转换。

使用 VAR 模型进行宏观预测

我们将扩展使用单一时间序列的示例，其中包括工业生产的月度数据以及消费者情绪的月度时间序列，这两者都由美联储的数据服务提供。我们将使用熟悉的 pandas-datareader 库从 1970 年到 2017 年检索数据：

df = web.DataReader(['UMCSENT', 'IPGMFN'],
                     'fred', '1970', '2017-12').dropna()
df.columns = ['sentiment', 'ip'] 

对工业生产系列进行对数变换，并使用 12 个滞后期进行季节性差分，得到稳定的结果：

df_transformed = pd.DataFrame({'ip': np.log(df.ip).diff(12),
                              'sentiment': df.sentiment.diff(12)}).dropna()
test_unit_root(df_transformed) # see notebook for details and additional plots
          p-value
ip          0.0003
sentiment   0.0000 

这给我们留下了以下系列：

图 9.10：转换后的时间序列：工业生产和消费者情绪

为了限制输出大小，我们将仅使用 statsmodels 的VARMAX实现（允许使用可选的外生变量）以及使用前 480 个观测来估计一个带有常数趋势的 VAR(1)模型：

model = VARMAX(df_transformed.loc[:'2017'], order=(1,1),
               trend='c').fit(maxiter=1000) 

这产生了以下摘要：

图 9.11：VAR(1)模型结果

输出包含了两个时间序列方程的系数，如前述 VAR(1)示例所述。statsmodels 提供了诊断图来检查残差是否符合白噪声假设。在这个简单的例子中，这并不完全成立，因为方差似乎不是恒定的（左上角），而且量化图显示了分布的差异，即尾部较大（左下角）：

图 9.12：statsmodels VAR 模型诊断图

您可以按以下方式生成样本外预测：

preds = model.predict(start=480, end=len(df_transformed)-1) 

下面的实际值和预测值的可视化展示了预测滞后于实际值，并且不能很好地捕捉非线性、样本外模式：

图 9.13：VAR 模型预测与实际值

协整 - 具有共享趋势的时间序列

我们在前一节关于多变量时间序列模型中简要提到了协整。现在让我们更详细地解释这个概念以及如何诊断其存在，然后再利用它进行统计套利交易策略。

我们已经看到一个时间序列可以具有创建随机趋势的单位根，并使时间序列高度持久。当我们在线性回归模型中将这种整合的时间序列作为特征使用时，而不是作为差分形式使用，它与结果的关系通常会显得统计上显著，尽管实际上并不是。这种现象称为伪回归（详见 Wooldridge，2008 年的第十八章，用于金融时间序列和卫星图像的 CNNs）。因此，建议的解决方案是在将它们用于模型之前对时间序列进行差分，使它们变得平稳。

然而，在结果和一个或多个输入变量之间存在协整关系时，有一个例外。要理解协整的概念，让我们首先记住回归模型的残差是输入和输出系列的线性组合。

通常，对一个或多个这样的系列的整合时间序列的回归的残差产生的是非平稳的残差，这些残差也是整合的，因此表现出类似于随机游走的行为。然而，对于一些时间序列，情况并非如此：回归产生的系数会形成时间序列的线性组合，这种线性组合是平稳的，尽管各个系列不是。这种时间序列是协整的。

一个非技术示例是一个醉酒的人在随机行走，他的狗（被拴在绳索上）陪伴着他。两个轨迹都是非平稳的，但协整，因为狗偶尔会回到它的主人身边。在交易的背景下，套利约束意味着现货和期货价格之间存在协整关系。

换句话说，两个或多个协整系列的线性组合有一个稳定的均值，这个线性组合会回归到这个稳定的均值。当各个系列是高阶整合的时候，这也适用，线性组合会减少整体的整合次数。

协整与相关性不同：两个系列可以高度相关，但不一定是协整的。例如，如果两个增长系列是彼此的常数倍数，它们的相关性将很高，但任何线性组合也会增长而不是回归到一个稳定的均值。

协整非常有用：如果两个或更多资产价格系列趋向于回归到一个共同的均值，我们可以利用与趋势偏离，因为它们应该暗示着未来价格的相反方向的移动。协整背后的数学更加复杂，因此我们将只关注实践方面；有关深入处理，请参阅 Lütkepohl（2005）。

在本节中，我们将讨论如何识别具有这种长期平稳关系的配对，估计任何失衡纠正的预期时间，以及如何利用这些工具来实施和回测长短对交易策略。

有两种测试协整的方法：

• 恩格尔-格兰杰两步法

• 约翰逊测试

在展示它们如何帮助识别趋向于恢复到共同趋势的共整合证券之前，我们将依次讨论每一个。

恩格尔-格兰杰两步法

恩格尔-格兰杰方法用于确定两个系列之间的共整合关系。它涉及到以下两个方面：

1. 对另一个系列进行回归以估计长期稳定关系

2. 将 ADF 单位根测试应用于回归残差

零假设是残差具有单位根并且是集成的；如果我们可以拒绝它，那么我们假设残差是稳定的，因此系列是共整合的（Engle and Granger 1987）。

这种方法的一个关键好处是回归系数表示的是使组合稳定的乘数，即，均值回归。不幸的是，测试结果将取决于我们考虑哪个变量是独立的，因此我们尝试两种方式，然后选择具有更负的测试统计量且具有较低 p 值的关系。

这项测试的另一个缺点是，它仅限于成对关系。更复杂的约翰逊程序可以识别多达十几个时间序列之间的显著共整合。

约翰逊似然比检验

相比之下，约翰逊程序测试了上一节中讨论的共整合对 VAR 模型施加的限制。更具体地说，在从通用 VAR(p)模型的两边减去目标向量之后，我们得到了误差修正模型（ECM）的表述：

结果修改后的 VAR(p)方程只有一个向量项在水平（y[t][-1]）中，并且没有使用运算符作为差异。共整合的性质取决于系数矩阵的等级（约翰逊 1991）。

尽管这个方程在结构上与 ADF 测试设置相似，但现在有了多个系列的潜在公共趋势组合。为了确定共整合关系的数量，约翰逊测试连续测试增加的等级，从 0 开始（无共整合）。我们将在下一节探讨应用于两个系列的情况。

冈萨洛和李（1998）讨论了由于错误指定的模型动态和其他实施方面的挑战，包括如何结合我们将在下一节中依赖的两种测试程序的实际统计套利策略。

具有共整合关系的统计套利

统计套利指的是利用某些统计模型或方法来利用资产的相对定价错误，同时保持一定程度的市场中性的策略。

对冲交易 是一个概念上简单直接的策略，至少从 20 世纪八十年代中期起就被算法交易者采用（Gatev、Goetzmann 和 Rouwenhorst 2006）。其目标是找到历史上价格走势相似的两个资产，跟踪价差（即它们价格之间的差异），并一旦价差扩大，买入价格低于共同趋势的失败者并做空价格高于共同趋势的赢家。如果这种关系持续存在，做多和/或做空腿将随着价格趋同而获利，并且位置将被关闭。

该方法通过从多个证券形成篮子，并使一个资产对抗另一个篮子而扩展到多变量环境。

在实践中，该策略需要两个步骤：

1. 形成阶段：确定具有长期均值回归关系的证券。理想情况下，价差应具有较高的方差，以允许频繁的盈利交易，并可靠地回归到共同趋势。

2. 交易阶段：当价格变动导致价差分歧和趋同时，触发进出交易规则。

在这一领域越来越活跃的研究中，过去几年中出现了几种形成和交易阶段的方法，跨越多个资产类别。下一小节概述了主要区别，然后我们将深入介绍一个示例应用。

如何选择和交易共动资产对

最近对对冲交易策略的全面调查（Krauss 2017）确定了四种不同的方法，以及一些更近期的方法，包括基于 ML 的预测：

• 距离法：最古老且研究最多的方法使用诸如相关性之类的距离度量标识候选对，并使用诸如布林带之类的非参数阈值触发进出交易。其计算简单性允许进行大规模应用，自从 Gatev 等人（2006）以来，在不同市场和资产类别中已经表现出盈利性质，并持续了相当长的时间。然而，最近的表现已经有所衰退。

• 协整法：如前所述，该方法依赖于两个或更多变量之间的长期关系的计量模型，并允许进行统计测试，承诺比简单距离度量更可靠。该类别中的示例使用 Engle-Granger 和 Johansen 程序来识别证券对和篮子，以及旨在捕获概念的更简单的启发式方法（Vidyamurthy 2004）。交易规则通常类似于距离度量使用的简单阈值。

• 时间序列法：着眼于交易阶段，该类别中的策略旨在将价差建模为均值回归的随机过程，并相应地优化进出规则（Elliott、Hoek 和 Malcolm 2005）。它假定有希望的配对已经被确定。

• 随机控制方法：类似于时间序列方法，其目标是使用随机控制理论优化交易规则，以找到值函数和策略函数来得到最优组合（Liu 和 Timmermann，2013）。我们将在第二十一章，用于合成时间序列数据的生成对抗网络中讨论这种类型的方法。

• 其他方法：除了基于无监督学习（如主成分分析，参见第十三章，使用无监督学习进行数据驱动的风险因子和资产配置）和统计模型（如 copulas，Patton 2012）的配对识别外，近年来机器学习也变得流行，用于基于相对价格或收益预测识别配对（Huck，2019）。我们将在接下来的章节中介绍几种可用于此目的的 ML 算法，并且说明相应的多变量配对交易策略

这些不同方法的摘要只是略显一斑地展示了配对交易策略设计所带来的灵活性。除了有关配对选择和交易规则逻辑的高层次问题之外，还有许多参数需要我们定义以供实施。这些参数包括以下内容：

• 用于筛选潜在配对或篮子的投资范围

• 形成期的长度

• 用于选择可交易候选者的关系强度

• 与公共均值的偏离程度和收敛程度来触发进入或退出交易或根据价差波动调整现有仓位

实践中的配对交易

距离方法使用（标准化的）资产价格或其收益的相关性来识别配对，简单且计算成本远远低于协整检验。笔记本cointegration_test对具有 4 年每日数据的约 150 只股票样本进行了说明：计算与 ETF 收益的相关性大约需要 30ms，而进行一系列协整检验（使用 statsmodels）则需要 18 秒 - 慢了 600 倍。

速度优势尤为宝贵。这是因为潜在配对的数量是要考虑到每一方的候选人数的乘积，因此评估 100 只股票和 100 只 ETF 的组合需要比较 10,000 个测试（我们稍后将讨论多重测试偏差的挑战）。

另一方面，距离度量不一定选择最有利可图的配对：相关性在完美共同运动时最大化，这反过来消除了实际的交易机会。经验证实证研究表明，协整配对的价差波动性几乎是距离配对的两倍（Huck 和 Afawubo，2015）。

为了平衡计算成本与生成对质量之间的权衡，Krauss（2017）根据他的文献综述建议采用一种结合了两种方法的程序：

1. 选择传播稳定且漂移小的对，以减少候选项的数量

2. 对剩余的传播方差最高的对进行协整检验

此过程旨在选择具有较低偏离风险的协整对，同时确保更具波动性的传播，进而产生更多的利润机会。

大量的测试引入了数据窥探偏差，如第六章，机器学习过程中所讨论的：多重检验可能会增加错误地拒绝无协整假设的假阳性数量。虽然统计显著性对于盈利交易可能不是必需的（Chan 2008），但商品对研究（Cummins and Bucca 2012）表明，控制家族内误差率以提高检验的功效，根据 Romano and Wolf（2010）的说法，可以带来更好的性能。

在接下来的小节中，我们将更详细地研究各种资产价格共同变动程度的预测能力对协整检验结果的影响。

示例代码使用了 172 只股票和 138 只 ETFs，在 2010 年至 2019 年间由 Stooq 提供的每日数据，这些股票和 ETFs 在纽约证券交易所和纳斯达克交易。

这些证券代表了其所属类别在样本期内的最大平均美元成交量；高度相关且平稳的资产已被移除。有关如何获取数据的说明，请参见 GitHub 仓库的data文件夹中的笔记本create_datasets，以及有关代码和额外预处理和探索细节的笔记本cointegration_tests。

基于距离的启发式方法来寻找协整对

compute_pair_metrics() 计算了超过 23,000 对股票和交易所交易基金（ETFs）在 2010-14 年和 2015-19 年的以下距离度量：

• 传播的漂移，定义为传播的时间趋势对传播的线性回归

• 传播的波动性

• 标准化价格系列之间和它们的回报之间的相关性

低漂移和波动，以及高相关性，是协整的简单代理。

为了评估这些启发式规则的预测能力，我们还使用 statsmodels 对前述对运行恩格尔-格兰杰和约翰逊协整检验。这发生在compute_pair_metrics()的后半部分的循环中。

我们首先估计我们需要为约翰逊检验指定的滞后数的最佳数量。对于两个测试，我们假设协整系列（传播）可能有一个不同于零的截距，但没有趋势：

def compute_pair_metrics(security, candidates):
    security = security.div(security.iloc[0])
    ticker = security.name
    candidates = candidates.div(candidates.iloc[0])
    # compute heuristics
    spreads = candidates.sub(security, axis=0)
    n, m = spreads.shape
    X = np.ones(shape=(n, 2))
    X[:, 1] = np.arange(1, n + 1)
    drift = ((np.linalg.inv(X.T @ X) @ X.T @ spreads).iloc[1]
             .to_frame('drift'))
    vol = spreads.std().to_frame('vol')
    corr_ret = (candidates.pct_change()
                .corrwith(security.pct_change())
                .to_frame('corr_ret'))
    corr = candidates.corrwith(security).to_frame('corr')
    metrics = drift.join(vol).join(corr).join(corr_ret).assign(n=n)
    tests = []
    # compute cointegration tests
    for candidate, prices in candidates.items():
        df = pd.DataFrame({'s1': security, 's2': prices})
        var = VAR(df)
        lags = var.select_order() # select VAR order
        k_ar_diff = lags.selected_orders['aic']
        # Johansen Test with constant Term and estd. lag order
        cj0 = coint_johansen(df, det_order=0, k_ar_diff=k_ar_diff)
        # Engle-Granger Tests
        t1, p1 = coint(security, prices, trend='c')[:2]
        t2, p2 = coint(prices, security, trend='c')[:2]
        tests.append([ticker, candidate, t1, p1, t2, p2, 
                      k_ar_diff, *cj0.lr1])

    return metrics.join(tests) 

为了检验协整检验的显著性，我们将约翰逊迹统计量对于秩 0 和 1 的临界值进行比较，并得到恩格尔-格兰杰 p 值。

我们遵循上一节末尾提到的 Gonzalo 和 Lee (1998) 的建议，即同时应用两个测试，并接受它们达成一致意见的对。作者建议在存在分歧时进行额外的尽职调查，但我们将跳过这一步：

spreads['trace_sig'] = ((spreads.trace0 > trace0_cv) &
                        (spreads.trace1 > trace1_cv)).astype(int)
spreads['eg_sig'] = (spreads.p < .05).astype(int) 

对于两个样本期间的超过 46,000 对，约翰逊测试将 3.2% 的关系视为显著，而恩格尔-格兰杰测试将 6.5% 视为显著。它们对 366 对（0.79%）达成一致意见。

启发式方法能多好地预测显著的协整？

当我们比较那些根据两个测试协整的系列和其余系列的启发式方法的分布时，波动性和漂移确实较低（绝对值）。图 9.14 显示这两个相关性指标的情况不太明确：

图 9.14：启发式方法的分布，按照两个协整测试的显著性来分解

为了评估启发式方法的预测准确性，我们首先运行一个 logistic 回归模型，使用这些特征来预测显著的协整。它达到了曲线下面积 (AUC) 交叉验证分数为 0.815；排除相关性指标后，它仍然得分 0.804。决策树在 AUC=0.821 时表现稍好，无论是否包含相关性特征。

由于强烈的类别不平衡，存在大量的假阳性：正确识别 366 个协整对的 80% 意味着有超过 16,500 个假阳性，但也消除了几乎 30,000 个候选对。有关更多详细信息，请参阅笔记本 cointegration_tests。

关键要点 是距离启发式方法可以帮助更高效地筛选大范围的内容，但这样做的代价是可能会错过一些共整对，并且仍然需要大量的测试。

准备策略回测

在本节中，我们将基于股票和 ETF 的样本以及 2017-2019 年期间实施基于协整的统计套利策略。为了简化演示，某些方面进行了简化。有关代码示例和额外细节，请参阅笔记本 statistical_arbitrage_with_cointegrated_pairs。

我们首先生成并存储所有候选对的协整测试及其产生的交易信号，然后，鉴于该过程的计算强度，我们对基于这些信号的策略进行回测。

预计算协整测试

首先，我们在一个两年的回溯期内对每一对可能的 23,000 个对进行季度协整测试，然后，我们选择那些约翰逊（Johansen）和恩格尔-格兰杰（Engle-Granger）测试均同意进行交易的对。我们应该在回溯期内排除静止的资产，但我们排除了整个期间都是静止的资产，所以我们跳过了这一步以简化流程。

此过程遵循先前概述的步骤；请参阅笔记本以获取详细信息。

图 9.15显示了选定用于交易的两个不同配对的原始股票和 ETF 系列；请注意样本期间共同趋势的明显存在：

图 9.15：样本期间两个选定配对的价格系列

获取入场和出场交易

现在，我们可以根据滚动对冲比率计算每个候选配对的价差。我们还计算布林带，因为我们将考虑价差大于其移动平均值两个滚动标准偏差的移动，作为长和短的入场信号，并将移动平均值的交叉点反向作为退出信号。

使用卡尔曼滤波器平滑价格

为此，我们首先应用滚动卡尔曼滤波器（KF）来消除一些噪音，正如在第四章，金融特征工程 - 如何研究 Alpha 因子中所示：

def KFSmoother(prices):
    """Estimate rolling mean"""

    kf = KalmanFilter(transition_matrices=np.eye(1),
                      observation_matrices=np.eye(1),
                      initial_state_mean=0,
                      initial_state_covariance=1,
                      observation_covariance=1,
                      transition_covariance=.05)
    state_means, _ = kf.filter(prices.values)
    return pd.Series(state_means.flatten(),
                     index=prices.index) 

使用卡尔曼滤波器计算滚动对冲比率

为了获得动态对冲比率，我们使用 KF 进行滚动线性回归，如下所示：

def KFHedgeRatio(x, y):
    """Estimate Hedge Ratio"""
    delta = 1e-3
    trans_cov = delta / (1 - delta) * np.eye(2)
    obs_mat = np.expand_dims(np.vstack([[x], [np.ones(len(x))]]).T, axis=1)
    kf = KalmanFilter(n_dim_obs=1, n_dim_state=2,
                      initial_state_mean=[0, 0],
                      initial_state_covariance=np.ones((2, 2)),
                      transition_matrices=np.eye(2),
                      observation_matrices=obs_mat,
                      observation_covariance=2,
                      transition_covariance=trans_cov)
    state_means, _ = kf.filter(y.values)
    return -state_means 

估计均值回归的半衰期

如果我们将价差视为一个连续时间的均值回归随机过程，我们可以将其建模为奥恩斯坦-乌伦贝克过程。这种视角的好处在于，我们得到了一个均值回归半衰期的公式，作为偏离后价差再次收敛所需时间的近似值（有关详细信息，请参阅 Chan 2013 年的第二章，市场和基本数据 - 来源和技术）：

def estimate_half_life(spread):
    X = spread.shift().iloc[1:].to_frame().assign(const=1)
    y = spread.diff().iloc[1:]
    beta = (np.linalg.inv(X.T@X)@X.T@y).iloc[0]
    halflife = int(round(-np.log(2) / beta, 0))
    return max(halflife, 1) 

计算价差和布林带

以下函数组织了前述计算，并将价差表达为 z 分数，该分数捕捉了相对于滚动标准偏差的移动平均值的偏差，窗口大小等于两个半衰期：

def get_spread(candidates, prices):
    pairs, half_lives = [], []
    periods = pd.DatetimeIndex(sorted(candidates.test_end.unique()))
    start = time()
    for p, test_end in enumerate(periods, 1):
        start_iteration = time()
        period_candidates = candidates.loc[candidates.test_end == test_end, 
                                          ['y', 'x']]
        trading_start = test_end + pd.DateOffset(days=1)
        t = trading_start - pd.DateOffset(years=2)
        T = trading_start + pd.DateOffset(months=6) - pd.DateOffset(days=1)
        max_window = len(prices.loc[t: test_end].index)
        print(test_end.date(), len(period_candidates))
        for i, (y, x) in enumerate(zip(period_candidates.y, 
                                       period_candidates.x), 1):
            pair = prices.loc[t: T, [y, x]]
            pair['hedge_ratio'] = KFHedgeRatio(
                y=KFSmoother(prices.loc[t: T, y]),
                x=KFSmoother(prices.loc[t: T, x]))[:, 0]
            pair['spread'] = pair[y].add(pair[x].mul(pair.hedge_ratio))
            half_life = estimate_half_life(pair.spread.loc[t: test_end])
            spread = pair.spread.rolling(window=min(2 * half_life, 
                                                    max_window))
            pair['z_score'] = pair.spread.sub(spread.mean()).div(spread.
std())
            pairs.append(pair.loc[trading_start: T].assign(s1=y, s2=x, period=p, pair=i).drop([x, y], axis=1))
            half_lives.append([test_end, y, x, half_life])
    return pairs, half_lives 

获取长和短头寸的入场和出场日期

最后，我们使用一组 z 分数来推导交易信号：

1. 如果 z 分数低于（高于）两个，我们进入长（空）头寸，这意味着价差已经移动了两个滚动标准偏差低于（高于）移动平均值

2. 当价差再次穿过移动平均线时，我们退出交易

我们每季度制定一组规则，用于通过协整测试的一组配对，这些测试是在先前的回溯期间进行的，但允许在随后的 3 个月内退出配对。

我们再次简化了这个过程，通过删除在这个 6 个月期间没有收盘的配对。或者，我们可以通过我们在策略中包含的止损风险管理来处理这个问题（请参阅关于回测的下一节）：

def get_trades(data):
    pair_trades = []
    for i, ((period, s1, s2), pair) in enumerate(
             data.groupby(['period', 's1', 's2']), 1):
        if i % 100 == 0:
            print(i)
        first3m = pair.first('3M').index
        last3m = pair.last('3M').index
        entry = pair.z_score.abs() > 2
        entry = ((entry.shift() != entry)
                 .mul(np.sign(pair.z_score))
                 .fillna(0)
                 .astype(int)
                 .sub(2))
        exit = (np.sign(pair.z_score.shift().fillna(method='bfill'))
                != np.sign(pair.z_score)).astype(int) - 1
        trades = (entry[entry != -2].append(exit[exit == 0])
                  .to_frame('side')
                  .sort_values(['date', 'side'])
                  .squeeze())
        trades.loc[trades < 0] += 2
        trades = trades[trades.abs().shift() != trades.abs()]
        window = trades.loc[first3m.min():first3m.max()]
        extra = trades.loc[last3m.min():last3m.max()]
        n = len(trades)
        if window.iloc[0] == 0:
            if n > 1:
                print('shift')
                window = window.iloc[1:]
        if window.iloc[-1] != 0:
            extra_exits = extra[extra == 0].head(1)
            if extra_exits.empty:
                continue
            else:
                window = window.append(extra_exits)
        trades = (pair[['s1', 's2', 'hedge_ratio', 'period', 'pair']]
                  .join(window. to_frame('side'), how='right'))
        trades.loc[trades.side == 0, 'hedge_ratio'] = np.nan
        trades.hedge_ratio = trades.hedge_ratio.ffill()
        pair_trades.append(trades)
    return pair_trades 

使用 backtrader 进行策略回测

现在，我们准备在我们的回测平台上制定我们的策略，执行它，并评估结果。为此，除了跟踪我们的配对以外，我们还需要跟踪单独的投资组合头寸，并监视活跃和非活跃配对的价差，以应用我们的交易规则。

使用自定义 DataClass 跟踪配对

为了考虑活跃的配对，我们定义了一个dataclass（在 Python 3.7 中引入 - 详见 Python 文档以获取详情）。这个数据结构，称为Pair，允许我们存储配对组件、它们的股数和对冲比率，并计算当前价差和收益等内容。请参见以下代码中的简化版本：

@dataclass
class Pair:
    period: int
    s1: str
    s2: str
    size1: float
    size2: float
    long: bool
    hr: float
    p1: float
    p2: float
    entry_date: date = None
    exit_date: date = None
    entry_spread: float = np.nan
    exit_spread: float = np.nan
    def compute_spread(self, p1, p2):
        return p1 * self.size1 + p2 * self.size2
    def compute_spread_return(self, p1, p2):
        current_spread = self.compute_spread(p1, p2)
        delta = self.entry_spread - current_spread
        return (delta / (np.sign(self.entry_spread) *
                         self.entry_spread)) 

运行和评估策略

关键的实施方面包括：

• 每天退出已触发退出规则或超过给定负回报的配对

• 对触发入场信号的价差的新多头和空头仓位的开仓

• 此外，我们调整头寸以考虑不同数量的配对

策略本身的代码在这里展示太占用空间；详见笔记本pairs_trading_backtest获取详情。

图 9.16显示，至少在 2017-2019 年期间，这个简化的策略有其时机（请注意，我们利用了一些前瞻性偏见并忽略了交易成本）。

在这些宽松的假设下，它在期间的开始和结束时表现不及标普 500 指数，并且在其他时候大致处于相同水平（左侧面板）。它产生 0.08 的 alpha 和-0.14 的负 beta（右侧面板），平均夏普比率为 0.75，Sortino 比率为 1.05（中央面板）：

图 9.16：策略绩效指标

尽管我们应该对这些绩效指标持谨慎态度，但这个策略展示了一种基于协整的统计套利的解剖学，以配对交易的形式呈现。让我们看看您可以采取哪些步骤来在此框架上进行改进以获得更好的表现。

扩展 - 如何做得更好

协整是一个非常有用的概念，可以识别出倾向于同步运动的股票配对或组合。与协整的统计复杂性相比，我们使用的是非常简单和静态的交易规则；季度基础上的计算也扭曲了策略，因为长期和短期持有的模式显示（请参见笔记本）。

要取得成功，你至少需要筛选更大的资产组合，并优化其中的几个参数，包括交易规则。此外，风险管理应考虑到当某些资产经常出现在同一交易对的同一侧时产生的集中仓位。

你也可以操作篮子而不是单个配对；然而，为了解决不断增长的候选数量，你可能需要限制篮子的组成。

正如在配对交易 - 具有协整的统计套利部分中所述，存在着旨在预测价格走势的替代方案。在接下来的章节中，我们将探讨各种机器学习模型，这些模型旨在预测给定投资范围和时间跨度的价格走势的绝对大小或方向。将这些预测用作多头和空头进入信号是对我们在本节中学习的配对交易框架的自然扩展或替代方法。

总结

在本章中，我们探讨了用于单个系列的单变量线性时间序列模型，以及用于多个交互系列的多变量模型。我们遇到了用于预测宏观基本面的应用程序，用于风险管理中广泛使用的预测资产或组合波动性的模型，以及捕捉多个宏观系列动态的多元 VAR 模型。我们还研究了协整概念，这是支撑流行的配对交易策略的基础。

类似于第七章，线性模型 - 从风险因素到收益预测，我们看到线性模型施加了很多结构，也就是说，它们做出了强烈的假设，可能需要进行转换和广泛的测试来验证这些假设是否成立。如果确实如此，模型的训练和解释就会很直接，而且这些模型提供了一个很好的基准，更复杂的模型可能能够改进。在接下来的两章中，我们将看到两个示例，即随机森林和梯度提升模型，并且我们将在第四部分中遇到更多示例，该部分是关于深度学习的。

第十章：贝叶斯 ML – 动态夏普比率和对冲交易

在本章中，我们将介绍贝叶斯方法在机器学习（ML）中的应用以及它们在开发和评估交易策略时对不确定性的不同视角所带来的价值。

贝叶斯统计学允许我们量化对未来事件的不确定性，并以原则性的方式在新信息到达时改进我们的估计。这种动态方法在金融市场的演变性质下适应良好。当相关数据较少且我们需要系统地整合先前知识或假设时，它特别有用。

我们将看到，贝叶斯方法对机器学习的应用允许更丰富地了解统计指标、参数估计和预测周围的不确定性。这些应用范围从更精细的风险管理到动态更新预测模型，其中包括市场环境的变化。资产配置的黑-利特曼方法（参见第五章，组合优化和绩效评估）可以被解释为贝叶斯模型。它计算资产的预期回报作为市场均衡和投资者观点的加权平均值，每个资产的波动性，跨资产相关性和对每个预测的信心。

更具体地说，在本章中，我们将涵盖：

• 贝叶斯统计学如何应用于 ML

• 使用 PyMC3 进行概率编程

• 使用 PyMC3 定义和训练 ML 模型

• 如何运行最先进的抽样方法进行近似推理

• 贝叶斯 ML 应用于计算动态夏普比率、动态对冲交易对冲比率和估计随机波动性

您可以在 GitHub 存储库的相应目录中找到本章的代码示例和附加资源的链接。笔记本包括图像的彩色版本。

贝叶斯机器学习的工作原理

经典统计学被认为遵循频率主义方法，因为它将概率解释为长期观察大量试验后事件发生的相对频率。在概率的背景下，事件是一个或多个试验的基本结果的组合，例如掷两个骰子的六个等可能结果中的任何一个，或者在某一天资产价格下跌了 10%或更多。

相比之下，贝叶斯统计学将概率视为事件发生的信心或信念的度量。因此，贝叶斯观点为主观看法和意见差留下了更多的空间，而不是频率主义解释。这种差异在那些不经常发生以至于无法得出长期频率客观度量的事件中最为显著。

换句话说，频率统计假设数据是从总体中随机抽样的，并旨在识别生成数据的固定参数。贝叶斯统计则将数据视为给定，并认为参数是随数据可推断的随机变量的分布。因此，频率主义方法至少需要与要估计的参数数量相同数量的数据点。另一方面，贝叶斯方法适用于较小的数据集，并且非常适合从一次一个样本进行在线学习。

贝叶斯观点对于许多在某些重要方面罕见或独特的现实事件非常有用。例如，下次选举的结果或者市场是否会在 3 个月内崩溃的问题。在每种情况下，都有相关的历史数据以及随着事件逼近而展开的独特环境。

我们首先会介绍贝叶斯定理，该定理通过将先验假设与新的经验证据相结合来阐明通过更新信念的概念，并将结果参数估计与其频率论对应物进行比较。然后我们将展示两种贝叶斯统计推断方法，即共轭先验和近似推断，这些方法揭示了潜在参数（即未观察到的参数）的后验分布，如期望值：

• 共轭先验 通过提供闭合形式的解决方案来促进更新过程，从而使我们能够精确计算解决方案。然而，这种精确的分析方法并不总是可用。

• 近似推断 模拟由假设和数据结合而成的分布，并使用该分布的样本来计算统计洞见。

如何从经验证据中更新假设

"当事实改变时，我改变我的想法。您怎么做，先生？"

–约翰·梅纳德·凯恩斯

250 多年前，托马斯·贝叶斯牧师提出的定理使用基本的概率理论来规定概率或信念应该如何随着相关新信息的到来而改变。前面的凯恩斯引言捕捉了那种精神。它依赖于条件和总概率以及链式法则；有关介绍和更多信息，请参见 Bishop（2006）和 Gelman 等人（2013）。

概率信念涉及单个参数或参数向量（也称为假设）。每个参数可以是离散的或连续的。可以是一维统计量，例如（离散的）分类变量的众数或（连续的）均值，或者是更高维的值集，例如协方差矩阵或深度神经网络的权重。

与频率统计学的关键区别在于，贝叶斯假设是以概率分布而不是参数值表示的。因此，频率派推断侧重于点估计，而贝叶斯推断产生概率分布。

贝叶斯定理通过计算从以下输入计算后验概率分布来更新对感兴趣参数的信念，如图 10.1 所示：

• 先验分布表示我们考虑每个可能假设的可能性有多大。

• 似然函数输出在给定参数 值时观察到数据集的概率，即针对特定假设。

• 证据衡量了给定所有可能假设的观察数据的可能性。因此，对于所有参数值，它是相同的，并用于规范化分子。

图 10.1：证据如何更新先验到后验概率分布

后验是先验和似然的乘积，除以证据。因此，它反映了假设的概率分布，通过考虑先前假设和数据更新。从不同角度看，后验概率来自应用链式规则，这反过来将数据和参数的联合分布分解为因子。

对于更高维度、连续变量，该公式变得更加复杂，涉及（多个）积分。此外，另一种表述使用赔率将后验赔率表达为先验赔率乘以似然比（参见 Gelman 等人 2013 年）。

精确推断 - 最大后验估计

将贝叶斯定理应用于精确计算后验概率的实际应用非常有限。这是因为计算分母中的证据项相当具有挑战性。证据反映了观察数据在所有可能参数值上的概率。它也被称为边际似然，因为它需要通过添加或对其分布进行积分来“边缘化”参数的分布。这通常仅在具有少量假设值的少量离散参数的简单情况下才可能。

最大后验概率（MAP）估计利用了证据是一个将后验缩放到满足概率分布要求的常数因子的事实。由于证据不依赖于 ，后验分布与似然和先验的乘积成正比。因此，MAP 估计选择使后验最大化的值，给定观察到的数据和先验信念，即后验的模式。

MAP 方法与定义概率分布的参数的最大似然估计（MLE）相对比。MLE 选择最大化观察训练数据的似然函数的参数值！。

通过定义可以看出，MAP 与 MLE 的不同之处在于包括先验分布。换句话说，除非先验是一个常数，否则 MAP 估计将与其 MLE 对应物不同：

MLE 解决方案往往反映了频率学观念，即概率估计应反映观察到的比率。另一方面，先验对 MAP 估计的影响通常对应于将反映先验假设的数据添加到 MLE 中。例如，一个关于硬币偏向的强先验可以通过添加倾斜的试验数据来纳入 MLE 环境中。

先验分布是贝叶斯模型的关键组成部分。我们现在将介绍一些便利的选择，以促进分析推断。

如何选择先验

先验应该反映对参数分布的了解，因为它影响着 MAP 估计。如果先验不确定，我们需要做出选择，通常从几个合理的选项中选择。一般来说，证明先验并通过测试替代是否导致相同结论是一种良好的做法。

有几种类型的先验：

• 客观先验将最大化数据对后验的影响。如果参数分布未知，我们可以选择一个无信息的先验，如参数值的相关范围上的均匀分布，也称为平坦先验。

• 相比之下，主观先验旨在将模型外的信息纳入估计中。在 Black-Litterman 环境中，投资者对资产未来收益的信念将是主观先验的一个例子。

• 一个经验先验结合了贝叶斯和频率主义方法，并利用历史数据消除了主观性，例如，通过估计各种时刻来拟合标准分布。使用一些历史上的日收益率平均值而不是对未来收益的信念，就是一个简单经验先验的例子。

在 ML 模型的上下文中，先验可以被视为正则化器，因为它限制了后验可以假设的值。例如，具有零先验概率的参数不是后验分布的一部分。通常，更多的好数据允许得出更强的结论，并减少先验的影响。

如何保持推断简单 - 共轭先验

当得到的后验与先验相同类或家族的分布相同时，先验分布与似然函数共轭。例如，当先验和似然都是正态分布时，后验也是正态分布的。

先验和似然的共轭意味着后验的封闭形式解，有利于更新过程，并避免使用数值方法来近似后验。此外，得到的后验可以用作下一次更新步骤的先验。

让我们用股票价格变动的二元分类例子来说明这个过程。

资产价格变动的动态概率估计

当数据由具有某一成功概率的二元伯努利随机变量组成时，重复试验中的成功次数遵循二项分布。共轭先验是支持区间[0,1]上的贝塔分布，并具有两个形状参数，以对成功概率上的任意先验分布进行建模。因此，后验分布也是一个贝塔分布，我们可以通过直接更新参数来得到。

我们将收集不同大小的二值化的每日标普 500 指数收益率样本，其中正面结果是价格上涨。从一个不具信息的先验开始，在区间[0,1]中为每个可能的成功概率分配相等的概率，我们计算不同证据样本的后验概率。

下面的代码示例显示，更新仅涉及将观察到的成功和失败的数量添加到先验分布的参数中以获得后验：

n_days = [0, 1, 3, 5, 10, 25, 50, 100, 500]
outcomes = sp500_binary.sample(n_days[-1])
p = np.linspace(0, 1, 100)
# uniform (uninformative) prior
a = b = 1
for i, days in enumerate(n_days):
    up = outcomes.iloc[:days].sum()
    down = days - up
    update = stats.beta.pdf(p, a + up , b + down) 

结果后验分布已绘制在下图中。它们说明了从将所有成功概率视为等可能的均匀先验到越来越尖峰分布的演变。

经过 500 个样本后，从 2010 年到 2017 年，正向移动的概率集中在实际概率的 54.7％附近。它还显示了 MLE 和 MAP 估计之间的小差异，后者倾向于被稍微拉向均匀先验的期望值：

图 10.2：经过最多 500 次更新后的标普 500 指数第二天上涨概率的后验分布

在实践中，使用共轭先验的情况受到低维情况的限制。此外，简化的 MAP 方法避免了计算证据项，但有一个关键缺点，即使在可用时也是如此：它不返回分布，因此我们无法得出不确定性的度量或将其用作先验。因此，我们需要使用数值方法和随机模拟进行近似推断，而不是精确推断，接下来我们将介绍。

确定性和随机近似推断

对于大多数具有实际相关性的模型，将无法解析地推导出精确的后验分布并计算潜在参数的期望值。模型可能具有太多参数，或者后验分布可能对于解析解而言过于复杂：

• 对于连续变量，积分可能没有封闭形式的解，而空间的维数和被积函数的复杂性可能会阻止数值积分。

• 对于离散变量，边缘化涉及对所有可能的隐藏变量配置求和，尽管这在原则上总是可能的，但实践中我们经常发现可能存在指数多个隐藏状态，使得这种计算代价过高。

尽管对于一些应用来说，未观察参数的后验分布可能很重要，但通常主要是需要评估期望，例如进行预测。在这种情况下，我们可以依赖于近似推断，其中包括随机和确定性方法：

• 基于马尔可夫链蒙特卡罗(MCMC)抽样的随机技术已经在许多领域推广了贝叶斯方法的使用。它们通常具有收敛到精确结果的属性。在实践中，抽样方法可能计算成本高，通常仅限于小规模问题。

• 称为变分推断或变分贝叶斯的确定性方法基于对后验分布的解析近似，并且可以很好地扩展到大型应用程序。它们进行简化假设，例如，后验分解为特定方式或具有特定参数形式，例如高斯分布。因此，它们不生成精确结果，并且可以用作抽样方法的补充。

我们将在接下来的两个部分概述这两种方法。

马尔可夫链蒙特卡罗抽样

抽样是指从给定分布p(x)中绘制样本X=(x[1], …, x[n])。假设样本是独立的，大数定律确保对于增加的样本数量，样本中给定实例x[i]的比例（对于离散情况）对应于其概率p(x=x[i])。在连续情况下，类似的推理适用于样本空间的给定区域。因此，对样本的平均值可以用作分布参数的期望值的无偏估计量。

实际挑战之一是确保独立抽样，因为分布是未知的。依赖样本可能仍然是无偏的，但往往会增加估计的方差，因此需要更多的样本来获得与独立样本相同的精确估计。

从多元分布中抽样在状态数量随维度数呈指数增长时计算成本很高。许多算法简化了这个过程；我们将在这里介绍几种基于 MCMC 的流行变体。

马尔可夫链是描述一组由转移概率连接的状态上的随机游走的动态随机模型。马尔可夫属性规定了过程没有记忆，并且下一步仅取决于当前状态。换句话说，这取决于当前、过去和未来是否独立，也就是说，过去状态的信息不会帮助我们预测未来超出我们从现在所知道的内容。

蒙特卡罗方法依赖于重复随机抽样来近似可能是确定性的但不允许精确解析解的结果。它在曼哈顿计划期间开发，用于估算原子级别的能量，并获得了其持久的代码名称以确保保密性。

许多算法将蒙特卡罗方法应用于马尔可夫链，并通常按以下方式进行：

1. 从当前位置开始。

2. 从提议分布中绘制一个新位置。

3. 根据数据和先验分布评估新位置的概率。

   1. 如果足够可能，移动到新位置。

   2. 否则，保持在当前位置。

4. 从步骤 1重复。

5. 经过一定数量的迭代后，返回所有被接受的位置。

MCMC 方法旨在识别和探索后验中集中显著概率密度的有趣区域。当它连续地移动到后验的附近高概率状态时，被称为无记忆过程的过程收敛。一个关键挑战是平衡对样本空间随机探索的需求与降低接受率的风险。

过程的初始步骤可能更多地反映了起始位置而不是后验，并且通常被丢弃作为burn-in 样本。 MCMC 的一个关键属性是在一定（但未知）数量的迭代之后，该过程应“忘记”其初始位置。

剩下的样本被称为该过程的迹。假设收敛，则样本的相对频率近似后验，并且可以根据大数定律计算期望值。

如前所述，估计的精度取决于由随机游走收集的样本的串行相关性，其中每个样本设计上仅依赖于前一个状态。更高的相关性限制了对后验的有效探索，并且需要进行诊断测试。

设计这样一个马尔可夫链的一般技术包括 Gibbs 抽样、Metropolis-Hastings 算法和更近期的哈密顿 MCMC 方法，后者通常表现更好。

Gibbs 抽样

Gibbs 抽样将多变量抽样简化为一系列一维绘制。从某个起始点开始，它在抽样第n个变量时迭代地保持n-1 个变量不变。它包含了这个样本并重复它。

该算法非常简单易实现，但产生高度相关的样本，减慢了收敛速度。顺序性质也阻止了并行化。详见 Casella 和 George (1992) 中对其的详细描述和解释。

Metropolis-Hastings 取样

Metropolis-Hastings 算法根据当前状态随机提出新的位置。它这样做是为了有效地探索样本空间并减少相对于 Gibbs 采样的样本相关性。为了确保它从后验中取样，它使用先验和似然的乘积来评估提议，这与后验成比例。它接受的概率取决于结果相对于当前样本的相应值。

提议评估方法的一个关键优点是，它使用与后验的精确评估不同的比例评估。但是，它可能需要很长时间才能收敛。这是因为与后验无关的随机移动可能会降低接受率，从而使大量步骤只产生少量（可能相关的）样本。接受率可以通过减小提议分布的方差来调节，但由此产生的较小步骤意味着较少的探索。详见 Chib 和 Greenberg (1995) 中对该算法的详细介绍。

汉密尔顿蒙特卡洛 – 进入 NUTS

汉密尔顿蒙特卡洛 (HMC) 是一种混合方法，利用了似然梯度的一阶导数信息。借助此信息，它提出了新的状态以进行探索，并克服了一些 MCMC 的挑战。此外，它结合了动量以有效地在后验中“跳跃”。因此，它比简单的随机行走 Metropolis 或 Gibbs 采样更快地收敛到高维目标分布。详见 Betancourt (2018) 中全面的概念介绍。

无 U-Turn 采样器 (NUTS, Hoffman 和 Gelman 2011) 是一种自调整的 HMC 扩展，它在选择提议之前自适应地调节后验周围的移动的大小和数量。它在高维和复杂的后验分布上表现良好，并允许拟合许多复杂模型而无需关于拟合算法本身的专业知识。正如我们将在下一节中看到的，它是 PyMC3 中的默认采样器。

变分推断和自动微分

变分推断 (VI) 是一种通过优化来近似概率密度的机器学习方法。在贝叶斯背景下，它近似后验分布，如下所示：

1. 选择一个参数化的概率分布族

2. 找到该族中距离目标最近的成员，其度量为 Kullback-Leibler 散度

与 MCMC 相比，变分贝叶斯倾向于更快地收敛并且对大数据具有更好的可扩展性。虽然 MCMC 通过来自链条的样本逼近后验概率，这些样本最终将收敛到目标任意接近的地方，但变分算法通过优化的结果来逼近后验概率，不能保证与目标重合。

变分推断更适用于大型数据集，例如，数亿个文本文档，因此我们可以快速探索许多模型。相反，当时间和计算资源的约束较少时，MCMC 将在较小的数据集上提供更精确的结果。例如，如果您花了 20 年时间收集了一个小但昂贵的数据集，并且确信您的模型是适当的，并且您需要精确的推断，则 MCMC 将是一个不错的选择。有关更详细的比较，请参见 Salimans、Kingma 和 Welling（2015）。

变分推断的缺点是需要模型特定的导出和定制优化例程的实现，这减慢了广泛采用的速度。

最近的自动微分变分推断（ADVI）算法自动化了这个过程，使得用户只需指定模型，表达为一个程序，ADVI 会自动生成相应的变分算法（详见 GitHub 上的参考资料以了解实现细节）。

我们将看到 PyMC3 支持各种变分推断技术，包括 ADVI。

使用 PyMC3 的概率编程

概率编程提供了一种描述和拟合概率分布的语言，以便我们可以设计、编码并自动估计和评估复杂的模型。它旨在抽象出一些计算和分析复杂性，以便我们可以专注于贝叶斯推理和推断的概念上更直观、更简单的方面。

自 Uber 开源了 Pyro（基于 PyTorch）后，该领域变得非常动态。最近，Google 为 TensorFlow 添加了一个概率模块

因此，贝叶斯方法在机器学习中的实际相关性和使用可能会增加，以生成对不确定性有洞察力的见解，特别是对于需要透明而不是黑匣子模型的用例。

在本节中，我们将介绍流行的PyMC3库，该库使用 Python 为 ML 模型实现了高级 MCMC 抽样和变分推断。与 Stan（以蒙特卡罗方法发明者 Stanislaw Ulam 命名，由哥伦比亚大学的 Andrew Gelman 自 2012 年以来开发）一起，PyMC3 是最受欢迎的概率编程语言。

使用 Theano 的贝叶斯机器学习

PyMC3 于 2017 年 1 月发布，以将 Hamiltonian MC 方法添加到 PyMC2（2012 年发布）中使用的 Metropolis-Hastings 采样器中。PyMC3 使用 Theano 作为其计算后端，用于动态 C 编译和自动微分。Theano 是由 Yoshua Bengio 的蒙特利尔学习算法研究所（MILA）开发的一个以矩阵为重点且支持 GPU 的优化库，这启发了 TensorFlow。由于更新的深度学习库的成功，MILA 最近停止进一步开发 Theano（详见第十六章，收益电话和 SEC 文件的单词嵌入）。

PyMC4 于 2019 年 12 月发布 alpha 版，使用 TensorFlow 代替 Theano，并旨在限制对 API 的影响（请参阅 GitHub 上的存储库链接）。

PyMC3 工作流程 - 预测衰退

PyMC3 旨在实现直观、易读且功能强大的语法，反映了统计学家描述模型的方式。建模过程通常遵循以下三个步骤：

1. 通过定义编码概率模型：

   1. 量化潜变量的知识和不确定性的先验分布

   2. 将条件参数化为观察数据的似然函数

2. 使用前一节描述的选项之一分析后验：

   1. 使用 MAP 推断获取一个点估计

   2. 使用 MCMC 方法从后验中采样

   3. 使用变分贝叶斯近似后验

3. 使用各种诊断工具检查您的模型

4. 生成预测

生成的模型可用于推理，以获得有关参数值的详细见解，以及对新数据点的预测结果。

我们将使用简单的 logistic 回归来说明这个工作流程，以模拟衰退的预测（请参阅笔记本pymc3_workflow）。随后，我们将使用 PyMC3 来计算和比较贝叶斯夏普比率，估计动态配对交易比率，并实现贝叶斯线性时间序列模型。

数据 - 领先的衰退指标

我们将使用一个小型和简单的数据集，以便我们可以专注于工作流程。我们将使用美联储经济数据（FRED）服务（详见第二章，市场和基本数据 - 来源和技术）下载由美国经济研究局（NBER）定义的美国衰退日期。我们还将获取四个通常用于预测衰退发生的变量（Kelley 2019），这些变量可以通过 FRED 获得，即：

• 长期国债收益率曲线的传播，定义为 10 年期和 3 个月期国债收益率之间的差异

• 密歇根大学的消费者情绪指数

• 全国金融状况指数（NFCI）

• NFCI 非金融杠杆子指数

衰退日期以季度为基础确定；我们将重新采样所有系列的频率为每月频率，以从 1982 年至 2019 年获得大约 457 个观测值。如果一个季度被标记为衰退，我们将考虑该季度中的所有月份都是衰退的。

我们将构建一个意图回答问题的模型：未来 x 个月美国经济是否会陷入衰退？ 换句话说，我们不只关注预测衰退的第一个月；这限制了 48 个衰退月的不平衡性。

为此，我们需要选择一个先导时间；已经进行了大量研究，以找到各种领先指标的适当时间范围：收益率曲线倾向于在衰退之前提前至多 24 个月发送信号；NFCI 指标倾向于具有较短的先导时间（参见 Kelley, 2019）。

以下表格在很大程度上证实了这种经验：它显示了二元衰退变量与四个领先指标之间的互信息（参见第六章，机器学习流程）在 1-24 个月的时间段内的情况：

图 10.3：衰退与领先指标之间的互信息，时间段为 1-24 个月

为了在 NFCI 指标的较短时间范围和收益率曲线之间取得平衡，我们将选择 12 个月作为我们的预测时间范围。 以下图表显示了每个指标的分布，按衰退状态划分：

图 10.4：按衰退状态划分的领先指标分布

这表明，经济衰退往往与国债收益率曲线的长期差价负相关，也被称为倒挂收益率曲线，当短期利率上升到长期利率之上时。 NFCI 指标的表现符合我们的预期；情绪指标似乎具有最弱的关联性。

模型定义 - 贝叶斯逻辑回归

如第六章，机器学习流程中讨论的那样，逻辑回归估计了一组特征与二元结果之间的线性关系，通过一个 Sigmoid 函数进行中介，以确保模型产生概率。 频率主义方法导致了参数的点估计，这些参数测量了每个特征对数据点属于正类的概率的影响，置信区间基于参数分布的假设。

相比之下，贝叶斯逻辑回归估计参数本身的后验分布。 后验允许更健壮地估计每个参数的所谓贝叶斯可信区间，其好处在于更透明地了解模型的不确定性。

概率程序由观察和未观察的随机变量（RVs）组成。 正如之前讨论的，我们通过似然分布定义观察到的 RVs，通过先验分布定义未观察到的 RVs。 PyMC3 包含许多概率分布，用于此目的。

PyMC3 库使得执行逻辑回归的近似贝叶斯推断非常简单。 逻辑回归模型基于左侧图中所述的 k 个特征，来建模经济在 i 个月后是否会陷入衰退的概率：

图 10.5：贝叶斯逻辑回归

我们将使用上下文管理器with来定义一个我们稍后可以引用的manual_logistic_model作为概率模型：

1. 未观测参数的 RVs，用于拦截和两个特征的被表达使用不含信息的先验，这些假设正态分布的均值为 0，标准差为 100。

2. 似然性根据逻辑回归的规范将参数与数据结合起来。

3. 结果被建模为具有由似然给定的成功概率的 Bernoulli RV：

   with pm.Model() as manual_logistic_model:
       # coefficients as rvs with uninformative priors
       intercept = pm.Normal('intercept', 0, sd=100)
       beta_1 = pm.Normal('beta_1', 0, sd=100)
       beta_2 = pm.Normal('beta_2', 0, sd=100)
       # Likelihood transforms rvs into probabilities p(y=1)
       # according to logistic regression model.
       likelihood = pm.invlogit(intercept +
                                beta_1 * data.yield_curve +
                                beta_2 * data.leverage)
       # Outcome as Bernoulli rv with success probability
       # given by sigmoid function conditioned on actual data
       pm.Bernoulli(name='logit',
                    p=likelihood,
                    observed=data.recession) 
   

模型可视化和板符号

命令pm.model_to_graphviz(manual_logistic_model)生成右侧图 10.5中显示的板符号。 它显示了未观察到的参数作为浅色椭圆，观察到的元素作为深色椭圆。 矩形表示在模型定义中包含的数据暗示的观察到的模型元素的重复次数。

广义线性模型模块

PyMC3 包括许多常见的模型，以便我们可以限制自定义应用程序的手动规范。

以下代码将相同的逻辑回归定义为广义线性模型 (GLM)家族的成员。它使用受统计语言 R 启发的公式格式，并由 patsy 库移植到 Python 中：

with pm.Model() as logistic_model:
    pm.glm.GLM.from_formula(recession ~ yield_curve + leverage,
                            data,
                            family=pm.glm.families.Binomial()) 

精确 MAP 推断

我们使用刚刚定义的模型的.find_MAP()方法获得三个参数的点 MAP 估计。 如预期的那样，较低的展开值增加了衰退的概率，较高的杠杆(但程度较小)也是如此：

with logistic_model:
    map_estimate = pm.find_MAP()
print_map(map_estimate)
Intercept     -4.892884
yield_curve   -3.032943
leverage       1.534055 

PyMC3 使用拟牛顿Broyden-Fletcher-Goldfarb-Shanno (BFGS)算法解决找到具有最高密度的后验点的优化问题，但提供了 SciPy 库提供的几种替代方案。

MAP 点估计与相应的statsmodels系数相同(请参阅笔记本pymc3_workflow)。

近似推断 – MCMC

如果我们只对模型参数的点估计感兴趣，那么对于这个简单模型来说，MAP 估计就足够了。 更复杂的自定义概率模型需要采样技术以获得参数的后验概率。

我们将使用所有变量的模型来说明 MCMC 推断：

formula = 'recession ~ yield_curve + leverage + financial_conditions + sentiment'
with pm.Model() as logistic_model:
    pm.glm.GLM.from_formula(formula=formula,
                            data=data,
                            family=pm.glm.families.Binomial())
# note that pymc3 uses y for the outcome
logistic_model.basic_RVs
[Intercept, yield_curve, leverage, financial_conditions, sentiment, y] 

请注意，测量在非常不同尺度上的变量可能会减慢采样过程。 因此，我们首先应用由 scikit-learn 提供的scale()函数来标准化所有特征。

一旦我们用新的公式定义了我们的模型，我们就可以执行推断以近似后验分布。 MCMC 采样算法可通过pm.sample()函数获得。

默认情况下，PyMC3 自动选择最有效的采样器并初始化采样过程以实现高效收敛。对于连续模型，PyMC3 选择在前一节中讨论的 NUTS 采样器。它还通过 ADVI 运行变分推断以找到采样器的良好起始参数。其中一种替代方法是使用 MAP 估计。

要了解收敛的情况，我们首先在调整采样器 1,000 次迭代后只绘制 100 个样本。这些将被丢弃。可以使用cores参数（除了使用 GPU 时）将采样过程并行化为多个链：

with logistic_model:
    trace = pm.sample(draws=100,
                      tune=1000,
                      init='adapt_diag',
                      chains=4,
                      cores=4,
                      random_seed=42) 

结果的trace包含每个 RV 的采样值。我们可以使用plot_traces()函数检查链的后验分布。

plot_traces(trace, burnin=0) 

图 10.6显示了前两个特征和截距的样本分布及其随时间的值（请参阅笔记本以获取完整输出）。此时，采样过程尚未收敛，因为对于每个特征，四条轨迹产生了完全不同的结果；左侧五个面板中垂直显示的数字是由四条轨迹生成的分布的模式的平均值：

图 10.6：100 个样本后的轨迹

我们可以通过提供先前运行的迹作为输入来继续采样。额外 20,000 个样本后，我们观察到一个非常不同的图像，如下图所示。这显示了采样过程现在更接近收敛。还要注意，初始系数点估计值与当前值相对较接近：

图 10.7：额外 50,000 个样本后的轨迹

我们可以计算可信区间，即贝叶斯区间估计的对应物，作为迹的百分位数。结果边界反映了我们对于给定概率阈值的参数值范围的置信度，而不是参数在大量试验中在此范围内的次数。图 10.8显示了变量的收益曲线和杠杆的可信区间，以系数值的指数函数（参见第七章，线性模型 - 从风险因素到收益预测）来表示。

请参阅笔记本pymc3_workflow的实现：

图 10.8：收益曲线和杠杆的可信区间

近似推断 - 变分贝叶斯

变分推断的接口与 MCMC 实现非常相似。我们只是使用fit()而不是sample()函数，可以选择包括一个早期停止的CheckParametersConvergence回调，如果分布拟合过程收敛到给定的容差：

with logistic_model:
    callback = CheckParametersConvergence(diff='absolute')
    approx = pm.fit(n=100000,
                    callbacks=[callback]) 

我们可以从近似分布中抽取样本以获得一个迹对象，就像我们之前为 MCMC 采样器所做的那样：

trace_advi = approx.sample(10000) 

检查跟踪摘要显示结果稍微不太准确。

模型诊断

贝叶斯模型诊断包括验证采样过程是否收敛并且是否一直从后验的高概率区域进行采样，以及确认模型是否很好地表示了数据。

收敛性

我们可以通过随时间变化的样本和它们的分布来可视化结果的质量。下图显示了初始 100 个样本和额外的 200,000 个样本后的后验分布，说明了收敛意味着多个链识别相同的分布：

图 10.9：400 个样本后和超过 200,000 个样本后的跟踪

PyMC3 为采样器生成各种摘要统计信息。这些信息可以作为统计模块中的单独函数提供，也可以通过将追踪传递给函数pm.summary()来获取。

下表包含（分别计算的）statsmodels 逻辑斯蒂回归系数作为第一列，以显示在这种简单情况下，两个模型略有一致，因为样本均值与系数不匹配。这很可能是由于准分离的高程度造成的：收益率曲线的高可预测性允许对 17%的数据点进行完美预测，这反过来导致逻辑回归的最大似然估计定义不清晰（有关更多信息，请参见笔记本中的 statsmodels 输出）。

参数	statsmodels	PyMC3
系数	均值	标准差	HPD 3%	HPD 97%	有效样本	R hat
截距	-5.22	-5.47	0.71	-6.82	-4.17	68,142	1.00
收益率曲线	-3.30	-3.47	0.51	-4.44	-2.55	70,479	1.00
杠杆	1.98	2.08	0.40	1.34	2.83	72,639	1.00
金融条件	-0.65	-0.70	0.33	-1.33	-0.07	91,104	1.00
情绪	-0.33	-0.34	0.26	-0.82	0.15	106,751	1.00

剩余列包含最高后验密度（HPD）估计的最小宽度可信区间，这是置信区间的贝叶斯版本，此处计算为 95%的水平。n_eff统计量总结了从绘制的有效（未被拒绝）样本数量。

R-hat，也被称为Gelman-Rubin 统计量，通过比较链间方差和链内方差来检查收敛性。如果采样器收敛了，这些方差应该相同，即链应该看起来相似。因此，该统计量应接近 1。

对于具有许多变量的高维模型，检查大量的跟踪变得很麻烦。在使用 NUTS 时，能量图帮助我们评估收敛问题。它总结了随机过程探索后验的效率。图表显示了能量和能量转换矩阵，它们应该匹配良好，就像下图右侧面板中显示的例子一样：

图 10.10: 森林和能量图

后验预测检查

后验预测检查 (PPCs) 对于检查模型与数据拟合程度非常有用。它们通过使用从后验抽取的参数从模型生成数据来实现这一点。我们使用函数 pm.sample_ppc 来实现此目的，并为每个观测值获取 n 个样本（GLM 模块自动将结果命名为 'y'）：

ppc = pm.sample_ppc(trace_NUTS, samples=500, model=logistic_model)
ppc['y'].shape
(500, 445) 

我们可以使用接收器操作特征曲线下面积（AUC，参见第六章，机器学习过程）得分来评估样本内适配度，例如，比较不同的模型：

roc_auc_score(y_score=np.mean(ppc['y'], axis=0),
              y_true=data.income)
0.9483627204030226 

结果相当高，几乎为 0.95。

如何生成预测

预测使用 Theano 的 共享变量 来在运行后验预测检查之前用测试数据替换训练数据。为了可视化和简化表述，我们使用收益曲线变量作为唯一的预测变量，并忽略我们数据的时间序列特性。

相反，我们使用 scikit-learn 的基本 train_test_split() 函数创建训练集和测试集，通过结果进行分层，以保持类别不平衡：

X = data[['yield_curve']]
labels = X.columns
y = data.recession
X_train, X_test, y_train, y_test = train_test_split(X, y,
                                                    test_size=0.2,
                                                    random_state=42,
                                                    stratify=y) 

然后我们为该训练集创建一个共享变量，然后在下一步中用测试集替换它。请注意，我们需要使用 NumPy 数组并提供列标签的列表：

X_shared = theano.shared(X_train.values)
with pm.Model() as logistic_model_pred:
    pm.glm.GLM(x=X_shared, labels=labels,
               y=y_train, family=pm.glm.families.Binomial()) 

然后我们像之前一样运行采样器：

with logistic_model_pred:
    pred_trace = pm.sample(draws=10000,
                           tune=1000,
                           chains=2,
                           cores=2,
                           init='adapt_diag') 

现在，我们将测试数据替换为共享变量上的训练数据，并将 pm.sample_ppc 函数应用于生成的 trace：

X_shared.set_value(X_test)
ppc = pm.sample_ppc(pred_trace,
                    model=logistic_model_pred,
                    samples=100)
y_score = np.mean(ppc['y'], axis=0)
roc_auc_score(y_score=np.mean(ppc['y'], axis=0),
              y_true=y_test)
0.8386 

这个简单模型的 AUC 分数是 0.86。显然，如果训练集已经包含了来自附近月份的衰退示例，那么对于另一个月份预测相同的衰退要容易得多。请记住，我们仅用于演示目的使用此模型。

图 10.11 绘制了从 100 个 Monte Carlo 链中抽样的预测及其周围的不确定性，以及实际的二元结果和对应于模型预测的 logistic 曲线：

图 10.11: 单变量模型预测

摘要和主要观点

我们建立了一个简单的 logistic 回归模型，用于预测美国经济在 12 个月内陷入衰退的概率，使用了四个领先指标。对于这个简单模型，我们可以得到精确的 MAP 估计系数值，然后使用这些系数值对模型进行参数化和预测。

但是，更复杂的自定义概率模型将不允许此捷径，而且 MAP 估计也不会生成有关点估计之外的后验分布的洞见。因此，我们演示了如何使用 PyMC3 进行近似推理。结果说明了我们如何了解每个模型参数的后验分布，但也显示出即使对于一个小模型，与 statsmodels MLE 估计相比，计算成本也会显著增加。尽管如此，对于复杂的概率模型，基于采样的解决方案是了解数据的唯一途径。

我们现在将说明如何将贝叶斯分析应用于一些与交易相关的用例。

交易的贝叶斯机器学习

既然我们熟悉了 ML 和用 PyMC3 进行概率编程的贝叶斯方法，让我们探讨一下几个相关的交易应用，即：

• 将夏普比率建模为更具洞察力的性能比较的概率模型

• 使用贝叶斯线性回归计算配对交易对冲比率

• 从贝叶斯角度分析线性时间序列模型

托马斯·韦基，PyMC3 的主要作者之一，也是 Quantopian 的数据科学主管，已经创建了几个示例，以下各节将继续并在其基础上建立。PyMC3 文档有许多额外的教程（请查看 GitHub 获取链接）。

用于性能比较的贝叶斯夏普比率

在本节中，我们将说明：

• 如何使用 PyMC3 定义夏普比率（SR）作为概率模型

• 如何比较不同回报序列的后验分布

为两个序列进行贝叶斯估计提供了非常丰富的见解，因为它提供了效应大小，组 SR 均值及其差异的可信值的完整分布，以及标准偏差及其差异的完整分布。Python 实现是由 Thomas Wiecki 完成的，并受到了 R 包 BEST（Meredith 和 Kruschke，2018）的启发。

使用贝叶斯 SR 的相关用例包括分析替代策略之间的差异，或者策略的样本内回报与样本外回报之间的差异（详见笔记本bayesian_sharpe_ratio）。贝叶斯 SR 也是 pyfolio 的贝叶斯信息单页的一部分。

定义一个自定义的概率模型

要将 SR 建模为概率模型，我们需要关于回报分布及其控制此分布的参数的先验。相对于低自由度（DF）的正态分布，学生 t 分布具有较厚的尾部，是捕获回报此方面的合理选择。

因此，我们需要对此分布的三个参数进行建模，即回报的均值和标准偏差，以及 DF。我们将假设对于均值和标准偏差，分别采用正态和均匀分布，对于 DF 采用指数分布，期望值足够低以确保尾部较厚。

这些概率输入是基于这些回报的，而年化 SR 来自标准计算，忽略无风险利率（使用每日回报）。我们将提供 2010-2018 年的 AMZN 股票回报作为输入（有关数据准备的更多信息，请参阅笔记本）：

mean_prior = data.stock.mean()
std_prior = data.stock.std()
std_low = std_prior / 1000
std_high = std_prior * 1000
with pm.Model() as sharpe_model:
    mean = pm.Normal('mean', mu=mean_prior, sd=std_prior)
    std = pm.Uniform('std', lower=std_low, upper=std_high)
    nu = pm.Exponential('nu_minus_two', 1 / 29, testval=4) + 2
    returns = pm.StudentT('returns', nu=nu, mu=mean, sd=std,
observed=data.stock)
    sharpe = returns.distribution.mean / returns.distribution.variance ** 
.5 * np.sqrt(252)
    pm.Deterministic('sharpe', sharpe) 

板符号，我们在关于 PyMC3 工作流程的上一节中引入的，可视化了三个参数及其关系，以及我们在以下图中提供的回报和观察次数：

图 10.12：贝叶斯 SR 在板符号中

然后我们运行我们在前一节介绍的 MCMC 抽样过程（请参阅笔记本bayesian_sharpe_ratio以了解随后的实现细节）。经过约 25000 次对四个链的抽样后，我们获得了模型参数的后验分布，结果显示在以下图中：

plot_posterior(data=trace); 

图 10.13：模型参数的后验分布

现在我们知道如何评估单个资产或投资组合的 SR，让我们看看如何使用贝叶斯 SR 来比较两个不同回报序列的性能。

比较两个回报序列的性能

为了比较两个回报序列的性能，我们将分别对每个组的 SR 进行建模，并将效应大小计算为波动率调整后回报之间的差异。以下图表显示的相应概率模型自然更大，因为它包括两个 SR 以及它们的差异：

图 10.14：在板符号中两个贝叶斯 SR 之间的差异

一旦我们定义了模型，我们就通过 MCMC 抽样过程运行它以获得其参数的后验分布。我们使用 2010-2018 年的 2,037 个每日回报来比较 AMZN 股票与同一时期的标普 500 指数回报。我们可以使用我们任何策略回测的回报，而不是 AMZN 的回报。

可视化轨迹揭示了对每个指标分布的精细性能洞察，如图 10.15中的各种图表所示：

图 10.15：两个贝叶斯 SR 之间的后验分布

最重要的度量标准是底部面板中两个 SR 之间的差异。给定完整的后验分布，可以直观地可视化或计算从 SR 角度来看一个回报序列优越的概率。

用于配对交易的贝叶斯滚动回归

在前一章中，我们介绍了配对交易作为一种流行的交易策略，它依赖于两个或更多资产的协整性。 鉴于这样的资产，我们需要估计对冲比率，以决定多头和空头位置的相对大小。 一种基本方法使用线性回归。 你可以在笔记本rolling_regression中找到此部分的代码，该笔记本遵循了托马斯·威克（Thomas Wiecki）的滚动回归示例（请参阅 GitHub 上 PyMC3 教程的链接）。

配对交易候选资产的一个流行示例是 ETF GLD，它反映了金价和像 GFI 这样的金矿股票。 我们使用 yfinance 获取 2004 年至 2020 年期间的收盘价数据。 图 10.16的左侧面板显示了历史价格序列，而右侧面板显示了历史价格的散点图，其中色调表示时间维度，以突出显示相关性的演变方式。 注意我们应该使用收益率，就像我们在第九章，波动率预测和统计套利的时间序列模型中所做的那样，来计算对冲比率; 然而，使用价格序列创建更引人注目的可视化效果。 建模过程本身保持不变：

图 10.16：两对交易候选资产的价格序列和随时间的相关性

我们想说明滚动贝叶斯线性回归如何随时间跟踪两个资产价格之间关系的变化。 主要思想是通过允许回归系数的变化将时间维度纳入线性回归中。 具体来说，我们将假设截距和斜率随时间呈随机游走：

我们使用 PyMC3 内置的pm.GaussianRandomWalk过程指定model_randomwalk。 它要求我们为截距 alpha 和斜率 beta 都定义标准差：

model_randomwalk = pm.Model()
with model_randomwalk:
    sigma_alpha = pm.Exponential('sigma_alpha', 50.)
    alpha = pm.GaussianRandomWalk('alpha', 
                                  sd=sigma_alpha,
                                  shape=len(prices))
    sigma_beta = pm.Exponential('sigma_beta', 50.)
    beta = pm.GaussianRandomWalk('beta', 
                                 sd=sigma_beta,
                                 shape=len(prices)) 

鉴于概率模型的规范，我们现在将定义回归并将其连接到输入数据：

with model_randomwalk:
    # Define regression
    regression = alpha + beta * prices_normed.GLD
    # Assume prices are normally distributed
    # Get mean from regression.
    sd = pm.HalfNormal('sd', sd=.1)
    likelihood = pm.Normal('y', 
                           mu=regression, 
                           sd=sd, 
                           observed=prices_normed.GFI) 

现在，我们可以运行 MCMC 采样器以生成模型参数的后验分布：

with model_randomwalk:
    trace_rw = pm.sample(tune=2000, 
                         cores=4, 
                         draws=200, 
                         nuts_kwargs=dict(target_accept=.9)) 

图 10.17描述了截距和斜率系数随时间的变化，强调了演变的相关性：

图 10.17：截距和斜率系数随时间的变化

使用动态回归系数，我们现在可以通过这种贝叶斯方法来可视化滚动回归建议的对冲比率随时间的变化，该方法将系数建模为随机游走。

下图将价格序列和回归线结合在一起，其中色调再次表示时间线（在笔记本中查看以获得彩色输出）：

图 10.18：滚动回归线和价格序列

对于我们的最后一个示例，我们将实现一个贝叶斯随机波动模型。

随机波动模型

正如前一章所讨论的，资产价格具有时变波动性。在某些时期，收益非常不稳定，而在其他时期则非常稳定。我们在 第九章，波动性预测和统计套利的时间序列模型 中从经典线性回归的角度探讨了 ARCH/GARCH 模型。

贝叶斯随机波动性模型通过潜在波动性变量捕捉这种波动现象，该变量被建模为随机过程。使用这种模型介绍了无 U 转向采样器 (Hoffman, et al. 2011)，笔记本 stochastic_volatility 使用了 S&P 500 每日数据来说明这种用例。图 10.19 显示了此期间的几个波动性集群：

图 10.19：每日 S&P 500 对数收益率

概率模型规定对数收益率遵循 t 分布，该分布具有厚尾，如资产收益通常观察到的情况。 t 分布受参数 ν 控制，表示自由度。它也被称为正态参数，因为 t 分布在 ν 增加时逼近正态分布。假定此参数具有参数 的指数分布。

此外，假定对数收益率的平均值为零，而标准差遵循具有指数分布的标准差随机游走：

我们将此模型实现为 PyMC3，以模仿其概率规范，使用对数收益率匹配模型：

prices = pd.read_hdf('../data/assets.h5', key='sp500/prices').loc['2000':,
                                                                  'Close']
log_returns = np.log(prices).diff().dropna()
with pm.Model() as model:
    step_size = pm.Exponential('sigma', 50.)
    s = GaussianRandomWalk('s', sd=step_size, 
                           shape=len(log_returns))
    nu = pm.Exponential('nu', .1)
    r = pm.StudentT('r', nu=nu, 
                    lam=pm.math.exp(-2*s), 
                    observed=log_returns) 

接下来，在经过 2,000 个样本的燃烧阶段后，我们绘制了 5,000 个 NUTS 样本，使用比默认值 0.8 更高的接受率，根据 PyMC3 文档建议，对于有问题的后验分布（请参阅 GitHub 上的适当链接）：

with model:
    trace = pm.sample(tune=2000, 
                      draws=5000,
                      nuts_kwargs=dict(target_accept=.9))
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [nu, s, sigma]
Sampling 4 chains, 0 divergences: 100%|██████████| 28000/28000 [27:46<00:00, 16.80draws/s]
The estimated number of effective samples is smaller than 200 for some parameters. 

在四个链总共 28,000 个样本之后，下图中的迹线图证实了采样过程已经收敛：

图 10.20：随机波动率模型的迹线图

当我们将样本绘制到 图 10.21 中的 S&P 500 收益率上时，我们可以看到这个简单的随机波动性模型相当好地跟踪了波动性集群：

图 10.21：模型

记住这代表了样本内拟合。作为下一步，你应该尝试评估预测准确性。我们在前面滚动线性回归的子章节中介绍了如何进行预测，并在前几章中使用了时间序列交叉验证，这为您提供了完成此目的所需的所有工具！

概要

在本章中，我们探讨了贝叶斯方法与机器学习的结合。我们看到它们有几个优点，包括能够编码先验知识或意见，更深入地了解模型估计和预测周围的不确定性，并适用于在线学习，其中每个训练样本逐渐影响模型的预测。

我们学习了如何使用 PyMC3 应用贝叶斯工作流程，从模型规范到估计、诊断和预测，并探讨了几个相关的应用场景。我们将在第十四章中遇到更多的贝叶斯模型，用于交易的文本数据 - 情感分析，在那里我们将讨论自然语言处理和主题建模，以及在第二十章中，用于条件风险因素和资产定价的自动编码器，在那里我们将介绍变分自动编码器。

下一章介绍了非线性的基于树的模型，即决策树，并展示了如何将多个模型组合成一组树的集合，以创建一个随机森林。

第十一章：随机森林 - 用于日本股票的长短策略

在本章中，我们将学习如何使用两种新的机器学习模型进行交易：决策树和随机森林。我们将看到决策树如何从数据中学习规则，这些规则编码了输入和输出变量之间的非线性关系。我们将说明如何训练决策树，并使用它进行回归和分类问题的预测，可视化并解释模型学到的规则，并调整模型的超参数以优化偏差-方差权衡，并防止过拟合。

决策树不仅是重要的独立模型，而且经常被用作其他模型的组成部分。在本章的第二部分中，我们将介绍集成模型，这些模型将多个单独模型组合起来，产生一个具有更低预测误差方差的单一聚合预测。

我们将阐述自助聚合，通常称为bagging，作为几种方法之一，用于随机化个体模型的构建，并减少由集合成员的预测误差造成的相关性。我们将说明如何通过 bagging 有效地减少方差，并学习如何配置、训练和调整随机森林。我们将看到随机森林作为一个（可能较大的）决策树集合，可以显著减少预测误差，但可能会牺牲一些解释性。

然后，我们将继续构建一个使用随机森林为过去 3 年大型日本股票生成盈利信号的长短交易策略。我们将获取并准备股价数据，调整随机森林模型的超参数，并根据模型的信号进行交易规则的回测。由此产生的长短策略使用机器学习而不是我们在第九章中看到的共整关系，以识别和交易在给定投资期限内价格可能朝相反方向移动的证券篮子。

简而言之，阅读完本章后，你将能够：

• 使用决策树进行回归和分类

• 从决策树中获得见解，并可视化从数据中学到的规则

• 理解为什么集成模型往往能够提供优越的结果

• 使用自助聚合来解决决策树过拟合的挑战

• 训练、调整和解释随机森林

• 使用随机森林设计和评估盈利交易策略

你可以在 GitHub 仓库的相应目录中找到本章的代码示例和其他资源的链接。笔记本包括图像的彩色版本。

决策树 - 从数据中学习规则

决策树是一种机器学习算法，它根据从数据中学习到的决策规则来预测目标变量的值。该算法可以通过改变管理树如何学习规则的目标函数应用于回归和分类问题。

我们将讨论决策树如何使用规则进行预测，如何训练它们来预测（连续的）收益以及（分类的）价格走势方向，以及如何有效地解释、可视化和调整它们。有关更多详细信息和背景信息，请参阅 Rokach 和 Maimon（2008）以及 Hastie、Tibshirani 和 Friedman（2009）。

树是如何学习和应用决策规则的

我们在第七章，线性模型 - 从风险因素到收益预测和第九章，时间序列模型用于波动率预测和统计套利中学习的线性模型通过学习一组参数来预测输出，可能在逻辑回归的情况下通过 S 形链接函数进行转换。

决策树采用不同的方法：它们学习并依次应用一组规则，将数据点分成子集，然后为每个子集做出一个预测。这些预测基于应用给定规则序列所导致的训练样本子集的结果值。分类树预测从相对类频率或最多类的值直接估计的概率，而回归树计算可用数据点的结果值均值的预测。

每一条规则都依赖于一个特定的特征，并使用一个阈值将样本分成两组，其中值要么低于要么高于该特征的阈值。二叉树自然地表示了模型的逻辑：根是所有样本的起点，节点代表决策规则的应用，数据沿着边移动，当它被分成更小的子集时，直到到达叶节点，模型进行预测。

对于线性模型，参数值允许解释输入变量对输出和模型预测的影响。相反，对于决策树，从根到叶的各种可能路径确定了特征及其值如何导致模型做出具体决策。因此，决策树能够捕捉线性模型无法“开箱即用”捕捉的特征之间的相互依赖。

以下图表突出显示了模型如何学习一条规则。在训练过程中，算法扫描特征，并且对于每个特征，它寻求找到一个分割数据以最小化由预测造成的损失的切分点。它使用将结果来自于分割的子集，按每个子集中的样本数量加权：

图 11.1：决策树如何从数据中学习规则

在训练期间构建整个树，学习算法重复此过程，即将特征空间，即p个输入变量X[1]、X[2]、...、X[p]的可能值集合，划分为互斥且集体穷尽的区域，每个区域由一个叶节点表示。不幸的是，由于特征空间的可能组合和阈值序列的爆炸性数量，算法将无法评估特征空间的每种可能分区。基于树的学习采用了一种自顶向下、贪婪的方法，称为递归二元分割，以克服这种计算限制。

此过程是递归的，因为它使用来自先前分割的数据子集。它是自顶向下的，因为它从树的根节点开始，所有观察仍然属于单个区域，然后通过向预测器空间添加一个以上的分割来连续创建树的两个新分支。它是贪婪的，因为算法根据对目标函数的直接影响选择最佳规则，而不是向前看并评估数步之后的损失。我们将在更具体的回归和分类树的上下文中返回分割逻辑，因为这代表了它们之间的主要差异。

训练样本数量随着递归分割向树中添加新节点而不断减少。如果规则将样本均匀分割，导致树完美平衡，每个节点都有相同数量的子节点，那么在第n级就会有 2^n 个节点，每个节点包含总观测数的相应部分。实际上，这是不太可能的，因此沿某些分支的样本数量可能会迅速减少，并且树倾向于沿不同路径生长到不同的深度。

递归分割将继续，直到每个叶节点仅包含单个样本，并且训练误差已经降低到零。我们将介绍几种方法来限制分割并防止决策树产生极端过拟合的自然倾向。

为了对新的观测值进行预测，模型使用在训练期间推断出的规则来决定数据点应分配到哪个叶节点，然后使用特征空间相应区域中训练观测的平均值（用于回归）或模式（用于分类）。特征空间中给定区域（即给定叶节点）中训练样本数量较少，会降低预测的置信度，并可能反映出过拟合。

实践中的决策树

在本节中，我们将说明如何使用基于树的模型来获得洞察并进行预测。为了演示回归树，我们预测收益，对于分类案例，我们回到了正向和负向资产价格变动的示例。本节的代码示例位于笔记本decision_trees中，除非另有说明。

数据 - 月度股票收益和特征

我们将选择涵盖 2006-2017 年期间的 Quandl 美国股票数据集的子集，并按照我们第一个特征工程示例中的过程进行操作，第四章，金融特征工程 - 如何研究阿尔法因子。我们将计算月度收益和基于它们的 5 年移动平均值的 500 种最常交易的股票的 25 个（希望是）预测性特征，产生 56,756 个观察值。这些特征包括：

• 过去 1、3、6 和 12 个月的历史收益。

• 动量指标将最近 1 个或 3 个月的收益与较长时间跨度的收益相关联。

• 设计用于捕捉波动性的技术指标，如（归一化的）平均真实范围（NATR 和 ATR）和像相对强弱指数（RSI）这样的动量指标。

• 根据滚动 OLS 回归的五个 Fama-French 因子的因子加载。

• 年份和月份以及部门的分类变量。

图 11.2显示了这些特征与我们用于回归的月度收益之间的互信息（左侧面板），以及它们的二值化分类对应物，代表了相同期间的正向或负向价格变动。它显示，从单变量的角度来看，无论是对于这些特征的哪一种结果，都存在着信号内容的显著差异。

更多细节可以在这一章的 GitHub 存储库中的data_prep笔记本中找到。本章的决策树模型不具备处理缺失或分类变量的能力，因此我们将放弃前者并对分类部门变量应用虚拟编码（参见第四章，金融特征工程 - 如何研究阿尔法因子和第六章，机器学习过程）：

图 11.2：特征与收益或价格变动方向的互信息

使用时间序列数据构建回归树

回归树根据分配给给定节点的训练样本的平均结果值进行预测，并且通常依靠平均平方误差在递归二进制分割过程中选择最佳规则。

给定一个训练集，算法在p个预测变量 X[1]、X[2]、...、X[p] 和n个可能的切分点 s[1]、s[2]、...、s[n] 上进行迭代，以找到最佳组合。最佳规则将特征空间分割成两个区域，{X|X[i] < s[j]} 和 {X|X[i] > s[j]}，其中 X[i] 特征的值要么在 s[j] 阈值以下，要么在 s[j] 阈值以上，以便基于训练子集的预测最大程度地减少相对于当前节点的平方残差。

让我们从一个简化的示例开始，以便进行可视化，并演示如何使用时间序列数据与决策树。我们将只使用 2 个月的滞后回报来预测以下月份，与前一章中的 AR(2) 模型类似：

使用 scikit-learn，配置和训练回归树非常简单：

from sklearn.tree import DecisionTreeRegressor
# configure regression tree
regression_tree = DecisionTreeRegressor(criterion='mse',      
                                        max_depth=6,         
                                        min_samples_leaf=50)
# Create training data
y = data.target
X = data.drop(target, axis=1)
X2 = X.loc[:, ['t-1', 't-2']]
# fit model
regression_tree.fit(X=X2, y=y)
# fit OLS model
ols_model = sm.OLS(endog=y, exog=sm.add_constant(X2)).fit() 

OLS 摘要和决策树前两个层级的可视化展示了模型之间的显著差异（见 图 11.3）。OLS 模型提供了三个参数，分别为截距和两个特征，符合该模型对函数的线性假设。

相反，回归树图表显示了前两个层级的每个节点使用的特征和阈值来拆分数据（请注意特征可以重复使用），以及均方误差（MSE）、样本数量和基于这些训练样本的预测值的当前值。此外，请注意，与线性回归的 66 微秒相比，训练决策树需要 58 毫秒。虽然两种模型在只有两个特征时运行速度很快，但差异是 1,000 倍：

图 11.3：OLS 结果和回归树

树状图还突出显示了节点间样本分布的不均匀性，因为在第一个分割之后，样本数量在 545 到 55,000 之间变化。

为了进一步说明输入变量和输出之间的不同假设关系的功能形式，我们可以将当前回报预测可视化为特征空间的函数，即基于滞后回报值的值范围的函数。下图显示了线性回归（左侧面板）和回归树的当前月回报与一段时间前回报之间的关系：

图 11.4：线性回归和回归树的决策表面

左侧的线性回归模型结果强调了滞后和当前回报之间关系的线性性，而右侧的回归树图表则说明了特征空间的递归分区中编码的非线性关系。

构建分类树

分类树的工作方式与回归版本相同，只是结果的分类性质需要不同的方法来进行预测和衡量损失。虽然回归树使用相关训练样本的平均结果来预测分配给叶节点的观测值的响应，但是分类树使用模式，即相关区域内训练样本中最常见的类别。分类树还可以基于相对类频率生成概率预测。

如何优化节点纯度

在构建分类树时，我们也使用递归二元拆分，但是我们不是使用减少均方误差来评估决策规则的质量，而是可以使用分类错误率，它简单地是给定（叶）节点中不属于最常见类别的训练样本的比例。

然而，更倾向于使用基尼不纯度或交叉熵等替代度量方法，因为它们对节点纯度的敏感性比分类错误率更高，正如您在图 11.5中所见。节点纯度指的是节点中单个类别占优势的程度。一个只包含属于单个类别结果的样本的节点是纯净的，并且意味着在特征空间的这个特定区域的成功分类。

看看如何计算具有K类别 0,1,...,K-1（在二进制情况下为K=2）的分类结果的这些度量值。对于给定的节点m，让p[mk]为来自k^(th)类的样本比例：

下图显示了当类别比例均匀时（在二进制情况下为 0.5），基尼不纯度和交叉熵度量在[0,1]区间内达到最大值。当类别比例接近零或一时，这些度量值会下降，并且由于拆分而导致的子节点趋向纯净。与此同时，它们对节点不纯度的惩罚比分类错误率更高：

图 11.5：分类损失函数

请注意，与基尼测度相比，交叉熵的计算时间几乎要长 20 倍（详见笔记本中的详情）。

如何训练分类树

现在我们将使用 80%的样本进行训练，预测剩余的 20%来训练、可视化和评估一个最多进行五次连续拆分的分类树。为了简化说明，我们将采用一种捷径并使用内置的train_test_split，它不会防止前瞻偏差，就像我们在第六章——机器学习过程中介绍的自定义MultipleTimeSeriesCV迭代器一样，稍后我们将在本章中使用。

树的配置意味着最多有 2⁵=32 个叶节点，平衡情况下平均会包含超过 1,400 个训练样本。看一下下面的代码：

# randomize train-test split
X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=0.2, random_state=42)
# configure & train tree learner
clf = DecisionTreeClassifier(criterion='gini',
                            max_depth=5,
                            random_state=42)
clf.fit(X=X_train, y=y_train)
# Output:
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=5,
            max_features=None, max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, presort=False, random_state=42,
            splitter='best') 

训练模型后的输出显示了所有DecisionTreeClassifier的参数。我们将在超参数调整部分详细讨论这些。

可视化决策树

您可以使用 Graphviz 库（请参阅 GitHub 安装说明）来可视化树，因为 scikit-learn 可以输出使用该库使用的 DOT 语言描述的树的描述。您可以配置输出以包括特征和类标签，并限制级别的数量以使图表可读，如下所示：

dot_data = export_graphviz(classifier,
                           out_file=None, # save to file and convert to png
                           feature_names=X.columns,
                           class_names=['Down', 'Up'],
                           max_depth=3,
                           filled=True,
                           rounded=True,
                           special_characters=True)
graphviz.Source(dot_data) 

下图显示了模型如何使用不同的特征，并指示了连续和分类（虚拟）变量的分裂规则。在每个节点的标签值下，图表显示了来自每个类的样本数量，并在类标签下显示了最常见的类（在样本期间上涨的月份更多）：

图 11.6：分类树的可视化

评估决策树预测

为了评估我们第一个分类树的预测准确度，我们将使用测试集生成预测的类概率，如下所示：

# only keep probabilities for pos. class
y_score = classifier.predict_proba(X=X_test)[:, 1] 

.predict_proba()方法为每个类别生成一个概率。在二元类别中，这些概率是互补的并总和为 1，因此我们只需要正类的值。为了评估泛化误差，我们将使用基于接收器操作特征的曲线下面积，我们在第六章，机器学习过程中介绍过。结果表明，相对于随机预测的基准值 0.5，有了显著的改进（但请记住，这里的交叉验证方法不考虑数据的时间序列性）：

roc_auc_score(y_score=y_score, y_true=y_test)
0.6341 

过拟合和正则化

决策树在数据集相对于样本数量具有大量特征时很容易过拟合。正如前几章讨论的那样，过拟合会增加预测误差，因为模型不仅学习了训练数据中包含的信号，还学习了噪音。

有多种方法可以解决过拟合的风险，包括：

• 降维（请参阅第十三章，使用无监督学习的数据驱动风险因子和资产配置）通过用更少、更具信息性和更少噪声的特征表示现有特征来改善特征与样本的比率。

• 集成模型，例如随机森林，结合了多个树，同时随机化树的构建，我们将在本章的第二部分中看到。

决策树提供了几个正则化超参数来限制树的增长和相关复杂度。虽然每次分裂都会增加节点数，但也会减少每个节点可用于支持预测的样本数量。对于每个额外的层级，需要两倍数量的样本才能使新节点以相同的样本密度填充。

树修剪是减少树复杂性的另一种工具。它通过消除添加了很少价值但增加了模型方差的节点或整个树的部分来实现。例如，成本复杂度修剪从一个大树开始，并通过将节点替换为叶子来递归地减小其大小，本质上是反向运行树构建。各种步骤产生一系列树，然后可以使用交叉验证来选择理想的大小。

如何规范化决策树

以下表格列出了 scikit-learn 决策树实现中用于此目的的关键参数。在介绍了最重要的参数后，我们将说明如何使用交叉验证来优化超参数设置，以便在偏差-方差权衡和降低预测误差方面：

参数	描述	默认	选项
max_depth	最大级别数：分割节点直到达到max_depth。所有叶子都是纯的，或者包含的样本少于min_samples_split。	None	int
max_features	考虑用于分割的特征数量。	None	None：所有特征 int：＃ 特征float：分数auto，sqrt：sqrt（n_features）log2：log2（n_features）
max_leaf_nodes	分割节点直到创建这么多叶子。	None	None：无限 int
min_impurity_decrease	如果不纯度减少至少这个值，则分割节点。	0	float
min_samples_leaf	只有在左右分支的每个中至少有min_samples_leaf训练样本时，才会考虑分割。	1	int;float（作为百分比N）
min_samples_split	分割内部节点所需的最小样本数。	2	int；float（百分之N）
min_weight_fraction_leaf	在叶子节点上需要的所有样本权重的最小加权分数。除非在拟合方法中提供了sample_weight，否则样本具有相同的权重。	0

max_depth 参数对连续分割的次数施加了硬限制，并代表了限制树生长的最直接方式。

min_samples_split 和 min_samples_leaf 参数是另一种基于数据的方法来限制树的生长。与对连续分割的次数施加硬限制不同，这些参数控制进一步分割数据所需的最小样本数。后者保证了每个叶子的一定样本数量，而前者在分割导致非常不均匀的分布时可能会创建非常小的叶子。较小的参数值有助于过拟合，而较高的数字可能会阻止树学习数据中的信号。默认值通常相当低，您应该使用交叉验证来探索一系列潜在值。您还可以使用浮点数来表示百分比，而不是绝对数字。

scikit-learn 文档中包含有关如何在不同用例中使用各种参数的其他详细信息；有关更多信息，请参阅 GitHub 上链接的资源。

决策树修剪

递归二元分裂很可能会在训练集上产生良好的预测结果，但往往会导致数据过度拟合，产生较差的泛化性能。这是因为它导致了过于复杂的树，这在大量叶节点或特征空间的划分中反映出来。较少的分裂和叶节点意味着总体较小的树，并且通常会导致更好的预测性能，以及可解释性。

限制叶节点数量的一种方法是除非它们产生目标度量的显着改善，否则避免进一步分裂。然而，这种策略的缺点是，有时候，导致小幅改善的分裂在样本组成不断变化时会使后续更有价值的分裂变得更加困难。

相比之下，树修剪首先通过生长一个非常大的树，然后移除或修剪节点，以将大树减小为一个较少复杂且过度拟合的子树。成本复杂度修剪通过对向树模型添加叶节点增加惩罚并引入正则化参数来生成一系列子树，类似于套索和岭线性回归模型，调节惩罚的影响。应用于大树，增加的惩罚将自动产生一系列子树。可以使用正则化参数的交叉验证来识别最佳的修剪子树。

这种方法是在 scikit-learn 版本 0.22 中引入的；有关各种方法的工作原理和性能，请参见 Esposito 等人(1997)的调查。

超参数调整

决策树提供了一系列超参数来控制和调整训练结果。交叉验证是获得对泛化误差的无偏估计的最重要工具，反过来又允许在各种配置选项之间做出明智选择。scikit-learn 提供了几个工具来简化交叉验证大量参数设置的过程，即我们将在下一节中介绍的GridSearchCV便利类。学习曲线还允许进行诊断，评估收集额外数据以减少泛化误差的潜在好处。

使用自定义度量标准的 GridsearchCV

如第六章，机器学习流程中所强调的，scikit-learn 提供了一种定义多个超参数值范围的方法。它自动化了交叉验证各种参数值组合的过程，以确定最佳配置。让我们逐步了解自动调整模型的过程。

第一步是实例化一个模型对象，并定义一个字典，其中关键词命名超参数，值列出要测试的参数设置：

reg_tree = DecisionTreeRegressor(random_state=42)
param_grid = {'max_depth': [2, 3, 4, 5, 6, 7, 8, 10, 12, 15],
              'min_samples_leaf': [5, 25, 50, 100],
              'max_features': ['sqrt', 'auto']} 

然后，实例化GridSearchCV对象，提供估算器对象和参数网格，以及评分方法和交叉验证选择，传递给初始化方法。

我们将自定义的MultipleTimeSeriesSplit类设置为对模型进行 60 个月或 5 年的数据训练，并使用随后的 6 个月验证性能，重复此过程 10 次以覆盖 5 年的样本外期间：

cv = MultipleTimeSeriesCV(n_splits=10,
                          train_period_length=60,
                          test_period_length=6,
                          lookahead=1) 

我们使用roc_auc指标对分类器进行评分，并使用 scikit-learn 的make_scorer函数为回归模型定义自定义信息系数（IC）指标：

def rank_correl(y, y_pred):
    return spearmanr(y, y_pred)[0]
ic = make_scorer(rank_correl) 

我们可以使用n_jobs参数并通过设置refit=True自动获取使用最佳超参数的训练模型。

所有设置都已就绪后，我们可以像任何其他模型一样拟合GridSearchCV：

gridsearch_reg = GridSearchCV(estimator=reg_tree,
                          param_grid=param_grid,
                          scoring=ic,
                          n_jobs=-1,
                          cv=cv,  # custom MultipleTimeSeriesSplit
                          refit=True,
                          return_train_score=True)
gridsearch_reg.fit(X=X, y=y) 

训练过程为我们的GridSearchCV对象生成了一些新的属性，最重要的是有关最佳设置和最佳交叉验证分数的信息（现在使用适当的设置，避免了前瞻偏差）。

以下表格分别列出了最佳回归模型和分类模型的参数和分数。具有更浅的树和更加正则化的叶节点，回归树的 IC 为 0.083，而分类器的 AUC 得分为 0.525：

参数	回归	分类
max_depth	6	12
max_features	sqrt	sqrt
min_samples_leaf	50	5
分数	0.0829	0.5250

自动化非常方便，但我们也希望检查性能如何随不同参数值的变化而变化。完成此过程后，GridSearchCV对象提供了详细的交叉验证结果，以便我们可以获得更多见解。

如何检查树结构

该笔记本还说明了如何手动运行交叉验证以获取自定义树属性，例如与某些超参数设置相关联的总节点数或叶节点数。以下函数访问内部的.tree_ 属性以检索有关总节点计数以及其中多少个节点是叶节点的信息：

def get_leaves_count(tree):
    t = tree.tree_
    n = t.node_count
    leaves = len([i for i in range(t.node_count) if t.children_left[i]== -1])
    return leaves 

我们可以将此信息与训练和测试分数结合起来，以获取有关模型在整个交叉验证过程中行为的详细知识，如下所示：

train_scores, val_scores, leaves = {}, {}, {}
for max_depth in range(1, 26):
    print(max_depth, end=' ', flush=True)
    clf = DecisionTreeClassifier(criterion='gini', 
                                 max_depth=max_depth,
                                 min_samples_leaf=10,
                                 max_features='auto',
                                 random_state=42)
    train_scores[max_depth], val_scores[max_depth] = [], [] 
    leaves[max_depth] = []
    for train_idx, test_idx in cv.split(X):
        X_train, = X.iloc[train_idx], 
        y_train  = y_binary.iloc[train_ idx]
        X_test, y_test = X.iloc[test_idx], y_binary.iloc[test_idx]
        clf.fit(X=X_train, y=y_train)
        train_pred = clf.predict_proba(X=X_train)[:, 1]
        train_score = roc_auc_score(y_score=train_pred, y_true=y_train)
        train_scores[max_depth].append(train_score)
        test_pred = clf.predict_proba(X=X_test)[:, 1]
        val_score = roc_auc_score(y_score=test_pred, y_true=y_test)
        val_scores[max_depth].append(val_score)    
        leaves[max_depth].append(get_leaves_count(clf)) 

以下图表显示了叶节点数随树深度增加而增加的情况。由于每个交叉验证折叠的样本大小为 60 个月，每个样本约包含 500 个数据点，因此将min_samples_leaf限制为 10 个样本时，叶节点数限制在约 3,000 个：

图 11.7：分类树的可视化

比较回归和分类的性能

为了更仔细地观察模型的性能，我们将展示各种深度的交叉验证性能，同时保持产生最佳网格搜索结果的其他参数设置。图 11.8显示了训练和验证分数，并突出了更深的树的过拟合程度。这是因为训练分数稳步增加，而验证性能保持不变或下降。

请注意，对于分类树，网格搜索建议使用 12 个级别以获得最佳预测准确性。然而，图表显示较简单的树（具有三个或七个级别）的 AUC 得分相似。我们更倾向于一个更浅的树，它承诺具有可比较的泛化性能，同时减少了过拟合的风险：

图 11.8：两个模型的训练和验证分数

使用学习曲线诊断训练集规模

学习曲线是一种有用的工具，显示了随着训练样本数量增加，验证和训练分数如何演变。

学习曲线的目的是找出模型在训练过程中是否以及在多大程度上会受益于使用更多数据。它还有助于诊断模型的泛化误差更可能是由偏差还是方差驱动。

如果训练分数符合预期的性能，并且随着训练样本数量的增加，验证分数表现出显著改善，那么训练更长的回溯期或获取更多数据可能会增加价值。另一方面，如果尽管训练集大小增加，但验证分数和训练分数都收敛到类似较差的值，则错误更可能是由于偏差，并且额外的训练数据不太可能有所帮助。

以下图片显示了最佳回归和分类模型的学习曲线：

图 11.9：每个模型最佳版本的学习曲线

特别是对于回归模型，随着训练集规模的扩大，验证性能有所提高。这表明延长训练周期可能会产生更好的结果。你可以自己试一试，看看是否奏效！

从特征重要性中获得洞察力

决策树不仅可以可视化以检查给定特征的决策路径，还可以总结每个特征对模型学习以拟合训练数据的规则的贡献。

特征重要性捕捉到每个特征产生的分裂如何帮助优化模型用于评估分裂质量的度量，我们的情况下是基尼不纯度。特征的重要性是计算为这个度量的（标准化的）总减少，并考虑到受分裂影响的样本数量。因此，在树的较早节点使用的特征，其中节点往往包含更多的样本，通常被认为是更重要的。

图 11.10 显示了每个模型前 15 个特征的特征重要性的图表。注意特征的顺序如何与本节开头给出的互信息分数的单变量评估不同。显然，决策树捕获时间段与其他特征之间的相互依赖关系的能力可能改变每个特征的价值：

图 11.10：最佳回归和分类模型的特征重要性

决策树的优势和劣势

回归树和分类树在预测方面与我们在前几章中探讨过的线性模型有着非常不同的方法。如何决定哪种模型更适合当前的问题？考虑以下内容：

• 如果结果与特征之间的关系大致是线性的（或可以相应地转换），那么线性回归可能会胜过更复杂的方法，例如不利用这种线性结构的决策树。

• 如果关系呈现高度非线性和更复杂，决策树可能会胜过传统模型。请记住，关系的复杂性需要是系统的或“真实的”，而不是由噪声驱动，这会导致更复杂的模型过拟合。

决策树具有几个优点：

• 它们相当容易理解和解释，尤其是因为它们可以很容易地可视化，因此更容易让非技术人员理解。决策树也被称为白盒模型，因为它们在如何得出预测方面具有高度透明性。黑盒模型，例如集成和神经网络，可能会提供更好的预测精度，但是决策逻辑往往更难理解和解释。

• 决策树比那些对数据做出更强假设或对异常值更敏感且需要数据标准化的模型需要更少的数据准备（例如正则化回归）。

• 一些决策树实现可以处理分类输入，不需要创建虚拟变量（提高内存效率），并且可以处理缺失值，正如我们将在 第十二章，提升您的交易策略 中看到的，但这不适用于 scikit-learn。

• 预测速度很快，因为它与叶节点的数量呈对数关系（除非树变得极不平衡）。

• 可以使用统计测试验证模型并考虑其可靠性（有关更多细节，请参阅参考文献）。

决策树还具有几个关键劣势：

• 决策树内置了对训练集的过度拟合的倾向，并产生了高的泛化误差。解决这一弱点的关键步骤是修剪和正则化，使用限制树生长的早停准则，如本节所述。

• 决策树还对不平衡类权重敏感，可能会产生偏向某一类的树。一种选择是对少数类进行过采样或对更频繁出现的类进行欠采样。通常更好的选择是使用类权重并直接调整目标函数。

• 决策树的高方差与它们密切适应训练集的能力有关。因此，数据中的细微变化可能会导致树结构和模型预测的广泛波动。一个关键的预防机制是使用一组具有低偏差且产生不相关预测误差的随机决策树的集成。

• 决策树学习的贪婪方法优化了减少当前节点预测误差的局部标准，并不保证全局最优结果。同样，由随机树组成的集成有助于缓解这个问题。

我们现在将转向缓解决策树使用时固有的过拟合风险的集成方法。

随机森林——使树更可靠

决策树不仅因其透明度和解释性而有用。它们也是更强大的集成模型的基本构建模块，这些模型组合了许多个体树，同时随机变化其设计，以解决我们刚刚讨论的过拟合问题。

为什么集成模型表现更好

集成学习涉及将几个机器学习模型组合成一个新的模型，旨在比任何单个模型做出更好的预测。更具体地说，集成整合了几个基本估计器的预测，这些估计器使用一个或多个学习算法进行训练，以减少它们自己产生的泛化误差。

为了使集成学习达到这个目标，个体模型必须是：

• 准确: 胜过一个简单的基线（例如样本平均值或类比例）

• 独立: 预测是通过不同的方式产生的，以产生不同的误差

集成方法是最成功的机器学习算法之一，特别适用于标准的数值数据。大型集成在机器学习竞赛中非常成功，可能由许多不同的个体模型组成，这些模型已经手工组合或使用另一个机器学习算法组合起来。

将不同模型的预测组合起来存在几个缺点。这些包括降低了解释性以及训练、预测和模型维护的复杂性和成本。因此，在实践中（不考虑竞赛），从大规模集成中获得的准确度微小增益可能不值得额外的成本。

根据它们如何优化组成模型然后将结果集成为单个集成预测，通常可以区分两组集成方法：

• 平均方法独立地训练多个基本估计器，然后对它们的预测进行平均。如果基本模型没有偏差，并且产生不高度相关的不同预测误差，那么合并的预测可能具有较低的方差，更可靠。这类似于从具有不相关回报的资产构建组合以减少波动性而不牺牲回报。

• 提升方法相反，按顺序训练基本估计器，其特定目标是减少组合估计器的偏差。其动机是将几个弱模型组合成一个强大的集成。

在本章的剩余部分我们将专注于自动平均方法，而在第十二章《提升你的交易策略》中关注提升方法。

自助聚合

我们发现决策树很可能由于高方差而做出糟糕的预测，这意味着树结构对可用训练样本非常敏感。我们还看到，方差较低的模型，如线性回归，产生类似的估计值，尽管训练样本不同，只要给定特征数足够多。

对于给定一组独立观察值，每个观察值的方差为 ，样本均值的标准误差由 给出。换句话说，对更大的观察集进行平均会减小方差。因此，减少模型方差和泛化误差的一种自然方法是从总体中收集许多训练集，在每个数据集上训练不同的模型，然后对结果预测进行平均。

在实践中，我们通常没有许多不同的训练集。这就是装袋法的用武之地，即自助聚合。装袋法是一种通用方法，用于降低机器学习模型的方差，当应用于决策树时特别有用且受欢迎。

我们将首先解释这种技术如何缓解过拟合，然后展示如何将其应用于决策树。

如何通过装袋法降低模型的方差

装袋指的是自助采样的聚合，这些自助采样是带替换的随机采样。这样的随机样本与原始数据集具有相同数量的观察值，但由于替换可能包含重复观察值。

装袋法提高了预测准确性，但降低了模型的可解释性，因为不再能够可视化树以了解每个特征的重要性。作为一个集成算法，装袋方法在这些自助采样样本上训练给定数量的基本估计器，然后将它们的预测聚合成最终的集成预测。

装袋法通过以下方式降低基本估计器的方差，从而降低它们的泛化误差：

1. 随机化每棵树的生长方式

2. 对它们的预测进行平均

这通常是一种简单的方法来改进给定模型，而无需更改基础算法。这种技术对于具有低偏差和高方差的复杂模型特别有效，例如深层次的决策树，因为其目标是限制过拟合。相比之下，提升方法在弱模型（例如浅决策树）上效果最佳。

有几种 bagging 方法的不同之处在于它们应用于训练集的随机抽样过程：

• 粘贴从训练数据中无替换地随机抽取样本，而 bagging 则使用替换抽样。

• 随机子空间随机从特征（即列）中抽取样本，不进行替换。

• 随机补丁通过随机抽样观察值和特征来训练基本估算器。

Bagged 决策树

要将 bagging 应用于决策树，我们通过反复抽样替换来创建从训练数据中生成 bootstrap 样本。然后，我们在每个样本上训练一个决策树，并通过对不同树的预测进行平均来创建一个集成预测。您可以在笔记本bagged_decision_trees中找到此示例的代码，除非另有说明。

Bagged 决策树通常生长较大，即它们具有许多层和叶子节点，并且不进行修剪，以使每棵树的偏差低而方差高。然后，对其预测进行平均的效果旨在减少其方差。通过构建在 bootstrap 样本上训练的数百甚至数千棵树的集成，已经证明了 bagging 能够显著提高预测性能。

要说明 bagging 对回归树方差的影响，我们可以使用 scikit-learn 提供的BaggingRegressor元估计器。它基于指定抽样策略的参数来训练用户定义的基础估算器。

• max_samples和max_features分别控制从行和列中绘制的子集的大小。

• bootstrap和bootstrap_features确定是否使用或不使用替换来绘制这些样本。

以下示例使用指数函数生成单个DecisionTreeRegressor和由 10 棵树组成的BaggingRegressor集成的训练样本，每棵树都生长了 10 层深度。这两个模型都在随机样本上进行训练，并为添加了噪声的实际函数预测结果。

由于我们知道真实函数，因此我们可以将均方误差分解为偏差、方差和噪声，并根据以下分解比较这些组件的相对大小：

我们将分别从训练集和测试集中随机抽取 100 组样本，每组样本包括 250 个训练观察和 500 个测试观察，以训练每个模型并收集预测结果：

noise = .5  # noise relative to std(y)
noise = y.std() * noise
X_test = choice(x, size=test_size, replace=False)
max_depth = 10
n_estimators=10
tree = DecisionTreeRegressor(max_depth=max_depth)
bagged_tree = BaggingRegressor(base_estimator=tree, n_estimators=n_estimators)
learners = {'Decision Tree': tree, 'Bagging Regressor': bagged_tree}
predictions = {k: pd.DataFrame() for k, v in learners.items()}
for i in range(reps):
    X_train = choice(x, train_size)
    y_train = f(X_train) + normal(scale=noise, size=train_size)
    for label, learner in learners.items():
        learner.fit(X=X_train.reshape(-1, 1), y=y_train)
        preds = pd.DataFrame({i: learner.predict(X_test.reshape(-1, 1))},
                              index=X_test)
        predictions[label] = pd.concat([predictions[label], preds], axis=1) 

对于每个模型，在图 11.11中显示的绘图如下：

• 平均预测值和平均值周围两个标准偏差的带状区间（上部面板）。

• 根据真实函数的值（底部面板），进行偏差-方差-噪声分解。

我们发现，个体决策树预测的方差（左侧）几乎是基于自助样本的小集成的两倍高：

图 11.11：个体和袋装决策树的偏差-方差分解

请参阅笔记本 bagged_decision_trees 以获取实现细节。

如何构建随机森林

随机森林算法基于装袋引入的随机化来进一步减少方差并提高预测性能。

除了在自助训练数据上训练每个集合成员外，随机森林还会随机从模型中使用的特征中进行抽样（不重复）。根据实现方式，随机样本可以为每棵树或每次分裂进行抽取。因此，算法在学习新规则时面临不同选择，无论是在树的级别还是在每次分裂时。

特征的样本大小在回归和分类树之间有所不同：

• 对于分类，样本大小通常是特征数量的平方根。

• 对于回归，特征数量可以是从三分之一到所有特征，并且应基于交叉验证进行选择。

以下图表说明了随机森林如何随机化个体树的训练，然后将它们的预测汇总成集成预测：

图 11.12：随机森林如何生长个体树

除了训练观察结果外，随机化特征的目标是进一步使个体树的预测误差失相关。所有特征并非同等重要，一小部分高度相关的特征将在树构建过程中更频繁且更早地被选择，使得整个集成中的决策树更加相似。然而，个体树的泛化误差之间的相关性越小，整体方差就会减小得越多。

如何训练和调整随机森林

关键的配置参数包括介绍的各个决策树的超参数，在如何调整超参数部分。以下表格列出了两个RandomForest类的其他选项：

关键字	默认值	描述
bootstrap	TRUE	训练期间启用自助采样
n_estimators	10	森林中的树的数量
oob_score	FALSE	使用袋外样本来估计未见数据上的 R2

bootstrap参数激活了刚才描述的装袋算法。装袋又可以启用计算袋外得分（oob_score），该分数估计了未包括在用于训练给定树的自助样本中的样本的泛化准确性（请参阅袋外测试部分）。

参数n_estimators定义了作为森林一部分生长的树的数量。更大的森林表现更好，但也需要更长的时间来构建。重要的是随着基学习者数量的增加来监控交叉验证误差。目标是确定在训练额外的树的成本上升超过减少验证误差的好处的时候，或者后者开始再次增加的时候。

参数max_features控制在学习新决策规则并分割节点时可用的随机选择特征子集的大小。较低的值会降低树之间的相关性，从而减少集成的方差，但也可能增加偏差。正如本节开头所指出的，良好的起始值是回归问题的训练特征数量以及分类问题的这个数字的平方根，但会取决于特征之间的关系，并且应该使用交叉验证进行优化。

随机森林被设计为包含深度完全生长的树，可以使用max_depth=None和min_samples_split=2来创建。然而，这些值不一定是最优的，特别是对于具有许多样本和因此可能非常深的树的高维数据，可以变得非常计算和内存密集。

由 scikit-learn 提供的RandomForest类支持并行训练和预测，通过将n_jobs参数设置为要在不同核心上运行的k个作业数。值-1使用所有可用的核心。进程间通信的开销可能会限制速度提升呈线性增长，因此k个作业可能需要超过单个作业的 1/k时间。尽管如此，对于数据大、分裂评估变得昂贵的大型森林或深层个体树，速度提升通常是相当显著的。

始终应该使用交叉验证来识别最佳的参数配置。以下步骤说明了这个过程。本示例的代码在笔记本random_forest_tuning中。

我们将使用GridSearchCV来识别一组最佳参数，用于分类树的集成：

rf_clf = RandomForestClassifier(n_estimators=100,
                                criterion='gini',
                                max_depth=None,
                                min_samples_split=2,
                                min_samples_leaf=1,
                                min_weight_fraction_leaf=0.0,
                                max_features='auto',
                                max_leaf_nodes=None,
                                min_impurity_decrease=0.0,
                                min_impurity_split=None,
                                bootstrap=True, oob_score=False,
                                n_jobs=-1, random_state=42) 

我们使用与以前决策树示例中相同的 10 折自定义交叉验证，并使用关键配置设置的值填充参数网格：

cv = MultipleTimeSeriesCV(n_splits=10, train_period_length=60,
                          test_period_length=6, lookahead=1)
clf = RandomForestClassifier(random_state=42, n_jobs=-1)
param_grid = {'n_estimators': [50, 100, 250],
              'max_depth': [5, 15, None],
              'min_samples_leaf': [5, 25, 100]} 

使用上述内容配置GridSearchCV：

gridsearch_clf = GridSearchCV(estimator=clf,
                              param_grid=param_grid,
                              scoring='roc_auc',
                              n_jobs=-1,
                              cv=cv,
                              refit=True,
                              return_train_score=True,
                              verbose=1) 

我们像以前一样运行网格搜索，并找到最佳的回归和分类模型的以下结果。与分类器相比，随机森林回归模型对较浅的树表现更好，但在其他方面使用相同的设置：

Parameter	Regression	Classification
max_depth	5	15
min_samples_leaf	5	5
n_estimators	100	100
Score	0.0435	0.5205

然而，这两个模型都不及它们各自的单独决策树模型表现好，这凸显了在数据噪声较大且过拟合风险较高时，更复杂的模型并不一定优于更简单的方法。

随机森林的特征重要性

随机森林集成可能包含数百棵单独的树，但仍然可以从袋装模型中获得特征重要性的总体摘要指标。

对于给定的特征，重要性评分是由于对该特征进行拆分而导致的目标函数值的总减少量，并且在所有树上进行了平均。由于目标函数考虑了拆分影响的特征数量，因此在树的顶部附近使用的特征将获得较高的分数，因为更小数量的可用节点中包含了更多的观测数据。通过对以随机方式生长的许多树进行平均，特征重要性估计失去了一些方差，并变得更加准确。

分数以回归树的均方误差和分类树的基尼不纯度或熵来衡量。 scikit-learn 进一步规范化特征重要性，使其总和为 1。因此，计算得到的特征重要性也很受欢迎，作为特征选择的替代方法，与我们在第六章“机器学习流程”中看到的互信息度量相比（参见 sklearn.feature_selection 模块中的 SelectFromModel）。

图 11.13显示了两个模型的前 15 个特征的值。回归模型更依赖于时间段，而性能更好的决策树：

图 11.13：随机森林特征重要性

袋外测试

随机森林提供了内置交叉验证的好处，因为个别树是在训练数据的自助版本上训练的。因此，每棵树平均只使用了可用观测数据的三分之二。为了理解原因，考虑到自助样本的大小与原始样本相同，每个观测数据被抽取的概率是相同的，即 1/n。因此，根本不进入自助样本的概率是(1-1/n)n，它迅速收敛到 1/e，或者大约三分之一。

这剩余的三分之一未包含在训练集中的观测数据用于生成一个称为袋外 (OOB) 观测数据的装袋树，并可以作为验证集。正如交叉验证一样，我们为每棵树构建的不包含此观测数据的装袋样本预测响应，然后对每个 OOB 样本的单个集成预测进行平均（如果回归是目标，则为预测响应的平均值，如果分类是目标，则为多数投票或预测概率）。这些预测产生了对泛化误差的无偏估计，方便在训练期间计算。

由此产生的 OOB 误差是对该观测的泛化误差的有效估计。这是因为预测是使用在缺乏该观测时学到的决策规则产生的。一旦随机森林足够大，OOB 误差会接近于留一交叉验证误差。对于大型数据集，OOB 方法估计测试错误非常高效，而交叉验证可能计算成本很高。

然而，与交叉验证相同的警告适用：您需要注意避免“超前展望偏差”，如果 OOB 观测可以无序选择，则会发生这种偏差。在实践中，这使得在时间序列数据中使用 OOB 测试非常困难，因为验证集需要根据数据的顺序性选择。

随机森林的优缺点

袋装集成模型既有优点又有缺点。

随机森林的优点包括：

• 根据用例，随机森林可以与最佳监督学习算法表现相当。

• 随机森林提供可靠的特征重要性估计。

• 它们提供了测试错误的有效估计，而不会产生与交叉验证相关的重复模型训练成本。

另一方面，随机森林的缺点包括：

• 集成模型本质上比单个决策树不太可解释。

• 训练大量深树可能会产生高计算成本（但可以并行化）并使用大量内存。

• 预测速度较慢，这可能对需要低延迟的应用程序产生挑战。

现在让我们看看如何使用随机森林进行交易策略。

日本股票的多空信号

在第九章，用于波动率预测和统计套利的时间序列模型中，我们使用协整检验来识别具有长期均衡关系的股票对，其价格呈现共同趋势，并回归到该趋势。

在本章中，我们将利用机器学习模型的预测结果来识别可能上涨或下跌的资产，以便我们可以相应地进行市场中性的做多和做空操作。这种方法类似于我们最初的交易策略，该策略在第七章，线性模型 - 从风险因素到收益预测，以及第八章，ML4T 工作流程 - 从模型到策略回测中使用了线性回归。

我们将使用 LightGBM 软件包而不是 scikit-learn 随机森林实现，后者主要设计用于梯度增强。 LightGBM 的几个优点之一是其能够有效地将分类变量编码为数值特征，而不是使用独热编码（Fisher 1958）。我们将在下一章中提供更详细的介绍，但是代码示例应该很容易理解，因为逻辑与 scikit-learn 版本相似。

数据 - 日本股票

我们将设计一个关于日本股票宇宙的策略，使用由 Stooq 提供的数据，这是一家波兰数据提供商，目前提供了各种资产类别、市场和频率的有趣数据集，我们在 第九章，用于波动率预测和统计套利的时间序列模型 中也依赖了这些数据。

尽管关于数据的来源和质量几乎没有透明度，但它目前的免费使用具有强大的优势。换句话说，我们可以使用每日、每小时和每 5 分钟的股票、债券、商品和外汇数据进行实验，但应该对结果持谨慎态度。

本书 GitHub 存储库的 data 目录中的 create_datasets 笔记本包含了下载数据并将其存储为 HDF5 格式的说明。对于本例子，我们使用了 2010 年至 2019 年间约 3000 支日本股票的价格数据。最后的两年将作为样本外测试期，而之前的几年将作为我们的模型选择的交叉验证样本。

请参阅笔记本 japanese_equity_features 中的代码示例。我们移除了连续五个以上缺失值的股票代码，并仅保留了 1000 支交易量最高的股票。

特征 - 滞后回报和技术指标

我们会保持相对简单，将历史回报与 1、5、10、21 和 63 个交易日的几个技术指标结合起来，这些指标由 TA-Lib 提供（见第四章，金融特征工程 - 如何研究阿尔法因子）。

具体来说，我们为每支股票计算：

• 百分比价格振荡器（PPO）：是移动平均线收敛/发散（MACD）指标的归一化版本，用于衡量 14 天和 26 天指数移动平均线之间的差异，以捕捉不同资产间的动量差异。

• 归一化平均真实波幅（NATR）：以一种可比较各资产的方式来衡量价格波动。

• 相对强度指数（RSI）：另一个流行的动量指标（详见第四章，金融特征工程 - 如何研究阿尔法因子）。

• 布林带：移动平均线与移动标准差的比率，用于识别均值回归的机会。

我们还将包括标记以表示年、月和星期，以及根据每个交易日的六个间隔的最新回报对股票进行排名，排名从 1 到 20。

结果 - 不同视角的前瞻回报

为了测试随机森林在给定这些特征情况下的预测能力，我们生成了相同间隔的前瞻回报，最多到 21 个交易日（1 个月）。

历史和前瞻回报所暗示的滞后和领先导致了一些数据的丢失，该丢失随着投资视角的增加而增加。我们最终得到了 941 支股票的 18 个特征和 4 个结果的 230 万个观察值。

基于 LightGBM 的 ML4T 工作流程

现在我们将选择一个能产生可交易信号的随机森林模型。有几项研究已经成功地做到了，例如 Krauss、Do 和 Huck（2017）以及 Rasekhschaffe 和 Jones（2019）等，以及其中引用的资源。

我们将使用微软开源的快速且内存高效的 LightGBM 实现，它是最受欢迎的梯度提升方法，也是下一章节的主题，我们将更详细地了解各种 LightGBM 特性。

我们将从讨论关键实验设计决策开始，然后构建和评估一个预测模型，其信号将驱动我们将在最后一步设计和评估的交易策略。除非另有说明，请参考笔记本 random_forest_return_signals 中的代码示例。

从宇宙选择到超参数调整

要开发使用机器学习模型的交易策略，我们需要就模型的范围和设计做出几项决策，包括：

• 回溯期：用于训练的历史交易日数

• 前瞻期：预测收益的未来天数

• 测试期间：使用相同模型进行连续预测的天数

• 超参数：要评估的参数和配置

• 集成：是否依赖于单一模型或多个模型的组合

为了评估感兴趣的选项，我们还需要选择交叉验证的宇宙和时间段，以及一个样本外测试期和宇宙。更具体地说，我们在我们的日本股票样本的子集上交叉验证了截至 2017 年的几个选项。

一旦我们确定了一个模型，我们将定义交易规则，并在完整宇宙上过去两年的样本外使用我们模型的信号回测策略，以验证其性能。

对于时间序列交叉验证，我们将依赖于我们在第七章《线性模型 - 从风险因素到收益预测》中开发的 MultipleTimeSeriesCV，以对训练和测试期的长度进行参数化，同时避免前瞻性偏差。这个自定义的 CV 类允许我们：

• 对每个 ticker 进行连续包含 train_length 天的样本训练模型。

• 在后续包含 test_length 天和 lookahead 天数的期间内验证其性能，除了训练期外，以避免数据泄漏。

• 在每次将训练和验证期向前滚动 test_length 天时，重复 n_splits 次。

在本节中，我们将进行模型选择步骤，并在接下来的一节中进行策略回测。

对股票进行抽样以加快交叉验证的速度

训练随机森林所需的时间比线性回归要长得多，这取决于配置，树的数量和深度是主要驱动因素。

为了使我们的实验易于管理，我们将选择 2010-17 年期间交易量最大的 250 支股票，以评估不同结果和模型配置的性能，如下所示：

DATA_DIR = Path('..', 'data')
prices = (pd.read_hdf(DATA_DIR / 'assets.h5', 'stooq/jp/tse/stocks/prices')
          .loc[idx[:, '2010': '2017'], :])
dollar_vol = prices.close.mul(prices.volume)
dollar_vol_rank = dollar_vol.groupby(level='date').rank(ascending=False)
universe = dollar_vol_rank.groupby(level='symbol').mean().nsmallest(250).index 

定义回顾、展望和向前滚动期间

运行我们的策略需要以滚动方式训练模型，使用一定数量的交易日（回顾期）从我们的宇宙中学习模型参数，并预测一定数量的未来日子的结果。在我们的示例中，我们将考虑在每次迭代期间进行 63、126、252、756 和 1,260 个交易日的训练，并向前滚动和预测 5、21 或 63 天。

我们将把参数配对成一个列表，以便轻松迭代和可选的采样和/或洗牌，如下所示：

train_lengths = [1260, 756, 252, 126, 63]
test_lengths = [5, 21, 63]
test_params = list(product(train_lengths, test_lengths))
n = len(test_params)
test_param_sample = np.random.choice(list(range(n)), 
                                     size=int(n), 
                                     replace=False)
test_params = [test_params[i] for i in test_param_sample] 

使用 LightGBM 进行超参数调整

LightGBM 模型接受大量参数，如文档详细说明（参见lightgbm.readthedocs.io/和下一章）。对于我们的目的，我们只需要通过定义boosting_type启用随机森林算法，将bagging_freq设置为正数，并将objective设置为regression：

base_params = dict(boosting_type='rf',
                   objective='regression',
                   bagging_freq=1) 

接下来，我们选择最有可能影响预测准确性的超参数，即：

• 为模型生长的树的数量（num_boost_round）

• 用于装袋的行（bagging_fraction）和列（feature_fraction）的份额

• 在叶子中需要的最小样本数（min_data_in_leaf）以控制过拟合

LightGBM 的另一个好处是，我们可以评估训练模型的子集（或在一定数量的评估之后继续训练），这使我们能够在单个训练会话中测试多个num_iteration值。

或者，您可以启用early_stopping在验证集的损失指标不再改善时中断训练。但是，由于模型使用了在实际情况下不可用的结果信息，交叉验证性能估计将被偏向上。

我们将使用以下值来控制装袋方法和树生长的超参数：

bagging_fraction_opts = [.5, .75, .95]
feature_fraction_opts = [.75, .95]
min_data_in_leaf_opts = [250, 500, 1000]
cv_params = list(product(bagging_fraction_opts,
                         feature_fraction_opts,
                         min_data_in_leaf_opts))
n_cv_params = len(cv_params) 

在各种视野上交叉验证信号

要评估给定一组超参数的模型，我们将使用回顾、展望和向前滚动期间生成预测。

首先，我们将识别分类变量，因为 LightGBM 不需要独热编码；而是根据结果对类别进行排序，这样对回归树的效果更好，根据费舍尔（1958 年）的说法。我们将创建变量来识别不同的时期：

categoricals = ['year', 'weekday', 'month']
for feature in categoricals:
    data[feature] = pd.factorize(data[feature], sort=True)[0] 

为此，我们将创建二进制的 LightGBM 数据集，并使用给定的train_length和test_length配置MultipleTimeSeriesCV，这确定了我们 2 年验证期间的拆分数量：

for train_length, test_length in test_params:
    n_splits = int(2 * YEAR / test_length)
    cv = MultipleTimeSeriesCV(n_splits=n_splits,
                              test_period_length=test_length,
                              lookahead=lookahead,
                              train_period_length=train_length)
    label = label_dict[lookahead]
    outcome_data = data.loc[:, features + [label]].dropna()
    lgb_data = lgb.Dataset(data=outcome_data.drop(label, axis=1),
                           label=outcome_data[label],
                           categorical_feature=categoricals,
                           free_raw_data=False) 

接下来，我们采取以下步骤：

1. 为这次迭代选择超参数。

2. 将我们刚刚创建的二进制 LightGM 数据集切分为训练集和测试集。

3. 训练模型。

4. 为一系列num_iteration设置生成验证集的预测：

   for p, (bagging_fraction, feature_fraction, min_data_in_leaf) \
           in enumerate(cv_params_):
       params = base_params.copy()
       params.update(dict(bagging_fraction=bagging_fraction,
                          feature_fraction=feature_fraction,
                          min_data_in_leaf=min_data_in_leaf))
       start = time()
       cv_preds, nrounds = [], []
       for i, (train_idx, test_idx) in \
               enumerate(cv.split(X=outcome_data)):
           lgb_train = lgb_data.subset(train_idx.tolist()).construct()
           lgb_test = lgb_data.subset(test_idx.tolist()).construct()
           model = lgb.train(params=params,
                             train_set=lgb_train,
                             num_boost_round=num_boost_round,
                             verbose_eval=False)
           test_set = outcome_data.iloc[test_idx, :]
           X_test = test_set.loc[:, model.feature_name()]
           y_test = test_set.loc[:, label]
           y_pred = {str(n): model.predict(X_test, num_iteration=n)
                     for n in num_iterations}
           cv_preds.append(y_test.to_frame('y_test')
                           .assign(**y_pred).assign(i=i))
           nrounds.append(model.best_iteration) 
   

5. 为了评估验证性能，我们计算了完整预测集的 IC，以及每日的 IC，对一系列迭代次数进行了评估：

   df = [by_day.apply(lambda x: spearmanr(x.y_test,
                                          x[str(n)])[0]).to_frame(n)
         for n in num_iterations]
   ic_by_day = pd.concat(df, axis=1)
   daily_ic.append(ic_by_day.assign(bagging_fraction=bagging_fraction,
                                    feature_fraction=feature_fraction,
                                    min_data_in_leaf=min_data_in_leaf))
   cv_ic = [spearmanr(cv_preds.y_test, cv_preds[str(n)])[0]
            for n in num_iterations]
   ic.append([bagging_fraction, feature_fraction,
              min_data_in_leaf, lookahead] + cv_ic) 
   

现在，我们需要评估预测的信号内容以选择我们的交易策略的模型。

分析交叉验证性能

首先，我们将查看各种训练和测试窗口以及各种超参数设置在所有超参数设置中的 IC 分布。然后，我们将更详细地研究超参数设置对模型预测准确性的影响。

不同回溯、滚动向前和展望期的 IC

下图说明了每日平均 IC 的分布和分位数，针对四个预测时段和五个训练窗口，以及最佳的 21 天测试窗口。不幸的是，它并没有给出关于较短或较长窗口哪个效果更好的确切见解，而是说明了由于我们测试的模型配置范围以及由此产生的结果缺乏一致性而导致的数据噪声程度：

图 11.14：各种模型配置的每日平均信息系数的分布

使用 OLS 回归分析随机森林配置参数

为了更详细地了解我们实验的参数如何影响结果，我们可以对这些参数在每日平均 IC 上进行 OLS 回归分析。图 11.15显示了 1 天和 5 天展望期的系数和置信区间。

所有变量都经过 one-hot 编码，并且可以相对于每个变量的最小类别进行解释，该最小类别由常数捕获。结果在不同时段间有所不同；对于 1 天的预测，最长的训练期效果最好，但对于 5 天的预测效果最差，没有明显的模式。较长的训练似乎在一定程度上改善了 1 天模型，但对于 5 天模型来说这一点不太明确。唯一比较一致的结果似乎表明较低的装袋分数和较高的最小样本设置：

图 11.15：各种随机森林配置参数的 OLS 系数和置信区间

集成预测 - 使用 Alphalens 进行信号分析

最终，我们关心模型预测关于我们的投资范围和持有期的信号内容。为此，我们将使用 Alphalens 评估等权重组合投资于预测收益不同分位数产生的收益差异。

如第四章，金融特征工程 - 如何研究 Alpha 因子中讨论的那样，Alphalens 计算并可视化了各种摘要统计量，总结了 Alpha 因子的预测性能。笔记本alphalens_signals_quality说明了如何使用实用函数get_clean_factor_and_forward_returns将模型预测与价格数据以适当的格式结合起来。

为了解决 CV 预测中固有噪声的一些问题，我们根据其平均每日 IC 选择了前三个 1 天模型，并平均了它们的结果。

当我们将结果信号提供给 Alphalens 时，我们发现了以下结果，持续 1 天：

• 年化 alpha 为 0.081，beta 为 0.083

• 前五分位数收益之间的平均期间差距为 5.16 基点

以下图像可视化了按因子五分位数的平均期间收益和与每个分位数中的股票相关的累积每日正向收益：

图 11.16：Alphalens 因子信号评估

上述图像显示，1 天预测似乎包含了基于前五分位数的收益差异的有用交易信号。我们现在将继续开发并回测一种策略，该策略使用了在验证期间产生了这些结果的前十个 1 天前瞻模型的预测。

策略 - 使用 Zipline 进行回测

要设计并使用 Zipline 回测交易策略，我们需要为测试期间的我们的股票池生成预测，摄入日本股票数据并将信号加载到 Zipline 中，设置一个 pipeline，并定义重新平衡规则以相应地触发交易。

将日本股票纳入 Zipline

我们遵循第八章，ML4T 工作流程 - 从模型到策略回测中描述的流程，将我们的 Stooq 股票 OHLCV 数据转换为 Zipline bundle。目录custom_bundle包含了创建资产 ID 和元数据的预处理模块，定义了一个执行繁重任务的摄入函数，并使用扩展注册了 bundle。

该文件夹包含一个带有额外说明的 README。

运行样本内外策略的回测

笔记本 random_forest_return_signals 显示了如何选择产生最佳验证 IC 性能的超参数，并相应地生成预测。

我们将使用我们的 1 天模型预测，并应用一些简单的逻辑：对于具有最高正预期收益和最低负预期收益的 25 个资产，我们将建立多头和空头头寸。只要每边至少有 15 个候选资产，我们每天都会交易，并清除所有不在当前最佳预测中的头寸。

这次，我们还将包括每股 $0.05 的小额交易佣金，但不会使用滑点，因为我们正在以相对较小的资本基础交易日本最流动的股票。

结果 - 使用 pyfolio 进行评估

图 11.17中的左侧面板显示了该策略相对于日经 225 指数的样本内（2016-17 年）和样本外（2018-19 年）表现，而日经 225 指数在整个时期基本持平。

该策略在样本内年化基础上获得了 10.4%的收益率，在样本外年化基础上获得了 5.5%的收益率。

右侧面板显示了 3 个月滚动夏普比率，在样本内达到 0.96，在样本外达到 0.61：

图 11.17：Pyfolio 策略评估

总体表现统计数据突显了交易成本（每股 0.05 美分）后的 36.6%累计收益率，意味着样本外 Alpha 为 0.06，Beta 为 0.08（相对于日经 225 指数）。样本内最大回撤为 11.0%，样本外最大回撤为 8.7%。

	全部	样本内	样本外
# 月份	48	25	23
年化收益率	8.00%	10.40%	5.50%
累积收益率	36.60%	22.80%	11.20%
年化波动率	10.20%	10.90%	9.60%
夏普比率	0.8	0.96	0.61
Calmar 比率	0.72	0.94	0.63
稳定性	0.82	0.82	0.64
最大回撤	-11.00%	-11.00%	-8.70%
Sortino 比率	1.26	1.53	0.95
日度风险价值	-1.30%	-1.30%	-1.20%
Alpha	0.08	0.11	0.06
Beta	0.06	0.04	0.08

pyfolio tearsheets 包含许多有关暴露、风险配置和其他方面的额外细节。

摘要

在本章中，我们了解了一种能够捕捉非线性关系的新型模型，与我们迄今为止探索的经典线性模型形成对比。我们看到了决策树如何学习规则来将特征空间划分为产生预测的区域，从而将输入数据分段成特定区域。

决策树非常有用，因为它们提供了关于特征和目标变量之间关系的独特见解，我们看到了如何可视化树结构中编码的决策规则序列。

不幸的是，决策树容易过拟合。我们了解到集成模型和自举聚合方法成功地克服了决策树的一些缺点，并将它们作为更强大的复合模型的组成部分变得有用。

在下一章中，我们将探讨另一个集成模型，即提升（boosting），它已被认为是最重要的机器学习算法之一。