期权专题9：雪球期权（一）（BSM+蒙特卡罗定价）

本文主要是介绍期权专题9：雪球期权（一）（BSM+蒙特卡罗定价），希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

1. 前期介绍

1.1 定价思路

1.2 雪球要素

1.3 标的定价

2. 代码复现

2.1 价格路径

2.1.1 单次路径

2.1.2 多次路径

2.2 雪球定价

2.2.1 加入要素

2.2.2 估值定价

3. 完整代码

雪球期权是近两年来备受青睐的衍生品投资标的之一，几个月前就想写一篇关于它的文章。但是因为其定价模型确实比较复杂，无论在思路和代码复现上都存在一定的难度，直到今日方才鼓起勇气一试。由于自我认知上的局限，分享的内容可能不够准确或者全面，望读者见谅，同时对有疑惑的地方欢迎探讨交流。

1. 前期介绍

1.1 定价思路

雪球期权本质上是障碍期权的组合，关于其一些基础性的介绍此处省略(网上相关介绍太多)。本文选用蒙特卡罗模拟法来实现雪球定价的净值化，定价的核心思路主要是以下三点：

1. 使用蒙特卡罗模拟出雪球标的未来可能出现的路径。

2. 根据雪球要素，计算每一条路径的折现净值。

3. 计算所有路径下折现净值的均值，即作为雪球定价的净值。

同时，还需要格外注意的一点是，雪球定价存在2种状态：一种是发生敲入状态下的，一种是未发生敲入状态下的。

1.2 雪球要素

本文选用经典雪球作为案例，其要素如下：

雪球要素
挂钩标的	中证500指数（399905.SZ）
存续期	12个月
封闭期	3
敲入线	1.03
敲出线	0.75
敲入观察频率	每日
敲出观察频率	每月
票息（年化）	20%

同时，为了简化定价模型的复杂程度，本文做出以下假设：

1. 假设波动率的数值为常数，值为22%（近五年指数年化波动率均值）。

2. 假设无风险利率为2.5%，一年的交易日为252天。

3. 假设不考虑股指升贴水，展期带来的损益。

4. 假设不考虑交易过程中所有的摩擦成本。

1.3 标的定价

假设标的满足对数正态分布的条件，将微分方程引入BSM模型，可以推导出雪球所挂钩标的对应的定价公式。

使用蒙特卡罗的思路，对上述公式进行逐逐步的模拟迭代，可以获得未来标的可能出现的价格路径。

2. 代码复现

为了计算的速度以及整体过程的简洁性，考虑将数据抽象出来，使用矩阵的思路来复现整个过程。考虑到文章的适读性，本章节将诸多步骤进行细化，尽可能将过程解释清楚，读者可根据自身理解，对各个小节的内容进行选择性的详略阅读。

2.1 价格路径

2.1.1 单次路径

此小节为步骤的细化分享，可选择性阅读。

首先设置基本的参数，得到初始的路径矩阵。

import numpy as np# 设置参数，分别是无风险利率，波动率和剩余期限
r = 0.025
vol = 0.22
residual_day = 252
# 设置初始值,需要模拟的路径数量
St = 1
path_num = residual_day
# 创建矩阵,并把初始值设为1，将时间差折算成年
s_path = np.zeros((path_num, 1))
s_path[0] = 1
Tt = 1/252 # 折算成年的时间差，此处以交易日为单位

得到的结果如下：

s_path是一个252行，1列的矩阵。接下来，根据初始值1，迭代矩阵的第二个数据，考虑到随机数可能产生异常的数值，因此需要对迭代的数据进行涨跌停的限制。

# 由初始值迭代第二个数据
N = np.random.standard_normal(1)  # 标准正态分布随机数
ST = s_path[0] * np.exp((r - 0.5 * vol ** 2) * Tt + vol * np.sqrt(Tt) * N)
# 进行涨跌停的限制
max_st = s_path[0] * 1.1
min_st = s_path[0] * 0.9
new_st = np.where(ST < min_st,  min_st, ST)
new_st = np.where(ST > max_st, max_st, new_st)
s_path[1] = new_st

对应得到的结果是：

通过上述代码进行循环迭代，得到单次路径的所有数据：

