雪球产品期权价值蒙特卡洛模拟（2）

先导入所需的模块

import os
import sys
import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
from tqdm import tqdm
import math
import scipy.stats as stats
from pyecharts.globals import CurrentConfig, NotebookType
CurrentConfig.NOTEBOOK_TYPE = NotebookType.JUPYTER_LAB
import pyecharts.options as opts
from pyecharts.charts import Scatter
import pyecharts.options as opts
from pyecharts.charts import Bar, Grid, Line, Liquid, Page, Pie, Timeline
from pyecharts.faker import Faker
from pyecharts import options as opts
from pyecharts.commons.utils import JsCode
from pyecharts.components import Table
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import time
import tushare as ts
import datetime as dt
pro = ts.pro_api('tushare接口token')

接着从tushare读取并整理沪深300指数的所需数据

index_dict = {'000300.SH':'沪深300' }
for i in index_dict:
    df = pro.index_daily(ts_code=i)
    df['trade_date'] = pd.to_datetime(df['trade_date'],format='%Y-%m-%d')

df = df.sort_values('trade_date')  #df按日期从前往后排序
df['trade_date'] = df['trade_date'].astype(str)
df['year'] = [i[0:4] for i in df['trade_date']]  #增加年份列
df

数据如下

基于沪深300指数数据，计算年化收益率及年化波动率：

days=(pd.to_datetime(df.trade_date[0])-pd.to_datetime(df.trade_date[4422])+dt.timedelta(1)).days
years=days/365.25
#totaltimes=df.close[0]/df.close[4422]*(df.pct_chg[4422]/100+1)
#annualreturngeo = math.pow(totaltimes,1/years)-1   #计算hs300指数的年化收益率（几何收益率）
annualvol=df.pct_chg.std()*math.sqrt(252)/100  #计算hs300指数的平均波动率
annualreturn=df.pct_chg.sum()/years/100
round(years, 2), round(annualreturn, 2), round(annualvol, 2)

 结果如下

years=18.2, annualreturn=0.11, annualvol=0.26

根据Black-Scholes期权定价模型定义Simulation_StockPrice()函数，用于计算从S0到S1。并通过循环模拟10000次从期初到期末的标的价格走势。结果得到一个10000*252大小的价格矩阵price_matrix。

def Simulation_StockPrice(S0, mu, Std, dt):
    S1 = S0 * np.exp((mu - 0.5 * Std**2)*dt + Std*np.sqrt(dt)*np.random.standard_normal())
    return S1    #根据S0计算S1
price_matrix = []   #先定义一个空数组
mu = 0.14132335622963155
Std = 0.26935858659655915
dt = 1/252. 
S0 = 5255.29

for i in range(10000):     #执行10000次模拟
    price_simulator = []   #定义一个空数组，或者将上一次循环的数组结果清零：这个数组放的是一次模拟的结果，即未来252天的股价表现
    S1 = S0
    for j in range(252):   #执行一个252天的循环，循环内每执行一步，就会模拟出第二天的股价
        S1 = Simulation_StockPrice(S1, mu, Std, dt)   #模拟第二天的股价，并赋值给S0
        price_simulator.append(S1)       #将S1放入数组里面，循环执行完毕后，会包括252个数，即未来252天的股价
    price_matrix.append(price_simulator)  #将上一个循环结果（即模拟的一次252天的路径放到这个新的数组里，如果10000条模拟结束之后，该数组应该是一个10000列*252行的二维数组）
print(price_matrix)

根据计算结果绘制10000次蒙特卡洛模拟的标的价格曲线如下

按照上一篇文章中列举的流行的一款雪球产品进行价值分析

o 敲出 – 100%本金+持有期内敲出票息

o 敲入但未敲出 – Min（期末价/期初价，100%）

o 未敲出也未敲入 – 100%本金+持有期内红利票息

 先针对其中一次（price_matrix[0]，第一次）价格路径进行分析：

count = 1
a = [price_matrix[0][j] for j in  [i*21-1 for i in range(1,13)]]
knock_out=0
for i in range(len(a)):
    if a[i]> price_matrix[0][0]*1.05:
        knock_out = a.index(a[i])
        break
    else:
        knock_out = 0
