Rで再現する仕事の都合で仕方なく状態空間モデルについて勉強していたのだけれど(なぜ私がこんな目に)、仕事で使うためには自分で計算できるようにならなければならない。(以下「CK本」)の著者らは、すべての事例についてデータとプログラムを公開している。ありがたいことであります。ssfpackという耳慣れないソフトを使わなければならない。新しいソフトの使い方を覚えるのは大変に面倒だ。できれば普段使っているソフトで済ませたい。汗と涙の甲斐あって、すべての章についていちおう再現できたので、ここに載せておくことにする。実践的な状態空間モデリング感想Rプログラム紹介Rのプログラムはこちらに置いておきますが、全ての計算例を再現しているので、それなりに長い。雰囲気をご覧頂くために一例だけご紹介します。この本を通して用いられている、英国ドライバー死傷者数(対数)の時系列についての分析である。この時系列に、いわゆるローカル・レベル・モデルをあてはめる。状態変数はレベル(死傷者数の潜在的な大きさ)と季節要素(12ヶ月)、ともに確率的に変動させる。\sigma^2_\omega)$4.2節は、Rプログラムのexercise.Chapter4.2()という関数で再現している。ライブラリを読み込み、定数を定義。gCKLogLikはCK本に掲載されている最大対数尤度で、あとで比較表を作るために使う。今回の方針として、用いる時系列データは要不要を問わず必ずtsクラスに放り込むことにしている。tsDataが英国ドライバー死傷者数の対数の時系列。nが時点数(=192)。length(tsData)nStateVarとnHyperParamは、それぞれ状態変数の数とハイパーパラメータの数。分析結果のまとめの際に用いる。季節が11個のダミー変数で表現されるので、状態変数はレベルと季節で計12個。\sigma^2_\omega$の3個である。agInitはハイパーパラメータの最尤推定の際に与える初期値。ここまでが前置きです。cat("\n")dlmパッケージによるモデリング。モデルを表現する関数funModel()を定義する。m1はローカル・レベル・モデルを表している。$\sigma^2_\xi$である。どちらも未知パラメータとして指定している。m2はこれに付け加える季節変動モデル。dWは撹乱要素の分散行列の対角要素で、季節は11個のダミー変数で表現するので長さ11のベクトルになる。0,0,0,0,0,0,0,0,0,0)$となる。季節要素も確率的に変動すると考え、撹乱要素の分散を未知パラメータとして指定する。不規則要素の分散はさっき指定したので、ここのdVは0でよい。ここがdlmでのもっとも大事な部分なので、誤りがないか念入りに確認する。funModel()の挙動を調べたコードが残ってますね、ははは。print(funModel(log(c(123,456,789))))パラメータを最尤推定。funModel()に与え、モデルを完成する(oFitted.DLM)。funModel(oMLE.DLM$par)このモデルを使ってtsDataをフィルタ化する(oFiltered.DLM)。lPredErr.DLM)。モデルを使ってtsDataを平滑化する(oSmoothed.DLM)。あとで使うために、平滑化状態系列を抽出する(tsSmostate.DLM)。dlmSmoothの返し値のsは多変量時系列のtsオブジェクトになっているので便利なのだが、観察時点の1期前から値が入っている点に注意。dlmパッケージに入っているdropFirst()で最初の時点を削除している。2種類求めておく(詳細は後述)。推定したハイパー・パラメータと初期状態を表示する。この部分の出力は下記の通り。順に、不規則要素の分散$\sigma^2_\epsilon$、レベル撹乱要素の分散$\sigma^2_\xi$、季節撹乱要素の分散$\sigma^2_\omega$、時点1のレベル$\mu_1$と季節要素$\gamma_1$、の推定値である。0.017273おつぎはKFASパッケージによるモデリング。モデルを表現するオブジェクトoModel.KFASを定義する。KFASでは、未知パラメータであることをNAで表現する。右辺のSSMtrend()はローカル・レベル・モデルを表している。Qが不規則要素の分散行列のリストで、ここではリストの要素はひとつ(状態変数がレベルだけだから)。その行列のサイズは1x1、中身は$\sigma^2_\xi$である。SSMseasonal()はこれに付け加える季節要素を表している。Qをサイズ1x1の行列にしていることに注意。$\sigma^2_\omega$となる。この2つを未知パラメータに指定。H行列は不規則要素の分散行列、つまり$\sigma^2_\epsilon$である。これも未知パラメータに指定。KFASにとってこの部分が肝になるので、念入りに調べる必要がある。oModel.KFASをprint()するといろいろ教えてくれるので、メモをとっている。季節ダミー2,...)の順に並んでいる。レベル撹乱要素の分散行列はサイズ2x2になっている。(NA,NA)パラメータを最尤推定。うまくいったら、この推定されたモデルをKFS()に与える。KFS()はデフォルトでフィルタ化と平滑化の両方を行ってくれる。あとで使うために予測誤差を抽出する(agStdPredErr.KFAS)。推定したハイパー・パラメータと初期状態を表示する。この部分の出力は下記の通り。dlmとほぼ同じ推定結果が得られている。0.017272dlmとKFASから得られた最大対数尤度と、CK本に載っている最大対数尤度を比較する。比較表を表示する関数sub.ShowLogLik()は別途に定義してある。この部分の出力は下記の通り。一行目はCK本に基づきdlmの予測誤差系列とその分散の系列から求めた値。時点数で割ると-1.721。二行目はdlmが出力した対数尤度で、なんと99.193も小さい(小さくなる理由は後述)。三行目と四行目はKFASでの結果。最後の行ではCK本に載っている値(=ssfpackで求めた値)を示している。(-1.718)-----------------------------------------------------------------残差診断の結果を表示する。関数sub.ShowDiagnostics()は別途に定義してある。この部分の出力は下記の通り。CK本とほぼ同じ結果が出力されるが、例題によっては均一分散性の臨界値がCK本と異なることがある(その理由は後述)。+------------------------------------------------------最後に、図の作成。それぞれの関数は別途定義してある(どれもごく単純な関数である)。)}最初に出力される図4.6はこんな感じ。赤線が時系列の背後にあると推定されるレベルである。意味ありげな減少が2回生じているのがわかる。