CPU-GPU 异构计算平台 Shark 应用：蒙特卡罗模拟定价

雪球产品是一类典型的自动赎回结构化产品，当标的资产的价格在特定日期超过某约定的水平时，产品触发自动赎回。因为雪球产品通常能够提供较高的票息，所以近年来得到了广大投资者的青睐，市场规模不断扩张。雪球期权作为雪球产品的重要组成部分，其期权结构本质是投资者卖出了带有复杂障碍条件的看跌期权：当标的资产上涨的时候获得有限制的收益，但当标的资产下跌的时候，可能承担和资产相同的跌幅。

面向对于雪球期权了解不多的读者，本文以一个经典的雪球期权合约为例，讲解雪球期权合约与传统多头或指数基金的损益差别，以便更好地理解后续内容。下表展示了一个典型的雪球期权合约结构以及具体参数：

对于以上典型的雪球期权合约，假设作为雪球期权的买方，最终的损益情况可以分为三种：

• 提前敲出：说明中证指数在3个月后的月频敲出观察日中某一个的收盘价 ≥ 期初价格 * 103% ，此时触发敲出条件，按照敲出票息进行结算，获得本金+25%票息
• 未敲入且未敲出：说明中证500指数在3个月后的所有月频敲出观察日的收盘价 < 期初价格 * 103%，同时日频敲入观察日的收盘价均 ≥ 期初价格 * 80%，此时既没有触发敲入条件，也没有触发敲出条件，按照红利票息进行结算，依然获得本金+25%票息
• 敲入后未敲出，即中证500指数3个月后在某一个日频敲入观察日的收盘价 < 期初价格 * 80%，同时月频敲出观察日的收盘价均 < 期初价格 * 103%，此时仅触发敲入条件而没有触发敲出条件，雪球期权亏损为 min((期末价格/期初价格-1),0) * 名义本金。假设期末价格为 St，期初价格为 S0，敲入且未敲出时合约到期的收益为 -max(S0-St，0)，相当于投资者卖出了一份行权价为 S0 的看跌期权的损益。

因而，将不同的雪球期权损益情形总结起来，就是下图流程图展示的损益结算规则：

同时，对于雪球期权的发行方而言，其损益与上述买方损益恰好相反。雪球期权发行方的目标是在给定目标利润率与其他参数的情况下，为其发行的雪球期权产品制定合适的票息率。

对于期权的数值定价方法，主要有 BS 公式法，有限差分法与蒙特卡洛法三种方法。BS 方法是期权定价的经典解析方法，其核心思想是通过构建无风险对冲组合，推导出期权价格满足的 PDE ，并求得闭式解。有限差分即是对期权 PDE 中的偏导用有限差分替代，进而通过构造网格求解三对角方程组的方式求解出期权本身的价值，但其计算方法相对复杂。蒙特卡洛则是通过假设标的资产价格为某个随机过程并大量模拟该随机过程下的轨道，进而得到期权的数值解。其优点一是实现起来非常简单，二是由于误差收敛率不依赖于问题的维数，从而非常适用于高维期权定价，常用于路径依赖 (path-dependent) 的各类奇异期权，以及多资产期权 (multi-asset) 模型。

本文所介绍的雪球期权属于奇异期权，其定价通常依赖于蒙特卡洛模拟方法，蒙特卡洛方法实现期权定价的基本原理如下：

假设风险中性测度下，资产价格服从 GBM (Geometric Brownian Motion, 几何布朗运动)，则满足：

其中：St 为 t 时刻的资产价格，r 为无风险利率，σ 为资产价格的波动率，Wt 为维纳过程。在离散时间形式下，雪球期权对应标的资产价格路径可以通过以下公式生成：

在模拟出价格路径的情况下，雪球期权买方仅需根据模拟出的价格路径，按照雪球期权损益规则计算期权现值，对不同价格路径对应的现值取均值即可得到雪球期权的价值：

