読者です 読者をやめる 読者になる 読者になる

StatsFragments

Python, R, Rust, 統計, 機械学習とか

R で 状態空間モデル: {dlm} で単変量モデルの状態空間表現を利用する

時系列分析 R

こちらの続きで、状態空間時系列分析入門の第8章の内容。

sinhrks.hatenablog.com

状態空間時系列分析入門

状態空間時系列分析入門

  • 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
  • 出版社/メーカー: シーエーピー出版
  • 発売日: 2008/09
  • メディア: 単行本
  • 購入: 2人 クリック: 4回
  • この商品を含むブログを見る

状態空間モデルでは、トレンド、季節効果、説明変数効果などさまざまな要因を考慮することができる。また これらの要因を単一の状態空間表現で統一的に扱える。時刻 { t } における観測値 { y_t } と状態 { a_t } は以下の式で書ける (前回と同じ)。

{ y_t = z'_t a_t + \epsilon_{t} \ \ \ \ \ \ \epsilon_t \sim NID(0, \sigma^{2}_{ \epsilon } ) }

{ a_{ t+1 } = T_t a_t + R_t \eta_t \ \ \ \ \ \ \eta_t \sim NID(0, Q_t) }

この式であらわされたモデルを {dlm} で利用する際にどうすればよいのかをまとめたい。

vignette にあるとおり、{dlm} ではモデルを dlm::dlmModPoly などのモデル記述関数の組み合わせで定義できる。これらの関数は便利だが、上の式との関係は外側からはよくわからない。

上式から {dlm} のモデルを作成する際には dlm::dlm を直接使うのがよさそうだ。dlm::dlm の引数は上式にほぼ対応しているため、式表現されたモデルをそのまま実装できる。また、モデル記述関数も内部的には dlm::dlm を使っているため、これら関数の理解にも役立つ。dlm::dlm の引数と上式の対応は前回書いたとおり。

  • FF: { z'_t } に対応。
  • V: { \sigma^{2} } に対応。
  • GG: { T_t } に対応。
  • W: { R_t Q_t } に対応。
  • m0: { a_1 } の事前分布の平均。
  • C0: { a_1 } の事前分布の分散。
  • J...: 各パラメータが時間変化するとき利用。

参考 サイト 「状態空間時系列分析入門」をRで再現する はこのテキストを {dlm}, {KFAS} で再現している。こちらは dlm::dlm ではなく dlm::dlmModPoly などの関数を使ってモデル作成している。

準備

まずデータを読み込む。データは著者のサポートサイト、An Introduction to State Space Time Series Analysis から入手できる。

library(dlm)
# install_github('sinhrks/ggfortify')
library(ggfortify)

ukdrivers <- read.table('UKdriversKSI.txt', skip = 1)
ukdrivers <- ts(ukdrivers[[1]], start = c(1969, 1), frequency = 12)
ukdrivers <- log(ukdrivers)

以降、テキストに記載のモデル順に dlm::dlm から利用する場合のプログラムを列挙する。値 / 推定するパラメータを変更したい場合は buildObj 内の対応する要素を適当に修正すればよい。

1. ローカルレベルモデル (確定的レベル)

それぞれのモデルについて 状態 { a_t } として何が考慮されているかを書いておく。

  • { a_t = u }: { u } は確定的レベル (経時変化しない)。

補足 それぞれのモデルの式表現を書こうとしたのだが、はてなマークダウン記法だと TeX 記法中 の "&" がうまく処理できないようなので諦めた。詳細はテキストと見比べてください。

y <- ukdrivers

parm<- log(var(y))
buildObj <- function(parm) {
  dlm(FF = matrix(1, 1, 1),
      V = matrix(exp(parm), 1, 1),
      GG = matrix(1, 1, 1),
      W = matrix(0, 1, 1),
      m0 = 0,
      C0 = matrix(1e+07, 1, 1))
}
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1]
# [1,]    1
# 
# $V
#            [,1]
# [1,] 0.02935256
# 
# $GG
#      [,1]
# [1,]    1
# 
# $W
#      [,1]
# [1,]    0
# 
# $m0
# [1] 0
# 
# $C0
#       [,1]
# [1,] 1e+07

p <- autoplot(dlmFilter(y, dlmModObj), fitted.colour = 'red', fitted.linetype = 'dashed')
autoplot(dlmSmooth(y, dlmModObj), colour = 'blue', p = p)

f:id:sinhrks:20150524153631p:plain

赤破線がカルマンフィルタ、青破線が平滑化した値。以降もプロットするが、確率的状態を使うと 実データにほぼフィットするため、グラフからでは差がよくわからない。

2. ローカル線形トレンドモデル (確率的レベルと確率的傾き)

  • { a_t = ( u_t, v_t )^T }: { u_t } は確率的レベル (経時変化する)、{ v_t } は確率的傾き(経時変化する)。
parm <- log(c(var(y), 0.001, 0.001))
buildObj <- function(parm) {
  dlm(FF = matrix(c(1, 0), 1, 2),
      V = matrix(exp(parm[1]), 1, 1),
      GG = matrix(c(1, 1,
                    0, 1), 2, 2, byrow = TRUE),
      W = matrix(c(exp(parm[2]), 0,
                   0, exp(parm[3])), 2, 2, byrow = TRUE),
      m0 = c(0, 0),
      C0 = matrix(c(1e+07, 0,
                    0, 1e+07), 2, 2, byrow = TRUE))
}
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1] [,2]
# [1,]    1    0
# 
# $V
#             [,1]
# [1,] 0.002118549
# 
# $GG
#      [,1] [,2]
# [1,]    1    1
# [2,]    0    1
# 
# $W
#            [,1]        [,2]
# [1,] 0.01212741 0.00000e+00
# [2,] 0.00000000 1.92431e-10
# 
# $m0
# [1] 0 0
# 
# $C0
#       [,1]  [,2]
# [1,] 1e+07 0e+00
# [2,] 0e+00 1e+07

