【BSM模型】用實際市場數據計算隱含波動率並驗證波動率微笑

在Black-Scholes期權定價模型中,不能直接觀察到的參數只有股票價格的波動率。波動率可以由歷史數據進行估計,這是歷史波動率。隱含波動率也是交易員非常關心的,隱含波動率是期權的市場價格中所包含的波動率,即由期權價格和期權定價公式反推的波動率。隱含波動率和歷史波動率作比較,可以指導投資者的操作。投資者可以直接買賣波動率,或者參考波動率確定買賣時機。

我們可以通過期權定價公式寫出隱含波動率的方程,但是直接解方程非常困難,因為這個方程不存在閉合解。既然是用程序求解,當然可以用計算機求方程解的神器-數值計算。牛頓迭代法和二分法是求隱含波動率常用的兩個方法。相比二分法,牛頓迭代法是更通用的近似求解方程的方法。

由於國內沒有場內個股期權,曲曲菜用上證50ETF期權做分析。首先從新浪財經的網站獲得期權的行情信息,並存入csv文件。到期時間我選了16天,51天,76天三種,分別存成三個文件。

新浪財經的期權行情數據(16天到期)

期權數據文件(16天到期)

然後就可以計算隱含波動率了,計算隱含波動率的python程序如下。

一.BSM模型

1.引入所用到的庫

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from math import log,sqrt,exp
from scipy import stats

2. 定價公式的程序實現

#根據公式計算期權價值
def bsm_call_value(s0,k,t,r,sigma):
d1 = ( log( s0/k ) + ( r + 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
d2 = ( log( s0/k ) + ( r - 0.5*sigma**2 )*t )/( sigma*sqrt(t) )
value = ( s0*stats.norm.cdf( d1,0.,1. ) - k*exp( -r*t )*stats.norm.cdf( d2,0.,1 ))
#print(cvalue,value)
return value

二.牛頓迭代法

1.介紹

設r是f(x)=0的根,選取x0作為r的初始近似值,過點(x0,f(x0))做曲線y=f(x)

的切線L ,

,則L與x軸交點的橫坐標

,稱x1為r的一次近似值。過點(x1,f(x1))做曲線y=f(x)的切線,並求該切線與x軸交點的橫坐標

,稱x2為r的二次近似值。重複以上程,得r的近似值序列,其中,

稱為r的n+1次近似值,上式稱為牛頓迭代公式。

隱含波動率的計算中,f(x)是BSM定價公式得到的價格和實際價格的差,x是隱含波動率,f(x)中只有隱含波動率是未知數,其他是已知數。由定義可知,x的導數是vega。

2.程序實現

#求vega
def bsm_vega(s0,k,t,r,sigma):
d1 = log( s0/k ) + ( r + 0.5*sigma**2 )*t/( sigma*sqrt(t) )
vega = s0*stats.norm.cdf(d1,0.,1.)*sqrt(t)
#print(vega,vega)
return vega

#牛頓迭代法求隱含波動率,迭代次數設為100
def bsm_call_imp_vol_newton(s0, k, t, r, c0, sigma_est, it = 100):
for i in range(it):
sigma_est -= ((bsm_call_value(s0, k, t, r, sigma_est) - c0)/
bsm_vega(s0, k, t, r, sigma_est))
return sigma_est

三.二分法

1.介紹

二分法求根的思想和二分查找相同,只不過二分查找比較的是未知數和目標,而二分法比較的是未知數的函數和目標函數值。二分法求根需要知道根所在的區間,將上下限分別設為區間的上下邊界點,初始值設為上下限的均值。通過迭代不斷更新並逼近方程的解。

2.程序實現

#二分法求隱含波動率
def bsm_call_imp_vol_dichotomy(s0,k,t,r,c):
c_est = 0
top = 3 #波動率上限
floor = 0 #波動率下限
sigma = ( floor + top )/2 #波動率初始值

while abs( c - c_est ) > 1e-8:
c_est = bsm_call_value(s0,k,t,r,sigma)
#根據價格判斷波動率是被低估還是高估,並對波動率做修正
if c - c_est > 0: #f(x)>0
floor = sigma
sigma = ( sigma + top )/2
else:
top = sigma
sigma = ( sigma + floor )/2
return sigma

四.讀取行情數據並初始化參數

1.程序實現

#讀取行情數據
call_50etf_16 = pd.read_csv(D:\data\options\call_50etf_16.csv,encoding=utf-8)
call_50etf_51 = pd.read_csv(D:\data\options\call_50etf_51.csv,encoding=utf-8)
call_50etf_79 = pd.read_csv(D:\data\options\call_50etf_79.csv,encoding=utf-8)

#取期權價格和行權價
price_16 = call_50etf_16[最新價]
k_16 = call_50etf_16[行權價]
price_51 = call_50etf_51[最新價]
k_51 = call_50etf_51[行權價]
price_79 = call_50etf_79[最新價]
k_79 = call_50etf_79[行權價]

#到期時間年化
t_16 = 16/365
t_51 = 51/365
t_79 = 79/365

#標的初始價格
s0 = 2.533
#無風險利率,用shibor近似
rf = 0.025

五.計算隱含波動率

1.程序實現

#用兩種方法分別求三隻期權的隱含波動率,並列印
sigma_init=1
sigma_16_newton=[]
sigma_16_dichotomy=[]
for i in range(call_50etf_16.shape[0]):
sigma_16_newton.append(bsm_call_imp_vol_newton(s0,k_16[i],t_16,rf,price_16[i],sigma_init))
sigma_16_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_16[i],t_16,rf,price_16[i]))

