ブラックショールズモデルにおけるオプショングリークスの計算式とプログラムを載せるだけ
何故か日本語で詳しく計算を追ってる記事が見つからなかったため

グリークスとは

株式に配当が含まれている場合の、ブラックショールズ方程式におけるヨーロピアンコール/プットオプション価格( $V_c,V_p$ )は現在時刻t、原資産価格 $S_t$ 、ボラティリティ $\sigma$ 、行使価格K、無リスク金利r、配当利回りq、満期Tまでの時間 $\tau (= T- t)$ によって決まる。具体的には以下の式

$\begin{align} V_c &= S_t e^{-q \tau} N(d_1) - K e^{-r \tau} N(d_2) \\ V_p &= -S_t e^{-q \tau} N(-d_1) + K e^{-r \tau} N(-d_2) \\ \end{align}$

ただし、

$\begin{align} N(x) &= \int^x_{-\infty} \phi(y) dy \\ \phi(y) &= \frac{1}{\sqrt{2 \pi}} e^{-\frac{y^2}{2}} \\ d_{1,2} &= \frac{log(\frac{S_t}{K}) + (r - q \pm \frac{\sigma^2}{2} \tau)}{\sigma \sqrt{\tau}} \end{align}$

グリークスとはこの式を各変数で偏微分したものであり、つまりは各変数に対するオプション価格の感応度を表したものである。具体的には以下の通り¹

$\begin{align} \Delta &= \frac{\partial V}{\partial S_t} \\ \gamma &= \frac{\partial^2 V}{\partial S_t^2} \\ \nu &= \frac{\partial V}{\partial \sigma}\\ \rho &= \frac{\partial V}{\partial r}\\ \epsilon &= \frac{\partial V}{\partial q}\\ \theta &= \frac{\partial V}{\partial t}\\ \end{align}$

計算

準備

この先度々出てくる式を先に計算する。

１.（微積分の交換）

$\begin{align} \frac{\partial N(x)}{\partial x} &= \frac{\partial}{\partial x} \int^{x}_{- \infty} \phi (z) dz \\ &= \phi(x) \tag{1} \end{align}$

２.

$\begin{align} \phi(d_2) &= \phi(d_1 - \sigma \sqrt{\tau}) \\ &= \frac{1}{\sqrt{2 \pi}} e^{-\frac{(d_1 - \sigma \tau)^2}{2}}\\ &= \left( \frac{1}{\sqrt{2 \pi}} e^{-\frac{d^2_1}{2}} \right) e^{ d_1 \sigma \sqrt{\tau} - \frac{\sigma^2 \tau}{2}}\\ &= \phi(d_1) exp [\left( log(\frac{S_t}{K}) + (r - q + \frac{\sigma^2}{2}) \tau \right) - \frac{\sigma^2 \tau}{2} ] \\ &= \phi(d_1) \frac{S_t}{K} e^{(r - q) \tau} \tag{2} \end{align}$

３.

$\begin{align} d_1 - d_2 &= \frac{log(\frac{S_t}{K}) + (r - q + \frac{\sigma^2}{2} \tau)}{\sigma \sqrt{\tau}} - \frac{log(\frac{S_t}{K}) + (r - q - \frac{\sigma^2}{2} \tau)}{\sigma \sqrt{\tau}} \\ &= \sigma \sqrt{\tau} \tag{3} \end{align}$

４.（プットコールパリティ）

$\begin{align} V_C - V_p &= S_t e^{-q \tau} N(d_1) - K e^{-r \tau} N(d_2) - (- S_t e^{-q \tau} N(-d_1) + K e^{-r \tau} N(-d_2))\\ &= S_t e^{-q \tau} (N(d_1) + N(-d_1)) - K e^{-r \tau} (N(d_2) + N(-d_2))\\ &= S_t e^{-q \tau} - K e^{-r \tau} \tag{4} \end{align}$

デルタ（Δ）

コール/プットオプションのデルタをそれぞれ $\Delta_c, \Delta_p$ と表す。

コールのデルタ

$\begin{align} \Delta_c &= \frac{\partial V_c}{\partial S_t} \\ &= e^{-q \tau} N(d_1) + S_t e^{-q \tau} \frac{\partial N(d_1)}{\partial S_t}- Ke^{-r \tau} \frac{\partial N(d_2)}{\partial S_t}\\ &= e^{-q \tau} N(d_1) + S_t e^{-q \tau} \frac{\partial N(d_1)}{\partial d_1} \frac{\partial d_1}{\partial S_t}- Ke^{-r \tau} \frac{\partial N(d_2)}{\partial d_2} \frac{\partial d_2}{\partial S_t}\\ \end{align}$

