［R］トレンドのあるモデルの比較（水準成分と傾き成分の推移）（「カルマンフィルタ」（共立出版）pp.76-79）

公式ページからダウンロードした入力ファイル「Weight.dat」をカレントディレクトリに置いておくこと。

> library(KFAS)
> dtf <- read.table("Weight.dat", header = FALSE)
> yt <- dtf$V1
> # ローカルレベルモデル 灰実線
> mod1tr <- SSModel(yt ~ SSMtrend(1, Q = NA), H = NA)
> fit1tr <- fitSSM(mod1tr, numeric(2), method = "BFGS")
> kfs1tr <- KFS(fit1tr$model)
> aht1tr <- kfs1tr$alphahat[, "level"]
> # 2次のトレンドモデル 黒実線
> mod2tr <- SSModel(yt ~ SSMtrend(2, Q = c(list(0), list(NA))), H = NA)
> fit2tr <- fitSSM(mod2tr, numeric(2), method = "BFGS")
> kfs2tr <- KFS(fit2tr$model)
> aht2tr <- kfs2tr$alphahat[, "level"]
> ahts2tr <- kfs2tr$alphahat[, "slope"]
> # ローカル線形トレンドモデル 黒破線
> modllt <- SSModel(yt ~ SSMtrend(2, Q = c(list(NA), list(NA))), H = NA)
> fitllt <- fitSSM(modllt, numeric(3), method = "BFGS")
> kfsllt <- KFS(fitllt$model)
> ahtllt <- kfsllt$alphahat[, "level"]
> ahtsllt <- kfsllt$alphahat[, "slope"]
> #
> png("kalmanfil_p078_1.png", width = 768, height = 512)
> plot(yt, lty = 3, type = "o", xlab = "経過日数", ylab = "水準成分")
> lines(aht1tr, lwd = 2, col = 8)  # ローカルレベルモデル(灰実線)
> lines(aht2tr, lwd = 2)           # 2次のトレンドモデル(黒実線)
> lines(ahtllt, lwd = 2, lty = 2)  # ローカル線形トレンドモデル(黒破線)
> dev.off()
null device 
          1 
> #
> png("kalmanfil_p078_2.png", width = 768, height = 512)
> plot(ahts2tr, lwd = 2, xlab = "経過日数", ylab = "傾き成分")  # 2次のトレンドモデル(黒実線)
> lines(ahtsllt, lwd = 2, lty = 2)  # ローカル線形トレンドモデル(黒破線)
> dev.off()
null device 
          1 
> #
> r <- qv <- seta12 <- seta22 <- seps2 <- mll <- aic <- mse <- double(3)
> r <- c(2, 2, 3)
> qv <- c(1, 2, 2)
> seta12[1] <- as.vector(fit1tr$model$Q[, , 1])
> seta12[2] <- as.vector(diag(fit2tr$model$Q[, , 1])[1])
> seta22[2] <- as.vector(diag(fit2tr$model$Q[, , 1])[2])
> seta12[3] <- as.vector(diag(fitllt$model$Q[, , 1])[1])
> seta22[3] <- as.vector(diag(fitllt$model$Q[, , 1])[2])
> seps2[1] <- as.vector(fit1tr$model$H[, , 1])
> seps2[2] <- as.vector(fit2tr$model$H[, , 1])
> seps2[3] <- as.vector(fitllt$model$H[, , 1])
> mll[1] <- kfs1tr$logLik - sum(kfs1tr$Finf > 0) * log(2 * pi) / 2
> mll[2] <- kfs2tr$logLik - sum(kfs2tr$Finf > 0) * log(2 * pi) / 2
> mll[3] <- kfsllt$logLik - sum(kfsllt$Finf > 0) * log(2 * pi) / 2
> aic <- -2 * mll + 2 * (r + qv)
> idx <- 3:60
> n <- length(idx)
> mse[1] <- sum(kfs1tr$v[idx] ^ 2) / n
> mse[2] <- sum(kfs2tr$v[idx] ^ 2) / n
> mse[3] <- sum(kfsllt$v[idx] ^ 2) / n
> dtf <- data.frame(r, q = qv, seta12, seta22, seps2, mll, aic, mse)
> # 表3.2
> print(dtf)
  r q     seta12       seta22     seps2       mll      aic       mse
1 2 1 0.07078040 0.000000e+00 0.1491052 -48.58750 103.1750 0.2984077
2 2 2 0.00000000 7.097290e-04 0.2432792 -54.89776 117.7955 0.3581985
3 3 2 0.08655912 5.260979e-07 0.1371066 -51.74947 113.4989 0.3234503