展开讲讲：隐含波动率（IV）、恐慌指数（VIX）、广义恐慌 ...

0.预览

先上结果。书接上回，本文将在上次的基础上，介绍恐慌指数（vix）相关内容。具体地，本文完成了以下几个方面的内容：

梳理了恐慌指数（VIX）的定义及计算方法，讲述了VIX指数计算方法的来源；
基于第一篇文章提取的数据，给出了相应的Python计算代码。
从数据层面阐述了隐含波动率（IV）、恐慌指数（VIX）二者之间的区别。

图0：上证50ETF 期权的VIX指数

本文将用到以下数据，数据提取方法见上文。如果需要实际计算，请前往GitHub下载（将用到 opt_basic.csv，opt_price.csv，interest_rate.csv）：
本文将用到以下三个Python第三方库，可以使用如下命令安装：
$ pip install numpy
$ pip install pandas
$ pip install matplotlib1.VIX指数定义及计算方法

1.1.VIX指数定义

我们来看 Investopedia 给VIX指数下的定义：

The CBOE Volatility Index, or VIX, is a real-time market index representing the market’s expectations for volatility over the coming 30 days. Investors use the VIX to measure the level of risk, fear, or stress in the market when making investment decisions.
CBOE波动性指数，即VIX，是一个实时市场指数，代表市场对未来30天波动性的预期。投资者在做投资决定时，使用VIX指数来衡量市场的风险、恐惧或压力水平。

这个定义相当清楚了，VIX指数就是一个指数，衡量市场未来30天的波动率预期，也反映市场恐惧、压力水平，所以也称为恐慌指数。
VIX指数的含义已经有相当多的朋友写文章探讨过了，写的也非常清晰详细。所以，这并非本文的讨论重点。本文的重点放在VIX指数的计算实现过程。请看下一小节。
1.2.VIX指数计算方法

芝加哥商品交易所CBOE于1993年推出VIX指数，我国上交所也曾推出过中国波指（IVIX.SH，后来停止了）。二者都给出了计算方法（为保证网页链接长期有效，本文直接将链接指向github仓库）：

CBOE给出了 CBOE VIX white paper，提供了计算的手把手教学；
上交所给出了上证 50 ETF波动率指数编制方案，给出了中国波指的计算方法。

上述链接点进去就是相应的文件，这两个阐述计算方法的文件都指向了同一个计算方法：

其中：

：近月波动率估计；
：期权合约剩余到期时间；
：时长为的无风险利率；
：无套利的远期价格，计算步骤如下。
- 获取时长为的无风险利率；
- 获取认购期权价格与认沽期权价格相差最小的执行价，记为；（注意这里不是标的资产价格）
- 计算无套利远期价格，，是执行价为S的认购期权价格。
：小于但最靠近的执行价；
：从小到大所有近月合约的执行价；（）
：第个执行价所对应的执行价间隔，一般为 ;
：执行价为的虚值期权的价格。（即有可能为认购期权，也有可能是认沽期权，保证它为虚值期权！）至于为什么是？quote、报价，比较合理吧？
：代表时间间隔数，分别代表近月（下一个月），次近月（下下个月）。

1.3.计算公式的来源

如果不想看公式推导，可以直接跳过这个小结。
看到这个公式，不禁要问，这个公式是怎么来的？为什么它能够计算出未来的波动率预期？带着这个疑问，我看到了这篇文章：《Does VIX Truly Measure Return Volatility?》。文章详细信息如下：

Chow, Victor and Jiang, Wanjun and Li, Jingrui, Does VIX Truly Measure Return Volatility? (January 22, 2018). Handbook of Financial Econometrics, Mathematics, Statistics, and Machine Learning, Forthcoming

文章详细推导了VIX计算公式，同时提出了一种更为准确的波动率度量方法——GVIX（generalized volatility index），下篇文章将会提到GVIX的Python实现。欢迎各位持续关注！
文章一开始从股价的几何布朗运动(geometric Brownian motion)和伊藤公式(Ito Calculus)开始：

公式（3）与公式（4）联立，将（3）-（4），可以解出波动率 :

然后，对时间求积分，即可将波动率从上述公式中剥离：

所谓VIX指数，就是对未来一个月的波动率预期：

（7）式我们先放一边，以下给出持有期收益（Holding-Period Return）和对数收益（Log-Return）的定义：