ここで、

$\begin{align} \frac{\partial d_1}{\partial S_t} = \frac{\partial d_2}{\partial S_t} \end{align}$

と(1)式より

$\begin{align} \Delta_c &= e^{-q \tau} N(d_1) + \frac{\partial d_1}{\partial S_t} (S_t e^{-q \tau} \phi (d_1)- Ke^{-r \tau} \phi (d_2))\\ \end{align}$

が得られる。更に二項目は(2)式を用いることで、

$\begin{align} \Delta_c &= S_t e^{-q \tau} \phi (d_1)- Ke^{-r \tau} \phi (d_2)\\ &= S_t e^{-q \tau} \phi (d_1)- Ke^{-r \tau} \left( \phi (d_1) \frac{S_t}{K} e^{(r - q) \tau} \right) \\ &= \phi(d_1) \left( S_t e^{-q \tau}- S_t e^{-q \tau} \right) \\ &= 0 \end{align}$

であることから０である。従って、

$\begin{align} \Delta_c &= e^{-q \tau} N(d_1) \end{align}$

プットのデルタ

(4)式より、

$\begin{align} \frac{\partial V_c}{\partial S_t} - \frac{\partial V_p}{\partial S_t} &= e^{ - qt} \end{align}$

よって、

$\begin{align} \Delta_p &= \Delta_c - e^{ - qt} \\ &= e^{-qt}(N(d_1) - 1) \\ &= -e^{-qt}(N(- d_1)) \end{align}$

ガンマ（γ）

コール/プットオプションのデルタをそれぞれ $\gamma_c, \gamma_p$ と表す。

コールのガンマ

デルタの計算結果を用いて以下のように計算することが出来る。

$\begin{align} \gamma_c &= \frac{\partial^2 V_c}{\partial S_t^2}\\ &= \frac{\partial \Delta_c}{\partial S_t}\\ &= e^{-q\tau} \frac{\partial N(d_1)} {\partial S_t}\\ &= e^{-q\tau} \frac{\partial N(d_1)} {\partial d_1} \frac{\partial d_1} {\partial S_t}\\ \end{align}$

(1)式より、

$\begin{align} \gamma_c &= e^{-q\tau} \phi(d_1) \left( \frac{1}{\sigma \sqrt{\tau}} \frac{1}{S_t} \right) \\ &= \frac{1} {\sigma S_t \sqrt{2 \pi \tau}} e^{- \frac{d_1^2}{2} - q \tau} \end{align}$

プットのガンマ

(4)式より、

$\begin{align} \frac{\partial^2 V_c}{\partial S_t^2} - \frac{\partial^2 V_p}{\partial S_t^2} &= 0 \end{align}$

よって、

$\begin{align} \gamma_p &= \gamma_c \\ &= \frac{1} {\sigma S_t \sqrt{2 \pi \tau}} e^{- \frac{d_1^2}{2} - q \tau} \end{align}$

ベガ（ν）

コール/プットオプションのベガをそれぞれ $\nu_c, \nu_p$ と表す。

コールのベガ

$\begin{align} \nu_c &= \frac{\partial V_c}{\partial \sigma} \\ &= S_t e^{-qt} \frac{\partial N(d_1)}{\partial \sigma} - K e^{-r \tau} \frac{\partial N(d_1)}{\partial \sigma}\\ &= S_t e^{-qt} \frac{\partial N(d_1)}{\partial d_1} \frac{\partial d_1}{\partial \sigma} - K e^{-r \tau} \frac{\partial N(d_1)}{\partial d_2} \frac{\partial d_2}{\partial \sigma}\\ \end{align}$

ここで、(1)式より

$\begin{align} \nu_c &= S e^{-qt} \phi(d_1) \frac{\partial d_1}{\partial \sigma} - K e^{-r \tau} \phi(d_2) \frac{\partial d_2}{\partial \sigma}\\ \end{align}$

また、(3)式より

$\begin{align} \frac{\partial d_1}{\partial \sigma} - \frac{\partial d_2}{\partial \sigma} &= \sqrt{\tau}\\ \frac{\partial d_1}{\partial \sigma} &= \frac{\partial d_2}{\partial \sigma} + \sqrt{\tau}\\ \end{align}$

