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

!!!「co2」データを使った状態空間表現の例
*投稿者： みゅ
*カテゴリ： なし
*優先度： 普通
*状態： 完了
*日時： 2008年06月12日 16時07分17秒

!!内容
*Rに含まれている「co2」データを使って、モデル化を提示
!元の「co2」データ
 plot(co2)
{{ref_image co2_1.png}}
*季節的な周期成分があるように見える．これらは後で考慮するとして、今は一番単純なモデル化を行う
 観測方程式 ： yt = at + e1
 状態方程式 ： at+1 = at + e2
 
 library(sspir)
 ssm_co2 <- SS(y=matrix(as.numeric(co2)),
     Fmat=function(tt,x,phi) return(matrix(1)),
     Vmat= function(tt,x,phi) return(matrix(1)),
     Gmat= function(tt,x,phi) return(matrix(1)),
     Wmat= function(tt,x,phi) return(matrix(0.1)),
     phi= NULL,
     m0=matrix(1),C0=matrix(1)
 )
*ノーテーションが文献によって異なるが、ここでは
**Fmat ： 観測方程式を形成する行列
**Vmat ： 観測方程式の分散行列
**Gmat ： 状態方程式を形成する行列
**Wmat ： 状態方程式の分散行列
**m0/C0は、状態ベクトルaの初期の分布
***これらの値のいくつかは未知であるが、このような場合に系列を開始する方法は、別の機会に触れる（フィルタリングの初期化）
***今の場合は初期の状態ベクトルが平均「１」、分散「１」の分布に従うとしている
*上の状態方程式からわかるように、状態ベクトルは１階のランダムウォークに従うとしている
!!!フィルタリング
 ssm_co2 <- SS(y=matrix(as.numeric(co2)),
     Fmat=function(tt,x,phi) return(matrix(1)),
     Vmat= function(tt,x,phi) return(matrix(1)),
     Gmat= function(tt,x,phi) return(matrix(1)),
     Wmat= function(tt,x,phi) return(matrix(0.1)),
     phi= NULL,
     m0=matrix(1),C0=matrix(1)
 )
 ssm_co2_f <- kfilter(ssm_co2)
 plot(ssm_co2_f$y, ty="l")
 lines(ssm_co2_f$m, ty="l", col="red")
{{ref_image co2_2.png}}
!!!平滑化
 ssm_co2_fs <- smoother(ssm_co2_f)
 plot(ssm_co2_f$y, ty="l")
 lines(ssm_co2_f$m, ty="l", col="red")
 lines(ssm_co2_fs$m, ty="l", col="blue")
{{ref_image co2_3.png}}
*ブルーのラインが平滑化後の系列
!!!ちなみに
*初期状態ベクトルの平均を以下のように変更すると
 ssm_co2 <- SS(y=matrix(as.numeric(co2)),
     Fmat=function(tt,x,phi) return(matrix(1)),
     Vmat= function(tt,x,phi) return(matrix(1)),
     Gmat= function(tt,x,phi) return(matrix(1)),
     Wmat= function(tt,x,phi) return(matrix(0.1)),
     phi= NULL,
     m0=matrix(as.numeric(co2)[1]),C0=matrix(1)
 )
 ssm_co2_f <- kfilter(ssm_co2)
 ssm_co2_fs <- smoother(ssm_co2_f)
 plot(ssm_co2_f$y, ty="l")
 lines(ssm_co2_f$m, ty="l", col="red")
 lines(ssm_co2_fs$m, ty="l", col="blue")
{{ref_image co2_4.png}}
!!!ちなみに（２）
*状態方程式の分散行列を１０倍にすると（Wmatを０．１→１に）
 ssm_co2 <- SS(y=matrix(as.numeric(co2)),
     Fmat=function(tt,x,phi) return(matrix(1)),
     Vmat= function(tt,x,phi) return(matrix(1)),
     Gmat= function(tt,x,phi) return(matrix(1)),
     Wmat= function(tt,x,phi) return(matrix(1)),
     phi= NULL,
     m0=matrix(as.numeric(co2)[1]),C0=matrix(1)
 )
 ssm_co2_f <- kfilter(ssm_co2)
 ssm_co2_fs <- smoother(ssm_co2_f)
 plot(ssm_co2_f$y, ty="l")
 lines(ssm_co2_f$m, ty="l", col="red")
 lines(ssm_co2_fs$m, ty="l", col="blue")
{{ref_image co2_5.png}}
!コメント
//{{comment}}