p <- autoplot(dlmFilter(y, dlmModObj), fitted.colour = 'red', fitted.linetype = 'dashed')
autoplot(dlmSmooth(y, dlmModObj), colour = 'blue', p = p)

f:id:sinhrks:20150524153644p:plain

3. 季節要素のあるローカルレベルモデル (確率的レベルと確率的季節)

  • { a_t = ( u_t, \gamma_{ 1, t }, \gamma_{ 2, t }, \gamma_{ 3, t } )^T }: { u_t } は確率的レベル、{ \gamma } は 四半期単位の確率的季節のダミー項。
ukinflation <- read.table('UKinflation.txt', skip = 1)
ukinflation <- ts(ukinflation[[1]], start = c(1950, 1), frequency = 4)
y <- ukinflation

parm <- log(c(var(y), 0.00001, 0.00001))
buildObj <- function(parm) {
  dlm(FF = matrix(c(1, 1, 0, 0), 1, 4),
      V = matrix(c(exp(parm[1])), 1, 1),
      GG = matrix(c(1,  0,  0,  0,
                    0, -1, -1, -1,
                    0,  1,  0,  0,
                    0,  0,  1,  0),
                  4, 4, byrow = TRUE),
      W = matrix(c(exp(parm[2]),            0, 0, 0,
                   0,            exp(parm[3]), 0, 0,
                   0,                       0, 0, 0,
                   0,                       0, 0, 0),
                 4, 4, byrow = TRUE),
      m0 = c(0, 0, 0, 0),
      C0 = matrix(c(1e+07,     0,     0,     0,
                    0,     1e+07,     0,     0,
                    0,         0, 1e+07,     0,
                    0,         0,     0, 1e+07),
                  4, 4, byrow = TRUE))
}
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1] [,2] [,3] [,4]
# [1,]    1    1    0    0
# 
# $V
#             [,1]
# [1,] 3.37127e-05
# 
# $GG
#      [,1] [,2] [,3] [,4]
# [1,]    1    0    0    0
# [2,]    0   -1   -1   -1
# [3,]    0    1    0    0
# [4,]    0    0    1    0
# 
# $W
#              [,1]         [,2] [,3] [,4]
# [1,] 2.124158e-05 0.000000e+00    0    0
# [2,] 0.000000e+00 4.345176e-07    0    0
# [3,] 0.000000e+00 0.000000e+00    0    0
# [4,] 0.000000e+00 0.000000e+00    0    0
# 
# $m0
# [1] 0 0 0 0
# 
# $C0
#       [,1]  [,2]  [,3]  [,4]
# [1,] 1e+07 0e+00 0e+00 0e+00
# [2,] 0e+00 1e+07 0e+00 0e+00
# [3,] 0e+00 0e+00 1e+07 0e+00
# [4,] 0e+00 0e+00 0e+00 1e+07

p <- autoplot(dlmFilter(y, dlmModObj), fitted.colour = 'red', fitted.linetype = 'dashed')
autoplot(dlmSmooth(y, dlmModObj), colour = 'blue', p = p)

f:id:sinhrks:20150524153704p:plain

4. 説明変数のあるローカルレベルモデル (確率的レベル)

  • { a_t = ( u_t, \beta_t )^T }: { u_t } は確率的レベル、{ \beta_t } は説明変数の係数。

説明変数は dlm::dlm へ引数 X として渡す。

ukpetrol <- read.table('logUKpetrolprice.txt', skip = 1)
ukpetrol <- ts(ukpetrol, start = start(ukdrivers), frequency = frequency(ukdrivers))