であるから

$\begin{align} \nu_c &= S_t e^{-qt} \phi(d_1) (\frac{\partial d_2}{\partial \sigma} + \sqrt{\tau})- K e^{-r \tau} \phi(d_2) \frac{\partial d_2}{\partial \sigma}\\ &= \frac{\partial d_2}{\partial \sigma} (S_t e^{-qt} \phi(d_1)- K e^{-r \tau} \phi(d_2) ) + S_t e^{-qt} \phi(d_1) \sqrt{\tau}\\ \end{align}$

ここで、コールのガンマを求める際と同様に、第一項は0になる。従って

$\begin{align} \nu_c & = S_t e^{-qt} \phi(d_1) \sqrt{\tau}\\ &= \left( S_t \sqrt{ \frac{\tau}{2 \pi} } e^{-\frac{d^2_1}{2} - q\tau} \right) \end{align}$

プットのベガ

(4)式より、

$\begin{align} \frac{\partial V_c}{\partial \sigma} - \frac{\partial V_p}{\partial \sigma} &= 0\\ \nu_c &= \nu_p \end{align}$

従って、

$\begin{align} \nu_p &= S_t e^{-qt} \phi(d_1) \sqrt{\tau}\\ &= \left( S_t \sqrt{ \frac{\tau}{2 \pi} } e^{-\frac{d^2_1}{2} - q\tau} \right) \end{align}$

ロー（ρ）

コールのロー

$\begin{align} \rho_c &= \frac{\partial C}{\partial r} \\ &= S_t e^{-q\tau} \frac{\partial N(d_1)}{\partial r} - K(-\tau) e^{-r \tau} N(d_2) - Ke^{-r \tau} \frac{\partial N(d_2)}{\partial r} \\ &= K \tau e^{-r \tau} N(d_2) + S_t e^{-q\tau} \frac{\partial N(d_1)}{\partial d_1} \frac{\partial d_1}{\partial r} - Ke^{-r \tau} \frac{\partial N(d_2)}{\partial d_2} \frac{\partial d_2}{\partial r} \end{align}$

ここで、

$\begin{align} \frac{\partial d_1}{\partial r} &= \frac{\partial d_2}{\partial r} \end{align}$

と(1)式より

$\begin{align} \rho_c &= K\tau e^{-r \tau} N(d_2) + \frac{\partial d_1}{\partial r} \left( S_t e^{-q\tau} \phi(d_1) - Ke^{-r \tau} \phi(d_2) \right) \\ \end{align}$

コールのガンマを求める際と同様に、第二項は０であるので

$\begin{align} \rho_c &= K\tau e^{-r \tau} N(d_2) \end{align}$

プットのロー

(4)式より、

$\begin{align} \frac{\partial V_c}{\partial r} - \frac{\partial V_p}{\partial r} &= -\frac{\partial K e^{-r \tau}}{\partial r}\\ \rho_p = \rho_c \end{align}$

従って、

$\begin{align} \rho_p &= K\tau e^{-r \tau} N(d_2) \end{align}$

イプシロン（ε）

コールのイプシロン

$\begin{align} \epsilon_c &= \frac{\partial C}{\partial q} \\ &= S_t (-\tau) e^{-q\tau} N(d_1) +S_t e^{-q\tau} \frac{\partial N(d_1)}{\partial q} - Ke^{-r \tau} \frac{\partial N(d_2)}{\partial q} \\ &= -S_t \tau e^{-q\tau} N(d_1) +S_t e^{-q\tau} \frac{\partial N(d_1)}{\partial d_1} \frac{\partial d_1}{\partial q} - Ke^{-r \tau} \frac{\partial N(d_2)}{\partial d_2} \frac{\partial d_2}{\partial q} \\ \end{align}$

ここで、

$\begin{align} \frac{\partial d_1}{\partial q} &= \frac{\partial d_2}{\partial q} \end{align}$

と(1)式より

$\begin{align} \epsilon_c &= -S_t \tau e^{-q \tau} N(d_1) +\frac{\partial d_2}{\partial q} \left( S_t e^{-q \tau} \phi (d_1) - Ke^{-r \tau} \phi (d_2) \right) \\ &= -S_t \tau e^{-q \tau} N(d_1) +\frac{\partial d_2}{\partial q} \left( S_t e^{-q \tau} \phi (d_1) - Ke^{-r \tau} \phi (d_2) \right) \\ \end{align}$

