して求めた2種類のフィルタ化シグナルとその推定誤差分散を比べてみたところ、dlmの推定誤差分散が滅茶苦茶大きくなってしまった。下図に1975年以降を示す(それ以前はもっと差が大きい)。上記の推定誤差分散の求め方が間違っているのだろう。いうわけで、dlmを使うのはあきらめ、KFASを使って図8.15を描きなおしてみた。観察期間の末期から描き始めている。...図8.16はこれまでのチャートの合成なので、省略する。英国ドライバー死傷者数、ただし途中の2期間を欠損にする。状態変数はレベル(確率的)と季節(確定的)。4.3節と同じモデルである。ここで大変困ったのだが、図8.19、季節要素の推定誤差分散の系列が再現できない。ここはdlmとKFASの両方で描いてみた。どちらにせよ、CK本のような図にはならない。そもそもこのモデルは確定的季節要素なのに、平滑化季節要素の推定誤差分散の系列が、欠損の期間だけこんな風に変化したりするものだろうか?著者らのssfpackプログラムをよく見ると、季節要素の推定誤差分散の系列を描いているというよりも、「シグナルの推定誤差分散からレベルの推定誤差分散を引いた系列」を描いているようである。KFASはシグナルの推定誤差分散を出力してくれるので、シグナルの推定誤差分散からレベルの推定誤差分散を引いた系列を描いてみたら、CK本とそっくりの図が得られた。うむむむむ...考えられる可能性は以下の3つだ。このモデルでは、シグナルの推定誤差分散とはレベルの推定誤差分散と季節要素の推定誤差分散の和だ。dlmとKFASの出力において季節の推定誤差分散だと私が思っているものは、実は季節要素の推定誤差分散ではない。このモデルでは、シグナルの推定誤差分散はレベルの推定誤差分散と季節要素の推定誤差分散の和ではない。8.19の説明は誤っており、本当は「シグナルの推定誤差分散からレベルの推定誤差分散を引いた系列」と書くべきだ。私が箸にも棒にもかからないくらいになにもかも誤っている。こんなことに時間を費やすのはやめて、いますぐに寝るべきだ。ひとことでいえばこういうことだ。例CK本のこの章は駆け足で進んでいく感じなので、私も駆け足でやりたいと思う。後席死傷者数についてのモデルをつくる。状態変数は、レベル(確率的)、季節要素(確定的)、石油価格の回帰係数(確定的)、旅行キロ数の回帰係数(確定的)、シートベルト法の回帰係数(確定的)。「制御系列」と呼んでいる。実験研究でいうところの「処理条件」と「統制条件」にあたる言葉である。シートベルト法の効果を調べるという観点から、後席は統制条件とみなすことができる、ということであろう。シートベルト法の回帰係数を前席だけに設ける。レベルの撹乱要素の分散行列をランク1に制約し、レベルを比例関係にする。ここの箇所、大変難しかったのだが、KFASのマニュアルに載っているコード例を参考になんとか再現できた。嬉しいのでコードを載せておく。)ポイントは以下のとおりである。SSModel()では、レベル以外のすべての要素は通常通りに記述する。レベルはSSMtrend()ではなくSSMcustom()で記述する。状態変数は前席のレベルと後席のレベルの2つ。遷移行列も2次の単位行列、P1infも2次の単位行列(初期状態が未知であることを表す)。レベル撹乱要素をひとつしか作らず、その分散を1に固定してしまう。従ってQ行列は1次の単位行列。その係数行列(というか、因子負荷行列)であるR行列は2行1列となる。Rの要素の値はなんでもいいのだが、とりあえず1にしておく。SSMcustomで追加したR行列は、モデル自体のR行列の28行目と29行目になる。custom2と呼ぶ。fitSSM()のなかで、updatefunにユーザ定義関数を指定する。ユーザ定義関数を使うのは、SSModel()で指定できる未知パラメータがQ行列とH行列に限定されているからだ。