for i in range(1,path_num):N = np.random.standard_normal(1)  # 标准正态分布随机数ST = s_path[i-1] * np.exp((r - 0.5 * vol ** 2) * Tt + vol * np.sqrt(Tt) * N)# 进行涨跌停的限制max_st = s_path[i-1] * 1.1min_st = s_path[i-1] * 0.9new_st = np.where(ST < min_st, min_st, ST)new_st = np.where(ST > max_st, max_st, new_st)s_path[i] = new_st

当前即可得到对应单次的完整路径，对应的结果为：

将s_path对应得到的数据进行可视化：

# 单次路径可视化
import matplotlib.pyplot as pltplt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
plt.rcParams['axes.unicode_minus'] = False
x = [x+1 for x in range(path_num)]
plt.plot(x,s_path,label='单次路径')
plt.legend()
plt.show()

对应等到的结果为：

2.1.2 多次路径

由于单次的路径存在较大的随机性，因此需要进行多次的模拟，此处暂时设定模拟次数为5000次。如果在单次路径的基础上添加循环继续模拟其他的路径，速度相对会比较慢，因此考虑还是使用矩阵来进行迭代。

首先设立基本参数，得到基础的路径矩阵。

import numpy as np# 设置参数，分别是无风险利率，波动率和剩余期限
r = 0.025
vol = 0.22
residual_day = 252
# 设置初始值,需要模拟的路径数量
St = 1
path_num = residual_day
times = 5000 # 模拟次数
# 创建矩阵,并把初始值设为1，将时间差折算成年
s_path = np.zeros((path_num, times))
s_path[0] = 1
Tt = 1/252 # 折算成年的时间差，此处以交易日为单位

对应得到的结果为：

此处的s_path是一个252行（代表日期），5000列（代表路径数量）。接下来，演示一下第二列价格的生成方式。

N = np.random.standard_normal(times)  # 标准正态分布随机数
ST = s_path[0] * np.exp((r - 0.5 * vol ** 2) * Tt + vol * np.sqrt(Tt) * N)
# 进行涨跌停的限制
max_st = s_path[0] * 1.1
min_st = s_path[0] * 0.9
new_st = np.where(ST < min_st,  min_st, ST)
new_st = np.where(ST > max_st, max_st, new_st)
s_path[1] = new_st

对应得到的解惑为：

到此处，个人感觉基本将路径生成的思路说清楚了。后续不再进行赘述，开始封装一个函数，来实现路径模拟的功能。

import numpy as np
import matplotlib.pyplot as pltdef get_s_path(St, r, vol, residual_day, times):'''模拟标的未来可能出现的价格路径:param St: int,标的初始价格，例如1:param r: int, 无风险利率，例如0.02:param vol: int, 波动率，例如0.02:param residual_day: int, 剩余期限(单位为天)，例如252:param times: int, 模拟次数，例如5000:return:  价格路径对应的矩阵'''path_num = residual_days_path = np.zeros((path_num, times))s_path[0] = StTt = 1 / 252for i in range(1, path_num):N = np.random.standard_normal(times)  # 标准正态分布随机数ST = s_path[i - 1] * np.exp((r - 0.5 * vol ** 2) * Tt + vol * np.sqrt(Tt) * N)# 进行涨跌停的限制max_st = s_path[i - 1] * 1.1min_st = s_path[i - 1] * 0.9new_st = np.where(ST < min_st, min_st, ST)new_st = np.where(ST > max_st, max_st, new_st)s_path[i] = new_streturn s_pathif __name__ == '__main__':r = 0.025vol = 0.22residual_day = 252St = 1times = 5000price_path = get_s_path(St, r, vol, residual_day, times)

对应得到的结果是：

将得到的结果进行可视化：

    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']plt.rcParams['axes.unicode_minus'] = Falsex = [x+1 for x in range(residual_day)]for i in range(times):plt.plot(x,price_path[:,i])plt.title('标的价格模拟路径'+'('+'模拟次数:'+str(times)+')')plt.show()

2.2 雪球定价

2.2.1 加入要素

