ブラックショールズモデルに基づきオプション価格を計算します。
ブラックショールズ方程式の解
ブラック–ショールズ方程式の解をそのまま使います。
コールオプションの価格\(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))
シミュレーションも理論値と同じような値になっています。