コールのガンマを求める際と同様に、第二項は０であるので

$\begin{align} \epsilon_c &= -S_t \tau e^{-q \tau} N(d_1) \end{align}$

プットのイプシロン

(4)式より

$\begin{align} \frac{\partial V_c}{\partial q} - \frac{\partial V_p}{\partial q} &= -\frac{\partial S_t e^{-q \tau}}{\partial q}\\ \end{align}$

従って、

$\begin{align} \epsilon_p &= \epsilon_c + S_t \tau e^{-q \tau} \\ &= -S_t \tau e^{-q \tau} N(d_1) + S_t \tau e^{-q \tau} \\ &= S_t \tau e^{-q \tau} \left( 1 - N(d_1) \right) \\ &= S_t \tau e^{-q \tau} N(- d_1) \end{align}$

シータ（Θ）

コールのシータ

$\begin{align} \theta_c &= \frac{\partial V_c}{\partial t}\\ &= S_t \frac{\partial e^{-q \tau} }{\partial t} N(d_1) + S_t e^{-q \tau} \frac{\partial N(d_1)}{\partial t} - K \frac{\partial e^{-r \tau} }{\partial t} N(d_2) - K e^{-r \tau} \frac{\partial N(d_2)}{\partial t}\\ &= S_t \frac{\partial e^{-q \tau} }{\partial \tau} \frac{\partial \tau }{\partial t} N(d_1) + S_t e^{-q \tau} \frac{\partial N(d_1)}{\partial d_1} \frac{\partial d_1}{\partial \tau} \frac{\partial \tau}{\partial t} - K \frac{\partial e^{-r \tau} }{\partial \tau} \frac{\partial \tau}{\partial t} N(d_2) - K e^{-r \tau} \frac{\partial N(d_2)}{\partial d_2} \frac{\partial d_2}{\partial \tau} \frac{\partial \tau}{\partial t}\\ \end{align}$

(1)式より

$\begin{align} \theta_c &= S_t (-q e^{-q \tau}) (-1) N(d_1) + S_t e^{-q \tau} \phi(d_1) \frac{\partial d_1}{\partial \tau} (-1) - K (-r e^{-r \tau} ) (-1) N(d_2) - K e^{-r \tau} \phi(d_2) \frac{\partial d_2}{\partial \tau} (-1) \\ &= q S_t e^{-q \tau} N(d_1) - S_t e^{-q \tau} \phi(d_1) \frac{\partial d_1}{\partial \tau} -r K e^{-r \tau} N(d_2) - K e^{-r \tau} \phi(d_2) \frac{\partial d_2}{\partial \tau}\\ &= q S_t e^{-q \tau} N(d_1) -r K e^{-r \tau} N(d_2) - S_t e^{-q \tau} \phi(d_1) \frac{\partial d_1}{\partial \tau} - K e^{-r \tau} \phi(d_2) \frac{\partial d_2}{\partial \tau}\\ \end{align}$

ここで、(3)式より

$\begin{align} \frac{\partial d_1}{\partial \tau} - \frac{\partial d_2}{\partial \tau} &= \frac{1}{2} \frac{\sigma}{\sqrt{\tau}} \\ \frac{\partial d_2}{\partial \tau} &= \frac{\partial d_1}{\partial \tau} - \frac{1}{2} \frac{\sigma}{\sqrt{\tau}} \end{align}$

と(2)式より

$\begin{align} \theta_c &= q S_t e^{-q \tau} N(d_1) -r K e^{-r \tau} N(d_2) - S_t e^{-q \tau} \phi(d_1) \frac{\partial d_1}{\partial \tau} - S_t e^{-q \tau} \phi(d_1) ( \frac{\partial d_1}{\partial \tau} - \frac {\sigma}{2 \sqrt{\tau}}) \\ &= q S_t e^{-q \tau} N(d_1) -r K e^{-r \tau} N(d_2) - S_t e^{-q \tau} \phi(d_1) \left( \frac{\partial d_1}{\partial \tau} - ( \frac{\partial d_1}{\partial \tau} - \frac {\sigma}{2 \sqrt{\tau}} ) \right ) \\ &= q S_t e^{-q \tau} N(d_1) -r K e^{-r \tau} N(d_2) - \frac {\sigma S_t}{2 \sqrt{\tau} } e^{-q \tau} \phi(d_1) \\ \end{align}$

