雪球结构产品定价
'''
Author: dkl
Description: MC雪球定价
Date: 2023-08-29 08:34:49
'''
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from tqdm import tqdm
mpl.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号
mpl.rcParams['font.sans-serif'] = ['SimHei']   # 用来正常显示中文标签
plt.rcParams['font.sans-serif'] = ['SimHei']   # 用来正常显示中文标签
np.random.seed(0)                              # fix random seed


def stock_monte_carlo(start_price, T, r, q, sigma, steps, simulations):
    '''
    Desctription
    ----------
    蒙特卡洛模拟股票价格

    Parameters
    ----------
    start_price : float
        模拟起始股价.
    T : float
        模拟时间长度，如T=1表示持续1年
    r : float
        股票收益率（年化）.
    sigma : float
        股票波动率（年化）.
    steps : float
        往后预测步数.
    simulations : float
        蒙特卡洛模拟次数.

    Returns
    -------
    price : numpy.array
        模拟股价记录矩阵.

    '''

    dt = T/steps
    # Define a price array
    price = np.zeros((steps + 1, simulations))
    price[0] = start_price
    for i in range(1, steps + 1):
        z = np.random.standard_normal(simulations)
        price[i] = price[i - 1] * \
            np.exp(((r-q) - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z)
    return price.T


def snowball_valuation(prices, B1, B2, strike_out_payment, div_payment, rf):
    """
    Description
    ----------
    雪球结构定价

    Parameters
    ----------
    prices: numpy.ndarray
        蒙特卡洛所得的一条模拟价格路径.
    B1: float
        雪球的敲入价格.
    B2: float
        雪球的敲出价格.
    strike_out_payment: float
        敲出票息.
    div_payment: float
        红利票息.
    rf: float
        无风险收益率（年化）.

    Returns
    -------
    value: float
        结构定价.
    ret: float
        期间收益.
    down_flag: bool
        是否敲入.
    up_flag: bool
        是否敲出.
    end_day: int
        存续天数.

    """
    value = None
    down_flag = False                        # 是否已经触发敲入
    up_flag = False                          # 是否触发敲出
    days = np.arange(0, prices.shape[0], 1)  # 包括起始日一共有253天
    up_B2 = (prices >= B2)                   # 高于敲出价格的日期
    KO_observe = (days % 21 == 0) & (days > 0)          # 敲出观察日
    KO_days = days[up_B2*KO_observe]         # 敲出日
    KI_days = days[prices < B1]              # 敲入日
    up_flag = (len(KO_days) > 0)
    down_flag = (len(KI_days) > 0)
    if len(KO_days) > 0:                     # 曾敲出
        end_day = KO_days[0]                 # 敲出日期
        value = strike_out_payment*(end_day/252)*np.exp(-rf*end_day/252)
    elif len(KO_days) == 0 and len(KI_days) > 0:  # 曾敲入
        end_day = 252                        # 持有至到期
        if prices[-1] < prices[0]:             # 如果到期标的小于期初价格
            value = (prices[-1]-prices[0])/prices[0] * np.exp(-rf*end_day/252)
        else:                                # 如果到期标的高于期初价格，获益为0
            value = 0
    elif len(KO_days) == 0 and len(KI_days) == 0:  # 未敲入也未敲出
        end_day = 252
        value = div_payment*(end_day/252)*np.exp(-rf*end_day/252)
    else:
        raise ValueError
    return value


def main(S0=1):
    # 用于MC的常数
    steps = 252
    T = 1
    r = 0.03
    q = 0
    sigma = 0.13
    simulation_time = 5000
    # 雪球结构参数
    B2 = 1.03  # 敲出界限：103% 每个月观察敲出
    B1 = 0.85  # 敲入界限：85% 每天观察敲入
    strike_out_payment = 0.2  # 敲出票息=20%（年化）
    div_payment = 0.2  # 红利票息=20%（年化）
    rf = 0.03
    record = stock_monte_carlo(S0, T, r, q, sigma, steps, simulation_time)
    value_lst = []
    for i in tqdm(range(len(record)), leave=False):
        prices = record[i]
        value = snowball_valuation(
            prices, B1, B2, strike_out_payment, div_payment, rf)
        value_lst.append(value)
    # 组合价值
    return value_lst

if __name__ == '__main__':
    # 雪球定价
    value_lst = main()
    print('MC下的雪球结构定价:', np.mean(value_lst))
    # 雪球收益分布
    plt.figure(figsize=(10, 5))
    plt.hist(value_lst, bins=30)
    plt.title('雪球收益分布')
    plt.savefig('MC下的雪球收益分布.png')
    plt.show()
    # 雪球delta计算
    start_price = 1.03*1.5
    end_price = 0.85*0.7
    S_arr = np.linspace(start_price, end_price, 250)
    fvalue_arr = np.zeros(250)
    for i in tqdm(range(len(S_arr))):
        S0 = S_arr[i]
        value_lst = main(S0)
        fvalue_arr[i] = np.mean(value_lst)
    delta_arr = (fvalue_arr[2:] - fvalue_arr[:-2])/(S_arr[2:]-S_arr[:-2])
    S_arr = S_arr[1:-1]
    plt.figure(figsize=(10, 5))
    plt.plot(S_arr, delta_arr)
    plt.title('MC下的雪球delta')
    plt.savefig('MC下的雪球delta.png')
    plt.show()

2_雪球蒙特卡洛定价/snowball_mc.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: 基本有限差分类
Date: 2023-04-22 08:33:21
'''
import numpy as np
from tqdm import tqdm


class BasicFDM(object):
    '''
    基本有限差分类, 包含显示、隐式、Crank-Nicolson三种定价方法
    '''

    def __init__(self,
                 S0,
                 r,
                 q,
                 T,
                 sigma,
                 S_max,
                 S_min,
                 M,
                 N,
                 method='Crank-Nicolson'):
        '''
        Description
        ----------
        构造函数

        Parameters
        ----------
        S0: float. 标的资产初始价格
        r: float. 无风险利率
        q: float. 分红率
        T: float. 时间
        sigma: float. 标的资产波动率
        S_max: float. 标的资产最大值
        S_min: float. 标的资产最小值
        M: int. 标的价格网格的分段数
        N: int. 时间网格的分段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        # 期权基本参数
        self.S0 = S0
        self.r = r
        self.q = q
        self.T = T
        self.sigma = sigma
        # 网格参数
        self.S_max = S_max
        self.S_min = S_min
        self.M = M
        self.N = N
        self.delta_S = (S_max - S_min) / M
        self.delta_T = T / N
        # 有限差分计算系数
        self.coef = None
        # 期权网格
        # 第一个维度是时间，第二个维度是标的价格
        self.f_mat = np.zeros((N + 1, M + 1))
        # 是否已设置边界条件
        self._flag_set_boundary = False
        if method not in ['Explicit', 'Implicit', 'Crank-Nicolson']:
            err_message = 'method must be in Explicit, Implicit, Crank-Nicolson'
            raise ValueError(err_message)
        self.method = method

    def set_boundary(self):
        '''
        Description
        ----------
        设置边界

        Parameters
        ----------
        None.

        Return
        ----------
        None
        '''
        self._flag_set_boundary = True
        return

    def get_abc(self):
        '''
        Description
        ----------
        计算abc

        Parameters
        ----------
        None.

        Return
        ----------
        tuple[numpy.array]. 每个点对应的是a1,a2到aM-1
        '''
        # Smin到Smax共M+1个节点，标号为0到M
        arr = np.linspace(self.S_min, self.S_max, self.M + 1) / self.delta_S
        # 但我们只需要1到M-1
        arr = arr[1: self.M]
        # 计算参数
        rjt = (self.r - self.q) * arr * self.delta_T
        sigma2j2t = (self.sigma**2) * (arr**2) * self.delta_T
        if self.method == 'Explicit':
            discount = 1 / (1 + (self.r-self.q) * self.delta_T)
            a = discount * (-0.5 * rjt + 0.5 * sigma2j2t)
            b = discount * (1 - sigma2j2t)
            c = discount * (0.5 * rjt + 0.5 * sigma2j2t)
        elif self.method == 'Implicit':
            a = 0.5 * rjt - 0.5 * sigma2j2t
            b = 1 + (self.r-self.q)*self.delta_T + sigma2j2t
            c = -0.5 * rjt - 0.5 * sigma2j2t
        else:
            a = -0.25 * rjt + 0.25 * sigma2j2t
            b = -0.5 * (self.r-self.q) * self.delta_T - 0.5 * sigma2j2t
            c = 0.25 * rjt + 0.25 * sigma2j2t
        return a, b, c

    def calc_coef(self):
        '''
        Description
        ----------
        计算有限差分需要的系数

        Parameters
        ----------
        None.

        Return
        ----------
        None.
        '''
        if self.method == 'Explicit':
            self.explicit_calc_coef()
        elif self.method == 'Implicit':
            self.implicit_calc_coef()
        else:
            self.crank_nicolson_calc_coef()

    def explicit_calc_coef(self):
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵P
        P = np.zeros((self.M - 1, self.M - 1))
        # 计算第一行
        a, b, c = self.get_abc()
        P[0, 0] = b[0]
        P[0, 1] = c[0]
        # 第2行到第M-2行
        for i in range(2, self.M - 1):
            P[i - 1, i - 2] = a[i - 1]
            P[i - 1, i - 1] = b[i - 1]
            P[i - 1, i] = c[i - 1]
        # 计算第M-1行
        P[self.M - 2, self.M - 3] = a[self.M - 2]
        P[self.M - 2, self.M - 2] = b[self.M - 2]
        # 计算b
        tempb1 = (a[0] * self.f_mat[:, 0]).reshape(-1, 1)
        tempb2 = np.zeros((self.N + 1, self.M - 3))
        tempb3 = (c[self.M - 2] * self.f_mat[:, self.M]).reshape(-1, 1)
        Q = np.hstack((tempb1, tempb2, tempb3))
        self.coef = (P, Q)
        return

    def implicit_calc_coef(self):
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵P
        P = np.zeros((self.M - 1, self.M - 1))
        # 计算abc,为a1-M-1
        a, b, c = self.get_abc()
        # 第一行为b1,c1
        P[0, 0] = b[0]
        P[0, 1] = c[0]
        # 第2行到第M-2行为0, ai, bi, ci, 0
        for i in range(2, self.M - 1):
            P[i - 1, i - 2] = a[i - 1]
            P[i - 1, i - 1] = b[i - 1]
            P[i - 1, i] = c[i - 1]
        # 计算第M-1行
        P[self.M - 2, self.M - 3] = a[self.M - 2]
        P[self.M - 2, self.M - 2] = b[self.M - 2]
        # 计算Q
        tempb1 = (a[0] * self.f_mat[:, 0]).reshape(-1, 1)
        tempb2 = np.zeros((self.N + 1, self.M - 3))
        tempb3 = (c[self.M - 2] * self.f_mat[:, self.M]).reshape(-1, 1)
        Q = np.hstack((tempb1, tempb2, tempb3))
        Pinv = np.linalg.inv(P)
        self.coef = (Pinv, Q)
        return

    def crank_nicolson_calc_coef(self):
        # 计算abc
        a, b, c = self.get_abc()
        # 计算P1
        P1 = np.zeros((self.M - 1, self.M - 1))
        # 第一行为1-b1,-c1,0,…,0
        P1[0, 0] = 1 - b[0]
        P1[0, 1] = -c[0]
        # 第二行及以后为(i-2)个0，-ai,1-bi,-ci,…,0
        for i in range(2, self.M - 1):
            P1[i - 1, i - 2] = -a[i - 1]
            P1[i - 1, i - 1] = 1 - b[i - 1]
            P1[i - 1, i] = -c[i - 1]
        P1[self.M - 2, self.M - 3] = -a[self.M - 2]
        P1[self.M - 2, self.M - 2] = 1 - b[self.M - 2]

        # 计算P2
        P2 = np.zeros((self.M - 1, self.M - 1))
        # 第一行为1+b1,c1,0,…,0
        P2[0, 0] = 1 + b[0]
        P2[0, 1] = c[0]
        # 第二行及以后为(i-2)个0，ai,bi,ci,…,0
        for i in range(2, self.M - 1):
            P2[i - 1, i - 2] = a[i - 1]
            P2[i - 1, i - 1] = 1 + b[i - 1]
            P2[i - 1, i] = c[i - 1]
        P2[self.M - 2, self.M - 3] = a[self.M - 2]
        P2[self.M - 2, self.M - 2] = 1 + b[self.M - 2]
        # 计算Q
        tempb1 = (a[0] * self.f_mat[:, 0]).reshape(-1, 1)
        tempb2 = np.zeros((self.N + 1, self.M - 3))
        tempb3 = (c[self.M - 2] * self.f_mat[:, self.M]).reshape(-1, 1)
        Q = np.hstack((tempb1, tempb2, tempb3))
        P1inv = np.linalg.inv(P1)
        self.coef = (P1inv, P2, Q)
        return

    def get_fi_1(self, i):
        '''
        Description
        ----------
        从第i期开始倒推i-1期的期权价值

        Parameters
        ----------
        i: int. 时间编号

        Return
        ----------
        np.array.第i-1期的期权价值
        '''
        # 计算相应的系数，按照时间进行倒推
        fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
        if self.method == 'Explicit':
            P, Q = self.coef
            Qi = Q[i, :].reshape(-1, 1)
            fi_1 = P @ fi + Qi
        elif self.method == 'Implicit':
            Pinv, Q = self.coef
            Qi_1 = Q[i-1, :].reshape(-1, 1)
            fi_1 = Pinv @ (fi - Qi_1)
        else:
            P1inv, P2, Q = self.coef
            Qi = Q[i, :].reshape(-1, 1)
            Qi_1 = Q[i-1, :].reshape(-1, 1)
            fi_1 = P1inv@(P2 @ fi + Qi + Qi_1)
        fi_1 = fi_1.reshape(-1)
        return fi_1

    def calc_f_mat(self):
        '''
        Description
        ----------
        计算期权价格矩阵

        Parameters
        ----------
        None.

        Return
        ----------
        None.
        '''
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算相应的系数，按照时间进行倒推
        self.calc_coef()
        for i in tqdm(range(self.N, 0, -1)):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            # 赋值给下个fi
            self.f_mat[i - 1, 1:self.M] = fi_1
        return

    def interpolate_f_value(self):
        '''
        Description
        ----------
        找到S0相邻的网格，进行插值计算.

        Parameters
        ----------
        None.

        Return
        ----------
        float.S0对应的期初期权价格
        '''
        idx = int((self.S0 - self.S_min) / self.delta_S)
        f_upper = self.f_mat[0, idx + 1]
        f_lower = self.f_mat[0, idx]
        s_upper = self.S_min + self.delta_S * (idx + 1)
        s_lower = self.S_min + self.delta_S * idx
        df_ds = (f_upper - f_lower) / (s_upper - s_lower)
        f_value = f_lower + df_ds * (self.S0 - s_lower)
        return f_value

    def get_f_value(self):
        '''
        Description
        ----------
        计算期初期权价值.在初始化后即可直接调用.

        Parameters
        ----------
        None.

        Return
        ----------
        float.S0对应的期初期权价格
        '''
        # 设置边界条件
        self.set_boundary()
        self.calc_f_mat()
        # 线性插值获得期权价值
        f_value = self.interpolate_f_value()
        return f_value

    def get_f_delta(self, S, flag_calc_f_mat=False):
        '''
        Description
        ----------
        计算期权在t时刻和S价格位置对应的delta值
        
        Parameters
        ----------
        t: 时刻
        S: 标的价格
        flag_calc_f_mat: 是否已经计算期权价格矩阵

        Return
        ----------
        float. delta值
        '''
        if not flag_calc_f_mat:
            self.set_boundary()
            # 计算边界价格
            self.calc_f_mat()
        idx_S = int((S - self.S_min) / self.delta_S)
        f2 = self.f_mat[0, idx_S+1]
        f1 = self.f_mat[0, idx_S-1]
        delta = (f2 - f1)/self.delta_S
        return delta
    
    def get_f_theta(self, t, flag_calc_f_mat=False):
        if not flag_calc_f_mat:
            self.set_boundary()
            # 计算边界价格
            self.calc_f_mat()
        idx_S = int((self.S0-self.S_min)/self.delta_S)
        idx_T = int(t/self.delta_T)
        f1 = self.f_mat[idx_T, idx_S]
        f2 = self.f_mat[idx_T+1, idx_S]
        theta = (f2 - f1)/self.delta_T
        return theta
3_雪球有限差分定价/basic_fdm.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: BS公式
Date: 2023-04-18 15:24:01
'''
import numpy as np
from scipy.stats import norm


def euro_bs_call(S0, K, r, T, sigma):
    d1 = (np.log(S0 / K) + (r + 0.5 * sigma**2) * T) / (sigma * np.sqrt(T))
    d2 = d1 - sigma * np.sqrt(T)
    return S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)


def euro_bs_put(S0, K, r, T, sigma):
    call = euro_bs_call(S0, K, r, T, sigma)
    put = call + K*np.exp(-r*T) - S0
    return put

3_雪球有限差分定价/bs.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: FDM雪球定价
Date: 2023-08-28 22:54:39
'''
import numpy as np
from tqdm import tqdm
from basic_fdm import BasicFDM
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号


class AutocallFDM(BasicFDM):
    '''
    Autocall期权定价
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min时，期权价格为0
        self.f_mat[:, 0] = 0
        # 边界条件2：到期时间T时，在敲出部分获得收益
        idx = int((self.barrier_ko - self.S_min)/self.delta_S)
        self.f_mat[self.N, 0:idx+1] = 0
        self.f_mat[self.N, idx+1:] = self.c * self.T
        # 边界条件3：敲出时才会获得收益
        for i in range(21, self.N+1, 21):
            # 敲出日
            self.f_mat[i, idx:self.M+1] = self.c*i*self.delta_T
            # 问题
            t_arr = np.arange(i-21, i)*self.delta_T
            # 敲出部分: 为下期收益的贴现
            self.f_mat[i-21:i, self.M] = self.c * (i * self.delta_T) * np.exp(-self.r*(i*self.delta_T-t_arr))
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算相应的系数，按照时间进行倒推
        self.calc_coef()
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            if ((i - 1) % 21 == 0) & (i != 1):
                idx = int((self.barrier_ko - self.S_min)/self.delta_S)
                self.f_mat[i - 1, 1:idx+1] = fi_1[0:idx]
            else:
                # 赋值给下个fi
                self.f_mat[i - 1, 1:self.M] = fi_1
        return


class DNTFDM(BasicFDM):
    '''
    双边触碰失效期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 barrier_ki,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        barrier_ki: float. 敲入价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ki = barrier_ki
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min或者敲入时，期权价格为0
        idx_in = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        self.f_mat[:, 0] = 0
        self.f_mat[:, 0:idx_in+1] = 0
        # 边界条件2：到期时间T时，仅在非敲入敲出部分获得全部票息收益
        self.f_mat[self.N, idx_in+1:idx_out+1] = self.c * self.T
        self.f_mat[self.N, 0:idx_in+1] = 0
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者敲出时，期权价格为0
        for i in range(21, self.N+1, 21):
            # 股票价格为Smax
            self.f_mat[i-21:i, self.M] = 0
            # 敲出部分
            self.f_mat[i, idx_out+1:self.M+1] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_ki = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_ko = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            # 如果是敲出日，范围应该是敲入到敲出
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, idx_ki+1:idx_ko+1] = fi_1[idx_ki:idx_ko]
            else:
                # 非敲出日范围应该是敲入到最后
                self.f_mat[i - 1, idx_ki+1:self.M] = fi_1[idx_ki:self.M-1]
        return


