R で 状態空間モデル: {dlm} で単変量モデルの状態空間表現を利用する
こちらの続きで、状態空間時系列分析入門の第8章の内容。
- 作者: J.J.F.コマンダー,S.J.クープマン,Jacques J.F. Commandeur,Sime Jan Koopman,和合肇
- 出版社/メーカー: シーエーピー出版
- 発売日: 2008/09
- メディア: 単行本
- 購入: 2人 クリック: 4回
- この商品を含むブログを見る
状態空間モデルでは、トレンド、季節効果、説明変数効果などさまざまな要因を考慮することができる。また これらの要因を単一の状態空間表現で統一的に扱える。時刻 における観測値 と状態 は以下の式で書ける (前回と同じ)。
この式であらわされたモデルを {dlm}
で利用する際にどうすればよいのかをまとめたい。
vignette にあるとおり、{dlm}
ではモデルを dlm::dlmModPoly
などのモデル記述関数の組み合わせで定義できる。これらの関数は便利だが、上の式との関係は外側からはよくわからない。
上式から {dlm}
のモデルを作成する際には dlm::dlm
を直接使うのがよさそうだ。dlm::dlm
の引数は上式にほぼ対応しているため、式表現されたモデルをそのまま実装できる。また、モデル記述関数も内部的には dlm::dlm
を使っているため、これら関数の理解にも役立つ。dlm::dlm
の引数と上式の対応は前回書いたとおり。
FF
: に対応。V
: に対応。GG
: に対応。W
: に対応。m0
: の事前分布の平均。C0
: の事前分布の分散。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. ローカルレベルモデル (確定的レベル)
それぞれのモデルについて 状態 として何が考慮されているかを書いておく。
- : は確定的レベル (経時変化しない)。
補足 それぞれのモデルの式表現を書こうとしたのだが、はてなマークダウン記法だと 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)
赤破線がカルマンフィルタ、青破線が平滑化した値。以降もプロットするが、確率的状態を使うと 実データにほぼフィットするため、グラフからでは差がよくわからない。
2. ローカル線形トレンドモデル (確率的レベルと確率的傾き)
- : は確率的レベル (経時変化する)、 は確率的傾き(経時変化する)。
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)
3. 季節要素のあるローカルレベルモデル (確率的レベルと確率的季節)
- : は確率的レベル、 は 四半期単位の確率的季節のダミー項。
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)
4. 説明変数のあるローカルレベルモデル (確率的レベル)
- : は確率的レベル、 は説明変数の係数。
説明変数は 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)
5. 干渉変数のあるローカルレベルモデル (確率的レベル)
- : は確率的レベル、 は干渉変数の係数。
変数名は異なるが、式は説明変数のあるローカルレベルモデルと同一。
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)
6. 季節要素と干渉変数のあるローカルレベルモデル (確率的レベルと確率的季節)
- : は確率的レベル、 は確率的季節、 は干渉変数の係数。
補足 テキスト記載の数式は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)
モデルの作り方
とはいえ、最初から上式のような形でモデルが定義できていることはないと思う。まずは簡単なモデルを 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