今週も特にありません

進捗どうですか?

Vasicekモデルでパスの発生方法を比較

Vasicekモデルの債券価格債券オプションについて書いたので、次はパスの発生方法を比較してみる

はじめはEuler–Maruyama method - Wikipedia, the free encyclopediaに書かれている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")

f:id:masaqol:20150217001831p:plain

時間経過とともに長期の平均回帰水準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パッケージを使ったものはほとんど同じという結果

次は、モンテカルロシミュレーションでそれぞれのパスから債券価格を求めてみます