温泉気分で♪

twitterしています。 https://twitter.com/nozawa_onsen

全体表示

[ リスト ]

カルマン講座のゆくえ

イメージ 1
結局、経済的なグループと、ちょっと工学よりをわけました。

で、工学よりは、少し進んでいますが経済よりは進んでいません。
相棒の都合が付かず遅れています。

そもそも経済予測で期待値を得るのって、単純に
カルマンフィルタを使えば、バリバリに予測出来るって
ものでもありませんからね。

とは言え何もやってないのも何なのでトレンドモデルのベース
だけ作ってみました。

まだ初期値やノイズなどは、調整していませんが参考まで。

F = |1,1|
    |0,1|

G = F

H = |1,0|

と言うモデルです。

tn = tn-1+Δtn+vn

トレンド(大まかな値)は、一期前のトレンドとトレンドの勾配に
ノイズが乗ったものと言うモデルです。

これでノイズの尤度を上げれば、それっぽくなるかな?ってね。

つか、どういう風なノイズを載せれば、トレンドって綺麗に
浮かび上がるんだろ?って話。

さて、どうなんでしょうかね?



library(signal)
Fs <- 1000;N <-150; time<-1:N; N2<-2; N3 <-2
a <-matrix(c(1.4 ,-0.98),nrow=1,ncol=2)
a2<-matrix(c(-0.98, 1.4),nrow=2,ncol=1)
aa<-c(1 , -a)

Nrn<-rnorm(N)
yobs<-matrix(filter(1,aa,1*Nrn),1)

par(mfrow=c(3,1));#画面を3×1 に分割
plot(time,yobs,type="l",ylim=c(-20,20));axis=c(0, N,-20, 20)

# State Space Representation
#
# x(t+1) = F x(t) + G w(t) # cov{w(t) w(t)} = {Q 0}
# y(t) = H x(t) + v(t) # {v(t),v(t)} = {0 R}
F <- matrix(c(1,0,1,1),2,2)
G <- matrix(c(1,0,1,1),2,2)
H <- matrix(c(1,0),1,2)
Q <- matrix(c(1,0,0,1),2,2)
R <- matrix(c(10),1,1)
# ----------------------------------
# Kalman Filter
# ----------------------------------
#
# initial value
x <- matrix(0,N2,1); # x(0|0)
P <- diag(N2); # P(0|0)
ypre <- matrix(numeric(N),1)
Spre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
for(t in 3:N){
#
# prediction
x <- F %*% x; # x(t+1|t) <- x(t|t)
P <- F %*% P %*% t(F) + G %*% Q %*% t(G); # P(t+1|t) <- P(t|t)
#H<-matrix(c(yobs[,t-1],yobs[,t-2]),1,2)
ypre[t] <- H %*% x
Spre[t] <- sqrt(H %*% P %*% t(H)+R)
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R); # Kalman Gain
x <- x + K * (yobs[,t] - ypre[t]); # x(t|t) <- x(t|t-1)
P <- P - K %*% H %*% P; # P(t|t) <- P(t|t-1)
}

plot(time,yobs,type="l",ylim=c(-20,20))
par(new=T)
plot(time,ypre,type="l",col="red",xlab="",ylab="",ylim=c(-20,20))
plot(time,Spre,type="l")
plot(time,Spre,type="l")


この記事に

閉じる コメント(0)

コメント投稿

顔アイコン

顔アイコン・表示画像の選択

名前パスワードブログ
絵文字
×
  • オリジナル
  • SoftBank1
  • SoftBank2
  • SoftBank3
  • SoftBank4
  • docomo1
  • docomo2
  • au1
  • au2
  • au3
  • au4
投稿

.


みんなの更新記事