print(imp_vol_newton_16:)
print(sigma_16_newton)
print(imp_vol_dichotomy_16:)
print(sigma_16_dichotomy)
print(-------------------------------------------------------------------------------------------------------------------------------------)
#
sigma_51_newton=[]
sigma_51_dichotomy=[]
for i in range(call_50etf_51.shape[0]):
sigma_51_newton.append(bsm_call_imp_vol_newton(s0,k_51[i],t_51,rf,price_51[i],sigma_init))
sigma_51_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_51[i],t_51,rf,price_51[i]))

print(imp_vol_newton_51:)
print(sigma_51_newton)
print(imp_vol_dichotomy_51:)
print(sigma_51_dichotomy)
print(-------------------------------------------------------------------------------------------------------------------------------------)

#
sigma_79_newton=[]
sigma_79_dichotomy=[]
for i in range(call_50etf_79.shape[0]):
sigma_79_newton.append(bsm_call_imp_vol_newton(s0,k_79[i],t_79,rf,price_79[i],sigma_init))
sigma_79_dichotomy.append(bsm_call_imp_vol_dichotomy(s0,k_79[i],t_79,rf,price_79[i]))

print(imp_vol_newton_79:)
print(sigma_79_newton)
print(imp_vol_dichotomy_79:)
print(sigma_79_dichotomy)
print(-------------------------------------------------------------------------------------------------------------------------------------)

2.隱含波動率列印結果

imp_vol_newton_16:
[0.3638689603157307, 0.3647065048442553, 0.33796965226626546, 0.3087508280418194, 0.29010352112059345, 0.27704783086782625, 0.2678504434116786, 0.27958298121452063, 0.2815159982582914, 0.2892477867510613, 0.2948778490238011, 0.3001603015476018, 0.3079952199130588]
imp_vol_dichotomy_16:
[0.3638690114021301, 0.3647065758705139, 0.33796969056129456, 0.3087505102157593, 0.290103480219841, 0.2770475149154663, 0.2678511142730713, 0.2795829623937607, 0.28151603043079376, 0.2892477214336395, 0.29487812519073486, 0.3001604676246643, 0.30799537897109985]
-------------------------------------------------------------------------------------------------------------------------------------
imp_vol_newton_51:
[0.2663830055838, 0.2623349770104388, 0.25886692974504383, 0.25514765329363565, 0.25388257585478174, 0.25112748137747243, 0.25311127192458244, 0.25282431711916714, 0.25334527303228827]
imp_vol_dichotomy_51:
[0.26638298481702805, 0.26233501732349396, 0.25886694341897964, 0.2551477253437042, 0.25388259440660477, 0.2511274963617325, 0.25311121344566345, 0.2528250217437744, 0.25334523618221283]
-------------------------------------------------------------------------------------------------------------------------------------
imp_vol_newton_79:
[0.30525607059637466, 0.2900199629942286, 0.29426758619557364, 0.26917014044939297, 0.2708636644397624, 0.25993171220348904, 0.2610946078237793, 0.2527101923035752, 0.2523987108427234, 0.24752459839396573, 0.2483057133362068, 0.24353321749378323, 0.24308705976680656, 0.23851008170784155, 0.2406859629151288, 0.24016117575185575]
imp_vol_dichotomy_79:
[0.30525608360767365, 0.29002001881599426, 0.29426760971546173, 0.2691701799631119, 0.2708636373281479, 0.25993166863918304, 0.261094618588686, 0.2527102008461952, 0.25239837169647217, 0.2475244402885437, 0.24830570071935654, 0.2435331791639328, 0.243087038397789, 0.2385101616382599, 0.24068592488765717, 0.24016119539737701]
-------------------------------------------------------------------------------------------------------------------------------------