class UOPFDM(BasicFDM):
    '''
    上涨失效看跌期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min时，不会敲出，且看跌期权必定执行
        ti = np.linspace(0, self.T, self.N+1)
        self.f_mat[:, 0] = self.S0*np.exp(-self.r*(self.T-ti)) - self.S_min
        # 边界条件2：到期时间T时，不敲出收益为max((S0-ST)/S0, 0),敲出收益为0
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        si = np.linspace(self.S_min, self.S_max, self.M+1)
        self.f_mat[self.N, 0:idx_out + 1] = \
            np.maximum((self.S0 - si)/self.S0, 0)[0:idx_out+1]
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者敲出时，收益为0
        for i in range(21, self.N+1, 21):
            self.f_mat[i, idx_out+1:] = 0
            self.f_mat[i-21:i, self.M] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, 1:idx_out+1] = fi_1.reshape(-1)[0:idx_out]
            else:
                self.f_mat[i - 1, 1:self.M] = fi_1.reshape(-1)
        return


class DKOPFDM(BasicFDM):
    '''
    双边失效看跌期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 barrier_ki,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        barrier_ki: float. 敲入价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ki = barrier_ki
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        idx_in = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        # 边界条件1： 股票价格为S_min或者向下敲出时，期权价格为0
        self.f_mat[:, 0] = 0
        self.f_mat[:, 0:idx_in+1] = 0
        # 边界条件2：到期时间T时，仅在非敲入敲出部分获得看跌期权收益
        si = np.linspace(self.S_min, self.S_max, self.M+1)
        vanillaT = np.maximum((self.S0 - si)/self.S0, 0)
        self.f_mat[self.N, idx_in+1:idx_out+1] = vanillaT[idx_in+1:idx_out+1]
        self.f_mat[self.N, 0:idx_in+1] = 0
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者向上敲出时，期权价格为0
        for i in range(21, self.N+1, 21):
            # 股票价格为Smax
            self.f_mat[i-21:i, self.M] = 0
            # 敲出部分
            self.f_mat[i, idx_out+1:] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_ki = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_ko = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            # 如果是敲出日，范围应该是敲入到敲出
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, idx_ki+1:idx_ko+1] = fi_1[idx_ki:idx_ko]
            else:
                # 非敲出日范围应该是敲入到最后
                self.f_mat[i - 1, idx_ki+1:self.M] = fi_1[idx_ki:self.M-1]
        return


