MMM研修 VII. 状態空間モデルの基礎
Shigeru ONO
2020年4月9日
このページでは、株式会社インサイト・ファクトリー 2020年度社内研修「マーケティング・ミックス・モデリング」VII章「状態空間モデルの基礎」のRコードを公開しています。
準備
library(tidyverse)## -- Attaching packages --------------------------------------------------------- tidyverse 1.3.0 --## √ ggplot2 3.3.0 √ purrr 0.3.3
## √ tibble 3.0.0 √ dplyr 0.8.5
## √ tidyr 1.0.2 √ stringr 1.4.0
## √ readr 1.3.1 √ forcats 0.5.0## -- Conflicts ------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(assertr)
library(patchwork)
library(dlm)##
## Attaching package: 'dlm'## The following object is masked from 'package:ggplot2':
##
## %+%library(KFAS)
WORKDIR <- "."
# dfSeatbelts
dfSeatbelts <- tibble(
dMonth = seq(as.Date("1969-01-01"), as.Date("1984-12-01"), by = "month"),
gDriver = as.vector(Seatbelts[,"drivers"]),
gPetrol = as.vector(Seatbelts[,"PetrolPrice"]),
gLogDriver = log(gDriver),
gLogPetrol = log(gPetrol)
) %>%
mutate(nID = seq_along(dMonth))Code 1. チャート
dfPlot <- dfSeatbelts %>%
gather(sVar, gValue, c(gLogDriver, gLogPetrol)) %>%
mutate(
fVar = factor(
sVar,
levels = c("gLogDriver", "gLogPetrol"),
labels = c("運転者死傷数(対数)", "原油価格(対数)")
)
)
g <- ggplot(data=dfPlot, aes(x = dMonth, y = gValue, color = fVar))
g <- g + geom_line()
g <- g + labs(x = "年月", y = NULL)
g <- g + theme_bw()
g <- g + guides(color = FALSE)
g <- g + facet_grid(fVar ~ ., scales = "free")
print(g)
ggsave(paste0(WORKDIR, "/chap7_code1.png"), width = 20, height = 12, units = "cm")Code 2. ローカルレベルモデル, 確定的レベル
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(1), # Z
T = matrix(1), # T
R = matrix(1), # R
Q = matrix(0), # Q
P1inf = diag(1) # 散漫初期状態の共分散行列. 説明略
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
## ローカル・レベル・モデルなので、SSMtrend()を使って以下のように略記できる
# oModel <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(0))),
# H = matrix(NA)
# )
#
# 撹乱項の分散の初期値のベクトル。適当な値でよい
gInit <- var(dfSeatbelts$gLogDriver) * 0.5
# パラメータ推定
oFitted <- fitSSM(
oModel,
inits = log(gInit),
method = "BFGS" # 状態変数がすべて確定的な場合は、これを指定するとよいらしい
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.02935256 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat = as.vector(oEstimated$alphahat), # スムージングによる状態推定
gAtt = as.vector(oEstimated$att) , # フィルタリングによる状態推定
gError = gLogDriver - gAlphaHat
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g3 <- g3 + geom_line(color = "gray")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g <- g1 + g2 + g3
print(g)
ggsave(paste0(WORKDIR, "/chap7_code2.1.png"), width = 24, height = 12, units = "cm")
# スムージングとフィルタリングのちがい
g <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g <- g + geom_line()
g <- g + geom_line(aes(y = gAlphaHat), color = "red")
g <- g + geom_line(aes(y = gAtt), color = "blue")
g <- g + labs(x = "年月", y = NULL)
g <- g + theme_bw()
print(g)
ggsave(paste0(WORKDIR, "/chap7_code2.2.png"), width = 20, height = 12, units = "cm")code 3. ローカルレベルモデル, 確率的レベル
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(1), # Z
T = matrix(1), # T
R = matrix(1), # R
Q = matrix(NA), # Q. NAは推定対象であることを表す
P1inf = diag(1)
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# ローカル・レベル・モデルなので、SSMtrend()を使って以下のように略記できる
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))),
H = matrix(NA)
)
# 撹乱項の分散の初期値のベクトル。先頭は観察撹乱項。適当な値でよい
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001)
# パラメータ推定
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.002220796 cat("sigma^2_xi =", as.vector(oFitted$model$Q), "\n")## sigma^2_xi = 0.01186672 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat = as.vector(oEstimated$alphahat), # スムージングによる状態推定
gAtt = as.vector(oEstimated$att), # フィルタリングによる状態推定
gError = gLogDriver - gAlphaHat
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g3 <- g3 + geom_line(color = "gray")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g <- g1 + g2 + g3
print(g)
ggsave(paste0(WORKDIR, "/chap7_code3.1.png"), width = 24, height = 12, units = "cm")
# スムージングとフィルタリング
g <- ggplot(data=dfPlot %>% slice(1:12), aes(x = dMonth, y = gLogDriver))
g <- g + geom_line()
g <- g + geom_point()
g <- g + geom_line(aes(y = gAlphaHat), color = "red")
g <- g + geom_point(aes(y = gAlphaHat), color = "red")
g <- g + geom_line(aes(y = gAtt), color = "blue")
g <- g + geom_point(aes(y = gAtt), color = "blue")
g <- g + labs(x = NULL, y = NULL)
g <- g + theme_bw()
print(g)
ggsave(paste0(WORKDIR, "/chap7_code3.2.png"), width = 20, height = 12, units = "cm")code 4. ローカル線形トレンド・モデル (確定的レベルと確定的傾き)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(1, 0), nrow=1), # Z
T = matrix(c(1,0,1,1), nrow=2), # T
R = matrix(c(1,0,0,1), nrow=2), # R
Q = matrix(c(0,0,0,0), nrow=2), # Q
P1inf = diag(c(1, 1))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# 線形トレンド・モデルなので、SSMtrend()を使って以下のように略記できる
# oModel <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=2, Q=list(matrix(0), matrix(0))),
# H = matrix(NA)
# )
# パラメータ推定
gInit <- var(dfSeatbelts$gLogDriver) / 2
oFitted <- fitSSM(
oModel,
inits = log(gInit),
method = "BFGS"
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.02299787 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gAtt = as.vector(oEstimated$att[,1]),
gAtt2 = as.vector(oEstimated$att[,2]),
gError = gLogDriver - gAlphaHat
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g3 <- g3 + geom_line(color = "gray")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g <- g1 + g2 + g3
print(g)
ggsave(paste0(WORKDIR, "/chap7_code4.1.png"), width = 24, height = 12, units = "cm")
# スムージングとフィルタリング
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + geom_line(aes(y = gAlphaHat), color = "red")
g1 <- g1 + geom_line(aes(y = gAtt), color = "blue")
g1 <- g1 + labs(x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g2 <- g2 + geom_line(color="red")
g2 <- g2 + geom_line(aes(y = gAtt2), color="blue")
g2 <- g2 + labs(x = NULL, y = "v[t]")
g2 <- g2 + theme_bw()
g <- g1 + g2 + plot_layout(nrow = 2)
print(g)
ggsave(paste0(WORKDIR, "/chap7_code4.2.png"), width = 24, height = 12, units = "cm")code 5. ローカル線形トレンド・モデル(確率的レベルと確定的傾き)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(1, 0), nrow=1), # Z
T = matrix(c(1,0,1,1), nrow=2), # T
R = matrix(c(1,0,0,1), nrow=2), # R
Q = matrix(c(NA,0,0,0), nrow=2), # Q
P1inf = diag(c(1, 1))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# View(oModel)
# # 線形トレンド・モデルなので、SSMtrend()を使って以下のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=2, Q=list(matrix(NA), matrix(0))),
# H = matrix(NA)
# )
# # View(oModel2)
# # パラメータ推定
agInit <- c(var(dfSeatbelts$gLogDriver) / 2, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.002116863 cat("sigma^2_xi =", as.vector(oFitted$model$Q[1,1,1]), "\n")## sigma^2_xi = 0.01212854 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gAtt = as.vector(oEstimated$att[,1]),
gAtt2 = as.vector(oEstimated$att[,2]),
gError = gLogDriver - gAlphaHat
)
# View(dfPlot)
# stop()
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g3 <- g3 + geom_line(color = "gray")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g <- g1 + g2 + g3
print(g)
ggsave(paste0(WORKDIR, "/chap7_code5.1.png"), width = 24, height = 12, units = "cm")
# スムージングとフィルタリング
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + geom_line(aes(y = gAlphaHat), color = "red")
g1 <- g1 + geom_line(aes(y = gAtt), color = "blue")
g1 <- g1 + labs(x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g2 <- g2 + geom_line(color="red")
g2 <- g2 + geom_line(aes(y = gAtt2), color="blue")
g2 <- g2 + labs(x = NULL, y = "v[t]")
g2 <- g2 + theme_bw()
g <- g1 + g2 + plot_layout(nrow = 2)
print(g)
ggsave(paste0(WORKDIR, "/chap7_code5.2.png"), width = 24, height = 12, units = "cm")code 6. ローカル線形トレンド・モデル(確率的レベルと確率的傾き)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(1, 0), nrow=1), # Z
T = matrix(c(1,0,1,1), nrow=2), # T
R = matrix(c(1,0,0,1), nrow=2), # R
Q = matrix(c(NA,0,0,NA), nrow=2), # Q
P1inf = diag(c(1, 1))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# View(oModel)
# 線形トレンド・モデルなので、SSMtrend()を使って以下のように略記できる
oModel2 <- SSModel(
dfSeatbelts$gLogDriver ~ SSMtrend(degree=2, Q=list(matrix(NA), matrix(NA))),
H = matrix(NA)
)
# View(oModel2)
# パラメータ推定
agInit <- c(var(dfSeatbelts$gLogDriver) / 2, 0.001, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.00212052 cat("sigma^2_xi =", as.vector(oFitted$model$Q[1,1,1]), "\n")## sigma^2_xi = 0.01212232 cat("sigma^2_zeta =", as.vector(oFitted$model$Q[2,2,1]), "\n")## sigma^2_zeta = 1.311489e-10 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gAtt = as.vector(oEstimated$att[,1]),
gAtt2 = as.vector(oEstimated$att[,2]),
gError = gLogDriver - gAlphaHat
)
# 描画 ボツ
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g3 <- g3 + geom_line(color = "gray")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g <- g1 + g2 + g3
print(g)
ggsave(paste0(WORKDIR, "/chap7_code6.1.png"), width = 24, height = 12, units = "cm")code 7. 季節要素のあるローカル・レベル・モデル(確定的レベルと確定的季節)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
), nrow=1),
T = matrix(c(
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0
), byrow = T, nrow = 12),
R = matrix(c(1, rep(0, 11)), nrow = 12),
Q = 0,
P1inf = diag(rep(1, 12))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# SSMtrend()とSSMSeasonal()を使って以下のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(0))) +
# SSMseasonal(12, sea.type="dummy"),
# H = matrix(NA)
# )
# パラメータ推定
agInit <- c(var(dfSeatbelts$gLogDriver) / 2, 0.001, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.01758663 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gError = gLogDriver - gAlphaHat1 - gAlphaHat2
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat1))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "gamma[1,t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1.1/2, 1.1/2)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
print(g)
ggsave(paste0(WORKDIR, "/chap7_code7.1.png"), width = 24, height = 12, units = "cm")code 8. 季節要素のあるローカル・レベル・モデル(確率的レベルと確定的季節)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
), nrow=1),
T = matrix(c(
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0
), byrow = T, nrow = 12),
R = matrix(c(1, rep(0, 11)), nrow = 12),
Q = NA,
P1inf = diag(rep(1, 12))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
# # SSMtrend()とSSMSeasonal()を使って以下のように略記できる
# oModel <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
# SSMseasonal(12, sea.type="dummy"),
# H = matrix(NA)
# )
# パラメータ推定
agInit <- c(var(dfSeatbelts$gLogDriver) / 2, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.003513563 cat("sigma^2_xi =", as.vector(oFitted$model$Q[1,1,1]), "\n")## sigma^2_xi = 0.000946002 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gError = gLogDriver - gAlphaHat1 - gAlphaHat2
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat1))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "gamma[1,t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1.1/2, 1.1/2)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
print(g)
ggsave(paste0(WORKDIR, "/chap7_code8.1.png"), width = 24, height = 12, units = "cm")code 9. 季節要素のあるローカル・レベル・モデル(確率的レベルと確率的季節)
# モデルの定義
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom( # 切片項がないので-1と指定
Z = matrix(c(
1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
), nrow=1),
T = matrix(c(
1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0
), byrow = T, nrow = 12),
R = matrix(c(1, rep(0, 11), 0, 1, rep(0, 10)), nrow = 12),
Q = matrix(c(NA, 0, 0, NA), nrow=2),
P1inf = diag(rep(1, 12))
),
H = matrix(NA) # H. NAは推定対象であることを表す
)
## View(oModel)
# # SSMtrend()とSSMSeasonal()を使って以下のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
# SSMseasonal(
# 12, sea.type="dummy",
# Q = matrix(NA), n = nrow(dfSeatbelts)
# ),
# H = matrix(NA)
# )
# ## View(oModel2)
# パラメータ推定
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit)
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.003514111 cat("sigma^2_xi =", as.vector(oFitted$model$Q[1,1,1]), "\n")## sigma^2_xi = 0.0009456102 cat("sigma^2_omega =", as.vector(oFitted$model$Q[2,2,1]), "\n")## sigma^2_omega = 3.646278e-11 # 状態推定(フィルタリング, スムージング)
oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]), # スムージングによる状態推定
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gError = gLogDriver - gAlphaHat1 - gAlphaHat2
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat1))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(6.8, 7.9)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(-1.1/2, 1.1/2)
g3 <- g3 + labs(title = "gamma[1,t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1.1/2, 1.1/2)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
print(g)
ggsave(paste0(WORKDIR, "/chap7_code9.1.png"), width = 24, height = 12, units = "cm")code 10. ARIMA(1,1,1)
# 参考: forecastパッケージ
library(forecast)## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo oArima <- Arima(dfSeatbelts$gLogDriver, order = c(1, 1, 1))
print(oArima)## Series: dfSeatbelts$gLogDriver
## ARIMA(1,1,1)
##
## Coefficients:
## ar1 ma1
## 0.6456 -0.9627
## s.e. 0.0649 0.0223
##
## sigma^2 estimated as 0.01402: log likelihood=136.89
## AIC=-267.78 AICc=-267.65 BIC=-258.02 # 初期モデルの定義
# SSMarima()を使って以下のように略記できる. ひとまずARMAパラメータは0と置く
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ SSMarima(
ar = 0, # ARパラメータ
ma = 0, # MAパラメータ
d = 1, # 差分を1回とる
Q = NA # Q行列は推定対象
),
H = 0 # H行列は0
)
# fitSSMはデフォルトではT行列とR行列を推定できないので、
# 「パラメータを与えるとモデルを返す」関数を定義する
sub_update <- function(par, model) {
# par: 順に{状態撹乱項の分散の対数, ARパラメータ, MAパラメータ}
# model: 現在のモデル
# 与えられた引数に基づき, T行列, R行列, Q行列をつくる
tmp <- SSMarima(
ar = artransform(par[2]), ma = par[3], d = 1, Q = exp(par[1])
)
# モデルを更新する
model["T", states = "arima"] <- tmp$T
model["R", states = "arima"] <- tmp$R
model["Q", states = "arima"] <- tmp$Q
model["P1", states = "arima"] <- tmp$P1
return(model)
}
# パラメータ推定
agInit <- c(0.001, 0.001, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit),
# 「パラメータを与えるとモデルを返す」関数を指定する
updatefn = sub_update,
# 状態撹乱項の分散は、exp(-10)からexp(10)の間で推定する
# ARパラメータは, -1から1のあいだで推定する
# MAパラメータは, -10から10のあいだで推定する
lower = c(-10, -1, -10),
upper = c(+10, +1, +10),
method = "L-BFGS-B"
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$Q), "\n")## sigma^2_epsilon = 0.01285874 cat("phi =", as.vector(oFitted$model$T[2,2,1]), "\n")## phi = 0.6456211 cat("theta =", as.vector(oFitted$model$R[3,1,1]), "\n")## theta = -1.03881Codel 11. 説明変数のあるローカル・レベル・モデル (確定的レベル, 確定的係数)
# 参考: 単回帰
print(summary(lm(gLogDriver ~ gLogPetrol, data=dfSeatbelts)))##
## Call:
## lm(formula = gLogDriver ~ gLogPetrol, data = dfSeatbelts)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37612 -0.09896 -0.01594 0.09077 0.36594
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.87873 0.20889 28.142 < 2e-16 ***
## gLogPetrol -0.67166 0.09173 -7.322 6.74e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1517 on 190 degrees of freedom
## Multiple R-squared: 0.2201, Adjusted R-squared: 0.216
## F-statistic: 53.61 on 1 and 190 DF, p-value: 6.742e-12 mgZ <- array(dim = c(1,2,nrow(dfSeatbelts)))
mgZ[1,1,] <- as.vector(dfSeatbelts$gLogPetrol)
mgZ[1,2,] <- 1
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom(
Z = mgZ,
T = matrix(c(1,0,0,1), nrow = 2),
R = matrix(c(1,0,0,1), nrow = 2),
Q = matrix(c(0, 0, 0, 0), nrow = 2),
P1 = matrix(c(0, 0, 0, 0), nrow = 2),
P1inf = matrix(c(1, 0, 0, 1), nrow = 2)
),
H = matrix(NA)
)
# # 次のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(0))) +
# SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(0)),
# H=matrix(NA)
# )
gInit <- c(var(dfSeatbelts$gLogDriver)/2)
oFitted <- fitSSM(
oModel,
inits = log(gInit),
method = "BFGS"
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.02301349 oEstimated <- KFS(oFitted$model)
cat("beta =", as.vector(oEstimated$alphahat[1,1]), "\n")## beta = -0.6716644 cat("mu =", as.vector(oEstimated$alphahat[1,2]), "\n")## mu = 5.878731 # 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]),
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gBetaX = gAlphaHat1 * gLogPetrol,
gError = gLogDriver - gAlphaHat2 - gBetaX
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(5.3, 7.3)
g2 <- g2 + labs(title = "mu", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gBetaX))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(1.0, 3.0)
g3 <- g3 + labs(title = "beta * X[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1, 1)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
print(g)
ggsave(paste0(WORKDIR, "/chap7_code11.1.png"), width = 24, height = 12, units = "cm")Code 12. 説明変数のあるローカル・レベル・モデル (確率的レベル, 確定的係数)
mgZ <- array(dim = c(1,2,nrow(dfSeatbelts)))
mgZ[1,1,] <- as.vector(dfSeatbelts$gLogPetrol)
mgZ[1,2,] <- 1
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom(
Z = mgZ,
T = matrix(c(1,0,0,1), nrow = 2),
R = matrix(c(1,0,0,1), nrow = 2),
Q = matrix(c(0, 0, 0, NA), nrow = 2),
P1 = matrix(c(0, 0, 0, 0), nrow = 2),
P1inf = matrix(c(1, 0, 0, 1), nrow = 2)
),
H = matrix(NA)
)
# 次のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
# SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(0)),
# H=matrix(NA)
# )
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit),
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.002348964 cat("sigma^2_xi =", as.vector(oFitted$model$Q[2,2,1]), "\n")## sigma^2_xi = 0.01166641 oEstimated <- KFS(oFitted$model, smoothing = c("state", "mean", "disturbance"))
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]),
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gBetaX = gAlphaHat1 * gLogPetrol,
gError = gLogDriver - gAlphaHat2 - gBetaX,
gRstd = as.vector(rstandard(oEstimated, type = "recursive")),
gEtastd = as.vector(rstandard(oEstimated, type = "state")[,2]),
gEpsstd = as.vector(rstandard(oEstimated, type = "pearson"))
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(5.5, 7.5)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gBetaX))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(0.5, 0.7)
g3 <- g3 + labs(title = "beta * X[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1, 1)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
print(g)
ggsave(paste0(WORKDIR, "/chap7_code12.1.png"), width = 24, height = 12, units = "cm")
g <- ggplot(data=dfPlot, aes(x = dMonth, y = gRstd))
g <- g + geom_line()
g <- g + labs(title = "e[t]", x = NULL, y = NULL)
g <- g + theme_bw()
print(g)## Warning: Removed 2 row(s) containing missing values (geom_path).
ggsave(paste0(WORKDIR, "/chap7_code12.2.png"), width = 12, height = 12, units = "cm")## Warning: Removed 2 row(s) containing missing values (geom_path). png(paste0(WORKDIR, "/chap7_code12.3.png"), width = 12, height = 12, res = 72, units = "cm")
acf(dfPlot$gRstd[!is.na(dfPlot$gRstd)], main = NULL)
dev.off()## png
## 2 g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gEpsstd))
g1 <- g1 + geom_hline(yintercept = 0, color="gray")
g1 <- g1 + geom_line()
g1 <- g1 + labs(title = "standardized hat(epsilon)[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gEtastd))
g2 <- g2 + geom_hline(yintercept = 0, color="gray")
g2 <- g2 + geom_line()
g2 <- g2 + labs(title = "standardized hat(eta)[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g <- g1 + g2
print(g)
ggsave(paste0(WORKDIR, "/chap7_code12.4.png"), width = 24, height = 12, units = "cm")Code 13. 説明変数のあるローカル・レベル・モデル (確率的レベル, 確率的係数)
mgZ <- array(dim = c(1,2,nrow(dfSeatbelts)))
mgZ[1,1,] <- as.vector(dfSeatbelts$gLogPetrol)
mgZ[1,2,] <- 1
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ -1 + SSMcustom(
Z = mgZ,
T = matrix(c(1,0,0,1), nrow = 2),
R = matrix(c(1,0,0,1), nrow = 2),
Q = matrix(c(NA, 0, 0, NA), nrow = 2),
P1 = matrix(c(0, 0, 0, 0), nrow = 2),
P1inf = matrix(c(1, 0, 0, 1), nrow = 2)
),
H = matrix(NA)
)
# 次のように略記できる
# oModel2 <- SSModel(
# dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
# SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(NA)),
# H=matrix(NA)
# )
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001, 0.001)
oFitted <- fitSSM(
oModel,
inits = log(agInit),
)
cat("sigma^2_epsilon =", as.vector(oFitted$model$H), "\n")## sigma^2_epsilon = 0.002356589 cat("sigma^2_xi =", as.vector(oFitted$model$Q[2,2,1]), "\n")## sigma^2_xi = 0.01097267 cat("sigma^2_tau =", as.vector(oFitted$model$Q[1,1,1]), "\n")## sigma^2_tau = 0.0001308724 oEstimated <- KFS(oFitted$model)
# 描画準備
dfPlot <- dfSeatbelts %>%
mutate(
gAlphaHat1 = as.vector(oEstimated$alphahat[,1]),
gAlphaHat2 = as.vector(oEstimated$alphahat[,2]),
gBetaX = gAlphaHat1 * gLogPetrol,
gError = gLogDriver - gAlphaHat2 - gBetaX
)
# 描画
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogDriver))
g1 <- g1 + geom_line()
g1 <- g1 + ylim(6.8, 7.9)
g1 <- g1 + labs(title = "y[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat2))
g2 <- g2 + geom_line(color = "red")
g2 <- g2 + ylim(5.5, 7.5)
g2 <- g2 + labs(title = "mu[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
g3 <- ggplot(data=dfPlot, aes(x = dMonth, y = gBetaX))
g3 <- g3 + geom_line(color = "red")
g3 <- g3 + ylim(0.5, 0.7)
g3 <- g3 + labs(title = "beta[t] * X[t]", x = NULL, y = NULL)
g3 <- g3 + theme_bw()
g3
g4 <- ggplot(data=dfPlot, aes(x = dMonth, y = gError))
g4 <- g4 + geom_line(color = "gray")
g4 <- g4 + ylim(-1, 1)
g4 <- g4 + labs(title = "epsilon[t]", x = NULL, y = NULL)
g4 <- g4 + theme_bw()
g <- g1 + g2 + g3 + g4
ggsave(paste0(WORKDIR, "/chap7_code13.1.png"), width = 24, height = 12, units = "cm")
g1 <- ggplot(data=dfPlot, aes(x = dMonth, y = gLogPetrol))
g1 <- g1 + geom_line(color="red")
g1 <- g1 + labs(title = "X[t]", x = NULL, y = NULL)
g1 <- g1 + theme_bw()
print(g1)
g2 <- ggplot(data=dfPlot, aes(x = dMonth, y = gAlphaHat1))
g2 <- g2 + geom_line(color="red")
g2 <- g2 + labs(title = "beta[t]", x = NULL, y = NULL)
g2 <- g2 + theme_bw()
print(g2)
g <- g1 + g2
print(g)
ggsave(paste0(WORKDIR, "/chap7_code13.2.png"), width = 24, height = 12, units = "cm")Code 14. モデル比較
sub.AIC <- function(gLogLik, nStateVar, nHyperParam){
# purpose: get AIC
# See Section 2.1.
# args: gLogLik: (numeric scalar) log likelihood
# nStateVar: (integer scalar) number of state variables
# nHyperParam: (integer scalar) number of hyper params
# return: (numeric scalar) AIC
- 2 * gLogLik + 2 * ( nStateVar + nHyperParam )
}
# 説明変数のあるローカル・レベル・モデル (確定的レベル, 確定的係数)
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(0))) +
SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(0)),
H=matrix(NA)
)
gInit <- c(var(dfSeatbelts$gLogDriver)/2)
oFitted1 <- fitSSM(
oModel,
inits = log(gInit),
method = "BFGS"
)
# 説明変数のあるローカル・レベル・モデル (確率的レベル, 確定的係数)
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(0)),
H=matrix(NA)
)
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001)
oFitted2 <- fitSSM(
oModel,
inits = log(agInit),
)
# 説明変数のあるローカル・レベル・モデル (確率的レベル, 確率的係数)
oModel <- SSModel(
dfSeatbelts$gLogDriver ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
SSMregression(~ -1 + dfSeatbelts$gLogPetrol, Q=matrix(NA)),
H=matrix(NA)
)
agInit <- c(var(dfSeatbelts$gLogDriver)/2, 0.001, 0.001)
oFitted3 <- fitSSM(
oModel,
inits = log(agInit),
)
print(sub.AIC(logLik(oFitted1$model), 2, 1))## [1] -165.1565 print(sub.AIC(logLik(oFitted2$model), 2, 2))## [1] -239.9243 print(sub.AIC(logLik(oFitted3$model), 2, 3))## [1] -237.9271以上