ログインしてさらにmixiを楽しもう

コメントを投稿して情報交換!
更新通知を受け取って、最新情報をゲット!

Rで経済時系列コミュの野沢型の現在

  • mixiチェック
  • このエントリーをはてなブックマークに追加
これの1期遅れを修正するのと、独立した関数化で、いいのかどうか思案中。

観測ノイズにGARCHを使っている。
状態ノイズだが誤差の理由は全部速度同定がおかしいからだと押し付ける。
その他は単なる速度付きのトレンド回帰

library(tseries)
library(matlab)
library(signal)
Fs <- 1000;N <-150; time<-1:N; N2<-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)
R11 <-predict(garch(diff(yobs[1,])))[,1]^2

par(mfrow=c(3,1));#画面を3×1 に分割
plot(time,yobs,type="l");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),N2,N2)
G <- matrix(c(1,0,1,1),N2,N2)
H <- matrix(c(1,0),1,N2)
Q <- matrix(c(0,0,0,1),N2,N2)
R <- matrix(c(0),1,1)
Err <- 0
# ----------------------------------
# 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 10:N){
#
# prediction
x <- F %*% x # x(t+1|t) <- x(t|t)
Q[2,2] <- Err^2
P <- F %*% P %*% t(F) + G %*% Q %*% t(G) # P(t+1|t) <- P(t|t)
ypre[t] <- H %*% x
Spre[t] <- sqrt(H %*% P %*% t(H)+R11[t-1])
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R11[t-1]) # Kalman Gain
Err <- yobs[,t] - ypre[t]
x <- x + K * Err # x(t|t) <- x(t|t-1)
P <- P - K %*% H %*% P # P(t|t) <- P(t|t-1)
}

plot(time,yobs,type="l")
par(new=T)
plot(time,ypre,type="l",col="red",xlab="",ylab="",axes=F)
plot(time,Spre,type="l")

コメント(95)

とりあえず、モデリングだけする。

階差モデルで、T-1とT-2のノイズは、同じにして良いはずだが、取り合えず独立でやってみるべきか?

1)ノーマル

T = T-1 + S-1
S = S-1

X
|T|
|S|

F
|1,1|
|0,1|

G=F

H
|1,0|

2)トレンド階差、スピードノーマル

T = 2*T-1 - T-2 + S-1
T-1 = T-1
S = S-1

X
|T|
|T-1|
|S|

F
|2,-1,1|
|1, 0,0|
|0, 0,1|

G=F

H
|1,0,0|

3)トレンドノーマル、スピード階差

T = T-1 +2*S-1 - S-2
S = 2*S-1 - S-2
S-1 = S-1


X
|T |
|S |
|S-1|


F
|1,1, 0|
|0,2,-1|
|0,1, 0|

G=F

H
|1,0,0|

4)トレンド階差、スピード階差

T = 2* T-1 - T-2 + 2* S-1 - S-2
T-1 = T-1
S = 2* S-1 - S-2
S-1 = S-1

X
|T |
|T-1|
|S |
|S-1|

F
|2,-1,2, -1|
|1, 0,0, 0|
|0, 0,2, -1|
|0, 0,1, 0|

G=F

H
|1,0,0,0|

あら、Gを間違ったな。。

とりあえず、こうしてみます。


2)
G
|1,1,1|
|0,0,0|
|0,0,1|

3)
G
|1,1,1|
|0,1,1|
|0,0,0|

4)
G
|1,1,1,1|
|0,0,0,0|
|0,0,1,0|
|0,0,0,0|

今日はドライブで300キロ以上走ったので、眠くなったので寝ます。

お陰で、かなりぼんやりモデリングしている為、明日もう一度見直して
コーディングして見ることにしますね。

つか、まぁここまで出来る人は、コーディングそのものは容易い話。
>魚が釣れれば基準量なんか、どうでも良い話のようにも思う。

何を言ってるかと言うと、基準量的に考えれば、モデリング精度そのものは、当然、トレンドもスピードも
階差しないほうが有利に決まっている話なんです。(つか尤度そのものも落ちる筈)

けどまぁ考えてみれば、 それがトレンド継続のシグナルとして良いかどうかは、全く別の話のような
気もすると言うこと。

例えばGPSの位置検出も、階差モデルを使っているようだが、基準量的には、おかしな話。

まぁこれは、走っていてノイズで、ぴょんぴょん飛ばれる事のほうが、基準量よりも商品として
重要と言うことなんだろか?と思ったりしている。

