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

CIRモデルについてもVasicekモデルと同様にシンプルに最尤推定でモデルのパラメータを推定できる関数を用意しておく

# maximum likelihood estimation for CIR model
#
# rate      : interest rate
# param     : model parameter
# delta     : time step
# initparam : initial parameter
#
CIRloglik <- function(rate, param, delta) {
  n <- length(rate)
  c <- 4 * param[1] / (param[3] ^ 2 * (1 - exp(-param[1] * delta)))
  df <- 4 * param[1] * param[2] / (param[3] ^ 2)
  ncp <- c * rate[1:(n - 1)] * exp(-param[1] * delta)
  sum(dchisq(c * rate[2:n], df = df, ncp = ncp, log = TRUE) + log(c))
}

CIRMLE <- function(rate, initparam, delta) {
  opt <- optim(initparam, CIRloglik, 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)
}

パラメータを設定して、パスを発生させる

library(sde)

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

推定結果を確認する

> CIRMLE(rate, initparam, delta)
       kappa           mu        sigma       loglik  convergence 
2.789964e-01 6.105776e-02 8.106622e-02 2.376850e+03 0.000000e+00

初期値の与え方にもよるが、Vasicekモデルと比べて多少パラメータの推定値にバイアスが出ることになります。 特に平均回帰パラメータkappaについては上方に推定する傾向にあります