BugTrack-状態空間モデリング/3

状態空間モデリング - 記事一覧

「co2」データを使った状態空間表現の例 - 季節成分を抽出する

• 投稿者： みゅ
• カテゴリ： なし
• 優先度： 普通
• 状態： 完了
• 日時： 2008年06月12日 16時59分04秒

内容

「co2」データを使用し、季節成分を抽出する

概要

• 季節成分は１２ヶ月で１周期になっている．季節成分をStとすると

St+1 + St + St-1 + St-2 + St-2 + ... + St-10 = ε3

• である（つまり１２ヶ月の合計が平均０）．なので

観測方程式 ： yt+1 = at+1 + St+1 + ε1
状態方程式 ： at+1 = at + ε2
              St+1 = -(St + St-1 + St-2 + ... + St-10) + ε3
              St   =   St
              St-1 =        St-1
              ...

コード（フィルタリングまで）

library(sspir)
Fmat <- matrix(c(1.0, 1.0, rep(0, 10)))
Vmat <- matrix(0.0, 1, 1)
v <- c(1.00)
diag(Vmat) <- v
Gmat <- matrix(0.0, 12, 12)
Gmat[1,1] <- 1.0
Gmat[2,2:12] <- -1.0
Gmat[3:12,2:11] <- diag(1,10)
w <- c(u=0.5^2, s=1.0^2)
Wmat <- diag(0,12)
diag(Wmat)[1] <- w[1]
diag(Wmat)[2] <- w[2]
m0 <- matrix(double(1), 1, 12)
m0[1,1] <- as.numeric(co2)[1]

ssm_co2 <- SS(y=matrix(as.numeric(co2)),
    Fmat=function(tt,x,phi) return(Fmat),
    Vmat=
        function(tt,x,phi) {
            Vmat0 <- Vmat
            diag(Vmat0) <- phi$v
            return(Vmat0)
        },
    Gmat=
        function(tt,x,phi){
            Gmat0 <- Gmat
            Gmat0
        },
    Wmat=
        function(tt,x,phi){
            Wmat0 <- Wmat
            diag(Wmat0)[1] <- phi$w[1]
            diag(Wmat0)[2] <- phi$w[2]
            Wmat0
        },
    phi=list(v=v, w=w),
    m0=m0,
    C0=diag(1,12)
)
ssm_co2_f <- kfilter(ssm_co2)

par(mfrow=c(2,1))
plot(ssm_co2_f$y, ty="l")
lines(ssm_co2_f$m[,1], ty="l", col="red")
plot(ssm_co2_f$m[,2], ty="l")
par(mfrow=c(1,1))

• トレンド成分

• 季節成分

• 系列の最初のほうで、季節成分がうまく分離されていない．これは状態変数の時点０での期待値【ｍ０】と分散【ｃ０】が適切に設定されていないため．これらを適切に設定する方法、あるいはフィルタリングを開始することをフィルタリングの初期化という．
• 初期化についてはまたの機会に・・・

平滑化

• 状態方程式の分散行列の下の１０行が０になっているが、このままでは平滑化でエラーが出る．
• 小さな正の数を与えることで、平滑化ができるようになる

コメント