dlm_code

##

## dlmによる状態空間モデル

##

library(dlm)

library(ggplot2)

library(gridExtra)

font <- "IPAexGothic" # フォントは適宜指定してください

options(digits = 3)

## (1) dlmによるローカルレベルモデル

# 松浦(2016)のデータを使用

# カレントディレクトリにデータがなければ、GitHubから取得

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

# グラフ描画

df.data <- data.frame(Time = seq_along(Y),

Value = Y)

p0 <- ggplot(df.data) +

geom_line(aes(x = Time, y = Value)) +

labs(x = "時間（四半期）", y = "販売個数（千個）") +

theme_classic(base_family = font)

print(p0)

# ローカルレベルモデル

# 観測ノイズとシステムノイズの分散をともに1とした場合

dlm.mod1<- dlmModPoly(order = 1, dV = 1, dW = 1,

m0 = 0, C0 = 1e+07)

dlm.flt1 <- dlmFilter(Y, dlm.mod1)

# フィルタリングされた値の最初の部分を表示

head(dlm.flt1$m)

# 最初の要素を除く

head(dropFirst(dlm.flt1$m))

# 観測ノイズの分散を1、システムノイズの分散を5とした場合

dlm.mod2<- dlmModPoly(order = 1, dV = 1, dW = 5,

m0 = 0, C0 = 1e+07)

dlm.flt2 <- dlmFilter(Y, dlm.mod2)

dlm.flt2$sd <- sqrt(unlist(dlmSvd2var(dlm.flt2$U.C, dlm.flt2$D.C)))

# 最尤推定

dlm.build3 <- function(par) {

dlmModPoly(order = 1,

dV = exp(par[1]), dW = exp(par[2]))

}

dlm.mle3 <- dlmMLE(Y, parm = c(0, 0),

build = dlm.build3)

dlm.mod3 <- dlm.build3(dlm.mle3$par)

# ノイズの誤差分散の値

V(dlm.mod3)

W(dlm.mod3)

# フィルタリング

dlm.flt3 <- dlmFilter(Y, dlm.mod3)

head(dropFirst(dlm.flt3$m))

# 平滑化

dlm.smt3 <- dlmSmooth(Y, dlm.mod3)

head(dropFirst(dlm.smt3$s))

# 各モデルの対数尤度（の符号を反転させた値）

dlmLL(Y, dlm.mod1)

dlmLL(Y, dlm.mod2)

dlmLL(Y, dlm.mod3)

# 誤差標準偏差

dlm.flt1$sd <- sqrt(unlist(dlmSvd2var(dlm.flt1$U.C, dlm.flt1$D.C)))

dlm.flt3$sd <- sqrt(unlist(dlmSvd2var(dlm.flt3$U.C, dlm.flt3$D.C)))

dlm.smt3$sd <- sqrt(unlist(dlmSvd2var(dlm.smt3$U.S, dlm.smt3$D.S)))

# グラフ

df1 <- data.frame(Time = rep(seq_along(Y), 6),

Value = c(Y, dropFirst(dlm.flt1$m),

Y, dropFirst(dlm.flt3$m),

Y, dropFirst(dlm.smt3$s)),

Type = factor(c(rep("観測値", length(Y)),

rep("フィルタリングした値", length(Y)),

rep("観測値", length(Y)),

rep("フィルタリングした値", length(Y)),

rep("観測値", length(Y)),

rep("平滑化した値", length(Y))),

levels = c("観測値", "フィルタリングした値",

"平滑化した値")),

Facet = factor(c(rep("A", length(Y) * 2),

rep("B", length(Y) * 2),

rep("C", length(Y) * 2))))

df1.ci <- data.frame(Time = rep(seq_along(Y), 6),

Lower = c(dropFirst(qnorm(0.025, dlm.flt1$m, dlm.flt1$sd)),

dropFirst(qnorm(0.025, dlm.flt3$m, dlm.flt3$sd)),

dropFirst(qnorm(0.025, dlm.smt3$s, dlm.smt3$sd))),

Upper = c(dropFirst(qnorm(0.975, dlm.flt1$m, dlm.flt1$sd)),

dropFirst(qnorm(0.975, dlm.flt3$m, dlm.flt3$sd)),

dropFirst(qnorm(0.975, dlm.smt3$s, dlm.smt3$sd))),

Facet = factor(rep(c(rep("A", length(Y)),

rep("B", length(Y)),

rep("C", length(Y))), 2)))

p1 <- ggplot(df1) +

geom_line(aes(x = Time, y = Value, color = Type,

size = 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)) +

labs(x = "時間（四半期）", y = "販売個数（千個）") +

geom_ribbon(data = df1.ci,

aes(x = Time, ymin = Lower, ymax = Upper),

alpha = 0.25) +

theme_classic(base_family = font) +

facet_wrap(~ Facet, ncol = 3, scales = "free_x") +

theme(strip.background = element_blank(),

strip.text = element_text(size = 9, hjust = 0.1),

legend.title = element_blank(), legend.position = "bottom")

print(p1)

# 2階差分のモデル

dlm.build3.2 <- function(par) {

mod <- dlmModPoly(order = 2, dV = exp(par[1]),

dW = c(0, exp(par[2])))

}

