期权定价Python实现

期权定价公式

期权定价通常使用BSM公式，如下所示：
$$C(S_t,K,t,T,r,\sigma )=S_{t}N(d_1)-e^{-r(T-t)}KN(d_2)\text{\hspace{1em}(1.1)}$$
其中$N(d)$为标准正态分布$x\sim N(0,1)$时$x<d$时的概率，如下所示：
$$N(d)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^d{e}^{-\frac{1}{2}{x}^2}dx\text{\hspace{1em}(1.2)}$$
式中$d_1$和$d_2$为两个常数，定义为：
$$d_1=\frac{\log \frac{S_t}{K}+(r+\frac{{\sigma }^2}{2})(T-t)}{\sigma \sqrt{T-t}}\text{\hspace{1em}(1.3)}$$
$$d_2=\frac{\log \frac{S_t}{K}+(r-\frac{{\sigma }^2}{2})(T-t)}{\sigma \sqrt{T-t}}\text{\hspace{1em}(1.4)}$$
上式中的符号说明如下所示：

• C：期权价格；
• $S_t$：标的资产价格；
• K：期权行权价格；
• t：当前时间；
• T：期权到期时间；
• r：短期无风险利率；
• $\sigma$：波动率；

迭代法

我们在初始状态下，先设定$\sigma$为一个随机值${\sigma }_0^{imp}$，其中上标imp代表是隐含波动率，下标代表是第0步。将${\sigma }_0^{imp}$代入上面的公式，就可以求出一个期权的计算值$\hat{C}$，其与当前市场上的期权价格$C$的差值可以用来代表我们当初随机取的${\sigma }_0^{imp}$接近真实值的程度。
我们可以通过牛顿法来求解，我们首先求出$\frac{\partial C}{\partial \sigma }$：
$$vega=\frac{\partial C}{\partial \sigma }=S_t{N}^{\prime }(d_1))\sqrt{T-t}\text{\hspace{1em}(1.5)}$$
我们可以通过迭代法来找到更好的$\sigma$的值：
$${\sigma }_{n+1}^{imp}={\sigma }_n^{imp}-\frac{1}{vega}({\hat{C}}_n-C)\text{\hspace{1em}(1.6)}$$
上式中$n$代表迭代步数，C代表期权的真实价格，${\hat{C}}_n$为当前的计算值。

代码实现

我们首先定义期权估值的BS模型类，如下所示：

#
import math
from enum import Enum,unique
from scipy.stats import norm

@unique
class DMode(Enum):
    D1 = 1
    D2 = 2

class BsModel(object):
    def __init__(self):
        self.name = 'fak.option.bs_model.BsModel'
    
    @staticmethod
    def calculate_C(St, K, t, T, r, sigma):
        '''
        已知波动率计算期权价格的BS公式
        '''
        d1 = BsModel.calculate_d1(St, K, t, T, r, sigma)
        d2 = BsModel.calculate_d2(St, K, t, T, r, sigma)
        nd1 = BsModel.calculate_nd(d1)
        nd2 = BsModel.calculate_nd(d2)
        return St * nd1 - math.exp(-r*(T-t)) * K * nd2

    @staticmethod
    def calculate_vega(St, K, t, T, r, sigma):
        '''
        计算期权价格对d1的导数，即希腊字母vega的值
        '''
        d1 = BsModel.calculate_d1(St, K, t, T, r, sigma)
        nd1p = BsModel.calculate_ndp(d1)
        return St * nd1p * math.sqrt(T-t)

    @staticmethod
    def calculate_d1(St, K, t, T, r, sigma):
        return BsModel.calculate_d(St, K, t, T, r, sigma, DMode.D1)

    @staticmethod
    def calculate_d2(St, K, t, T, r, sigma):
        return BsModel.calculate_d(St, K, t, T, r, sigma, DMode.D2)

    @staticmethod
    def calculate_d(St, K, t, T, r, sigma, mode):
        if DMode.D1 == mode:
            numerator = math.log(St / K)  + (r + sigma*sigma / 2.0)*(T - t)
        else:
            numerator = math.log(St / K) + (r - sigma*sigma / 2.0)*(T - t)
        denominator = sigma * math.sqrt(T-t)
        return numerator / denominator
        
    @staticmethod
    def calculate_nd(x0, mu=0.0, sigma=1.0):
        '''
        计算由均值为mu，方差为sigma的正态分布下，x<=x0时的概率
        '''
        return norm.cdf(x0, loc=mu, scale=sigma)

    @staticmethod
    def calculate_ndp(x0, mu=0.0, sigma=1.0):
        return norm.pdf(x0, loc=mu, scale=sigma)

代码解读如下所示：

• 第6~9行：定义枚举类来表示BS公式中计算d1还是d2；
• 第16~24行：根据式$(1.1)$计算给定波动率条件下期权价格；
• 第27~33行：根据式$(1.5)$计算期权价格对波动率的导数，即希腊字母vega；
• 第36~50行：计算BS公式中的d1和d2的值，通过mode来区分；
• 第53~57行：计算式$(1.1)$中$N(d)$，该项表示正态分布下值小于d的概率，需要用norm.cdf来计算；
• 第60、61行：计算式$(1.5)$中N(d)的导数，我们用norm.pdf来计算；
  接下来，我们先假设我们知道波动率$\sigma$，然后可以得到对应的期权价格$C$，然后我们随机给一个$\sigma$的值，就可以计算出期权价格的计算值$C_t^{imp}$，接下来我们求出vega，然后根据式$(1.6)$求出新的波动率${\sigma }_t^{imp}$，循环上述过程，我们就可以求出波动率的真实值，在本例中就是我们初始时给定的波动率的值，代码如下所示：

    def startup(self):
        # 初始化变量
        St = 2100.0 # 标的资产价格
        K = 2200.0 # 行权价格
        t = 10 # 当前时间
        T = 30 # 到期时间
        r = 0.01 # 短期无风险利率
        sigma = 0.3 # 波动率
        # 假设已知波动率的情况下计算期权价格和vega
        C = BsModel.calculate_C(St, K, t, T, r, sigma)
        vega = BsModel.calculate_vega(St, K, t, T, r, sigma)
        # 随机给出一个波动率，用迭代法求出波动率真值
        sigma_t = 1200
        for n in range(10000):
            C_t = BsModel.calculate_C(St, K, t, T, r, sigma_t)
            vega_t = BsModel.calculate_vega(St, K, t, T, r, sigma_t)
            sigma_t = sigma_t - 1/vega * (C_t - C)
            print('C0={0}; vega={1}; sigma={2};'.format(C_t, vega_t, sigma_t))
        print('真实值：C={0}; vega={1};'.format(C, vega))

运行结果如下所示：

如上图所示，我们经过1000步迭代，$\sigma$从我们最初给定的随机值1200，逐步收敛到我们初始时设定的真实0.3，这正是我们想要的结果。