对于雪球期权的发行方而言，其目标是给定利润率计算雪球的票息率（敲出票息与红利票息），仍然遵循上述步骤模拟出价格路径并计算雪球现值。之后可利用二分法等数值方法对票息率进行逼近，从而求出满足目标收益率的票息率。

利用 seq 函数与 repmat 函数构造指定数量的时间序列矩阵，并利用 norm 函数生成正态分布矩阵，将时间序列矩阵与正态分布矩阵带入 Black Sholes 公式构造标的价格矩阵。进一步地，本案例设置了一定的涨跌幅限制（这里是日上下10%），从而更贴合实际交易环境。该部分的实例代码如下所示：

@gpu
def simulation(S, r, T, sigma, N, days_in_year){
    // 计算每日时间间隔, 折算成年的时间差
    dt = 1 \ days_in_year
    // 计算到期时间, 折算成日
    Tdays = T*days_in_year
    // 月交易日数 （每个月观察一次是否敲出）
    monthDay = int(days_in_year/12)     
    all_t = dt*(0..Tdays)
    // 一列是一个时间序列，有 N 列
    dtMatrix = repmat(matrix(all_t),1,N)    
    normMatrix = norm(0,1,(Tdays+1)*N)$(Tdays+1):N
    normMatrix = normMatrix * (all_t>0)
    // 一列是一个价格路径，有 N 列
    ST = S * exp((r-0.5*sigma*sigma)*dtMatrix + sigma*sqrt(dt)*cumsum(normMatrix)) 
    Spath = ST[1..Tdays, ]
    // 进行涨跌停限制
    prevPrice = ST[0..(Tdays-1), ]
    uplimit = prevPrice*1.1
    lowlimit = prevPrice*0.9
    Spath = iif(Spath > uplimit, uplimit, iif(Spath < lowlimit, lowlimit, Spath))
    return Spath
}

模拟出的价格路径矩阵 Spath 的可视化如图所示，这里只展示前 25 条标的价格路径。

如上文所介绍，用户对于仅使用 Shark Graph 支持数据格式与相关函数的函数加上 @gpu 注解，即可使用DolphinDB 脚本在 GPU 上加速计算。在获得标的价格矩阵后，可以利用 DolphinDB 的向量化编程，对于雪球期权每条价格路线的损益状态进行向量化判断，得到每一条价格路径的雪球期权现值，取均值即为最终雪球期权的定价结果。

以下是 Shark 雪球期权定价的主体函数，该函数会返回每一条雪球期权价值组成的向量 payoff，向下敲入次数 knock_in_times ，向上敲出次数 knock_out_times ，正常到期次数 existence_times 与亏损次数 lose_times ：