在获取标的价格的模拟路径后，需要在路径上添加雪球对应的要素结构。此处代码续接上一小节的price_path。

    # 雪球要素knock_out = 1.03  # 敲出价格knock_in = 0.75  # 敲入价格coupon = 0.2  # 年化票息lock_time = 3  # 锁定期,单位为月T = 1  # 期限，折算为年all_time = T * 252  # 期限，折算为天survival_day = all_time - residual_day  # 已存续期限# 利用要素对路径进行判断one_path = price_path[:, 0]# 价格低于敲入价格的交易日对应的序号knock_in_day = np.where(one_path <= knock_in)# 价格高于敲出价格的交易日对应的序号knock_out_day = np.where(one_path >= knock_out)# 上述每个交易日的序号需要加上雪球已存续期限new_out_day = knock_out_day[0] + survival_day# 判断上述日期是否是观察日,1年252个交易日，因此一个月设定为12obs_may_day = new_out_day[new_out_day % 21 == 0]# 由于有锁定期，需要剔除前三个月的敲出观察日obs_day = obs_may_day[obs_may_day / 21 > 3]

此处，需要注意两点：

1. 未使用具体日期来作为观察日的检测，由于数据本身是模拟的，同时模拟的次数较大。因此采取估算的思路来处理观察日与使用具体日期的差异不大，同时这样处理能有效提升运算的效率。

2. 由于敲出观察日事先是固定的，因此对于已存续一段时间的雪球，其敲出观察日的序号需要加上已存续的交易日，这样才能使得观察日的序号对称。举个例子：某雪球已经存续60个交易日，在对第61个交易日进行估值时，模拟路径对应的第一个交易日，实际是作为雪球的第61个交易日。

上述思路确实有点绕，难以用文字完全表述清楚，只能靠读者自身去理解感悟，有种‘名可名，非常名’的意思。

2.2.2 估值定价

定价计算需要区分雪球在已存续的时间里是否发生过敲入。首先，讨论常规的情况，即在当前估值节点前，雪球未发生过敲入。

    status = 'out'# 过去未发生敲入if status == 'out':# 情况1:未来发生过敲出if len(obs_day) > 0:t = obs_day[0]re = coupon * (t / 252) * np.exp(-r * t / 252)nav = 1 * (1 + re)  # 1表示期初的净值else:# 情况2:未来未发生敲入和敲出if len(knock_in_day[0]) == 0:re = coupon * np.exp(-r * T)nav = 1 + re# 情况3:未来未发敲出，但发生敲入else:# 期末净值大于等于1if one_path[-1] >= 1:re = 0nav = 1 * (1 + re)# 期末净值小于1elif one_path[-1] < 1:re = (one_path[-1] - 1) * np.exp(-r * T)nav = 1 * (1 + re)

接下来，讨论另外一种情况，在当前估值节点前，雪球发生过敲入。相比于未发生过敲入的雪球，发生过敲入的雪球不存在未敲入和敲出的情况，即不具备获取全额票息的可能。

    elif status == 'in':# 情况1:未来发生过敲出if len(obs_day) > 0:t = obs_day[0]re = coupon * (t / 252) * np.exp(-r * t / 252)nav = 1 * (1 + re)  # 1表示期初的净值else:# 情况2:未来未发生过敲出# 期末净值大于等于1if one_path[-1] >= 1:re = 0nav = 1 * (1 + re)# 期末净值小于1elif one_path[-1] < 1:re = (one_path[-1] - 1) * np.exp(-r * T)nav = 1 * (1 + re)

在上述代码的基础上，循环5000次，计算所有nav的均值，即对应当前节点雪球的预估净值。

3. 完整代码

前文讲的相对较细，也显得有些凌乱，希望整体上是说清楚了。此处，对以上代码进行封装处理，以便查阅代码的完整性。

