计算9.2.7的样本二阶原点矩方差怎么算三阶中心矩三阶偏度？
点击联系发帖人
时间：2023-09-25 14:05
在第7章”线性模型——从风险因子到资产收益预测“中，我们介绍了用于推理和预测的线性模型，从对产出有直接影响的横截面输入的同期关系的静态模型开始。我们介绍了普通最小二乘法（Ordinary Least Squares，OLS）学习算法，并看到它对一个正确指定的模型产生了无偏系数，其残差与输入变量不相关。加上残差具有恒定方差的假设，保证了普通最小二乘法在无偏估计中产生最小的平均平方预测误差。当我们了解到Fama-Macbeth回归是如何估计风险因子在不同时间和不同资产中的价值时，我们也遇到了既有横截面又有时间序列的面板数据。然而，不同时间的收益之间的关系通常是相当低的，所以这个程序在很大程度上可以忽略时间维度。此外，我们还介绍了正则化的岭回归和套索回归模型，这些模型产生有偏差的系数估计，但可以减少平均平方预测误差。这些预测模型从更动态的角度出发，将历史收益与其他输入结合起来，预测未来收益。在本章中，我们将建立动态线性模型来明确表示时间，并包括在特定间隔或滞后期观察的变量。时间序列数据的一个关键特征是它们的顺序：我们的数据不是像截面数据那样的单个观察的随机样本，而是一个随机过程的单一实现，我们不能重复。我们的目标是识别时间序列中的系统模式，帮助我们预测时间序列在未来的表现。更具体地说，我们将专注于从产出的历史序列中提取信号的模型，也可以选择其他同期或滞后的输入变量来预测产出的未来值。例如，我们可能会尝试利用过去的收益，结合基准或宏观经济变量的历史收益来预测一只股票的未来收益。在第4部分转向递归或卷积神经网络等非线性模型之前，我们将重点讨论线性时间序列模型。鉴于交易中固有的时间维度，时间序列模型非常流行。主要应用包括预测资产收益和波动率，以及识别资产价格序列的共同运动。随着越来越多的联网设备收集具有潜在信号内容的定期测量数据，时间序列数据可能会变得更加普遍。我们将首先介绍我们可以用来诊断时间序列特征和提取捕捉潜在模式的特征的工具。然后，我们将介绍如何诊断和实现时间序列的平稳性。接下来，我们将介绍单变量和多变量时间序列模型，并应用它们来预测宏观数据和波动模式。最后，我们将介绍协整的概念，以及如何应用它来制定一个配对交易策略。特别是，我们将涵盖以下主题：如何使用时间序列分析为建模过程做准备并提供信息估计和诊断单变量自回归和移动平均模型建立自回归条件异方差（Autoregressive Conditional Heteroskedasticity，ARCH）模型来预测波动率如何建立多元向量自回归模型使用协整法来制定配对交易策略你可以在GitHub资源库的相应目录中找到本章的代码样本和其他资源的链接。笔记本包括彩色版本的图片。关于从投资角度对本章主题的详尽介绍，请参见Tsay（2005）和Fabozzi, Focardi, and Kolm（2010）。9.1
诊断和特征提取工具时间序列是一个由离散的滞后期隔开的数值序列，这些间隔通常是均匀的（缺失的数值除外）。时间序列通常被建模为由一系列随机变量y(t_1),\cdots,y(t_T)组成的随机过程，每个时间点t_i,i=1,\cdots,T都有一个变量。单变量时间序列由每个时间点的单个值y组成，而多元时间序列由多个观测值组成，这些观测值可以用向量表示。期间数\Delta t=t_i-t_j，在不同的时间点t_i和t_j之间，被称为滞后期，每个时间序列都有T-1个不同的滞后期。正如在给定时间点上不同变量之间的关系是横截面模型的关键，被给定滞后期分开的数据点之间的关系是分析和利用时间序列中的模式的基础。对于横截面模型，我们区分输入变量和输出变量，或目标变量和预测变量，分别用y和x标记。在时间序列上下文中，结果y的部分或全部滞后期值y_{t_1},y_{t_2},\cdots,y_{t_T}在横截面上下文中扮演输入值或x值的角色。如果一个时间序列是一个独立同分布（Independent and Identically Distributed，IID）随机变量的序列，\epsilon_t，且其均值和方差都是有限的，那么该时间序列被称为白噪声。特别是，如果随机变量是正态分布，均值为零，方差为常数\sigma，则该序列被称为高斯白噪声。如果一个时间序列可以被写为过去扰动项\epsilon_t和序列的平均值\mu的加权和，那么它就是线性的。这些扰动项在这里被假定表示为白噪声：y_t=\mu+\sum_{i=0}^\infty a_i\epsilon_{t-i},\; a_0=1,\;\epsilon\sim i.i.d \\时间序列分析的一个关键目标是了解由系数a_i驱动的动态行为。时间序列的分析提供了针对这种类型的数据的方法，目的是提取有用的模式，反过来帮助我们建立预测模型。我们将为此介绍最重要的工具，包括分解为关键的系统要素，分析自相关，以及滚动窗口统计，如移动平均线。在本章的大部分例子中，我们将使用美联储提供的数据，你可以使用pandas-datareader访问这些数据，我们在第2章"市场和基本面数据——来源和技术 "中介绍了该数据。本节的代码例子可以在笔记本tsa_and_stationarity中找到。9.1.1
如何分解时间序列模式时间序列数据通常包含一个可以分解为几个组成部分的混合模式。特别是，一个时间序列往往结合了系统成分，如趋势、季节性和周期与非系统噪声。这些成分可以被建模为线性组合（例如，当波动不取决于序列的水平时）或以非线性、乘法的形式。基于模型的假设，它们也可以被自动拆分。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显示了附加组件的结果图表。假设趋势和季节性成分更具有确定性，可以进行简单的外推，那么剩余成分将是后续建模工作的重点：还有更复杂的基于模型的方法，例如，参见第6章：机器学习过程”，在Hyndman和Athanasopoulos(2018)。9.1.2
滚动窗口统计和移动平均鉴于时间序列数据的顺序排列，很自然地要计算熟悉的给定长度的周期的描述性统计。目标是检测该序列是否稳定或随时间变化，并获得一个平滑的表示，在过滤掉噪音的同时捕捉到系统方面。滚动窗口统计服务于这一过程：它们产生一个新的时间序列，其中每个数据点代表一个为原始数据的某个周期计算的汇总统计。移动平均数是最熟悉的例子。原始数据点可以以相等的权重进入计算，或者，例如，强调较近的数据点的权重。指数移动平均线递归地计算权重，对过去更远的数据点进行衰减。新的数据点通常是前面所有数据点的汇总，但它们也可以从周围的窗口中计算出来。pandas库包括滚动或扩展窗口，并允许各种权重分布。在第二步，你可以对窗口捕获的每一组数据进行计算。这些计算包括针对单个序列的内置函数，如平均值或总和，以及几个序列的相关性或协方差，还有用户定义的函数。例如，我们在第4章“金融特征工程——如何研究阿尔法因子 "和第7章”线性模型——从风险因子到收益预测“中使用了这种功能来设计特征。下一节中的移动平均数和指数平滑例子也将应用这些工具。早期的预测模型包括带有指数权重的移动平均模型，称为指数平滑模型。我们将再次遇到移动平均数作为线性时间序列的关键组成部分。依靠指数平滑方法的预测使用过去观测值的加权平均数，其中的权重随着观测值的老化而呈指数衰减。因此，一个较新的观测值会得到一个较高的相关权重。这些方法在没有非常复杂或突然的模式的时间序列中很受欢迎。9.1.3
如何测量自相关性自相关（也称为序列相关）使相关的概念适用于时间序列的上下文：正如相关系数衡量两个变量之间线性关系的强度一样，自相关系数\rho_k衡量由给定滞后期k的时间序列值之间的线性关系的程度：\rho_k=\frac{\sum_{t=k+1}^T(y_t-\bar{y})(y_{t-k}-\bar{y})}{\sum_{t=1}^T(y_t-\bar{y})^2} \\因此，我们可以为长度为T的时间序列中的每一个T-1滞后期计算一个自相关系数，自相关函数（Autocorrelation Function，ACF）计算相关系数作为滞后期的函数。滞后期大于1的自相关（即相隔超过一个时间步长的观测值之间）既反映了这些观测值之间的直接相关性，也反映了中间的数据点的间接影响。部分自相关消除了这种影响，只测量给定滞后期距离T的数据点之间的线性依赖。x_t和滞后期值x_{t-1},x_{t-2},\cdots,x_{T-2}作为特征（也被称为AR(T-1)模型，我们将在下一节讨论单变量时间序列模型）。如前所述，部分自相关函数（Partial Autocorrelation Function，PACF）提供了在较短的滞后期的相关性的影响被移除后所有的相关性。也有一些算法根据部分自相关函数和自相关函数之间的确切理论关系，从样本自相关中估计部分自相关。相关图只是一个连续滞后期的自相关函数或部分自相关k=0,1,\cdots,n的图。它使我们能够一目了然地检查各滞后期的相关结构（见图9.3的例子）。相关图的主要用途是检测去除确定性趋势或季节性后的任意自相关。自相关函数和部分自相关函数都是设计线性时间序列模型的关键诊断工具，我们将在下面关于时间序列转换的部分回顾自相关函数和部分自相关函数图的例子。9.2
如何诊断并达到平稳平稳的时间序列的统计属性，如平均数、方差或自相关，是独立于周期的，也就是说，它们不随时间变化。因此，平稳性意味着一个时间序列不存在趋势或季节性影响。此外，它要求描述性统计，如平均数或标准差，在计算不同的滚动窗口时，是恒定的或不随时间发生显著变化。一个平稳的时间序列会回复到它的平均值，方差差有一个恒定的振幅，而短期运动在统计意义上总是相似的。更正式地说，严格的平稳性要求任何时间序列观察的子集的联合分布在所有矩方面都与时间无关。因此，除了均值和方差之外，如偏度和峰度等更高的阶矩也需要保持不变，不管不同观测值之间的滞后期。在大多数应用中，比如本章中我们可以用来建立资产收益模型的大多数时间序列模型，我们将平稳性限制在一阶矩和二阶矩，这样时间序列就是具有恒定均值、方差和自相关的协方差平稳性。然而，在建立波动率模型时，我们放弃了这一假设，并明确地假设方差以可预测的方式随时间变化。请注意，我们特别考虑了不同滞后期的输出值之间存在的依赖关系，就像我们希望线性回归的输入数据与结果相关一样。平稳性意味着这些关系是稳定的。平稳性是经典统计模型的一个关键假设。下面两个小节介绍了可以帮助使时间序列平稳的变换，以及如何处理单位根引起的随机趋势这一特殊情况。9.2.1
变换时间序列以达到平稳为了满足许多时间序列模型的平稳性假设，我们需要对原始序列进行转换，通常分几个步骤。常见的转换包括（自然）对数，将指数增长模式转换为线性趋势并稳定方差。通缩意味着将一个时间序列除以另一个引起趋势行为的序列，例如，将一个名义序列除以一个价格指数，将其转换为一个实际的指标。如果一个序列恢复到一个稳定的长期线性趋势，那么它就是趋势平稳的。通常可以通过使用线性回归拟合一条趋势线并使用残差来使其成为平稳的。这意味着在线性回归模型中把时间指数作为一个独立变量，可能与记录或放空相结合。在许多情况下，去趋势化并不足以使序列平稳。相反，我们需要将原始数据转化为一系列的周期与周期和/或季节与季节的差异。换句话说，我们使用相邻的数据点或季节性滞后期的数值相互减去的结果。请注意，当这种差分被应用于对数转换的系列时，其结果代表了金融背景下的瞬时增长率或收益。如果一个单变量序列在经过d次差分后成为平稳的，我们就说它是d阶的积分，如果d=1，就说是简单的积分。这种行为是由于单位根引起的，我们接下来会解释。9.2.2
处理而不是如何处理单位根在确定将使时间序列平稳的转换构成了一个特殊的问题。我们将首先解释单位根的概念，然后讨论诊断测试和解决方案。9.2.2.1
单位根与随机游走时间序列通常被建模为以下自回归形式的随机过程，因此当前值是过去值的加权和，加上随机扰动：y_t=a_1y_{t-1}+a_2y_{t-2}+\cdots+a_py_{t-p}+\epsilon_t \\我们将在下一节关于单变量时间序列模型的内容中更详细地探讨这些模型，作为ARIMA模型的AR构建模块。这样一个过程有一个如下形式的特征方程：m^p-m^{p-1}a_1-m^{p-2}a_2-\cdots-a_p=0 \\如果这个多项式的（最多）p个根中的一个等于1，那么这个过程被称为有单位根。它将是非平稳的，但不一定有趋势。如果特征方程的其余根的绝对值小于1，那么该过程的1阶差分将是平稳的，**该过程是1阶积分或\mathbf{I}(1)**。如果其余的根的绝对值大于1，那么积分的阶数就会更高，需要进行额外的差分。在实践中，利率或资产价格的时间序列往往不是平稳的，因为没有一个价格水平可以让该序列还原。最突出的非平稳序列的例子是随机游走。给定一个价格的时间序列p_t，其起始价格p_0（例如，一只股票的IPO价格）和白噪声扰动\epsilon_t，那么随机游走统计符合以下自回归关系：p_t=p_{t-1}+\epsilon_t=\sum_{s=0}^t\epsilon_0+p_0 \\重复替换表明，当前的价格p_t是所有先前的扰动\epsilon_t和初始价格p_0的总和。如果该方程包括一个常数项，那么就说随机游走有漂移。因此，随机游走是一个自回归随机过程，其形式如下：y_t=a_1y_{t-1}+\epsilon_t,\;a_1=1 \\它的特征方程m-a_1=0，有单位根，既是非平稳的，又是1阶积分的。一方面，考虑到\epsilon的独立同分布性质，时间序列的方差等于t\sigma^2，这不是2阶平稳的，并意味着原则上序列可以随着时间的推移取任意值。另一方面，取1阶差分\Delta p_t=p_t-p_{t-1}，使得\Delta p_t=\epsilon_t在关于\epsilon的统计假设下是平稳的。具有单位根的非平稳序列的典型特征是长记忆：由于当前值是过去干扰的总和，大的创新比均值回复的平稳序列持续的时间长得多。9.2.2.2
如何诊断单位根统计单位根检验是客观地确定是否需要进行（额外的）差分的一种常见方法。这些是平稳性的统计假设检验，目的是确定是否需要进行差分。增强Dickey—Fuller检验（Augmented Dickey-Fuller Test，ADF test）评估了时间序列样本有单位根的无效假设和平稳性的替代假设。它将被差分的时间序列回归到时间趋势、1阶滞后期和所有滞后期差分上，并从滞后期时间序列值的系数中计算出检验统计量。 statsmodels库使其易于实现（参见笔记本tsa_and_stationarity）。从形式上看，一个时间序列的增强Dickey—Fuller检验运行线性回归，其中\alpha是一个常数，\beta是一个时间趋势的系数，p是指模型中使用的滞后期的数量：\Delta y_t=\alpha+\beta t+\gamma y_{t-1}+\delta_1\Delta y_{t-1}+\cdots+\delta_{p-1}\Delta y_{t-p+1}+\epsilon_t \\\alpha=\beta=0的约束意味着随机行走，而只有\beta=0意味着有漂移的随机行走。滞后期阶数通常使用第7章”线性模型——从风险因子到收益预测“中介绍的Akaike信息准则（AIC）和Bayesian信息准则（BIC）决定。增强Dickey—Fuller检验统计量使用样本系数 ，样本系数在单位根非平稳性的无效假设下，等于零，否则为负。本文的目的是说明，对于一个单整序列，其滞后期序列值不应提供预测滞后期差值以外和滞后差值以外的1阶段差值的有用信息。9.2.2.3
如何去除单位根并处理所产生的序列除了使用相邻数据点之间的差分来消除恒定的变化模式外，我们还可以应用季节性差分来消除季节性变化的模式。这涉及到在代表季节性模式长度的滞后期距离上取值的差分。对于月度数据，这通常涉及滞后期12的差分，而对于季度数据，这涉及时间价格4的差分，以消除季节性和线性趋势。确定正确的转换，特别是确定适当的差分数量和滞后期，并不总是一目了然。有人提出了一些启发式方法，可以归纳为以下几点：1期滞后期自相关接近零或负值，或自相关一般较小且无模式：不需要高阶差分。滞后期达10个以上的正自相关：该序列可能需要高阶差分。1期滞后期自相关<-0.5：该序列可能被过度差分了。可以用AR或MA项来纠正轻微的过度差分或不足（见下一节关于单变量时间序列模型）。一些作者建议将部分差分作为使单整序列平稳的更灵活的方法，并且可能比离散间隔的简单或季节性差分能够保留更多的信息或信号。例如，见Marcos Lopez de Prado（2018）中的第5章”投资组合优化和绩效评估“。9.2.3
实践中的时间序列变换图9.2中的图表显示了纳斯达克股票和工业生产指数30年来的时间序列的原始形式，以及分别应用对数和随后应用1阶和季节性差分（时间价格12）之后的转换版本。图表还显示了增强Dickey—Fuller检验的p值，它使我们能够拒绝两种情况下所有转换后的单位根非平稳性的假设：我们可以使用Q-Q图来进一步分析相关的时间序列特征，该图将时间序列观测值的分布量级与正态分布的量级以及基于自相关函数和部分自相关函数的相关图进行比较。对于图9.3中的纳斯达克图，我们可以看到，虽然没有趋势，但方差不是恒定的，而是在20世纪80年代末、2001年和2008年的市场动荡时期显示出集群式的高峰。Q-Q图突出了分布的肥尾，其极端值比正态分布所显示的更为频繁。自相关函数和部分自相关函数表现出相似的模式，在几个延迟时间的自相关似乎很显著：对于工业制造业生产的月度时间序列，我们可以看到2008年危机后的一个巨大的负离群点，以及Q-Q图中相应的倾斜（见图9.4）。自相关比纳斯达克的收益率要高得多，并且平稳地下降。部分自相关函数在滞后期1和13显示出明显的正自相关模式，在时间价格3和4显示出明显的负系数：9.3
单变量时间序列模型多重线性回归模型将感兴趣的变量表示为输入的线性组合，加上一个随机扰动。相比之下，单变量时间序列模型将时间序列的当前值与该序列的滞后值、当前噪声和可能的过去噪声项的线性组合联系起来。指数平滑模型是基于对数据中趋势和季节性的描述，而ARIMA模型旨在描述数据中的自相关。ARIMA(p, d, q)模型需要平稳性并利用两个构件：自回归（Autoregressive，AR）项由时间序列的p个滞后值组成移动平均（Moving Average，MA）项含由q个滞后扰动组成I代表差分，因为该模型可以通过微分d次来解释单位根非平稳性。术语自回归强调，ARIMA模型意味着时间序列对其自身值的回归。我们将介绍ARIMA构件，AR和MA模型，并解释如何在自回归移动平均（Autoregressive Moving-Average，ARMA）模型中结合它们，这些模型可以像ARIMA模型那样考虑序列整合，或者像AR(I)MAX模型那样包括外生变量。此外，我们将说明如何包括季节性的AR和MA项来扩展工具箱，使其也包括SARMAX模型。9.3.1
如何建立自回归模型p阶AR模型旨在捕捉不同滞后时间序列值之间的线性相关性，可以表示为：AR(p):\;y_t=\phi_0+\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\epsilon_t,\;\epsilon\sim i.i.d \\这与对y的滞后值进行多元线性回归非常相似。该模型具有如下特征方程：1-\phi_1x-\phi_2x^2-\cdots-\phi_px^p=0 \\这个p次多项式在x中的解的倒数是特征根，如果所有根的绝对值都小于1，AR(p)过程是平稳的，否则是不稳定的。对于平稳序列，多步预测将收敛于该序列的均值。我们可以用熟悉的最小二乘法估计模型参数，使用p+1,\cdots,T观测值，以确保每个滞后项和结果都有数据。9.3.1.1
如何识别滞后期的数量在实践中，挑战在于决定滞后项的适当阶数p。我们在 "如何测量自相关 "一节中讨论过的序列相关的时间序列分析工具，在做出这一决定时发挥了关键作用。更具体地说，对相关图的目视检查经常提供有用的线索：自相关函数估计了不同滞后期观测值之间的自相关性，这反过来又源于直接和间接的线性相关性。因此，如果一个k阶的AR模型是正确的模型，自相关函数将显示一个显著的序列相关性，直到滞后k，并且由于线性关系的间接影响造成的惯性，将延伸到随后的滞后期，直到它最终随着影响的减弱而减弱。反过来，局部自相关函数只测量相隔一定滞后期的观测值之间的直接线性关系，所以它不会反映出超过k滞后期的相关性。9.3.1.2
如何诊断模型拟合如果模型正确地捕捉了跨滞后期的线性相关性，那么残差应该类似于白噪声，而且自相关函数应该强调没有明显的自相关系数。除了残差图之外，Ljung-Box Q-统计量允许我们检验残差序列遵循白噪声的假设。无效假设是所有的m个序列相关系数都是零，而另一个假设是一些系数不是零。检验统计量是由不同滞后期k的样本自相关系数k计算出来的，并遵循X^2分布：Q(m)=T(T+2)\sum_{t=1}^m\frac{\rho_L^2}{T-l} \\正如我们将看到的，statsmodels库提供了关于不同滞后期的系数的显著性的信息，不显著的系数应该被删除。如果Q统计量拒绝了没有自相关的无效假设，你应该考虑额外的AR项。9.3.2
如何建立移动平均模型MA(q)模型在类似回归的模型中使用q个过去的扰动而不是时间序列的滞后值，如下所示：MA(q):\;y_t=c+\epsilon_t+\theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+\cdots+\theta_p\epsilon_{t-p}\;\epsilon\sim i.i.d \\由于我们没有观察到白噪声扰动值\epsilon_t，所以MA(q)不是像我们目前看到的那样的回归模型。MA(q)模型不是使用最小二乘法，而是使用最大似然法（Maximum Likelihood，MLE）进行估计，在序列开始时初始化或估计干扰值，然后递归和迭代地计算剩余部分。MA(q)模型的名字来自于将y_t的每个值表示为过去q个创新点的加权移动平均。换句话说，当前的估计值代表了相对于模型过去错误的修正。在MA(q)模型中使用移动平均不同于指数平滑或季节性时间序列成分的估计，因为MA(q)模型的目的是预测未来的数值，而不是去噪或估计过去数值的趋势周期。MA(q)过程总是平稳的，因为它们是白噪声变量的加权和，这些变量本身是平稳的。9.3.2.1
如何识别滞后期的数量一个由MA(q)过程产生的时间序列是由之前q模型预测的残差驱动的。因此，MA(q)过程的自相关函数将显示直到滞后q的值的显著系数，然后急剧下降，因为这是模型假设序列值产生的方式。请注意这与我们刚才描述的AR情况有什么不同，在AR情况下，部分自相关函数会显示类似的模式。9.3.2.2
AR模型和MA模型之间的关系一个AR(p)模型总是可以用重复替换的过程来表达。如在“如何处理单位根引起的随机趋势”一节中的随机游走例子。当MA(q)过程的系数满足一定的规模约束时，它也会变得可逆，可以表示为一个AR(\infty)过程（详见Tsay，2005）。9.3.3
如何构建ARIMA模型和扩展自回归移动平均数—ARIMA(p, d, q)—模型结合了AR(p)和MA(q)过程，利用这些构件的互补性，简化了模型的开发。它们采用更紧凑的形式，减少了参数的数量，从而降低了过拟合的风险。该模型还利用时间序列的第d次差分消除了单位根非平稳性。使用序列1阶差分的ARIMA(p,1,q)模型与ARMA(p,q)模型相同。用y'表示经过d次非季节性差分后的原始序列，ARIMA(p, d, q)模型很简单：\begin{align} ARIMA(p,d,q):\quad y'_t&=AR(P)+MA(q) \\ &=\phi_0+\phi_1y'_{t-1}+\cdots+\phi_py'_{t-p}+\epsilon_t+ \theta_1\epsilon_{t-1}+\cdots+\theta_q\epsilon_{t-p},\;\epsilon\sim i.i.d. \end{align} \\ARIMA模型也可以用极大似然估计法进行估计。根据实现的不同，高阶模型通常包含低阶模型。例如，到0.11版本为止，statsmodels库包括所有低阶p和q项，并且不允许删除低于最高值的滞后期的系数。在这种情况下，高阶模型总是能更好地拟合。要注意不要因为使用了太多的项而使你的模型与数据过度拟合。最近的版本，也就是本文写作时的0.11版本，增加了一个实验性的新ARIMA模型，具有更灵活的配置选项。9.3.3.1
如何对差分级数建模在使用数据时，也有一些设计单变量时间序列模型的指导准则：一个没有差分的模型假定原始序列是平稳的，包括均值回归的。它通常包括一个常数项，以允许非零均值的存在。一个有一阶差分的模型假设原始序列有一个恒定的趋势，因此应该包括一个常数项。二阶差分的模型假定原始序列有一个时间变化的趋势，且不应包括常数。9.3.3.2
如何识别AR和MA项的数量由于AR(p)和MA(q)项相互影响，自相关函数和部分自相关函数提供的信息不再可靠，只能作为一个起点。传统上，在选择模型设计时，Akaike信息准则和Bayesian信息准则被用来依赖样本内拟合。另外，我们可以依靠样本外测试来交叉验证多个参数选择。下面的总结为在单独考虑AR和MA模型的情况下如何选择模型阶数提供了一些指导：超过部分自相关函数中断的滞后期就是AR项的数量。如果差分序列的部分自相关函数显示出急剧下降和/或滞后1自相关系数为正，则增加一个或多个AR项。超过自相关函数中断的滞后期就是MA项的数量。如果差分序列的自相关函数显示出急剧的下降和/或滞后1自相关系数为负，考虑在模型中增加一个MA项。AR项和MA项可能会相互抵消，因此如果模型同时包含AR项和MA项，请尽量将它们的数量减少1，以避免过拟合，特别是当更复杂的模型需要超过10次迭代才能收敛时。如果AR系数之和接近1，表明模型的AR部分存在单位根，那么就消除一个AR项，并对模型进行一（更多）次差分。如果MA系数之和接近1，表明模型的MA部分存在单位根，则消除一个MA项，并将差分的阶数减少1。不稳定的长期预测表明在模型的AR或MA部分可能有一个单位根。9.3.3.3
增加特征—ARMAX带有外生输入的自回归移动平均模型（Autoregressive Moving-Average with Exogenous Inputs，ARMAX）模型在ARMA时间序列模型的右侧增加了输入变量或协变量（假设序列是平稳的，所以我们可以跳过差分）：\begin{align} ARMAX(p,d,q):\quad y_t&=\beta x_t+AR(P)+MA(q) \\ &=\beta x_t+\phi_0+\phi_1y_{t-1}+\cdots+\phi_py_{t-p}+\epsilon_t+ \theta_1\epsilon_{t-1}+\cdots+\theta_q\epsilon_{t-p},\;\epsilon\sim i.i.d. \end{align} \\这类似于线性回归模型，但解释起来相当困难。这是因为\beta对y_t的影响并不像线性回归中那样是x增加一个单位的影响。相反，方程右侧存在y_t的滞后值，这意味着系数只能根据响应变量的滞后值来解释，这是不直观的。9.3.3.4
增加季节差异—SARIMAX对于具有季节性效应的时间序列，我们可以引入AR和MA项来捕捉季节性的周期性。例如，当使用月度数据且季节性效应长度为1年时，季节性AR和MA项将反映这个特定的滞后长度。ARIMAX(p, d, q)模型就变成了SARIMAX(p, d, q)\times(P, D, Q)模型，写出来比较复杂，但statsmodels库文档（见GitHub上的链接）详细提供了这些信息。现在我们将使用宏数据建立一个季节性的ARMA模型来说明其实现。9.3.4
如何预测宏观基本面我们将为1988-2017年期间的工业生产时间序列的月度数据建立一个SARIMAX模型。正如关于分析工具的第一节中所说明的，数据 已经进行了对数转换，我们使用的是季节性（滞后12）差分。我们使用了10年的训练数据的滚动窗口对普通和传统的AR和MA参数范围的模型进行了估计，并评估超前一步预测的均方根误差（Root Mean Square Error，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)]我们还集成了Akaike信息准则和Bayesian信息准则，显示出非常高的秩相关系数0.94，Bayesian信息准则比Akaike信息准则更倾向于参数略少的模型。按均方根误差计算，最好的五个模型是：
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) \times (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())我们得到以下摘要：系数是显著的，而且Q统计量拒绝了进一步自相关的假设。相关图同样表明，我们已经成功地消除了该系列的自相关：9.3.5
如何使用时间序列模型预测波动率单变量时间序列模型在金融领域的一个特别重要的应用是对波动性的预测。这是因为随着时间的推移，波动率通常不是恒定的，而是聚集在一起的。方差的变化给使用假设平稳性的经典ARIMA模型进行时间序列预测带来了挑战。为了应对这一挑战，我们现在将对波动率进行建模，以便预测方差的变化。异方差是描述一个变量方差变化的技术术语。ARCH模型将误差项的方差表示为以前各期误差的函数。更具体地说，它假设误差方差遵循AR(p)模型。广义自回归条件异方差（Generalized Autoregressive Conditional Heteroskedasticity，GARCH）模型扩大了ARCH模型的范围，允许ARMA模型。时间序列的预测通常结合ARIMA模型的预期平均值和ARCH/GARCH模型的预期方差。2003年的诺贝尔经济学奖授予了Robert Engle和Clive Granger，以表彰他们开发了这一类模型。前者还管理着纽约大学斯特恩学院的波动率实验室（http://vlab.sten.nyu.edu），该实验室有许多关于我们将要讨论的模型的在线例子和工具。9.3.5.1
ARCH模型ARCH(p)模型只是一个应用于时间序列模型残差方差的AR(p)模型，它使t时刻的方差以方差的滞后观测值为条件。更具体地说，误差项\epsilon_t是线性模型（如ARIMA）在原始时间序列上的残差，分为随时间变化的标准差\sigma_t和扰动z_t，具体如下：\begin{align} ARCH(p):\quad var(x_t)&=\sigma_t^2 \\ &=\omega+\alpha_1\epsilon_{t-1}^2+\cdots+\alpha_p\epsilon_{t-p}^2 \\ \epsilon_t&=\sigma_tz_t \\ z_t&\sim i.i.d \end{align} \\ARCH(p)模型可以用普通最小二乘法来估计。Engle提出了一种使用拉格朗日乘数检验来识别适当的ARCH阶数的方法，它相当于对线性回归中所有系数为零的假设的F检验（见第7章“线性模型——从风险因子到收益预测”）。ARCH模型的一个主要优点是它产生的波动率估计值具有正的超额峰度，即相对于正态分布的肥尾，这反过来又与关于收益的经验观察相一致。缺点包括假设正负波动率冲击的影响相同，而资产价格的反应往往不同。它也不能解释波动率的变化，而且很可能过度预测波动率，因为它们对收益率序列的巨大、孤立的冲击反应缓慢。对于一个适当指定的ARCH模型，标准化的残差（除以模型对标准差周期的估计）应该类似于白噪声，可以进行Ljung-Box Q检验。9.3.5.2
推广ARCH——GARCH模型ARCH模型相对简单，但往往需要许多参数来捕捉资产收益序列的波动模式。GARCH模型适用于具有扰动\epsilon_t=r_t-\mu的对数收益序列r_t，该序列遵循GARCH(p,q)模型，如：\epsilon_t=\sigma_tz_t,\; \sigma_t^2=\omega+\sum_{i=1}^p\alpha_i\epsilon_{t-i}^2+\sum_{j=1}^q\beta_i\sigma_{t-j}^2\; z_t\sim i.i.d \\对于误差项的方差\epsilon_t，GARCH(p, q)模型假设为ARMA(p, q)模型。与ARCH模型类似，GARCH(1,1)过程的尾部分布比正态分布更重。该模型遇到了与ARCH模型相同的弱点。例如，它对正负冲击的反应相同。为了配置ARCH和GARCH模型的滞后阶数，使用训练好的时间序列的平方残差来预测原始序列的平均值。残差是以零为中心的，所以它们的平方也是方差。然后，检查平方残差的自相关函数和部分自相关函数图，确定时间序列方差的自相关模式。9.3.5.3
如何建立一个预测波动率的模型为资产收益序列建立波动率模型包括四个步骤：根据自相关函数和部分自相关函数所揭示的序列依赖性，为金融时间序列建立一个ARMA时间序列模型测试模型的残差是否具有ARCH/GARCH效应，同样依靠自相关函数和部分自相关函数对残差平方序列的测试如果序列相关效应显著，则指定一个波动率模型，并联合估计均值和波动率方程仔细检查拟合的模型，必要时对其进行改进在对收益率序列进行波动率预测时，可以对序列依赖性进行限制，这样就可以使用常数平均值来代替ARMA模型。arch库（见GitHub上的文档链接）提供了几个方法来估计波动率预测模型。你可以将预期均值建模为一个常数，也可以建模为AR(p)模型，如“如何建立自回归模型”一节中所讨论的，还可以建模为更近期的异质自回归过程（Heterogeneous Autoregressive Processes，HAR），它使用每日（1天）、每周（5天）和每月（22天）的滞后期来捕捉短期、中期和长期投资者的交易频率。均值模型可以用几个条件异方差模型共同定义和估计，这些模型除了ARCH和GARCH之外，还包括指数GARCH（Exponential GARCH，EGARCH）模型，它允许正负收益之间的不对称效应，以及异质ARCH（Heterogeneous，HARCH）模型，它是对HAR均值模型的补充。我们将使用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重新缩放的每日收益序列只表现出有限的自相关，但与平均值的平方偏差确实有大量的记忆，反映在缓慢衰减的自相关函数和部分自相关函数上，前两个延迟时间很高，只有在前六个滞后期后才截止： plot_correlogram(nasdaq_returns.sub(nasdaq_returns.mean()).pow(2), lags=120,
title='NASDAQ Daily Volatility')函数plot_correlogram的输出如下：因此，我们可以估计一个GARCH模型来捕捉过去波动率的线性关系。我们将使用10年的滚动窗口来估计GARCH(p, q)模型，p和q的范围是1-4，以产生一步样本外预测结果。然后我们比较预测波动率的均方根误差与实际收益率与其平均值的平方差，以确定最具预测性的模型。 我们使用winsorized数据来限制极端收益值反映在波动率的非常高的正偏度上的影响：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)模型实现了最低的均方根误差（与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())输出显示了最大化的对数似然，以及Akaike信息准则和Bayesian信息准则，这些标准在根据样本内表现选择模型时通常是最小化的（见第7章“线性模型——从风险因子到收益预测”）。它还显示了均值模型的结果，在本例子中，它只是一个常数估计，以及GARCH参数的常数omega，AR参数\alpha和MA参数\beta，所有这些都是统计显著的：现在让我们探索多时间序列模型和协整的概念，这将实现新的交易策略。9.4
多元时间序列模型多变量时间序列模型旨在同时捕捉多个时间序列的动态，并利用这些序列之间的依赖关系进行更可靠的预测。对这个问题最全面的介绍是Lütkepohl（2005）。9.4.1
等效系统方法单变量时间序列模型，如我们刚才讨论的ARMA方法，仅限于目标变量与其滞后值或滞后干扰和外生序列之间的统计关系，在ARMAX的情况下。相反，多变量时间序列模型还允许其他时间序列的滞后值影响目标变量。这种影响适用于所有序列，导致复杂的相互作用，如下图所示：除了潜在的更好的预测外，多变量时间序列还被用来深入了解跨序列的依赖性。例如，在经济学中，多变量时间序列被用来了解一个变量（如利率）的政策变化如何在不同时间范围内影响其他变量。多变量模型产生的脉冲响应函数就可以达到这个目的，它允许我们模拟一个变量如何对其他变量的突然变化做出反应。Granger因果关系的概念分析了一个变量是否有助于预测另一个变量（在最小二乘意义上）。此外，多变量时间序列模型允许对预测误差方差进行分解，以分析其他序列的贡献。9.4.2
向量自回归（VAR）模型我们将看到向量自回归VAR(p)模型如何通过建立一个由k个方程组成的系统，将AR(p)模型扩展到k个序列，其中每个方程包含所有k个序列的p个滞后值。在最简单的情况下，k=2的VAR(1)模型有如下形式：y_{1,t}=c_1+\alpha_{1,1}y_{1,t-1}+\alpha_{1,2}y_{2,t-2}+\epsilon_{1,t} \\ y_{2,t}=c_2+\alpha_{2,1}y_{2,t-1}+\alpha_{2,2}y_{2,t-2}+\epsilon_{2,t} \\这个模型可以用矩阵形式更简洁地表示：\left[\begin{aligned} y_{1,t} \\ y_{2,t} \end{aligned} \right] =
\left[\begin{aligned} c_1 \\ c_2 \end{aligned} \right] + \begin{bmatrix} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{bmatrix} \left[\begin{aligned} y_{1,t-1} \\ y_{2,t-2} \end{aligned} \right] + \left[\begin{aligned} \epsilon_{1,t} \\ \epsilon_{2,t} \end{aligned} \right] \\输出的滞后值的系数提供了关于序列本身的动态信息，而交叉变量系数则提供了关于序列之间相互作用的信息。这种表示方法可以扩展到k个p阶时间序列，如下所示：\begin{array}\pmb{y}_t \\ k\times 1\end{array}=\begin{array}\pmb{c} \\ k\times 1\end{array}+ \begin{array}\pmb{A}_t \\ k\times k\end{array}\begin{array}\pmb{y}_{t-1} \\ k\times 1\end{array}+ \cdots+\begin{array}\pmb{A}_p \\ k\times k\end{array}\begin{array}\pmb{y}_{t-p} \\ k\times 1\end{array} +\begin{array}\pmb{\epsilon}_t \\ k\times 1\end{array} \\VAR(p)模型也需要平稳性，以便单变量时间序列建模的初始步骤可以延续。首先，探索序列并确定必要的转换。然后，应用增强Dickey—Fuller检验来验证每个序列是否满足平稳性标准，否则应用进一步的转换。可以用以初始信息为条件的普通最小二乘法来估计，也可以用极大似然估计法来估计，对于正态分布的误差来说，与极大似然估计法是等价的，但在其他情况下则不然。如果k个序列中的一些或全部是单位根非平稳的，它们可能是协整的（见下一节）。这种将单位根概念扩展到多个时间序列的做法意味着两个或多个序列的线性组合是平稳的，因此可以实现均值回归。VAR模型无法在不进行差分的情况下处理这种情况，相应地使用向量误差修正模型（Vector Error Correction Model，VECM，Johansen和Juselius 1990）。我们将进一步探讨协整问题，因为如果存在协整并假定其持续存在，它可以被用于配对交易策略。滞后阶数的确定也是以每个序列的自相关函数和部分自相关函数为线索，但由于所有序列都具有相同的滞后阶数，因此滞后阶数的确定受到限制。在模型估计后，残差诊断也需要类似于白噪声的结果，模型选择可以使用样本内信息标准，如果目标是使用模型进行预测，则使用样本外预测性能来交叉验证其他模型设计。正如在单变量情况下所提到的，对原始时间序列的预测需要我们在训练模型之前将应用于使一个序列平稳的变换进行逆转。9.4.3
运用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这给我们留下了以下序列：为了限制输出的规模，我们将只是使用statsmodels库的VARMAX实现（允许可选的外生变量），用前480个观察值估计一个具有恒定趋势的VAR(1)模型： model = VARMAX(df_transformed.loc[:'2017'], order=(1,1), trend='c').fit(maxiter=1000)这会产生以下摘要：输出包含两个时间序列方程的系数，如前面的VAR(1)插图中所概述的。 statsmodels库提供了诊断图来检查残差是否符合白噪声假设。在这个简单的例子中，情况并非如此，因为方差似乎并不恒定（左上），而且分位数图显示了分布的差异，即肥尾（左下）：你可以像下面这样生成样本外预测：preds = model.predict(start=480, end=len(df_transformed)-1)下面可视化了实际值和预测值，显示了预测如何滞后于实际值，并且不能很好地捕捉非线性、样本外模式：9.5
协整—具有共同趋势的时间序列在上一节多元时间序列模型中，我们简要地提到了协整。现在让我们解释这个概念，以及如何在利用它进行统计套利交易策略之前更详细地诊断它的存在。我们已经看到了一个时间序列是如何具有单位根的，它创造了一个随机的趋势，并使时间序列具有高度的持久性。当我们用这样的集成时间序列的原始形式，而不是差分形式，作为线性回归模型的一个特征时，它与结果的关系往往会在统计上显得显著，尽管事实并非如此。这种现象被称为伪回归（详见第18章”金融时间序列和卫星图像的CNNs", Wooldridge, 2008）。因此，推荐的解决方案是对时间序列进行差分，使其在模型中使用之前成为平稳的。然而，当结果与一个或多个输入变量之间存在协整关系时，就有一个例外。为了理解协整的概念，首先让我们记住，回归模型的残差是输入和输出序列的线性组合。通常情况下，一个集成时间序列在一个或多个这样的序列上的回归残差也会产生非平稳的集成残差，因此表现为随机游走。然而，对于某些时间序列，情况并非如此：回归系数的线性组合的残差是平稳的，即使单个序列不是平稳的。这样的时间序列是协整的。一个非技术性的例子是，一个醉汉在他的狗（拴着的）陪同下随机散步。两条轨迹都是非平稳的，但却是协整的，因为狗偶尔会回到它的主人身边。在交易方面，套利约束意味着现货和期货价格之间的协整。换句话说，两个或多个协整序列的线性组合有一个稳定的平均值，这个线性组合会还原到这个平均值。这也适用于单个序列的高阶积分，线性组合降低了整体的积分阶。协整与相关不同：两个序列可以高度相关，但不一定是协整的。例如，如果两个增长的序列是彼此的恒定倍数，它们的相关性会很高，但任何线性组合也会增长，而不是恢复到一个稳定的平均值。协整是非常有用的：如果两个或更多的资产价格序列倾向于恢复到一个共同的平均值，我们可以利用偏离趋势的情况，因为它们应该意味着未来的价格移动方向相反。协整背后涉及的数学问题较多，所以我们将只关注实践方面；关于深入的处理，见Lütkepohl（2005）。在本节中，我们将讨论如何识别具有这种长期平稳关系的配对，估计任何不平衡修正的预期时间，以及如何利用这些工具来实施和回溯测试一个多空货币对交易策略。有两种方法可以检验协整关系：Engle—Granger两步法Johansen检验法我们将依次讨论每一种方法，然后展示它们如何帮助识别倾向于恢复到共同趋势的协整证券，我们可以利用这一事实来进行统计套利策略。9.5.1
Engle—Granger两步法Engle—Granger方法用于识别两个序列之间的协整关系。它涉及以下两个方面：用一个序列对另一个序列进行回归来估计平稳的长期关系对回归残差进行自相关函数单位根检验无效假设是残差有单位根，并且是集成的；如果我们能拒绝它，那么我们就假设残差是平稳的，因此，序列是协整的（Engle和Granger 1987）。这种方法的一个主要好处是，回归系数代表了使组合平稳的乘数，也就是说，均值可回归。不幸的是，检验结果会有所不同，这取决于我们认为哪个变量是独立的，因此，我们尝试两种方法，然后挑选p值较低的、检验统计量较较低的关系。另一个缺点是，这种检验仅限于成对的关系。更复杂的Johansen程序可以识别多达十几个时间序列之间的显著协整关系。9.5.2
Johansen似然比率检验而Johansen程序则是检验协整对VAR模型的限制，正如上一节所讨论的。更具体地说，在从一般的VAR(p)模型的两边减去目标向量后，我们得到误差修正模型（Error Correction Model，ECM）的表述：\Delta\pmb{y}_t=\pmb{c}+\pmb{\Pi}\pmb{y}_{t-1}+\pmb{\Gamma}_1\Delta y_{t-1}+\cdots +\pmb{\Gamma}_p\pmb{y}_{t-p}+\pmb{\epsilon}_t \\修改后的VAR(p)方程中只有一个行向量项（\pmb{y}_{t-1}）没有使用\Delta运算符表示差值。协整的性质取决于这一项系数矩阵的秩（Johansen 1991）。虽然这个方程在结构上看起来与自相关函数检验设置相似，但由于涉及多个序列，因此现在有几个共同趋势的潜在。为了确定协整关系的数量，Johansen检验从0（没有协整）开始，依次检验\Pi秩的增加。在下一节中，我们将探讨两个序列案例的应用。Gonzalo和Lee（1998）讨论了由于错误的模型动力学和其他实施方面的实际挑战，包括如何结合我们在下一节的样本统计套利策略中所依靠的两种检验程序。9.6
带有协整的统计套利统计套利是指采用某种统计模型或方法，利用看似相对错误的资产定价，同时保持一定程度的市场中性的策略。配对交易是一种概念简单的策略，至少从80年代中期开始就被算法交易者所采用（Gatev, Goetzmann，和Rouwenhorst 2006）。其目的是找到历史上价格一起变动的两种资产，跟踪价差（它们之间的价格差异），一旦价差扩大，买入跌破共同趋势的失败者，并做空胜利者。如果这种关系持续下去，随着价格的收敛和头寸的平仓，多头和/或空头将带来利润。这种方法扩展到多元环境，由多种证券组成篮子，并将一种资产与一个篮子或两个篮子彼此进行交易。在实践中，该策略需要两个步骤：形成阶段：识别具有长期均值回归关系的证券。理想情况下，价差应具有较高的方差，以便在可靠地回复到共同趋势的同时频繁地进行有利可图的交易。交易阶段：当价格变动导致价差背离和收敛时，触发进入和退出交易规则。在过去的几年里，在这个领域日益活跃的研究中，出现了几种形成和交易阶段的方法，涉及多个资产类别。下一小节概述了关键的差异，然后我们深入到一个应用实例中。9.6.1
如何选择和交易流动资产配对最近对成对交易策略的全面调查（Krauss 2017）确定了四种不同的方法，加上其他一些较新的方法，包括基于机器学习的预测：距离方法：这是最古老和研究最多的方法，用距离指标如相关性来确定候选配对，并使用非参数阈值如布林线来触发进入和退出交易。自Gatev等人（2006年）以来，其计算的简单性允许大规模的应用，并在整个市场和资产类别中表现出盈利能力。然而，最近的表现却有所下降。协整方法：如前所述，这种方法依赖于两个或多个变量之间长期关系的计量经济学模型。并允许进行统计检验，承诺比简单的距离度量更可靠。这一类的例子使用Engle-Granger和Johansen程序来识别证券配对和证券篮子，以及旨在捕捉概念的更简单的启发式方法（Vidyamurthy 2004）。交易规则通常类似于使用距离度量的简单阈值。时间序列方法：以交易阶段为重点，这类策略旨在将价差建模为均值回归的随机过程，并相应优化进入和退出规则（Elliott, Hoek，和Malcolm 2005）。它假定已经确定了有希望的配对。随机控制方法：与时间序列方法类似，目标是利用随机控制理论优化交易规则，找到价值和策略函数，以达到最佳的投资组合（Liu和Timmermann 2013）。我们将在第21章“合成时间序列数据的生成对抗网络”中讨论这种类型的方法。其他方法：除了基于主成分分析等无监督学习（见第13章“数据驱动的风险因子和无监督学习的资产配置”）和协同学等统计模型（Patton 2012）的配对识别，机器学习最近也开始流行，根据其相对价格或收益预测来识别配对（Huck 2019）。我们将在接下来的章节中介绍几种可用于此目的的机器学习算法，并说明相应的多变量交易策略。以上对各种方法的总结，仅仅是对配对交易策略设计所提供的灵活性的一瞥。除了关于配对选择和交易规则逻辑的更高层次的问题外，还有许多参数需要我们定义来实现。这些参数包括以下内容：筛选潜在配对或篮子的投资范围形成期的长度用于挑选可交易候选者的关系强度当价差波动时，对触发进入或退出交易或调整现有头寸的常用手段的偏离程度和趋同程度9.6.2
配对交易实践距离方法使用（归一化）资产价格或其收益的相关性来识别配对，它很简单，而且比协整测试的计算量要少很多。笔记本cointegration_test说明了这一点，它是以4年的每日数据为样本的150只股票为例：计算与ETF收益的相关性需要30ms，而一套协整检验（使用statsmodels库）需要18秒，慢了600倍。速度优势特别有价值。这是因为潜在配对的数量是两边要考虑的候选者数量的乘积，所以评估100只股票和100只ETF的组合需要比较10,000次测试（我们将在后面讨论多重测试偏差的挑战）。另一方面，距离指标不一定会选择最有利可图的配对：完美的联合移动会使相关性最大化，这反过来又会消除实际的交易机会。实证研究证实，协整配对的价差波动率几乎是距离配对价差波动率的两倍（Huck和Afawubo 2015）。为了平衡计算成本和所产生的配对质量之间的权衡，Krauss（2017）根据他的文献综述，推荐了一个结合两种方法的程序：选择分布稳定且漂移小的对，以减少候选集的数量检验剩下的具有最高价差的配对是否存在协整关系这个过程的目的是选择具有较低分化风险的协整对，同时确保更多的波动价差，反过来，产生更多的盈利机会。大量的测试会引入数据窥探的偏差，正如第6章“机器学习过程中”所讨论的：多次测试很可能会增加误判拒绝无协整性的无效假设的假阳性数量。虽然统计意义可能不是交易盈利的必要条件（Chan，2008），但对商品对的研究（Cummins和Bucca，2012）表明，根据Romano和Wolf（2010），控制家族误差率以提高测试的力量，可以带来更好的性能。在下一小节中，我们将更仔细地研究各种针对资产价格整合程度的启发式方法对协整检验结果的预测能力。该示例代码使用了在纽约证券交易所和纳斯达克交易所的172只股票和138只ETF的样本，并由Stooq提供了2010-2019年的每日数据。证券代表了其各自类别在样本期内最大的平均美元交易量；高度相关和固定的资产已被删除。关于如何获得数据的说明，请参见GitHub存储库数据文件夹中的笔记本create_datasets，关于相关代码和额外的预处理和探索细节的笔记本cointegration_tests。9.6.2.1
基于距离的启发式寻找协整配对compute_pair_metrics()计算了2010-14年和2015-19年超过23,000对股票和交易所交易基金（ETF）的以下距离指标：价差的漂移，定义为时间趋势对价差的线性回归价差的波动性归一化的价格序列之间以及它们的收益之间的相关性低漂移和波动性以及高相关性是协整的简单代理。为了评估这些启发式方法的预测能力，我们还使用statsmodels库对前面的配对进行Engle-Granger和Johansen协整检验。这发生在compute_pair_metrics()的后半部分的循环中。我们首先估计我们需要为Johansen检验指定的最佳滞后期数量。对于这两个检验，我们假设协整序列（价差）可能有一个不同于零的截距，但没有趋势。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的Johansen跟踪统计量与它们各自的临界值进行比较，得到Engle-Granger的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多配对关系，Johansen检验认为3.2%的关系是显著的，而Engle-Granger检验则认为6.5%。他们在366对关系（0.79%）上意见一致。9.6.2.2
启发式方法预测显著协整的效果如何？当我们把根据两种检验方法都是协整的序列的启发式方法的分布与不协整的其余序列进行比较时，波动率和漂移确实较低（以绝对值计算）。图9.14显示，两种相关度量的情况不太清楚：为了评估启发式方法的预测准确性，我们首先运行一个具有这些特征的逻辑回归模型来预测显著的协整关系。它取得了一个 曲线下面积（Area-Under-the-Curve，AUC）交叉验证分数为0.815；如果不考虑相关指标，它的分数仍为0.804。决策树的表现稍好，AUC=0.821，无论是否有相关性特征。特别是由于严重的类别不平衡，出现了大量的假阳性：正确识别366个协整对中的80%意味着超过16500个假阳性，但排除了近30000个候选者。更多细节见笔记本cointegration_ tests。关键的启示是，距离启发法可以帮助更有效地筛选一个大范围的数据，但这是以错过一些协整对为代价的，仍然需要大量的测试。9.6.3
准备策略回测在本节中，我们将对2017-2019年期间的股票和ETF的样本实施基于协整的统计套利策略。为了精简表述，有些方面被简化了。有关代码示例和其他细节，请参见笔记本 statistical_arbitrage_with_ cointegrated_pairs。我们首先生成并存储所有候选配对的协整检验以及由此产生的交易信号。然后，考虑到该过程的计算强度，我们对基于这些信号的策略进行回测。9.6.3.1
预计算协整检验首先，我们对23,000个潜在配对中的每一个进行2年回溯期的季度协整测试，然后，我们选择Johansen和Engle- Granger检验都一致的配对进行交易。我们应该排除那些在回溯期内静止的资产，但是我们排除了在整个时期内静止的资产，所以我们跳过这一步来简化它。这一程序遵循之前概述的步骤；详情请见笔记本。图9.15显示了为交易而选择的两个不同配对的原始股票和ETF序列；注意在样本期间明显存在着共同的趋势。9.6.3.2
进行入市和出市交易现在，我们可以根据滚动对冲比率来计算每个候选配对的价差。我们还计算了布林带，因为我们将把价差远离其移动平均线大于两个滚动标准差的移动视为多头和空头进入信号，而移动平均线的反向交叉则是退出信号。用卡尔曼率滤波平滑价格为此，我们首先应用滚动卡尔曼滤波器（Kalman Filter，KF）来去除一些噪音，如第4章“金融特征工程——如何研究阿尔法因素子”所展示的：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)使用卡尔曼滤波器计算滚动套期保值率为了获得动态套期保值比率，我们将卡尔曼滤波器应用于滚动线性回归，如下所示：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估算均值回归的半衰期如果我们将利差视为连续时间内的均值回归随机过程，我们可以将其建模为Ornstein-Uhlenbeck过程。这种观点的好处是，我们可以获得均值回归的半衰期公式，作为偏差后价差再次收敛所需时间的近似值（详见Chan 2013，第2章“市场和基本面数据——来源和技术”）：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分数来推导出交易信号：如果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_trades9.6.4
使用backtrader库对策略进行回测现在，我们已经准备好在我们的回测平台上制定我们的策略，执行它，并评估结果。为此，我们需要跟踪我们的配对，除了个人投资组合头寸外，监测活跃和不活跃配对的分布，以应用我们的交易规则。9.6.4.1
用一个自定义的数据类跟踪配对为了说明活跃的配对，我们定义了一个数据类（在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))9.6.4.2
运行和评估该战略关键的执行方面包括：每天从已触发退出规则或超过特定负收益的配对中退出为那些价差触发了进入信号的配对对开立新的多头和空头头寸此外，我们调整头寸以考虑到不同数量的配对策略本身的代码占用了太多空间，无法在此显示；请看笔记本上的pair_trading_backtest了解详情。图9-16表明，至少在2017-2019年期间，这种简化的策略曾有过辉煌时刻（请注意，我们利用了一些前瞻性偏差，忽略了交易成本）。在这些宽松的假设下，它在期初和期末的表现低于标准普尔500指数，其他方面则基本保持一致（左图）。它的阿尔法值为0.08，贝塔值为-0.14（右图），平均夏普比率为0.75，Sortino比率为1.05（中央图）：虽然我们应该对这些业绩指标持谨慎态度，但该策略展示了基于协整的统计套利的剖析，其形式是配对交易。让我们来看看你可以采取的几个步骤，在这个框架的基础上产生更好的业绩。9.6.5
扩展—如何做得更好协整是一个非常有用的概念，可以识别倾向于一致行动的股票对或股票组。与协整的统计复杂性相比，我们使用的是非常简单和静态的交易规则；按季度计算也扭曲了策略，正如多头和空头持有的模式所示（见笔记本）。要想获得成功，你至少需要筛选一个更大的集合，并优化几个参数，包括交易规则。此外，风险管理应考虑到当某些资产相对频繁地出现在一个交易对的同一侧时产生的集中头寸。你也可以用篮子来操作，而不是单独的配对；然而，为了解决越来越多的候选数量，你可能需要限制篮子的组成。正如在“配对交易—利用协整进行统计套利”一节中提到的，还有一些旨在预测价格变动的替代方法。在接下来的章节中，我们将探讨各种机器学习模型，旨在预测特定投资领域的绝对规模或方向。使用这些预测作为多头和空头入场信号，是我们在本节中研究的配对交易框架的自然延伸或替代方案。9.7
总结在本章中，我们探讨了单个序列的单变量情况下的线性时间序列模型，以及几个相互作用的序列的多变量模型。我们遇到了预测宏观基本面的应用，预测资产或投资组合波动性的模型，以及在风险管理中广泛使用的多变量VAR模型，以及捕捉多个宏观序列的动态变化。我们还研究了协整的概念，它是流行的配对交易策略的基础。与第7章"线性模型——从风险因子到收益预测“类似，我们看到了线性模型如何强加了大量结构，也就是说，它们做出了强有力的假设，有可能需要进行转换和大量测试来验证这些假设是否得到满足。如果符合，模型的训练和解释就很简单，而且这些模型可以提供一个良好的基线，更复杂的模型可以在此基础上进行改进。在接下来的两章中，我们将看到这方面的两个例子，即随机森林和梯度提升模型，我们将在第4部分遇到更多的例子，这部分是关于深度学习的。}
诊断和特征提取工具时间序列是一个由离散的滞后期隔开的数值序列，这些间隔通常是均匀的（缺失的数值除外）。时间序列通常被建模为由一系列随机变量y(t_1),\cdots,y(t_T)组成的随机过程，每个时间点t_i,i=1,\cdots,T都有一个变量。单变量时间序列由每个时间点的单个值y组成，而多元时间序列由多个观测值组成，这些观测值可以用向量表示。期间数\Delta t=t_i-t_j，在不同的时间点t_i和t_j之间，被称为滞后期，每个时间序列都有T-1个不同的滞后期。正如在给定时间点上不同变量之间的关系是横截面模型的关键，被给定滞后期分开的数据点之间的关系是分析和利用时间序列中的模式的基础。对于横截面模型，我们区分输入变量和输出变量，或目标变量和预测变量，分别用y和x标记。在时间序列上下文中，结果y的部分或全部滞后期值y_{t_1},y_{t_2},\cdots,y_{t_T}在横截面上下文中扮演输入值或x值的角色。如果一个时间序列是一个独立同分布（Independent and Identically Distributed，IID）随机变量的序列，\epsilon_t，且其均值和方差都是有限的，那么该时间序列被称为白噪声。特别是，如果随机变量是正态分布，均值为零，方差为常数\sigma，则该序列被称为高斯白噪声。如果一个时间序列可以被写为过去扰动项\epsilon_t和序列的平均值\mu的加权和，那么它就是线性的。这些扰动项在这里被假定表示为白噪声：y_t=\mu+\sum_{i=0}^\infty a_i\epsilon_{t-i},\; a_0=1,\;\epsilon\sim i.i.d \\时间序列分析的一个关键目标是了解由系数a_i驱动的动态行为。时间序列的分析提供了针对这种类型的数据的方法，目的是提取有用的模式，反过来帮助我们建立预测模型。我们将为此介绍最重要的工具，包括分解为关键的系统要素，分析自相关，以及滚动窗口统计，如移动平均线。在本章的大部分例子中，我们将使用美联储提供的数据，你可以使用pandas-datareader访问这些数据，我们在第2章"市场和基本面数据——来源和技术 "中介绍了该数据。本节的代码例子可以在笔记本tsa_and_stationarity中找到。9.1.1
如何分解时间序列模式时间序列数据通常包含一个可以分解为几个组成部分的混合模式。特别是，一个时间序列往往结合了系统成分，如趋势、季节性和周期与非系统噪声。这些成分可以被建模为线性组合（例如，当波动不取决于序列的水平时）或以非线性、乘法的形式。基于模型的假设，它们也可以被自动拆分。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显示了附加组件的结果图表。假设趋势和季节性成分更具有确定性，可以进行简单的外推，那么剩余成分将是后续建模工作的重点：还有更复杂的基于模型的方法，例如，参见第6章：机器学习过程”，在Hyndman和Athanasopoulos(2018)。9.1.2
滚动窗口统计和移动平均鉴于时间序列数据的顺序排列，很自然地要计算熟悉的给定长度的周期的描述性统计。目标是检测该序列是否稳定或随时间变化，并获得一个平滑的表示，在过滤掉噪音的同时捕捉到系统方面。滚动窗口统计服务于这一过程：它们产生一个新的时间序列，其中每个数据点代表一个为原始数据的某个周期计算的汇总统计。移动平均数是最熟悉的例子。原始数据点可以以相等的权重进入计算，或者，例如，强调较近的数据点的权重。指数移动平均线递归地计算权重，对过去更远的数据点进行衰减。新的数据点通常是前面所有数据点的汇总，但它们也可以从周围的窗口中计算出来。pandas库包括滚动或扩展窗口，并允许各种权重分布。在第二步，你可以对窗口捕获的每一组数据进行计算。这些计算包括针对单个序列的内置函数，如平均值或总和，以及几个序列的相关性或协方差，还有用户定义的函数。例如，我们在第4章“金融特征工程——如何研究阿尔法因子 "和第7章”线性模型——从风险因子到收益预测“中使用了这种功能来设计特征。下一节中的移动平均数和指数平滑例子也将应用这些工具。早期的预测模型包括带有指数权重的移动平均模型，称为指数平滑模型。我们将再次遇到移动平均数作为线性时间序列的关键组成部分。依靠指数平滑方法的预测使用过去观测值的加权平均数，其中的权重随着观测值的老化而呈指数衰减。因此，一个较新的观测值会得到一个较高的相关权重。这些方法在没有非常复杂或突然的模式的时间序列中很受欢迎。9.1.3
如何测量自相关性自相关（也称为序列相关）使相关的概念适用于时间序列的上下文：正如相关系数衡量两个变量之间线性关系的强度一样，自相关系数\rho_k衡量由给定滞后期k的时间序列值之间的线性关系的程度：\rho_k=\frac{\sum_{t=k+1}^T(y_t-\bar{y})(y_{t-k}-\bar{y})}{\sum_{t=1}^T(y_t-\bar{y})^2} \\因此，我们可以为长度为T的时间序列中的每一个T-1滞后期计算一个自相关系数，自相关函数（Autocorrelation Function，ACF）计算相关系数作为滞后期的函数。滞后期大于1的自相关（即相隔超过一个时间步长的观测值之间）既反映了这些观测值之间的直接相关性，也反映了中间的数据点的间接影响。部分自相关消除了这种影响，只测量给定滞后期距离T的数据点之间的线性依赖。x_t和滞后期值x_{t-1},x_{t-2},\cdots,x_{T-2}作为特征（也被称为AR(T-1)模型，我们将在下一节讨论单变量时间序列模型）。如前所述，部分自相关函数（Partial Autocorrelation Function，PACF）提供了在较短的滞后期的相关性的影响被移除后所有的相关性。也有一些算法根据部分自相关函数和自相关函数之间的确切理论关系，从样本自相关中估计部分自相关。相关图只是一个连续滞后期的自相关函数或部分自相关k=0,1,\cdots,n的图。它使我们能够一目了然地检查各滞后期的相关结构（见图9.3的例子）。相关图的主要用途是检测去除确定性趋势或季节性后的任意自相关。自相关函数和部分自相关函数都是设计线性时间序列模型的关键诊断工具，我们将在下面关于时间序列转换的部分回顾自相关函数和部分自相关函数图的例子。9.2
如何诊断并达到平稳平稳的时间序列的统计属性，如平均数、方差或自相关，是独立于周期的，也就是说，它们不随时间变化。因此，平稳性意味着一个时间序列不存在趋势或季节性影响。此外，它要求描述性统计，如平均数或标准差，在计算不同的滚动窗口时，是恒定的或不随时间发生显著变化。一个平稳的时间序列会回复到它的平均值，方差差有一个恒定的振幅，而短期运动在统计意义上总是相似的。更正式地说，严格的平稳性要求任何时间序列观察的子集的联合分布在所有矩方面都与时间无关。因此，除了均值和方差之外，如偏度和峰度等更高的阶矩也需要保持不变，不管不同观测值之间的滞后期。在大多数应用中，比如本章中我们可以用来建立资产收益模型的大多数时间序列模型，我们将平稳性限制在一阶矩和二阶矩，这样时间序列就是具有恒定均值、方差和自相关的协方差平稳性。然而，在建立波动率模型时，我们放弃了这一假设，并明确地假设方差以可预测的方式随时间变化。请注意，我们特别考虑了不同滞后期的输出值之间存在的依赖关系，就像我们希望线性回归的输入数据与结果相关一样。平稳性意味着这些关系是稳定的。平稳性是经典统计模型的一个关键假设。下面两个小节介绍了可以帮助使时间序列平稳的变换，以及如何处理单位根引起的随机趋势这一特殊情况。9.2.1
变换时间序列以达到平稳为了满足许多时间序列模型的平稳性假设，我们需要对原始序列进行转换，通常分几个步骤。常见的转换包括（自然）对数，将指数增长模式转换为线性趋势并稳定方差。通缩意味着将一个时间序列除以另一个引起趋势行为的序列，例如，将一个名义序列除以一个价格指数，将其转换为一个实际的指标。如果一个序列恢复到一个稳定的长期线性趋势，那么它就是趋势平稳的。通常可以通过使用线性回归拟合一条趋势线并使用残差来使其成为平稳的。这意味着在线性回归模型中把时间指数作为一个独立变量，可能与记录或放空相结合。在许多情况下，去趋势化并不足以使序列平稳。相反，我们需要将原始数据转化为一系列的周期与周期和/或季节与季节的差异。换句话说，我们使用相邻的数据点或季节性滞后期的数值相互减去的结果。请注意，当这种差分被应用于对数转换的系列时，其结果代表了金融背景下的瞬时增长率或收益。如果一个单变量序列在经过d次差分后成为平稳的，我们就说它是d阶的积分，如果d=1，就说是简单的积分。这种行为是由于单位根引起的，我们接下来会解释。9.2.2
处理而不是如何处理单位根在确定将使时间序列平稳的转换构成了一个特殊的问题。我们将首先解释单位根的概念，然后讨论诊断测试和解决方案。9.2.2.1
单位根与随机游走时间序列通常被建模为以下自回归形式的随机过程，因此当前值是过去值的加权和，加上随机扰动：y_t=a_1y_{t-1}+a_2y_{t-2}+\cdots+a_py_{t-p}+\epsilon_t \\我们将在下一节关于单变量时间序列模型的内容中更详细地探讨这些模型，作为ARIMA模型的AR构建模块。这样一个过程有一个如下形式的特征方程：m^p-m^{p-1}a_1-m^{p-2}a_2-\cdots-a_p=0 \\如果这个多项式的（最多）p个根中的一个等于1，那么这个过程被称为有单位根。它将是非平稳的，但不一定有趋势。如果特征方程的其余根的绝对值小于1，那么该过程的1阶差分将是平稳的，**该过程是1阶积分或\mathbf{I}(1)**。如果其余的根的绝对值大于1，那么积分的阶数就会更高，需要进行额外的差分。在实践中，利率或资产价格的时间序列往往不是平稳的，因为没有一个价格水平可以让该序列还原。最突出的非平稳序列的例子是随机游走。给定一个价格的时间序列p_t，其起始价格p_0（例如，一只股票的IPO价格）和白噪声扰动\epsilon_t，那么随机游走统计符合以下自回归关系：p_t=p_{t-1}+\epsilon_t=\sum_{s=0}^t\epsilon_0+p_0 \\重复替换表明，当前的价格p_t是所有先前的扰动\epsilon_t和初始价格p_0的总和。如果该方程包括一个常数项，那么就说随机游走有漂移。因此，随机游走是一个自回归随机过程，其形式如下：y_t=a_1y_{t-1}+\epsilon_t,\;a_1=1 \\它的特征方程m-a_1=0，有单位根，既是非平稳的，又是1阶积分的。一方面，考虑到\epsilon的独立同分布性质，时间序列的方差等于t\sigma^2，这不是2阶平稳的，并意味着原则上序列可以随着时间的推移取任意值。另一方面，取1阶差分\Delta p_t=p_t-p_{t-1}，使得\Delta p_t=\epsilon_t在关于\epsilon的统计假设下是平稳的。具有单位根的非平稳序列的典型特征是长记忆：由于当前值是过去干扰的总和，大的创新比均值回复的平稳序列持续的时间长得多。9.2.2.2
如何诊断单位根统计单位根检验是客观地确定是否需要进行（额外的）差分的一种常见方法。这些是平稳性的统计假设检验，目的是确定是否需要进行差分。增强Dickey—Fuller检验（Augmented Dickey-Fuller Test，ADF test）评估了时间序列样本有单位根的无效假设和平稳性的替代假设。它将被差分的时间序列回归到时间趋势、1阶滞后期和所有滞后期差分上，并从滞后期时间序列值的系数中计算出检验统计量。 statsmodels库使其易于实现（参见笔记本tsa_and_stationarity）。从形式上看，一个时间序列的增强Dickey—Fuller检验运行线性回归，其中\alpha是一个常数，\beta是一个时间趋势的系数，p是指模型中使用的滞后期的数量：\Delta y_t=\alpha+\beta t+\gamma y_{t-1}+\delta_1\Delta y_{t-1}+\cdots+\delta_{p-1}\Delta y_{t-p+1}+\epsilon_t \\\alpha=\beta=0的约束意味着随机行走，而只有\beta=0意味着有漂移的随机行走。滞后期阶数通常使用第7章”线性模型——从风险因子到收益预测“中介绍的Akaike信息准则（AIC）和Bayesian信息准则（BIC）决定。增强Dickey—Fuller检验统计量使用样本系数 ，样本系数在单位根非平稳性的无效假设下，等于零，否则为负。本文的目的是说明，对于一个单整序列，其滞后期序列值不应提供预测滞后期差值以外和滞后差值以外的1阶段差值的有用信息。9.2.2.3
如何去除单位根并处理所产生的序列除了使用相邻数据点之间的差分来消除恒定的变化模式外，我们还可以应用季节性差分来消除季节性变化的模式。这涉及到在代表季节性模式长度的滞后期距离上取值的差分。对于月度数据，这通常涉及滞后期12的差分，而对于季度数据，这涉及时间价格4的差分，以消除季节性和线性趋势。确定正确的转换，特别是确定适当的差分数量和滞后期，并不总是一目了然。有人提出了一些启发式方法，可以归纳为以下几点：1期滞后期自相关接近零或负值，或自相关一般较小且无模式：不需要高阶差分。滞后期达10个以上的正自相关：该序列可能需要高阶差分。1期滞后期自相关<-0.5：该序列可能被过度差分了。可以用AR或MA项来纠正轻微的过度差分或不足（见下一节关于单变量时间序列模型）。一些作者建议将部分差分作为使单整序列平稳的更灵活的方法，并且可能比离散间隔的简单或季节性差分能够保留更多的信息或信号。例如，见Marcos Lopez de Prado（2018）中的第5章”投资组合优化和绩效评估“。9.2.3
实践中的时间序列变换图9.2中的图表显示了纳斯达克股票和工业生产指数30年来的时间序列的原始形式，以及分别应用对数和随后应用1阶和季节性差分（时间价格12）之后的转换版本。图表还显示了增强Dickey—Fuller检验的p值，它使我们能够拒绝两种情况下所有转换后的单位根非平稳性的假设：我们可以使用Q-Q图来进一步分析相关的时间序列特征，该图将时间序列观测值的分布量级与正态分布的量级以及基于自相关函数和部分自相关函数的相关图进行比较。对于图9.3中的纳斯达克图，我们可以看到，虽然没有趋势，但方差不是恒定的，而是在20世纪80年代末、2001年和2008年的市场动荡时期显示出集群式的高峰。Q-Q图突出了分布的肥尾，其极端值比正态分布所显示的更为频繁。自相关函数和部分自相关函数表现出相似的模式，在几个延迟时间的自相关似乎很显著：对于工业制造业生产的月度时间序列，我们可以看到2008年危机后的一个巨大的负离群点，以及Q-Q图中相应的倾斜（见图9.4）。自相关比纳斯达克的收益率要高得多，并且平稳地下降。部分自相关函数在滞后期1和13显示出明显的正自相关模式，在时间价格3和4显示出明显的负系数：9.3
单变量时间序列模型多重线性回归模型将感兴趣的变量表示为输入的线性组合，加上一个随机扰动。相比之下，单变量时间序列模型将时间序列的当前值与该序列的滞后值、当前噪声和可能的过去噪声项的线性组合联系起来。指数平滑模型是基于对数据中趋势和季节性的描述，而ARIMA模型旨在描述数据中的自相关。ARIMA(p, d, q)模型需要平稳性并利用两个构件：自回归（Autoregressive，AR）项由时间序列的p个滞后值组成移动平均（Moving Average，MA）项含由q个滞后扰动组成I代表差分，因为该模型可以通过微分d次来解释单位根非平稳性。术语自回归强调，ARIMA模型意味着时间序列对其自身值的回归。我们将介绍ARIMA构件，AR和MA模型，并解释如何在自回归移动平均（Autoregressive Moving-Average，ARMA）模型中结合它们，这些模型可以像ARIMA模型那样考虑序列整合，或者像AR(I)MAX模型那样包括外生变量。此外，我们将说明如何包括季节性的AR和MA项来扩展工具箱，使其也包括SARMAX模型。9.3.1
如何建立自回归模型p阶AR模型旨在捕捉不同滞后时间序列值之间的线性相关性，可以表示为：AR(p):\;y_t=\phi_0+\phi_1y_{t-1}+\phi_2y_{t-2}+\cdots+\phi_py_{t-p}+\epsilon_t,\;\epsilon\sim i.i.d \\这与对y的滞后值进行多元线性回归非常相似。该模型具有如下特征方程：1-\phi_1x-\phi_2x^2-\cdots-\phi_px^p=0 \\这个p次多项式在x中的解的倒数是特征根，如果所有根的绝对值都小于1，AR(p)过程是平稳的，否则是不稳定的。对于平稳序列，多步预测将收敛于该序列的均值。我们可以用熟悉的最小二乘法估计模型参数，使用p+1,\cdots,T观测值，以确保每个滞后项和结果都有数据。9.3.1.1
如何识别滞后期的数量在实践中，挑战在于决定滞后项的适当阶数p。我们在 "如何测量自相关 "一节中讨论过的序列相关的时间序列分析工具，在做出这一决定时发挥了关键作用。更具体地说，对相关图的目视检查经常提供有用的线索：自相关函数估计了不同滞后期观测值之间的自相关性，这反过来又源于直接和间接的线性相关性。因此，如果一个k阶的AR模型是正确的模型，自相关函数将显示一个显著的序列相关性，直到滞后k，并且由于线性关系的间接影响造成的惯性，将延伸到随后的滞后期，直到它最终随着影响的减弱而减弱。反过来，局部自相关函数只测量相隔一定滞后期的观测值之间的直接线性关系，所以它不会反映出超过k滞后期的相关性。9.3.1.2
如何诊断模型拟合如果模型正确地捕捉了跨滞后期的线性相关性，那么残差应该类似于白噪声，而且自相关函数应该强调没有明显的自相关系数。除了残差图之外，Ljung-Box Q-统计量允许我们检验残差序列遵循白噪声的假设。无效假设是所有的m个序列相关系数都是零，而另一个假设是一些系数不是零。检验统计量是由不同滞后期k的样本自相关系数k计算出来的，并遵循X^2分布：Q(m)=T(T+2)\sum_{t=1}^m\frac{\rho_L^2}{T-l} \\正如我们将看到的，statsmodels库提供了关于不同滞后期的系数的显著性的信息，不显著的系数应该被删除。如果Q统计量拒绝了没有自相关的无效假设，你应该考虑额外的AR项。9.3.2
如何建立移动平均模型MA(q)模型在类似回归的模型中使用q个过去的扰动而不是时间序列的滞后值，如下所示：MA(q):\;y_t=c+\epsilon_t+\theta_1\epsilon_{t-1}+\theta_2\epsilon_{t-2}+\cdots+\theta_p\epsilon_{t-p}\;\epsilon\sim i.i.d \\由于我们没有观察到白噪声扰动值\epsilon_t，所以MA(q)不是像我们目前看到的那样的回归模型。MA(q)模型不是使用最小二乘法，而是使用最大似然法（Maximum Likelihood，MLE）进行估计，在序列开始时初始化或估计干扰值，然后递归和迭代地计算剩余部分。MA(q)模型的名字来自于将y_t的每个值表示为过去q个创新点的加权移动平均。换句话说，当前的估计值代表了相对于模型过去错误的修正。在MA(q)模型中使用移动平均不同于指数平滑或季节性时间序列成分的估计，因为MA(q)模型的目的是预测未来的数值，而不是去噪或估计过去数值的趋势周期。MA(q)过程总是平稳的，因为它们是白噪声变量的加权和，这些变量本身是平稳的。9.3.2.1
如何识别滞后期的数量一个由MA(q)过程产生的时间序列是由之前q模型预测的残差驱动的。因此，MA(q)过程的自相关函数将显示直到滞后q的值的显著系数，然后急剧下降，因为这是模型假设序列值产生的方式。请注意这与我们刚才描述的AR情况有什么不同，在AR情况下，部分自相关函数会显示类似的模式。9.3.2.2
AR模型和MA模型之间的关系一个AR(p)模型总是可以用重复替换的过程来表达。如在“如何处理单位根引起的随机趋势”一节中的随机游走例子。当MA(q)过程的系数满足一定的规模约束时，它也会变得可逆，可以表示为一个AR(\infty)过程（详见Tsay，2005）。9.3.3
如何构建ARIMA模型和扩展自回归移动平均数—ARIMA(p, d, q)—模型结合了AR(p)和MA(q)过程，利用这些构件的互补性，简化了模型的开发。它们采用更紧凑的形式，减少了参数的数量，从而降低了过拟合的风险。该模型还利用时间序列的第d次差分消除了单位根非平稳性。使用序列1阶差分的ARIMA(p,1,q)模型与ARMA(p,q)模型相同。用y'表示经过d次非季节性差分后的原始序列，ARIMA(p, d, q)模型很简单：\begin{align} ARIMA(p,d,q):\quad y'_t&=AR(P)+MA(q) \\ &=\phi_0+\phi_1y'_{t-1}+\cdots+\phi_py'_{t-p}+\epsilon_t+ \theta_1\epsilon_{t-1}+\cdots+\theta_q\epsilon_{t-p},\;\epsilon\sim i.i.d. \end{align} \\ARIMA模型也可以用极大似然估计法进行估计。根据实现的不同，高阶模型通常包含低阶模型。例如，到0.11版本为止，statsmodels库包括所有低阶p和q项，并且不允许删除低于最高值的滞后期的系数。在这种情况下，高阶模型总是能更好地拟合。要注意不要因为使用了太多的项而使你的模型与数据过度拟合。最近的版本，也就是本文写作时的0.11版本，增加了一个实验性的新ARIMA模型，具有更灵活的配置选项。9.3.3.1
如何对差分级数建模在使用数据时，也有一些设计单变量时间序列模型的指导准则：一个没有差分的模型假定原始序列是平稳的，包括均值回归的。它通常包括一个常数项，以允许非零均值的存在。一个有一阶差分的模型假设原始序列有一个恒定的趋势，因此应该包括一个常数项。二阶差分的模型假定原始序列有一个时间变化的趋势，且不应包括常数。9.3.3.2
如何识别AR和MA项的数量由于AR(p)和MA(q)项相互影响，自相关函数和部分自相关函数提供的信息不再可靠，只能作为一个起点。传统上，在选择模型设计时，Akaike信息准则和Bayesian信息准则被用来依赖样本内拟合。另外，我们可以依靠样本外测试来交叉验证多个参数选择。下面的总结为在单独考虑AR和MA模型的情况下如何选择模型阶数提供了一些指导：超过部分自相关函数中断的滞后期就是AR项的数量。如果差分序列的部分自相关函数显示出急剧下降和/或滞后1自相关系数为正，则增加一个或多个AR项。超过自相关函数中断的滞后期就是MA项的数量。如果差分序列的自相关函数显示出急剧的下降和/或滞后1自相关系数为负，考虑在模型中增加一个MA项。AR项和MA项可能会相互抵消，因此如果模型同时包含AR项和MA项，请尽量将它们的数量减少1，以避免过拟合，特别是当更复杂的模型需要超过10次迭代才能收敛时。如果AR系数之和接近1，表明模型的AR部分存在单位根，那么就消除一个AR项，并对模型进行一（更多）次差分。如果MA系数之和接近1，表明模型的MA部分存在单位根，则消除一个MA项，并将差分的阶数减少1。不稳定的长期预测表明在模型的AR或MA部分可能有一个单位根。9.3.3.3
增加特征—ARMAX带有外生输入的自回归移动平均模型（Autoregressive Moving-Average with Exogenous Inputs，ARMAX）模型在ARMA时间序列模型的右侧增加了输入变量或协变量（假设序列是平稳的，所以我们可以跳过差分）：\begin{align} ARMAX(p,d,q):\quad y_t&=\beta x_t+AR(P)+MA(q) \\ &=\beta x_t+\phi_0+\phi_1y_{t-1}+\cdots+\phi_py_{t-p}+\epsilon_t+ \theta_1\epsilon_{t-1}+\cdots+\theta_q\epsilon_{t-p},\;\epsilon\sim i.i.d. \end{align} \\这类似于线性回归模型，但解释起来相当困难。这是因为\beta对y_t的影响并不像线性回归中那样是x增加一个单位的影响。相反，方程右侧存在y_t的滞后值，这意味着系数只能根据响应变量的滞后值来解释，这是不直观的。9.3.3.4
增加季节差异—SARIMAX对于具有季节性效应的时间序列，我们可以引入AR和MA项来捕捉季节性的周期性。例如，当使用月度数据且季节性效应长度为1年时，季节性AR和MA项将反映这个特定的滞后长度。ARIMAX(p, d, q)模型就变成了SARIMAX(p, d, q)\times(P, D, Q)模型，写出来比较复杂，但statsmodels库文档（见GitHub上的链接）详细提供了这些信息。现在我们将使用宏数据建立一个季节性的ARMA模型来说明其实现。9.3.4
如何预测宏观基本面我们将为1988-2017年期间的工业生产时间序列的月度数据建立一个SARIMAX模型。正如关于分析工具的第一节中所说明的，数据 已经进行了对数转换，我们使用的是季节性（滞后12）差分。我们使用了10年的训练数据的滚动窗口对普通和传统的AR和MA参数范围的模型进行了估计，并评估超前一步预测的均方根误差（Root Mean Square Error，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)]我们还集成了Akaike信息准则和Bayesian信息准则，显示出非常高的秩相关系数0.94，Bayesian信息准则比Akaike信息准则更倾向于参数略少的模型。按均方根误差计算，最好的五个模型是：
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) \times (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())我们得到以下摘要：系数是显著的，而且Q统计量拒绝了进一步自相关的假设。相关图同样表明，我们已经成功地消除了该系列的自相关：9.3.5
如何使用时间序列模型预测波动率单变量时间序列模型在金融领域的一个特别重要的应用是对波动性的预测。这是因为随着时间的推移，波动率通常不是恒定的，而是聚集在一起的。方差的变化给使用假设平稳性的经典ARIMA模型进行时间序列预测带来了挑战。为了应对这一挑战，我们现在将对波动率进行建模，以便预测方差的变化。异方差是描述一个变量方差变化的技术术语。ARCH模型将误差项的方差表示为以前各期误差的函数。更具体地说，它假设误差方差遵循AR(p)模型。广义自回归条件异方差（Generalized Autoregressive Conditional Heteroskedasticity，GARCH）模型扩大了ARCH模型的范围，允许ARMA模型。时间序列的预测通常结合ARIMA模型的预期平均值和ARCH/GARCH模型的预期方差。2003年的诺贝尔经济学奖授予了Robert Engle和Clive Granger，以表彰他们开发了这一类模型。前者还管理着纽约大学斯特恩学院的波动率实验室（http://vlab.sten.nyu.edu），该实验室有许多关于我们将要讨论的模型的在线例子和工具。9.3.5.1
ARCH模型ARCH(p)模型只是一个应用于时间序列模型残差方差的AR(p)模型，它使t时刻的方差以方差的滞后观测值为条件。更具体地说，误差项\epsilon_t是线性模型（如ARIMA）在原始时间序列上的残差，分为随时间变化的标准差\sigma_t和扰动z_t，具体如下：\begin{align} ARCH(p):\quad var(x_t)&=\sigma_t^2 \\ &=\omega+\alpha_1\epsilon_{t-1}^2+\cdots+\alpha_p\epsilon_{t-p}^2 \\ \epsilon_t&=\sigma_tz_t \\ z_t&\sim i.i.d \end{align} \\ARCH(p)模型可以用普通最小二乘法来估计。Engle提出了一种使用拉格朗日乘数检验来识别适当的ARCH阶数的方法，它相当于对线性回归中所有系数为零的假设的F检验（见第7章“线性模型——从风险因子到收益预测”）。ARCH模型的一个主要优点是它产生的波动率估计值具有正的超额峰度，即相对于正态分布的肥尾，这反过来又与关于收益的经验观察相一致。缺点包括假设正负波动率冲击的影响相同，而资产价格的反应往往不同。它也不能解释波动率的变化，而且很可能过度预测波动率，因为它们对收益率序列的巨大、孤立的冲击反应缓慢。对于一个适当指定的ARCH模型，标准化的残差（除以模型对标准差周期的估计）应该类似于白噪声，可以进行Ljung-Box Q检验。9.3.5.2
推广ARCH——GARCH模型ARCH模型相对简单，但往往需要许多参数来捕捉资产收益序列的波动模式。GARCH模型适用于具有扰动\epsilon_t=r_t-\mu的对数收益序列r_t，该序列遵循GARCH(p,q)模型，如：\epsilon_t=\sigma_tz_t,\; \sigma_t^2=\omega+\sum_{i=1}^p\alpha_i\epsilon_{t-i}^2+\sum_{j=1}^q\beta_i\sigma_{t-j}^2\; z_t\sim i.i.d \\对于误差项的方差\epsilon_t，GARCH(p, q)模型假设为ARMA(p, q)模型。与ARCH模型类似，GARCH(1,1)过程的尾部分布比正态分布更重。该模型遇到了与ARCH模型相同的弱点。例如，它对正负冲击的反应相同。为了配置ARCH和GARCH模型的滞后阶数，使用训练好的时间序列的平方残差来预测原始序列的平均值。残差是以零为中心的，所以它们的平方也是方差。然后，检查平方残差的自相关函数和部分自相关函数图，确定时间序列方差的自相关模式。9.3.5.3
如何建立一个预测波动率的模型为资产收益序列建立波动率模型包括四个步骤：根据自相关函数和部分自相关函数所揭示的序列依赖性，为金融时间序列建立一个ARMA时间序列模型测试模型的残差是否具有ARCH/GARCH效应，同样依靠自相关函数和部分自相关函数对残差平方序列的测试如果序列相关效应显著，则指定一个波动率模型，并联合估计均值和波动率方程仔细检查拟合的模型，必要时对其进行改进在对收益率序列进行波动率预测时，可以对序列依赖性进行限制，这样就可以使用常数平均值来代替ARMA模型。arch库（见GitHub上的文档链接）提供了几个方法来估计波动率预测模型。你可以将预期均值建模为一个常数，也可以建模为AR(p)模型，如“如何建立自回归模型”一节中所讨论的，还可以建模为更近期的异质自回归过程（Heterogeneous Autoregressive Processes，HAR），它使用每日（1天）、每周（5天）和每月（22天）的滞后期来捕捉短期、中期和长期投资者的交易频率。均值模型可以用几个条件异方差模型共同定义和估计，这些模型除了ARCH和GARCH之外，还包括指数GARCH（Exponential GARCH，EGARCH）模型，它允许正负收益之间的不对称效应，以及异质ARCH（Heterogeneous，HARCH）模型，它是对HAR均值模型的补充。我们将使用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重新缩放的每日收益序列只表现出有限的自相关，但与平均值的平方偏差确实有大量的记忆，反映在缓慢衰减的自相关函数和部分自相关函数上，前两个延迟时间很高，只有在前六个滞后期后才截止： plot_correlogram(nasdaq_returns.sub(nasdaq_returns.mean()).pow(2), lags=120,
title='NASDAQ Daily Volatility')函数plot_correlogram的输出如下：因此，我们可以估计一个GARCH模型来捕捉过去波动率的线性关系。我们将使用10年的滚动窗口来估计GARCH(p, q)模型，p和q的范围是1-4，以产生一步样本外预测结果。然后我们比较预测波动率的均方根误差与实际收益率与其平均值的平方差，以确定最具预测性的模型。 我们使用winsorized数据来限制极端收益值反映在波动率的非常高的正偏度上的影响：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)模型实现了最低的均方根误差（与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())输出显示了最大化的对数似然，以及Akaike信息准则和Bayesian信息准则，这些标准在根据样本内表现选择模型时通常是最小化的（见第7章“线性模型——从风险因子到收益预测”）。它还显示了均值模型的结果，在本例子中，它只是一个常数估计，以及GARCH参数的常数omega，AR参数\alpha和MA参数\beta，所有这些都是统计显著的：现在让我们探索多时间序列模型和协整的概念，这将实现新的交易策略。9.4
多元时间序列模型多变量时间序列模型旨在同时捕捉多个时间序列的动态，并利用这些序列之间的依赖关系进行更可靠的预测。对这个问题最全面的介绍是Lütkepohl（2005）。9.4.1
等效系统方法单变量时间序列模型，如我们刚才讨论的ARMA方法，仅限于目标变量与其滞后值或滞后干扰和外生序列之间的统计关系，在ARMAX的情况下。相反，多变量时间序列模型还允许其他时间序列的滞后值影响目标变量。这种影响适用于所有序列，导致复杂的相互作用，如下图所示：除了潜在的更好的预测外，多变量时间序列还被用来深入了解跨序列的依赖性。例如，在经济学中，多变量时间序列被用来了解一个变量（如利率）的政策变化如何在不同时间范围内影响其他变量。多变量模型产生的脉冲响应函数就可以达到这个目的，它允许我们模拟一个变量如何对其他变量的突然变化做出反应。Granger因果关系的概念分析了一个变量是否有助于预测另一个变量（在最小二乘意义上）。此外，多变量时间序列模型允许对预测误差方差进行分解，以分析其他序列的贡献。9.4.2
向量自回归（VAR）模型我们将看到向量自回归VAR(p)模型如何通过建立一个由k个方程组成的系统，将AR(p)模型扩展到k个序列，其中每个方程包含所有k个序列的p个滞后值。在最简单的情况下，k=2的VAR(1)模型有如下形式：y_{1,t}=c_1+\alpha_{1,1}y_{1,t-1}+\alpha_{1,2}y_{2,t-2}+\epsilon_{1,t} \\ y_{2,t}=c_2+\alpha_{2,1}y_{2,t-1}+\alpha_{2,2}y_{2,t-2}+\epsilon_{2,t} \\这个模型可以用矩阵形式更简洁地表示：\left[\begin{aligned} y_{1,t} \\ y_{2,t} \end{aligned} \right] =
\left[\begin{aligned} c_1 \\ c_2 \end{aligned} \right] + \begin{bmatrix} a_{1,1} & a_{1,2} \\ a_{2,1} & a_{2,2} \end{bmatrix} \left[\begin{aligned} y_{1,t-1} \\ y_{2,t-2} \end{aligned} \right] + \left[\begin{aligned} \epsilon_{1,t} \\ \epsilon_{2,t} \end{aligned} \right] \\输出的滞后值的系数提供了关于序列本身的动态信息，而交叉变量系数则提供了关于序列之间相互作用的信息。这种表示方法可以扩展到k个p阶时间序列，如下所示：\begin{array}\pmb{y}_t \\ k\times 1\end{array}=\begin{array}\pmb{c} \\ k\times 1\end{array}+ \begin{array}\pmb{A}_t \\ k\times k\end{array}\begin{array}\pmb{y}_{t-1} \\ k\times 1\end{array}+ \cdots+\begin{array}\pmb{A}_p \\ k\times k\end{array}\begin{array}\pmb{y}_{t-p} \\ k\times 1\end{array} +\begin{array}\pmb{\epsilon}_t \\ k\times 1\end{array} \\VAR(p)模型也需要平稳性，以便单变量时间序列建模的初始步骤可以延续。首先，探索序列并确定必要的转换。然后，应用增强Dickey—Fuller检验来验证每个序列是否满足平稳性标准，否则应用进一步的转换。可以用以初始信息为条件的普通最小二乘法来估计，也可以用极大似然估计法来估计，对于正态分布的误差来说，与极大似然估计法是等价的，但在其他情况下则不然。如果k个序列中的一些或全部是单位根非平稳的，它们可能是协整的（见下一节）。这种将单位根概念扩展到多个时间序列的做法意味着两个或多个序列的线性组合是平稳的，因此可以实现均值回归。VAR模型无法在不进行差分的情况下处理这种情况，相应地使用向量误差修正模型（Vector Error Correction Model，VECM，Johansen和Juselius 1990）。我们将进一步探讨协整问题，因为如果存在协整并假定其持续存在，它可以被用于配对交易策略。滞后阶数的确定也是以每个序列的自相关函数和部分自相关函数为线索，但由于所有序列都具有相同的滞后阶数，因此滞后阶数的确定受到限制。在模型估计后，残差诊断也需要类似于白噪声的结果，模型选择可以使用样本内信息标准，如果目标是使用模型进行预测，则使用样本外预测性能来交叉验证其他模型设计。正如在单变量情况下所提到的，对原始时间序列的预测需要我们在训练模型之前将应用于使一个序列平稳的变换进行逆转。9.4.3
运用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这给我们留下了以下序列：为了限制输出的规模，我们将只是使用statsmodels库的VARMAX实现（允许可选的外生变量），用前480个观察值估计一个具有恒定趋势的VAR(1)模型： model = VARMAX(df_transformed.loc[:'2017'], order=(1,1), trend='c').fit(maxiter=1000)这会产生以下摘要：输出包含两个时间序列方程的系数，如前面的VAR(1)插图中所概述的。 statsmodels库提供了诊断图来检查残差是否符合白噪声假设。在这个简单的例子中，情况并非如此，因为方差似乎并不恒定（左上），而且分位数图显示了分布的差异，即肥尾（左下）：你可以像下面这样生成样本外预测：preds = model.predict(start=480, end=len(df_transformed)-1)下面可视化了实际值和预测值，显示了预测如何滞后于实际值，并且不能很好地捕捉非线性、样本外模式：9.5
协整—具有共同趋势的时间序列在上一节多元时间序列模型中，我们简要地提到了协整。现在让我们解释这个概念，以及如何在利用它进行统计套利交易策略之前更详细地诊断它的存在。我们已经看到了一个时间序列是如何具有单位根的，它创造了一个随机的趋势，并使时间序列具有高度的持久性。当我们用这样的集成时间序列的原始形式，而不是差分形式，作为线性回归模型的一个特征时，它与结果的关系往往会在统计上显得显著，尽管事实并非如此。这种现象被称为伪回归（详见第18章”金融时间序列和卫星图像的CNNs", Wooldridge, 2008）。因此，推荐的解决方案是对时间序列进行差分，使其在模型中使用之前成为平稳的。然而，当结果与一个或多个输入变量之间存在协整关系时，就有一个例外。为了理解协整的概念，首先让我们记住，回归模型的残差是输入和输出序列的线性组合。通常情况下，一个集成时间序列在一个或多个这样的序列上的回归残差也会产生非平稳的集成残差，因此表现为随机游走。然而，对于某些时间序列，情况并非如此：回归系数的线性组合的残差是平稳的，即使单个序列不是平稳的。这样的时间序列是协整的。一个非技术性的例子是，一个醉汉在他的狗（拴着的）陪同下随机散步。两条轨迹都是非平稳的，但却是协整的，因为狗偶尔会回到它的主人身边。在交易方面，套利约束意味着现货和期货价格之间的协整。换句话说，两个或多个协整序列的线性组合有一个稳定的平均值，这个线性组合会还原到这个平均值。这也适用于单个序列的高阶积分，线性组合降低了整体的积分阶。协整与相关不同：两个序列可以高度相关，但不一定是协整的。例如，如果两个增长的序列是彼此的恒定倍数，它们的相关性会很高，但任何线性组合也会增长，而不是恢复到一个稳定的平均值。协整是非常有用的：如果两个或更多的资产价格序列倾向于恢复到一个共同的平均值，我们可以利用偏离趋势的情况，因为它们应该意味着未来的价格移动方向相反。协整背后涉及的数学问题较多，所以我们将只关注实践方面；关于深入的处理，见Lütkepohl（2005）。在本节中，我们将讨论如何识别具有这种长期平稳关系的配对，估计任何不平衡修正的预期时间，以及如何利用这些工具来实施和回溯测试一个多空货币对交易策略。有两种方法可以检验协整关系：Engle—Granger两步法Johansen检验法我们将依次讨论每一种方法，然后展示它们如何帮助识别倾向于恢复到共同趋势的协整证券，我们可以利用这一事实来进行统计套利策略。9.5.1
Engle—Granger两步法Engle—Granger方法用于识别两个序列之间的协整关系。它涉及以下两个方面：用一个序列对另一个序列进行回归来估计平稳的长期关系对回归残差进行自相关函数单位根检验无效假设是残差有单位根，并且是集成的；如果我们能拒绝它，那么我们就假设残差是平稳的，因此，序列是协整的（Engle和Granger 1987）。这种方法的一个主要好处是，回归系数代表了使组合平稳的乘数，也就是说，均值可回归。不幸的是，检验结果会有所不同，这取决于我们认为哪个变量是独立的，因此，我们尝试两种方法，然后挑选p值较低的、检验统计量较较低的关系。另一个缺点是，这种检验仅限于成对的关系。更复杂的Johansen程序可以识别多达十几个时间序列之间的显著协整关系。9.5.2
Johansen似然比率检验而Johansen程序则是检验协整对VAR模型的限制，正如上一节所讨论的。更具体地说，在从一般的VAR(p)模型的两边减去目标向量后，我们得到误差修正模型（Error Correction Model，ECM）的表述：\Delta\pmb{y}_t=\pmb{c}+\pmb{\Pi}\pmb{y}_{t-1}+\pmb{\Gamma}_1\Delta y_{t-1}+\cdots +\pmb{\Gamma}_p\pmb{y}_{t-p}+\pmb{\epsilon}_t \\修改后的VAR(p)方程中只有一个行向量项（\pmb{y}_{t-1}）没有使用\Delta运算符表示差值。协整的性质取决于这一项系数矩阵的秩（Johansen 1991）。虽然这个方程在结构上看起来与自相关函数检验设置相似，但由于涉及多个序列，因此现在有几个共同趋势的潜在。为了确定协整关系的数量，Johansen检验从0（没有协整）开始，依次检验\Pi秩的增加。在下一节中，我们将探讨两个序列案例的应用。Gonzalo和Lee（1998）讨论了由于错误的模型动力学和其他实施方面的实际挑战，包括如何结合我们在下一节的样本统计套利策略中所依靠的两种检验程序。9.6
带有协整的统计套利统计套利是指采用某种统计模型或方法，利用看似相对错误的资产定价，同时保持一定程度的市场中性的策略。配对交易是一种概念简单的策略，至少从80年代中期开始就被算法交易者所采用（Gatev, Goetzmann，和Rouwenhorst 2006）。其目的是找到历史上价格一起变动的两种资产，跟踪价差（它们之间的价格差异），一旦价差扩大，买入跌破共同趋势的失败者，并做空胜利者。如果这种关系持续下去，随着价格的收敛和头寸的平仓，多头和/或空头将带来利润。这种方法扩展到多元环境，由多种证券组成篮子，并将一种资产与一个篮子或两个篮子彼此进行交易。在实践中，该策略需要两个步骤：形成阶段：识别具有长期均值回归关系的证券。理想情况下，价差应具有较高的方差，以便在可靠地回复到共同趋势的同时频繁地进行有利可图的交易。交易阶段：当价格变动导致价差背离和收敛时，触发进入和退出交易规则。在过去的几年里，在这个领域日益活跃的研究中，出现了几种形成和交易阶段的方法，涉及多个资产类别。下一小节概述了关键的差异，然后我们深入到一个应用实例中。9.6.1
如何选择和交易流动资产配对最近对成对交易策略的全面调查（Krauss 2017）确定了四种不同的方法，加上其他一些较新的方法，包括基于机器学习的预测：距离方法：这是最古老和研究最多的方法，用距离指标如相关性来确定候选配对，并使用非参数阈值如布林线来触发进入和退出交易。自Gatev等人（2006年）以来，其计算的简单性允许大规模的应用，并在整个市场和资产类别中表现出盈利能力。然而，最近的表现却有所下降。协整方法：如前所述，这种方法依赖于两个或多个变量之间长期关系的计量经济学模型。并允许进行统计检验，承诺比简单的距离度量更可靠。这一类的例子使用Engle-Granger和Johansen程序来识别证券配对和证券篮子，以及旨在捕捉概念的更简单的启发式方法（Vidyamurthy 2004）。交易规则通常类似于使用距离度量的简单阈值。时间序列方法：以交易阶段为重点，这类策略旨在将价差建模为均值回归的随机过程，并相应优化进入和退出规则（Elliott, Hoek，和Malcolm 2005）。它假定已经确定了有希望的配对。随机控制方法：与时间序列方法类似，目标是利用随机控制理论优化交易规则，找到价值和策略函数，以达到最佳的投资组合（Liu和Timmermann 2013）。我们将在第21章“合成时间序列数据的生成对抗网络”中讨论这种类型的方法。其他方法：除了基于主成分分析等无监督学习（见第13章“数据驱动的风险因子和无监督学习的资产配置”）和协同学等统计模型（Patton 2012）的配对识别，机器学习最近也开始流行，根据其相对价格或收益预测来识别配对（Huck 2019）。我们将在接下来的章节中介绍几种可用于此目的的机器学习算法，并说明相应的多变量交易策略。以上对各种方法的总结，仅仅是对配对交易策略设计所提供的灵活性的一瞥。除了关于配对选择和交易规则逻辑的更高层次的问题外，还有许多参数需要我们定义来实现。这些参数包括以下内容：筛选潜在配对或篮子的投资范围形成期的长度用于挑选可交易候选者的关系强度当价差波动时，对触发进入或退出交易或调整现有头寸的常用手段的偏离程度和趋同程度9.6.2
配对交易实践距离方法使用（归一化）资产价格或其收益的相关性来识别配对，它很简单，而且比协整测试的计算量要少很多。笔记本cointegration_test说明了这一点，它是以4年的每日数据为样本的150只股票为例：计算与ETF收益的相关性需要30ms，而一套协整检验（使用statsmodels库）需要18秒，慢了600倍。速度优势特别有价值。这是因为潜在配对的数量是两边要考虑的候选者数量的乘积，所以评估100只股票和100只ETF的组合需要比较10,000次测试（我们将在后面讨论多重测试偏差的挑战）。另一方面，距离指标不一定会选择最有利可图的配对：完美的联合移动会使相关性最大化，这反过来又会消除实际的交易机会。实证研究证实，协整配对的价差波动率几乎是距离配对价差波动率的两倍（Huck和Afawubo 2015）。为了平衡计算成本和所产生的配对质量之间的权衡，Krauss（2017）根据他的文献综述，推荐了一个结合两种方法的程序：选择分布稳定且漂移小的对，以减少候选集的数量检验剩下的具有最高价差的配对是否存在协整关系这个过程的目的是选择具有较低分化风险的协整对，同时确保更多的波动价差，反过来，产生更多的盈利机会。大量的测试会引入数据窥探的偏差，正如第6章“机器学习过程中”所讨论的：多次测试很可能会增加误判拒绝无协整性的无效假设的假阳性数量。虽然统计意义可能不是交易盈利的必要条件（Chan，2008），但对商品对的研究（Cummins和Bucca，2012）表明，根据Romano和Wolf（2010），控制家族误差率以提高测试的力量，可以带来更好的性能。在下一小节中，我们将更仔细地研究各种针对资产价格整合程度的启发式方法对协整检验结果的预测能力。该示例代码使用了在纽约证券交易所和纳斯达克交易所的172只股票和138只ETF的样本，并由Stooq提供了2010-2019年的每日数据。证券代表了其各自类别在样本期内最大的平均美元交易量；高度相关和固定的资产已被删除。关于如何获得数据的说明，请参见GitHub存储库数据文件夹中的笔记本create_datasets，关于相关代码和额外的预处理和探索细节的笔记本cointegration_tests。9.6.2.1
基于距离的启发式寻找协整配对compute_pair_metrics()计算了2010-14年和2015-19年超过23,000对股票和交易所交易基金（ETF）的以下距离指标：价差的漂移，定义为时间趋势对价差的线性回归价差的波动性归一化的价格序列之间以及它们的收益之间的相关性低漂移和波动性以及高相关性是协整的简单代理。为了评估这些启发式方法的预测能力，我们还使用statsmodels库对前面的配对进行Engle-Granger和Johansen协整检验。这发生在compute_pair_metrics()的后半部分的循环中。我们首先估计我们需要为Johansen检验指定的最佳滞后期数量。对于这两个检验，我们假设协整序列（价差）可能有一个不同于零的截距，但没有趋势。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的Johansen跟踪统计量与它们各自的临界值进行比较，得到Engle-Granger的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多配对关系，Johansen检验认为3.2%的关系是显著的，而Engle-Granger检验则认为6.5%。他们在366对关系（0.79%）上意见一致。9.6.2.2
启发式方法预测显著协整的效果如何？当我们把根据两种检验方法都是协整的序列的启发式方法的分布与不协整的其余序列进行比较时，波动率和漂移确实较低（以绝对值计算）。图9.14显示，两种相关度量的情况不太清楚：为了评估启发式方法的预测准确性，我们首先运行一个具有这些特征的逻辑回归模型来预测显著的协整关系。它取得了一个 曲线下面积（Area-Under-the-Curve，AUC）交叉验证分数为0.815；如果不考虑相关指标，它的分数仍为0.804。决策树的表现稍好，AUC=0.821，无论是否有相关性特征。特别是由于严重的类别不平衡，出现了大量的假阳性：正确识别366个协整对中的80%意味着超过16500个假阳性，但排除了近30000个候选者。更多细节见笔记本cointegration_ tests。关键的启示是，距离启发法可以帮助更有效地筛选一个大范围的数据，但这是以错过一些协整对为代价的，仍然需要大量的测试。9.6.3
准备策略回测在本节中，我们将对2017-2019年期间的股票和ETF的样本实施基于协整的统计套利策略。为了精简表述，有些方面被简化了。有关代码示例和其他细节，请参见笔记本 statistical_arbitrage_with_ cointegrated_pairs。我们首先生成并存储所有候选配对的协整检验以及由此产生的交易信号。然后，考虑到该过程的计算强度，我们对基于这些信号的策略进行回测。9.6.3.1
预计算协整检验首先，我们对23,000个潜在配对中的每一个进行2年回溯期的季度协整测试，然后，我们选择Johansen和Engle- Granger检验都一致的配对进行交易。我们应该排除那些在回溯期内静止的资产，但是我们排除了在整个时期内静止的资产，所以我们跳过这一步来简化它。这一程序遵循之前概述的步骤；详情请见笔记本。图9.15显示了为交易而选择的两个不同配对的原始股票和ETF序列；注意在样本期间明显存在着共同的趋势。9.6.3.2
进行入市和出市交易现在，我们可以根据滚动对冲比率来计算每个候选配对的价差。我们还计算了布林带，因为我们将把价差远离其移动平均线大于两个滚动标准差的移动视为多头和空头进入信号，而移动平均线的反向交叉则是退出信号。用卡尔曼率滤波平滑价格为此，我们首先应用滚动卡尔曼滤波器（Kalman Filter，KF）来去除一些噪音，如第4章“金融特征工程——如何研究阿尔法因素子”所展示的：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)使用卡尔曼滤波器计算滚动套期保值率为了获得动态套期保值比率，我们将卡尔曼滤波器应用于滚动线性回归，如下所示：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估算均值回归的半衰期如果我们将利差视为连续时间内的均值回归随机过程，我们可以将其建模为Ornstein-Uhlenbeck过程。这种观点的好处是，我们可以获得均值回归的半衰期公式，作为偏差后价差再次收敛所需时间的近似值（详见Chan 2013，第2章“市场和基本面数据——来源和技术”）：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分数来推导出交易信号：如果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_trades9.6.4
使用backtrader库对策略进行回测现在，我们已经准备好在我们的回测平台上制定我们的策略，执行它，并评估结果。为此，我们需要跟踪我们的配对，除了个人投资组合头寸外，监测活跃和不活跃配对的分布，以应用我们的交易规则。9.6.4.1
用一个自定义的数据类跟踪配对为了说明活跃的配对，我们定义了一个数据类（在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))9.6.4.2
运行和评估该战略关键的执行方面包括：每天从已触发退出规则或超过特定负收益的配对中退出为那些价差触发了进入信号的配对对开立新的多头和空头头寸此外，我们调整头寸以考虑到不同数量的配对策略本身的代码占用了太多空间，无法在此显示；请看笔记本上的pair_trading_backtest了解详情。图9-16表明，至少在2017-2019年期间，这种简化的策略曾有过辉煌时刻（请注意，我们利用了一些前瞻性偏差，忽略了交易成本）。在这些宽松的假设下，它在期初和期末的表现低于标准普尔500指数，其他方面则基本保持一致（左图）。它的阿尔法值为0.08，贝塔值为-0.14（右图），平均夏普比率为0.75，Sortino比率为1.05（中央图）：虽然我们应该对这些业绩指标持谨慎态度，但该策略展示了基于协整的统计套利的剖析，其形式是配对交易。让我们来看看你可以采取的几个步骤，在这个框架的基础上产生更好的业绩。9.6.5
扩展—如何做得更好协整是一个非常有用的概念，可以识别倾向于一致行动的股票对或股票组。与协整的统计复杂性相比，我们使用的是非常简单和静态的交易规则；按季度计算也扭曲了策略，正如多头和空头持有的模式所示（见笔记本）。要想获得成功，你至少需要筛选一个更大的集合，并优化几个参数，包括交易规则。此外，风险管理应考虑到当某些资产相对频繁地出现在一个交易对的同一侧时产生的集中头寸。你也可以用篮子来操作，而不是单独的配对；然而，为了解决越来越多的候选数量，你可能需要限制篮子的组成。正如在“配对交易—利用协整进行统计套利”一节中提到的，还有一些旨在预测价格变动的替代方法。在接下来的章节中，我们将探讨各种机器学习模型，旨在预测特定投资领域的绝对规模或方向。使用这些预测作为多头和空头入场信号，是我们在本节中研究的配对交易框架的自然延伸或替代方案。9.7
总结在本章中，我们探讨了单个序列的单变量情况下的线性时间序列模型，以及几个相互作用的序列的多变量模型。我们遇到了预测宏观基本面的应用，预测资产或投资组合波动性的模型，以及在风险管理中广泛使用的多变量VAR模型，以及捕捉多个宏观序列的动态变化。我们还研究了协整的概念，它是流行的配对交易策略的基础。与第7章"线性模型——从风险因子到收益预测“类似，我们看到了线性模型如何强加了大量结构，也就是说，它们做出了强有力的假设，有可能需要进行转换和大量测试来验证这些假设是否得到满足。如果符合，模型的训练和解释就很简单，而且这些模型可以提供一个良好的基线，更复杂的模型可以在此基础上进行改进。在接下来的两章中，我们将看到这方面的两个例子，即随机森林和梯度提升模型，我们将在第4部分遇到更多的例子，这部分是关于深度学习的。}
背景描述数据分析中离不开对数据的相关性分析，并且需要把这些相关性进行可视化（绘图），以方便人们对各种特征属性之间呈现出来的相关性有更直接、清晰的感知和理解，提升数据的价值和数据挖掘的效益。本文以“鸢尾花数据集”为基础，主要关注于各种关系图的绘制，以及统计分析的数据可视化，提供和展示了16种关系图及6种统计分析图和回归图的方法（详见以下目录）。由于从sklearn中获取的“鸢尾花”数据集中，目标值（iris.target)是“0”和“1”，这种类型的数据方便实现“机器学习”的建模，但是在数据绘图中不利于理解，因此我们会把数据集中的目标值与鸢尾花的品种(species)进行关联，转换为新的数据集，以获取更好的可视化结果。目录第一部分：鸢尾花数据集的获取、转换并保存为csv文件第二部分：关系图
1. 箱型图
2. 小提琴图 -- violinplot(）
3. 分簇散点图
4. 散点图 与 relplot()函数 (附时间序列图示例）
5. 散点矩阵图
6. 直方图矩阵
7. 密度图
8. 直方密度线图
9. 热力图及半角热力图（半三角）
10. 平行坐标图
11. 多变量联合分布图 -- pairplot() 函数
12. 多组分类重叠密度图 -- joy plot() 函数
13. 联合分布图 -- jointplot() 函数
第三部分：统计分析图表及回归图
1. 检验是否符合正态分布 -- p_test,
Skewness(), Kurtosis()的计算
2. 正态概率分布图
3. 正态分布图 -- norm.pdf() 函数的计算与绘制
4. 回归图 -- lmplot() 函数
5. 回归图 -- regplot() 函数
6. 回归图 -- lmplot() 函数
示例代码（Python代码）：第一部分：鸢尾花数据集的获取与转换**鸢尾花数据集（Iris数据集）**是一类多重变量分析的数据集。数据集包含150个数据样本，分为3类，每类50个数据，每个数据包含4个属性。通过花萼长度，花萼宽度，花瓣长度，花瓣宽度4个属性可以预测鸢尾花卉属于（Setosa，Versicolour，Virginica）三个种类中的哪一类。原数据集的目标值（iris.target）是“0”和“1”，但在数据可视化展示的时候，为了更清晰容易地理解数据的分类，我们需要将数据集中的目标值与鸢尾花的品种连接起来，即：把数据集中的编号"iris.target"与鸢尾花的品种连接，转换为新的数据集。import numpy as np
import pandas as pd
import sklearn
from sklearn import datasets
# 从sklearn中获取鸢尾花数据集，并转换为DataFrame
iris = datasets.load_iris()
dataset = pd.DataFrame(iris.data, columns=iris.feature_names)
# 或：更新特征列的名称，如“sepal_length”
# columns = [x.strip("(cm)").strip().replace(" ", "_") for x in iris.feature_names]
# dataset = pd.DataFrame(iris.data, columns=columns)
dataset["target"]= iris.target
# 对数据集中的目标值进行鸢尾花品种的转换
dict_species = dict(zip(np.array([0, 1, 2]), iris.target_names,))
# dict_species
dataset["speices"] = dataset["target"].map(dict_species)
# 将整理好的数据集以csv文件的形式保存下来，也可以保存为Excel 文件
outputfile = r"d:/iris.csv"
dataset.to_csv(outputfile)
dataset.info()
保存下来的csv文件第二部分：关系图的绘制1. 箱型图# 1.1 绘制箱型图 -- 根据不同类别的数据绘制箱型图
data = dataset.drop(columns=["species"])
plt.boxplot(data, labels=data.columns, showmeans=True)
plt.show()
# 1. 2 绘制箱型图 -- 按鸢尾花的品种展示花萼长度和花瓣长度
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
sns.boxplot(x="species", y="sepal_length", data=dataset, ax=axes[0])
sns.boxplot(x="species", y="petal_length", data=dataset, ax=axes[1])
plt.show()
2. 小提琴图 – violinplot() 函数# 绘制小提琴图
import seaborn as sns
sns.violinplot(x="species", y="petal_length", data=dataset)
plt.show()
3. 分簇散点图# 分簇散点图
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
sns.swarmplot(x="species", y="sepal_length", data=dataset, ax=axes[0])
sns.swarmplot(x="species", y="sepal_width", data=dataset, ax=axes[1])
plt.show()
4. 散点图# 4.1 绘制散点图
fig = plt.subplots(1, 1, figsize=(8, 4))
sns.scatterplot(x="sepal_length", y="petal_length", hue="species",
data=dataset)
plt.show()
# 4.2 散点图 -- 以颜色条区分类别
plt.style.use("ggplot")
plt.figure(figsize=(7, 5))
cmap = plt.cm.get_cmap("RdBu")
sc = plt.scatter(x=dataset["sepal_length"], y=dataset["sepal_width"], c=iris.target, s=20, cmap=cmap)
bar = plt.colorbar(sc)
bar.set_label("species")
plt.xlabel("sepal_length")
plt.ylabel("sepal_width")
plt.show()
# 4.3 relplot()函数 -- 设置kind="line"
import seaborn as sns
sns.set(style="ticks", palette="colorblind", color_codes=True)
sns.relplot(x = "sepal_length", y="sepal_width", data=dataset, kind="line")
plt.title("speal_length vs. speal_width")
plt.show()
# 4.4 relplot()函数 -- 设置kind="scatter"
sns.relplot(x ="sepal_length", y="sepal_width", data=dataset)
plt.title("speal_length vs. speal_width (scatter)")
plt.show()
# 4.5 relplot()函数 -- 设置kind="scatter",并通过hue参数做类别区分
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.relplot(x="sepal_length", y="sepal_width", hue="species",
size="species", sizes=(50, 120), style="species", data=dataset)
plt.title("speal_length vs. speal_width (by species)")
plt.show()
# 4.6 replot()函数绘制时间序列图的示例（本示例是样品的序列号）
sns.set(style="whitegrid", palette="muted", color_codes=True)
sns.relplot(x=dataset.index, y="sepal_length", data=dataset,
hue="species", style="species", kind="line")
plt.title("sepal_length by index")
plt.show()
5. 散点矩阵图# 5.1 绘制散点图矩阵示例
import seaborn as sns
sns.set(style="ticks")
sns.pairplot(dataset, vars=["sepal_length", "petal_length"])
plt.show()
# 5.2 设定hue参数以区分不同的种类
import seaborn as sns
sns.set(style="ticks")
sns.pairplot(dataset, vars=["sepal_length", "petal_length"], hue="species")
plt.show()
6. 直方图矩阵# 绘制直方图矩阵 -- 对数据集中的所有特征属性绘制直方图
dataset.hist(sharex=True)
plt.show()
7. 绘制密度图# 绘制密度图 -- 整个数据集
dataset.plot(kind="kde")
plt.show()
8. 直方密度线图#
绘制直方密度线图 -- 并按不同品种分别呈现sepal_length数据的分布情况
import seaborn as sns
# kde: 是否显示数据分布曲线，默认为False
sns.displot(x="sepal_length", data=dataset, bins=20, kde=True, hue="species")
plt.show()
9. 热力图以及半角热力图# 9.1 绘制热力图
data = dataset.drop(columns=["target"])
corrmat = data.corr()
k = 4
# 排序：根据相关性程度从大到小进行排序（选定1个特征属性作为对照）
cols = corrmat.nlargest(k, "sepal_length")["sepal_length"].index
cm = np.corrcoef(data[cols].values.T)
sns.set(font_scale=1.25)
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt=".2f", annot_kws={"size": 10},
yticklabels=cols.values, xticklabels=cols.values)
plt.show()
# 9.2 构建半角热力图：方法一
corrmat = data.corr()
mask = np.ones_like(corrmat)
num = 4
mask[np.tril_indices(num)] = 0.
# 上三角被屏蔽
# mask[np.triu_indices(num)] = 0.
# 下三角被屏蔽
sns.heatmap(corrmat, cbar=True, square=True, fmt=".2f", mask=mask,
annot=True, yticklabels=cols.values, xticklabels=cols.values)
plt.show()
# 9.3 构建半角热力图：方法二
corrmat = data.corr()
mask = np.zeros_like(corrmat)
# 将mask的对角线及以上设置为True
mask[np.triu_indices_from(mask)] = True
sns.heatmap(corrmat, mask=mask, square=True, annot=True, fmt="0.2f")
plt.show()
10. 平行坐标图平行坐标图（Paralllel Coordinates Plot) 是对于具有多个属性的一种可视化方法，可解决在维度增加时，散点矩阵变得不太有效的问题。在平行坐标图中，数据集的一行数据在平行坐标图中用一条折线表示，纵向是属性，横向是属性类别。安装方法：pip install pyecharts# 绘制数据的平行坐标图示例
from pyecharts.charts import Parallel
import pyecharts.options as opts
import seaborn as sns
import numpy as np
data_ = np.array(dataset[["sepal_length", "sepal_width", "petal_length", "petal_width"]]).tolist()
parallel_axis = [{"dim": 0, "name": "萼片长度"},
{"dim": 1, "name": "萼片宽度"},
{"dim": 2, "name": "花瓣长度"},
{"dim": 3, "name": "花瓣宽度"}]
parallel = Parallel(init_opts=opts.InitOpts(width="600px", height="400px"))
parallel.add_schema(schema=parallel_axis)
parallel.add("鸢尾花（iris）的平行图", data=data_, linestyle_opts=opts.LineStyleOpts(width=4, opacity=0.5))
parallel.render_notebook()
11. 多变量联合分布图 – pairplot() 函数11.1 以散点图的形式循环展示数据属性之间的相关性# 11.1：绘制以散点图的形式展示数据属性之间的相关性
import seaborn as sns
plt.figure(figsize=(10, 8), dpi=80)
plot_setting = dict(s=80, edgecolor="white", linewidth=2.5)
sns.pairplot(dataset, kind="scatter", hue="species",
plot_kws=plot_setting)
plt.show()
11.2 以回归线的方式循环展示数据属性之间的相关性# 11.2：绘制以回归线的方式循环展示数据属性之间的相关性
plt.figure(figsize=(10, 8), dpi=80)
sns.pairplot(dataset, kind="reg", hue="species")
plt.show()
# 11.3 绘制2个特征属性的关系，并以颜色和符合区分不同的类别
plt.style.use("ggplot")
sns.pairplot(data=dataset[["petal_length", "petal_width", "species"]],
hue="species", markers=["o", "*", "^"])
plt.show()
12. 多组分类重叠密度图 – Joyplot() 函数多组分类重叠密度图（Joy plot）又称为“峰峦图”，是一种可视化大量分组数据的方法，通过部分堆叠、重叠的密度图来展示不同类别的密度曲线折叠状况，直观地在一个维度上呈现和比较不同组别数据的分布。安装方法：pip install joyplot# Joy Plot
import joypy
plt.figure(figsize=(10, 6), dpi=80)
fig, axes = joypy.joyplot(dataset, column=["sepal_length", "sepal_width"], by="species", figsize=(10, 6),
grid=True, title="Sepal_length vs. Sepal_width")
plt.show()
13. 联合分布图 – jointplot() 函数# 13.1 联合分布图 -- jointplot()
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.jointplot(x="sepal_length", y="sepal_width", data=dataset)
plt.show()
# 13.2 联合分布图 使用hue参数来区分不同的类别
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.jointplot(x="petal_length", y="petal_width", hue="species", data=dataset)
plt.show()
第三部分：统计分析图表及回归图1. 检验是否符合正态分布 – p_test, Skewness(), Kurtosis()的计算# 1.1 检查是否属于正态分布及偏度（Skewness）和峰度(Kurtosis)
print("偏度(Skewness): %f" % dataset["sepal_length"].skew())
print("峰度(Kurtosis): %f" % dataset["sepal_length"].kurt())
# 在统计学中，峰度（Kurtosis）衡量实数随机变量概率分布的峰态。
# 峰度高就意味着方差增大是由低频度的大于或小于平均值的极端差值引起的。
# 1.2检验数据的正态分布及p_test 的计算 -- 分析petal_length是否符合正态分布：
import scipy.stats as ss
p_test= np.array(dataset["sepal_length"].T)
print(ss.normaltest(p_test))
from matplotlib import pyplot as plt
p_test = pd.Series(p_test)
p_test.plot(kind="kde")
# 2.2 分析petal_length是否符合正态分布：
import scipy.stats as ss
data_ = dataset[dataset["species"] == "setosa"]
p_test= np.array(data_["sepal_length"].T)
print(ss.normaltest(p_test))
from matplotlib import pyplot as plt
p_test = pd.Series(p_test)
p_test.plot(kind="kde")
2. 正态概率分布图# 正态概率分布图：-- histogram and normal probability plot
from scipy import stats
sns.distplot(dataset["sepal_length"], fit=norm)
fig = plt.figure()
res = stats.probplot(dataset["sepal_length"], plot=plt)
3. 正态分布图 – norm.pdf() 函数的计算与绘制# 3.1 绘制正态分布图
from scipy.stats import norm
fig, axes = plt.subplots()
sigma = dataset["sepal_length"].std()
mu = dataset["sepal_length"].mean()
num_bins = 20
x = dataset["sepal_length"]
n, bins, patches = axes.hist(x, num_bins, density=1)
# 计算正态分布概率密度函数
y = norm.pdf(bins, mu, sigma)
axes.plot(bins, y, "r--")
axes.set_title("sepal_length(cm)的正态分布图")
fig.tight_layout()
plt.show()
# 3.2 绘制正态分布图 -- 关注于某种品种的花瓣长度的数据分析
mpl.rcParams["font.family"] = "SimHei"
data_ = dataset[dataset["species"]=="setosa"]
from scipy.stats import norm
fig, axes = plt.subplots()
sigma = data_["sepal_length"].std()
mu = data_["sepal_length"].mean()
num_bins = 20
x = data_["sepal_length"]
n, bins, patches = axes.hist(x, num_bins, density=1,color="g")
y = norm.pdf(bins, mu, sigma)
axes.plot(bins, y, "r--")
axes.set_title("Setosa: sepal_length(cm)的正态分布图")
fig.tight_layout()
plt.show()
4. 回归图 – lmplot()函数# 绘制回归图
import seaborn as sns
sns.lmplot(x="sepal_length", y="petal_length", hue="species", data=dataset)
plt.show()
5. 回归图 – 使用regplot()函数# 绘制线性回归图 -- 分别绘制
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
sns.regplot(x="petal_length", y="petal_width", data=dataset, color="g", ax=axes[0])
sns.regplot(x="sepal_length", y="sepal_width", data=dataset, color="orange", ax=axes[1])
plt.show()
6. 回归图 – 使用lmplot()函数# 使用lmplot()函数绘制和呈现不同类别的两两特征属性的相关情况
plt.style.use("ggplot")
sns.lmplot(x="sepal_length", y="sepal_width", hue="species",
data=dataset, fit_reg=False, size=5)
sns.lmplot(x="sepal_length", y="sepal_width", col="species", hue="species",
data=dataset, fit_reg=True, size=3, aspect=0.9, col_wrap=3)
plt.show()
附录：1. Joyplot() 函数的介绍‘data：绘制数据集’‘column’：使用data的中的有限列进行绘图‘by=None’:分组列‘gird=false：添加网格线‘xlabelsize=none x轴标签的大小‘ylabelsize=none y轴标签的大小‘xrot=none x轴刻度线标签旋转角度‘yrot=none y轴刻度线标签旋转角度‘hist=flase显示直方图‘fade=flase如果设定的是true，则显示渐变色‘ylim’=‘max共享y轴的刻度ll=‘true 曲线下的填充颜色linecolor=‘None；曲线的颜色blackground=none:背景颜色overlap=1:控制重叠程度‘title’=none 添加图表的标题‘colormap=none 色谱}
1. 箱型图
2. 小提琴图 -- violinplot(）
3. 分簇散点图
4. 散点图 与 relplot()函数 (附时间序列图示例）
5. 散点矩阵图
6. 直方图矩阵
7. 密度图
8. 直方密度线图
9. 热力图及半角热力图（半三角）
10. 平行坐标图
11. 多变量联合分布图 -- pairplot() 函数
12. 多组分类重叠密度图 -- joy plot() 函数
13. 联合分布图 -- jointplot() 函数
第三部分：统计分析图表及回归图
1. 检验是否符合正态分布 -- p_test,
Skewness(), Kurtosis()的计算
2. 正态概率分布图
3. 正态分布图 -- norm.pdf() 函数的计算与绘制
4. 回归图 -- lmplot() 函数
5. 回归图 -- regplot() 函数
6. 回归图 -- lmplot() 函数
示例代码（Python代码）：第一部分：鸢尾花数据集的获取与转换**鸢尾花数据集（Iris数据集）**是一类多重变量分析的数据集。数据集包含150个数据样本，分为3类，每类50个数据，每个数据包含4个属性。通过花萼长度，花萼宽度，花瓣长度，花瓣宽度4个属性可以预测鸢尾花卉属于（Setosa，Versicolour，Virginica）三个种类中的哪一类。原数据集的目标值（iris.target）是“0”和“1”，但在数据可视化展示的时候，为了更清晰容易地理解数据的分类，我们需要将数据集中的目标值与鸢尾花的品种连接起来，即：把数据集中的编号"iris.target"与鸢尾花的品种连接，转换为新的数据集。import numpy as np
import pandas as pd
import sklearn
from sklearn import datasets
# 从sklearn中获取鸢尾花数据集，并转换为DataFrame
iris = datasets.load_iris()
dataset = pd.DataFrame(iris.data, columns=iris.feature_names)
# 或：更新特征列的名称，如“sepal_length”
# columns = [x.strip("(cm)").strip().replace(" ", "_") for x in iris.feature_names]
# dataset = pd.DataFrame(iris.data, columns=columns)
dataset["target"]= iris.target
# 对数据集中的目标值进行鸢尾花品种的转换
dict_species = dict(zip(np.array([0, 1, 2]), iris.target_names,))
# dict_species
dataset["speices"] = dataset["target"].map(dict_species)
# 将整理好的数据集以csv文件的形式保存下来，也可以保存为Excel 文件
outputfile = r"d:/iris.csv"
dataset.to_csv(outputfile)
dataset.info()
保存下来的csv文件第二部分：关系图的绘制1. 箱型图# 1.1 绘制箱型图 -- 根据不同类别的数据绘制箱型图
data = dataset.drop(columns=["species"])
plt.boxplot(data, labels=data.columns, showmeans=True)
plt.show()
# 1. 2 绘制箱型图 -- 按鸢尾花的品种展示花萼长度和花瓣长度
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
sns.boxplot(x="species", y="sepal_length", data=dataset, ax=axes[0])
sns.boxplot(x="species", y="petal_length", data=dataset, ax=axes[1])
plt.show()
2. 小提琴图 – violinplot() 函数# 绘制小提琴图
import seaborn as sns
sns.violinplot(x="species", y="petal_length", data=dataset)
plt.show()
3. 分簇散点图# 分簇散点图
fig, axes = plt.subplots(1, 2, figsize=(10, 5))
sns.swarmplot(x="species", y="sepal_length", data=dataset, ax=axes[0])
sns.swarmplot(x="species", y="sepal_width", data=dataset, ax=axes[1])
plt.show()
4. 散点图# 4.1 绘制散点图
fig = plt.subplots(1, 1, figsize=(8, 4))
sns.scatterplot(x="sepal_length", y="petal_length", hue="species",
data=dataset)
plt.show()
# 4.2 散点图 -- 以颜色条区分类别
plt.style.use("ggplot")
plt.figure(figsize=(7, 5))
cmap = plt.cm.get_cmap("RdBu")
sc = plt.scatter(x=dataset["sepal_length"], y=dataset["sepal_width"], c=iris.target, s=20, cmap=cmap)
bar = plt.colorbar(sc)
bar.set_label("species")
plt.xlabel("sepal_length")
plt.ylabel("sepal_width")
plt.show()
# 4.3 relplot()函数 -- 设置kind="line"
import seaborn as sns
sns.set(style="ticks", palette="colorblind", color_codes=True)
sns.relplot(x = "sepal_length", y="sepal_width", data=dataset, kind="line")
plt.title("speal_length vs. speal_width")
plt.show()
# 4.4 relplot()函数 -- 设置kind="scatter"
sns.relplot(x ="sepal_length", y="sepal_width", data=dataset)
plt.title("speal_length vs. speal_width (scatter)")
plt.show()
# 4.5 relplot()函数 -- 设置kind="scatter",并通过hue参数做类别区分
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.relplot(x="sepal_length", y="sepal_width", hue="species",
size="species", sizes=(50, 120), style="species", data=dataset)
plt.title("speal_length vs. speal_width (by species)")
plt.show()
# 4.6 replot()函数绘制时间序列图的示例（本示例是样品的序列号）
sns.set(style="whitegrid", palette="muted", color_codes=True)
sns.relplot(x=dataset.index, y="sepal_length", data=dataset,
hue="species", style="species", kind="line")
plt.title("sepal_length by index")
plt.show()
5. 散点矩阵图# 5.1 绘制散点图矩阵示例
import seaborn as sns
sns.set(style="ticks")
sns.pairplot(dataset, vars=["sepal_length", "petal_length"])
plt.show()
# 5.2 设定hue参数以区分不同的种类
import seaborn as sns
sns.set(style="ticks")
sns.pairplot(dataset, vars=["sepal_length", "petal_length"], hue="species")
plt.show()
6. 直方图矩阵# 绘制直方图矩阵 -- 对数据集中的所有特征属性绘制直方图
dataset.hist(sharex=True)
plt.show()
7. 绘制密度图# 绘制密度图 -- 整个数据集
dataset.plot(kind="kde")
plt.show()
8. 直方密度线图#
绘制直方密度线图 -- 并按不同品种分别呈现sepal_length数据的分布情况
import seaborn as sns
# kde: 是否显示数据分布曲线，默认为False
sns.displot(x="sepal_length", data=dataset, bins=20, kde=True, hue="species")
plt.show()
9. 热力图以及半角热力图# 9.1 绘制热力图
data = dataset.drop(columns=["target"])
corrmat = data.corr()
k = 4
# 排序：根据相关性程度从大到小进行排序（选定1个特征属性作为对照）
cols = corrmat.nlargest(k, "sepal_length")["sepal_length"].index
cm = np.corrcoef(data[cols].values.T)
sns.set(font_scale=1.25)
hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt=".2f", annot_kws={"size": 10},
yticklabels=cols.values, xticklabels=cols.values)
plt.show()
# 9.2 构建半角热力图：方法一
corrmat = data.corr()
mask = np.ones_like(corrmat)
num = 4
mask[np.tril_indices(num)] = 0.
# 上三角被屏蔽
# mask[np.triu_indices(num)] = 0.
# 下三角被屏蔽
sns.heatmap(corrmat, cbar=True, square=True, fmt=".2f", mask=mask,
annot=True, yticklabels=cols.values, xticklabels=cols.values)
plt.show()
# 9.3 构建半角热力图：方法二
corrmat = data.corr()
mask = np.zeros_like(corrmat)
# 将mask的对角线及以上设置为True
mask[np.triu_indices_from(mask)] = True
sns.heatmap(corrmat, mask=mask, square=True, annot=True, fmt="0.2f")
plt.show()
10. 平行坐标图平行坐标图（Paralllel Coordinates Plot) 是对于具有多个属性的一种可视化方法，可解决在维度增加时，散点矩阵变得不太有效的问题。在平行坐标图中，数据集的一行数据在平行坐标图中用一条折线表示，纵向是属性，横向是属性类别。安装方法：pip install pyecharts# 绘制数据的平行坐标图示例
from pyecharts.charts import Parallel
import pyecharts.options as opts
import seaborn as sns
import numpy as np
data_ = np.array(dataset[["sepal_length", "sepal_width", "petal_length", "petal_width"]]).tolist()
parallel_axis = [{"dim": 0, "name": "萼片长度"},
{"dim": 1, "name": "萼片宽度"},
{"dim": 2, "name": "花瓣长度"},
{"dim": 3, "name": "花瓣宽度"}]
parallel = Parallel(init_opts=opts.InitOpts(width="600px", height="400px"))
parallel.add_schema(schema=parallel_axis)
parallel.add("鸢尾花（iris）的平行图", data=data_, linestyle_opts=opts.LineStyleOpts(width=4, opacity=0.5))
parallel.render_notebook()
11. 多变量联合分布图 – pairplot() 函数11.1 以散点图的形式循环展示数据属性之间的相关性# 11.1：绘制以散点图的形式展示数据属性之间的相关性
import seaborn as sns
plt.figure(figsize=(10, 8), dpi=80)
plot_setting = dict(s=80, edgecolor="white", linewidth=2.5)
sns.pairplot(dataset, kind="scatter", hue="species",
plot_kws=plot_setting)
plt.show()
11.2 以回归线的方式循环展示数据属性之间的相关性# 11.2：绘制以回归线的方式循环展示数据属性之间的相关性
plt.figure(figsize=(10, 8), dpi=80)
sns.pairplot(dataset, kind="reg", hue="species")
plt.show()
# 11.3 绘制2个特征属性的关系，并以颜色和符合区分不同的类别
plt.style.use("ggplot")
sns.pairplot(data=dataset[["petal_length", "petal_width", "species"]],
hue="species", markers=["o", "*", "^"])
plt.show()
12. 多组分类重叠密度图 – Joyplot() 函数多组分类重叠密度图（Joy plot）又称为“峰峦图”，是一种可视化大量分组数据的方法，通过部分堆叠、重叠的密度图来展示不同类别的密度曲线折叠状况，直观地在一个维度上呈现和比较不同组别数据的分布。安装方法：pip install joyplot# Joy Plot
import joypy
plt.figure(figsize=(10, 6), dpi=80)
fig, axes = joypy.joyplot(dataset, column=["sepal_length", "sepal_width"], by="species", figsize=(10, 6),
grid=True, title="Sepal_length vs. Sepal_width")
plt.show()
13. 联合分布图 – jointplot() 函数# 13.1 联合分布图 -- jointplot()
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.jointplot(x="sepal_length", y="sepal_width", data=dataset)
plt.show()
# 13.2 联合分布图 使用hue参数来区分不同的类别
sns.set(style="darkgrid", palette="muted", color_codes=True)
sns.jointplot(x="petal_length", y="petal_width", hue="species", data=dataset)
plt.show()
第三部分：统计分析图表及回归图1. 检验是否符合正态分布 – p_test, Skewness(), Kurtosis()的计算# 1.1 检查是否属于正态分布及偏度（Skewness）和峰度(Kurtosis)
print("偏度(Skewness): %f" % dataset["sepal_length"].skew())
print("峰度(Kurtosis): %f" % dataset["sepal_length"].kurt())
# 在统计学中，峰度（Kurtosis）衡量实数随机变量概率分布的峰态。
# 峰度高就意味着方差增大是由低频度的大于或小于平均值的极端差值引起的。
# 1.2检验数据的正态分布及p_test 的计算 -- 分析petal_length是否符合正态分布：
import scipy.stats as ss
p_test= np.array(dataset["sepal_length"].T)
print(ss.normaltest(p_test))
from matplotlib import pyplot as plt
p_test = pd.Series(p_test)
p_test.plot(kind="kde")
# 2.2 分析petal_length是否符合正态分布：
import scipy.stats as ss
data_ = dataset[dataset["species"] == "setosa"]
p_test= np.array(data_["sepal_length"].T)
print(ss.normaltest(p_test))
from matplotlib import pyplot as plt
p_test = pd.Series(p_test)
p_test.plot(kind="kde")
2. 正态概率分布图# 正态概率分布图：-- histogram and normal probability plot
from scipy import stats
sns.distplot(dataset["sepal_length"], fit=norm)
fig = plt.figure()
res = stats.probplot(dataset["sepal_length"], plot=plt)
3. 正态分布图 – norm.pdf() 函数的计算与绘制# 3.1 绘制正态分布图
from scipy.stats import norm
fig, axes = plt.subplots()
sigma = dataset["sepal_length"].std()
mu = dataset["sepal_length"].mean()
num_bins = 20
x = dataset["sepal_length"]
n, bins, patches = axes.hist(x, num_bins, density=1)
# 计算正态分布概率密度函数
y = norm.pdf(bins, mu, sigma)
axes.plot(bins, y, "r--")
axes.set_title("sepal_length(cm)的正态分布图")
fig.tight_layout()
plt.show()
# 3.2 绘制正态分布图 -- 关注于某种品种的花瓣长度的数据分析
mpl.rcParams["font.family"] = "SimHei"
data_ = dataset[dataset["species"]=="setosa"]
from scipy.stats import norm
fig, axes = plt.subplots()
sigma = data_["sepal_length"].std()
mu = data_["sepal_length"].mean()
num_bins = 20
x = data_["sepal_length"]
n, bins, patches = axes.hist(x, num_bins, density=1,color="g")
y = norm.pdf(bins, mu, sigma)
axes.plot(bins, y, "r--")
axes.set_title("Setosa: sepal_length(cm)的正态分布图")
fig.tight_layout()
plt.show()
4. 回归图 – lmplot()函数# 绘制回归图
import seaborn as sns
sns.lmplot(x="sepal_length", y="petal_length", hue="species", data=dataset)
plt.show()
5. 回归图 – 使用regplot()函数# 绘制线性回归图 -- 分别绘制
fig, axes = plt.subplots(1, 2, figsize=(8, 4))
sns.regplot(x="petal_length", y="petal_width", data=dataset, color="g", ax=axes[0])
sns.regplot(x="sepal_length", y="sepal_width", data=dataset, color="orange", ax=axes[1])
plt.show()
6. 回归图 – 使用lmplot()函数# 使用lmplot()函数绘制和呈现不同类别的两两特征属性的相关情况
plt.style.use("ggplot")
sns.lmplot(x="sepal_length", y="sepal_width", hue="species",
data=dataset, fit_reg=False, size=5)
sns.lmplot(x="sepal_length", y="sepal_width", col="species", hue="species",
data=dataset, fit_reg=True, size=3, aspect=0.9, col_wrap=3)
plt.show()
附录：1. Joyplot() 函数的介绍‘data：绘制数据集’‘column’：使用data的中的有限列进行绘图‘by=None’:分组列‘gird=false：添加网格线‘xlabelsize=none x轴标签的大小‘ylabelsize=none y轴标签的大小‘xrot=none x轴刻度线标签旋转角度‘yrot=none y轴刻度线标签旋转角度‘hist=flase显示直方图‘fade=flase如果设定的是true，则显示渐变色‘ylim’=‘max共享y轴的刻度ll=‘true 曲线下的填充颜色linecolor=‘None；曲线的颜色blackground=none:背景颜色overlap=1:控制重叠程度‘title’=none 添加图表的标题‘colormap=none 色谱}
我要回帖
更多关于 二阶原点矩的数学期望和方差 的文章
更多推荐
- ·怎么在excel重复值中展示+号/-号？
- ·看100一光年要飞多久之外的天体需要多久，可望不可及还是不渴望不可及！
- ·新疆乌鲁木齐2月份平均气温5月份初的气温是多少度?
- ·沁县oppo手机自带天气软件下载山西沁县oppo手机自带天气软件下载预报一周oppo手机自带天气软件下载
- ·72V电动车72v充电器不充电怎么修电路板这个位置应该是一个什么IC芯片呢？
- ·pat+—(have)at+at the school和at school的区别?
- ·等于1的1是带分数还是假分数，分数，小数，假分数，整数？
- ·积相等的三两个数和一定,差越小积越大为什么3两个数和一定,差越小积越大越接近,它们的和是最小的？
- ·请问国产半导体还有机会吗发展吗？
- ·沈桃沈康高速连接线图编号是什么？
- ·二分之一+2的二次方之一+2+³之一++……++2(-2)的2019次方+(-2)的2020次方分一=?
- ·关于在智能制造中,你对数据安全的理解的硬件设备的安全性，请问怎么保证？
- ·计算9.2.7的样本二阶原点矩方差怎么算三阶中心矩三阶偏度？
- ·一氧化碳是不是刺激性气体属 ()气体A)刺激性B 窒息性C)以上皆有？
- ·华英证券怎么样销售能力如何？
- ·华英证券排名的业务团队可靠吗？
- ·华英证券排名的债券业绩如何？
- ·华英华菁证券是小券商吗规模大吗？
- ·华英证券排名的债券业绩如何？
- ·秋田满满是哪个国家的里的企业？
