モンテカルロ法で条件付期待値計算をする際の試行錯誤−１

問題設定

ここではある時点 $t$ における、適当な確率変数の条件付期待値を計算する事を考えたい。例えば適当な確率空間 $(\Omega, \mathcal{F}_t, \mathbb{P})$ を設定し、以下のような条件付き期待値の計算を行うとする*1。
$V(S_t, t) = E\left[\left. e^{-r(T-t)} max(S_T - K, 0) \right| \mathcal{F_t} \right]$ ・・・(A)
ここで $S_t$ は確率測度 $\mathbb{P}$ の下で以下の確率微分方程式に従うものとしている。
$\frac{dS_t}{S_t} = r {\rm d}t + \sigma {\rm d}W_t$
ただし、 $r$ , $\sigma$ はそれぞれ一定値とし、 ${\rm d}W_t$ は標準ブラウン運動（の差分）とする。この確率微分方程式をモンテカルロシミュレーションさせる際には、数値差分法に起因する余分な計算誤差を減らすために一度確率微分方程式解くことで
$S_t = S_0 \exp^{(r-0.5\sigma^2)t+\sigma W_t}$
という形にしてから実行することにした。

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

この問題設定は俗にいうブラックショールズ方程式の枠組みが当てはまるものとなっているので、その解は
$V(S_t, t) = S_t \Phi(d1) - K \exp^{-r(T-t)}\Phi(d2)$
$d1 = \frac{\log(\frac{S_t}{K})+(r+0.5\sigma^2)(T-t)}{\sigma \sqrt{T-t}}$
$d2 = d1 - \sigma \sqrt{T-t}$
$\Phi(x)$ :標準正規分布の累積分布関数
と表現する事ができる。これを基準にモンテカルロ法での評価値の良しあしを考える。

素直なモンテカルロ法での実装

(A)を素直にモンテカルロ法で評価してやる場合には、モンテカルロパス数を $N$ として
$V(S_t, t) \sim \frac{1}{N} \sum_{i=1}^{N} e^{-r(T-t)} max(S_T^{(i)} - K, 0)$ ・・・(B)
というある種の標本平均値をモンテカルロシミュレーションを通して計算した ${S_T^{(i)}, i=1 \cdots N}$ から計算してやればよい。
これはその辺の本に載ってる本当に素直なモンテカルロのやり方。

塔定理を使った変形と実装方法

一方、条件付期待値に関する俗に言う塔定理（tower property)を使うと
$V(S_t, t) = E\left[\left. e^{-r(T-t)} max(S_T - K, 0) \right| \mathcal{F_t} \right] = E\left[\left. E\left[\left. e^{-r(T-t)} max(S_T - K, 0) \right| \mathcal{F_s} \right] \right| \mathcal{F_t} \right]$
$= E\left[\left. e^{-r(s-t)} E\left[\left. e^{-r(T-s)} max(S_T - K, 0) \right| \mathcal{F_s} \right] \right| \mathcal{F_t} \right] = E\left[\left. e^{-r(s-t)} V(S_s, s) \right| \mathcal{F_t} \right]$ （ただし $s > t$ ）
と変形することができるので、これを使って解析解とモンテカルロシミュレーションの組み合わせで条件付期待値を評価してやりたい、つまり
$V(S_t, t) \sim \frac{1}{N} \sum_{i=1}^{N} e^{-r(s-t)} V(S_s^{(i)}, s)$ ・・・(C)
という量を計算することで評価値を算出する。ここで(C)のやり方の場合、ブラックショールズ方程式の解をモンテカルロパス数分計算しなければならない点に注意。

teramonagi氏が本当にやりたかったこと

ここで私がとりあえずやりたい事は
「(B)と(C)の計算方法においてどちらがより大きくor小さくモンテカルロ誤差（標準誤差）が出てくるのか？」
ということなんで、とりあえずやってみた。特に(C)の場合、塔定理を適用する際の時点 $s$ に関して任意性( $0 <= s <= T$ )があるので、そこはパラメーターとして複数のsに対して計算を実施。結果は以下のようになった。

