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

潘大师爱好者

已于 2023-09-11 21:03:28 修改

阅读量5.9k

点赞数 13

分类专栏：期权文章标签：算法 python

于 2023-09-06 16:46:50 首次发布

本文链接：https://blog.csdn.net/mhmh123789/article/details/132663291

版权

期权专栏收录该内容

12 篇文章 20 订阅

订阅专栏

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

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.1
    min_st = s_path[i-1] * 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[i] = new_st

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

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

# 单次路径可视化
import matplotlib.pyplot as plt

plt.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 plt


def 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_day
    s_path = np.zeros((path_num, times))
    s_path[0] = St
    Tt = 1 / 252
    for 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.1
        min_st = s_path[i - 1] * 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[i] = new_st
    return s_path


if __name__ == '__main__':
    r = 0.025
    vol = 0.22
    residual_day = 252
    St = 1
    times = 5000
    price_path = get_s_path(St, r, vol, residual_day, times)

对应得到的结果是：

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

    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
    plt.rcParams['axes.unicode_minus'] = False
    x = [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个交易日，因此一个月设定为12
    obs_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:
                # 期末净值大于等于1
                if one_path[-1] >= 1:
                    re = 0
                    nav = 1 * (1 + re)
                # 期末净值小于1
                elif 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:未来未发生过敲出
            # 期末净值大于等于1
            if one_path[-1] >= 1:
                re = 0
                nav = 1 * (1 + re)
            # 期末净值小于1
            elif one_path[-1] < 1:
                re = (one_path[-1] - 1) * np.exp(-r * T)
                nav = 1 * (1 + re)

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

3. 完整代码

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

import numpy as np


def 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_day
    s_path = np.zeros((path_num, times))
    s_path[0] = St
    Tt = 1 / 252
    for 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.1
        min_st = s_path[i - 1] * 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[i] = new_st
    return s_path


def 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个交易日，因此一个月设定为12
    obs_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:
                # 期末净值大于等于1
                if one_path[-1] >= 1:
                    re = 0
                    nav = 1 * (1 + re)
                # 期末净值小于1
                elif 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:未来未发生过敲出
            # 期末净值大于等于1
            if one_path[-1] >= 1:
                re = 0
                nav = 1 * (1 + re)
            # 期末净值小于1
            elif one_path[-1] < 1:
                re = (one_path[-1] - 1) * np.exp(-r * T)
                nav = 1 * (1 + re)
    else:
        nav = 0
    new_value = nav
    return new_value


def 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.025
    vol = 0.22
    residual_day = 252
    St = 1
    times = 5000
    price_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)

对应得到的结果为：