Vasicekモデルのパラメータ推定

SMFI5パッケージでもVasicekモデルのパラメータ推定はできるが、参考程度にパラメータを推定する際、シンプルに最尤推定できるものを用意しておいて、簡単に使い回せるようにしておきたい

# maximum likelihood estimation for Vasicek model
#
# rate      : interest rate
# param     : Vasicek model parameter
# delta     : time step
# initparam : initial parameter
#
Vasicekloglik <- function(rate, param, delta) {
  n <- length(rate)
  m <- param[2] + (rate[1:(n-1)] - param[2]) * exp(-param[1] * delta)
  v <- param[3]^2 * (1 - exp(-2 * param[1] * delta)) / (2 * param[1])
  sum(dnorm(rate[2:n], mean = m, sd = sqrt(v), log = TRUE))
}
  
VasicekMLE <- function(rate, initparam, delta) {
  opt <- optim(initparam, Vasicekloglik, rate = rate, delta = delta, control = list(fnscale = -1))
  c(kappa = opt$par[1], mu = opt$par[2], sigma = opt$par[3], loglik = opt$value, convergence = opt$convergence)
}

optimのhessian = TRUEを足したり、２つの関数を１つにしてしまったりということもできる

sdeパッケージを使ってパスを発生させ、パラメータ推定する

library(sde)

r0 <- 0.01
n <- 500
delta <- 1/52
kappa <- 0.2
mu <- 0.03
sigma <- 0.01
initparam <- c(kappa, mu, sigma)
set.seed(123)
rate <- sde.sim(X0 = r0, N = n, delta = delta, theta = c(kappa * mu, kappa, sigma), model = "OU")

結果を確認する

> VasicekMLE(rate, initparam, delta)
       kappa           mu        sigma       loglik  convergence 
5.670395e-01 3.114274e-02 9.740629e-03 2.596687e+03 0.000000e+00 

最適化の初期値に真の値を入れていることもあり、良い値に推定できています

金利の平均回帰速度kappaはだいぶ真の値よりも上方に推定されていますが、推定バイアスを改善させましょうというような論文がいくつかありますので、この辺をじっくり読むというのも大変面白いです