六.繪製波動率曲線

1.繪製牛頓法曲線的程序實現

#繪製牛頓法求解的波動率曲線
plt.plot(k_16,sigma_16_newton,label=16-days,lw=1.5,)
plt.plot(k_16,sigma_16_newton,r.)
plt.plot(k_51,sigma_51_newton,label=51-days,lw=1.5)
plt.plot(k_51,sigma_51_newton,g.)
plt.plot(k_79,sigma_79_newton,label=79-days,lw=1.5)
plt.plot(k_79,sigma_79_newton,k.)

plt.grid(True)
plt.xlabel(Strike price)
plt.ylabel(Implied volatility)
plt.legend()
plt.show()

2.上一步繪製出的圖形

3.繪製二分法曲線的程序實現

#繪製二分法求解的波動率曲線
plt.plot(k_16,sigma_16_dichotomy,label=16-days,lw=1.5,)
plt.plot(k_16,sigma_16_dichotomy,r.)
plt.plot(k_51,sigma_51_dichotomy,label=51-days,lw=1.5)
plt.plot(k_51,sigma_51_dichotomy,g.)
plt.plot(k_79,sigma_79_dichotomy,label=79-days,lw=1.5)
plt.plot(k_79,sigma_79_dichotomy,k.)

plt.grid(True)
plt.xlabel(Strike price)
plt.ylabel(Implied volatility)
plt.legend()
plt.show()

4.上一步繪製出的圖形

從圖形可以看出,51天到期和79天到期的隱含波動率隨執行價格的遞增,呈現遞減趨勢,這就是股票期權的波動率微笑(volatility smile)。16天到期的隱含波動率是隨執行價格遞增先是遞減,至標的價格附近後,開始緩慢遞增,這也是波動率微笑,雖然對股票期權來說,這個微笑不是很標準。從圖形還可以看出,距離到期時間越近,隱含波動率越大。

波動率微笑反映了隱含波動率和執行價格的關係。外匯期權的波動率是對稱的微笑, 股票期權的波動率微笑是不對稱的,更準確的叫法是volatility skew(波動率傾斜),或者volatility smirk(波動率假笑)。反映到圖形上,就是左高有低,隱含波動率隨執行價格的遞增而遞減。

股票期權波動率微笑的原因,常見的解釋是槓桿效應和恐慌情緒,但是也有人認為這就是一個市場的反應,沒有特別的原因。(如需要本文的源代碼和數據文件可以聯繫我)

代碼在我的GitHub:

ququcai/volatility_smile?

github.com圖標

參考資料

[1] 約翰 赫爾.期權、期貨及其他衍生品

[2] Yves Hilpsch. Python for Finance: Analyze Big Financial Data

本文作者:曲曲菜(微信公眾號:曲曲菜)

知乎專欄:AI和金融模型

原創作品,未標明作者不得轉載。

推薦閱讀:

TAG:隱含波動率 | BlackScholesMerton模型 | 上證指數 |