class SnowballFDM(object):
    def __init__(self,
                 S0=1,
                 barrier_ko=1.03,
                 barrier_ki=0.85,
                 c=0.2,
                 r=0.03,
                 q=0,
                 T=1,
                 sigma=0.13,
                 S_min=0,
                 S_max=4,
                 M=1000,
                 N=252,
                 method='Crank-Nicolson'
                 ):
        self.S0 = S0
        self.barrier_ko = barrier_ko
        self.barrier_ki = barrier_ki
        self.c = c
        self.r = r
        self.q = q
        self.T = T
        self.sigma = sigma
        self.S_min = S_min
        self.S_max = S_max
        self.M = M
        self.N = N
        self.delta_S = (S_max - S_min) / M
        self.delta_T = T / N
        self.method = method

    def get_snowball_fvalue(self):
        fdm1 = AutocallFDM(self.S0, self.barrier_ko, self.c, self.r, self.q,
                           self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm2 = DNTFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, self.r, self.q,
                      self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm3 = UOPFDM(self.S0, self.barrier_ko, self.c, self.r, self.q, self.T,
                      self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm4 = DKOPFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, self.r, self.q,
                       self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        f1 = fdm1.get_f_value()
        f2 = fdm2.get_f_value()
        f3 = fdm3.get_f_value()
        f4 = fdm4.get_f_value()
        fvalue = f1 + f2 - f3 + f4
        return fvalue

    def get_snowball_f_mat(self, sigma=None, r=None):
        if sigma is None:
            sigma = self.sigma
        if r is None:
            r = self.r
        fdm1 = AutocallFDM(self.S0, self.barrier_ko, self.c, r, self.q,
                           self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm2 = DNTFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, r, self.q,
                      self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm3 = UOPFDM(self.S0, self.barrier_ko, self.c, r, self.q, self.T,
                      sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm4 = DKOPFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, r, self.q,
                       self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm1.set_boundary()
        fdm2.set_boundary()
        fdm3.set_boundary()
        fdm4.set_boundary()
        fdm1.calc_f_mat()
        fdm2.calc_f_mat()
        fdm3.calc_f_mat()
        fdm4.calc_f_mat()
        f_mat = fdm1.f_mat + fdm2.f_mat - fdm3.f_mat + fdm4.f_mat
        return f_mat

    def plot_snowball_delta(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算delta矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            delta_mat = (f_mat[:, 2:] - f_mat[:, :-2]) / np.tile((S_arr[2:] - S_arr[:-2]).reshape(1, -1), (self.N+1, 1))
            S_arr = S_arr[1:-1]
            # 选取Smin到Smax的值
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            delta_mat = delta_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, delta_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (delta_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, delta_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('delta(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球delta_3D.png')
            plt.show()
        else:
            fvalue_arr = f_mat[0, :]
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # delta_arr对应的是Smin+deltaS到Smax-deltaS之间的值
            delta_arr = (fvalue_arr[2:] - fvalue_arr[:-2])/(S_arr[2:]-S_arr[:-2])
            S_arr = S_arr[1:-1]
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            delta_arr = delta_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('delta(FDM)')
            plt.plot(S_arr, delta_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球delta.png')
            plt.show()
        return

    def plot_snowball_gamma(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算gamma矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            d2f = f_mat[:, 2:] + f_mat[:, :-2] - 2*f_mat[:, 1:-1]
            ds2 = ((S_arr[2:]-S_arr[:-2])/2)**2
            ds2 = np.tile((ds2).reshape(1, -1), (self.N+1, 1))
            gamma_mat = d2f/ds2
            S_arr = S_arr[1:-1]
            # 选取Smin到Smax的值
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            gamma_mat = gamma_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, gamma_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (gamma_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, gamma_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('gamma(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球gamma_3D.png')
            plt.show()
        else:
            # 绘制gamma折线图
            fvalue_arr = f_mat[0, :]
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            d2f = fvalue_arr[2:] + fvalue_arr[:-2] - 2*fvalue_arr[1:-1]
            ds2 = ((S_arr[2:]-S_arr[:-2])/2)**2
            gamma_arr = d2f/ds2
            S_arr = S_arr[1:-1]
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            gamma_arr = gamma_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('gamma(FDM)')
            plt.plot(S_arr, gamma_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.yticks(np.arange(min(gamma_arr)//50*50, (max(gamma_arr)//50+1)*50, 50))
            plt.savefig('FDM下的雪球gamma.png')
            plt.show()
        return

    def plot_snowball_theta(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算theta矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            theta_mat = (f_mat[1:, :] - f_mat[:-1, :]) / np.tile((T_arr[1:] - T_arr[:-1]).reshape(-1, 1), (1, self.M+1))
            T_arr = T_arr[1:]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            theta_mat = theta_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, theta_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (theta_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, theta_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('theta(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            ax.set_zlim3d(np.quantile(theta_mat, 0.01), np.quantile(theta_mat, 0.99))
            plt.savefig('FDM下的雪球theta_3D.png')
            plt.show()
        else:
            # 绘制theta折线图
            theta_arr = (f_mat[1, :] - f_mat[0, :]) / self.delta_T
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            theta_arr = theta_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('theta(FDM)')
            plt.plot(S_arr, theta_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球theta.png')
            plt.show()
        return

    def plot_snowball_vega(self, plot_3d=False):
        start_sigma = 0.01
        end_sigma = 0.4
        delta_sigma = 0.04
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        sigma_arr = np.arange(start_sigma, end_sigma, delta_sigma)
        f_mat = np.zeros((len(sigma_arr), self.M+1))
        for i in tqdm(range(len(sigma_arr))):
            temp_f_mat = self.get_snowball_f_mat(sigma=sigma_arr[i])
            f_mat[i, :] = temp_f_mat[0, :]
        if plot_3d:
            # 计算vega矩阵
            vega_mat = (f_mat[2:, :] - f_mat[:-2, :]) / np.tile((sigma_arr[2:] - sigma_arr[:-2]).reshape(-1, 1), (1, self.M+1))
            sigma_arr = sigma_arr[1:-1]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            vega_mat = vega_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            sigma_mat = np.tile(sigma_arr.reshape(-1, 1), (1, vega_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (vega_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(sigma_mat, S_mat, vega_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('vega(FDM)')
            ax.set_xlabel('sigma')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球vega_3D.png')
            plt.show()
        else:
            # 绘制vega折线图
            idx = int((self.sigma - start_sigma)/delta_sigma)
            vega_arr = (f_mat[idx+1, :] - f_mat[idx-1, :]) / (2*delta_sigma)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            vega_arr = vega_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('vega(FDM)')
            plt.plot(S_arr, vega_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球vega.png')
            plt.show()
        return

    def plot_snowball_rho(self, plot_3d=False):
        start_r = 0.005
        end_r = 0.05
        delta_r = 0.005
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        r_arr = np.arange(start_r, end_r, delta_r)
        f_mat = np.zeros((len(r_arr), self.M+1))
        for i in tqdm(range(len(r_arr))):
            temp_f_mat = self.get_snowball_f_mat(r=r_arr[i])
            f_mat[i, :] = temp_f_mat[0, :]
        if plot_3d:
            # 计算rho矩阵
            rho_mat = (f_mat[2:, :] - f_mat[:-2, :]) / np.tile((r_arr[2:] - r_arr[:-2]).reshape(-1, 1), (1, self.M+1))
            r_arr = r_arr[1:-1]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            rho_mat = rho_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            r_mat = np.tile(r_arr.reshape(-1, 1), (1, rho_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (rho_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(r_mat, S_mat, rho_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('rho(FDM)')
            ax.set_xlabel('r')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球rho_3D.png')
            plt.show()
        else:
            # 绘制rho折线图
            idx = int((self.r - start_r)/delta_r)
            rho_arr = (f_mat[idx+1, :] - f_mat[idx-1, :]) / (2*delta_r)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            rho_arr = rho_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('rho(FDM)')
            plt.plot(S_arr, rho_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球rho.png')
            plt.show()
        return


if __name__ == '__main__':
    fdm = SnowballFDM()
    print('FDM下的雪球结构定价:', fdm.get_snowball_fvalue())
    plot_3d = True
    fdm.plot_snowball_delta(plot_3d=plot_3d)
    fdm.plot_snowball_gamma(plot_3d=plot_3d)
    fdm.plot_snowball_theta(plot_3d=plot_3d)
    fdm.plot_snowball_vega(plot_3d=plot_3d)
    fdm.plot_snowball_rho(plot_3d=plot_3d)

3_雪球有限差分定价/snowball_fdm.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: FDM雪球定价
Date: 2023-08-28 22:54:39
'''
import numpy as np
import pandas as pd
from tqdm import tqdm
from basic_fdm import BasicFDM
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号


class AutocallFDM(BasicFDM):
    '''
    Autocall期权定价
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min时，期权价格为0
        self.f_mat[:, 0] = 0
        # 边界条件2：到期时间T时，在敲出部分获得收益
        idx = int((self.barrier_ko - self.S_min)/self.delta_S)
        self.f_mat[self.N, 0:idx+1] = 0
        self.f_mat[self.N, idx+1:] = self.c * self.T
        # 边界条件3：敲出时才会获得收益
        for i in range(21, self.N+1, 21):
            # 敲出日
            self.f_mat[i, idx:self.M+1] = self.c*i*self.delta_T
            # 问题
            t_arr = np.arange(i-21, i)*self.delta_T
            # 敲出部分: 为下期收益的贴现
            self.f_mat[i-21:i, self.M] = self.c * (i * self.delta_T) * np.exp(-self.r*(i*self.delta_T-t_arr))
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算相应的系数，按照时间进行倒推
        self.calc_coef()
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            if ((i - 1) % 21 == 0) & (i != 1):
                idx = int((self.barrier_ko - self.S_min)/self.delta_S)
                self.f_mat[i - 1, 1:idx+1] = fi_1[0:idx]
            else:
                # 赋值给下个fi
                self.f_mat[i - 1, 1:self.M] = fi_1
        return


class DNTFDM(BasicFDM):
    '''
    双边触碰失效期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 barrier_ki,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        barrier_ki: float. 敲入价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ki = barrier_ki
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min或者敲入时，期权价格为0
        idx_in = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        self.f_mat[:, 0] = 0
        self.f_mat[:, 0:idx_in+1] = 0
        # 边界条件2：到期时间T时，仅在非敲入敲出部分获得全部票息收益
        self.f_mat[self.N, idx_in+1:idx_out+1] = self.c * self.T
        self.f_mat[self.N, 0:idx_in+1] = 0
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者敲出时，期权价格为0
        for i in range(21, self.N+1, 21):
            # 股票价格为Smax
            self.f_mat[i-21:i, self.M] = 0
            # 敲出部分
            self.f_mat[i, idx_out+1:self.M+1] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_ki = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_ko = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            # 如果是敲出日，范围应该是敲入到敲出
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, idx_ki+1:idx_ko+1] = fi_1[idx_ki:idx_ko]
            else:
                # 非敲出日范围应该是敲入到最后
                self.f_mat[i - 1, idx_ki+1:self.M] = fi_1[idx_ki:self.M-1]
        return


class UOPFDM(BasicFDM):
    '''
    上涨失效看跌期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        # 边界条件1： 股票价格为S_min时，不会敲出，且看跌期权必定执行
        ti = np.linspace(0, self.T, self.N+1)
        self.f_mat[:, 0] = self.S0*np.exp(-self.r*(self.T-ti)) - self.S_min
        # 边界条件2：到期时间T时，不敲出收益为max((S0-ST)/S0, 0),敲出收益为0
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        si = np.linspace(self.S_min, self.S_max, self.M+1)
        self.f_mat[self.N, 0:idx_out + 1] = \
            np.maximum((self.S0 - si)/self.S0, 0)[0:idx_out+1]
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者敲出时，收益为0
        for i in range(21, self.N+1, 21):
            self.f_mat[i, idx_out+1:] = 0
            self.f_mat[i-21:i, self.M] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, 1:idx_out+1] = fi_1.reshape(-1)[0:idx_out]
            else:
                self.f_mat[i - 1, 1:self.M] = fi_1.reshape(-1)
        return


class DKOPFDM(BasicFDM):
    '''
    双边失效看跌期权有限差分代码
    '''

    def __init__(self,
                 S0,
                 barrier_ko,
                 barrier_ki,
                 c,
                 r,
                 q,
                 T,
                 sigma,
                 S_min,
                 S_max,
                 M,
                 N,
                 method
                 ):
        '''
        Description
        ----------
        构造函数，输入期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        barrier_ko: float. 敲出价格
        barrier_ki: float. 敲入价格
        c: float. 票面利率
        r: float. 无风险利率
        q: float. 红利率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson
        '''
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N)
        self.c = c
        self.barrier_ki = barrier_ki
        self.barrier_ko = barrier_ko

    def set_boundary(self):
        idx_in = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_out = int((self.barrier_ko - self.S_min)/self.delta_S)
        # 边界条件1： 股票价格为S_min或者向下敲出时，期权价格为0
        self.f_mat[:, 0] = 0
        self.f_mat[:, 0:idx_in+1] = 0
        # 边界条件2：到期时间T时，仅在非敲入敲出部分获得看跌期权收益
        si = np.linspace(self.S_min, self.S_max, self.M+1)
        vanillaT = np.maximum((self.S0 - si)/self.S0, 0)
        self.f_mat[self.N, idx_in+1:idx_out+1] = vanillaT[idx_in+1:idx_out+1]
        self.f_mat[self.N, 0:idx_in+1] = 0
        self.f_mat[self.N, idx_out+1:] = 0
        # 边界条件3: 股票价格为Smax或者向上敲出时，期权价格为0
        for i in range(21, self.N+1, 21):
            # 股票价格为Smax
            self.f_mat[i-21:i, self.M] = 0
            # 敲出部分
            self.f_mat[i, idx_out+1:] = 0
        self._flag_set_boundary = True
        return

    def calc_f_mat(self):
        # 初始化网格矩阵
        if not self._flag_set_boundary:
            raise ValueError('please set f_mat boundary first')
        # 计算系数矩阵
        self.calc_coef()
        idx_ki = int((self.barrier_ki - self.S_min)/self.delta_S)
        idx_ko = int((self.barrier_ko - self.S_min)/self.delta_S)
        for i in tqdm(range(self.N, 0, -1), leave=False):
            fi = (self.f_mat[i, 1:self.M]).reshape(-1, 1)
            # 计算fi-1
            fi_1 = self.get_fi_1(i)
            # 如果是敲出日，范围应该是敲入到敲出
            if ((i - 1) % 21 == 0) & (i != 1):
                self.f_mat[i - 1, idx_ki+1:idx_ko+1] = fi_1[idx_ki:idx_ko]
            else:
                # 非敲出日范围应该是敲入到最后
                self.f_mat[i - 1, idx_ki+1:self.M] = fi_1[idx_ki:self.M-1]
        return


class SnowballFDM(object):
    def __init__(self,
                 S0=1,
                 barrier_ko=1.03,
                 barrier_ki=0.85,
                 c=0.2,
                 r=0.03,
                 q=0,
                 T=1,
                 sigma=0.13,
                 S_min=0,
                 S_max=4,
                 M=1000,
                 N=252,
                 method='Crank-Nicolson'
                 ):
        self.S0 = S0
        self.barrier_ko = barrier_ko
        self.barrier_ki = barrier_ki
        self.c = c
        self.r = r
        self.q = q
        self.T = T
        self.sigma = sigma
        self.S_min = S_min
        self.S_max = S_max
        self.M = M
        self.N = N
        self.delta_S = (S_max - S_min) / M
        self.delta_T = T / N
        self.method = method

    def get_snowball_fvalue(self):
        fdm1 = AutocallFDM(self.S0, self.barrier_ko, self.c, self.r, self.q,
                           self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm2 = DNTFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, self.r, self.q,
                      self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm3 = UOPFDM(self.S0, self.barrier_ko, self.c, self.r, self.q, self.T,
                      self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm4 = DKOPFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, self.r, self.q,
                       self.T, self.sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        f1 = fdm1.get_f_value()
        f2 = fdm2.get_f_value()
        f3 = fdm3.get_f_value()
        f4 = fdm4.get_f_value()
        fvalue = f1 + f2 - f3 + f4
        return fvalue

    def get_snowball_f_mat(self, sigma=None, r=None):
        if sigma is None:
            sigma = self.sigma
        if r is None:
            r = self.r
        fdm1 = AutocallFDM(self.S0, self.barrier_ko, self.c, r, self.q,
                           self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm2 = DNTFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, r, self.q,
                      self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm3 = UOPFDM(self.S0, self.barrier_ko, self.c, r, self.q, self.T,
                      sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm4 = DKOPFDM(self.S0, self.barrier_ko, self.barrier_ki, self.c, r, self.q,
                       self.T, sigma, self.S_min, self.S_max, self.M, self.N, self.method)
        fdm1.set_boundary()
        fdm2.set_boundary()
        fdm3.set_boundary()
        fdm4.set_boundary()
        fdm1.calc_f_mat()
        fdm2.calc_f_mat()
        fdm3.calc_f_mat()
        fdm4.calc_f_mat()
        f_mat = fdm1.f_mat + fdm2.f_mat - fdm3.f_mat + fdm4.f_mat
        return f_mat

    def plot_snowball_delta(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算delta矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            delta_mat = (f_mat[:, 2:] - f_mat[:, :-2]) / np.tile((S_arr[2:] - S_arr[:-2]).reshape(1, -1), (self.N+1, 1))
            S_arr = S_arr[1:-1]
            # 选取Smin到Smax的值
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            delta_mat = delta_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, delta_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (delta_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, delta_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('delta(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球delta_3D.png')
            plt.show()
        else:
            fvalue_arr = f_mat[0, :]
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # delta_arr对应的是Smin+deltaS到Smax-deltaS之间的值
            delta_arr = (fvalue_arr[2:] - fvalue_arr[:-2])/(S_arr[2:]-S_arr[:-2])
            S_arr = S_arr[1:-1]
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            delta_arr = delta_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('delta(FDM)')
            plt.plot(S_arr, delta_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球delta.png')
            plt.show()
        return

    def plot_snowball_gamma(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算gamma矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            d2f = f_mat[:, 2:] + f_mat[:, :-2] - 2*f_mat[:, 1:-1]
            ds2 = ((S_arr[2:]-S_arr[:-2])/2)**2
            ds2 = np.tile((ds2).reshape(1, -1), (self.N+1, 1))
            gamma_mat = d2f/ds2
            S_arr = S_arr[1:-1]
            # 选取Smin到Smax的值
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            gamma_mat = gamma_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, gamma_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (gamma_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, gamma_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('gamma(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球gamma_3D.png')
            plt.show()
        else:
            # 绘制gamma折线图
            fvalue_arr = f_mat[0, :]
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            d2f = fvalue_arr[2:] + fvalue_arr[:-2] - 2*fvalue_arr[1:-1]
            ds2 = ((S_arr[2:]-S_arr[:-2])/2)**2
            gamma_arr = d2f/ds2
            S_arr = S_arr[1:-1]
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            gamma_arr = gamma_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('gamma(FDM)')
            plt.plot(S_arr, gamma_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.yticks(np.arange(min(gamma_arr)//50*50, (max(gamma_arr)//50+1)*50, 50))
            plt.savefig('FDM下的雪球gamma.png')
            plt.show()
        return

    def plot_snowball_theta(self, plot_3d=False):
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        # 计算期权矩阵
        f_mat = self.get_snowball_f_mat()
        if plot_3d:
            # 计算theta矩阵
            T_arr = np.linspace(0, self.T, self.N+1)
            theta_mat = (f_mat[1:, :] - f_mat[:-1, :]) / np.tile((T_arr[1:] - T_arr[:-1]).reshape(-1, 1), (1, self.M+1))
            T_arr = T_arr[1:]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            theta_mat = theta_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            T_mat = np.tile(T_arr.reshape(-1, 1), (1, theta_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (theta_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(T_mat, S_mat, theta_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('theta(FDM)')
            ax.set_xlabel('T')
            ax.set_ylabel('S')
            ax.set_zlim3d(np.quantile(theta_mat, 0.01), np.quantile(theta_mat, 0.99))
            plt.savefig('FDM下的雪球theta_3D.png')
            plt.show()
        else:
            # 绘制theta折线图
            theta_arr = (f_mat[1, :] - f_mat[0, :]) / self.delta_T
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            theta_arr = theta_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('theta(FDM)')
            plt.plot(S_arr, theta_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球theta.png')
            plt.show()
        return

    def plot_snowball_vega(self, plot_3d=False):
        start_sigma = 0.01
        end_sigma = 0.4
        delta_sigma = 0.04
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        sigma_arr = np.arange(start_sigma, end_sigma, delta_sigma)
        f_mat = np.zeros((len(sigma_arr), self.M+1))
        for i in tqdm(range(len(sigma_arr))):
            temp_f_mat = self.get_snowball_f_mat(sigma=sigma_arr[i])
            f_mat[i, :] = temp_f_mat[0, :]
        if plot_3d:
            # 计算vega矩阵
            vega_mat = (f_mat[2:, :] - f_mat[:-2, :]) / np.tile((sigma_arr[2:] - sigma_arr[:-2]).reshape(-1, 1), (1, self.M+1))
            sigma_arr = sigma_arr[1:-1]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            vega_mat = vega_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            sigma_mat = np.tile(sigma_arr.reshape(-1, 1), (1, vega_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (vega_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(sigma_mat, S_mat, vega_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('vega(FDM)')
            ax.set_xlabel('sigma')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球vega_3D.png')
            plt.show()
        else:
            # 绘制vega折线图
            idx = int((self.sigma - start_sigma)/delta_sigma)
            vega_arr = (f_mat[idx+1, :] - f_mat[idx-1, :]) / (2*delta_sigma)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            vega_arr = vega_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('vega(FDM)')
            plt.plot(S_arr, vega_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球vega.png')
            plt.show()
        return

    def plot_snowball_rho(self, plot_3d=False):
        start_r = 0.005
        end_r = 0.05
        delta_r = 0.005
        start_price = round(self.barrier_ki*0.7, 1)
        end_price = round(self.barrier_ko*1.5, 1)
        r_arr = np.arange(start_r, end_r, delta_r)
        f_mat = np.zeros((len(r_arr), self.M+1))
        for i in tqdm(range(len(r_arr))):
            temp_f_mat = self.get_snowball_f_mat(r=r_arr[i])
            f_mat[i, :] = temp_f_mat[0, :]
        if plot_3d:
            # 计算rho矩阵
            rho_mat = (f_mat[2:, :] - f_mat[:-2, :]) / np.tile((r_arr[2:] - r_arr[:-2]).reshape(-1, 1), (1, self.M+1))
            r_arr = r_arr[1:-1]
            # 选取Smin到Smax的值
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))[0]
            rho_mat = rho_mat[:, cond]
            S_arr = S_arr[cond]
            # 绘制曲面
            r_mat = np.tile(r_arr.reshape(-1, 1), (1, rho_mat.shape[1]))
            S_mat = np.tile(S_arr.reshape(1, -1), (rho_mat.shape[0], 1))
            fig = plt.figure(figsize=(10, 5))
            ax = fig.add_subplot(projection='3d')
            surf = ax.plot_surface(r_mat, S_mat, rho_mat, cmap='rainbow', linewidth=0, antialiased=True)
            fig.colorbar(surf)
            ax.set_title('rho(FDM)')
            ax.set_xlabel('r')
            ax.set_ylabel('S')
            plt.savefig('FDM下的雪球rho_3D.png')
            plt.show()
        else:
            # 绘制rho折线图
            idx = int((self.r - start_r)/delta_r)
            rho_arr = (f_mat[idx+1, :] - f_mat[idx-1, :]) / (2*delta_r)
            S_arr = np.linspace(self.S_min, self.S_max, self.M + 1)
            # 选出从start_price到end_price之间的delta
            cond = np.where((S_arr >= start_price) & (S_arr <= end_price))
            rho_arr = rho_arr[cond]
            S_arr = S_arr[cond]
            # 绘制曲线
            plt.figure(figsize=(10, 5))
            plt.title('rho(FDM)')
            plt.plot(S_arr, rho_arr)
            plt.xticks(np.arange(start_price, end_price, 0.1))
            plt.savefig('FDM下的雪球rho.png')
            plt.show()
        return


if __name__ == '__main__':
    df = pd.DataFrame({'q': np.arange(0, 0.2, 0.01)})
    for idx in tqdm(range(len(df))):
        q = df.loc[idx, 'q']
        fdm = SnowballFDM(q=q)
        df.loc[idx, 'price'] = fdm.get_snowball_fvalue()
    df.to_excel('红利率变化下的雪球价格.xlsx', index=False)

3_雪球有限差分定价/snowball_fdm_q.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: FDM欧式香草期权定价
Date: 2023-04-18 15:23:33
'''
import numpy as np
from basic_fdm import BasicFDM
from bs import euro_bs_call, euro_bs_put


class VanillaFDM(BasicFDM):
    '''
    欧式期权有限差分代码
    '''

    def __init__(self,
                 S0=1,
                 K=1,
                 r=0.03,
                 q=0,
                 T=1,
                 sigma=0.13,
                 S_min=0,
                 S_max=2,
                 M=100,
                 N=252,
                 method='Crank-Nicolson',
                 opt_type='call'):
        '''
        Description
        ----------
        构造函数，输入欧式期权参数

        Parameters
        ----------
        S0: float. 当前标的价格
        K: float. 执行价格-----用于测试欧式期权定价，后面删除
        r: float. 无风险收益率
        T: float. 到期时间
        sigma: float. 标的波动率
        S_min: float. 网格中最小标的价格
        S_max: float. 网格中最大标的价格
        M: int. 价格维度上网格段数
        N: int. 时间维度上网格段数
        method: str. 定价方法，有Explicit, Implicit, Crank-Nicolson，默认为Crank-Nicolson
        opt_type: str. 是'call'或者'put'，默认为'call'
        '''
        # 期权参数
        self.K = K  # 期权执行价
        super().__init__(S0, r, q, T, sigma, S_max, S_min, M, N, method)
        if opt_type not in ['call', 'put']:
            raise ValueError('opt_type must be call or put!')
        self.opt_type = opt_type

    def set_boundary(self):
        if self.opt_type == 'call':
            self._call_set_boundary()
        else:
            self._put_set_boundary()
        self._flag_set_boundary = True

    def _call_set_boundary(self):
        # 边界条件1： 股票价格为S_min时，期权价格为0
        self.f_mat[:, 0] = 0
        # 边界条件2：到期时间T时，期权价格f=max(ST-K, 0)
        for j in range(self.M + 1):
            ST = self.S_min + j * self.delta_S
            self.f_mat[self.N, j] = np.maximum(ST - self.K, 0)
        # 边界条件3：S=S_max时，call=S_max-Kexp(-r(T-t))
        for i in range(self.N + 1):
            ti = self.T - i * self.delta_T
            self.f_mat[i, self.M] = self.S_max - self.K * np.exp(-self.r * ti)
        self._init_f_mat_flag = True

    def _put_set_boundary(self):
        # 边界条件1： 初始价格为S_min时，期权价格为f=Kexp(-r(T-t))-Smin
        for i in range(self.N + 1):
            ti = self.T - i * self.delta_T
            self.f_mat[i, 0] = self.K * np.exp(-self.r * ti) - self.S_min
        # 边界条件2：到期时间T时，期权价格f=max(K-ST, 0)
        for j in range(self.M + 1):
            ST = self.S_min + j * self.delta_S
            self.f_mat[self.N, j] = np.maximum(self.K - ST, 0)
        # 边界条件3：S=S_max时，put=0
        self.f_mat[:, self.M] = 0
        self._init_f_mat_flag = True


if __name__ == '__main__':
    # 欧式期权参数
    S0 = 1
    K = 1
    r = 0.03
    q = 0
    T = 1
    sigma = 0.2
    S_max = 4
    S_min = 0
    M = 1000
    N = 252
    fdm_call = VanillaFDM(S0=S0,
                          K=K,
                          r=r,
                          q=q,
                          T=T,
                          sigma=sigma,
                          S_max=S_max,
                          S_min=S_min,
                          M=M,
                          N=N,
                          method='Crank-Nicolson',
                          opt_type='call')
    vanilla_call = fdm_call.get_f_value()
    print('有限差分计算call价值', vanilla_call)
    print('BS计算call价值', euro_bs_call(S0, K, r, T, sigma))
    fdm_put = VanillaFDM(S0=S0,
                         K=K,
                         r=r,
                         T=T,
                         sigma=sigma,
                         S_max=S_max,
                         S_min=S_min,
                         M=M,
                         N=N,
                         opt_type='put')
    vanilla_put = fdm_put.get_f_value()
    print('有限差分计算put价值', vanilla_put)
    print('BS计算put价值', euro_bs_put(S0, K, r, T, sigma))

3_雪球有限差分定价/vanilla_fdm.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: MC情景分析
Date: 2023-08-29 08:32:53
'''
import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from functools import wraps
from datetime import datetime
from tqdm import tqdm
mpl.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号
mpl.rcParams['font.sans-serif'] = ['SimHei']   # 用来正常显示中文标签
plt.rcParams['font.sans-serif'] = ['SimHei']   # 用来正常显示中文标签
np.random.seed(0)                              # fix random seed


def timer(func):
    # 计时器，在函数上一行加上@timer即可
    @wraps(func)
    def wrapper(*args, **kwargs):
        print("Started running function: %s" % func.__name__)
        start_time = datetime.now()
        res = func(*args, **kwargs)
        end_time = datetime.now()
        cost_time = (end_time - start_time).total_seconds()
        print("Finished running function: %s" % func.__name__)
        print("----------Time cost: {}s----------".format(cost_time))
        print("")
        return res
    return wrapper


def stock_monte_carlo(start_price, T, r, q, sigma, steps, simulations):
    '''
    Desctription
    ----------
    蒙特卡洛模拟股票价格

    Parameters
    ----------
    start_price : float
        模拟起始股价.
    T : float
        模拟时间长度，如T=1表示持续1年
    r : float
        股票收益率（年化）.
    sigma : float
        股票波动率（年化）.
    steps : float
        往后预测步数.
    simulations : float
        蒙特卡洛模拟次数.

    Returns
    -------
    price : numpy.array
        模拟股价记录矩阵.

    '''

    dt = T/steps
    # Define a price array
    price = np.zeros((steps + 1, simulations))
    price[0] = start_price
    for i in range(1, steps + 1):
        z = np.random.standard_normal(simulations)
        price[i] = price[i - 1] * np.exp(((r-q) - 0.5 * sigma ** 2) * dt + sigma * np.sqrt(dt) * z)
    return price.T


def snowball_valuation(prices, B1, B2, strike_out_payment, div_payment, rf):
    """
    Description
    ----------
    雪球结构定价

    Parameters
    ----------
    prices : numpy.ndarray
        蒙特卡洛所得的一条模拟价格路径.
    B1 : float
        雪球的敲入价格.
    B2 : float
        雪球的敲出价格.
    strike_out_payment : float
        敲出票息.
    div_payment : float
        红利票息.
    rf : float
        无风险收益率（年化）.

    Returns
    -------
    value : float
        结构定价.
    ret : float
        期间收益.
    down_flag : bool
        是否敲入.
    up_flag : bool
        是否敲出.
    end_day : int
        存续天数.

    """

    value = None
    down_flag = False                        # 是否已经触发敲入
    up_flag = False                          # 是否触发敲出
    days = np.arange(0, prices.shape[0], 1)  # 包括起始日一共有253天
    up_B2 = (prices >= B2)                   # 高于敲出价格的日期
    KO_observe = (days % 21 == 0)            # 敲出观察日
    KO_days = days[up_B2*KO_observe]         # 敲出日
    KI_days = days[prices < B1]              # 敲入日
    up_flag = (len(KO_days) > 0)
    down_flag = (len(KI_days) > 0)
    if len(KO_days) > 0:                     # 曾敲出
        end_day = KO_days[0]                 # 敲出日期
        value = strike_out_payment*(end_day/252)*np.exp(-rf*end_day/252)
    elif len(KO_days) == 0 and len(KI_days) > 0:  # 曾敲入
        end_day = 252                        # 持有至到期
        if prices[-1] < prices[0]:             # 如果到期标的小于期初价格
            value = (prices[-1]-prices[0]) * np.exp(-rf*end_day/252)
        else:                                # 如果到期标的高于期初价格，获益为0
            value = 0
    else:  # 未敲入也未敲出
        end_day = 252
        value = div_payment*(end_day/252)*np.exp(-rf*end_day/252)
    return value, down_flag, up_flag, end_day


def main():
    # 用于MC的常数
    steps = 252
    T = 1
    r = 0.03
    q = 0
    sigma = 0.13
    simulation_time = 300000
    # 雪球结构参数
    S0 = 1.0
    B2 = 1.03*S0  # 敲出界限：103% 每个月观察敲出
    B1 = 0.85*S0  # 敲入界限：85% 每天观察敲入
    strike_out_payment = 0.2  # 敲出票息=20%（年化）
    div_payment = 0.2  # 红利票息=20%（年化）
    rf = 0.03
    record = stock_monte_carlo(S0, T, r, q, sigma, steps, simulation_time)
    result_df = pd.DataFrame(index=range(len(record)), columns=["value", "ret", "down_KI", "up_KO", "exist_days"])
    for i in tqdm(range(len(record))):
        prices = record[i]
        value, down_flag, up_flag, exist_days = snowball_valuation(prices, B1, B2, strike_out_payment, div_payment, rf)
        result_df.loc[i]["value"] = value
        result_df.loc[i]["down_KI"] = down_flag
        result_df.loc[i]["up_KO"] = up_flag
        result_df.loc[i]["exist_days"] = exist_days

    snowball_df = pd.DataFrame(index=["价值", "最大损失", "存续交易日数", "概率"],
                               columns=["雪球结构整体", "敲出", "敲入未敲出", "未敲入也未敲出"])
    # 组合价值
    snowball_df.loc["价值"]["雪球结构整体"] = result_df["value"].mean()
    snowball_df.loc["最大损失"]["雪球结构整体"] = result_df["value"].min()
    snowball_df.loc["存续交易日数"]["雪球结构整体"] = result_df["exist_days"].mean()

    # 敲出
    snowball_df.loc["价值"]["敲出"] = result_df.query("up_KO==True")["value"].mean()
    snowball_df.loc["最大损失"]["敲出"] = result_df.query("up_KO==True")["value"].min()
    snowball_df.loc["存续交易日数"]["敲出"] = result_df.query("up_KO==True")["exist_days"].mean()
    snowball_df.loc["概率"]["敲出"] = result_df.query("up_KO==True").shape[0]/result_df.shape[0]
    # 敲入
    snowball_df.loc["价值"]["敲入未敲出"] = result_df.query("down_KI==True and up_KO==False")["value"].mean()
    snowball_df.loc["最大损失"]["敲入未敲出"] = result_df.query("down_KI==True and up_KO==False")["value"].min()
    snowball_df.loc["存续交易日数"]["敲入未敲出"] = result_df.query("down_KI==True and up_KO==False")["exist_days"].mean()
    snowball_df.loc["概率"]["敲入未敲出"] = result_df.query("down_KI==True and up_KO==False").shape[0]/result_df.shape[0]
    # 未敲入也未敲出
    snowball_df.loc["价值"]["未敲入也未敲出"] = result_df.query("down_KI==False and up_KO==False")["value"].mean()
    snowball_df.loc["最大损失"]["未敲入也未敲出"] = result_df.query("down_KI==False and up_KO==False")["value"].min()
    snowball_df.loc["存续交易日数"]["未敲入也未敲出"] = result_df.query("down_KI==False and up_KO==False")["exist_days"].mean()
    snowball_df.loc["概率"]["未敲入也未敲出"] = result_df.query("down_KI==False and up_KO==False").shape[0]/result_df.shape[0]
    snowball_df.to_excel("MC下的雪球情景分析.xlsx")


if __name__ == '__main__':
    main()

4_雪球蒙特卡洛情景分析/snowball_mc_scenario.py
dkl0707/SnowballPricing
'''
Author: dkl
Description: 欧式看跌期权theta
Date: 2023-09-20 11:54:57
'''
import numpy as np
import matplotlib.pyplot as plt
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']
mpl.rcParams['axes.unicode_minus'] = False     # 用来正常显示负号


def theta_option(S, K, sigma, r, T, optype):
    import numpy as np
    from scipy.stats import norm
    d1 = (np.log(S/K)+(r+sigma**2/2)*T)/(sigma*np.sqrt(T))
    d2 = d1 - sigma*np.sqrt(T)
    theta_call = -(S*sigma*np.exp(-d1**2/2)) / \
        (2*np.sqrt(2*np.pi*T)) - r*K*np.exp(-r*T)*norm.cdf(d2)
    if optype == 'call':
        theta = theta_call
    else:
        theta = theta_call + r*K*np.exp(-r*T)
    return theta


# 标的价格与Theta的关系
S_list = np.linspace(0.5, 1.5, 100)
theta_list = theta_option(S=S_list, K=1, sigma=0.13, r=0.03, T=1, optype='put')


plt.figure(figsize=(12, 6))
plt.plot(S_list, -theta_list, label='看跌期权')
plt.grid('True')
plt.savefig('欧式看跌期权空头theta.png')
plt.show()

5_欧式看跌期权theta绘制/欧式看跌期权空头theta.py
dkl0707/SnowballPricing