//このページが表示される方は、URLから「action=SOURCE&」を削除してみてください
[[状態空間モデリング - 記事一覧]]

!!!欠測値の扱い
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2008年06月20日 11時20分51秒
//{{bugstate}}
!!内容
R言語の状態空間を扱える【sspir】パッケージで、欠測値を扱う
!!!概要
*参考書「[状態空間モデリングによる時系列分析入門|http://www.amazon.co.jp/%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AA%E3%83%B3%E3%82%B0%E3%81%AB%E3%82%88%E3%82%8B%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90%E5%85%A5%E9%96%80-J-%E3%83%80%E3%83%BC%E3%83%93%E3%83%B3/dp/4916092457/ref=sr_1_1?ie=UTF8&s=books&qid=1213841649&sr=8-1]」によると、状態空間表現で欠測値を扱うのは非常に容易とのこと
*方法は
 v = ｙ - a
*の箇所で、ｖを０とおけばよいだけだから
!!!sspir
*sspirでは欠測値の取り扱いが現状では出来ない
*が、元の関数をちょっといじってやるだけですぐにできる
 ｙがNAだったら、参考書で言うｖを０にする
*これだけ
*実際のコードと、「ｃｏ２」データにおける例は後日・・・
!!!コード
*sspirパッケージの関数定義で【kfilter.SS】という定義がある．この中で時系列方向に、フィルタリングの逐次計算を行っているのが【filterstep】関数
 > filterstep
 function (y, Fmat, Gmat, Vt, Wt, mx, Cx)
 {
     a <- Gmat %*% t(mx)
     R <- Gmat %*% Cx %*% t(Gmat) + Wt
     f <- t(Fmat) %*% a
     Q <- t(Fmat) %*% R %*% Fmat + Vt
     e <- y - f
     A <- R %*% Fmat %*% mysolve(Q)
     m <- a + A %*% e
     C <- R - A %*% Q %*% t(A)
     if (length(y) > 1)
         loglikterm <- log(dmvnorm(as.numeric(y), as.numeric(f),
             Q))
     else loglikterm <- -0.5 * (log(2 * pi) + log(Q) + (y - f)^2/Q)
     list(m = m, C = C, loglikterm = loglikterm)
 }
 <environment: namespace:sspir>
*この関数の【ｅ】変数が参考書で言う「ｖ」．なので該当箇所を
       if( any(is.na(y)) ) {
         e <- f
         e[] <- 0
       }else {
         e <- y - f
       }
         A <- R %*% Fmat %*% mysolve_c(Q)
       if( is.na(y) ) {
         A[] <- 0                       --- ①
       }
*のように書き換えてあげればまったく問題なく．
*この場合、欠測値をＮＡでｙに与えておく
**いま気がついたのだけれども、なぜ①のような記述をしたのか思い出せない・・・
**【ｅ】が０なんだから【Ａ】がどんな値をとろうが
     m <- a + A %*% e
**には、影響がないように思うのだが・・・．それとも、次の行の
     C <- R - A %*% Q %*% t(A)
**に関係あるんだったか・・・．条件のところが
       if( any(is.na(y)) ) {
**でなく
       if( is.na(y) ) {
**になってるから、また別のケースを考慮してるのだろう、きっと・・・
!!!例（[「co2」データを使った状態空間表現の例|http://r-project.matrix.jp/R/wiki.cgi?page=BugTrack%2D%BE%F5%C2%D6%B6%F5%B4%D6%A5%E2%A5%C7%A5%EA%A5%F3%A5%B0%2F2]より）
*フィルタリングまで
{{ref_image co2_6_1.png}}
↑状態変数が１階のランダムウォークなので、欠測値があっても直近の値で補完されるだけ
*平滑化
{{ref_image co2_6_2.png}}
!!!季節調整込み
*フィルタリング
**もとの系列
{{ref_image co2_6_3.png}}
**トレンド ＋ 季節成分
{{ref_image co2_6_4.png}}
↑季節調整成分を入れると、このように季節成分で補完される
!!コメント
//{{comment}}