仕事の都合で仕方なく状態空間モデルについて勉強していたのだけれど(なぜ私がこんな目に)、仕事で使うためには自分で計算できるようにならなければならない。
参考にしているCommandeur & Koopman 「状態空間時系列分析入門」（以下「CK本」）の著者らは、すべての事例についてデータとプログラムを公開している。ありがたいことであります。しかし、ssfpackという耳慣れないソフトを使わなければならない。わざわざ新しいソフトの使い方を覚えるのは大変に面倒だ。できれば普段使っているソフトで済ませたい。
というわけで、勉強かたがた、CK本に出てくる計算例を片っ端から R で再現してみた。汗と涙の甲斐あって、すべての章についていちおう再現できたので、ここに載せておくことにする。
もくじ:
Rのプログラムはこちらに置いておきますが、全ての計算例を再現しているので、それなりに長い。雰囲気をご覧頂くために一例だけご紹介します。
CK本の4.2節。この本を通して用いられている、英国ドライバー死傷者数(対数)の時系列についての分析である。
この時系列に、いわゆるローカル・レベル・モデルをあてはめる。状態変数はレベル(死傷者数の潜在的な大きさ)と季節要素(12ヶ月)、ともに確率的に変動させる。
時点$t$の観察値を$y_t$, 不規則要素を$\epsilon_t$, レベルを$\mu_t$, レベル撹乱要素を$\xi_t$, 季節要素を$\gamma_t$, 季節撹乱要素を$\omega_t$として、
$y_t = \mu_t + \gamma_t + \epsilon_t$
$\mu_{t+1} = \mu_t + \xi_t$
$t=11, \ldots, n$ に対して $ \gamma_{t+1} = - (\gamma_t + \gamma_{t-1} + \cdots + \gamma_{t-10} ) + \omega_t$
$\epsilon_t \sim NID(0, \sigma^2_\epsilon)$
$\xi_t \sim NID(0, \sigma^2_\xi)$
$\omega_t \sim NID(0, \sigma^2_\omega)$
4.2節は、Rプログラムのexercise.Chapter4.2()という関数で再現している。
まずライブラリを読み込み、定数を定義。gCKLogLikはCK本に掲載されている最大対数尤度で、あとで比較表を作るために使う。
exercise.Chapter4.2 <- function(){
library(dlm)
library(KFAS)
# constants - - - - - -
sFILENAME <- "UKdriversKSI.txt"
gCKLogLik <- 0.9369063
# - - - - - - - - - - -
データを読み込む。今回の方針として、用いる時系列データは要不要を問わず必ずtsクラスに放り込むことにしている。tsDataが英国ドライバー死傷者数の対数の時系列。nが時点数(=192)。
# read data
dfData <- read.table(sFILENAME, skip=1)
tsData <- ts(log(dfData[,1]), frequency=12, start=c(1969, 1))
n <- length(tsData)
nStateVarとnHyperParamは、それぞれ状態変数の数とハイパーパラメータの数。分析結果のまとめの際に用いる。季節が11個のダミー変数で表現されるので、状態変数はレベルと季節で計12個。ハイパー・パラメータは$\sigma^2_\epsilon, \sigma^2_\xi, \sigma^2_\omega$の3個である。
agInitはハイパーパラメータの最尤推定の際に与える初期値。適当に決めている。
ここまでが前置きです。
# model specs
nStateVar <- 12
nHyperParam <- 3
# init value for variances (epsilon, xi, omega)
agInit <- c(var(tsData), 0.001, 0.001)
cat("\n")
dlmパッケージによるモデリング。
まずモデルを表現する関数funModel()を定義する。
m1はローカル・レベル・モデルを表している。dVは不規則要素の分散 $\sigma^2_\epsilon$ 、dWはレベル撹乱要素の分散(正確にいうと、レベル撹乱要素の分散行列の対角要素) $\sigma^2_\xi$である。どちらも未知パラメータとして指定している。
m2はこれに付け加える季節変動モデル。dWは撹乱要素の分散行列の対角要素で、季節は11個のダミー変数で表現するので長さ11のベクトルになる。撹乱要素をいれるのは最初のダミー変数だけでよいので、dWは結局$(\sigma^2_\omega, 0,0,0,0,0,0,0,0,0,0)$となる。季節要素も確率的に変動すると考え、撹乱要素の分散を未知パラメータとして指定する。なお、不規則要素の分散はさっき指定したので、ここのdVは0でよい。
ここがdlmでのもっとも大事な部分なので、誤りがないか念入りに確認する。funModel()の挙動を調べたコードが残ってますね、ははは。
cat("estimated by dlm package:\n")
funModel <- function(parm){
m1 <- dlmModPoly(order = 1, dV=exp(parm[1]), dW=exp(parm[2]))
m2 <- dlmModSeas(12, dV=0, dW=c(exp(parm[3]), rep(0,10)))
return( m1 + m2 )
}
## print(funModel(log(c(123,456,789))))
パラメータを最尤推定。うまくいったら、この推定されたパラメータを改めて funModel()に与え、モデルを完成する(oFitted.DLM)。
# Fitting
oMLE.DLM <- dlmMLE(
tsData,
parm = log(agInit),
build = funModel
)
stopifnot(oMLE.DLM$convergence == 0)
oFitted.DLM <- funModel(oMLE.DLM$par)
このモデルを使ってtsDataをフィルタ化する(oFiltered.DLM)。あとで使うために、予測誤差と標準化予測誤差を抽出する(agStdPredErr.DLM, lPredErr.DLM)。
また、モデルを使ってtsDataを平滑化する(oSmoothed.DLM)。あとで使うために、平滑化状態系列を抽出する(tsSmostate.DLM)。dlmSmoothの返し値のsは多変量時系列のtsオブジェクトになっているので便利なのだが、観察時点の1期前から値が入っている点に注意。dlmパッケージに入っているdropFirst()で最初の時点を削除している。
# Filtering
oFiltered.DLM <- dlmFilter(tsData, oFitted.DLM)
agStdPredErr.DLM <- residuals(oFiltered.DLM, type = c("standardized"), sd=FALSE)
lPredErr.DLM <- residuals(oFiltered.DLM, type = c("raw"), sd=TRUE)
# Smoothing
oSmoothed.DLM <- dlmSmooth(tsData, oFitted.DLM)
tsSmoState.DLM <- dropFirst(oSmoothed.DLM$s)
ここからは結果の表示。最大対数尤度を2種類求めておく(詳細は後述)。また、推定したハイパー・パラメータと初期状態を表示する。
# LL
gLLbyFun.DLM <- -oMLE.DLM$value - 0.5*n*log(2*pi)
gLLbyErr.DLM <- sub.LogLik(lPredErr.DLM$res, lPredErr.DLM$sd^2, nStateVar)
# report
agW <- diag(oFitted.DLM$W)
cat(sprintf("hat(sigma^2_epsilon) = %f\n", drop(oFitted.DLM$V)))
cat(sprintf("hat(sigma^2_xi) = %e\n", agW[1]))
cat(sprintf("hat(sigma^2_omega) = %e\n", agW[2]))
cat(sprintf("hat(mu_1) = %f\n", oSmoothed.DLM$s[2,1]))
cat(sprintf("hat(gamma_1) = %f\n", oSmoothed.DLM$s[2,2]))
cat("\n")
なお、この部分の出力は下記の通り。順に、不規則要素の分散$\sigma^2_\epsilon$、レベル撹乱要素の分散$\sigma^2_\xi$、季節撹乱要素の分散$\sigma^2_\omega$、時点1のレベル$\mu_1$と季節要素$\gamma_1$、の推定値である。
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 9.458618e-04
hat(sigma^2_omega) = 4.368203e-10
hat(mu_1) = 7.411847
hat(gamma_1) = 0.017273
おつぎはKFASパッケージによるモデリング。
モデルを表現するオブジェクトoModel.KFASを定義する。KFASでは、未知パラメータであることをNAで表現する。
まずはモデル式で表現する。右辺のSSMtrend()はローカル・レベル・モデルを表している。Qが不規則要素の分散行列のリストで、ここではリストの要素はひとつ(状態変数がレベルだけだから)。その行列のサイズは1x1、中身は$\sigma^2_\xi$である。SSMseasonal()はこれに付け加える季節要素を表している。Qをサイズ1x1の行列にしていることに注意。こうするとQは、最初のダミー変数にのみ効く攪乱要素の分散、すなわち $\sigma^2_\omega$となる。この2つを未知パラメータに指定。
H行列は不規則要素の分散行列、つまり$\sigma^2_\epsilon$である。これも未知パラメータに指定。
KFASにとってこの部分が肝になるので、念入りに調べる必要がある。oModel.KFASをprint()するといろいろ教えてくれるので、メモをとっている。状態変数は(レベル, 季節ダミー1, 季節ダミー2,...)の順に並んでいる。レベル撹乱要素の分散行列はサイズ2x2になっている。
cat("estimated by KFAS package:\n")
oModel.KFAS <- SSModel(
tsData ~ SSMtrend(degree=1, Q=list(matrix(NA))) +
SSMseasonal(12, sea.type="dummy", Q=matrix(NA), n=n),
H=matrix(NA)
)
## print(oModel.KFAS) -> states are (level, seasondummy1, seasondummy2, ...)
## print(oModel.KFAS$Q) -> Q is (2,2) matrix where diag is (NA,NA)
パラメータを最尤推定。うまくいったら、この推定されたモデルをKFS()に与える。KFS()はデフォルトでフィルタ化と平滑化の両方を行ってくれる。あとで使うために予測誤差を抽出する(agStdPredErr.KFAS)。
# Fitting
oFitted.KFAS <- fitSSM(
oModel.KFAS,
inits = log(agInit)
)
stopifnot(oFitted.KFAS$optim.out$convergence == 0)
# Filtering & Smoothing
oEstimated.KFAS <- KFS(oFitted.KFAS$model)
agStdPredErr.KFAS <- rstandard(oEstimated.KFAS, type="recursive")
ここからは結果の表示。最大対数尤度を2種類求めておく。また、推定したハイパー・パラメータと初期状態を表示する。
# LL
gLLbyFun.KFAS <- oEstimated.KFAS$logLik
gLLbyErr.KFAS <- sub.LogLik(oEstimated.KFAS$v[,1], t(oEstimated.KFAS$F), nStateVar)
# report
agQ <- diag(oFitted.KFAS$model$Q[,,1])
cat(sprintf("hat(sigma^2_epsilon) = %f\n", drop(oFitted.KFAS$model$H)))
cat(sprintf("hat(sigma^2_xi) = %e\n", agQ[1]))
cat(sprintf("hat(sigma^2_omega) = %e\n", agQ[2]))
cat(sprintf("hat(mu_1) = %f\n", coef(oEstimated.KFAS)[1,1]))
cat(sprintf("hat(gamma_1) = %f\n", coef(oEstimated.KFAS)[1,2]))
cat("\n")
この部分の出力は下記の通り。dlmとほぼ同じ推定結果が得られている。
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 9.455296e-04
hat(sigma^2_omega) = 5.481093e-11
hat(mu_1) = 7.411848
hat(gamma_1) = 0.017272
dlmとKFASから得られた最大対数尤度と、CK本に載っている最大対数尤度を比較する。比較表を表示する関数sub.ShowLogLik()は別途に定義してある。
cat("max loglikelihood and AIC:\n")
sub.ShowLogLik (
n, nStateVar, nHyperParam,
gLLbyFun.DLM, gLLbyErr.DLM, gLLbyFun.KFAS, gLLbyErr.KFAS, gCKLogLik*n
)
cat("\n")
この部分の出力は下記の通り。一行目はCK本に基づきdlmの予測誤差系列とその分散の系列から求めた値。対数尤度は180.1930、時点数で割ると0.939, AICは-330.39, 時点数で割ると-1.721。二行目はdlmが出力した対数尤度で、なんと99.193も小さい(小さくなる理由は後述)。三行目と四行目はKFASでの結果。最後の行ではCK本に載っている値(=ssfpackで求めた値)を示している。このように、最大対数尤度はソフトによってちがうんですよ、奥さん...
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 180.1930 ( 0.939) -330.39 (-1.721)
(by output) 80.9995 [ -99.193] ( 0.422) -132.00 (-0.687)
KFAS (by pred.err) 180.1930 [ +0.000] ( 0.939) -330.39 (-1.721)
(by output) 177.7081 [ -2.485] ( 0.926) -325.42 (-1.695)
The CK Book 179.8860 [ -0.307] ( 0.937) -329.77 (-1.718)
-----------------------------------------------------------------
残差診断の結果を表示する。関数sub.ShowDiagnostics()は別途に定義してある。
cat("diagnostics of standardized prediction error (estimated by dlm):\n")
sub.ShowDiagnostics (
agStdPredErr.DLM, nStateVar, nHyperParam, 15, c(1,12)
)
cat("\n")
この部分の出力は下記の通り。CK本とほぼ同じ結果が出力されるが、例題によっては均一分散性の臨界値がCK本と異なることがある(その理由は後述)。
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 14.369 22.36 +
r ( 1) 0.040 +-0.14 +
r (12) 0.033 +-0.14 +
homoscedasticity H (60) 1.093 1.67 +
normality N 5.156 5.99 +
------------------------------------------------------
最後に、図の作成。それぞれの関数は別途定義してある(どれもごく単純な関数である)。
# plot
sub.plotFigure_twoseries(
tsData,
tsSmoState.DLM[,1],
"Figure 4.6",
c("log UK drivers KSI", "stochastic level")
)
sub.plotFigure_twoseries(
NULL,
tsSmoState.DLM[,2],
"Figure 4.7",
"stochastic seasonal"
)
sub.plotFigure_twoseries(
NULL,
ts(window(tsSmoState.DLM[,2], start=c(1969,1), end=c(1969,12)), frequency = 1),
"Figure 4.8",
"seasonal 1969"
)
sub.plotFigure_residual(
tsData - (tsSmoState.DLM[,1]+tsSmoState.DLM[,2]),
"Figure 4.9",
"irregular"
)
}
最初に出力される図4.6はこんな感じ。赤線が時系列の背後にあると推定されるレベルである。意味ありげな減少が2回生じているのがわかる。
主に参照するのは訳書(2008, 第1刷)のほう。訳者の先生に感謝。ときどき原著(2007)もチェックした。
なお、この本の計算例をRで再現しようという酔狂な人は私だけではないようで... すでに大変充実した解説を公開されている方がいる。私も何箇所か参考にさせていただきました。感謝しております。私とは再現の範囲が違うし、ところどころ間違えているところもあると思うので(もちろん私の方も)、見比べていただければと思います。
Rのプログラムはこちら。実行にあたっては、著者らが配布しているデータとプログラムに含まれている、以下のデータファイルを作業ディレクトリに置く必要がある。
話をわかりやすくするために、すべての例題について同じ変数名・関数名をつかうように心がけた。たとえばdlmのdlmFilter()の返し値は、どの例題でも変数 oFiltered.DLM に格納している。以下、このように私が命名した変数名なり関数名を、ブラケットをつけて[oFiltered.DLM]というように表記する。
m1 <- dlmModPoly(order=1, dW=123)
m2 <- dlmModSeas(4, dW=c(456, 0, 0))
library(KFAS)
X <- matrix(rnorm(100), ncol=2)
m <- SSMregression (~ 1 + X, Q = matrix(123))
英国ドライバー死傷者数を使った例示。本文中の計算例・チャートを再現したアウトプットは以下のとおり。
> exercise.Chapter1 ()
Linear regression:
Call:
lm(formula = tsData ~ time(tsData))
Coefficients:
(Intercept) time(tsData)
7.545843 -0.001448
plotting Figure 1.1 ...
plotting Figure 1.2 ...
plotting Figure 1.3 ...
plotting Figure 1.4 ...
plotting Figure 1.5 ...
英国ドライバー死傷者数。状態変数はレベル(確定的)。実質的には平均と分散を推定するだけである。
> exercise.Chapter2.1 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.029353
hat(mu_1) = 7.406108
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.029353
hat(mu_1) = 7.406108
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 62.3949 ( 0.325) -120.79 (-0.629)
(by output) 54.3359 [ -8.059] ( 0.283) -104.67 (-0.545)
KFAS (by pred.err) 62.3949 [ -0.000] ( 0.325) -120.79 (-0.629)
(by output) 62.3949 [ -0.000] ( 0.325) -120.79 (-0.629)
The CK Book 63.3139 [ +0.919] ( 0.330) -122.63 (-0.639)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 415.212 25.00 -
r ( 1) 0.699 +-0.14 -
r (12) 0.677 +-0.14 -
homoscedasticity H (64) 2.058 1.64 -
normality N 0.733 5.99 +
------------------------------------------------------
plotting Figure 2.1 ...
plotting Figure 2.2 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)。
> exercise.Chapter2.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.002221
hat(sigma^2_xi) = 0.011866
hat(mu_1) = 7.414954
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.002221
hat(sigma^2_xi) = 0.011867
hat(mu_1) = 7.414958
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 122.9587 ( 0.640) -239.92 (-1.250)
(by output) 114.8996 [ -8.059] ( 0.598) -223.80 (-1.166)
KFAS (by pred.err) 122.9587 [ -0.000] ( 0.640) -239.92 (-1.250)
(by output) 122.9587 [ -0.000] ( 0.640) -239.92 (-1.250)
The CK Book 123.8776 [ +0.919] ( 0.645) -241.76 (-1.259)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 105.382 23.68 -
r ( 1) 0.009 +-0.14 +
r (12) 0.537 +-0.14 -
homoscedasticity H (64) 1.064 1.64 +
normality N 13.242 5.99 -
------------------------------------------------------
plotting Figure 2.3 ...
plotting Figure 2.4 ...
ノルウェイ道路交通事故数。状態変数はレベル(確率的)。ローカル･レベル･モデルでうまく記述できる例である。
> exercise.Chapter2.3 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003268
hat(sigma^2_xi) = 0.004703
hat(mu_1) = 6.304804
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003268
hat(sigma^2_xi) = 0.004703
hat(mu_1) = 6.304802
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 27.8744 ( 0.820) -49.75 (-1.463)
(by output) 19.8153 [ -8.059] ( 0.583) -33.63 (-0.989)
KFAS (by pred.err) 27.8744 [ -0.000] ( 0.820) -49.75 (-1.463)
(by output) 27.8744 [ -0.000] ( 0.820) -49.75 (-1.463)
The CK Book 28.7933 [ +0.919] ( 0.847) -51.59 (-1.517)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm)
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (10) 6.228 16.92 +
r ( 1) -0.127 +-0.34 +
r ( 4) -0.105 +-0.34 +
homoscedasticity H (11) 1.746 3.47 +
normality N 1.191 5.99 +
------------------------------------------------------
plotting Figure 2.5 ...
plotting Figure 2.6 ...
plotting Figure 8.5 ...
plotting Figure 8.7a ...
plotting Figure 8.7b ...
plotting Figure 8.13 ...
8.4節の図8.5, 図8.7もここで出力する。フィルタ化レベル系列と平滑化レベル系列の比較、予測誤差系列とその分散の系列。
8.6節の図8.13もここで出力する。フィルタ化レベルならびに観測期間後予測とその信頼区間。
英国ドライバー死傷者数。状態変数はレベル(確定的)と傾き(確定的)。実質的には単回帰である。
> exercise.Chapter3.1 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.022998
hat(mu_1) = 7.544395
hat(nu_1) = -0.001448
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.022998
hat(mu_1) = 7.544395
hat(nu_1) = -0.001448
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 77.6641 ( 0.405) -149.33 (-0.778)
(by output) 61.5460 [ -16.118] ( 0.321) -117.09 (-0.610)
KFAS (by pred.err) 77.6641 [ -0.000] ( 0.405) -149.33 (-0.778)
(by output) 77.6641 [ -0.000] ( 0.405) -149.33 (-0.778)
The CK Book 79.5020 [ +1.838] ( 0.414) -153.00 (-0.797)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 305.682 25.00 -
r ( 1) 0.610 +-0.14 -
r (12) 0.631 +-0.14 -
homoscedasticity H (63) 1.360 1.65 +
normality N 1.790 5.99 +
------------------------------------------------------
英国ドライバー死傷者数。状態変数はレベル(確率的)と傾き(確定的)。
> exercise.Chapter3.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.002119
hat(sigma^2_xi) = 0.012127
hat(sigma^2_zeta) = 0.000000
hat(mu_1) = 7.415732
hat(nu_1) = 0.000289
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.002117
hat(sigma^2_xi) = 0.012130
hat(sigma^2_zeta) = 0.000000
hat(mu_1) = 7.415741
hat(nu_1) = 0.000289
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 118.1224 ( 0.615) -226.24 (-1.178)
(by output) 102.0043 [ -16.118] ( 0.531) -194.01 (-1.010)
KFAS (by pred.err) 118.1225 [ +0.000] ( 0.615) -226.24 (-1.178)
(by output) 118.1225 [ +0.000] ( 0.615) -226.24 (-1.178)
The CK Book 119.9604 [ +1.838] ( 0.625) -229.92 (-1.198)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 100.615 22.36 -
r ( 1) 0.005 +-0.14 +
r (12) 0.532 +-0.14 -
homoscedasticity H (63) 1.058 1.65 +
normality N 14.946 5.99 -
------------------------------------------------------
plotting Figure 3.1 ...
plotting Figure 3.2 ...
plotting Figure 3.3 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)と傾き(確定的)。
収束すると、対数尤度関数の値は0.6247935である。不規則要素[observation disturbances]の分散の最尤推定値は $\hat\sigma^2_\epsilon = 0.00211869$、レベル撹乱要素の分散の最尤推定値は $\hat\sigma^2=0.0121271$ である。系列の最初の時点におけるレベルと傾きの最尤推定値はそれぞれ$\hat\mu_2 = 7.4157$と$\hat\nu_1 = 0.00028897$である。
> exercise.Chapter3.3 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.002118
hat(sigma^2_xi) = 0.012128
hat(mu_1) = 7.415736
hat(nu_1) = 0.000289
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.002119
hat(sigma^2_xi) = 0.012127
hat(mu_1) = 7.415731
hat(nu_1) = 0.000289
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 118.1225 ( 0.615) -228.24 (-1.189)
(by output) 102.0044 [ -16.118] ( 0.531) -196.01 (-1.021)
KFAS (by pred.err) 118.1225 [ -0.000] ( 0.615) -228.24 (-1.189)
(by output) 118.1225 [ -0.000] ( 0.615) -228.24 (-1.189)
The CK Book 119.9604 [ +1.838] ( 0.625) -231.92 (-1.208)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 100.609 23.68 -
r ( 1) 0.005 +-0.14 +
r (12) 0.532 +-0.14 -
homoscedasticity H (63) 1.058 1.65 +
normality N 14.946 5.99 -
------------------------------------------------------
plotting Figure 3.4 ...
フィンランドの道路交通事故数。状態変数はレベルと傾き。ローカル線形トレンドモデルがうまくいく場合である。まずは最初の段落の、レベルも傾きも確率的な場合。
> exercise.Chapter3.4a ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003201
hat(sigma^2_xi) = 3.931371e-09
hat(sigma^2_zeta) = 0.001533
hat(mu_1) = 7.013342
hat(nu_1) = 0.006849
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003201
hat(sigma^2_xi) = 9.179120e-10
hat(sigma^2_zeta) = 0.001533
hat(mu_1) = 7.013343
hat(nu_1) = 0.006848
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 24.9023 ( 0.732) -39.80 (-1.171)
(by output) 8.7842 [ -16.118] ( 0.258) -7.57 (-0.223)
KFAS (by pred.err) 24.9023 [ +0.000] ( 0.732) -39.80 (-1.171)
(by output) 24.9023 [ +0.000] ( 0.732) -39.80 (-1.171)
The CK Book 26.7401 [ +1.838] ( 0.786) -43.48 (-1.279)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 9.873 22.36 +
r ( 1) -0.028 +-0.34 +
r (12) -0.117 +-0.34 +
homoscedasticity 1/H(11) 1.348 3.47 +
normality N 0.644 5.99 +
------------------------------------------------------
レベルを確定的にする。
収束すると、対数尤度関数の値は0.7864746である。不規則要素と傾きの分散[observation and slope disturbances]の最尤推定値は、それぞれ $\hat\sigma^2_\epsilon = 0.00320083$、$\hat\sigma^2_\zeta = 0.00153314$ である[2個目の$\sigma$の添え字は$\xi$ではなく$\zeta$]。系列の最初の時点におけるレベルと傾きの最尤推定値はそれぞれ$\hat\mu_2 = 7.0133$と$\hat\nu_1 = 0.0068482$である。
> exercise.Chapter3.4b ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003201
hat(sigma^2_zeta) = 0.001533
hat(mu_1) = 7.013343
hat(nu_1) = 0.006848
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003201
hat(sigma^2_zeta) = 0.001533
hat(mu_1) = 7.013344
hat(nu_1) = 0.006847
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 24.9023 ( 0.732) -41.80 (-1.230)
(by output) 8.7842 [ -16.118] ( 0.258) -9.57 (-0.281)
KFAS (by pred.err) 24.9023 [ +0.000] ( 0.732) -41.80 (-1.230)
(by output) 24.9023 [ +0.000] ( 0.732) -41.80 (-1.230)
The CK Book 26.7401 [ +1.838] ( 0.786) -45.48 (-1.338)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (10) 7.044 16.92 +
r ( 1) -0.028 +-0.34 +
r ( 4) -0.094 +-0.34 +
homoscedasticity 1/H(11) 1.348 3.47 +
normality N 0.644 5.99 +
------------------------------------------------------
plotting Figure 3.5a ...
plotting Figure 3.5b ...
plotting Figure 3.6 ...
plotting Figure 8.14 ...
ここで8.6節の図8.14も出力する。観測期間後予測とその信頼区間。
英国ドライバー死傷者数。状態変数はレベル(確定的)と季節(確定的)。
> exercise.Chapter4.1 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.017589
hat(mu_1) = 7.406107
hat(gamma_1}) = 0.022155
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.017589
hat(mu_1) = 7.406108
hat(gamma_1) = 0.022155
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 80.5740 ( 0.420) -135.15 (-0.704)
(by output) -18.6195 [ -99.193] (-0.097) 63.24 ( 0.329)
KFAS (by pred.err) 80.5740 [ -0.000] ( 0.420) -135.15 (-0.704)
(by output) 78.0891 [ -2.485] ( 0.407) -130.18 (-0.678)
The CK Book 80.1576 [ -0.416] ( 0.417) -134.32 (-0.700)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 751.575 25.00 -
r ( 1) 0.724 +-0.14 -
r (12) 0.431 +-0.14 -
homoscedasticity H (60) 3.400 1.67 -
normality N 1.971 5.99 +
------------------------------------------------------
plotting Figure 4.2 ...
plotting Figure 4.3 ...
plotting Figure 4.4 ...
plotting Figure 4.5 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)と季節(確率的)。
> exercise.Chapter4.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 9.458618e-04
hat(sigma^2_omega) = 4.368203e-10
hat(mu_1) = 7.411847
hat(gamma_1) = 0.017273
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 9.455296e-04
hat(sigma^2_omega) = 5.481093e-11
hat(mu_1) = 7.411848
hat(gamma_1) = 0.017272
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 180.1930 ( 0.939) -330.39 (-1.721)
(by output) 80.9995 [ -99.193] ( 0.422) -132.00 (-0.687)
KFAS (by pred.err) 180.1930 [ +0.000] ( 0.939) -330.39 (-1.721)
(by output) 177.7081 [ -2.485] ( 0.926) -325.42 (-1.695)
The CK Book 179.8860 [ -0.307] ( 0.937) -329.77 (-1.718)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 14.369 22.36 +
r ( 1) 0.040 +-0.14 +
r (12) 0.033 +-0.14 +
homoscedasticity H (60) 1.093 1.67 +
normality N 5.156 5.99 +
------------------------------------------------------
plotting Figure 4.6 ...
plotting Figure 4.7 ...
plotting Figure 4.8 ...
plotting Figure 4.9 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)と季節(確定的)。
> exercise.Chapter4.3 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 0.000946
hat(mu_1) = 7.411848
hat(gamma_1) = 0.017272
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003514
hat(sigma^2_xi) = 0.000946
hat(mu_1) = 7.411848
hat(gamma_1) = 0.017272
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 180.1930 ( 0.939) -332.39 (-1.731)
(by output) 80.9995 [ -99.193] ( 0.422) -134.00 (-0.698)
KFAS (by pred.err) 180.1930 [ -0.000] ( 0.939) -332.39 (-1.731)
(by output) 177.7081 [ -2.485] ( 0.926) -327.42 (-1.705)
The CK Book 179.7765 [ -0.416] ( 0.936) -331.55 (-1.727)
-----------------------------------------------------------------
plotting Figure 8.1 ...
plotting Figure 8.2 ...
plotting Figure 8.3 ...
plotting Figure 8.4 (by dlm) ...
plotting Figure 8.4 (by KFAS) ...
plotting Figure 8.11a ...
plotting Figure 8.11b ...
ここで8.3節の図8.1, 図8.2, 図8.3, 図8.4も出力する。平滑化レベルの推定誤差分散の系列、平滑化レベル系列とその信頼区間、平滑化季節要素系列とその信頼区間、平滑化シグナル系列とその信頼区間。
ところで、平滑化シグナルの推定誤差分散の求め方がわからない...
8.5節の図8.11節も出力する。補助残差(標準化した平滑化レベル撹乱要素系列と標準化した平滑化不規則要素系列)。
英国インフレーション。状態変数はレベル(確率的)と季節(確率的)。4.2節と同じタイプのモデルである。
> exercise.Chapter4.4 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 3.371270e-05
hat(sigma^2_xi) = 2.124158e-05
hat(sigma^2_gamma) = 4.345176e-07
hat(mu_1) = 0.006043
hat(gamma_1) = 0.001931
hat(mu_208) = 0.002049
estimated by KFAS package:
hat(sigma^2_epsilon) = 3.372868e-05
hat(sigma^2_xi) = 2.123500e-05
hat(sigma^2_gamma) = 4.347504e-07
hat(mu_1) = 0.006045
hat(gamma_1) = 0.001931
hat(gamma_208) = 0.002050
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 663.5979 ( 3.190) -1313.20 (-6.313)
(by output) 629.9754 [ -33.622] ( 3.029) -1245.95 (-5.990)
KFAS (by pred.err) 663.5979 [ -0.000] ( 3.190) -1313.20 (-6.313)
(by output) 662.2116 [ -1.386] ( 3.184) -1310.42 (-6.300)
The CK Book 665.2805 [ +1.683] ( 3.198) -1316.56 (-6.330)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (10) 7.561 15.51 +
r ( 1) 0.049 +-0.14 +
r ( 4) -0.061 +-0.14 +
homoscedasticity 1/H(68) 2.737 1.62 -
normality N 172.198 5.99 -
------------------------------------------------------
plotting Figure 4.10a ...
plotting Figure 4.10b ...
plotting Figure 4.10c ...
英国ドライバー死傷者数。状態変数はレベル(確定的)と回帰係数(確定的)。まずは時点番号を説明変数にする。
> exercise.Chapter5.1a ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.022998
hat(mu_1) = 7.545843
hat(beta_1) = -0.001448
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.022998
hat(mu_1) = 7.545843
hat(beta_1) = -0.001448
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 77.6641 ( 0.405) -149.33 (-0.778)
(by output) 61.5460 [ -16.118] ( 0.321) -117.09 (-0.610)
KFAS (by pred.err) 77.6641 [ -0.000] ( 0.405) -149.33 (-0.778)
(by output) 77.6641 [ -0.000] ( 0.405) -149.33 (-0.778)
The CK Book 79.5020 [ +1.838] ( 0.414) -153.00 (-0.797)
-----------------------------------------------------------------
説明変数を石油価格に変更。
> exercise.Chapter5.1b ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.023014
hat(mu_1) = 5.878731
hat(beta_1) = -0.671664
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.023014
hat(mu_1) = 5.878731
hat(beta_1) = -0.671664
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 78.6132 ( 0.409) -151.23 (-0.788)
(by output) 67.6223 [ -10.991] ( 0.352) -129.24 (-0.673)
KFAS (by pred.err) 78.6126 [ -0.001] ( 0.409) -151.23 (-0.788)
(by output) 83.7404 [ +5.127] ( 0.436) -161.48 (-0.841)
The CK Book 85.5783 [ +6.965] ( 0.446) -165.16 (-0.860)
-----------------------------------------------------------------
plotting Figure 5.1 ...
plotting Figure 5.2 ...
plotting Figure 5.3 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)と回帰係数(確定的)。
> exercise.Chapter5.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.002348
hat(sigma^2_xi) = 0.011667
hat(mu_1) = 6.820359
hat(beta_1) = -0.261050
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.002348
hat(sigma^2_xi) = 0.011667
hat(mu_1) = 6.820355
hat(beta_1) = -0.261052
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 116.9968 ( 0.609) -225.99 (-1.177)
(by output) 106.0062 [ -10.991] ( 0.552) -204.01 (-1.063)
KFAS (by pred.err) 116.9965 [ -0.000] ( 0.609) -225.99 (-1.177)
(by output) 122.1243 [ +5.127] ( 0.636) -236.25 (-1.230)
The CK Book 123.9621 [ +6.965] ( 0.646) -239.92 (-1.250)
-----------------------------------------------------------------
plotting Figure 5.4 ...
plotting Figure 5.5 ...
英国ドライバー死傷者数。時点160まで0, 時点170から1になる干渉変数を投入する。状態変数はレベル(確定的)と回帰係数(確定的)。
> exercise.Chapter6.1 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.022243
hat(mu_1) = 7.437386
hat(lambda_1) = -0.261108
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.022243
hat(mu_1) = 7.437386
hat(lambda_1) = -0.261108
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 76.5029 ( 0.398) -147.01 (-0.766)
(by output) 69.8587 [ -6.644] ( 0.364) -133.72 (-0.696)
KFAS (by pred.err) 81.4368 [ +4.934] ( 0.424) -156.87 (-0.817)
(by output) 85.9768 [ +9.474] ( 0.448) -165.95 (-0.864)
The CK Book 87.8147 [ +11.312] ( 0.457) -169.63 (-0.883)
-----------------------------------------------------------------
plotting Figure 6.1 ...
plotting Figure 6.2 ...
plotting Figure 6.3 ...
making additional chart ...
ここまでの例題では、最大対数尤度はソフトによってそれぞれ異なっていたものの、予測誤差系列とその分散の系列は dlm でも KFAS でも一致した。ところが、6.1節以降の干渉変数を含むモデルでは、予測誤差の分散の系列がずれる。従って予測誤差系列とその分散の系列から再計算した最大対数尤度もずれる。
6.1節の場合の様子を下図に示す。分散の値が大きく異なるのは次の2時点である。ひとつは最初の時点で、これは最大対数尤度の再計算には影響しない。もうひとつは1983年2月、すなわち干渉変数がはじめて1になる時点で、KFASでは前月とほぼ同じ値 (0.0223) なのに対し、dlmでは10,000,000と出力される(図では欠損として扱っている)。
詳しいことはわからないが、おそらくカルマン・フィルタの初期化の仕組みがちょっとちがうのであろう。
英国ドライバー死傷者数。状態変数はレベル(確率的)と回帰係数(確定的)。
> exercise.Chapter6.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.002693
hat(sigma^2_xi) = 0.010412
hat(mu_1) = 7.410664
hat(lambda_1) = -0.378495
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.002692
hat(sigma^2_xi) = 0.010411
hat(mu_1) = 7.410664
hat(lambda_1) = -0.378495
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 115.7397 ( 0.603) -223.48 (-1.164)
(by output) 109.3564 [ -6.383] ( 0.570) -210.71 (-1.097)
KFAS (by pred.err) 120.6065 [ +4.867] ( 0.628) -233.21 (-1.215)
(by output) 125.4745 [ +9.735] ( 0.654) -242.95 (-1.265)
The CK Book 127.3123 [ +11.573] ( 0.663) -246.62 (-1.285)
-----------------------------------------------------------------
plotting Figure 6.4 ...
plotting Figure 6.5 ...
いよいよ複雑になってまいりました。英国ドライバー死傷者数。説明変数(石油価格)と干渉変数(英国シートベルト法)をいれる。状態変数はレベル(確定的)、季節(確定的)、石油価格の回帰係数(確定的)、シートベルト法の回帰係数(確定的)。
> exercise.Chapter7.1 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.007402
hat(mu_1) = 6.401570
hat(gamma_1) = 0.006701
hat(beta_1) = -0.452131
hat(lambda_1) = -0.197140
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.007402
hat(mu_1) = 6.401571
hat(gamma_1) = 0.006700
hat(beta_1) = -0.452130
hat(lambda_1) = -0.197139
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 138.9588 ( 0.724) -247.92 (-1.291)
(by output) 37.3235 [-101.635] ( 0.194) -44.65 (-0.233)
KFAS (by pred.err) 142.8375 [ +3.879] ( 0.744) -255.67 (-1.332)
(by output) 150.1502 [ +11.191] ( 0.782) -270.30 (-1.408)
The CK Book 154.0565 [ +15.098] ( 0.802) -278.11 (-1.449)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 147.020 25.00 -
r ( 1) 0.426 +-0.14 -
r (12) 0.198 +-0.14 -
homoscedasticity 1/H(59) 1.110 1.67 +
normality N 0.560 5.99 +
------------------------------------------------------
plotting Figure 7.1 ...
plotting Figure 7.5 ...
7.3節の図7.5もここで出力する。平滑化残差のコレログラム。
英国ドライバー死傷者数。状態変数はレベル(確率的)、季節(確率的)、石油価格の回帰係数(確定的)、シートベルト法の回帰係数(確定的)。
> exercise.Chapter7.2 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.004034
hat(sigma^2_xi) = 2.681767e-04
hat(sigma^2_omega) = 1.572947e-09
hat(mu_1) = 6.781404
hat(gamma_1) = 0.008544
hat(beta_1) = -0.276738
hat(lambda_1) = -0.237591
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.004034
hat(sigma^2_xi) = 2.679640e-04
hat(sigma^2_omega) = 8.033134e-11
hat(mu_1) = 6.781380
hat(gamma_1) = 0.008543
hat(beta_1) = -0.276751
hat(lambda_1) = -0.237584
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 172.9417 ( 0.901) -311.88 (-1.624)
(by output) 71.4011 [-101.541] ( 0.372) -108.80 (-0.567)
KFAS (by pred.err) 174.2003 [ +1.259] ( 0.907) -314.40 (-1.638)
(by output) 184.2277 [ +11.286] ( 0.960) -334.46 (-1.742)
The CK Book 188.6443 [ +15.703] ( 0.983) -343.29 (-1.788)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (15) 18.675 22.36 +
r ( 1) 0.078 +-0.14 +
r (12) 0.068 +-0.14 +
homoscedasticity 1/H(59) 1.025 1.67 +
normality N 1.444 5.99 +
------------------------------------------------------
plotting Figure 7.2 ...
plotting Figure 7.3 ...
plotting Figure 7.4 ...
英国ドライバー死傷者数。状態変数はレベル(確率的)、季節(確定的)、石油価格の回帰係数(確定的)、シートベルト法の回帰係数(確定的)。
> exercise.Chapter7.3 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.004034
hat(sigma^2_xi) = 2.680726e-04
hat(mu_1) = 6.781414
hat(gamma_1) = 0.008542
hat(beta_1) = -0.276736
hat(lambda_1) = -0.237587
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.004035
hat(sigma^2_xi) = 2.680042e-04
hat(mu_1) = 6.781378
hat(gamma_1) = 0.008543
hat(beta_1) = -0.276752
hat(lambda_1) = -0.237584
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 172.9417 ( 0.901) -313.88 (-1.635)
(by output) 71.4011 [-101.541] ( 0.372) -110.80 (-0.577)
KFAS (by pred.err) 174.2019 [ +1.260] ( 0.907) -316.40 (-1.648)
(by output) 184.2277 [ +11.286] ( 0.960) -336.46 (-1.752)
The CK Book 188.1341 [ +15.192] ( 0.980) -344.27 (-1.793)
-----------------------------------------------------------------
plotting Figure 7.6 ...
plotting Figure 8.8 ...
plotting Figure 8.9 ...
making Figure 8.10 ...
plotting Figure 8.12a ...
plotting Figure 8.12b ...
8.5節の図8.8, 図8.9, 図8.10, 図8.12もここで出力する。標準化予測誤差系列、そのコレログラム、ヒストグラム、補助残差(標準化した平滑化レベル撹乱要素系列と標準化した平滑化不規則要素系列)。
英国インフレーション。4.4節のモデルに、2つのパルス干渉変数を投入する。状態変数はレベル(確率的)、季節(確率的)、2つのパルス干渉変数の回帰係数(確定的)。
> exercise.Chapter7.4 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 2.202838e-05
hat(sigma^2_xi) = 1.863874e-05
hat(sigma^2_gamma) = 4.332568e-07
hat(mu_1) = 0.005652
hat(gamma_1) = 0.002000
hat(lambda1_1) = 0.033350
hat(lambda2_1) = 0.042357
estimated by KFAS package:
hat(sigma^2_epsilon) = 2.201617e-05
hat(sigma^2_xi) = 1.864427e-05
hat(sigma^2_omega) = 4.334952e-07
hat(mu_1) = 0.005651
hat(gamma_1) = 0.002000
hat(lambda1_1) = 0.033339
hat(lambda1_1) = 0.042363
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 663.7023 ( 3.191) -1309.40 (-6.295)
(by output) 634.1082 [ -29.594] ( 3.049) -1250.22 (-6.011)
KFAS (by pred.err) 663.4805 [ -0.222] ( 3.190) -1308.96 (-6.293)
(by output) 682.4624 [ +18.760] ( 3.281) -1346.92 (-6.476)
The CK Book 687.4448 [ +23.742] ( 3.305) -1356.89 (-6.524)
-----------------------------------------------------------------
diagnostics of standardized prediction error (estimated by dlm):
------------------------------------------------------
stat value critical satisfied
------------------------------------------------------
independence Q (10) 11.786 15.51 +
r ( 1) 0.036 +-0.14 +
r ( 4) -0.068 +-0.14 +
homoscedasticity 1/H(67) 2.500 1.62 -
normality N 0.097 5.99 +
------------------------------------------------------
plotting Figure 7.7a ...
plotting Figure 7.7b ...
plotting Figure 7.7c ...
訳書p.85上部の数式に誤りがある(原著に対する著者らの訂正が反映されていない)。正しくは、
$y_t = t_1 + (t-1) v_t + \beta_1 \sum_{i=1}^{t-2} (t-i-1) x_i + \varepsilon_t$
ここで $\sum_{i=1}^{t-2} (t-i-1) x_i = 0$ ただし $t=1,2$のとき
4.3節で出力した。
2.3節で出力した。
図8.8, 図8.9, 図8.10, 図8.12は7.3節で、図8.11は4.3節で出力した。
図8.13は2.3節で、図8.14は3.4節で出力した。
訳書p.105以降は、英国ドライバー死傷者数の最初の169時点の分析。状態変数は、レベル(確率的)、季節(確率的)、説明変数(石油価格の対数)の回帰係数(確定的)。
最後の例として、英国ドライバー死傷者数の対数を、シートベルト法が施行される前の169時点だけについて分析する。この分析から得た予測を、1983年2月に英国でシートベルト法が施行された後のドライバー死傷者の実際の動きと比較する。(後略)
> exercise.Chapter8.6a ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.004140
hat(sigma^2_xi) = 2.524502e-04
hat(sigma^2_omega) = 2.367462e-09
hat(mu_1) = 6.746456
hat(gamma_1) = 0.005978
hat(beta_1) = -0.292133
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.004140
hat(sigma^2_xi) = 2.526106e-04
hat(sigma^2_omega) = 3.739316e-12
hat(mu_1) = 6.746495
hat(gamma_1) = 0.005964
hat(beta_1) = -0.292121
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 155.2805 ( 0.919) -276.56 (-1.636)
(by output) 53.7384 [-101.542] ( 0.318) -73.48 (-0.435)
KFAS (by pred.err) 155.2805 [ +0.000] ( 0.919) -276.56 (-1.636)
(by output) 158.5060 [ +3.226] ( 0.938) -283.01 (-1.675)
The CK Book 161.5061 [ +6.226] ( 0.956) -289.01 (-1.710)
-----------------------------------------------------------------
季節を確定的にする。さらに時点170以降について、石油価格の対数を使って予測する。
> exercise.Chapter8.6b ()
estimated by dlm package:
LL by dlm (to compare with exercise.Chapter8.6c()): -209.039
hat(sigma^2_epsilon) = 0.004140
hat(sigma^2_xi) = 2.526023e-04
hat(mu_1) = 6.746507
hat(gamma_1) = 0.005968
hat(beta_1) = -0.292114
estimated by KFAS package:
LL by KFAS (to compare with exercise.Chapter8.6c()): 158.506
hat(sigma^2_epsilon) = 0.004140
hat(sigma^2_xi) = 2.525567e-04
hat(mu_1) = 6.746491
hat(gamma_1) = 0.005964
hat(beta_1) = -0.292123
max loglikelihood and AIC:
-----------------------------------------------------------------
LogLik [diff] (/n) AIC (/n)
-----------------------------------------------------------------
dlm (by pred.err) 156.8938 ( 0.928) -283.79 (-1.679)
(by output) 53.7384 [-103.155] ( 0.318) -77.48 (-0.458)
KFAS (by pred.err) 156.8938 [ -0.000] ( 0.928) -283.79 (-1.679)
(by output) 158.5060 [ +1.612] ( 0.938) -287.01 (-1.698)
The CK Book 161.4934 [ +4.600] ( 0.956) -292.99 (-1.734)
-----------------------------------------------------------------
plotting Figure 8.15 (by predict() of KFAS) ...
ついでに、観察期間後について予測する問題ではなく観察期間に欠損があるデータについて予測する問題だと考え、やり直してみたのだけれど、フィルタ化シグナルの推定誤差分散の求め方がわからない...
> exercise.Chapter8.6c ()
LL by dlm (to compare with exercise.Chapter8.6b()): -209.039
LL by KFAS (to compare with exercise.Chapter8.6b()): 158.506
plotting Figure 8.15 (filtering by KFAS) ...
図8.16はこれまでのチャートの合成なので、省略する。
英国ドライバー死傷者数、ただし途中の2期間を欠損にする。状態変数はレベル(確率的)と季節(確定的)。4.3節と同じモデルである。
ここで大変困ったのだが、図8.19、季節要素の推定誤差分散の系列が再現できない。
考えられる可能性は以下の3つだ。
ひとことでいえばこういうことだ。平滑化シグナルの推定誤差分散とは、平滑化した状態変数の推定誤差分散の和なのでしょうか？ どなたか教えて下さい...
> exercise.Chapter8.7 ()
estimated by dlm package:
hat(sigma^2_epsilon) = 0.003862
hat(sigma^2_xi) = 0.000751
hat(mu_1) = 7.411790
hat(gamma_1) = 0.017583
estimated by KFAS package:
hat(sigma^2_epsilon) = 0.003863
hat(sigma^2_xi) = 0.000751
hat(mu_1) = 7.411789
hat(gamma_1) = 0.017584
plotting Figure 8.17 (by dlm) ...
plotting Figure 8.18 (by dlm) ...
plotting Figure 8.19 (var(seasonal) by dlm) ...
plotting Figure 8.19 (var(seasonal) by KFAS) ...
plotting Figure 8.19 (var(signal)-var(level) by KFAS) ...
plotting Figure 8.20 (by dlm) ...
plotting Figure 8.21 (by dlm) ...
CK本のこの章は駆け足で進んでいく感じなので、私も駆け足でやりたいと思う。dlmは省略し、KFASだけで試します。
英国死傷者数。前席死傷者数と後席死傷者数についてのモデルをつくる。状態変数は、レベル(確率的)、季節要素(確定的)、石油価格の回帰係数(確定的)、旅行キロ数の回帰係数(確定的)、シートベルト法の回帰係数(確定的)。
> exercise.Chapter9.4a ()
estimated by KFAS package:
estimated H matrix:
[,1] [,2]
[1,] 0.005401354 0.004448926
[2,] 0.004448926 0.008566328
estimated Q matrix:
[,1] [,2]
[1,] 0.0002560403 0.0002250350
[2,] 0.0002250350 0.0002321929
correlation between level disturbance:
[1] 0.9752466
regresssion coefficients:
[,1]
agKilo.agFront 0.1518435928
agPetrol.agFront -0.3075254484
abSeatbelt.agFront -0.3370416745
agKilo.agRear 0.5503489353
agPetrol.agRear -0.0856779645
abSeatbelt.agRear 0.0008966818
making Figure 9.1 ...
making Figure 9.2 ...
making Figure 9.3 ...
making Figure 9.4 ...
訳書p,123から。シートベルト法の回帰係数を前席だけに設ける。また、レベルの撹乱要素の分散行列をランク1に制約し、レベルを比例関係にする。
ここの箇所、大変難しかったのだが、KFASのマニュアルに載っているコード例を参考になんとか再現できた。嬉しいのでコードを載せておく。
oModel.KFAS <- SSModel(
cbind(agFront,agRear) ~ -1
+ agPetrol
+ agKilo
+ SSMregression(~-1+abSeatbelt, data=tsData, index=1)
+ SSMseasonal(period=12, sea.type='dummy')
+ SSMcustom(
Z = diag(2),
T = diag(2),
R = matrix(1,nrow=2,ncol=1),
Q = matrix(1),
P1inf = diag(2)
)
,
data = tsData,
H = matrix(NA, ncol=2, nrow=2)
)
# see manual of KFS()
funUpdate <- function(pars, model, ...){
model$H[,,1] <- exp(0.5*pars[1:2])
model$H[1,2,1] <- model$H[2,1,1] <- tanh(pars[3])*prod(sqrt(exp(0.5*pars[1:2])))
model$R[28:29] <- exp(pars[4:5])
## print(model$H[,,1])
## print(model$R[28:29])
model
}
# fitting
oFitted.KFAS <- fitSSM(
oModel.KFAS,
inits = c(-7, -7, 1, -1, -3)
updatefn = funUpdate,
method = "BFGS"
)
ポイントは以下のとおりである。
なるほどねー。頭いいなあ。 このやり方だとレベル撹乱要素はひとつしかないので、図9.5は省略する。
> exercise.Chapter9.4b ()
estimated by KFAS package:
estimated H matrix:
[,1] [,2]
[1,] 0.005444613 0.004381946
[2,] 0.004381946 0.008869446
estimated R matrix for custom disturbance:
custom1 custom2
0.01520784 0.01445783
shown as Q matrix:
custom1 custom2
[1,] 0.0002312785 0.0002198725
[2,] 0.0002198725 0.0002090290
regression coefficients:
[,1]
agPetrol.agFront -0.31498815
agKilo.agFront 0.17830894
agPetrol.agRear -0.08717566
agKilo.agRear 0.41825762
abSeatbelt.agFront -0.33639794
making Figure 9.6 ...
making Figure 9.7 ...
making Figure 9.8 ...
making Figure 9.9 ...
この章は計算例なし。さしたる意味もないが、いちおう図を再現しておく。
> exercise.Chapter10 ()
plotting Figure10.1 ...
plotting Figure 10.2 ...
plotting Figure10.3 ...
plotting Figure 10.4 ...
plotting Figure10.5 ...
plotting Figure 10.6 ...
plotting Figure10.7 ...
plotting Figure 10.8 ...
plotting Figure10.9 ...
plotting Figure 10.10 ...
ssfpackの説明なので、再現する必要がないだろう。なお、訳書p.153の(11.7)式に誤りがある(原著に対する著者らの訂正が反映されていない)。正しくは以下の通り (総和記号の下添字が t=d+1)。
$dLogLik = log L(y|\Psi) = -\frac{np}{2} log(2 \pi) -\frac{1}{2} \sum_{t=d+1}^n (log|F_t| + v' F^{-1}_t v_t)$
全体を振り返って思うに、状態空間モデルを使う際は、システム行列がどうなっているのかが頭に入っていないと大変に苦労する。文系ユーザであっても、行列表現の形できちんと勉強しておいた方が良い。以前出席した統計数理研究所のセミナーで、講師の偉い先生が「中味はただの行列演算だから、とにかく一度全部自分でプログラミングしてみなさい、パッケージに頼らずに」と述べていて、そんな無茶なと思ったものだが、いまにして思えばあのアドバイスは正しかった。
dlmとKFASを比べると、最初はdlmのほうが親切だと思ったのだが、章が進むにつれてKFASのほうに魅力を感じるようになった。dlmがクラスの人気者なら、KFASはとっつきにくいけど実は心の暖かい眼鏡っ娘、という感じである。(なにをいっているんだか)