Vasicekモデルでパスの発生方法を比較
Vasicekモデルの債券価格と債券オプションについて書いたので、次はパスの発生方法を比較してみる
はじめはEuler–Maruyama method - Wikipediaに書かれているPythonコードをそのままRに書き直したもの。愚直にforループを回してみる
Vasicek_euler <- function(kappa, mu, sigma, r0, dt, n) { r <- numeric(n) r[1] <- r0 w <- rnorm(n) for (i in 2:n) { r[i] <- r[i-1] + kappa * (mu - r[i-1]) * dt + sigma * sqrt(dt) * w[i-1] } r }
週次データで10年までを想定してパラメータをセット
kappa <- 0.2 mu <- 0.015 sigma <- 0.003 r0 <- 0.001 dt <- 1/52 n <- 520
パスを発生させる
> set.seed(123) > p1 <- Vasicek_euler(kappa, mu, sigma, r0, dt, n) > head(p1) [1] 0.0010000000 0.0008206742 0.0007794504 0.0014826067 0.0015639300 0.0016693941
複数パスを発生させて、プロット
library(dplyr) library(tidyr) library(ggplot2) time <- 1:n * dt paths <- replicate(30, Vasicek_euler(kappa, mu, sigma, r0, dt, n)) cbind(time, paths) %>% as.data.frame() %>% gather(path, rate, -time) %>% ggplot(aes(x = time, y = rate, group = path)) + geom_line(aes(colour = path)) + ggtitle("Short Rate Paths of Vasicek Model") + theme(legend.position = "none")
時間経過とともに長期の平均回帰水準0.015付近で推移していくようになる。今ではものすごくおかしなことではないのかもしれないが、Vasicekモデルでは金利がマイナスの状態になることもある
VasicekモデルはOU過程であり、その平均と分散も明確にわかるので、より厳密に書くことができる
Vasicek_exact <- function(kappa, mu, sigma, r0, dt, n) { r <- numeric(n) r[1] <- r0 w <- rnorm(n) e <- exp(-kappa * dt) s <- sigma * sqrt((1 - e^2) / (2 * kappa)) for (i in 2:n) { r[i] <- e * r[i-1] + mu * (1 - e) + s * w[i-1] } r }
こちらも同じパラメータでパスを発生させて、前者と比較する
> set.seed(123) > p2 <- Vasicek_exact(kappa, mu, sigma, r0, dt, n) > head(p2) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p1 - p2) [1] 0.000000e+00 -3.442722e-07 -4.220666e-07 9.296463e-07 1.082250e-06 1.280618e-06
多少ではあるが、誤差が出ることがわかる
パッケージを使ったものとも比較してみる。sdeパッケージのsde.simは異なる方法でパスを発生させることできる。初期値からnステップ発生するので、n+1番目は抜いておく
library(sde) set.seed(123) p3 <- sde.sim(X0 = r0, N = n, delta = dt, theta = c(kappa * mu, kappa, sigma), model = "OU")[-(n+1)] set.seed(123) p4 <- sde.sim(X0 = r0, N = n, delta = dt, method = "cdist", rcdist = rcOU, theta = c(kappa * mu, kappa, sigma))[-(n+1)] set.seed(123) p5 <- sde.sim(X0 = r0, N = n, delta = dt, drift = expression(0.2 * (0.015 - x)), sigma = expression(0.003))[-(n+1)]
それぞれ、モデルを直接指定する方法、条件付き分布から発生させる方法、ドリフト項とディフュージョン項を指定する方法
こちらも一部を確認
> head(p3) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p4) [1] 0.0010000000 0.0008210185 0.0007798725 0.0014816771 0.0015628477 0.0016681135 > head(p5) [1] 0.0010000000 0.0008210190 0.0007798732 0.0014816757 0.0015628461 0.0016681116
当然だが、自作関数のものとほとんど同じ値になることを確認できる
yuimaパッケージはさらに複雑なものもシミュレーションすることができるようになっている。ここでは、一番簡単そうな方法を一つだけやってみる
library(yuima) set.seed(123) sampling <- setSampling(n = n, delta = dt) model <- setModel(drift = "kappa * (mu - x)", diffusion = "sigma", state.var = "x", time.var = "t", solve.var = "x", xinit = r0) p6 <- simulate(model, true.param = list(kappa = kappa, mu = mu, sigma = sigma), sampling = sampling)@data@original.data[-(n+1)]
一部を確認
> head(p6) [1] 0.0010000000 0.0008206742 0.0007794504 0.0014826067 0.0015639300 0.0016693941
結局、どれくらいの違いがあるのか？
> sum((p1 - p2)^2) [1] 5.358932e-08 > sum((p1 - p3)^2) [1] 5.358932e-08 > sum((p1 - p4)^2) [1] 5.358932e-08 > sum((p1 - p5)^2) [1] 5.373786e-08 > sum((p1 - p6)^2) [1] 1.982772e-32 > sum((p2 - p3)^2) [1] 3.662136e-30 > sum((p2 - p4)^2) [1] 3.662136e-30 > sum((p2 - p5)^2) [1] 1.03119e-13 > sum((p2 - p6)^2) [1] 5.358932e-08 > sum((p3 - p4)^2) [1] 0 > sum((p3 - p5)^2) [1] 1.03119e-13 > sum((p3 - p6)^2) [1] 5.358932e-08 > sum((p4 - p5)^2) [1] 1.03119e-13 > sum((p4 - p6)^2) [1] 5.358932e-08 > sum((p5 - p6)^2) [1] 5.373786e-08
多少の違いはあるようですが、sdeパッケージを使った3番目と4番目は全く同じ。2番目のものとsdeパッケージを使ったものとはほとんど同じ。1番はじめのEuler法のものとyuimaパッケージを使ったものはほとんど同じという結果
次は、モンテカルロシミュレーションでそれぞれのパスから債券価格を求めてみます
