精通 Python 金融第二版（三）
原文：
zh.annas-archive.org/md5/8b046e39ce2c1a10ac13fd89834aaadc译者：飞龙
第六章：时间序列数据的统计分析
在金融投资组合中，其组成资产的回报取决于许多因素，如宏观和微观经济条件以及各种金融变量。随着因素数量的增加，建模投资组合行为所涉及的复杂性也在增加。鉴于计算资源是有限的，再加上时间限制，为新因素进行额外计算只会增加投资组合建模计算的瓶颈。一种用于降维的线性技术是主成分分析（PCA）。正如其名称所示，PCA 将投资组合资产价格的变动分解为其主要成分或共同因素，以进行进一步的统计分析。不能解释投资组合资产变动很多的共同因素在其因素中获得较少的权重，并且通常被忽略。通过保留最有用的因素，可以大大简化投资组合分析，而不会影响计算时间和空间成本。
在时间序列数据的统计分析中，数据保持平稳对于避免虚假回归是很重要的。非平稳数据可能由受趋势影响的基础过程、季节效应、单位根的存在或三者的组合产生。非平稳数据的统计特性，如均值和方差，会随时间变化。非平稳数据需要转换为平稳数据，以便进行统计分析以产生一致和可靠的结果。这可以通过去除趋势和季节性成分来实现。然后可以使用平稳数据进行预测或预测。
在本章中，我们将涵盖以下主题：
-
对道琼斯及其 30 个成分进行主成分分析
-
重建道琼斯指数
-
理解平稳和非平稳数据之间的区别
-
检查数据的平稳性
-
平稳和非平稳过程的类型
-
使用增广迪基-富勒检验来检验单位根的存在
-
通过去趋势化、差分和季节性分解制作平稳数据
-
使用自回归积分移动平均法进行时间序列预测和预测
道琼斯工业平均指数及其 30 个成分
道琼斯工业平均指数（DJIA）是由 30 家最大的美国公司组成的股票市场指数。通常被称为道琼斯，它由 S＆P 道琼斯指数有限责任公司拥有，并以价格加权的方式计算（有关道琼斯的更多信息，请参见us.spindices.com/index-family/us-equity/dow-jones-averages）。
本节涉及将道琼斯及其成分的数据集下载到pandas DataFrame 对象中，以供本章后续部分使用。
从 Quandl 下载道琼斯成分数据集
以下代码从 Quandl 检索道琼斯成分数据集。我们将使用的数据提供者是 WIKI Prices，这是一个由公众成员组成的社区，向公众免费提供数据集。这些数据并非没有错误，因此请谨慎使用。在撰写本文时，该数据源不再受到 Quandl 社区的积极支持，尽管过去的数据集仍可供使用。我们将下载 2017 年的历史每日收盘价：
In [ ]:
import quandl
QUANDL_API_KEY = 'BCzkk3NDWt7H9yjzx-DY' # Your own Quandl key here
quandl.ApiConfig.api_key = QUANDL_API_KEY
SYMBOLS = [
'AAPL','MMM', 'AXP', 'BA', 'CAT',
'CVX', 'CSCO', 'KO', 'DD', 'XOM',
'GS', 'HD', 'IBM', 'INTC', 'JNJ',
'JPM', 'MCD', 'MRK', 'MSFT', 'NKE',
'PFE', 'PG', 'UNH', 'UTX', 'TRV',
'VZ', 'V', 'WMT', 'WBA', 'DIS',
]
wiki_symbols = ['WIKI/%s'%symbol for symbol in SYMBOLS]
df_components = quandl.get(
wiki_symbols,
start_date='2017-01-01',
end_date='2017-12-31',
column_index=11)
df_components.columns = SYMBOLS # Renaming the columns
wiki_symbols变量包含我们用于下载的 Quandl 代码列表。请注意，在quandl.get()的参数中，我们指定了column_index=11。这告诉 Quandl 仅下载每个数据集的第 11 列，这与调整后的每日收盘价相符。数据集以单个pandas DataFrame 对象的形式下载到我们的df_components变量中。
让我们在分析之前对数据集进行归一化处理：
In [ ]:
filled_df_components = df_components.fillna(method='ffill')
daily_df_components = filled_df_components.resample('24h').ffill()
daily_df_components = daily_df_components.fillna(method='bfill')
如果您检查这个数据源中的每个值，您会注意到NaN值或缺失数据。由于我们使用的是容易出错的数据，并且为了快速研究 PCA，我们可以通过传播先前观察到的值临时填充这些未知变量。fillna(method='ffill')方法有助于执行此操作，并将结果存储在filled_df_components变量中。
标准化的另一个步骤是以固定间隔重新取样时间序列，并将其与我们稍后将要下载的道琼斯时间序列数据集完全匹配。daily_df_components变量存储了按日重新取样时间序列的结果，重新取样期间的任何缺失值都使用向前填充方法传播。最后，为了解决起始数据不完整的问题，我们将简单地使用fillna(method='bfill')对值进行回填。
为了展示 PCA 的目的，我们必须使用免费的低质量数据集。如果您需要高质量的数据集，请考虑订阅数据发布商。
Quandl 不提供道琼斯工业平均指数的免费数据集。在下一节中，我们将探索另一个名为 Alpha Vantage 的数据提供商，作为下载数据集的替代方法。
关于 Alpha Vantage
Alpha Vantage (www.alphavantage.co)是一个数据提供商，提供股票、外汇和加密货币的实时和历史数据。与 Quandl 类似，您可以获得 Alpha Vantage REST API 接口的 Python 包装器，并直接将免费数据集下载到pandas DataFrame 中。
获取 Alpha Vantage API 密钥
从您的网络浏览器访问www.alphavantage.co，并从主页点击立即获取您的免费 API 密钥。您将被带到注册页面。填写关于您自己的基本信息并提交表单。您的 API 密钥将显示在同一页。复制此 API 密钥以在下一节中使用：
安装 Alpha Vantage Python 包装器
从您的终端窗口，输入以下命令以安装 Alpha Vantage 的 Python 模块：
$ pip install alpha_vantage
从 Alpha Vantage 下载道琼斯数据集
以下代码连接到 Alpha Vantage 并下载道琼斯数据集，股票代码为^DJI。用您自己的 API 密钥替换常量变量ALPHA_VANTAGE_API_KEY的值：
In [ ]:
"""
Download the all-time DJIA dataset
"""
from alpha_vantage.timeseries import TimeSeries
# Update your Alpha Vantage API key here...
ALPHA_VANTAGE_API_KEY = 'PZ2ISG9CYY379KLI'
ts = TimeSeries(key=ALPHA_VANTAGE_API_KEY, output_format='pandas')
df, meta_data = ts.get_daily_adjusted(symbol='^DJI', outputsize='full')
alpha_vantage.timeseries模块的TimeSeries类是用 API 密钥实例化的，并指定数据集自动下载为pandas DataFrame 对象。get_daily_adjusted()方法使用outputsize='full'参数下载给定股票符号的整个可用每日调整价格，并将其存储在df变量中作为DataFrame对象。
让我们使用info()命令检查一下这个 DataFrame：
In [ ]:
df.info()
Out[ ]:
<class 'pandas.core.frame.DataFrame'>
Index: 4760 entries, 2000-01-03 to 2018-11-30
Data columns (total 8 columns):
1\. open 4760 non-null float64
2\. high 4760 non-null float64
3\. low 4760 non-null float64
4\. close 4760 non-null float64
5\. adjusted close 4760 non-null float64
6\. volume 4760 non-null float64
7\. dividend amount 4760 non-null float64
8\. split coefficient 4760 non-null float64
dtypes: float64(8)
memory usage: 316.1+ KB
我们从 Alpha Vantage 下载的道琼斯数据集提供了从最近可用交易日期一直回到 2000 年的完整时间序列数据。它包含几列给我们额外的信息。
让我们也检查一下这个 DataFrame 的索引：
In [ ]:
df.index
Out[ ]:
Index(['2000-01-03', '2000-01-04', '2000-01-05', '2000-01-06', '2000-01-07',
'2000-01-10', '2000-01-11', '2000-01-12', '2000-01-13', '2000-01-14',
...
'2018-08-17', '2018-08-20', '2018-08-21', '2018-08-22', '2018-08-23',
'2018-08-24', '2018-08-27', '2018-08-28', '2018-08-29', '2018-08-30'],
dtype='object', name='date', length=4696)
输出表明索引值由字符串类型的对象组成。让我们将这个 DataFrame 转换为适合我们分析的形式：
In [ ]:
import pandas as pd
# Prepare the dataframe
df_dji = pd.DataFrame(df['5\. adjusted close'])
df_dji.columns = ['DJIA']
df_dji.index = pd.to_datetime(df_dji.index)
# Trim the new dataframe and resample
djia_2017 = pd.DataFrame(df_dji.loc['2017-01-01':'2017-12-31'])
djia_2017 = djia_2017.resample('24h').ffill()
在这里，我们正在获取 2017 年道琼斯的调整收盘价，按日重新取样。结果的 DataFrame 对象存储在djia_2017中，我们可以用它来应用 PCA。
应用核 PCA
在本节中，我们将执行核 PCA 以找到特征向量和特征值，以便我们可以重建道琼斯指数。
寻找特征向量和特征值
我们可以使用 Python 的sklearn.decomposition模块的KernelPCA类执行核 PCA。默认的核方法是线性的。在 PCA 中使用的数据集需要被标准化，我们可以使用 z-scoring 来实现。以下代码执行此操作：
In [ ]:
from sklearn.decomposition import KernelPCA
fn_z_score = lambda x: (x - x.mean()) / x.std()
df_z_components = daily_df_components.apply(fn_z_score)
fitted_pca = KernelPCA().fit(df_z_components)
fn_z_score变量是一个内联函数，用于对pandas DataFrame 执行 z 得分，该函数使用apply()方法应用。这些归一化的数据集可以使用fit()方法拟合到核 PCA 中。每日道琼斯成分价格的拟合结果存储在fitted_pca变量中，该变量是相同的KernelPCA对象。
PCA 的两个主要输出是特征向量和特征值。特征向量是包含主成分线方向的向量，当应用线性变换时不会改变。特征值是标量值，指示数据在特定特征向量方向上的方差量。实际上，具有最高特征值的特征向量形成主成分。
KernelPCA对象的alphas_和lambdas_属性返回中心化核矩阵数据集的特征向量和特征值。当我们绘制特征值时，我们得到以下结果：
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (12,8)
plt.plot(fitted_pca.lambdas_)
plt.ylabel('Eigenvalues')
plt.show();
然后我们应该得到以下输出：
我们可以看到，前几个特征值解释了数据中的大部分方差，并且在后面的成分中变得更加忽略。获取前五个特征值，让我们看看每个特征值给我们提供了多少解释：
In [ ]:
fn_weighted_avg = lambda x: x / x.sum()
weighted_values = fn_weighted_avg(fitted_pca.lambdas_)[:5]
In [ ]:
print(weighted_values)
Out[ ]:
array([0.64863002, 0.13966718, 0.05558246, 0.05461861, 0.02313883])
我们可以看到，第一个成分解释了数据方差的 65%，第二个成分解释了 14%，依此类推。将这些值相加，我们得到以下结果：
In [ ]:
weighted_values.sum()
Out[ ]:
0.9216371041932268
前五个特征值将解释数据集方差的 92%。
使用 PCA 重建道琼斯指数
默认情况下，KernelPCA实例化时使用n_components=None参数，这将构建一个具有非零成分的核 PCA。我们还可以创建一个具有五个成分的 PCA 指数：
In [ ]:
import numpy as np
kernel_pca = KernelPCA(n_components=5).fit(df_z_components)
pca_5 = kernel_pca.transform(-daily_df_components)
weights = fn_weighted_avg(kernel_pca.lambdas_)
reconstructed_values = np.dot(pca_5, weights)
# Combine DJIA and PCA index for comparison
df_combined = djia_2017.copy()
df_combined['pca_5'] = reconstructed_values
df_combined = df_combined.apply(fn_z_score)
df_combined.plot(figsize=(12, 8));
使用fit()方法，我们使用具有五个成分的线性核 PCA 函数拟合了归一化数据集。transform()方法使用核 PCA 转换原始数据集。这些值使用由特征向量指示的权重进行归一化，通过点矩阵乘法计算。然后，我们使用copy()方法创建了道琼斯时间序列pandas DataFrame 的副本，并将其与df_combined DataFrame 中的重建值组合在一起。
新的 DataFrame 通过 z 得分进行归一化，并绘制出来，以查看重建的 PCA 指数跟踪原始道琼斯运动的情况。这给我们以下输出：
上面的图显示了 2017 年原始道琼斯指数与重建的道琼斯指数相比，使用了五个主成分。
平稳和非平稳时间序列
对于进行统计分析的时间序列数据，重要的是数据是平稳的，以便正确进行统计建模，因为这样的用途可能是用于预测和预测。本节介绍了时间序列数据中的平稳性和非平稳性的概念。
平稳性和非平稳性
在经验时间序列研究中，观察到价格变动向某些长期均值漂移，要么向上，要么向下。平稳时间序列是其统计特性（如均值、方差和自相关）随时间保持恒定的时间序列。相反，非平稳时间序列数据的观察结果其统计特性随时间变化，很可能是由于趋势、季节性、存在单位根或三者的组合。
在时间序列分析中，假设基础过程的数据是平稳的。否则，对非平稳数据进行建模可能会产生不可预测的结果。这将导致一种称为伪回归的情况。伪回归是指产生误导性的统计证据，表明独立的非平稳变量之间存在关系的回归。为了获得一致和可靠的结果，非平稳数据需要转换为平稳数据。
检查平稳性
有多种方法可以检查时间序列数据是平稳还是非平稳：
-
通过可视化：您可以查看时间序列图，以明显指示趋势或季节性。
-
通过统计摘要：您可以查看数据的统计摘要，寻找显著差异。例如，您可以对时间序列数据进行分组，并比较每组的均值和方差。
-
通过统计检验：您可以使用统计检验，如增广迪基-富勒检验，来检查是否满足或违反了平稳性期望。
非平稳过程的类型
以下几点有助于识别时间序列数据中的非平稳行为，以便考虑转换为平稳数据：
-
纯随机游走：具有单位根或随机趋势的过程。这是一个非均值回归的过程，其方差随时间演变并趋于无穷大。
-
带漂移的随机游走：具有随机游走和恒定漂移的过程。
-
确定性趋势：均值围绕着固定的趋势增长的过程，该趋势是恒定的且与时间无关。
-
带漂移和确定性趋势的随机游走：将随机游走与漂移分量和确定性趋势结合的过程。
平稳过程的类型
以下是时间序列研究中可能遇到的平稳性定义：
-
平稳过程：生成平稳观测序列的过程。
-
趋势平稳：不呈现趋势的过程。
-
季节性平稳：不呈现季节性的过程。
-
严格平稳：也称为强平稳。当随机变量的无条件联合概率分布在时间（或x轴上）移动时不发生变化的过程。
-
弱平稳：也称为协方差平稳或二阶平稳。当随机变量的均值、方差和相关性在时间移动时不发生变化的过程。
增广迪基-富勒检验
增广迪基-富勒检验（ADF）是一种统计检验，用于确定时间序列数据中是否存在单位根。单位根可能会导致时间序列分析中的不可预测结果。对单位根检验形成零假设，以确定时间序列数据受趋势影响的程度。通过接受零假设，我们接受时间序列数据是非平稳的证据。通过拒绝零假设，或接受备择假设，我们接受时间序列数据是由平稳过程生成的证据。这个过程也被称为趋势平稳。增广迪基-富勒检验统计量的值为负数。较低的 ADF 值表示更强烈地拒绝零假设。
以下是用于 ADF 测试的一些基本自回归模型：
- 没有常数和趋势：
- 没有常数和趋势：
- 带有常数和趋势：
这里，α是漂移常数，β是时间趋势的系数，γ是我们的假设系数，p是一阶差分自回归过程的滞后阶数，ϵ[t]是独立同分布的残差项。当α=0和β=0时，模型是一个随机游走过程。当β=0时，模型是一个带漂移的随机游走过程。滞后阶数p的选择应使得残差不具有序列相关性。一些选择滞后阶数的信息准则的方法包括最小化阿卡信息准则（AIC）、贝叶斯信息准则（BIC）和汉南-奎恩信息准则。
然后可以将假设表述如下：
-
零假设，H[0]：如果未能被拒绝，表明时间序列包含单位根并且是非平稳的
-
备择假设，H[1]：如果拒绝H[0]，则表明时间序列不包含单位根并且是平稳的
为了接受或拒绝零假设，我们使用 p 值。如果 p 值低于 5%甚至 1%的阈值，我们拒绝零假设。如果 p 值高于此阈值，我们可能未能拒绝零假设，并将时间序列视为非平稳的。换句话说，如果我们的阈值为 5%或 0.05，请注意以下内容：
-
p 值> 0.05：我们未能拒绝零假设H[0]，并得出结论，数据具有单位根并且是非平稳的
-
p 值≤0.05：我们拒绝零假设H[0]，并得出结论，数据具有单位根并且是非平稳的
statsmodels库提供了实现此测试的adfuller()函数。
分析具有趋势的时间序列
让我们检查一个时间序列数据集。例如，考虑在芝加哥商品交易所交易的黄金期货价格。在 Quandl 上，可以通过以下代码下载黄金期货连续合约：CHRIS/CME_GC1。这些数据由维基连续期货社区小组策划，仅考虑了最近的合约。数据集的第六列包含了结算价格。以下代码从 2000 年开始下载数据集：
In [ ]:
import quandl
QUANDL_API_KEY = 'BCzkk3NDWt7H9yjzx-DY' # Your Quandl key here
quandl.ApiConfig.api_key = QUANDL_API_KEY
df = quandl.get(
'CHRIS/CME_GC1',
column_index=6,
collapse='monthly',
start_date='2000-01-01')
使用以下命令检查数据集的头部：
In [ ]:
df.head()
我们得到以下表格：
| Settle | Date |
|---|---|
| 2000-01-31 | 283.2 |
| 2000-02-29 | 294.2 |
| 2000-03-31 | 278.4 |
| 2000-04-30 | 274.7 |
| 2000-05-31 | 271.7 |
将滚动均值和标准差计算到df_mean和df_std变量中，窗口期为一年：
In [ ] :
df_settle = df['Settle'].resample('MS').ffill().dropna()
df_rolling = df_settle.rolling(12)
df_mean = df_rolling.mean()
df_std = df_rolling.std()
resample()方法有助于确保数据在月度基础上平滑，并且ffill()方法向前填充任何缺失值。
可以在pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html找到用于指定resample()方法的常见有用时间序列频率列表.
让我们可视化滚动均值与原始时间序列的图表：
In [ ] :
plt.figure(figsize=(12, 8))
plt.plot(df_settle, label='Original')
plt.plot(df_mean, label='Mean')
plt.legend();
我们获得以下输出：
将滚动标准差可视化分开，我们得到以下结果：
In [ ] :
df_std.plot(figsize=(12, 8));
我们获得以下输出：
使用statsmodels模块，用adfuller()方法对我们的数据集进行 ADF 单位根检验：
In [ ]:
from statsmodels.tsa.stattools import adfuller
result = adfuller(df_settle)
print('ADF statistic: ', result[0])
print('p-value:', result[1])
critical_values = result[4]
for key, value in critical_values.items():
print('Critical value (%s): %.3f' % (key, value))
Out[ ]:
ADF statistic: -1.4017828015895548
p-value: 0.5814211232134314
Critical value (1%): -3.461
Critical value (5%): -2.875
Critical value (10%): -2.574
adfuller()方法返回一个包含七个值的元组。特别地，我们对第一个、第二个和第五个值感兴趣，它们分别给出了检验统计量、p 值和临界值字典。
从图表中可以观察到，均值和标准差随时间波动，均值呈现总体上升趋势。ADF 检验统计值大于临界值（特别是在 5%时），p-value大于 0.05。基于这些结果，我们无法拒绝存在单位根的原假设，并认为我们的数据是非平稳的。
使时间序列平稳
非平稳时间序列数据可能受到趋势或季节性的影响。趋势性时间序列数据的均值随时间不断变化。受季节性影响的数据在特定时间间隔内有变化。在使时间序列数据平稳时，必须去除趋势和季节性影响。去趋势、差分和分解就是这样的方法。然后得到的平稳数据适合进行统计预测。
让我们详细看看所有三种方法。
去趋势
从非平稳数据中去除趋势线的过程称为去趋势。这涉及一个将大值归一化为小值的转换步骤。例如可以是对数函数、平方根函数，甚至是立方根。进一步的步骤是从移动平均值中减去转换值。
让我们对相同的数据集df_settle执行去趋势，使用对数变换并从两个周期的移动平均值中减去，如下 Python 代码所示：
In [ ]:
import numpy as np
df_log = np.log(df_settle)
In [ ]:
df_log_ma= df_log.rolling(2).mean()
df_detrend = df_log - df_log_ma
df_detrend.dropna(inplace=True)
# Mean and standard deviation of detrended data
df_detrend_rolling = df_detrend.rolling(12)
df_detrend_ma = df_detrend_rolling.mean()
df_detrend_std = df_detrend_rolling.std()
# Plot
plt.figure(figsize=(12, 8))
plt.plot(df_detrend, label='Detrended')
plt.plot(df_detrend_ma, label='Mean')
plt.plot(df_detrend_std, label='Std')
plt.legend(loc='upper right');
df_log变量是我们使用numpy模块的对数函数转换的pandas DataFrame，df_detrend变量包含去趋势数据。我们绘制这些去趋势数据，以可视化其在滚动一年期间的均值和标准差。
我们得到以下输出：
观察到均值和标准差没有表现出长期趋势。
观察去趋势数据的 ADF 检验统计量，我们得到以下结果：
In [ ]:
from statsmodels.tsa.stattools import adfuller
result = adfuller(df_detrend)
print('ADF statistic: ', result[0])
print('p-value: %.5f' % result[1])
critical_values = result[4]
for key, value in critical_values.items():
print('Critical value (%s): %.3f' % (key, value))
Out[ ]:
ADF statistic: -17.04239232215001
p-value: 0.00000
Critical value (1%): -3.460
Critical value (5%): -2.874
Critical value (10%): -2.574
这个去趋势数据的p-value小于 0.05。我们的 ADF 检验统计量低于所有临界值。我们可以拒绝原假设，并说这个数据是平稳的。
通过差分去除趋势
差分涉及将时间序列值与时间滞后进行差分。时间序列的一阶差分由以下公式给出：
我们可以重复使用前一节中的df_log变量作为我们的对数转换时间序列，并利用 NumPy 模块的diff()和shift()方法进行差分，代码如下：
In [ ]:
df_log_diff = df_log.diff(periods=3).dropna()
# Mean and standard deviation of differenced data
df_diff_rolling = df_log_diff.rolling(12)
df_diff_ma = df_diff_rolling.mean()
df_diff_std = df_diff_rolling.std()
# Plot the stationary data
plt.figure(figsize=(12, 8))
plt.plot(df_log_diff, label='Differenced')
plt.plot(df_diff_ma, label='Mean')
plt.plot(df_diff_std, label='Std')
plt.legend(loc='upper right');
diff()的参数periods=3表示在计算差异时数据集向后移动三个周期。
这提供了以下输出：
从图表中可以观察到，滚动均值和标准差随时间变化很少。
观察我们的 ADF 检验统计量，我们得到以下结果：
In [ ]:
from statsmodels.tsa.stattools import adfuller
result = adfuller(df_log_diff)
print('ADF statistic:', result[0])
print('p-value: %.5f' % result[1])
critical_values = result[4]
for key, value in critical_values.items():
print('Critical value (%s): %.3f' % (key, value))
Out[ ]:
ADF statistic: -2.931684356800213
p-value: 0.04179
Critical value (1%): -3.462
Critical value (5%): -2.875
Critical value (10%): -2.574
从 ADF 检验中，此数据的p-value小于 0.05。我们的 ADF 检验统计量低于 5%的临界值，表明此数据在 95%的置信水平下是平稳的。我们可以拒绝原假设，并说这个数据是平稳的。
季节性分解
分解涉及对趋势和季节性进行建模，然后将它们移除。我们可以使用statsmodel.tsa.seasonal模块来使用移动平均模型非平稳时间序列数据，并移除其趋势和季节性成分。
通过重复使用包含先前部分数据集对数的df_log变量，我们得到以下结果：
In [ ]:
from statsmodels.tsa.seasonal import seasonal_decompose
decompose_result = seasonal_decompose(df_log.dropna(), freq=12)
df_trend = decompose_result.trend
df_season = decompose_result.seasonal
df_residual = decompose_result.resid
statsmodels.tsa.seasonal的seasonal_decompose()方法需要一个参数freq，它是一个整数值，指定每个季节周期的周期数。由于我们使用的是月度数据，我们期望每个季节年有 12 个周期。该方法返回一个对象，主要包括趋势和季节分量，以及最终的pandas系列数据，其趋势和季节分量已被移除。
有关statsmodels.tsa.seasonal模块的seasonal_decompose()方法的更多信息可以在www.statsmodels.org/dev/generated/statsmodels.tsa.seasonal.seasonal_decompose.html找到。
通过运行以下 Python 代码来可视化不同的图表：
In [ ]:
plt.rcParams['figure.figsize'] = (12, 8)
fig = decompose_result.plot()
我们得到以下图表：
在这里，我们可以看到单独的趋势和季节分量从数据集中被移除并绘制，残差在底部绘制。让我们可视化一下我们的残差的统计特性：
In [ ]:
df_log_diff = df_residual.diff().dropna()
# Mean and standard deviation of differenced data
df_diff_rolling = df_log_diff.rolling(12)
df_diff_ma = df_diff_rolling.mean()
df_diff_std = df_diff_rolling.std()
# Plot the stationary data
plt.figure(figsize=(12, 8))
plt.plot(df_log_diff, label='Differenced')
plt.plot(df_diff_ma, label='Mean')
plt.plot(df_diff_std, label='Std')
plt.legend();
我们得到以下图表：
从图表中观察到，滚动均值和标准差随时间变化很少。
通过检查我们的残差数据的平稳性，我们得到以下结果：
In [ ]:
from statsmodels.tsa.stattools import adfuller
result = adfuller(df_residual.dropna())
print('ADF statistic:', result[0])
print('p-value: %.5f' % result[1])
critical_values = result[4]
for key, value in critical_values.items():
print('Critical value (%s): %.3f' % (key, value))
Out[ ]:
ADF statistic: -6.468683205304995
p-value: 0.00000
Critical value (1%): -3.463
Critical value (5%): -2.876
Critical value (10%): -2.574
从 ADF 测试中，这些数据的p-value小于 0.05。我们的 ADF 测试统计量低于所有临界值。我们可以拒绝零假设，并说这些数据是平稳的。
ADF 测试的缺点
在使用 ADF 测试可靠检查非平稳数据时，有一些考虑事项：
-
ADF 测试不能真正区分纯单元根生成过程和非单元根生成过程。在长期移动平均过程中，ADF 测试在拒绝零假设方面存在偏差。其他检验平稳性的方法，如Kwiatkowski–Phillips–Schmidt–Shin（KPSS）检验和Phillips-Perron检验，采用了不同的方法来处理单位根的存在。
-
在确定滞后长度p时没有固定的方法。如果p太小，剩余误差中的序列相关性可能会影响检验的大小。如果p太大，检验的能力将会下降。对于这个滞后阶数，还需要进行额外的考虑。
-
随着确定性项被添加到测试回归中，单位根测试的能力减弱。
时间序列的预测和预测
在上一节中，我们确定了时间序列数据中的非平稳性，并讨论了使时间序列数据平稳的技术。有了平稳的数据，我们可以进行统计建模，如预测和预测。预测涉及生成样本内数据的最佳估计。预测涉及生成样本外数据的最佳估计。预测未来值是基于先前观察到的值。其中一个常用的方法是自回归积分移动平均法。
关于自回归积分移动平均
自回归积分移动平均（ARIMA）是基于线性回归的平稳时间序列的预测模型。顾名思义，它基于三个组件：
-
自回归（AR）：使用观察和滞后值之间的依赖关系的模型
-
积分（I）：使用差分观察和以前时间戳的观察来使时间序列平稳
-
移动平均（MA）：使用观察误差项和先前误差项的组合之间的依赖关系的模型，e[t]
ARIMA 模型的标记是ARIMA(p, d, q)，对应于三个组件的参数。可以通过改变p、d和q的值来指定非季节性 ARIMA 模型，如下所示：
- ARIMA(p,0,0): 一阶自回归模型，用AR(p)表示。p是滞后阶数，表示模型中滞后观察值的数量。例如，ARIMA(2,0,0)是AR(2)，表示如下：
在这里，ϕ[1]和ϕ[2]是模型的参数。
- ARIMA(0,d,0): 整合分量中的一阶差分，也称为随机游走，用I(d)表示。d是差分的程度，表示数据被减去过去值的次数。例如，ARIMA(0,1,0)是I(1)，表示如下：
在这里，μ是季节差分的均值。
- ARIMA(0,0,q)：移动平均分量，用MA(q)表示。阶数q决定了模型中要包括的项数：
通过网格搜索找到模型参数
网格搜索，也称为超参数优化方法，可用于迭代地探索不同的参数组合，以拟合我们的 ARIMA 模型。我们可以在每次迭代中使用statsmodels模块的SARIMAX()函数拟合季节性 ARIMA 模型，返回一个MLEResults类的对象。MLEResults对象具有一个aic属性，用于返回 AIC 值。具有最低 AIC 值的模型为我们提供了最佳拟合模型，确定了我们的p、d和q参数。有关 SARIMAX 的更多信息，请访问www.statsmodels.org/dev/generated/statsmodels.tsa.statespace.sarimax.SARIMAX.html。
我们将网格搜索过程定义为arima_grid_search()函数，如下所示：
In [ ]:
import itertools
import warnings
from statsmodels.tsa.statespace.sarimax import SARIMAX
warnings.filterwarnings("ignore")
def arima_grid_search(dataframe, s):
p = d = q = range(2)
param_combinations = list(itertools.product(p, d, q))
lowest_aic, pdq, pdqs = None, None, None
total_iterations = 0
for order in param_combinations:
for (p, q, d) in param_combinations:
seasonal_order = (p, q, d, s)
total_iterations += 1
try:
model = SARIMAX(df_settle, order=order,
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False,
disp=False
)
model_result = model.fit(maxiter=200, disp=False)
if not lowest_aic or model_result.aic < lowest_aic:
lowest_aic = model_result.aic
pdq, pdqs = order, seasonal_order
except Exception as ex:
continue
return lowest_aic, pdq, pdqs
我们的变量df_settle保存了我们在上一节中下载的期货数据的月度价格。在SARIMAX（具有外生回归器的季节性自回归整合移动平均模型）函数中，我们提供了seasonal_order参数，这是ARIMA(p,d,q,s)季节性分量，其中s是数据集中一个季节的周期数。由于我们使用的是月度数据，我们使用 12 个周期来定义季节模式。enforce_stationarity=False参数不会将 AR 参数转换为强制模型的 AR 分量。enforce_invertibility=False参数不会将 MA 参数转换为强制模型的 MA 分量。disp=False参数在拟合模型时抑制输出信息。
定义了网格函数后，我们现在可以使用我们的月度数据调用它，并打印出具有最低 AIC 值的模型参数：
In [ ]:
lowest_aic, order, seasonal_order = arima_grid_search(df_settle, 12)
In [ ]:
print('ARIMA{}x{}'.format(order, seasonal_order))
print('Lowest AIC: %.3f'%lowest_aic)
Out[ ]:
ARIMA(0, 1, 1)x(0, 1, 1, 12)
Lowest AIC: 2149.636
ARIMA(0,1,1,12)季节性分量模型将在 2149.636 的 AIC 值处得到最低值。我们将使用这些参数在下一节中拟合我们的 SARIMAX 模型。
拟合 SARIMAX 模型
获得最佳模型参数后，使用summary()方法检查拟合结果的模型属性，以查看详细的统计信息：
In [ ]:
model = SARIMAX(
df_settle,
order=order,
seasonal_order=seasonal_order,
enforce_stationarity=False,
enforce_invertibility=False,
disp=False
)
model_results = model.fit(maxiter=200, disp=False)
print(model_results.summary())
这给出了以下输出：
Statespace Model Results
==========================================================================================
Dep. Variable: Settle No. Observations: 226
Model: SARIMAX(0, 1, 1)x(0, 1, 1, 12) Log Likelihood -1087.247
Date: Sun, 02 Dec 2018 AIC 2180.495
Time: 17:38:32 BIC 2190.375
Sample: 02-01-2000 HQIC 2184.494
- 11-01-2018
Covariance Type: opg
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
ma.L1 -0.1716 0.044 -3.872 0.000 -0.258 -0.085
ma.S.L12 -1.0000 447.710 -0.002 0.998 -878.496 876.496
sigma2 2854.6342 1.28e+06 0.002 0.998 -2.5e+06 2.51e+06
===================================================================================
Ljung-Box (Q): 67.93 Jarque-Bera (JB): 52.74
Prob(Q): 0.00 Prob(JB): 0.00
Heteroskedasticity (H): 6.98 Skew: -0.34
Prob(H) (two-sided): 0.00 Kurtosis: 5.43
===================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).
重要的是要运行模型诊断，以调查模型假设是否被违反：
In [ ]:
model_results.plot_diagnostics(figsize=(12, 8));
我们得到以下输出：
右上角的图显示了标准化残差的核密度估计（KDE），这表明误差服从均值接近于零的高斯分布。让我们看一下残差的更准确的统计信息：
In [ ] :
model_results.resid.describe()
Out[ ]:
count 223.000000
mean 0.353088
std 57.734027
min -196.799109
25% -22.036234
50% 3.500942
75% 22.872743
max 283.200000
dtype: float64
从残差的描述中，非零均值表明预测可能存在正向偏差。
预测和预测 SARIMAX 模型
model_results变量是statsmodel模块的SARIMAXResults对象，代表 SARIMAX 模型的输出。它包含一个get_prediction()方法，用于执行样本内预测和样本外预测。它还包含一个conf_int()方法，返回预测的置信区间，包括拟合参数的下限和上限，默认情况下为 95%置信区间。让我们应用这些方法：
In [ ]:
n = len(df_settle.index)
prediction = model_results.get_prediction(
start=n-12*5,
end=n+5
)
prediction_ci = prediction.conf_int()
get_prediction()方法中的start参数表示我们正在对最近五年的价格进行样本内预测。同时，使用end参数，我们正在对接下来的五个月进行样本外预测。
通过检查前三个预测的置信区间值，我们得到以下结果：
In [ ]:
print(prediction_ci.head(3))
Out[ ]:
lower Settle upper Settle
2017-09-01 1180.143917 1396.583325
2017-10-01 1204.307842 1420.747250
2017-11-01 1176.828881 1393.268289
让我们将预测和预测的价格与我们的原始数据集从 2008 年开始进行对比：
In [ ]:
plt.figure(figsize=(12, 6))
ax = df_settle['2008':].plot(label='actual')
prediction_ci.plot(
ax=ax, style=['--', '--'],
label='predicted/forecasted')
ci_index = prediction_ci.index
lower_ci = prediction_ci.iloc[:, 0]
upper_ci = prediction_ci.iloc[:, 1]
ax.fill_between(ci_index, lower_ci, upper_ci,
color='r', alpha=.1)
ax.set_xlabel('Time (years)')
ax.set_ylabel('Prices')
plt.legend()
plt.show()
这给我们以下输出：
实线图显示了观察值，而虚线图显示了五年滚动预测，紧密跟随并受到阴影区域内的置信区间的限制。请注意，随着接下来五个月的预测进入未来，置信区间扩大以反映对前景的不确定性。
摘要
在本章中，我们介绍了 PCA 作为投资组合建模中的降维技术。通过将投资组合资产价格的波动分解为其主要成分或共同因素，可以保留最有用的因素，并且可以大大简化投资组合分析，而不会影响计算时间和空间复杂性。通过使用sklearn.decomposition模块的KernelPCA函数将 PCA 应用于道琼指数及其 30 个成分，我们获得了特征向量和特征值，用于用五个成分重构道琼指数。
在时间序列数据的统计分析中，数据被视为平稳或非平稳。平稳的时间序列数据是其统计特性随时间保持不变的数据。非平稳的时间序列数据其统计特性随时间变化，很可能是由于趋势、季节性、存在单位根或三者的组合。从非平稳数据建模可能产生虚假回归。为了获得一致和可靠的结果，非平稳数据需要转换为平稳数据。
我们使用统计检验，如 ADF，来检查是否满足或违反了平稳性期望。statsmodels.tsa.stattools模块的adfuller方法提供了检验统计量、p 值和临界值，从中我们可以拒绝零假设，即数据具有单位根且是非平稳的。
我们通过去趋势化、差分和季节性分解将非平稳数据转换为平稳数据。通过使用 ARIMA，我们使用statsmodels.tsa.statespace.sarimax模块的SARIMAX函数拟合模型，以找到通过迭代网格搜索过程给出最低 AIC 值的合适模型参数。拟合结果用于预测和预测。
在下一章中，我们将使用 VIX 进行交互式金融分析。
第三部分：实践方法
在本部分中，我们将应用第 1 部分Python 入门和第 2 部分金融概念中涵盖的理论概念，构建完全功能的工作系统。
本部分将包括以下章节：
-
[第七章]，使用 VIX 进行交互式金融分析
-
[第八章]，构建算法交易平台
-
[第九章]，实施回测系统
-
[第十章]，金融的机器学习
-
[第十一章]，金融的深度学习
第七章：与 VIX 一起进行交互式金融分析
投资者使用波动率衍生品来分散和对冲他们在股票和信用组合中的风险。由于股票基金的长期投资者面临下行风险，波动率可以用作尾部风险的对冲和看涨期权的替代品。在美国，芝加哥期权交易所（CBOE）的波动率指数（VIX），或简称 VIX，衡量了具有平均到期日为 30 天的标准普尔 500 股票指数期权隐含的短期波动率。世界各地许多人使用 VIX 来衡量未来 30 天的股票市场波动性。在欧洲，等效的波动率对应指标是 EURO STOXX 50 波动率（VSTOXX）市场指数。对于利用 S&P 500 指数的基准策略，其与 VIX 的负相关性的性质提供了一种可行的方式来避免基准再平衡成本。波动性的统计性质允许交易员执行均值回归策略、离散交易和波动性价差交易等策略。
在本章中，我们将看看如何对 VIX 和 S&P 500 指数进行数据分析。使用 S&P 500 指数期权，我们可以重建 VIX 并将其与观察值进行比较。这里呈现的代码在 Jupyter Notebook 上运行，这是 Python 的交互式组件，可以帮助我们可视化数据并研究它们之间的关系。
在本章中，我们将讨论以下主题：
-
介绍 EURO STOXX 50 指数、VSTOXX 和 VIX
-
对 S&P 500 指数和 VIX 进行金融分析
-
根据 CBOE VIX 白皮书逐步重建 VIX 指数
-
寻找 VIX 指数的近期和次期期权
-
确定期权数据集的行权价边界
-
通过行权价对 VIX 的贡献进行制表
-
计算近期和次期期权的远期水平
-
计算近期和次期期权的波动率值
-
同时计算多个 VIX 指数
-
将计算出的指数结果与实际的标准普尔 500 指数进行比较
波动率衍生品
全球最受欢迎的两个波动率指数是 VIX 和 VSTOXX，分别在美国和欧洲可用。VIX 基于标准普尔 500 指数，在 CBOE 上发布。虽然 VIX 本身不进行交易，但 VIX 的衍生产品，如期权、期货、交易所交易基金和一系列基于波动性的证券可供投资者选择。CBOE 网站提供了许多期权和市场指数的全面信息，如标准普尔 500 标准和周期期权，以及我们可以分析的 VIX。在本章的后面部分，我们将首先了解这些产品的背景，然后进行金融分析。
STOXX 和 Eurex
在美国，标准普尔 500 指数是最广泛关注的股票市场指数之一，由标准普尔道琼斯指数创建。在欧洲，STOXX 有限公司是这样一家公司。
成立于 1997 年，STOXX 有限公司总部位于瑞士苏黎世，在全球计算大约 7000 个指数。作为一个指数提供商，它开发、维护、分发和推广一系列严格基于规则和透明的指数。
STOXX 在这些类别提供了许多股票指数：基准指数、蓝筹股指数、股息指数、规模指数、行业指数、风格指数、优化指数、策略指数、主题指数、可持续性指数、信仰指数、智能贝塔指数和计算产品。
Eurex 交易所是德国法兰克福的衍生品交易所，提供超过 1900 种产品，包括股票指数、期货、期权、交易所交易基金、股息、债券和回购。STOXX 的许多产品和衍生品在 Eurex 上交易。
EURO STOXX 50 指数
由 STOXX 有限公司设计，EURO STOXX 50 指数是全球最流动的股票指数之一，服务于 Eurex 上列出的许多指数产品。它于 1998 年 2 月 26 日推出，由来自奥地利、比利时、芬兰、法国、德国、希腊、爱尔兰、意大利、卢森堡、荷兰、葡萄牙和西班牙的 50 只蓝筹股组成。EURO STOXX 50 指数期货和期权合约可在 Eurex 交易所上交易。指数每 15 秒基于实时价格重新计算一次。
EURO STOXX 50 指数的股票代码是 SX5E。EURO STOXX 50 指数期权的股票代码是 OESX。
VSTOXX
VSTOXX 或 EURO STOXX 50 波动率是由 Eurex 交易所提供服务的一类波动率衍生品。VSTOXX 市场指数基于一篮子 OESX 报价的平价或虚价。它衡量了未来 30 天在 EURO STOXX 50 指数上的隐含市场波动率。
投资者利用波动率衍生品进行基准策略，利用 EURO STOXX 50 指数的负相关性，可以避免基准再平衡成本。波动率的统计性质使交易员能够执行均值回归策略、离散交易和波动率价差交易等。指数每 5 秒重新计算一次。
VSTOXX 的股票代码是 V2TX。基于 VSTOXX 指数的 VSTOXX 期权和 VSTOXX 迷你期货在 Eurex 交易所交易。
标普 500 指数
标普 500 指数（SPX）的历史可以追溯到 1923 年，当时它被称为综合指数。最初，它跟踪了少量股票。1957 年，跟踪的股票数量扩大到 500 只，并成为了 SPX。
构成 SPX 的股票在纽约证券交易所（NYSE）或全国证券经纪人自动报价系统（NASDAQ）上公开上市。该指数被认为是美国经济的主要代表，通过大市值普通股。指数每 15 秒重新计算一次，并由路透美国控股公司分发。
交易所使用的常见股票代码是 SPX 和 INX，一些网站上是^GSPC。
SPX 期权
芝加哥期权交易所提供各种期权合约进行交易，包括标普 500 指数等股票指数期权。SPX 指数期权产品具有不同的到期日。标准或传统的 SPX 期权每月第三个星期五到期，并在交易日开始结算。SPX 周期（SPXW）期权产品可能每周到期，分别在星期一、星期三和星期五，或每月在月末最后一个交易日到期。如果到期日落在交易所假日，到期日将提前至前一个交易日。其他 SPX 期权是迷你期权，交易量为名义规模的十分之一，以及标普 500 指数存托凭证交易基金（SPDR ETF）。大多数 SPX 指数期权是欧式风格，除了 SPDR ETF 是美式风格。
VIX
与 STOXX 一样，芝加哥期权交易所 VIX 衡量了标普 500 股票指数期权价格隐含的短期波动率。芝加哥期权交易所 VIX 始于 1993 年，基于标普 100 指数，于 2003 年更新为基于 SPX，并于 2014 年再次更新以包括 SPXW 期权。全球许多人认为 VIX 是未来 30 天股市波动的流行测量工具。VIX 每 15 秒重新计算一次，并由芝加哥期权交易所分发。
VIX 期权和 VIX 期货基于 VIX，在芝加哥期权交易所交易。
标普 500 和 VIX 的金融分析
在本节中，我们将研究 VIX 与标普 500 市场指数之间的关系。
收集数据
我们将使用 Alpha Vantage 作为我们的数据提供商。让我们按以下步骤下载 SPX 和 VIX 数据集：
- 查询具有股票代码
^GSPC的全时 S&P 500 历史数据：
In [ ]:
from alpha_vantage.timeseries import TimeSeries
# Update your Alpha Vantage API key here...
ALPHA_VANTAGE_API_KEY = 'PZ2ISG9CYY379KLI'
ts = TimeSeries(key=ALPHA_VANTAGE_API_KEY, output_format='pandas')
df_spx_data, meta_data = ts.get_daily_adjusted(
symbol='^GSPC', outputsize='full')
- 对于具有股票代码
^VIX的 VIX 指数也做同样的操作：
In [ ]:
df_vix_data, meta_data = ts.get_daily_adjusted(
symbol='^VIX', outputsize='full')
- 检查 DataFrame 对象
df_spx_data的内容：
In [ ]:
df_spx_data.info()
Out[ ]:
<class 'pandas.core.frame.DataFrame'>
Index: 4774 entries, 2000-01-03 to 2018-12-21
Data columns (total 8 columns):
1\. open 4774 non-null float64
2\. high 4774 non-null float64
3\. low 4774 non-null float64
4\. close 4774 non-null float64
5\. adjusted close 4774 non-null float64
6\. volume 4774 non-null float64
7\. dividend amount 4774 non-null float64
8\. split coefficient 4774 non-null float64
dtypes: float64(8)
memory usage: 317.0+ KB
- 检查 DataFrame 对象
df_vix_data的内容：
In [ ]:
df_vix_data.info()
Out[ ]:
<class 'pandas.core.frame.DataFrame'>
Index: 4774 entries, 2000-01-03 to 2018-12-21
Data columns (total 8 columns):
1\. open 4774 non-null float64
2\. high 4774 non-null float64
3\. low 4774 non-null float64
4\. close 4774 non-null float64
5\. adjusted close 4774 non-null float64
6\. volume 4774 non-null float64
7\. dividend amount 4774 non-null float64
8\. split coefficient 4774 non-null float64
dtypes: float64(8)
memory usage: 317.0+ KB
- 注意，两个数据集的开始日期都是从 2000 年 1 月 3 日开始的，第五列标记为
5\. adjusted close包含我们感兴趣的值。提取这两列并将它们合并成一个pandasDataFrame：
In [ ]:
import pandas as pd
df = pd.DataFrame({
'SPX': df_spx_data['5\. adjusted close'],
'VIX': df_vix_data['5\. adjusted close']
})
df.index = pd.to_datetime(df.index)
pandas的to_datetime()方法将作为字符串对象给出的交易日期转换为 pandas 的DatetimeIndex对象。检查我们最终的 DataFrame 对象df的头部，得到以下结果：
In [ ]:
df.head(3)
这给我们以下表格：
| 日期 | SPX | VIX |
|---|---|---|
| 2000-01-03 | 1455.22 | 24.21 |
| 2000-01-04 | 1399.42 | 27.01 |
| 2000-01-05 | 1402.11 | 26.41 |
查看我们格式化的指数，得到以下结果：
In [ ]:
df.index
Out[ ]:
DatetimeIndex(['2000-01-03', '2000-01-04', '2000-01-05', '2000-01-06',
'2000-01-07', '2000-01-10', '2000-01-11', '2000-01-12',
'2000-01-13', '2000-01-14',
...
'2018-10-11', '2018-10-12', '2018-10-15', '2018-10-16',
'2018-10-17', '2018-10-18', '2018-10-19', '2018-10-22',
'2018-10-23', '2018-10-24'],
dtype='datetime64[ns]', name='date', length=4734, freq=None)
有了正确格式的pandas DataFrame，让我们继续处理这个数据集。
执行分析
pandas的describe()方法给出了 DataFrame 对象中每列的摘要统计和值的分布：
In [ ]:
df.describe()
这给我们以下表格：
| SPX | |
|---|---|
| count | 4734.000000 |
| mean | 1493.538998 |
| std | 500.541938 |
| min | 676.530000 |
| 25% | 1140.650000 |
| 50% | 1332.730000 |
| 75% | 1840.515000 |
| max | 2930.750000 |
另一个相关的方法info()，之前使用过，给我们 DataFrame 的技术摘要，如指数范围和内存使用情况：
In [ ]:
df.info()
Out[ ]:
<class 'pandas.core.frame.DataFrame'>
DatetimeIndex: 4734 entries, 2000-01-03 to 2018-10-24
Data columns (total 2 columns):
SPX 4734 non-null float64
VIX 4734 non-null float64
dtypes: float64(2)
memory usage: 111.0 KB
让我们绘制 S&P 500 和 VIX，看看它们从 2010 年开始是什么样子：
In [ ]:
%matplotlib inline
import matplotlib.pyplot as plt
plt.figure(figsize = (12, 8))
ax_spx = df['SPX'].plot()
ax_vix = df['VIX'].plot(secondary_y=True)
ax_spx.legend(loc=1)
ax_vix.legend(loc=2)
plt.show();
这给我们以下图表：
注意，当 S&P 500 上涨时，VIX 似乎下降，表现出负相关关系。我们需要进行更多的统计分析来确保。
也许我们对两个指数的日回报感兴趣。diff()方法返回前一期值之间的差异集。直方图可用于在 100 个 bin 间隔上给出数据密度估计的大致感觉：
In [ ]:
df.diff().hist(
figsize=(10, 5),
color='blue',
bins=100);
hist()方法给出了以下直方图：
使用pct_change()命令也可以实现相同的效果，它给出了前一期值的百分比变化：
In [ ]:
df.pct_change().hist(
figsize=(10, 5),
color='blue',
bins=100);
我们得到了相同的直方图，以百分比变化为单位：
对于收益的定量分析，我们对每日收益的对数感兴趣。为什么使用对数收益而不是简单收益？有几个原因，但最重要的是归一化，这避免了负价格的问题。
我们可以使用pandas的shift()函数将值向前移动一定数量的期间。dropna()方法删除对数计算转换末尾未使用的值。NumPy 的log()方法帮助计算 DataFrame 对象中所有值的对数，并将其存储在log_returns变量中作为 DataFrame 对象。然后可以绘制对数值，得到每日对数收益的图表。以下是绘制对数值的代码：
In [ ]:
import numpy as np
log_returns = np.log(df / df.shift(1)).dropna()
log_returns.plot(
subplots=True,
figsize=(10, 8),
color='blue',
grid=True
);
for ax in plt.gcf().axes:
ax.legend(loc='upper left')
我们得到以下输出：
顶部和底部图表分别显示了 SPX 和 VIX 的对数收益，从 2000 年到现在的期间。
SPX 和 VIX 之间的相关性
我们可以使用corr()方法来推导pandas DataFrame 对象中每列值之间的相关值，如以下 Python 示例所示：
In [ ]:
log_returns.corr()
这给我们以下相关性表：
| SPX | VIX | |
|---|---|---|
| SPX | 1.000000 | -0.733161 |
| VIX | -0.733161 | 1.000000 |
在-0.731433 处，SPX 与 VIX 呈负相关。为了更好地可视化这种关系，我们可以将每组日对数收益值绘制成散点图。使用statsmodels.api模块来获得散点数据之间的普通最小二乘回归线：
In [ ]:
import statsmodels.api as sm
log_returns.plot(
figsize=(10,8),
x="SPX",
y="VIX",
kind='scatter')
ols_fit = sm.OLS(log_returns['VIX'].values,
log_returns['SPX'].values).fit()
plt.plot(log_returns['SPX'], ols_fit.fittedvalues, 'r');
我们得到以下输出：
如前图所示的向下倾斜回归线证实了标普 500 和 VIX 指数之间的负相关关系。
pandas的rolling().corr()方法计算两个时间序列之间的滚动窗口相关性。我们使用252的值来表示移动窗口中的交易日数，以计算年度滚动相关性，使用以下命令：
In [ ]:
plt.ylabel('Rolling Annual Correlation')
df_corr = df['SPX'].rolling(252).corr(other=df['VIX'])
df_corr.plot(figsize=(12,8));
我们得到以下输出：
从前面的图表可以看出，SPX 和 VIX 呈负相关，在大部分时间内在 0.0 和-0.9 之间波动，每年使用 252 个交易日。
计算 VIX 指数
在本节中，我们将逐步复制 VIX 指数。VIX 指数的计算在芝加哥期权交易所网站上有记录。您可以在www.cboe.com/micro/vix/vixwhite.pdf获取芝加哥期权交易所 VIX 白皮书的副本。
导入 SPX 期权数据
假设您从经纪人那里收集了 SPX 期权数据，或者从 CBOE 网站等外部来源购买了历史数据。为了本章的目的，我们观察了 2018 年 10 月 15 日（星期一）到 2018 年 10 月 19 日（星期五）的期末 SPX 期权链价格，并将其保存为逗号分隔值（CSV）文件。这些文件的示例副本提供在源代码存储库的文件夹下。
在下面的示例中，编写一个名为read_file()的函数，它接受文件路径作为第一个参数，指示 CSV 文件的位置，并返回元数据列表和期权链数据列表的元组：
In [ ]:
import csv
META_DATA_ROWS = 3 # Header data starts at line 4
COLS = 7 # Each option data occupy 7 columns
def read_file(filepath):
meta_rows = []
calls_and_puts = []
with open(filepath, 'r') as file:
reader = csv.reader(file)
for row, cells in enumerate(reader):
if row < META_DATA_ROWS:
meta_rows.append(cells)
else:
call = cells[:COLS]
put = cells[COLS:-1]
calls_and_puts.append((call, put))
return (meta_rows, calls_and_puts)
请注意，您自己的期权数据结构可能与此示例不同。请谨慎检查并相应修改此函数。导入数据集后，我们可以继续解析和提取有用信息。
解析 SPX 期权数据
在这个例子中，我们假设 CSV 文件的前三行包含元信息，后面的期权链价格从第四行开始。对于每一行期权定价数据，前七列包含看涨合同的买入和卖出报价，接下来的七列是看跌合同的。每七列的第一列包含描述到期日、行权价和合同代码的字符串。按照以下步骤从我们的 CSV 文件中解析信息：
- 每行元信息都被添加到名为
meta_data的列表变量中，而每行期权数据都被添加到名为calls_and_puts的列表变量中。使用此函数读取单个文件会得到以下结果：
In [ ]:
(meta_rows, calls_and_puts) = \
read_file('files/chapter07/SPX_EOD_2018_10_15.csv')
- 打印每行元数据提供以下内容：
In [ ]:
for line in meta_rows:
print(line)
Out[ ]:
['SPX (S&P 500 INDEX)', '2750.79', '-16.34']
['Oct 15 2018 @ 20:00 ET']
['Calls', 'Last Sale', 'Net', 'Bid', 'Ask', 'Vol', 'Open Int', 'Puts', 'Last Sale', 'Net', 'Bid', 'Ask', 'Vol', 'Open Int']
- 期权报价的当前时间可以在我们的元数据的第二行找到。由于东部时间比格林威治标准时间（GMT）晚 5 小时，我们将替换
ET字符串并将整个字符串解析为datetime对象。以下函数get_dt_current()演示了这一点：
In [ ]:
from dateutil import parser
def get_dt_current(meta_rows):
"""
Extracts time information.
:param meta_rows: 2D array
:return: parsed datetime object
"""
# First cell of second row contains time info
date_time_row = meta_rows[1][0]
# Format text as ET time string
current_time = date_time_row.strip()\
.replace('@ ', '')\
.replace('ET', '-05:00')\
.replace(',', '')
dt_current = parser.parse(current_time)
return dt_current
- 从我们的期权数据的元信息中提取日期和时间信息作为芝加哥当地时间：
In [ ]:
dt_current = get_dt_current(meta_rows)
print(dt_current)
Out[ ]:
2018-10-15 20:00:00-05:00
- 现在，让我们看一下我们的期权报价数据的前两行：
In [ ]:
for line in calls_and_puts[:2]:
print(line)
Out[ ]:
(['2018 Oct 15 1700.00 (SPXW1815J1700)', '0.0', '0.0', '1039.30', '1063.00', '0', '0'], ['2018 Oct 15 1700.00 (SPXW1815V1700)', '0.15', '0.0', ' ', '0.05', '0'])
(['2018 Oct 15 1800.00 (SPXW1815J1800)', '0.0', '0.0', '939.40', '963.00', '0', '0'], ['2018 Oct 15 1800.00 (SPXW1815V1800)', '0.10', '0.0', ' ', '0.05', '0'])
列表中的每个项目都包含两个对象的元组，每个对象都包含一个看涨期权和一个看跌期权定价数据的列表，这些数据具有相同的行权价。参考我们打印的标题，每个期权价格列表数据的七个项目包含合同代码和到期日、最后成交价、价格净变动、买价、卖价、成交量和未平仓量。
让我们编写一个函数来解析每个 SPX 期权数据集的描述：
In [ ]:
from decimal import Decimal
def parse_expiry_and_strike(text):
"""
Extracts information about the contract data.
:param text: the string to parse.
:return: a tuple of expiry date and strike price
"""
# SPXW should expire at 3PM Chicago time.
[year, month, day, strike, option_code] = text.split(' ')
expiry = '%s %s %s 3:00PM -05:00' % (year, month, day)
dt_object = parser.parse(expiry)
"""
Third friday SPX standard options expire at start of trading
8.30 A.M. Chicago time.
"""
if is_third_friday(dt_object):
dt_object = dt_object.replace(hour=8, minute=30)
strike = Decimal(strike)
return (dt_object, strike)
实用函数parse_expiry_and_strike()返回一个到期日期对象的元组，以及一个Decimal对象作为行权价。
每个合同数据都是一个字符串，包含到期年、月、日和行权价，后跟合同代码，所有用空格分隔。我们取日期组件并重构一个日期和时间字符串，可以轻松地由之前导入的dateutil解析函数解析。周期期权在纽约时间下午 4 点到期，或芝加哥时间下午 3 点到期。标准的第三个星期五期权是上午结算的，将在交易日开始时上午 8 点 30 分到期。我们根据执行is_third_friday()检查的需要替换到期时间，实现如下：
In [ ]:
def is_third_friday(dt_object):
return dt_object.weekday() == 4 and 15 <= dt_object.day <= 21
使用一个简单的合同代码数据测试我们的函数并打印结果。
In [ ]:
test_contract_code = '2018 Sep 26 1800.00 (*)'
(expiry, strike) = parse_expiry_and_strike(test_contract_code)
In [ ]:
print('Expiry:', expiry)
print('Strike price:', strike)
Out[ ]:
Expiry: 2018-09-26 15:00:00-05:00
Strike price: 1800.00
自 2018 年 9 月 26 日起，星期三，SPXW 期权将在芝加哥当地时间下午 3 点到期。
这一次，让我们使用一个落在第三个星期五的合同代码数据来测试我们的函数：
In [ ]:
test_contract_code = '2018 Oct 19 2555.00 (*)'
(expiry, strike) = parse_expiry_and_strike(test_contract_code)
In [ ]:
print('Expiry:', expiry)
print('Strike price:', strike)
Out[ ]:
Expiry: 2018-10-19 08:30:00-05:00
Strike price: 2555.00
我们使用的测试合同代码数据是 2018 年 10 月 19 日，这是 10 月的第三个星期五。这是一个标准的 SPX 期权，将在交易日开始时，在芝加哥时间上午 8 点 30 分结算。
有了我们的实用函数，我们现在可以继续解析单个看涨或看跌期权价格条目，并返回我们可以使用的有用信息：
In [ ]:
def format_option_data(option_data):
[desc, _, _, bid_str, ask_str] = option_data[:5]
bid = Decimal(bid_str.strip() or '0')
ask = Decimal(ask_str.strip() or '0')
mid = (bid+ask) / Decimal(2)
(expiry, strike) = parse_expiry_and_strike(desc)
return (expiry, strike, bid, ask, mid)
实用函数format_option_data()以option_data作为其参数，其中包含我们之前看到的数据列表。索引零处的描述性数据包含合同代码数据，我们可以使用parse_expiry_and_strike()函数进行解析。索引三和四包含买价和卖价，用于计算中间价。中间价是买价和卖价的平均值。该函数返回期权到期日的元组，以及买价、卖价和中间价作为Decimal对象。
寻找近期和次近期期权
VIX 指数使用 24 天到 36 天到期的看涨和看跌期权的市场报价来衡量 SPX 的 30 天预期波动率。在这些日期之间，将有两个 SPX 期权合同到期日。最接近到期的期权被称为近期期权，而稍后到期的期权被称为次近期期权。每周发生一次，当期权到期日超出 24 到 36 天的范围时，新的合同到期日将被选择为新的近期和次近期期权。
为了帮助我们找到近期和次近期期权，让我们按到期日对看涨和看跌期权数据进行组织，每个期权数据都有一个以行权价为索引的pandas DataFrame。我们需要以下 DataFrame 列定义：
In [ ]:
CALL_COLS = ['call_bid', 'call_ask', 'call_mid']
PUT_COLS = ['put_bid', 'put_ask', 'put_mid']
COLUMNS = CALL_COLS + PUT_COLS + ['diff']
以下函数generate_options_chain()将我们的列表数据集calls_and_puts组织成一个单一的字典变量chain：
In [ ]:
import pandas as pd
def generate_options_chain(calls_and_puts):
chain = {}
for row in calls_and_puts:
(call, put) = row
(call_expiry, call_strike, call_bid, call_ask, call_mid) = \
format_option_data(call)
(put_expiry, put_strike, put_bid, put_ask, put_mid) = \
format_option_data(put)
# Ensure each line contains the same put and call maturity
assert(call_expiry == put_expiry)
# Get or create the DataFrame at the expiry
df = chain.get(call_expiry, pd.DataFrame(columns=COLUMNS))
df.loc[call_strike, CALL_COLS] = \
[call_bid, call_ask, call_mid]
df.loc[call_strike, PUT_COLS] = \
[put_bid, put_ask, put_mid]
df.loc[call_strike, 'diff'] = abs(put_mid-call_mid)
chain[call_expiry] = df
return chain
In [ ]:
chain = generate_options_chain(calls_and_puts)
chain变量的键是期权到期日，每个键都引用一个pandas DataFrame 对象。对format_option_data()函数进行两次调用，以获取感兴趣的看涨和看跌数据。assert关键字确保了我们的看涨和看跌到期日的完整性，基于我们数据集中的每行都指向相同的到期日的假设。否则，将抛出异常并要求我们检查数据集是否存在任何损坏的迹象。
loc关键字为特定行权价分配列值，对于看涨期权和看跌期权数据。此外，diff列包含看涨和看跌报价的中间价格的绝对差异，我们稍后将使用。
让我们查看我们的chain字典中的前两个和最后两个键：
In [ ]:
chain_keys = list(chain.keys())
for row in chain_keys[:2]:
print(row)
print('...')
for row in chain_keys[-2:]:
print(row)
Out[ ]:
2018-10-15 15:00:00-05:00
2018-10-17 15:00:00-05:00
...
2020-06-19 08:30:00-05:00
2020-12-18 08:30:00-05:00
我们的数据集包含未来两年内到期的期权价格。从中，我们使用以下函数选择我们的近期和下期到期日：
In [ ]:
def find_option_terms(chain, dt_current):
"""
Find the near-term and next-term dates from
the given indexes of the dictionary.
:param chain: dictionary object
:param dt_current: DateTime object of option quotes
:return: tuple of 2 datetime objects
"""
dt_near = None
dt_next = None
for dt_object in chain.keys():
delta = dt_object - dt_current
if delta.days > 23:
# Skip non-fridays
if dt_object.weekday() != 4:
continue
# Save the near term date
if dt_near is None:
dt_near = dt_object
continue
# Save the next term date
if dt_next is None:
dt_next = dt_object
break
return (dt_near, dt_next)
Out[ ]:
(dt_near, dt_next) = find_option_terms(chain, dt_current)
在这里，我们只是选择了到数据集时间后 23 天内到期的前两个期权。这两个期权到期日如下：
In [ ]:
print('Found near-term maturity', dt_near,
'with', dt_near-dt_current, 'to expiry')
print('Found next-term maturity', dt_next,
'with', dt_next-dt_current, 'to expiry')
Out[ ]:
Found near-term maturity 2018-11-09 15:00:00-05:00 with 24 days, 19:00:00 to expiry
Found next-term maturity 2018-11-16 08:30:00-05:00 with 31 days, 12:30:00 to expiry
近期到期日为 2018 年 11 月 9 日，下期到期日为 2018 年 11 月 16 日。
计算所需的分钟数
计算 VIX 的公式如下：
在这里，适用以下规定：
-
T[1]是到近期期权结算的年数
-
T[2]是到下期期权结算的年数
-
N[T1]是到近期期权结算的分钟数
-
N[T2]是到下期期权结算的分钟数
-
N[30]是 30 天内的分钟数
-
N[365]是一年 365 天的分钟数
让我们在 Python 中找出这些值：
In [ ]:
dt_start_year = dt_current.replace(
month=1, day=1, hour=0, minute=0, second=0)
dt_end_year = dt_start_year.replace(year=dt_current.year+1)
N_t1 = Decimal((dt_near-dt_current).total_seconds() // 60)
N_t2 = Decimal((dt_next-dt_current).total_seconds() // 60)
N_30 = Decimal(30 * 24 * 60)
N_365 = Decimal((dt_end_year-dt_start_year).total_seconds() // 60)
两个datetime对象的差异返回一个timedelta对象，其“total_seconds（）”方法以秒为单位给出差异。将秒数除以六十即可得到分钟数。一年的分钟数是通过取下一年的开始和当前年的开始之间的差异来找到的，而一个月的分钟数简单地是三十天内的秒数之和。
获得的值如下：
In [ ]:
print('N_365:', N_365)
print('N_30:', N_30)
print('N_t1:', N_t1)
print('N_t2:', N_t2)
Out[ ]:
N_365: 525600
N_30: 43200
N_t1: 35700
N_t2: 45390
计算 T 的一般公式如下：
在这里，适用以下规定：
-
M[当前日]是直到当天午夜剩余的分钟数
-
M[其他日]是当前日和到期日之间的分钟总和
-
M[结算日]是从到期日的午夜到到期时间的分钟数
有了这些，我们可以找到 T[1]和 T[2]，即近期和下期期权每年剩余的时间：
In [ ]:
t1 = N_t1 / N_365
t2 = N_t2 / N_365
In [ ]:
print('t1:%.5f'%t1)
print('t2:%.5f'%t2)
Out[ ]:
t1:0.06792
t2:0.08636
近期期权到期日为 0.6792 年，下期期权到期日为 0.08636 年。
计算前向 SPX 指数水平
对于每个合同月，前向 SPX 水平F如下所示：
在这里，选择绝对差异最小的行权价。请注意，对于 VIX 指数的计算，不考虑出价为零的期权。这表明随着 SPX 和期权的波动性变化，出价可能变为零，并且用于计算 VIX 指数的期权数量可能在任何时刻发生变化！
我们可以用“determine_forward_level（）”函数表示前向指数水平的计算，如下面的代码所示：
In [ ]:
import math
def determine_forward_level(df, r, t):
"""
Calculate the forward SPX Index level.
:param df: pandas DataFrame for a single option chain
:param r: risk-free interest rate for t
:param t: time to settlement in years
:return: Decimal object
"""
min_diff = min(df['diff'])
pd_k = df[df['diff'] == min_diff]
k = pd_k.index.values[0]
call_price = pd_k.loc[k, 'call_mid']
put_price = pd_k.loc[k, 'put_mid']
return k + Decimal(math.exp(r*t))*(call_price-put_price
df参数是包含近期或下期期权价格的数据框。 min_diff变量包含在先前的差异列中计算的所有绝对价格差异的最小值。 pd_k变量包含我们将选择的 DataFrame，其中我们将选择具有最小绝对价格差异的行权价。
请注意，出于简单起见，我们假设两个期权链的利率均为 2.17%。在实践中，近期和次期期权的利率基于美国国债收益率曲线利率的三次样条计算，或者恒定到期国债收益率（CMTs）。美国国债收益率曲线利率可从美国财政部网站www.treasury.gov/resource-center/data-chart-center/interest-rates/Pages/TextView.aspx?data=yieldYear&year=2018获取。
让我们计算近期期权的前向 SPX 水平为f1：
In [ ]:
r = Decimal(2.17/100)
In [ ]:
df_near = chain.get(dt_near)
f1 = determine_forward_level(df_near, r, t1)
In [ ]:
print('f1:', f1)
Out[ ]:
f1: 2747.596459994546094129930225
我们将使用前向 SPX 水平F作为 2747.596。
寻找所需的前向行权价格
前向行权价格是紧挨着前向 SPX 水平的行权价格，用k0表示，并由以下find_k0()函数确定：
In [ ]:
def find_k0(df, f):
return df[df.index<f].tail(1).index.values[0]
近期期权的k0值可以通过以下函数调用简单找到：
In [ ]:
k0_near = find_k0(df_near, f1)
In [ ]:
print('k0_near:', k0_near)
Out[ ]:
k0_near: 2745.00
近期前向行权价格被确定为 2745。
确定行权价格边界
在选择用于 VIX 指数计算的期权时，忽略买价为零的认购和认沽期权。对于远虚值（OTM）认沽期权，其行权价格低于k0，下限价格边界在遇到两个连续的零买价时终止。同样，对于行权价格高于k0的远虚值认购期权，上限价格边界在遇到两个连续的零买价时终止。
以下函数find_lower_and_upper_bounds()说明了在 Python 代码中找到下限和上限的过程：
In [ ]:
def find_lower_and_upper_bounds(df, k0):
"""
Find the lower and upper boundary strike prices.
:param df: the pandas DataFrame of option chain
:param k0: the forward strike price
:return: a tuple of two Decimal objects
"""
# Find lower bound
otm_puts = df[df.index<k0].filter(['put_bid', 'put_ask'])
k_lower = 0
for i, k in enumerate(otm_puts.index[::-1][:-2]):
k_lower = k
put_bid_t1 = otm_puts.iloc[-i-1-1]['put_bid']
put_bid_t2 = otm_puts.iloc[-i-1-2]['put_bid']
if put_bid_t1 == 0 and put_bid_t2 == 0:
break
if put_bid_t2 == 0:
k_lower = otm_puts.index[-i-1-1]
# Find upper bound
otm_calls = df[df.index>k0].filter(['call_bid', 'call_ask'])
k_upper = 0
for i, k in enumerate(otm_calls.index[:-2]):
call_bid_t1 = otm_calls.iloc[i+1]['call_bid']
call_bid_t2 = otm_calls.iloc[i+2]['call_bid']
if call_bid_t1 == 0 and call_bid_t2 == 0:
k_upper = k
break
return (k_lower, k_upper)
df参数是期权价格的pandas DataFrame。otm_puts变量包含虚值认沽期权数据，并通过for循环按降序迭代。在每次迭代时，k_lower变量存储当前行权价格，同时我们在循环中向前查看两个报价。当for循环由于遇到两个零报价而终止，或者到达列表末尾时，k_lower将包含下限行权价格。
在寻找上限行权价格时采用相同的方法。由于虚值认购期权的行权价格已经按降序排列，我们只需使用iloc命令上的前向索引引用来读取价格。
当我们将近期期权链数据提供给这个函数时，下限和上限行权价格可以从k_lower和k_upper变量中获得，如下面的代码所示：
In [ ]:
(k_lower_near, k_upper_near) = \
find_lower_and_upper_bounds(df_near, k0_near)
In [ ]:
print(k_lower_near, k_upper_near
Out[ ]:
1250.00 3040.00
将使用行权价格从 1,500 到 3,200 的近期期权来计算 VIX 指数。
按行权价格制表
由于 VIX 指数由平均到期日为 30 天的认购和认沽期权的价格组成，所以所选到期日的每个期权都会对 VIX 指数的计算产生一定的影响。这个影响量可以用以下一般公式表示：
在这里，T是期权到期时间，R是期权到期时的无风险利率，K[i]是第i个虚值期权的行权价格，△K[i]是K[i]两侧的半差，使得△K[i]=0.5(K[i+1]-K[i-1])。
我们可以用以下calculate_contrib_by_strike()函数来表示这个公式：
In [ ]:
def calculate_contrib_by_strike(delta_k, k, r, t, q):
return (delta_k / k**2)*Decimal(math.exp(r*t))*q
在计算△K[i]=0.5(K[i+1]-K[i-1])时，我们使用实用函数find_prev_k()来寻找K[i-1]，如下所示：
In [ ]:
def find_prev_k(k, i, k_lower, df, bid_column):
"""
Finds the strike price immediately below k
with non-zero bid.
:param k: current strike price at i
:param i: current index of df
:param k_lower: lower strike price boundary of df
:param bid_column: The column name that reads the bid price.
Can be 'put_bid' or 'call_bid'.
:return: strike price as Decimal object.
"""
if k <= k_lower:
k_prev = df.index[i-1]
return k_prev
# Iterate backwards to find put bids
k_prev = 0
prev_bid = 0
steps = 1
while prev_bid == 0:
k_prev = df.index[i-steps]
prev_bid = df.loc[k_prev][bid_column]
steps += 1
return k_prev
类似地，我们使用相同的程序来寻找K[i+1]，使用实用函数find_next_k()，如下所示：
In [ ]:
def find_next_k(k, i, k_upper, df, bid_column):
"""
Finds the strike price immediately above k
with non-zero bid.
:param k: current strike price at i
:param i: current index of df
:param k_upper: upper strike price boundary of df
:param bid_column: The column name that reads the bid price.
Can be 'put_bid' or 'call_bid'.
:return: strike price as Decimal object.
"""
if k >= k_upper:
k_next = df.index[i+1]
return k_next
k_next = 0
next_bid = 0
steps = 1
while next_bid == 0:
k_next = df.index[i+steps]
next_bid = df.loc[k_next][bid_column]
steps += 1
return k_next
有了前面的实用函数，我们现在可以创建一个名为tabulate_contrib_by_strike()的函数，使用迭代过程来计算pandas DataFrame df中可用的每个行权价的期权的贡献，返回一个包含用于计算 VIX 指数的最终数据集的新 DataFrame：
In [ ]:
import pandas as pd
def tabulate_contrib_by_strike(df, k0, k_lower, k_upper, r, t):
"""
Computes the contribution to the VIX index
for every strike price in df.
:param df: pandas DataFrame containing the option dataset
:param k0: forward strike price index level
:param k_lower: lower boundary strike price
:param k_upper: upper boundary strike price
:param r: the risk-free interest rate
:param t: the time to expiry, in years
:return: new pandas DataFrame with contributions by strike price
"""
COLUMNS = ['Option Type', 'mid', 'contrib']
pd_contrib = pd.DataFrame(columns=COLUMNS)
for i, k in enumerate(df.index):
mid, bid, bid_column = 0, 0, ''
if k_lower <= k < k0:
option_type = 'Put'
bid_column = 'put_bid'
mid = df.loc[k]['put_mid']
bid = df.loc[k][bid_column]
elif k == k0:
option_type = 'atm'
elif k0 < k <= k_upper:
option_type = 'Call'
bid_column = 'call_bid'
mid = df.loc[k]['call_mid']
bid = df.loc[k][bid_column]
else:
continue # skip out-of-range strike prices
if bid == 0:
continue # skip zero bids
k_prev = find_prev_k(k, i, k_lower, df, bid_column)
k_next = find_next_k(k, i, k_upper, df, bid_column)
delta_k = Decimal((k_next-k_prev)/2)
contrib = calculate_contrib_by_strike(delta_k, k, r, t, mid)
pd_contrib.loc[k, COLUMNS] = [option_type, mid, contrib]
return pd_contrib
生成的 DataFrame 以行权价为索引，包含三列——期权类型为看涨或看跌，买卖价差的平均值，以及对 VIX 指数的贡献。
列出我们近期期权的贡献给出以下结果：
In [ ]:
pd_contrib_near = tabulate_contrib_by_strike(
df_near, k0_near, k_lower_near, k_upper_near, r, t1)
查看结果的头部提供以下信息：
In [ ]:
pd_contrib_near.head()
这给出以下表格：
| 期权类型 | 中间值 | 贡献 | |
|---|---|---|---|
| 1250.00 | 看跌期权 | 0.10 | 0.000003204720007271874493426366826 |
| 1300.00 | 看跌期权 | 0.125 | 0.000003703679742131881579865901010 |
| 1350.00 | 看跌期权 | 0.15 | 0.000004121296305647986745661479970 |
| 1400.00 | 看跌期权 | 0.20 | 0.000005109566338124799893855814454 |
| 1450.00 | 看跌期权 | 0.20 | 0.000004763258036967708819004706934 |
查看结果的尾部提供以下信息：
In [ ]:
pd_contrib_near.tail()
这也给我们提供了以下表格：
| 期权类型 | 中间值 | 贡献 | |
|---|---|---|---|
| 3020.00 | 看涨期权 | 0.175 | 9.608028452572290489411343569E-8 |
| 3025.00 | 看涨期权 | 0.225 | 1.231237623174939828257858985E-7 |
| 3030.00 | 看涨期权 | 0.175 | 9.544713775211615220689389699E-8 |
| 3035.00 | 看涨期权 | 0.20 | 1.087233242345573774601901086E-7 |
| 3040.00 | 看涨期权 | 0.15 | 8.127448187590304540304760266E-8 |
pd_contrib_near变量包含了单个 DataFrame 中包含的近期看涨和看跌虚值期权。
计算波动性
所选期权的波动性计算如下：
由于我们已经计算了求和项的贡献，这个公式可以简单地在 Python 中写成calculate_volatility()函数：
In [ ]:
def calculate_volatility(pd_contrib, t, f, k0):
"""
Calculate the volatility for a single-term option
:param pd_contrib: pandas DataFrame
containing contributions by strike
:param t: time to settlement of the option
:param f: forward index level
:param k0: immediate strike price below the forward level
:return: volatility as Decimal object
"""
term_1 = Decimal(2/t)*pd_contrib['contrib'].sum()
term_2 = Decimal(1/t)*(f/k0 - 1)**2
return term_1 - term_2
计算近期期权的波动性给出以下结果：
In [ ]:
volatility_near = calculate_volatility(
pd_contrib_near, t1, f1, k0_near)
In [ ]:
print('volatility_near:', volatility_near)
Out[ ]:
volatility_near: 0.04891704334249740486501736967
近期期权的波动性为 0.04891。
计算下一个期权
就像我们对近期期权所做的那样，使用已经定义好的函数进行下一个期权的计算是非常简单的：
In [ ] :
df_next = chain.get(dt_next)
f2 = determine_forward_level(df_next, r, t2)
k0_next = find_k0(df_next, f2)
(k_lower_next, k_upper_next) = \
find_lower_and_upper_bounds(df_next, k0_next)
pd_contrib_next = tabulate_contrib_by_strike(
df_next, k0_next, k_lower_next, k_upper_next, r, t2)
volatility_next = calculate_volatility(
pd_contrib_next, t2, f2, k0_next)
In [ ]:
print('volatility_next:', volatility_next)
Out[ ]:
volatility_next: 0.04524308316212813982254693873
由于dt_next是我们的下一个到期日，调用chain.get()从期权链存储中检索下一个到期期权的价格。有了这些数据，我们确定了下一个到期期权的前向 SPX 水平f2，找到了它的前向行权价k0_next，并找到了它的下限和上限行权价边界。接下来，我们列出了在行权价边界内计算 VIX 指数的每个期权的贡献，从中我们使用calculate_volatility()函数计算了下一个期权的波动性。
下一个期权的波动性为 0.0452。
计算 VIX 指数
最后，30 天加权平均的 VIX 指数写成如下形式：
在 Python 代码中表示这个公式给出以下结果：
In [ ]:
def calculate_vix_index(t1, volatility_1, t2,
volatility_2, N_t1, N_t2, N_30, N_365):
inner_term_1 = t1*Decimal(volatility_1)*(N_t2-N_30)/(N_t2-N_t1)
inner_term_2 = t2*Decimal(volatility_2)*(N_30-N_t1)/(N_t2-N_t1)
sqrt_terms = math.sqrt((inner_term_1+inner_term_2)*N_365/N_30)
return 100 * sqrt_terms
用近期和下一个期权的值进行替换得到以下结果：
In [ ]:
vix = calculate_vix_index(
t1, volatility_near, t2,
volatility_next, N_t1, N_t2,
N_30, N_365)
In [ ]:
print('At', dt_current, 'the VIX is', vix)
Out[ ]:
At 2018-10-15 20:00:00-05:00 the VIX is 21.431114075693934
我们得到了 2018 年 10 月 15 日收盘时的 VIX 指数为 21.43。
计算多个 VIX 指数
对于特定交易日计算出的单个 VIX 值，我们可以重复使用定义的函数来计算一段时间内的 VIX 值。
让我们编写一个名为process_file()的函数，来处理单个文件路径，并返回计算出的 VIX 指数：
In [ ]:
def process_file(filepath):
"""
Reads the filepath and calculates the VIX index.
:param filepath: path the options chain file
:return: VIX index value
"""
headers, calls_and_puts = read_file(filepath)
dt_current = get_dt_current(headers)
chain = generate_options_chain(calls_and_puts)
(dt_near, dt_next) = find_option_terms(chain, dt_current)
N_t1 = Decimal((dt_near-dt_current).total_seconds() // 60)
N_t2 = Decimal((dt_next-dt_current).total_seconds() // 60)
t1 = N_t1 / N_365
t2 = N_t2 / N_365
# Process near-term options
df_near = chain.get(dt_near)
f1 = determine_forward_level(df_near, r, t1)
k0_near = find_k0(df_near, f1)
(k_lower_near, k_upper_near) = find_lower_and_upper_bounds(
df_near, k0_near)
pd_contrib_near = tabulate_contrib_by_strike(
df_near, k0_near, k_lower_near, k_upper_near, r, t1)
volatility_near = calculate_volatility(
pd_contrib_near, t1, f1, k0_near)
# Process next-term options
df_next = chain.get(dt_next)
f2 = determine_forward_level(df_next, r, t2)
k0_next = find_k0(df_next, f2)
(k_lower_next, k_upper_next) = find_lower_and_upper_bounds(
df_next, k0_next)
pd_contrib_next = tabulate_contrib_by_strike(
df_next, k0_next, k_lower_next, k_upper_next, r, t2)
volatility_next = calculate_volatility(
pd_contrib_next, t2, f2, k0_next)
vix = calculate_vix_index(
t1, volatility_near, t2,
volatility_next, N_t1, N_t2,
N_30, N_365)
return vix
假设我们观察了期权链数据，并将其收集到 2018 年 10 月 15 日至 19 日的 CSV 文件中。我们可以将文件名和文件路径模式定义为常量变量：
In [ ]:
FILE_DATES = [
'2018_10_15',
'2018_10_16',
'2018_10_17',
'2018_10_18',
'2018_10_19',
]
FILE_PATH_PATTERN = 'files/chapter07/SPX_EOD_%s.csv'
通过日期进行迭代，并将计算出的 VIX 值设置到一个名为'VIX'的pandas DataFrame 列中，得到以下结果：
In [ ] :
pd_calcs = pd.DataFrame(columns=['VIX'])
for file_date in FILE_DATES:
filepath = FILE_PATH_PATTERN % file_date
vix = process_file(filepath)
date_obj = parser.parse(file_date.replace('_', '-'))
pd_calcs.loc[date_obj, 'VIX'] = vix
使用head()命令观察我们的数据提供了以下结果：
In [ ]:
pd_calcs.head(5)
这给我们提供了以下表格，其中包含了 VIX 在 5 天内的数值：
| VIX | |
|---|---|
| 2018-10-15 | 21.4311 |
| 2018-10-16 | 17.7384 |
| 2018-10-17 | 17.4741 |
| 2018-10-18 | 20.0477 |
| 2018-10-19 | 19.9196 |
比较结果
让我们通过重用在之前的部分中下载的 DataFrame df_vix_data VIX 指数，提取出 2018 年 10 月 15 日至 19 日对应周的相关数值，比较计算出的 VIX 值与实际 VIX 值：
In [ ]:
df_vix = df_vix_data['2018-10-14':'2018-10-21']['5\. adjusted close']
该时期的实际 VIX 收盘价如下：
In [ ]:
df_vix.head(5)
Out [ ]:
date
2018-10-15 21.30
2018-10-16 17.62
2018-10-17 17.40
2018-10-18 20.06
2018-10-19 19.89
Name: 5\. adjusted close, dtype: float64
让我们将实际的 VIX 值和计算出的值合并到一个 DataFrame 中，并绘制它们：
In [ ]:
df_merged = pd.DataFrame({
'Calculated': pd_calcs['VIX'],
'Actual': df_vix,
})
df_merged.plot(figsize=(10, 6), grid=True, style=['b', 'ro']);
这给我们提供了以下输出：
红点中的计算值似乎非常接近实际的 VIX 值。
总结
在本章中，我们研究了波动率衍生品及其在投资者中的用途，以实现在股票和信用组合中的多样化和对冲风险。由于股票基金的长期投资者面临下行风险，波动率可以用作尾部风险的对冲工具，并替代认购期权。在美国，芝加哥期权交易所 VIX 衡量了由 SPX 期权价格隐含的短期波动率。在欧洲，VSTOXX 市场指数基于 OESX 一篮子的市场价格，并衡量了下一个 30 天内欧洲 STOXX 50 指数的隐含市场波动率。世界各地的许多人使用 VIX 作为下一个 30 天股票市场波动率的流行测量工具。为了帮助我们更好地理解 VIX 指数是如何计算的，我们研究了它的组成部分和确定其价值的公式。
为了帮助我们确定 SPX 和 VIX 之间的关系，我们下载了这些数据并进行了各种金融分析，得出它们之间存在负相关的结论。这种关系提供了一种可行的方式，通过基于基准的交易策略来避免频繁的再平衡成本。波动性的统计性质使波动率衍生品交易者能够通过利用均值回归策略、离散交易和波动率价差交易等方式获得回报。
在研究基于 VIX 的交易策略时，我们复制了单个时间段的 VIX 指数。由于 VIX 指数是对未来 30 天波动性展望的一种情绪，它由两个 SPX 期权链组成，到期日在 24 至 36 天之间。随着 SPX 的涨跌，SPX 期权的波动性也会发生变化，期权买价可能会变为零。用于计算 VIX 指数的期权数量可能会因此而改变。为了简化本章中对 VIX 计算的分解，我们假设包括的期权数量是静态的。我们还假设了在 5 天内 CMT 是恒定的。实际上，期权价格和无风险利率是不断变化的，VIX 指数大约每 15 秒重新计算一次。
在下一节中，我们将建立一个算法交易平台。
标签：精通,Python,数据,第二,next,df,2018,VIX,期权 From： https://www.cnblogs.com/apachecn/p/18140508