y <- ukdrivers
x <- ukpetrol

parm <- c(var(y), 0.001)  

buildObj <- function(parm) {
  dlm(FF = matrix(c(1, 0), 1, 2),
      V = matrix(exp(parm[1]), 1, 1),
      GG = matrix(c(1, 0,
                    0, 1),
                  2, 2, byrow = TRUE),
      W = matrix(c(exp(parm[2]), 0,
                   0,            0),
                 2, 2, byrow = TRUE),
      m0 = c(0, 0),
      C0 = matrix(c(1e+07,     0,
                    0,     1e+07),
                  2, 2, byrow = TRUE),
      JFF = matrix(c(0, 1), 1, 2),
      X = as.matrix(x))
}
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1] [,2]
# [1,]    1    1
# 
# $V
#             [,1]
# [1,] 0.002347965
# 
# $GG
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1
# 
# $W
#            [,1] [,2]
# [1,] 0.01166743    0
# [2,] 0.00000000    0
# 
# $JFF
#      [,1] [,2]
# [1,]    0    1
# 
# $X
#      V1    
# [1,] -2.273
# [2,] -2.279
# [3,] ...   
# 
# $m0
# [1] 0 0
# 
# $C0
#       [,1]  [,2]
# [1,] 1e+07 0e+00
# [2,] 0e+00 1e+07

filtered <- dlmFilter(y, dlmModObj)$m
smoothed <- dlmSmooth(y, dlmModObj)$s
p <- autoplot(y)
p <- autoplot(filtered[, 1] + filtered[, 2] * x, colour = 'red', linetype = 'dashed', p = p)
autoplot(smoothed[, 1] + smoothed[, 2] * x, colour = 'blue', p = p)

f:id:sinhrks:20150524153715p:plain

5. 干渉変数のあるローカルレベルモデル (確率的レベル)

  • { a_t = ( u_t, \lambda_t )^T }: { u_t } は確率的レベル、{ \lambda_t } は干渉変数の係数。

変数名は異なるが、式は説明変数のあるローカルレベルモデルと同一。

ukseats <- c(rep(0, (1982 - 1968) * 12 + 1), rep(1, (1984 - 1982) * 12 - 1))
ukseats <- ts(ukseats, start = start(ukdrivers), frequency = frequency(ukdrivers))

x <- ukseats
parm <- c(var(y), 0.001)  

buildObj <- function(parm) {
  dlm(FF = matrix(c(1, 0), 1, 2),
      V = matrix(exp(parm[1]), 1, 1),
      GG = matrix(c(1, 0,
                    0, 1),
                  2, 2, byrow = TRUE),
      W = matrix(c(exp(parm[2]), 0,
                   0,            0),
                 2, 2, byrow = TRUE),
      m0 = c(0, 0),
      C0 = matrix(c(1e+07,     0,
                    0,     1e+07),
                  2, 2, byrow = TRUE),
      JFF = matrix(c(0, 1), 1, 2),
      X = as.matrix(x))
}

dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1] [,2]
# [1,]    1    0
# 
# $V
#             [,1]
# [1,] 0.002692686
# 
# $GG
#      [,1] [,2]
# [1,]    1    0
# [2,]    0    1
# 
# $W
#            [,1] [,2]
# [1,] 0.01041175    0
# [2,] 0.00000000    0
# 
# $JFF
#      [,1] [,2]
# [1,]    0    1
# 
# $X
#      [,1]
# [1,] 0   
# [2,] 0   
# [3,] ... 
# 
# $m0
# [1] 0 0
# 
# $C0
#       [,1]  [,2]
# [1,] 1e+07 0e+00
# [2,] 0e+00 1e+07

filtered <- dlmFilter(y, dlmModObj)$m
smoothed <- dlmSmooth(y, dlmModObj)$s
p <- autoplot(y)
p <- autoplot(filtered[, 1] + filtered[, 2] * x, colour = 'red', linetype = 'dashed', p = p)
autoplot(smoothed[, 1] + smoothed[, 2] * x, colour = 'blue', p = p)

f:id:sinhrks:20150524153729p:plain

6. 季節要素と干渉変数のあるローカルレベルモデル (確率的レベルと確率的季節)

  • { a_t = ( u_t, \gamma_{ 1, t }, \gamma_{ 2, t }, \gamma_{ 3, t }, \lambda_t )^T }: { u_t } は確率的レベル、{ \gamma } は確率的季節、{ \lambda_t } は干渉変数の係数。

補足 テキスト記載の数式は3項の季節要素、説明変数、干渉変数を含む 6 状態が記載されている。季節要素が3項のため英国インフレーションモデルについて記載していると思われるが、このモデルには説明変数がないため状態数があわない。そのため、テキスト中の式を英国インフレーションモデルに合うよう書き換えている。

