yuima は確率過程のモデリングとその統計推測を可能にするフレームワークです．広範なクラスの確率微分方程式のシミュレーションが可能です．今回はそのような確率過程の汎函数の漸近展開に基づく計算方法を紹介します．確率変数の期待値を近似するのに Monte Carlo 法は普遍的な方法ですが，漸近展開が用いられる場合，その計算時間は比較にならないほど速くなります．
A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
例えば，\(d\)-次元拡散過程 \(X=(X^{(\epsilon)}_t)_{t\in[0,T]},\epsilon\in(0,1]\) を次のように定める： \[ X^{(\epsilon)}_t=x_0+\int^t_0a(X_s^{(\epsilon)},\epsilon)ds+\int^t_0b(X_s^{(\epsilon)},\epsilon)dW_s. \]
ただし，\(W_t\) は \(r\)-次元 Wiener 過程とした．そして，その汎函数 \[ F^{(\epsilon)}=\sum_{\alpha=0}^r\int^T_0f_\alpha(X_t^{(\epsilon)},\epsilon)dW_t^\alpha+F(X_T^{(\epsilon)},\epsilon),\quad W^0_t=t \] を推定したいとする．
このように式で表せても，特に Asian option についてはこの線型な設定においてさえ数値計算法が要請される．
ここでは，\(\epsilon\searrow0\) の極限で系が決定論的であるとする．すなわち， \[ b(-,0)=0 \] \[ f_\alpha(-,0)=0\quad(\alpha\in[r]) \] とする．すると \(X_t^{(0)}\) は次の常微分方程式 \[ \frac{d X_t^{(0)}}{d t}=a(X_t^{(0)},0),\quad X_0^{(0)}=x_0 \] の解であるから，\(F^{(0)}\) も定数 \[ F^{(0)}=\int^T_0f_0(X_t^{(0)},0)dt+F(X_T^{(0)},0) \] で与えられる．
さらに，\(a,b,f_\alpha,F\) がしかるべき正則性条件を満たすとき，汎函数 \(F^{(\epsilon)}\) にはある版が存在して \(\epsilon\in[0,1)\) に関して殆ど確実に滑らかである．特に， \[ \widetilde{F}^{(\epsilon)}:=\frac{F^{(\epsilon)}-F^{(0)}}{\epsilon} \] は次の確率展開を持つ： \[ \widetilde{F}^{(\epsilon)}\sim\widetilde{F}^{[0]}+\epsilon\widetilde{F}^{[1]}+\epsilon^2\widetilde{F}^{[2]}+\cdots\quad(\epsilon\searrow0) \]
この展開は，Malliavin 解析の Sobolev 空間において厳密に成り立つ．これを導くのが (Watanabe, 1987) の理論である．
これに基づき，汎函数 \(F^{(\epsilon)}\) の近似を構成する機能が yuima に実装されている．
この漸近展開をオプションの価格付けに応用したのが (Yoshida, 1992) である．
setFunctionalコンストラクタsetFunctional(model, F, f, xinit, e)Black-Scholes モデル \[ dX_t^{(\epsilon)}=\mu X_t^{(\epsilon)}dt+\epsilon X_t^{(\epsilon)}dW_t \] が定める幾何 Brown 運動 \((X_t)\) のパラメータが \(\mu=1,x_0=1\) を満たす場合において，Asian call option の価格は，汎函数 \[ g(x):=\max\left(F^{(0)}-K+\epsilon x,0\right) \] を計算すれば良い．
\(F^{(\epsilon)}\) の極限 \(F^{(0)}\) の値は，関数 F0 を yuima.functional スロットが埋まった yuima オブジェクトに適用することで得られる．
library(yuima)
model <- setModel(drift="x", diffusion=matrix("x*e", 1, 1))
K <- 100
yuima <- setYuima(model=model, sampling=setSampling(Terminal=1, n=1000))
yuima <- setFunctional(yuima, f=list(expression(x/T), expression(0)), F=0, xinit=150, e=0.5)
F0 <- F0(yuima)print(F0)[1] 257.6134g <- function(x) {
tmp <- (F0 - 100) + (0.5*x)
tmp[(0.5*x) < (100-F0)] <- 0
tmp
}
asymp <- asymptotic_term(yuima, block=10, expression(0), g)str(asymp)List of 3
$ d0: num 159
$ d1: num -3.93
$ d2: num [1, 1] 4これを適切な和をとれば良い．
e = 0.5
asy1 <- asymp$d0 + e*asymp$d1
asy1[1] 156.608asy2 <- asymp$d0 + e*asymp$d1 + e^2 * asymp$d2
asy2 [,1]
[1,] 157.6082は Asian call price の，それぞれ１次と２次の漸近展開を与える．
この設定では対数正規分布に対する Edgeworth 展開によっても計算ができる (Levy, 1992)．
\(X_t\) が幾何 Brown 運動の場合にしか (Levy, 1992) の近似は用いることはできないが，yuima のアプローチでは可能である
例えば利子率の期間構造のモデルである CIR 模型 (Cox et al., 1985) \[ dX_t=(\alpha-\beta X_t)\,dt+\sqrt{\gamma X_t}\,dW_t \] がこの境界例になっている．
\(X\) が \[ dX_t=0.9X_tdt+\epsilon\sqrt{X_t}dW_t,\quad X_0=1, \] である場合の \(K=10,\epsilon=0.4\) における European call option の価格を考える．
a <- 0.9; e <- 0.4; Terminal <- 3; xinit <- 1; K <- 10
drift <- "a * x"; diffusion <- "e * sqrt(x)"
model <- setModel(drift = drift, diffusion = diffusion)
n <- 1000 * Terminal
yuima <- setYuima(model = model, sampling = setSampling(Terminal = Terminal, n = n))
f <- list(c(expression(0)), c(expression(0)))
F <- expression(x)
yuima.ae <- setFunctional(yuima, f = f, F = F, xinit = xinit, e = e)
rho <- expression(0)
F1 <- F0(yuima.ae)
get_ge <- function(x, epsilon, K, F0) {
tmp <- (F0 - K) + (epsilon * x[1])
tmp[(epsilon * x[1]) > (K - F0)] <- 0
return(-tmp)
}
g <- function(x) {
return(get_ge(x, e, K, F1))
}
time1 <- proc.time()
asymp <- asymptotic_term(yuima.ae, block = 100, rho, g)
time2 <- proc.time()ae.value0 <- asymp$d0
ae.value0[1] 0.7219652ae.value1 <- asymp$d0 + e * asymp$d1
ae.value1[1] 0.5787545ae.value2 <- as.numeric(asymp$d0 + e * asymp$d1 + e^2 * asymp$d2)
ae.value2[1] 0.5617722ae.time <- time2 - time1
ae.time ユーザ システム 経過
0.537 0.016 0.561 この状態での European call option の価格は，２次までの漸近展開が与える値が 100 万データ数による Monte Carlo 推定量の精度に匹敵し，当然計算量は漸近展開の方が圧倒的に少ない．