使用python计算期权隐含波动率并作图

## 前言

这其实是我们一次课程作业，以上证50ETF期权为例说明波动率微笑现象。按习惯我先上网搜了一下看有没有前辈写过这样的代码，毕竟重复造轮子不好嘛。没想到真的有，原文链接：[https://www.jianshu.com/p/e73f538859df](https://www.jianshu.com/p/e73f538859df)。但是这份代码有个问题，就是需要自己手动搜集数据，而且输出的数据不是标准的DataFrame。趁着做作业的机会，我借鉴并改写了作者的代码，主要实现了以下改进：

* 使用plotly作图，生成可交互式图像。
* 利用tushare自动拉取数据，计算某日所有50ETF期权的隐含波动率，同时可以生成标准的DataFrame数据。
* 借鉴[这篇文章](https://www.jianshu.com/p/34515ec65a53)的思路，利用看涨看跌期权平价关系计算隐含分红率
* 可生成波动率微笑图和波动率曲面图。

详细代码见：[https://github.com/caly5144/shu-s-project/tree/master/options](https://github.com/caly5144/shu-s-project/tree/master/options)

## 理论

从期权定价模型本身来说，公式中的波动率指的是未来的波动率数据，这使历史波动率始终存在着较大的缺陷。为了回避这一缺陷，一些学者将目光转向隐含波动率。隐含波动率具体计算思路是：不论用何种理论的理论模型对市场上交易的期权用模型确定其合理的期权价值，一般需要用到6个标准参数：执行价格（K）、到期时间（t）、标的资产价格（S）、利率（rf）、标的资产息票率（q）和波动率。前五个参数可以从市场中和期权合约中获取，只有波动率是未知的参数。因此，将期权市场价格以及除波动率之外的5个参数代入期权定价公式后推导出波动率，由此计算出的波动率称为隐含波动率。因此，隐含波动率是根据期权的市场价格反推出来的波动率，也是市场自身对未来波动率的预期值，是市场参与者通过期权交易的买卖报价对未来波动率达成的共识。

> 上述说明来自蒋祥林《金融衍生工具》

数据的话，直接使用[tushare](https://tushare.pro/)即可，这是一个提供各类金融数据的API，我以前也做过介绍，见：[https://www.sitstars.com/archives/37/](https://www.sitstars.com/archives/37/)

| 参数 | tushare函数 | tushare参数 | 说明 |
| – | – | – | – |
| 执行价格（K） | opt_basic | exercise_price | |
| 到期时间（t） | opt_basic | delist_date | 需要进行处理 |
| 标的资产价格（S） | fund_daily | close | 上证50ETF代码为510050.SH |
| 利率（rf） | shibor | 1w | 其实用7日银行间国债回购利率更好，但是tushare不提供，只能用一周shibor利率代替 |
| 期权最新价格（c） | opt_daily | settle | |

分红率用看涨看跌期权平价关系求出：

$$
cleft(tright)+Ke^{-rleft(T-tright)}=pleft(tright)+e^{-qleft(T-tright)}Sleft(tright)

$$

将上述参数代入Black-Scholes定价公式，即可求出相应隐含波动率

$$
d_1=frac{ln{left(Sleft(tright)/Kright)}+left(r-q+sigma^2/2right)left(T-tright)}{sigmasqrt{T-t}}

$$

$$
d_2=frac{ln{left(Sleft(tright)/Kright)}+left(r-q-sigma^2/2right)left(T-tright)}{sigmasqrt{T-t}}=d_1-sigmasqrt{T-t}

$$

看涨期权：

$$
cleft(tright)=Sleft(tright)e^{-qleft(T-tright)}Nleft(d_1right)-Ke^{-rleft(T-tright)}Nleft(d_2right)

$$

看跌期权：

$$
p=Ke^{-rleft(T-tright)}Nleft(-d_2right)-Se^{-qleft(T-tright)}Nleft(-d_1right)

$$

## 代码实现

### 导入库

“`python
import pandas as pd
import datetime
import tushare as ts
import numpy as np
from math import log,sqrt,exp
from scipy import stats
import plotly.graph_objects as go
import plotly
import plotly.express as px

pro = ts.pro_api()
plotly.offline.init_notebook_mode(connected=True)
“`

### 提取数据并清洗

“`python
def extra_data(date): # 提取数据
# 提取50ETF合约基础信息
df_basic = pro.opt_basic(exchange=’SSE’, fields=’ts_code,name,call_put,exercise_price,list_date,delist_date’)
df_basic = df_basic.loc[df_basic[‘name’].str.contains(’50ETF’)]
df_basic = df_basic[(df_basic.list_date<=date)&(df_basic.delist_date>date)] # 提取当天市场上交易的期权合约
df_basic = df_basic.drop([‘name’,’list_date’],axis=1)
df_basic[‘date’] = date

# 提取日线行情数据
df_cal = pro.trade_cal(exchange=’SSE’, cal_date=date, fields = ‘cal_date,is_open,pretrade_date’)
if df_cal.iloc[0, 1] == 0:
date = df_cal.iloc[0, 2] # 判断当天是否为交易日，若否则选择前一个交易日

opt_list = df_basic[‘ts_code’].tolist() # 获取50ETF期权合约列表
df_daily = pro.opt_daily(trade_date=date,exchange = ‘SSE’,fields=’ts_code,trade_date,settle’)
df_daily = df_daily[df_daily[‘ts_code’].isin(opt_list)]

# 提取50etf指数数据
df_50etf = pro.fund_daily(ts_code=’510050.SH’, trade_date = date,fields = ‘close’)
s = df_50etf.iloc[0, 0]

# 提取无风险利率数据（用一周shibor利率表示）
df_shibor = pro.shibor(date = date,fields = ‘1w’)
rf = df_shibor.iloc[0,0]/100

# 数据合并
df = pd.merge(df_basic,df_daily,how=’left’,on=[‘ts_code’])
df[‘s’] = s
df[‘r’] = rf
df = df.rename(columns={‘exercise_price’:’k’, ‘settle’:’c’})
#print(df)
return df

def data_clear(df): # 数据清洗
def days(df): # 计算期权到期时间
start_date = datetime.datetime.strptime(df.date,”%Y%m%d”)
end_date = datetime.datetime.strptime(df.delist_date,”%Y%m%d”)
delta = end_date – start_date
return int(delta.days)/365

def iq(df): # 计算隐含分红率
#q = -log((df.settle+df.exercise_price*exp(-df.interest*df.delta)-df.settle_p)/(df.s0))/df.delta
q = -log((df.c+df.k*exp(-df.r*df.t)-df.c_p)/(df.s))/df.t
return q

df[‘t’] = df.apply(days,axis = 1)
df = df.drop([‘delist_date’,’date’],axis = 1)

# 计算隐含分红率
df_c = df[df[‘call_put’]==’C’]
df_p = df[df[‘call_put’]==’P’]
df_p = df_p.rename(columns={‘c’:’c_p’,’ts_code’:’ts_code_p’,
‘call_put’:’call_put_p’})
df = pd.merge(df_c,df_p,how=’left’,on=[‘trade_date’,’k’,’t’,’r’,’s’])

df[‘q’] = df.apply(iq,axis = 1)
c_list = [x for x in range(8)]+[11]

df_c = df.iloc[:,c_list]
df_p = df[[‘ts_code_p’,’trade_date’,’c_p’,’call_put_p’,’k’,’t’,’r’,’s’,’q’]]
df_p = df_p.rename(columns={‘c_p’:’c’,’ts_code_p’:’ts_code’,
‘call_put_p’:’call_put’})
df_c = df_c.append(df_p)
#print(df_c)
return df_c
“`

### 计算隐含波动率

因为公式中有正态分布存在，所以我们没办法直接求出解，只能用二分法求出数值解。

“`python
#根据BS公式计算期权价值
def bsm_value(s,k,t,r,sigma,q,option_type):
d1 = ( log( s/k ) + ( r -q + 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
d2 = ( log( s/k ) + ( r -q – 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
if option_type.lower() == ‘c’:
value = (s*exp(-q*t)*stats.norm.cdf( d1) – k*exp( -r*t )*stats.norm.cdf( d2))
else:
value = k * exp(-r * t) * stats.norm.cdf(-d2) – s*exp(-q*t) * stats.norm.cdf(-d1)
return value

#二分法求隐含波动率
def bsm_imp_vol(s,k,t,r,c,q,option_type):
c_est = 0 # 期权价格估计值
top = 1 #波动率上限
floor = 0 #波动率下限
sigma = ( floor + top )/2 #波动率初始值
count = 0 # 计数器
while abs( c – c_est ) > 0.000001:
c_est = bsm_value(s,k,t,r,sigma,q,option_type)
#根据价格判断波动率是被低估还是高估，并对波动率做修正
count += 1
if count > 100: # 时间价值为0的期权是算不出隐含波动率的，因此迭代到一定次数就不再迭代了
sigma = 0
break

if c – c_est > 0: #f(x)>0
floor = sigma
sigma = ( sigma + top )/2
else:
top = sigma
sigma = ( sigma + floor )/2
return sigma

def cal_iv(df): # 计算主程序
option_list = df.ts_code.tolist()

df = df.set_index(‘ts_code’)
alist = []

for option in option_list:
s = df.loc[option,’s’]
k = df.loc[option,’k’]
t = df.loc[option,’t’]
r = df.loc[option,’r’]
c = df.loc[option,’c’]
q = df.loc[option,’q’]
option_type = df.loc[option,’call_put’]
sigma = bsm_imp_vol(s,k,t,r,c,q,option_type)
alist.append(sigma)
df[‘iv’] = alist
return df
“`

### 多项式拟合填充缺失值

由于部分期限的50etf期权品种较少，所以使用多项式拟合的方法补全这部分数据，由波动率微笑可知，这是一个二次函数，因此令`degree = 2`

此外，用plotly做3D曲面图时，需要数据格式如下：


因此，用pandas的`pivot_table`函数进行整理。

“`python
def data_pivot(df): # 数据透视
df = df.reset_index()
option_type = ‘C’ # 具有相同执行价格、相同剩余到期时间的看涨看跌期权隐含波动率相等，因此算一个就够了
df = df[df[‘call_put’]==option_type]
df = df.drop([‘ts_code’,’trade_date’,’c’,’s’,’r’,’call_put’,’q’],axis=1)
df[‘t’] = df[‘t’]*365
df[‘t’] = df[‘t’].astype(int)
df = df.pivot_table(index=[“k”],columns=[“t”],values=[“iv”])
df.columns = df.columns.droplevel(0)
df.index.name = None
df = df.reset_index()
df = df.rename(columns={‘index’:’k’})

return df

def fitting(df): # 多项式拟合
col_list = df.columns
for i in range(df.shape[1]-1):
x_col = col_list[0]
y_col = col_list[i+1]
df1 = df.dropna(subset=[y_col])

x = df1.iloc[:,0]
y = df1.iloc[:,i+1]

degree = 2

weights = np.polyfit(x, y, degree)
model = np.poly1d(weights)
predict = np.poly1d(model)
x_given_list = df[pd.isnull(df[y_col]) == True][x_col].tolist()
# 所有空值对应的k组成列表
for x_given in x_given_list:
y_predict = predict(x_given)
df.loc[df[x_col]==x_given, y_col] = y_predict
return df
“`

### 作图

“`python
def im_surface(): # 波动率曲面作图
df = plot_df()
df = fitting(df)
#df.to_excel(‘iv_fitting.xlsx’)
df = df.set_index(‘k’)

y = np.array(df.index)
x = np.array(df.columns)
fig = go.Figure(data=[go.Surface(z=df.values, x=x, y=y)])

fig.update_layout(scene = dict(
xaxis_title=’剩余期限’,
yaxis_title=’执行价格’,
zaxis_title=’隐含波动率’),
width=1400,
margin=dict(r=20, b=10, l=10, t=10))
#fig.write_image(“fig1.jpg”)
plotly.offline.plot(fig)

def smile_plot(): # 波动率微笑作图
df = plot_df()
df = df.set_index(‘k’)
df = df.stack().reset_index()
df.columns = [‘k’, ‘days’, ‘iv’]
fig = px.line(df, x=”k”, y=”iv”, color=”days”,line_shape=”spline”)
plotly.offline.plot(fig)
“`

### 综合

“`python
def main():
date = ‘20200515’ # 可任意更改日期
df = extra_data(date) # 提取数据
df = data_clear(df) # 数据清洗
df = cal_iv(df) # 计算隐含波动率
data_pivot(df) # 数据透视表
smile_plot() # 波动率微笑
im_surface() # 波动率曲面

if __name__ == ‘__main__’:
main()
“`

## 效果


Read More