とかくギザギザな相場データは、やっぱ階差モデルのほうが良いような気もすると言うこと。
さて、起きてきました。

まず、筋書きを決めておきます。

1)とっとと実装する。
2)データを数種類使いシャープレシオを出す。
3)この時にシステムノイズは、最適化された値を使う。

まずは、標準系(トレンド、スピード、共に階差なし)を、Nile意外でも
計測出来るようにして、それから他を実装します。
とりあえず標準型の計測用。

ナイルは赤字?のようです。

library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q22=1)
{
N <- length(yobs)
N2 <- 2
# 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),N2,N2)
G <- matrix(c(1,0, 1,1),N2,N2)
H <- matrix(c(1,0),1,N2)
Q <- matrix(c(Q11,0, 0,Q22),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[2]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*data[L]
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
#plot(Nile)
#par(new=T)
#plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)
レシオ計算が間違っていました。

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

これによると

> sratio(yobs,result$speed,30)
[1] -0.3980177

と、逆張りすれば、良いとの、お告げのようです。

昼寝も終わったんで、続きを頑張ります。
オマケ:

> dn<-diff(Nile)
> sum(diff(dn)^2)/sum(dn^2)
[1] 2.801595

と、こういう処理をした結果が、2を越える銘柄は逆張りが有効。

経済学では、DW(ダービン、ワトソン)比と言うらしいが、要するに
乱数拡散の場合は、2丁度になる。

で、

2の場合:相関なし
2以上:負の相関
2以下:正の相関
 
先ほどの計算は、diffで速度を取り出して、それのDW比を計算させて
いるので、ようするに逆張りが有利と言う、さっきの結果と一致する。
あちこち見て回ると、結構むつかしい説明がしてあるが

Vt= V1√T

なのが、ランダムウオーク。

で、このイコールの部分が崩れるのが相関が発生したと言うこと。

これは、標準偏差での話で、これを分散の形にすれば2になるだけ。

と言えば、ランジュバン風?
今度は、こんな結果が出てきた。

> sratio(yobs,result$speed,30)
[1] -0.4771214

2)番を書いてるつもり。

トレンドノイズはTn-1のもTn-2のも同等で扱った。
(たぶんこれで良いと思う)

同定されるスピードが極端に小さくなったが、間違いがあるのか、トレンドに食われた結果なのか
判断が付かない。

どうやって見分ければ良いのだろう?


---------------------------------------------------------------------------------------

library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q33=1)
{
N <- length(yobs)
N2 <- 3
# 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(2,1,0, -1,0,0, 1,0,1),N2,N2)
G <- matrix(c(1,0,0, 1,1,0, 1,1,0),N2,N2)
H <- matrix(c(1,0,0),1,N2)
Q <- matrix(c(Q11,0,0, 0,Q11,0, 0,0,Q33),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[3]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
#plot(Nile)
#par(new=T)
#plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)


まぁ相手はNileだし、70回ぶんなので、いいかと思っているけれど、えらい成績がいいね。
(無茶悪いと言う見方も出きるが、そういう場合は逆さまにみる。つまり逆張り)

たまたま、なんだろうか?



ははは。

間違い発見。使うノイズを、2番用じゃなく3番用を使っていました。

G <- matrix(c(1,0,0, 1,0,0, 1,0,1),N2,N2)

これが2番用

成績は、やや落ちて

> sratio(yobs,result$speed,30)
[1] 0.4370424

こんな感じ。

ただ気になるのが、だんだんとスピードが最後になって大きくなっているところ。

初期化がうまく行ってないんだろうか?
悪い所は後で見直すとして、これが3番

成績は、こんな感じ。

> sratio(yobs,result$speed,30)
[1] -0.3773699

確かにスピードの変化は、ゆったりしてきている。

library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q22=1)
{
N <- length(yobs)
N2 <- 3
# 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,0, 1,2,1, 0,-1,0),N2,N2)
G <- matrix(c(1,0,0, 1,1,0, 1,1,0),N2,N2)
H <- matrix(c(1,0,0),1,N2)
Q <- matrix(c(Q11,0,0, 0,Q22,0, 0,0,Q22),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[3]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
#plot(Nile)
#par(new=T)
#plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)
で、これが4番

> sratio(yobs,result$speed,30)
[1] -0.08480667

途中で、speedが、すっぽ抜けたような所が出ている。

なんか実装を間違っているんだろうか?

library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q33=1)
{
N <- length(yobs)
N2 <- 4
# 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(2,1,0,0, -1,0,0,0, 2,0,2,1, -1,0,-1,0),N2,N2)
G <- matrix(c(1,0,0,0, 1,0,0,0, 1,0,1,0, 1,0,0,0),N2,N2)
H <- matrix(c(1,0,0,0),1,N2)
Q <- matrix(c(Q11,0,0,0, 0,Q11,0,0, 0,0,Q33,0, 0,0,0,Q33),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[3]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
#plot(Nile)
#par(new=T)
#plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)

呼び方を、改めます。

ノーマル  -> T0S0
トレンド階差 -> T1S0
スピード階差 -> T0S1
両方 -> T1S1
で、T0S1ですが

result$speed[t] <- x[3]

を以下の用に修正しました。

result$speed[t] <- x[2]

成績ですが

> sratio(yobs,result$speed,30)
[1] -0.3773699

こんな感じ。


ノイズの接続パターンを変えてみる。

T0S1のGとQを、こんな風に変更

G <- matrix(c(1,0,0, 0,1,0, 0,0,0),N2,N2)
Q <- matrix(c(Q11,0,0, 0,Q22,0, 0,0,0),N2,N2)

> sratio(yobs,result$speed,30)
[1] -0.5331954

と、こんな結果。

一方で、T1S0を

G <- matrix (c(1,0,0, 0,0,0, 0,0,1),N2,N2)
Q <- matrix(c(Q11,0,0, 0,0,0, 0,0,Q33),N2,N2)

としてみると

> sratio(yobs,result$speed,30)
[1] -0.03410078

T1S1を

G <- matrix(c(1,0,0,0, 0,0,0,0, 0,0,1,0, 0,0,0,0),N2,N2)
Q <- matrix(c(Q11,0,0,0, 0,0,0,0, 0,0,Q33,0, 0,0,0,0),N2,N2)

とすると

> sratio(yobs,result$speed,30)
[1] -0.1119054



釣り道具?として使う場合、魅力があるとしたら、スピードをシグナルにする場合は、
S1なタイプだと思う。

もちろんS1にすれば良く釣れるのか?と言えば、当然違うとは思うが。

ようするに、S1で釣りたい場合は、S1に見合った養殖が必要だと言うことだと思うが
共和分などは、そのままでもS1向けだと思う。

成績を示したほうがいいが、番組規程に掛かると思うので、割愛することにします。

ただし、初期化ですが、今回やってるようなのは、とても実用性が低いと思う。

なんか、妥当な方法は無いものか?
T0S1の接続だが、このほうが自然か?

G <- matrix(c(1,0,0, 1,1,0, 0,0,0),N2,N2)

成績は

> sratio(yobs,result$speed,30)
[1] -0.4313325

と、落ちるが、フィットしたトレンドが自然な感じ。

たとえば、家電品のカルマンのシステムノイズは、どうやって決定するのかと言えば
要するに実験で求めた、最適化した初期値が、そのまま入っているだけ。

ノイズの精度に、鈍感なカルマンフィルタは、まぁこれでも実用になるが、相場でも
同じように一度ノイズを計っておけば、ダメと言う訳でもない。

まぁ、ださいから嫌なんだが。。
連休の成果として以下を上げておきます。

T0S1の、ノイズ接続いじりした、いまのとこの最終番。

観測ノイズに、garchを使い、システムノイズは、最尤してるわけで、まぁ経済時系列で
使うぶんには、これで文句言われる筋合いは、無いと思う。

トレンドはそのままで、速度には一階差を取ってあり、ゆっくり目の速度を拾い出すもの。

--------------------------

library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q22=1)
{
N <- length(yobs)
N2 <- 3
# 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,0, 1,2,1, 0,-1,0),N2,N2)
G <- matrix(c(1,0,0, 1,1,0, 0,0,0),N2,N2)
H <- matrix(c(1,0,0),1,N2)
Q <- matrix(c(Q11,0,0, 0,Q22,0, 0,0,0),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[2]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
plot(Nile)
par(new=T)
plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)
実は、ちょい本業でカルマントラックする用事があり、これを持っていって、ひどい目にあって
レビィさんに、ロケットのトラッキングのプログラムを、探して貰ったと。

どうも、尤度が、合わないのが釈然としなかったんですが、間違っていました。

これで試して、いいようなので、お仕事用にしました。

つか、今日は仕事で、一日カルマンフィルタいじりをしてました。

-----------------------------------
library(tseries)

T_S_G_GARCH <- function(yobs,Q11=1,Q22=1)
{
Q11 <-Q11^2
Q22 <-Q22^2
N <- length(yobs)
N2 <- 3
# 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,0, 1,2,1, 0,-1,0),N2,N2)
G <- matrix(c(1,0,0, 1,1,0, 0,0,0),N2,N2)
H <- matrix(c(1,0,0),1,N2)
Q <- matrix(c(Q11,0,0, 0,Q22,0, 0,0,0),N2,N2)
sink("garch-suppress.txt")
R <- predict(garch(diff(yobs)))[,1]^2
sink()
#
err2 <- 0
result <- data.frame(position=rep(0, N),speed=rep(0, N),err2=0)
# --------------
# Kalman Filter
# --------------
#
# initial value
x <- matrix(0,N2,1) # x(0|0)
P <- diag(N2) # P(0|0)
ypre <- matrix(numeric(N),1)
K <- matrix(0,N2,1)
#
d.yobs <- diff(yobs)
for(t in 1:10){
result$position[t] <- yobs[t]
result$speed[t] <- d.yobs[t]
}
for(t in 11: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)
ypre[t] <- H %*% x
#
# filtering
K <- P %*% t(H) %*% solve(H %*% P%*%t (H)+R[t-1]) # 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)
#
err2 <- err2+ (yobs[t] - ypre[t])^2
result$position[t] <- x[1]
result$speed[t] <- x[2]
}
result$err2 <- err2
result;
}


min_var <-function(par)
{
result <-T_S_G_GARCH(yobs,par[1],par[2])
return(result$err2[1])
}

sratio <- function(data,signal,start=1)
{
ratio <- rep(NaN,length(data))
for(L in start:length(data)){
ratio[L] <- sign(signal[L-1])*(data[L]-data[L-1])
}
ratio <- ratio[!is.nan(ratio)]
return(mean(ratio)/sd(ratio))
}

yobs <<- Nile
par <-optim(c(1,1),min_var,method="BFGS")$par
result <- T_S_G_GARCH(Nile,par[1],par[2])
plot(Nile)
par(new=T)
plot(result$position,type="l",col="red",xlab="",ylab="",axes=F)
sratio(yobs,result$speed,30)
要するにノイズ調整が、正常になっただけです。
>77番
>ロケットのトラッキングのプログラム
>尤度が、合わないのが釈然としなかったんですが、間違っていました。

 たしか、昨年末の宿題でしたが、
 添削結果の合否は如何だったのでしょうか??。(笑)
 
#---------
>T_S_G_GARCH。。。のコード

 Nileデータでは、相似形となります.
 しかし、NK225先物では下記の通りです
> sratio(yobs,result$speed,30)
[1] 0.08002135
カルマンのノイズは、分散なので、こんな風にして、調整するがいいみたい。

SSAなんかも時間掛かりますが、同じ風にも出来ると思いますよ。

> 79番 NK225先物での相似形出力
> sratio(yobs,result$speed,30)
(訂正報告)
[1] 0.09213023 逐次分散累計→Qに逐次代入
[1] -0.01593006 optim関数で誤差最適化→Qに代入
上記、いづれも相似形作図となりますが
   SR値に相違が生じます。

> 『カルマンのノイズは分散』なので、
> こんな風にして、調整するがいいみたい。
min_var <-function(par)
{
  result <-T_S_G_GARCH(yobs,par[1],par[2])
  return(result$err2[1])
}
par <-optim(c(1,1),min_var,method="BFGS")$par

 optim関数を誤差最適化に使用されたとは少々ビックリです。 
 77番T0S1(3)スピード階差のコードは、大変参考になりました。
 有難うございました。

 ところで、
 次はどのような課題を考えていらっしゃるのでしょうか??
 固有値の各定理をカルマン展開したいと思いますが、小生では実力不足。
 ただいま、ランダム行列理論で充電中といった近況です。
>optim関数を誤差最適化に使用された

カルマンのノイズは分散なので、これを覚えておくと、予測の精度を
上げていくのには、いいと思います。

けど、とんでもなく遅いのが泣き所ですよね。

>次はどのような課題を考えていらっしゃるのでしょうか??

あれこれ思うものは、あるにはあるんですが、よし!これにしよう!ってのが
無い段階です。

まぁ、ゆっくり考えて、おもしろそうなものを探して行きますね。

>まぁ、ゆっくり考えて、おもしろそうなものを探して行きますね。

要するに、なにが面白いか?というところが難しいわけです。

たまにチラっと言っていますが、相場は予想し、その精度を上げて勝つって事そのものも、ひとつの勝ちパターンでしかないわけです。

HFTなんかが話題になって以降は、あれで勝つ為には、たとえば低レーテンシー化なんてのが、どこの機関にも共通する、とても重要テーマのようで、こういうのは、明らかに、予測しないで勝つって、連中の考えることなわけです。

どっち道、現在の相場で勝つ主流は、予測して建てる事では無いわけです。

でも、ここでは予想の話ししかしないのは、これが基礎トレーニングだと思っているから。

なんの事なのかと言うと、野球選手も、ボクシングの選手も、毎日走りこみをするわけですが、自分にとっては、それに相当する部分だと思っているわけです。

で、走りこみに限らず、だいたいそういうのって、どれもこれもが単調で辛く、面白くも無いのばっかりなわけですが、それはスポーツ選手に限らず、私なんかの領域でも全くそうで、それこそ面白くもないわけです。

私自身は、自分の事を、結局は技術屋だと看做していて、価値あるものを作り出す事が自分の使命でも、やりたい事でもあると、常に思っておりまして、その為には、ない頭も、それなりに回るようにしておく必要が強くあり、そういうの基礎トレーニングに、予測が適していると思っているので、予測ばかりを、しかも、人前でやっているわけです。

そういう視点で、やることを探しています。
どの商売も、結局は期待値のある事を、やり続けているので、存続しているわけですが、じゃ予測で期待値を作り出しているのか?と言われれば、そこは全く違うわけです。

実は相場もそうで、たとえば旨い予測の方法があったとして、それを見せて喜ぶ機関があれば、それだけで今は2流だと判断して全く差し支えないです。

どうしてかって、予測よりも期待値を出す方法が色々とあるし、そういう方法を彼らは知っているからです。

最近ではCBOEのVIXも下がってきたのに、日本市場の出来高は情けないし、投資家心理は痛んだままだと言う人もいますが、これは全くのオカドチガイ。

なぜ北米の流動性が上がって来てるか?って、トレードスタイルそのものが変化しているのと、結局は北米の取引所がそれに対応出来ているのが、主な理由で、投資家心理で上がっている訳では全くないと思うけどな。。(ぼそ

ましてや最近になって流動性が株価を支えるなんて言ってるのは、言ってる奴がどんな奴らかを見れば、ははぁんそう言うことか。。って判る筈。

相場の本質は、ケーザイ学では当然ないし、予測でさえも、当然ありません。

そう強く、そう思っているのですけどね??。。(ぼそぼそ
じゃ、どういうのが、最近のスタイルなのよ?ってのを、チラっと話すと、出来高的に見ると、とても面白い傾向があるようです。

大きくふたつあって、流動性に貢献しまくるタイプと、流動性に貢献しなさすぎるタイプがあります。

G社さんなんかが見てて面白いのは、結局これの両方をやってるようで、両方で批判されているところ。

どっちか片方に気が付く人は結構いるものらしく、G社は流動性を上げろ!と批判する人と、流動性を下げろ!と批判する人が両方いるんですね。。(ぼそ

でまぁ、しなさすぎのタイプは、日本市場でもうまく機能するけれど、しすぎのタイプは、さほど旨くも行ってないようで。。

W大の学生さんグループが捕まったのは、早い話が、しすぎのタイプなんですが、これが学生さんグループでも機能したのは、結局は、しょぼい銘柄でやっていたからで、同じ事をG社さんが、その銘柄で、やったんじゃ食っては行けないわけで。。

個人が本気で投資家を目指すなら、貢献しないほうを検討するのが最善かとは思いますが。。(ぼそぼそ



ひとつだけハッキリ言っておける事があります。

今の北米市場で通用しなくなったのは、小規模なシテ戦の類です。

これは北米市場が健全化した為ではなく、設備的な問題で、小規模なシテが
通用しなくなった訳です。

日本市場も、数年のオーダーで見れば流動性も上がるし、シテ戦も相当少なくなるだろうと予想されます。

まぁこれがまた健全化とは、全く関係無しにですけどね。。(ぼそ
ここに書いたら、色々お調べになったようで、アローヘッドに興味を持たれている方も
いらっしゃるようですね。

まぁアローヘッドは、カナメな部分ではあるわけですが、まだまだ周辺が整備されて
いない事や、アローヘッドで、結局どうなるか?と言うのは、まだ未知数なわけで。。

見てると、どこもかしこも周辺整備絡みの整備が、しばらくは大変そうな様子ではあります。

まぁ、マーカンタイルなんかと同じ事で、1カイ2ヤリな、クラシカルスタイルな
デイトレーダーの類は、ほぼ確実に撲滅状態になると想定されます。

ご注意くださいませ。。(ぼそ
>シテ戦も相当少なくなるだろう

こうなったら、時間外のシンガで、海賊に混じって戦うしかないのか??

まぁ、こんな蛮勇のある人も居ないとは思うが。。(ぼそぼそ
>1カイ2ヤリな、

そういえば、昔のチャット仲間に、この手がいっぱいいたな。。

今でも生き伸びているので、連絡しといてやるか。。

まぁ板を目で見て売買するって時代は、確実に終わりなのかも。。(ぼそ
私の朝の日記。
全員に公開としておいたので、あっちに詳しいこと書くと、
温泉さんに、変な虫が付くとイカンと思い、
あちらへの返答は、やめておきました。。
>ネットワークの性能を大別すると、
>1)帯域幅2)レーテンシー、に、分類される。
>日本でHFTが、たいしてうまく行かない理由は2)が、しょぼいから。

レーテンシーってどういうことなんでしょう?


>こっちは、それ絡みな話が忙しくなりそうなんですけどね。。(ぼそ

あ〜 そういうお仕事なんですね。
>レーテンシー

そうだねぇ。

たとえば1Mのファイルを送ったら1秒掛かったとします。
では次に2Mのファイルを送ったら2秒掛かった。

ここで0Mのファイルを送ったら0秒掛かる。。ってのは計る前から当たり前なんだけど、
では極力0Mに近い所だと、実際は何秒掛かるの?と言うのがレーテンシーです。

ようするにこれも、一次関数みたいに、ファイルの大きさに応じて重くなるぶんと
元からある切片みたいな捉え方をするんですが、レーテンシーは、この切片です。

これって実は、いわゆる回線の太さ(帯域)を、でかくしても、単純には改善されない
部分なわけですが、ここが改善しなかったので、機関も個人も1円抜きは同じ条件
だったわけです。

平たく言えば1円抜きは、板の多い方に動くとするわけですが、まぁ今までは目で見ても
機械に見させても大差が無かった所が、レーテンシーが改善されると、機械が圧倒的に
有利になります。

アローヘッドも、この切片部分が大きく改善されるわけですが、こうなると手動の個人は
太刀打ち出来なくなるわけです。

じゃ、個人も、そこに機械を導入すれば?と言う話にしても、早くなるのは東証と機関の間なので
個人サイドが機械を持っても、個人と機関の間のぶんが、以前にまして差が大きく見える
ようになるわけです。

ふむじゃ、1円抜きはヤメて、予測重視な日ばかりにすれば?と言うことにしても、更に相関が
見られる時間幅も、実は短くなるようで、今まで指数よりも株のほうが予想は、やりやすかった
わけですが、少なくとも指数並に相関の消える時間は短くなるものと思われます。

つうわけで、これからの個人は、もう養殖系?しか残って無いのかな?と思う次第。
>もう養殖系?しか残って無いのかな

ムギュ。

ちょっとヘコむ話ですね。
やっぱし、もう一度、養殖を勉強 頑張らねば。。
これだと個人に不平等すぎるから、HFTを禁止する規制を行うと言う説もありますが
こんなの、どうせ、おおうそ。

ホントなのだとしたら、とっくに北米で規制してる筈で。。

あと、機械化が進むと、日本市場は、ストップ高/安とか、すぐに付いちゃって、危険
極まりないような気もするんですけど、どうなるんでしょう?

まぁ、暫くは様子をみたほうが、いいでしょう。

けど実際問題は、後半年は、周辺システムが全部は、出揃わないから、完璧な
本領発揮とまでは行かないはず。

まぁメイン部分は、すぐに動くわけだけど、それまでに1円抜きを楽しんでおく??

ログインすると、残り55件のコメントが見れるよ

mixiユーザー
ログインしてコメントしよう!

Rで経済時系列 更新情報

Rで経済時系列のメンバーはこんなコミュニティにも参加しています

星印の数は、共通して参加しているメンバーが多いほど増えます。