0.预览
先上结果，本篇完成了隐含波动率相关内容（以上证vix.shtml" target="_blank" class="relatedlink">50ETF期权为例），具体为以下几个方面：
- 梳理了隐含波动率的定义，简要介绍了它的计算方法；
- 介绍了计算数据获取方法（上证50ETF行情数据、上证50ETF期权数据、无风险利率数据（SHIBOR）），并给出了数据处理相关代码。
- 给出了隐含波动率的Python数值计算方法，包括BSM模型实现，牛顿法实现；
- 给出了相应的数据可视化方案。
本来想一篇文章把这三个指标讲完，在跑通代码后发现太长了，所以打算写成一个系列，这篇是第一篇。如果对此话题感兴趣，欢迎各位持续关注。
![]()
图0：上证50ETF期权隐含波动率曲线
工欲善其事，必先利其器。在开始动手之前需要提前安装一些依赖包（其中jqdatasdk是聚宽的数据接口，用于提取相关数据。聚宽数据接口使用权限申请请点这里）：
$ pip install jqdatasdk
$ pip install pandas
$ pip install numpy
$ pip install scipy
$ pip install matplotlib
$ pip install tqdm完成相关依赖包安装后，在python中使用下列代码完成导入，后续代码将不再赘述导入过程：
from jqdatasdk import * # 数据接口
import pandas as pd
import numpy as np # 计算
from tqdm import tqdm as tq # 程序中加个进度条
import matplotlib.pyplot as plt # 绘图
from scipy.stats import norm # 正态分布
# 登录聚宽账号，便于后续数据提取
auth('UserID', 'Password')1.数据搜集与预处理
1.1.数据获取与保存
本步骤要获取4个数据文件， 时间跨度为2015-01-01至2022-03-30（我写代码的时候）：
- 历史SHIBOR数据，利用线性插值法获取每日对应的利率期限结构，并保存为 interest_rate.csv；
- 期权基本信息，包含期权代码、上市日期、行权日、执行价格等等。保存为 opt_basic.csv；
- 期权行情数据，包含期权代码、高开低收、结算价、交易日等等。保存为 opt_price.csv；
- 标的资产数据。包含上证50ETF行情数据。保存为 etf_price.csv。
用于本文演示的数据已上传github，有需要练习的朋友可前往下载（顺便点个赞啊谢谢！）：
聚宽数据接口非常好用，有条件的朋友可以支持一下。我们来看看如何提取并保存数据。
首先，提取利率数据。
start_date = '2015-01-01'
df = pd.DataFrame()
# 大致计算了以下，每次返回5000行数据，三次能把数据提取完毕
for _ in range(3):
q = query(macro.MAC_LEND_RATE)\
.filter(macro.MAC_LEND_RATE.market_id == '5')\
.filter(macro.MAC_LEND_RATE.day > start_date)
df_ = macro.run_query(q).sort_values(by='day').drop(columns=['id', 'currency_id', 'market_id', 'currency_name'])
df = pd.concat([df, df_])
start_date = df_['day'].iloc[-1]
# 代号对应天数
days = {1: 30, 3: 90, 6:180, 7: 7, 12: 360, 14: 14, 20: 1, 23: 270}
df['days'] = df['term_id'].map(days)
df
![]()
图1：代码运行结果。接下来要对其进行一些处理。
对利率数据进行处理，使用线性插值方法，得到每一天对应的利率期限结构：
df_r = pd.DataFrame()
for day, tmp in df.groupby('day'):
tmp = tmp[['days', 'interest_rate']].set_index('days').to_dict()['interest_rate']
tmp = pd.Series(tmp, index=range(1, 361), name=day).interpolate('linear')
df_r = pd.concat([df_r, tmp], axis=1)
df_r = df_r.T
# 这里得到我们第一个数据文件 interest_rate.csv
# df_r.to_csv('interest_rate.csv')
df_r
![]()
图2：代码运行结果，利率的期限结构。
然后，提取期权基本信息数据。我们以上证50ETF为例，所以只提取上证50ETF期权基本信息（规定底层资产代码为 510050.XSHG）：
q = query(opt.OPT_CONTRACT_INFO).filter(opt.OPT_CONTRACT_INFO.underlying_symbol=='510050.XSHG')
opt_basic = opt.run_query(q)
# 这里得到我们第二个数据文件 opt_basic.csv
# opt_basic.to_csv('opt_basic.csv', index=False)
opt_basic
![]()
图3：代码运行结果，期权的基本信息。
然后，基于期权的基本信息，按照期权代码提取每个期权的行情信息，并做数据表格的合并：
opt_price = pd.DataFrame()
for code in tq(opt_basic['code']):
tmp = opt.run_query(query(opt.OPT_DAILY_PRICE).filter(opt.OPT_DAILY_PRICE.code==code))
opt_price = pd.concat([opt_price, tmp])
opt_price.reset_index()
# 这里得到我们第三个数据文件 opt_price.csv
# opt_price.to_csv('opt_price.csv', index=False)
opt_price.reset_index(drop=True)
![]()
图4：代码运行结果，期权行情数据。
最后，我们来提取上证50ETF的行情数据并作保存：
etf_price = get_price('510050.XSHG', start_date='2015-01-01', end_date='2022-03-30', fq=None)
# 这里得到我们第四个数据文件 etf_price.csv
# etf_price.to_csv('etf_price.csv')
etf_price
![]()
图5：代码运行结果，etf行情数据。
至此，我们就将计算用的所有数据都提取完毕了。下节将展示数据预处理的过程。
1.2.数据预处理
首先，加载期权基本信息数据，从中获取需要的字段，形成字典，用于后续的map方法：
opt_basic = pd.read_csv('opt_basic.csv')
# 有用的字段
fields = ['code', 'trading_code', 'contract_type', 'exercise_price', 'list_date', 'delist_date','exercise_date']
# 将用于计算的字段转换成dict，用于map方法。
exe_date = opt_basic[fields].set_index('code')['exercise_date'].to_dict()
exe_price = opt_basic[fields].set_index('code')['exercise_price'].to_dict()
opt_type = opt_basic[fields].set_index('code')['contract_type'].apply(lambda x: 'call' if x=='CO' else 'put').to_dict()然后，读取期权行情数据，并将这些字段对应到行情数据之中：
opt_price = pd.read_csv('opt_price.csv')
# 提取计算需要的字段
fields = ['code', 'date', 'close']
opt_price = opt_price[fields].reset_index(drop=True)
# 时间格式处理
opt_price['exercise_date'] = opt_price['code'].map(exe_date)
opt_price['date'] = pd.to_datetime(opt_price['date'])
opt_price['exercise_date'] = pd.to_datetime(opt_price['exercise_date'])
opt_price['K'] = opt_price['code'].map(exe_price)
opt_price['opt_type'] = opt_price['code'].map(opt_type)
opt_price = opt_price.dropna()
# 计算距离行权日的时间
opt_price['T-days'] = (opt_price['exercise_date'] - opt_price['date']).apply(lambda x: x.days)
opt_price['T'] = opt_price['T-days'] / 360
opt_price = opt_price[opt_price['T-days'] > 0]
# 开始对应利率
opt_price['r'] = np.nan
opt_price我们来检查以下代码运行的中间结果：
![]()
图6：代码运行中间结果
然后，我们应当读取利率期限结构表格和ETF行情数据，在分别对应到上述表格之中。先看利率数据是如何对应的：
df_r = pd.read_csv('interest_rate.csv', index_col=0)
# 对上时间
for dt in tq(set(opt_price['date'])):
try:
df_r.loc[dt.strftime("%F"), :]
except:
df_r.loc[dt.strftime("%F"), :] = np.nan
# 其中有缺失值，做向下填充处理
df_r = df_r.sort_index().ffill()
df_r.columns = df_r.columns.astype(int)
for i, row in tq(opt_price.iterrows()):
dt, tdays = row['date'].strftime("%F"), row['T-days']
opt_price.loc[i, 'r'] = df_r.loc[dt, tdays] / 100ETF数据对应：
etf_price = pd.read_csv('etf_price.csv', index_col=0)
for dt in tq(set(opt_price['date'])):
try:
etf_price.loc[dt.strftime("%F"), :]
except:
etf_price.loc[dt.strftime("%F"), :] = np.nan
# etf_price = etf_price.sort_index().ffill()
etf_price.index = pd.DatetimeIndex(etf_price.index)
s = etf_price['close'].to_dict()
opt_price['S'] = opt_price['date'].map(s)
opt_price = opt_price.dropna()
opt_price数据对应的最终结果展示如下：
![]()
图7：代码运行结果，所有数据对应上了。
至此，已经搜集了计算隐含波动率所有所需的数据。隐含波动率的计算方法及代码，请接着阅读下一节。
2.隐含波动率（IV）
2.1.BSM公式
根据BSM模型 ，认购期权的理论价格是：
![]()
![]()
![]()
看涨看跌平价（Put-Call Parity）：
![]()
核心参数有：
![]() 和 ![]() ：看涨和看跌期权价格
![]() ：期权执行价
![]() ：无风险利率
![]() ：期权剩余时间
![]() ：标的资产价格
![]() ：标的资产波动率
上述核心参数中，唯有 ![]() 我们没有办法获得直接的数据。但是，可以通过解方程的方法把 ![]() 解出来。解出来的这个波动率就称为隐含波动率（IV）。即，有如下方程：
![]()
或者：
![]()
我们先着手实现BSM公式：
def N(x):
return norm._cdf(x)
def bsm_model(S, K, r, T, sigma, opt_type='call'):
"""BSM模型，参数符号和文中一致"""
d1 = (np.log(S/K)+(r+0.5*sigma**2)*T) / (sigma*np.sqrt(T))
d2 = d1 - sigma*np.sqrt(T)
Call = S*N(d1) - np.exp(-r*T)*K*N(d2)
if opt_type == 'call':
return Call
elif opt_type == 'put':
Put = Call + K*np.exp(-r*T) - S
return Put2.2.牛顿法迭代求解
如何解这样的超越方程呢？解析解是不可能有的，但我们可以求得它的数值解。这里我们采用牛顿法。
假设我们要解这样一个方程 ![]() 。对于任意 ![]() 满足：
![]() ，其中 ![]() 为区间 ![]() 。
- 零点
![]() 在区间 ![]() 中。
- 对于任意
![]() , ![]() 是连续的。
我们可以利用以下迭代来求解出零点（起码足够接近）：
![]()
在我们的例子中，只需令 ![]() ，求解方程 ![]() 即可。
如何求得数值导数？假设有一个连续可导的函数 ![]() 和一个足够小的正数 ![]() , 在某点 ![]() 附近的导数可以计算为：
![]()
来看看代码实现：
def get_derivative(x, func, epsilon=1e-8):
"""获取一个函数在某点的导数"""
return (func(x + epsilon) - func(x - epsilon)) / (2*epsilon)
def Newton_method(x, func, loss=1e-6, iter_limit=10):
"""
牛顿法迭代求解。
x: 给定函数初始值；
func: 求零点的函数;
loss: 迭代退出阈值，默认为1e-6；
iter_limit: 迭代次数上限，达到上限仍未收敛则返回np.nan。默认值为10。
return: 函数零点数值解。
"""
l = func(x)
count = 0
while np.abs(l) > loss:
derivative = get_derivative(x, func)
x -= func(x) / derivative
l = func(x)
count += 1
if count >= iter_limit:
return np.nan
return x
f = lambda sigma: bsm_model(3.516, 2.60, 0.023391, 0.155556, sigma) - 1
Newton_method(2, f)
# 0.7685290263718533运用上述函数，我们可以轻松实现隐含波动率的迭代求解。有人可能会问为什么不用结算价，改一下就可以，结果差不多。来看看迭代求解的代码：
data = opt_price.copy()
for i, row in tq(opt_price.iterrows()):
opt_price_, S_, K_, r_, T_, opt_type_ = row['close'], row['S'], row['K'], row['r'], row['T'], row['opt_type']
f = lambda sigma: bsm_model(S_, K_, r_, T_, sigma, opt_type=opt_type_) - opt_price_
iv = Newton_method(0.2, f)
data.loc[i, 'IV'] = iv
# 有些数值不收敛，就跑到无穷大了，或许是流动性原因？
data = data[data['IV'].apply(lambda x: ~np.isinf(x))]
data
![]()
图8：隐含波动率求解结果
至此，我们就完成了隐含波动率的计算。接下来把它们画出来看看？
2.3.数据可视化
先看看波动率会不会微笑。我就取数据的最后一天绘制出来看看？
# 认购期权的波动率微笑
data.query("(date=='2022-03-30')&(opt_type=='call')")[['K', 'IV']].groupby('K').mean().plot()
![]()
图9：认购期权的波动率耐克笑。这是对的，右边应该比左边高。
接着，把整体取平均，按日期把历史的隐含波动率绘制出来：
fig, ax = plt.subplots(figsize=(15, 6))
data[['date', 'IV']].groupby('date').mean().plot(ax=ax)
ax.grid()
# 保存图像
# plt.savefig('source/IV-历史曲线.png', bbox_inches='tight', dpi=200)
![]()
图10：历史隐含波动率
3.总结
本文实现了隐含波动率的原理阐述与数值计算，代码纯手码，难免有错漏之处。如果发现哪里有错误，请一定要告诉我！另外，这是一个系列文章的第一篇，如果各位喜欢这个系列，欢迎各位关注！
亲爱的读者，感谢你看到这里！不知不觉这是我在知乎写的第三篇博客了，如果你觉得喜欢、有趣，欢迎与我交流！（华工解封了，结束了3个星期的封锁，感激不尽！） |
|
转播
