实际波动率与隐含波动率的计算 python

计算第i日 后，n(20)日的实际波动率，以及期权的隐含波动率

# 实际波动率

import math
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd


def cal_hist_vol():
    df = pd.read_csv('FUND_MKT_Quotation.csv')
    df = df.drop(
        columns=['OpenPrice', 'HighPrice', 'LowPrice', 'Volume', 'Amount', 'AvgPrice', 'ChangeRatio', 'TurnoverRate',
                 'Covert', 'CovertRate'])
    df['ratio'] = df['ClosePrice'] / df['PreClosePrice']
    df['log_return'] = np.log(df['ratio'])
    df['avg_20_log_return'] = (df['log_return'] + df['log_return'].shift(-1) + df['log_return'].shift(-2) + df[
        'log_return'].shift(-3) + df['log_return'].shift(-4) + df['log_return'].shift(-5) + df['log_return'].shift(-6) +
                               df['log_return'].shift(-7) + df['log_return'].shift(-8) + df['log_return'].shift(-9)
                               + df['log_return'].shift(-10) + df['log_return'].shift(-11) + df['log_return'].shift(
                -12) + df['log_return'].shift(-13) + df['log_return'].shift(-14) + df['log_return'].shift(-15) + df[
                                   'log_return'].shift(-16) + df['log_return'].shift(-17) + df['log_return'].shift(
                -18) + df['log_return'].shift(-19)) / 20
    df['diff_r_i-r_bar'] = np.power((df['log_return'] - df['avg_20_log_return']), 2)
    df['sum_diff_power'] = df['diff_r_i-r_bar'] + df['diff_r_i-r_bar'].shift(-1) + df['diff_r_i-r_bar'].shift(-2) + df[
        'diff_r_i-r_bar'].shift(-3) + df['diff_r_i-r_bar'].shift(-4) + df['diff_r_i-r_bar'].shift(-5) + df[
                               'diff_r_i-r_bar'].shift(-6) + df['diff_r_i-r_bar'].shift(-7) + df['diff_r_i-r_bar'].shift(-8) + df['diff_r_i-r_bar'].shift(-9)
    + df['diff_r_i-r_bar'].shift(-10) + df['diff_r_i-r_bar'].shift(-11) + df['diff_r_i-r_bar'].shift(-12) + df['diff_r_i-r_bar'].shift(-13) + df['diff_r_i-r_bar'].shift(-14)
    + df['diff_r_i-r_bar'].shift(-15) + df['diff_r_i-r_bar'].shift(-16) + df['diff_r_i-r_bar'].shift(-17) + df['diff_r_i-r_bar'].shift(-18) + df['diff_r_i-r_bar'].shift(-19)

    df['RV'] = np.sqrt(242 * df['sum_diff_power'] / 19)

    df.to_csv('RV')

if __name__ == '__main__':
    cal_hist_vol()

# 隐含波动率 ——牛顿法、二分法

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


# 根据公式计算期权价值
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


# 求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))
    if sigma_est >= 0.5:
        sigma_est =None
    return sigma_est


# 二分法求隐含波动率
# 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
#

# 读取行情数据

call_300_18 = pd.read_csv(r'C:\Users\heluxiaofu\PycharmProjects\pythonProject\论文\波动率微笑\2.3_02_call.csv',
                          encoding='utf-8')
# put_300_18 = pd.read_csv(r'C:\Users\heluxiaofu\PycharmProjects\pythonProject\论文\波动率微笑\2.3_02_put.csv', encoding='utf-8')

call_300_17 = pd.read_csv(r'C:\Users\heluxiaofu\PycharmProjects\pythonProject\论文\波动率微笑\2.4_02_call.csv')
# 取期权价格和行权价
price_18 = call_300_18['ClosePrice']
k_18 = call_300_18['StrikePrice']
price_17 = call_300_17['ClosePrice']
k_17 = call_300_17['StrikePrice']
# 到期时间年化
t_18 = 0.0493  # 18天
t_17 = 0.0466
# 标的初始价格
s0_18 = 3688.358
s0_17 = 3785.638
# 无风险利率，用shibor近似
rf = 0.015

# 用两种方法分别求三只期权的隐含波动率，并打印
sigma_init = 1
sigma_18_newton = []
for i in range(call_300_18.shape[0]):
    sigma_18_newton.append(bsm_call_imp_vol_newton(s0_18, k_18[i], t_18, rf, price_18[i], sigma_init))

sigma_17_newton = []
for i in range(call_300_17.shape[0]):
    sigma_17_newton.append(bsm_call_imp_vol_newton(s0_17, k_17[i], t_17, rf, price_17[i], sigma_init))

print('imp_vol_newton_18:')
print(sigma_18_newton)
print('imp_vol_newton_17:')
print(sigma_17_newton)
# 绘制牛顿法求解的波动率曲线
plt.plot(k_18, sigma_18_newton, label='18toMaturity', lw=1.5, )
plt.plot(k_18, sigma_18_newton, 'r.')
plt.plot(k_17, sigma_17_newton, label='17toMaturity', lw=2, )
plt.plot(k_17, sigma_17_newton, 'r.')

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