import numpy as npdef get_s_path(St, r, vol, residual_day, times):'''模拟标的未来可能出现的价格路径:param St: int,标的初始价格，例如1:param r: int, 无风险利率，例如0.02:param vol: int, 波动率，例如0.02:param residual_day: int, 剩余期限(单位为天)，例如252:param times: int, 模拟次数，例如5000:return:  价格路径对应的矩阵'''path_num = residual_days_path = np.zeros((path_num, times))s_path[0] = StTt = 1 / 252for i in range(1, path_num):N = np.random.standard_normal(times)  # 标准正态分布随机数ST = s_path[i - 1] * np.exp((r - 0.5 * vol ** 2) * Tt + vol * np.sqrt(Tt) * N)# 进行涨跌停的限制max_st = s_path[i - 1] * 1.1min_st = s_path[i - 1] * 0.9new_st = np.where(ST < min_st, min_st, ST)new_st = np.where(ST > max_st, max_st, new_st)s_path[i] = new_streturn s_pathdef get_one_nav(status, coupon, T, r, residual_day, one_path, knock_in, knock_out, lock_time):'''获取雪球单次的模拟路径下，对应的净值:param status: str,雪球过去的状态，'out'表示过去未发生过敲入，'in'表示过去发生过敲入:param coupon: int,雪球对应的年化票息，例如0.2:param T: int,雪球对应的期限（折算成年），例如1:param r: int,无风险利率,例如0.025:param residual_day: int,剩余天数,例如252:param one_path: 矩阵，单次对应的模拟路径:param knock_in: int,雪球对应的敲入价格:param knock_out: int，雪球对应的敲出价格:return: int,预估的雪球净值'''all_time = T * 252  # 期限，折算为天survival_day = all_time - residual_day  # 已存续期限# 价格低于敲入价格的交易日对应的序号knock_in_day = np.where(one_path <= knock_in)# 价格高于敲出价格的交易日对应的序号knock_out_day = np.where(one_path >= knock_out)# 上述每个交易日的序号需要加上雪球已存续期限new_out_day = knock_out_day[0] + survival_day# 判断上述日期是否是观察日,1年252个交易日，因此一个月设定为12obs_may_day = new_out_day[new_out_day % 21 == 0]# 由于有锁定期，需要剔除前三个月的观察日obs_day = obs_may_day[obs_may_day / 21 > lock_time]# 过去未发生敲入if status == 'out':# 情况1:未来发生过敲出if len(obs_day) > 0:t = obs_day[0]re = coupon * (t / 252) * np.exp(-r * t / 252)nav = 1 * (1 + re)  # 1表示期初的净值else:# 情况2:未来未发生敲入和敲出if len(knock_in_day[0]) == 0:re = coupon * np.exp(-r * T)nav = 1 + re# 情况3:未来未发敲出，但发生敲入else:# 期末净值大于等于1if one_path[-1] >= 1:re = 0nav = 1 * (1 + re)# 期末净值小于1elif one_path[-1] < 1:re = (one_path[-1] - 1) * np.exp(-r * T)nav = 1 * (1 + re)elif status == 'in':# 情况1:未来发生过敲出if len(obs_day) > 0:t = obs_day[0]re = coupon * (t / 252) * np.exp(-r * t / 252)nav = 1 * (1 + re)  # 1表示期初的净值else:# 情况2:未来未发生过敲出# 期末净值大于等于1if one_path[-1] >= 1:re = 0nav = 1 * (1 + re)# 期末净值小于1elif one_path[-1] < 1:re = (one_path[-1] - 1) * np.exp(-r * T)nav = 1 * (1 + re)else:nav = 0new_value = navreturn new_valuedef get_snowball_nav(status, coupon, T, r, residual_day, knock_in, knock_out, lock_time, price_path):# 获取多次模拟路径下雪球的估值# price_path:模拟的标的路径# 其余参数和get_one_nav相同nav_list = []times = np.array([price_path]).shape[-1]  # 价格矩阵的列数for num in range(times):one_path = price_path[:, num]one_nav = get_one_nav(status, coupon, T, r, residual_day, one_path, knock_in, knock_out, lock_time)nav_list.append(one_nav)return np.mean(nav_list)if __name__ == '__main__':r = 0.025vol = 0.22residual_day = 252St = 1times = 5000price_path = get_s_path(St, r, vol, residual_day, times)# 雪球要素knock_out = 1.03  # 敲出价格knock_in = 0.75  # 敲入价格coupon = 0.2  # 年化票息lock_time = 3  # 锁定期,单位为月T = 1  # 期限，折算为年status_in = 'in'in_nav = get_snowball_nav(status_in, coupon, T, r, residual_day, knock_in, knock_out, lock_time, price_path)status_out = 'out'out_nav = get_snowball_nav(status_out, coupon, T, r, residual_day, knock_in, knock_out, lock_time, price_path)

对应得到的结果为：