if (knock_out != 0) :  
    pv = ((knock_out  / 12) * 0.2 * price_matrix[0][0]) / (1 + 0.015*(knock_out  / 12))
elif  all(i >= price_matrix[0][0]*0.8 for i in price_matrix[0]):
    pv = 0.2 * price_matrix[0][0]
else:
    if (price_matrix[0][-1] > price_matrix[0][0] ):
        pv = 0
    else: 
        pv = (price_matrix[0][-1] - price_matrix[0][0]) * 0.9852
knock_out, pv

结果是敲出的：

(1, 88.42933303060764)

类似的，对这10000次模拟的价格路径都进行测试：

test_pv_list = []
for count in tqdm(range(10000)):
    a = [price_metrix[count][j] for j in  [i*21-1 for i in range(1,13)]]
    knock_out=0
    for i in range(len(a)):
        if a[i]> price_metrix[count][0]*1.05:
            knock_out = a.index(a[i])
            break
        else:
            knock_out = 0
    if (knock_out != 0) :  
        pv = ((knock_out  / 12) * 0.2 * price_metrix[count][0]) / (1 + 0.015*(knock_out  / 12))
    elif  all(i >= price_metrix[count][0]*0.8 for i in price_metrix[count]):
        pv = 0.2 * price_metrix[count][0]
    else:
        if (price_metrix[count][-1] > price_metrix[count][0] ):
            pv = 0
        else: 
            pv = (price_metrix[count][-1] - price_metrix[count][0]) * 0.9852
    test_pv_list.append([knock_out, pv])

pv_list_df = pd.DataFrame(test_pv_list)
pv_list_df.columns = ['敲出日','PV']
pv_list_df['敲出日'] = pv_list_df['敲出日'] * 21

定义期权价值函数并根据沪深300雪球产品参数计算看跌期权的价值

def option_value(S,K,r,T,sigma,q,option_type):
    d1 = ( np.log(S/K) + ( r-q+0.5*sigma**2 )*T )/( sigma*math.sqrt(T) )
    d2 = ( np.log(S/K) + ( r-q-0.5*sigma**2 )*T )/( sigma*math.sqrt(T) )
    if option_type == 'call':
        value = (math.exp(-q*T)*S*stats.norm.cdf(d1) - K*math.exp(-r*T)*stats.norm.cdf(d2))
    elif option_type == 'put':
        value = math.exp(-r*T)*K*stats.norm.cdf(-d2) - S*math.exp(-q*T)*stats.norm.cdf(-d1)
    else:
        value = np.nan
    return value

option_value(5255.29,5255.29,0.015,0.2166,0.26935858659655917,0,'put')

结果为：

253.7916790717095

绘制标的资产波动率和雪球产品收益率的关系曲线

line = (
        Line()
        .add_xaxis([0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6])
        .add_yaxis("收益率", [10.12, 10.2, 9.27, 6.55, 4.05, 2.09, 1.22, -0.17, -1.41, -2.41, -3.75, -4.36])
        .set_global_opts(
     title_opts=opts.TitleOpts("标的资产波动率与雪球产品收益率"),
     legend_opts=opts.LegendOpts(pos_top="5%"),
            
        )
 )
line.render('1.html')

曲线如下

 继续绘制标的资产预期收益率和雪球产品收益率的关系曲线

x2  = [ str(i) for i in [0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4]]
y2 =   [-1.0406,1.2914,2.5072, 3.4441,4.4952,4.9936,5.3288,5.5710,5.6041]

c2 = (
    Line()
    .add_xaxis(x2)
    .add_yaxis("", y2, is_smooth=True )
    .set_global_opts(
        title_opts=opts.TitleOpts(title="标的资产预期收益与雪球产品收益率"),
        yaxis_opts=opts.AxisOpts(
            name='(%)',
            type_="value",
            axistick_opts=opts.AxisTickOpts(is_show=True),
            splitline_opts=opts.SplitLineOpts(is_show=True),
        ),
    )
)
c2.render('1.html')

 曲线如下

 计算雪球产品的期望收益率（irr），模拟10000次

n=10000  #模拟n次