y <- ukinflation

ukpulse <- rep(0, length.out = length(ukinflation))
ukpulse[4*(1975-1950)+2] <- 1
ukpulse[4*(1979-1950)+3] <- 1
ukpulse <- ts(ukpulse, start = start(ukinflation), frequency = frequency(ukinflation))
x <- ukpulse

parm <- log(c(var(y), 0.00001, 0.00001))
buildObj <- function(parm) {
  dlm(FF = matrix(c(1, 1, 0, 0, 1), 1, 5),
      V = matrix(c(exp(parm[1])), 1, 1),
      GG = matrix(c(1,  0,  0,  0, 0,
                    0, -1, -1, -1, 0, 
                    0,  1,  0,  0, 0, 
                    0,  0,  1,  0, 0,
                    0,  0,  0,  0, 1),
                  5, 5, byrow = TRUE),
      W = matrix(c(exp(parm[2]),            0, 0, 0, 0,
                   0,            exp(parm[3]), 0, 0, 0,
                   0,                       0, 0, 0, 0,
                   0,                       0, 0, 0, 0,
                   0,                       0, 0, 0, 0),
                 5, 5, byrow = TRUE),
      m0 = c(0, 0, 0, 0, 0),
      C0 = matrix(c(1e+07,     0,     0,     0,     0,
                    0,     1e+07,     0,     0,     0,
                    0,         0, 1e+07,     0,     0, 
                    0,         0,     0, 1e+07,     0,
                    0,         0,     0,     0, 1e+07),
                  5, 5, byrow = TRUE),
      JFF = matrix(c(0, 0, 0, 0, 1), 1, 5),
      X = as.matrix(x))
}
dlmFitObj <- dlmMLE(y, parm = parm, build = buildObj)
(dlmModObj <- buildObj(dlmFitObj$par))
# $FF
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    1    1    0    0    1
# 
# $V
#             [,1]
# [1,] 2.27569e-05
# 
# $GG
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    1    0    0    0    0
# [2,]    0   -1   -1   -1    0
# [3,]    0    1    0    0    0
# [4,]    0    0    1    0    0
# [5,]    0    0    0    0    1
# 
# $W
#              [,1]         [,2] [,3] [,4] [,5]
# [1,] 1.783797e-05 0.000000e+00    0    0    0
# [2,] 0.000000e+00 4.310206e-07    0    0    0
# [3,] 0.000000e+00 0.000000e+00    0    0    0
# [4,] 0.000000e+00 0.000000e+00    0    0    0
# [5,] 0.000000e+00 0.000000e+00    0    0    0
# 
# $JFF
#      [,1] [,2] [,3] [,4] [,5]
# [1,]    0    0    0    0    1
# 
# $X
#      [,1]
# [1,] 0   
# [2,] 0   
# [3,] ... 
# 
# $m0
# [1] 0 0 0 0 0
# 
# $C0
#       [,1]  [,2]  [,3]  [,4]  [,5]
# [1,] 1e+07 0e+00 0e+00 0e+00 0e+00
# [2,] 0e+00 1e+07 0e+00 0e+00 0e+00
# [3,] 0e+00 0e+00 1e+07 0e+00 0e+00
# [4,] 0e+00 0e+00 0e+00 1e+07 0e+00
# [5,] 0e+00 0e+00 0e+00 0e+00 1e+07

p <- autoplot(dlmFilter(y, dlmModObj), fitted.colour = 'red', fitted.linetype = 'dashed')
autoplot(dlmSmooth(y, dlmModObj), colour = 'blue', p = p)

f:id:sinhrks:20150524155321p:plain

モデルの作り方

とはいえ、最初から上式のような形でモデルが定義できていることはないと思う。まずは簡単なモデルを dlm::dlmModPoly などで作って、それをベースに書き換えていくのが良さそう。作ったモデルに適当なパラメータを与えて初期化すれば、モデルの次元 / 中身を確認できる。

y <- ukdrivers
parm<- log(c(var(y), 0.001))
buildFun <- function(parm){
  dlmModPoly(order = 2, dV = exp(parm[1]), dW = c(exp(parm[2]), 0))
}
buildFun(parm)
# $FF
#      [,1] [,2]
# [1,]    1    0
# 
# $V
#            [,1]
# [1,] 0.02935256
# 
# $GG
#      [,1] [,2]
# [1,]    1    1
# [2,]    0    1
# 
# $W
#       [,1] [,2]
# [1,] 0.001    0
# [2,] 0.000    0
# 
# $m0
# [1] 0 0
# 
# $C0
#       [,1]  [,2]
# [1,] 1e+07 0e+00
# [2,] 0e+00 1e+07