プットのシータ

(4)式より

$\begin{align} \frac{\partial V_c}{\partial t} - \frac{\partial V_p}{\partial t} &= S_t \frac{\partial e^{-q \tau}}{\partial \tau} \frac{\partial \tau}{\partial t} - K \frac{\partial e^{-r \tau}}{\partial \tau} \frac{\partial \tau}{\partial t}\\ &= q S_t e^{-q \tau} - r K e^{-r \tau} \end{align}$

よって

$\begin{align} \theta_p &= \theta_c - q S_t e^{-q \tau} + r K e^{-r \tau} \end{align}$

従って

$\begin{align} \theta_p &= \theta_c - q S_t e^{-q \tau} + r K e^{-r \tau}\\ &= q S_t e^{-q \tau} N(d_1) -r K e^{-r \tau} N(d_2) - \frac {\sigma S_t}{2 \sqrt{\tau} } e^{-q \tau} \phi(d_1) - q S_t e^{-q \tau} + r K e^{-r \tau} \\ &= -q S_t e^{-q \tau} N(-d_1) + r K e^{-r \tau} N(-d_2) - \frac {\sigma S_t}{2 \sqrt{\tau} } e^{-q \tau} \phi(d_1) \end{align}$

実装

pythonで各グリークスを実装する。

準備

共通する準備の部分

#グリークスを計算する。

#値の設定
K = 105;
r = 0.01;
q = 0.02;
sigma = 0.3;
tau = 0.1; #T - t

##d1,d2を計算する。
import math
from scipy.stats import norm
import matplotlib.pyplot as plt

def d1(S,K,r,q,sigma,tau):
    d1 = (math.log(S/K) + (r - q + sigma**2 / 2) * tau) / (sigma * math.sqrt(tau))
    
    return d1

def d2(S,K,r,q,sigma,tau):
    d1 = (math.log(S/K) - (r - q + sigma**2 / 2) * tau) / (sigma * math.sqrt(tau))
    
    return d1

#株価のリストを作成
S_list = list(range(80,120))
d1_list = []
d2_list = []

#d1を計算する。
for i in S_list:
    d1_list.append(d1(i,K,r,q,sigma,tau))

#d2を計算する。
for i in S_list:
    d2_list.append(d2(i,K,r,q,sigma,tau))

デルタ

#デルタを求める。
###コールのデルタ
def call_delta(q,T,d1):
    delta = math.exp(- q * T) * norm.cdf(d1, loc = 0, scale = 1)
    
    return delta
    
###プットのデルタ
def put_delta(q,T,d1):
    delta = - math.exp(- q * T) * norm.cdf(- d1, loc = 0, scale = 1)
    
    return delta


##コール、プットのデルタを計算し、格納する。
c_delta_list = []
p_delta_list = []

for i in range(len(S_list)):
    c_delta_list.append(call_delta(q,T,d1_list[i]))
    p_delta_list.append(put_delta(q,T,d1_list[i]))))

#デルタをプロットする。
import matplotlib.pyplot as plt
fig, ax = plt.subplots()

ax.plot(S_list, c_delta_list, label="CallDelta")
ax.plot(S_list, p_delta_list, label="PutDelta")
ax.set_title("Call and Put Delta")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Delta")

plt.show();

f:id:ara-statistic:20201228022355p:plain

ガンマ

#ガンマを求める
##ガンマを求める関数を定義
def gamma(S,sigma,tau,d):
    gamma = (1 / (S * sigma * math.sqrt(tau))) * (1 / math.sqrt(2 * math.pi)) * math.exp(- d ** 2 / 2) * math.exp(- q * tau)
    
    return gamma

##ガンマを実際に求め、値を格納する。
gamma_list = []
for i in range(len(S_list)):
    gamma_list.append(gamma(S_list[i],sigma,tau,d1_list[i]))

##ガンマをプロットする。
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
plt.plot(S_list, gamma_list, label="Gamma")
plt.title("Gamma")
plt.xlabel("underlying asset price")
plt.ylabel("Gamma")

plt.show()

f:id:ara-statistic:20201228022339p:plain

ベガ