def Simulation_StockPrice(S0, mu, Std, dt):
    S1 = S0 * np.exp((mu - 0.5 *Std**2)*dt + Std*np.sqrt(dt)*np.random.standard_normal())
    return S1

price_metrix = []
mu = 0.14132335622963155  #zz：HS300期望收益率（年化）
Std = 0.26935858659655915 #zz：HS300波动率（年化）
dt = 1/252. 
S0 = 5255.29  #HS指数期初值
for i in tqdm(range(n)):
    S0 = 5255.29
    price_simulator = [S0]
    for j in range(252):
        S0 = Simulation_StockPrice(S0, mu, Std, dt)
        price_simulator.append(S0)      
    price_metrix.append(price_simulator)

count = 1
a = [price_metrix[0][j] for j in  [i*21 for i in range(1,13)]]  #观察日HS300指数值（每月末）
knock_out=0

for i in range(len(a)):
    if a[i]> price_metrix[0][0]*1.05:
        knock_out = a.index(a[i])
        break
    else:
        knock_out = 0
#Case1_计算情形1&2发生情况下“收益”的pv（考虑贴现）
if (knock_out != 0) :  
    pv = ((knock_out  / 12) * 0.2 * price_metrix[0][0]) / (1 + 0.015*(knock_out  / 12))
#Case2_计算情形3发生情况下“收益”的pv（考虑贴现）
elif  all(i >= price_metrix[0][0]*0.8 for i in price_metrix[0]):
    pv = 0.2 * price_metrix[0][0] / (1 + 0.015)  #zz：加上贴现（1年）
#Case3_计算情形4&5发生情况下“收益”的pv（考虑贴现）
else:
    if (price_metrix[0][-1] > price_metrix[0][0] ):  #期末价格高于初始值
        pv = 0                                       #收益为0
    else: 
        pv = (price_metrix[0][-1] - price_metrix[0][0]) / (1 + 0.015)  #期末价格低于初始值，按照实际收益兑付（考虑1年贴现）

test_pv_list = []
for count in tqdm(range(n)):
    a = [price_metrix[count][j] for j in  [i*21 for i in range(1,13)]]  #a为观察日价格列表
    knock_out=0
    for i in range(len(a)):
        if a[i]> price_metrix[count][0]*1.05:
            knock_out = a.index(a[i])  #knockout为敲出日在a列表中的index，knock_out大于0
            break
        else:
            knock_out = 0  #不敲出则该循环的knock_out等于0
    if (knock_out != 0) :  
        pv = ((knock_out  / 12) * 0.2) / (1 + 0.015*(knock_out  / 12))
    elif  all(i >= price_metrix[count][0]*0.8 for i in price_metrix[count]):
        pv = 0.2/(1+0.015)
    else:
        if (price_metrix[count][-1] > price_metrix[count][0] ):
            pv = 0/(1+0.015)
        else: 
            pv = (price_metrix[count][-1]/5255.29 - 1)/(1+0.015)
    test_pv_list.append([knock_out, pv])
pv_list_df = pd.DataFrame(test_pv_list)
pv_list_df.columns = ['敲出日','PV']
pv_list_df['敲出日'] = pv_list_df['敲出日'] * 21
#pv_list_df.to_excel('12445.xlsx')

print('The annual return of snowball is','%.2f%%'%(pv_list_df['PV'].mean()*1.015 * 100) )       #+0.015，zz：应该不需要加0.015
print('The annual volatility of snowball is','%.2f%%'%(pv_list_df['PV'].std()*1.015 * 100) )    #求平均收益率
print('The sharp ratio of snowball is', round((pv_list_df['PV'].mean()*1.015-0.015)/(pv_list_df['PV'].std()*1.015),2)) #求夏普比率

计算结果如下：

The annual return of snowball is 5.28%
The annual volatility of snowball is 14.75%
The sharp ratio of snowball is 0.26

可见雪球产品的期望收益率并不高，不能被票面收益率所诱导，其波动率14.75%倒是比标的资产要低的，但是标的资产预期收益率也更高。最终看下sharpe ratio，其sharpe ratio只有0.26，可以说很一般，比标的资产沪深300指数也要差一些。