[Python] オプション価格の計算

ブラックショールズモデルに基づきオプション価格を計算します。

ブラックショールズ方程式の解

コールオプションの価格$C$、プットオプションの価格$P$は以下のようになります。

$${\displaystyle N(x)={\frac {1}{\sqrt {2\pi }}}\int _{-\infty }^{x}e^{-{\frac {y^{2}}{2}}}dy}$$

$${\displaystyle d_{1}={\frac {\log({\frac {S_{t}}{E}})+(r+{\frac {\sigma ^{2}}{2}})(T-t)}{\sigma {\sqrt {T-t}}}}}$$

$${\begin{align} \displaystyle d_{2} &={\frac {\log({\frac {S_{t}}{E}})+(r-{\frac {\sigma ^{2}}{2}})(T-t)}{\sigma {\sqrt {T-t}}}} \\&= d_1- \sigma {\sqrt {T-t}} \end{align}}$$

として、

$$ C = S(0)N(d_1)- E exp(-r(T-t))N(d_2) $$

$$ P = -S(0)N(-d_1) + E exp(-r(T-t))N(-d_2) $$

import numpy as np
from scipy import stats

class BlackScholes(object):

    def __init__(self, stock_price, strike_price, expiry, risk_free_rate, sigma):
        self.s = stock_price
        self.e = strike_price
        self.t = expiry # expiry 1 is 1_year=365_days 
        self.r = risk_free_rate
        self.sigma = sigma # volatility of the given stock
        self.d1 = (np.log(self.s/self.e) + (self.r + self.sigma**2/2)*self.t) / (self.sigma*self.t**(1/2))
        self.d2 = (np.log(self.s/self.e) + (self.r-self.sigma**2/2)*self.t) / (self.sigma*self.t**(1/2))

    def call_price(self):
        call_price = self.s*stats.norm.cdf(self.d1) - self.e*np.exp(-self.r*self.t)*stats.norm.cdf(self.d2) 
        return call_price

    def put_price(self):
        put_price = -self.s*stats.norm.cdf(-self.d1) + self.e*np.exp(-self.r*self.t)*stats.norm.cdf(-self.d2) 
        return put_price

if __name__ == '__main__':
    stock_price = 20000
    strike_price = 20000
    expiry = 1
    risk_free_rate = 0.02
    sigma = 0.2

    black_scholes = BlackScholes(stock_price, strike_price, expiry, risk_free_rate, sigma)
    # 1783.207455714506
    print(black_scholes.call_price())
    # 1387.1809218496146
    print(black_scholes.put_price())

以下のサイトと比較してみます。

ブラックショールズモデル計算機

公式にそのまま当てはめたので当然ながら一致します。

モンテカルロ法のシミュレーション

ブラックショールズモデルにより権利行使日での株価を計算します。

その株価に基づきペイオフを計算し、ペイオフを現在価値へ割り引くことで、現在のオプション価格を計算します。

株価

株価$S(t)$を求めます。

株価の微分方程式は以下のようになります。

$$ d S(t) = \mu S(t) dt + \sigma S(t) dW $$

株価は負にはならないので、$F(S) = \log S(t)$として、伊藤の補題の公式を適用します。

Mathematical formulation of Itô’s lemma

$$ d \log S(t) = (\mu – \frac{1}{2} \sigma^2)dt + \sigma dW$$

$$ log S(t) = log S(0) + (\mu – \frac{1}{2} \sigma^2) t + \sigma \int_0^t dW $$

$dW$はウィーナー過程の増分なので、

$$ \int_0^t dW = N(0, t) = \sqrt{t} N(0, 1) $$

また、以下より$\mu=r$とします。

“The drift of stock price becomes the risk-free interest rate” under RNP

$$ S(t) = S(0) exp ((r – \frac{1}{2}\sigma^2)t + \sigma \sqrt{t} N(0, 1)) $$

オプション価格

行使価格$E$としてペイオフを求め、さらに、 $e^{-rt}$をかけることで現在価値へと割り引き、オプション価格を計算します。

コールオプション

$$ max(S-E, 0) \cdot e^{-rt} $$

プットオプション

$$ max(E-S, 0) \cdot e^{-rt} $$

以上をコードにします。

import numpy as np
from scipy import stats

class BlackScholes(object):

    def __init__(self, stock_price, strike_price, expiry, risk_free_rate, sigma):
        self.s = stock_price
        self.e = strike_price
        self.t = expiry # expiry 1 is 1_year=365_days 
        self.r = risk_free_rate
        self.sigma = sigma # volatility of the given stock
        self.d1 = (np.log(self.s/self.e) + (self.r + self.sigma**2/2)*self.t) / (self.sigma*self.t**(1/2))
        self.d2 = (np.log(self.s/self.e) + (self.r-self.sigma**2/2)*self.t) / (self.sigma*self.t**(1/2))

    def call_price(self):
        call_price = self.s*stats.norm.cdf(self.d1) - self.e*np.exp(-self.r*self.t)*stats.norm.cdf(self.d2) 
        return call_price

    def put_price(self):
        put_price = -self.s*stats.norm.cdf(-self.d1) + self.e*np.exp(-self.r*self.t)*stats.norm.cdf(-self.d2) 
        return put_price

    def monte_carlo_call(self, iterations):
        np.random.seed(seed=0)

        option_data = np.zeros([iterations, 2])
        rands = np.random.normal(0, 1, [1, iterations])
        stock_price = self.s * np.exp(self.t*(self.r-0.5*self.sigma**2) + self.sigma*self.t**(1/2)*rands)
        option_data[:, 1] = stock_price - self.e 

        average = np.sum(np.amax(option_data, axis=1))/iterations

        return average * np.exp(-1*self.r*self.t)

    def monte_carlo_put(self, iterations):
        np.random.seed(seed=0)

        option_data = np.zeros([iterations, 2])
        rands = np.random.normal(0, 1, [1, iterations])
        stock_price = self.s * np.exp(self.t*(self.r-0.5*self.sigma**2) + self.sigma*self.t**(1/2)*rands)
        option_data[:, 1] = self.e - stock_price

        average = np.sum(np.amax(option_data, axis=1))/iterations

        return average * np.exp(-1*self.r*self.t)


if __name__ == '__main__':
    stock_price = 20000
    strike_price = 20000
    expiry = 1
    risk_free_rate = 0.02
    sigma = 0.2
    
    black_scholes = BlackScholes(stock_price, strike_price, expiry, risk_free_rate, sigma)
    # 1783.207455714506
    print(black_scholes.call_price())
    # 1387.1809218496146
    print(black_scholes.put_price())

    iteratoins = 10000000
    # 1783.4595785128108
    print(black_scholes.monte_carlo_call(iteratoins))
    # 1386.4309824340962
    print(black_scholes.monte_carlo_put(iteratoins))

シミュレーションも理論値と同じような値になっています。