# 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) }
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)
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)
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)
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)
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)
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)