@gpu 
def snowball_pricing(S,r,T,sigma,N,coupon,k_out,k_in,lock_period,days_in_year){ 
    // 计算到期时间, 折算成日
    Tdays = T*days_in_year         
    // 月交易日数 （每个月观察一次是否敲出）
    monthDay = int(days_in_year/12)
    Spath = simulation(S, r, T, sigma, N, days_in_year)
    // 构造需要观察是否敲入的时间序列（超过封闭期）
    observation_out_idx = monthDay*(1..(Tdays/monthDay))
    observation_out_idx = observation_out_idx[observation_out_idx>lock_period]-1
    // 判断每个月是否需要敲出
    knockOutMatrix = Spath[observation_out_idx, ]>=k_out
    // 获取敲出的日期
    knock_out_day = observation_out_idx*knockOutMatrix
    // 计算首次敲出的日期
    knock_out_day = min(knock_out_day[knock_out_day>0])
    // 判断是否需要敲入
    knock_in_day = sum(Spath < k_in)

    // 计算收益率：
    // 情景1：发生了向上敲出事件: 此时按照合约规定，该雪球类产品提前终止
    // 投资者收益=年化票息率*合约持有期限年化
    payoff1 = coupon*knock_out_day/days_in_year*exp(-r * knock_out_day/days_in_year)
    // 情景2：未敲出且未敲入: 此时按照合约规定， 该雪球类产品到期自动终止
    // 投资者收益=年化票息率*合约持有期限年化
    payoff2 = coupon * exp(-r * T)
    // 情景3：只发生向下敲入，不发生向上敲出, 且合约到期时，标的价格低于期初标的价格
    // 此时按照合约规定，该雪球类产品到期自动终止，此时敲入的看跌期权为实值，投资者需要承担标的下跌造成的损失。
    payoff3 = (last(Spath) - S)*exp(-r * T)
    // 判断每一个价格路径的收益状态
    payoff = iif(knock_out_day>0, payoff1,
             iif(knock_in_day==0, payoff2,
             iif(last(Spath)<S, payoff3, 0)))
  
    // 统计敲出次数：
    knock_out_times = sum(knock_out_day>0)
    // 统计未敲出且未敲入次数 
    existence_times = sum(knock_out_day<=0&&knock_in_day==0)
    // 统计敲入次数 
    knock_in_times = sum(knock_out_day<=0&&knock_in_day>0)
    // 统计亏损次数
    lose_times = sum(knock_out_day<=0&&knock_in_day>0&&last(Spath)<S)
    return payoff, knock_out_times, knock_in_times, existence_times,lose_times
}

蒙特卡洛方法进行雪球期权定价的模拟路径数量、标的初始价格、敲入边界价格、敲出边界价格、标的波动率等均可以通过参数进行配置，本案例所用到的雪球期权定价参数含义如下表所示：

参数配置与具体的运行代码如下所示：

N = 300000          // 模拟路径数量
S = 1.0             // 标的初始价格
k_in = 0.85         // 敲入边界
k_out = 1.03        // 敲出边界
T = 1               // 到期时间（年）
days_in_year = 252  // 一年的交易日数
sigma = 0.13        // 标的波动率（年化）
r = 0.03            // 无风险利率（年化）
coupon = 0.2        // 贴水利率（年化）
lock_period = 90     // 锁仓期（日）

timer payoff, knock_out_times, knock_in_times, existence_times, lose_times = snowball_pricing(S, r, T, sigma, N, coupon, k_out, k_in, lock_period, days_in_year)
print("雪球价值:", avg(payoff))
print("敲出概率:", knock_out_times \ N)
print("未敲入也未敲出概率:", existence_times \ N)
print("敲入且未敲出概率:", knock_in_times \ N)
print("亏损概率:", lose_times \ N)

使用 Shark Graph 加速运算的输出结果如下所示，可以看出，以上参数对应的雪球期权预期收益率为 5%，敲出概率为 73%，具备一定的投资价值，雪球期权定价的运行速度不到 40ms。

Time elapsed: 36.827 ms

雪球价值:
0.050549256728865
敲出概率:
0.737543333333333
未敲入也未敲出概率:
0.136706666666667
敲入且未敲出概率:
0.12575
亏损概率:
0.123776666666667

首先，不同于雪球期权定价，给定利润率进行雪球期权票息率的确定需要利用数值方法对相同的价格矩阵进行雪球价值计算，因而首先单独定义雪球价值计算函数，依然使用 @gpu 注解利用 Shark Graph 进行 GPU 加速计算：

@gpu
def calculate_pay_off(coupon, Spath, S, r, T, knock_out_day, knock_in_day, days_in_year){
    // 分三种情景计算收益率：
    payoff1 = coupon * knock_out_day / days_in_year * exp(-r * knock_out_day / days_in_year)
    payoff2 = coupon * exp(-r * T)
    payoff3 = (last(Spath) - S)*exp(-r * T)
    // 判断每一个价格路径的收益状态
    payoff = iif(knock_out_day>0, payoff1,
             iif(knock_in_day==0, payoff2,
             iif(last(Spath)<S, payoff3, 0)))
    return payoff
}