对期的对数股价做含有余项的泰勒展开，从这一步开始，就决定了VIX的计算必然基于虚值期权：

持有期收益与对数收益之差可以表示为下式，然后对其求期望：

我们得到了预期对数收益：

在无套利条件（No-arbitrage）下，我们容易有：

联立（7）、（11）、（12），得到VIX的计算表达式：

用下考研时学过的等价无穷小（武老师：不然你就慢了）：

最后，我们就得到了VIX的表达式：

将表达式离散化，就得到式（1）。这样，我们就弄懂了这个VIX计算表达来源。如果不懂的话我们就假装自己懂了（不懂原理不影响我们的计算），然后进入下一个章节。
2.VIX指数Python实现

使用如下命令导入需要用到的第三方库，后文将不再赘述导入过程：
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt为避免反复前后翻动，我将（1）（2）式搬到这里：

将下载好的数据（csv文件）放在&#34;data/&#34;文件夹下，然后着手实现VIX指数计算。
2.1.数据预处理

首先，我们提取计算需要使用的数据，并做一些整理。一是将期权合约信息和期权行情数据整合到一起，形成一个dataframe，二是处理一下时间方面的格式。具体代码见下方prepare_data函数：
def prepare_data():
&#34;&#34;&#34;提取VIX计算需要的数据&#34;&#34;&#34;
# 期权数据
opt_basic = pd.read_csv(&#39;data/opt_basic.csv&#39;)
opt_price = pd.read_csv(&#39;data/opt_price.csv&#39;)
opt_data = pd.merge(opt_price, opt_basic, on=&#39;code&#39;)
# 转换为日期格式
opt_data[&#39;date&#39;] = pd.to_datetime(opt_data[&#39;date&#39;])
opt_data[&#39;exercise_date&#39;] = pd.to_datetime(opt_data[&#39;exercise_date&#39;])
# 筛选需要的字段
use_fields = [&#39;date&#39;, &#39;close&#39;, &#39;contract_type&#39;, &#39;exercise_price&#39;, &#39;exercise_date&#39;]
opt_data = opt_data[use_fields]
# 重命名
opt_data[&#39;contract_type&#39;] = opt_data[&#39;contract_type&#39;].map({&#39;CO&#39;: &#39;call&#39;, &#39;PO&#39;: &#39;put&#39;})
opt_data[&#39;T-days&#39;] = (opt_data[&#39;exercise_date&#39;] - opt_data[&#39;date&#39;]).apply(lambda x: x.days)
opt_data[&#39;T&#39;] = opt_data[&#39;T-days&#39;] / 365
# 利率数据
interest_rate = pd.read_csv(&#39;data/interest_rate.csv&#39;, index_col=0)
# 转换格式
interest_rate.index = pd.DatetimeIndex(interest_rate.index)
interest_rate.columns = interest_rate.columns.astype(int)
# 填充遗漏天数
omit_date = list(set(opt_data[&#39;date&#39;]) - set(interest_rate.index))
for d in omit_date:
interest_rate.loc[d, :] = np.nan
interest_rate = interest_rate.sort_index().ffill()
return opt_data, interest_rate

opt_data, interest_rate = prepare_data()代码运行结果见下图：

图1：opt_data数据

图2：interest_rate数据

一种比较可行的思路是按照交易日逐日计算VIX指数。我们先取出一个交易日的数据，完整地计算下来，以此类推，就能够计算出VIX指数的时间序列。使用下列代码获取一个交易日的数据样本：
for dt, data in opt_data.groupby(&#39;date&#39;):
break
# dt == Timestamp(&#39;2015-02-09 00:00:00&#39;)
# data见下图
data.head(10)

图3：一个交易日数据的样例

我们先对合约进行筛选，选出近月（Front-month，简写为fm）、次近月（Sub-front-month，简写为sfm）的期权合约。其中T-days为合约的到期时间，我们需要筛选出，合约存续时间大于7天的，最小的两个存续时间，分别为近月和次近月。代码如下：
# 筛选近月、次近月的合约存续时间
all_T = data[data[&#39;T-days&#39;] > 7][&#39;T-days&#39;].unique()
all_T.sort()
# 以【日】为单位和以【年】为单位的存续时间
Tdays_fm, Tdays_sfm = all_T[:2]
T_fm, T_sfm = all_T[:2] / 365
# 以此筛选出用于计算的合约数据
data_fm, data_sfm = data[data[&#39;T&#39;] == T_fm], data[data[&#39;T&#39;] == T_sfm]
data_fm

图4：data_fm数据展示

2.2.参数计算

一通操作猛如虎，这里我们直接一次性求出所有计算参数。给定上述合约信息，我们可以求得所有计算当日VIX指数的参数。但其中，数据透视表的操作尤其重要：
# 数据透视表的使用示例
data_fm_ = data_fm.pivot_table(index=&#39;exercise_price&#39;, columns=&#39;contract_type&#39;, values=&#39;close&#39;)
data_fm_[&#39;call-put&#39;] = np.abs(data_fm_[&#39;call&#39;] - data_fm_[&#39;put&#39;])
data_sfm_ = data_sfm.pivot_table(index=&#39;exercise_price&#39;, columns=&#39;contract_type&#39;, values=&#39;close&#39;)
data_sfm_[&#39;call-put&#39;] = np.abs(data_sfm_[&#39;call&#39;] - data_sfm_[&#39;put&#39;])

图5：数据透视表样例

如前文所述，我们计算VIX指数需要如下参数，我们分别来讲讲计算思路：

：时长为的无风险利率；我们直接用合约剩余时长 和日期 dt 即可取得。
：无套利的远期价格，计算步骤如下。
获取时长为的无风险利率；
获取认购期权价格与认沽期权价格相差最小的执行价，记为；（注意这里不是标的资产价格）
计算无套利远期价格，，是执行价为S的认购期权价格。
：小于但最靠近的执行价；
：从小到大所有近月合约的执行价；（）
：第个执行价所对应的执行价间隔，一般为 ;但是，处于两端的数据没法计算，所以我们采取填充的方式。这个数据的计算方法为移动窗口法。
：执行价为的虚值期权的价格。（即有可能为认购期权，也有可能是认沽期权，保证它为虚值期权！）

def pivot_table_(data):
data = data.pivot_table(index=&#39;exercise_price&#39;, columns=&#39;contract_type&#39;, values=&#39;close&#39;)
data[&#39;call-put&#39;] = data[&#39;call&#39;] - data[&#39;put&#39;]
data[&#39;abs-diff&#39;] = np.abs(data[&#39;call-put&#39;])
return data

def get_params(data, dt, Tdays):
&#34;&#34;&#34;
参数计算函数.
data: 前面计算的数据透视表；
dt: 计算VIX指数的时间；
Tdays: 到期时间，单位为（日）；
return:
      S: 无套利标的资产价格
      R: 对应时间的无风险利率
      F: 无套利远期价格
      K0: 离F最近的，但小于F的执行价
      Q: pandas.DataFrame, 里面包含计算参数Q, K_i, delta_K_i
&#34;&#34;&#34;
S = data.sort_values(by=&#39;abs-diff&#39;).index[0]
R = interest_rate.loc[dt, Tdays] / 100
F = S + np.exp(Tdays / 365 * R) * data.loc[S, &#39;call-put&#39;]
K0 = np.max(data.index[data.index < F])
# 合约选择
Q = data[[&#39;call&#39;, &#39;put&#39;]].copy()
Q.loc[Q.index >= F, &#39;Q&#39;] = Q.loc[Q.index >= F, &#39;call&#39;]
Q.loc[Q.index < F, &#39;Q&#39;] = Q.loc[Q.index < F, &#39;put&#39;]
Q = Q[&#39;Q&#39;].reset_index()
# 一大波操作：重设K_i index，移动窗口方法计算delta_K_i
Q[&#39;delta_K&#39;] = Q[&#39;exercise_price&#39;]\
      .rolling(3, center=True)\
      .apply(lambda x: (x.iloc[-1] - x.iloc[0]) / 2)\
      .ffill().bfill() # 填充两端的缺失值
return S, R, F, K0, Q

data_fm_ = pivot_table_(data_fm.copy())
data_sfm_ = pivot_table_(data_sfm.copy())

S_fm, R_fm, F_fm, K0_fm, Q_fm = get_params(data_fm_, dt, Tdays_fm)
S_sfm, R_sfm, F_sfm, K0_sfm, Q_sfm = get_params(data_sfm_, dt, Tdays_sfm)

# S_fm, R_fm, F_fm, K0_fm = (2.3, 0.04996766666666667, 2.325754666938786, 2.3)
# S_sfm, R_sfm, F_sfm, K0_sfm = (2.3, 0.049402999999999996, 2.3140361211554956, 2.3)2.3.VIX指数计算

根据上述计算的参数，以及VIX的计算表达式，我们可以比较容易地计算出单日的VIX指数（如果算错一定要告诉我，谢谢谢谢！！！）：
def get_sigma(T, R, F, K0, Q):
&#34;&#34;&#34;
计算式（1）sigma的函数。至于为什么拆分成两部分，或许是“激动的心，颤抖的手”。
&#34;&#34;&#34;
sigma_part1 = (2/T) * np.sum(np.exp(R*T)*Q[&#39;Q&#39;]*Q[&#39;delta_K&#39;]/np.square(Q[&#39;exercise_price&#39;]))
sigma_part2 = (1/T) * (F/K0 - 1)**2
return sigma_part1 + sigma_part2

sigma1 = get_sigma(T_fm, R_fm, F_fm, K0_fm, Q_fm)
sigma2 = get_sigma(T_sfm, R_sfm, F_sfm, K0_sfm, Q_sfm)

# sigma1, sigma2 = (0.06465976520531855, 0.054218952375645016)

def get_VIX(sigma1, sigma2, T1, T2):
&#34;&#34;&#34;计算式（2）VIX计算函数，同样激动。&#34;&#34;&#34;
vix_part1 = (T1 * sigma1) * (T2 - 30 / 365) / (T2 - T1)
vix_part2 = (T2 * sigma2) * (30 / 365 - T1) / (T2 - T1)
return np.sqrt((365 / 30) * (vix_part1 + vix_part2))

vix = get_VIX(sigma1, sigma2, T_fm, T_sfm)
# vix = 0.2778286173181711至此，我们就能计算VIX指数了。我们需要将代码整合一下，计算所有时间的VIX指数。我们回到最初的那个for循环，开始编写我们的代码：
# 设置一个空字典来储存数据
VIX = {}
for dt, data in opt_data.groupby(&#39;date&#39;):
# 筛选近月、次近月的合约存续时间
all_T = data[data[&#39;T-days&#39;] > 7][&#39;T-days&#39;].unique()
all_T.sort()
T_fm, T_sfm = all_T[:2] / 365
Tdays_fm, Tdays_sfm = all_T[:2]
# 以此筛选出用于计算的合约数据
data_fm, data_sfm = data[data[&#39;T&#39;] == T_fm], data[data[&#39;T&#39;] == T_sfm]
# 数据透视表操作
data_fm_ = pivot_table_(data_fm.copy())
data_sfm_ = pivot_table_(data_sfm.copy())
# 获取计算参数
S_fm, R_fm, F_fm, K0_fm, Q_fm = get_params(data_fm_, dt, Tdays_fm)
S_sfm, R_sfm, F_sfm, K0_sfm, Q_sfm = get_params(data_sfm_, dt, Tdays_sfm)
# 计算sigma
sigma1 = get_sigma(T_fm, R_fm, F_fm, K0_fm, Q_fm)
sigma2 = get_sigma(T_sfm, R_sfm, F_sfm, K0_sfm, Q_sfm)
# 计算VIX指数
VIX[dt] = get_VIX(sigma1, sigma2, T_fm, T_sfm)

# 里面有一些天数算不出来（但不多），估计是合约移仓换月的时候无数可算。
VIX = pd.Series(VIX, name=&#39;VIX&#39;).dropna()

# 保存数据
VIX.to_csv(&#39;source/VIX.csv&#39;)我们来做一些简单的数据可视化。代码如下：
fig, ax = plt.subplots(figsize=(15, 6))
VIX.plot(ax=ax)
ax.grid()
# 保存图像
# plt.savefig(&#39;source/VIX.png&#39;, bbox_inches=&#39;tight&#39;, dpi=200)

图6：VIX指数

3.VIX指数与IV的联系与区别

在github页面下载计算好的IV、VIX序列，放在&#34;data/&#34;文件夹下，然后进行下列分析。首先，读取数据合并到同一个dataframe中：
IV = pd.read_csv(&#39;data/IV.csv&#39;, index_col=0)
IV = IV.groupby(&#39;date&#39;).mean()
IV.index.name = None
IV.index = pd.DatetimeIndex(IV.index)
VIX = pd.read_csv(&#39;data/VIX.csv&#39;, index_col=0)
VIX.index = pd.DatetimeIndex(VIX.index)
df = pd.concat([IV, VIX], axis=1)
df

图7：合并的dataframe

参照matplotlib官方文档，画一个VIX-IV散点图：
def scatter_hist(x, y, ax, ax_histx, ax_histy):
# no labels
ax_histx.tick_params(axis=&#34;x&#34;, labelbottom=False)
ax_histy.tick_params(axis=&#34;y&#34;, labelleft=False)

# the scatter plot:
ax.scatter(x, y, alpha=.5)

# now determine nice limits by hand:
binwidth = 0.01
xymax = max(np.max(np.abs(x)), np.max(np.abs(y)))
lim = (int(xymax/binwidth) + 1) * binwidth

bins = np.arange(0.1, lim + binwidth, binwidth)
ax_histx.hist(x, bins=bins)
ax_histy.hist(y, bins=bins, orientation=&#39;horizontal&#39;)
pass

# definitions for the axes
left, width = 0.1, 0.65
bottom, height = 0.1, 0.65
spacing = 0.0005

rect_scatter = [left, bottom, width, height]
rect_histx = [left, bottom + height + spacing, width, 0.2]
rect_histy = [left + width + spacing, bottom, 0.2, height]

# start with a square Figure
fig = plt.figure(figsize=(12, 12))

ax = fig.add_axes(rect_scatter)
ax_histx = fig.add_axes(rect_histx, sharex=ax)
ax_histy = fig.add_axes(rect_histy, sharey=ax)

x, y = df[&#39;IV&#39;], df[&#39;VIX&#39;]
# use the previously defined function
scatter_hist(x, y, ax, ax_histx, ax_histy)
ax.plot([0.1, 0.6], [0.1, 0.6], color=&#39;k&#39;, lw=2)
ax.set_xlabel(&#39;$IV$&#39;, fontsize=20)
ax.set_ylabel(&#39;$VIX$&#39;, fontsize=20)
# 保存图片
plt.savefig(&#39;source/scatter_plot.png&#39;, dpi=200, bbox_inches=&#39;tight&#39;)

图8：IV与VIX对比散点图

图中的黑色直线为。从图中不难看出：

VIX指数略低于IV指数，但整体围绕在45度线左右；
VIX指数与IV指数的分布十分相近。

我们再考察一下。计算并绘制出折线图：
fig, ax = plt.subplots(figsize=(15, 6))
df[&#39;spread&#39;] = df[&#39;IV&#39;] - df[&#39;VIX&#39;]
df[&#39;spread&#39;].plot()
ax.axhline(y=0, ls=&#39;--&#39;, c=&#39;0&#39;, alpha=0.8)
# 保存图像
plt.savefig(&#39;source/IV-VIX-spread.png&#39;, dpi=200, bbox_inches=&#39;tight&#39;)

图9：IV-VIX spread

从图中可以看出，在二者波动性较大的时间，二者之差波动性也较大。当然还有个刚刚已经得出的结论，就是IV往往比VIX要更大。
从数据上来看，IV与VIX相差并不大，其原因可能是他们两个的计算方法都源自于伊藤引理，并都基于收益率的正态假设，对对数股价的泰勒展开也只是止步于二阶。所以，从一定程度上看应该是同一个东西。但是区别可能在于波动率汇总的过程：IV序列是所有期权的IV按日取平均，但VIX指数则采用近月与次近月加权平均的方法。
上述只是我的一些想法，很有可能是错的，如果有正确的观点请告诉我！
4.总结

亲爱的读者，感谢你看到这里！本文阐述了VIX指数的由来、Python实现，以及VIX指数与IV在数值层面的区别。我也只能做到数值层面，也就是把它们都算出来，然后做比较。原理层面的区别还望大佬给出解释，我的想法可能不太靠谱。
这是本系列的第二篇文章了，和第一篇的间隔有些久，因为这对我来说确实是挑战了。如果文章中有错误，请一定要告诉我！最后，再次感谢大家的阅读！