#ベガを計算する。
##ベガを計算する関数を定義する
def vega(S,sigma,tau,d1):
    vega = S * math.sqrt(tau / (2 * math.pi)) * math.exp( -d1 ** 2 / 2 - q * tau)
    
    return vega
##ベガを実際に求める。
vega_list = []

for i in range(len(S_list)):
    vega_list.append(vega(S_list[i],sigma,tau,d1_list[i]))

##ベガをプロットする。
fig, ax = plt.subplots()
plt.plot(S_list, vega_list, label="Vega")
plt.title("Vega")
plt.xlabel("underlying asset price")
plt.ylabel("Vega")

plt.show()

f:id:ara-statistic:20201228022351p:plain

ロー

#ローを計算する。
##ローを計算する関数を定義する

###コールのロー
def call_rho(K,r,tau,d2):
    rho = K * tau * math.exp(-r * tau) * norm.cdf(d2, loc = 0, scale = 1)
    
    return rho

###プットのロー
def put_rho(K,r,tau,d2):
    rho = -K * tau * math.exp(-r * tau) * norm.cdf(-d2, loc = 0, scale = 1)
    
    return rho


##ローを実際に求める。
c_rho_list = []
p_rho_list = []

for i in range(len(S_list)):
    c_rho_list.append(call_rho(K,r,tau,d2_list[i]))
    p_rho_list.append(put_rho(K,r,tau,d2_list[i]))

#ローをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_rho_list, label="CallRho")
ax.plot(S_list, p_rho_list, label="PutRho")
ax.set_title("Call and Put Rho")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Rho")

plt.show();

f:id:ara-statistic:20201228022342p:plain

イプシロン

#ローエフを計算する。
##ローエフを計算する関数を定義する

###コールのローエフ
def call_rhof(S,q,tau,d1):
    rhof = - S * tau * math.exp(- q * tau) * norm.cdf(d1, loc = 0, scale = 1)
    
    return rhof

###プットのローエフ
def put_rhof(S,q,tau,d1):
    rhof = S * tau * math.exp(- q * tau) * norm.cdf(- d1, loc = 0, scale = 1)
    
    return rhof


##ローエフを実際に求める。
c_rhof_list = []
p_rhof_list = []

for i in range(len(S_list)):
    c_rhof_list.append(call_rhof(S_list[i],q,tau,d1_list[i]))
    p_rhof_list.append(put_rhof(S_list[i],q,tau,d1_list[i]))

#ローエフをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_rhof_list, label="CallRhof")
ax.plot(S_list, p_rhof_list, label="PutRhof")
ax.set_title("Call and Put Rhof")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Rhof")

plt.show();

f:id:ara-statistic:20201228022345p:plain

シータ

#シータを計算する。

##コールのシータ
def call_theta(S,K,r,q,sigma,tau,d1,d2):
    theta = q * S * math.exp(- q * tau) * norm.cdf(d1, loc = 0, scale = 1) - r * K * math.exp(- r * tau) * norm.cdf(d2, loc = 0, scale = 1) - S * sigma * math.exp(-q * tau) / (2 * math.sqrt(tau)) * (1 / math.sqrt(2 * math.pi)) * math.exp( - d1 ** 2 / 2)
    
    return theta

###プットのシータ
def put_theta(S,K,r,q,sigma,tau,d1,d2):
    theta = - q * S * math.exp(- q * tau) * norm.cdf(- d1, loc = 0, scale = 1) + r * K * math.exp(- r * tau) * norm.cdf(- d2, loc = 0, scale = 1) - S * sigma * math.exp(-q * tau) / (2 * math.sqrt(tau)) * (1 / math.sqrt(2 * math.pi)) * math.exp( - d1 ** 2 / 2)
    
    return theta


##シータを実際に求める。
c_theta_list = []
p_theta_list = []

for i in range(len(S_list)):
    c_theta_list.append(call_theta(S_list[i],K,r,q,sigma,tau,d1_list[i],d2_list[i]))
    p_theta_list.append(put_theta(S_list[i],K,r,q,sigma,tau,d1_list[i],d2_list[i]))    


#シータをプロットする。
fig, ax = plt.subplots()

ax.plot(S_list, c_theta_list, label="CallTheta")
ax.plot(S_list, p_theta_list, label="PutTheta")
ax.set_title("Call and Put Theta")
ax.legend(loc = 0)

ax.set_xlabel("underlying asset price")
ax.set_ylabel("Theta")

plt.show();