之后，这里使用了等距网格搜索，固定间隔生成雪球期权票息率的可能范围，间距为 0.001，在此基础上进行搜索票息率，直至前后两个票息率对应的现值刚好包含了发行方期望的现值：

@gpu
def snowball_pricing(S, r, T, sigma, N, k_out, k_in, lock_period, days_in_year, profit_rate){ 
    Tdays = T*days_in_year              
    monthDay = int(days_in_year/12)    
    Spath = simulation(S, r, T, sigma, N, days_in_year) // 模拟价格路径矩阵
    observation_out_idx = monthDay*(1..(Tdays/monthDay))
    observation_out_idx = observation_out_idx[observation_out_idx>lock_period]-1
    knockOutMatrix = Spath[observation_out_idx, ]>=k_out
    knock_out_day = observation_out_idx*knockOutMatrix
    knock_out_day = min(knock_out_day[knock_out_day>0])
    knock_in_day = sum(Spath < k_in) 
    knock_out_times = sum(knock_out_day>0) // 模拟价格路径发生敲入的数量
    existence_times = sum(knock_out_day<=0&&knock_in_day==0) // 模拟价格路径发生未敲入未敲出的数量
    knock_in_times = sum(knock_out_day<=0&&knock_in_day>0) // 模拟价格路径发生期末敲入未敲出的数量
    lose_times = sum(knock_out_day<=0&&knock_in_day>0&&last(Spath)<S) // 投资者亏损次数(发行方盈利次数)
    
    // 固定间隔生成票息率可能范围
    coupon_test_range = (0..1000)*0.001 // 精确到小数点后四位

    // for loop循环搜索
    lastvalue = 0
    lastcoupon = 0
    target = (1-profit_rate)
    for (coupon in coupon_test_range){
        // 计算路径对应的现值 
        value = 1.0 + avg(calculate_pay_off(coupon, Spath, S, r, T, knock_out_day, knock_in_day, days_in_year)) 
        if (value > target > lastvalue){
            break
        }
        lastvalue = value
        lastcoupon = coupon
    }
    // 获取区间边界
    coupon_lower_bound = lastcoupon // coupon上边界
    coupon_upper_bound = coupon // coupon下边界
    return 0.5*(coupon_lower_bound + coupon_upper_bound), knock_out_times, knock_in_times, existence_times, lose_times
}

下表列示了用于确定雪球发行方票息率的函数的相关参数，除了profit_rate 参数其余参数均与雪球定价函数的参数保持一致。

参数配置与具体的运行代码如下所示：

N = 300000          // 模拟路径条数
S = 1.0             // 标的初始价格
k_in = 0.85         // 敲入边界
k_out = 1.03        // 敲出边界
T = 1               // 到期时间（年）
days_in_year = 252  // 一年的交易日数
sigma = 0.13        // 标的波动率（年化）
r = 0.03            // 无风险利率（年化）
lock_period = 0     // 锁仓期（日）
profit_rate = 0.01  // 雪球期权发行方目标利润率

timer coupon, knock_out_times, knock_in_times, existence_times, lose_times = snowball_pricing(S, r, T, sigma, N, k_out, k_in, lock_period, days_in_year, profit_rate)
print("雪球票息率:", coupon)
print("敲出概率:", knock_out_times \ N)
print("未敲入也未敲出概率:", existence_times \ N)
print("敲入且未敲出概率:", knock_in_times \ N)
print("发行方盈利概率:", lose_times \ N)  // 发行方盈利概率=买方亏损概率

使用 Shark Graph 加速运算的输出结果如下所示，可以看出，以上参数对应的雪球期权票息率为 2.35%。

Time elapsed: 116.619 ms

雪球票息率:
0.0235
敲出概率:
0.7385
未敲入也未敲出概率:
0.135646666666667
敲入且未敲出概率:
0.125853333333333
亏损概率:
0.12393