pforeachで最尤推定した漸近分布を確認する

最尤推定で推定したパラメータの漸近分布をシミュレーションによって確認する
繰り返し推定を行う際にpforeach(hoxo-m/pforeach · GitHub)を使う

CIRモデルのパラメータ推定と関数もパラメータの設定も同じとする

r0 <- 0.01
n <- 1000
delta <- 1/52
kappa <- 0.2
mu <- 0.03
sigma <- 0.08
initparam <- c(kappa, mu, sigma)

パスを発生させ、パラメータ推定を繰り返す。収束しなかったものは取り除く

ggplot2で推定したパラメータの分布をプロットしたいので、最終的にdata.frameで返ってくるようにする

.c = data.frameとして、data.frameのまま結合してもパラメータの名前がついたdata.frameが返って来なかったため、rbindで結果を結合したあとにdata.frameに変換している（これが最適な方法なのかについては…）

library(sde)
library(dplyr)
library(pforeach)

est <- pforeach(i = 1:1000, .c = rbind)({
  sde.sim(X0 = r0, N = n, delta = delta, 
          theta = c(kappa * mu, kappa, sigma), model = "CIR") %>%
  CIRMLE(initparam, delta)
}) %>% as.data.frame %>% filter(convergence != 1)

foreachの少し込み入った使い方に面倒くささを感じていたが、pforeachを使えば簡単、簡潔に記述できる

replicateを使う場合には、データを発生させて、パラメータを推定してという一連の処理を書いた関数を書いておいて、その関数を繰り返すということをしなければならないが、pforeachではループの中にザッと処理を書いておくだけで繰り返し推定が行える。返り値もあるので、その後一度に加工して扱える。素晴らしい

推定したパラメータの分布をプロットする

library(ggplot2)

ggplot(est, aes(x = kappa)) + geom_histogram(binwidth = 0.1)
ggplot(est, aes(x = mu)) + geom_histogram(binwidth = 0.005)
ggplot(est, aes(x = sigma)) + geom_histogram(binwidth = 0.0001)

平均回帰パラメータkappaの分布のみここでは確認する。推定された値の平均を計算すると真値よりも大きい値となる。そして、プロットしても上方に推定するバイアスが生じることが確認できる（TRUE = 0.2）