このページでは、株式会社インサイト・ファクトリー 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)