KFAS
##
## KFASによる状態空間モデル
##
library(KFAS)
library(ggplot2)
library(gridExtra)
font <- "IPAexGothic" # フォントは適宜指定してください
options(digits = 3)
# (1) KFASによる季節調整モデル
# 松浦(2016)のデータを使用
data.file <- "data-season.txt"
if (file.exists(data.file)) {
data <- read.csv(data.file)
} else {
data.url <- "https://raw.githubusercontent.com/iwanami-datascience/vol1/master/matsuura/example2/input/data-season.txt"
data <- read.csv(url(data.url))
}
Y <- data$Y
kfas.mod1 <- SSModel(Y ~ SSMtrend(degree = 1, Q = NA,
a1 = 20) +
SSMseasonal(period = 4,
sea.type = "dummy",
Q = NA), H = NA)
# 係数行列
kfas.mod1$Z
kfas.mod1$T
# パラメータを最尤推定
kfas.fit1 <- fitSSM(kfas.mod1, inits = rep(0, 3))
# 観測ノイズの分散
kfas.fit1$model$H
# システムノイズの分散
kfas.fit1$model$Q
# 尤度
logLik(kfas.fit1$model)
# カルマンフィルタを適用
kfas.fil1 <- KFS(kfas.fit1$model, filtering = "state",
smoothing = c("state"))
# 状態の1期先予測
head(kfas.fil1$a, 3)
# 平滑化した状態の値
head(kfas.fil1$alphahat, 3)
# 予測
kfas.fct1 <- predict(kfas.fit1$model, n.ahead = 8,
interval = "prediction", level = 0.95,
type = "response", nsim = 1000)
obs.exp <- kfas.fct1[, 1]
lower <- kfas.fct1[, 2]
upper <- kfas.fct1[, 3]
# グラフ
df1 <- data.frame(Time = c(rep(seq_along(Y), 2),
length(Y) + seq_along(obs.exp)),
Value = c(Y, kfas.fil1$alphahat[, 1], obs.exp),
Lower = c(rep(NA, length(Y) * 2), lower),
Upper = c(rep(NA, length(Y) * 2), upper),
Type = factor(c(rep("観測値", length(Y)),
rep("平滑化した値", length(Y)),
rep("予測値", length(obs.exp))),
levels = c("観測値", "平滑化した値",
"予測値")))
p1 <- ggplot(df1) +
geom_line(aes(x = Time, y = Value, color = Type, alpha = Type)) +
scale_color_manual(values = c("black", "gray10", "gray60")) +
scale_size_manual(values = c(0.5, 1.5, 1.5)) +
scale_alpha_manual(values = c(1, 0.6, 0.6)) +
geom_ribbon(aes(x = Time, ymin = Lower, ymax = Upper),
fill = "gray50", alpha = 0.5) +
labs(title = "A", x = "時間（四半期）", y = "販売個数（千個）") +
theme_classic(base_family = font) +
theme(legend.title = element_blank(),
legend.justification = c(0, 1),
legend.position = c(0.1, 1))
print(p1)
# AR(3)モデル
mod.ar <- SSModel(Y ~ 0 +
SSMarima(ar = rep(0, 3), Q = NA),
H = 0)
update.ar <- function(pars, mod) {
a <- artransform(pars[1:3])
sigma2 <- exp(pars[4])
mod <- SSModel(Y ~ 0 +
SSMarima(ar = a, Q = sigma2),
H = 0)
return(mod)
}
fit.ar <- fitSSM(mod.ar, inits = rep(1, 4),
updatefn = update.ar,
method = "L-BFGS-B")
# 結果
print(fit.ar$model$T, digits = 3)
print(fit.ar$model$Q, digits = 3)
# (2) ポアソン分布のモデル
# 久保(2012)のデータを使用
# 1次元の空間データを時系列と見なして解析する
if (file.exists(data.file)) {
data.file <- "Y.RData"
load(data.file)
} else {
data.url <- "http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/spatial/Y.RData"
load(url(data.url))
}
df2 <- data.frame(Time = seq_along(Y), Y = Y)
p2 <- ggplot(df2) +
geom_line(aes(x = Time, y = Y)) +
theme_classic(base_family = font)
print(p2)
# ローカルレベルモデル
kfas.mod2 <- SSModel(Y ~ SSMtrend(degree = 1, Q = NA),
distribution = "poisson")
kfas.fit2 <- fitSSM(kfas.mod2, inits = c(0),
method = "BFGS")
# ノイズ分散
kfas.fit2$model$Q
kfas.fit2$model$H
# 観測値の予測
kfas.fil2 <- KFS(kfas.fit2$model, smoothing = "mean")
fitted(kfas.fil2)
# 平滑化した値
kfas.pre2 <- predict(kfas.fil2$model, interval = "confidence",
level = 0.95, type = "response")
df2 <- data.frame(Time = rep(seq_along(Y), 2),
Y = c(Y, kfas.pre2[, "fit"]),
Type = factor(rep(c("観測値", "平滑化した値"),
each = length(Y)),
levels = c("観測値", "平滑化した値")))
df2.ci <- data.frame(Time = seq_along(Y),
Lower = kfas.pre2[, "lwr"],
Upper = kfas.pre2[, "upr"])
p2 <- ggplot(df2) +
geom_line(aes(x = Time, y = Y, color = Type,
size = Type, alpha = Type)) +
scale_color_manual(values = c("black", "gray20")) +
scale_size_manual(values = c(0.5, 1.5)) +
scale_alpha_manual(values = c(1, 0.6)) +
geom_ribbon(data = df2.ci,
aes(x = Time, ymin = Lower, ymax = Upper),
fill = "gray50", alpha = 0.5) +
labs(title = "A") +
theme_classic(base_family = font) +
theme(legend.title = element_blank(),
legend.text = element_text(size = 8),
legend.justification = c(1, 1),
legend.position = c(1, 1.07),
legend.background = element_rect(fill = "transparent"))
print(p2)
# 2階差分のモデル
kfas.mod2.2 <- SSModel(Y ~ SSMtrend(degree = 2,
Q = list(0, NA)),
distribution = "poisson")
kfas.fit2.2 <- fitSSM(kfas.mod2.2, inits = c(0),
method = "BFGS")
# システムノイズの分散
kfas.fit2.2$model$Q
# 観測値の予測
kfas.flt22 <- KFS(kfas.fit2.2$model, smoothing = "mean")
fitted(kfas.flt2.2)
kfas.pre22 <- predict(kfas.flt22$model, interval = "confidence",
level = 0.95, type = "response")
df22 <- data.frame(Time = rep(seq_along(Y), 2),
Y = c(Y, kfas.pre22[, "fit"]),
Type = factor(rep(c("観測値", "平滑化した値"),
each = length(Y)),
levels = c("観測値", "平滑化した値")))
df22.ci <- data.frame(Time = seq_along(Y),
Lower = kfas.pre22[, "lwr"],
Upper = kfas.pre22[, "upr"])
p22 <- ggplot(df22) +
geom_line(aes(x = Time, y = Y, color = Type, size = Type)) +
scale_color_manual(values = c("black", "gray50")) +
geom_ribbon(data = df22.ci,
aes(x = Time, ymin = Lower, ymax = Upper),
fill = "gray50", alpha = 0.5) +
scale_size_manual(values = c(0.5, 1)) +
theme_classic(base_family = font) +
theme(legend.title = element_blank())
print(p22)