（横軸が $s$ 、縦軸が標準誤差、緑線が(B)、青線が(C)に対応。パス数は双方10,000パスで固定）
・・・つまり条件付のタイミングが早い( $s\sim0$ )の場合はブラックショールズ方程式の解析解をモンテカルロ数分だけ計算することになるので、ほぼモンテカルロ誤差がなくなる一方、条件付のタイミングが遅い（満期に近い、 $s \sim T$ ）の場合は満期時点でのペイオフ関数( $\max(x-K,0)$ )をモンテカルロパスごとに計算することになるので、結果としてはほぼ素直なモンテカルロ法の結果と変わらないということか。非常に素直な結論。逆にそんなのあるのかしらないけど、条件付期待値をある時点でだけ解析的に評価できるのならばモンテカルロ誤差は減らせるということでもある。

モンテカルロ誤差自体は
$(B) \sim \frac{\exp^{-r(T-t)}stdev(max(S_T-K,0))}{\sqrt{N}}$
$(C) \sim \frac{\exp^{-r(T-s)}stdev(V(S_s,s))}{\sqrt{N}}$
のように見積もれるので、 $max(S_T-K,0)$ や $V(S_s,s)$ の分布から考察するってのもありか。

以下、計算に使ったコード。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
#generate underlying(stock) path by geometric brownian motion
def generate_path(s0, T, r, vol, N, M):
    dt = T/M
    w = np.cumsum(np.reshape(np.random.standard_normal(N*M), (N,M)), 1) * (dt**0.5)
    t = np.cumsum(np.ones((N,M)), 1)*dt   
    return s0 * np.exp((r-0.5*vol**2)*t + vol*w)

#call option price by black-scholes model
def call_price_bs(s0, k, T, r, vol):
    d1 = (np.log(s0/k)+(r+0.5*(vol**2.0))*T) / (vol * np.sqrt(T))
    d2 = d1-vol*np.sqrt(T)
    return (s0*norm.cdf(d1) - k*np.exp(-r*T)*norm.cdf(d2))

#call option price by montecarlo    
def call_price_monte(s0, k, T, r, vol, N, M):
    s = generate_path(s0, T, r, vol, N, M)
    discount = np.exp(-r*T)
    callprice = map(lambda x: discount*max(x, 0), s[:, -1] - k)
    return (np.average(callprice), np.std(callprice)/(N**0.5))

#call option price by conditional montecarlo
def call_price_monte_conditional(s0, k, t, T, r, vol, N, M):
    s_t = generate_path(s0, t, r, vol, N, M)[:, -1]
    discount = np.exp(-r*t)
    callprice = [discount * call_price_bs(s, k, T-t, r, vol) for s in s_t]
    return (np.average(callprice), np.std(callprice)/(N**0.5))
    

if __name__ == "__main__":
    #initial stock price
    S0  = 100.0
    #strike    
    K   = 100.0
    #maturity
    T   = 3.0
    #interest rate    
    R   = 0.01
    #volatility    
    VOL = 0.2
    #number of montecarlo path
    N = 10000
    #number of time-grid    
    M = 200

    print "Exact solution"
    print call_price_bs(S0, K, T, R, VOL)
    print "Simple montecarlo "
    price_monte = call_price_monte(S0, K, T, R, VOL, N, M)
    print "montecarlo with tower property"
    ts = np.linspace(0, T, 10, endpoint=True)
    price_monte_conditional = [call_price_monte_conditional(S0, K, t, T, R, VOL, N, M) for t in ts]
    #plot montecarlo error   
    plt.plot(ts, map(lambda x:x[1], price_monte_conditional), ts, [price_monte[1] for i in ts])
    plt.show()

*1:金融工学でいうユーロピアンコールオプションの評価式。めんどくさいから金利は一定値に固定