dlm.mle3.2 <- dlmMLE(Y, parm = c(0, 0),

build = dlm.build3.2)

# 尤度

dlmLL(Y, dlm.build3.2(dlm.mle3.2$par))

## (2) dlmによる季節調整モデル

dlm.build4 <- function(par) {

mod <- dlmModPoly(order = 1) +

dlmModSeas(frequency = 4)

V(mod) <- exp(par[1])

diag(W(mod))[1:2] <- exp(par[2:3])

m0(mod) <- c(20, 0, 0, 0)

return(mod)

}

# 最尤推定

dlm.mle4 <- dlmMLE(Y, parm = rep(0, 3),

build = dlm.build4)

dlm.mod4 <- dlm.build4(dlm.mle4$par)

# 係数行列

FF(dlm.mod4)

GG(dlm.mod4)

# ノイズの誤差分散

V(dlm.mod4)

W(dlm.mod4)

# 対数尤度（の符号を反転させた値）

dlmLL(Y, dlm.mod4)

# 予測

dlm.fil4 <- dlmFilter(Y, dlm.mod4)

dlm.for4 <- dlmForecast(dlm.fil4, nAhead = 8)

# 予測値の最初の部分

head(dlm.for4$f, 3)

head(dlm.for4$a, 3)

head(unlist(dlm.for4$Q))

# 平滑化

dlm.smo4 <- dlmSmooth(dlm.fil4)

head(dropFirst(dlm.smo4$s), 3)

# 誤差の標準偏差

var.s <- dlmSvd2var(dlm.smo4$U.S, dlm.smo4$D.S)

sd.s <- sapply(1:length(var.s), function(i) sqrt(diag(var.s[[i]])[1:2]))

# グラフ描画

obs.exp <- dlm.for4$f

lev.exp <- dlm.for4$a[, 1]

ses.exp <- dlm.for4$a[, 2]

stdev <- sqrt(unlist(dlm.for4$Q))

df2 <- data.frame(Time = c(rep(seq_along(Y), 2),

length(Y) + 0:(nrow(obs.exp))),

Value = c(Y, dropFirst(dlm.smo4$s[, 1]), Y[length(Y)], obs.exp),

Type = factor(c(rep("観測値", length(Y)),

rep("水準成分", length(Y)),

rep("予測値", nrow(obs.exp) + 1)),

levels = c("観測値", "水準成分",

"予測値")))

df2.ci1 <- data.frame(Time = seq_along(Y),

Lower = dropFirst(qnorm(0.025, dlm.smo4$s[, 1], sd.s[1, ])),

Upper = dropFirst(qnorm(0.975, dlm.smo4$s[, 1], sd.s[1, ])))

df2.ci2 <- data.frame(Time = length(Y) + 1:nrow(obs.exp),

Lower = qnorm(0.025, obs.exp, stdev),

Upper = qnorm(0.975, obs.exp, stdev))

p2 <- ggplot(df2) +

geom_line(aes(x = Time, y = Value, color = Type,

size = Type, alpha = Type)) +

scale_color_manual(values = c("black", "gray10",

"gray40")) +

scale_size_manual(values = c(0.33, 1.0, 1.0)) +

scale_alpha_manual(values = c(1, 0.6, 0.6)) +

geom_ribbon(data = df2.ci1, aes(x = Time, ymin = Lower, ymax = Upper),

fill = "gray50", alpha = 0.5) +

geom_ribbon(data = df2.ci2, 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))

df3 <- data.frame(Time = seq_along(Y),

Value = dropFirst(dlm.smo4$s[, 2]))

df3.ci <- data.frame(Time = seq_along(Y),

Lower = dropFirst(qnorm(0.025, dlm.smo4$s[, 2], sd.s[2, ])),

Upper = dropFirst(qnorm(0.975, dlm.smo4$s[, 2], sd.s[2, ])))

p3 <- ggplot(df3) +

geom_line(aes(x = Time, y = Value), color = "gray25", size = 1) +

geom_ribbon(data = df3.ci,

aes(x = Time, ymin = Lower, ymax = Upper),

fill = "gray50", alpha = 0.5) +

labs(title = "B", x = "時間（四半期）", y = "販売個数（千個）") +

theme_classic(base_family = font)

grid.arrange(p2, p3, ncol = 2)

# ベイズ推定

# 季節調整モデル

set.seed(1)

dlm.fit5 <- dlmGibbsDIG(Y, dlmModPoly(order = 1) +

dlmModSeas(frequency = 4),

ind = 1:2,

shape.y = 0.01, rate.y = 0.01,

shape.theta = 0.01,

rate.theta = 0.01,

n.sample = 1010, thin = 29)

# 事後平均

mcmcMean(dlm.fit5$dV[-(1:10)])

mcmcMean(dlm.fit5$dW[-(1:10), ])

# 軌跡

p5 <- ggplot(data.frame(Iter = rep(1:1010, 3),

Value = c(dlm.fit5$dV, dlm.fit5$dW[, 1], dlm.fit5$dW[, 2]),

Var = rep(c("dV", "dW1", "dW2"), each = 1010))) +

geom_line(aes(x = Iter, y = Value)) +

facet_wrap(~ Var, ncol = 1, scales = "free_y") +

theme_classic(base_family = font)

print(p5)