In this section we show how to modify the extended kalman filter and
smoother to cater for a non-linear system equation. The resulting code
is applied on the advertisement data described in Section 2.7.2.
Lines that have been modified are marked with ###.
# The data
Yt <- round(c(40,41,31,40,45,44,39,50,32,42,33,24,25,32,28,
25,36,38,36,29,43,34,42,50,43,43,52,45,30,55,
33,32,39,32,30,44,27,44,30,32,30,NA,NA,NA,33,
48,40,44,40,34,37,37,23,30,21,23,22,25,23,14,
21,16,19,07,26,16,21,07,22,10,15,15,22,11,14)/100*66,0)
nt <- rep(66,75)
nt[is.na(Yt)]<-0
Xt <- matrix(
c( 5, 0, 20,780,610,515,110, 0, 0, 0, 0, 0, 0, 0, 0,
150,460,370,145,120,200,340,440,380,390,500, 10, 60,385,350,
315,330, 35, 0,280,290,340,220, 50, 0, 0, 10, 85,465,510,
550,230,460, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)/100,
ncol=1)
# Priors from West and Harrison (1997) page 542
m0 <- c(0.1,0.85,0.90,0.02,0.30)
C0 <- matrix(c( 6.25, 6.25, 0, 0, 0,
6.25, 406.25, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 2.25, 0,
0, 0, 0, 0, 100),
ncol=5, byrow=T)*0.0001
# Design matrix
Ft <- function(i,x,phi) c(66,0,0,0,66)
# Generate instance of class ssm
# Note that evolution variances and evolution transfer matrices
# are left unspecified
ss <- ssm(Ft=Ft,
m0=m0,
Xt=Xt,
C0=C0,
Yt=Yt,
nt=nt,
fam="binomial",
link="identity")
# Evolution transfer function
gt <- function(zt, xt){
# z is five dimensional vector
# xt is explanatory - here TWR
z1 <- zt[1]
z2 <- zt[2]
z3 <- zt[3]
z4 <- zt[4]
z5 <- zt[5]
result <- c(z1,z2,z3,z4,z2-z1-(z2-z1-z3*z5)*exp(-z4*xt))
result
}
# Derivative of gt
Gt <- function(zt, xt){
# zt is five dimensional vector
# xt is explanatory - here TWR
z1 <- zt[1]
z2 <- zt[2]
z3 <- zt[3]
z4 <- zt[4]
z5 <- zt[5]
result <- matrix(
c( 1, 0, 0, 0, exp(-z4*xt)-1,
0, 1, 0, 0, 1-exp(-z4*xt),
0, 0, 1, 0, z5*exp(-z4*xt),
0, 0, 0, 1, xt*(z2-z1-z3*z5)*exp(-z4*xt),
0, 0, 0, 0, z3*exp(-z4*xt)), ncol=5, byrow=F)
result
}
#####################################################
# MODIFIED KALMAN FILTER FOR BINOMIAL DATA
# IDENTITY LINK, NON-LINEAR SYSTEM EQUATION
#####################################################
modified.filter <- function(ssm, Gt, gt)
{
# require
# ssm: object of class ssm
# gt: state function # NEW
# Gt: Hessian of gt # NEW
# Returns an updated ssm object
# mt: Filtered states
# Ct: Filtered MSE
# Rt: Prior variances (needed in the Kalman smoother)
# llh: Log-likelihood
# Extract information from ssm object
Ft <- ssm$Ft
# Gt <- ssm$Gt # Not needed
# Wt <- ssm$Wt # Not needed
Xt <- ssm$Xt
phi <- ssm$phi
m0 <- ssm$m0
C0 <- ssm$C0
Yt <- ssm$Yt
nt <- ssm$nt
m.start <- ssm$smoothed$m.tilde
# Dimensions
n <- length(Yt)
latent.dim <- dim(C0)[1]
# initialize
mt <- matrix(NA, ncol=latent.dim, nrow=n)
Ct <- array(NA, dim=c(latent.dim,latent.dim,n))
Rt <- array(NA, dim=c(latent.dim,latent.dim,n))
llh <- -n*log(2*pi)/2
# Filter first observation
G <- Gt(m0, Xt[1,])
FF <- Ft(1, Xt[1,], phi)
### at <- G%*%m0
### REPLACED BY
at <- gt(m0, Xt[1])
### Rt[,,1] <- G%*%C0%*%t(G) + Wt(1,Xt[1,],phi)
### REPLACED BY
Rt[,,1] <- G%*%C0%*%t(G)*1.03
if (is.na(m.start[1]))
{
lambda.ping <- FF%*%at
}
else
{
lambda.ping <- FF%*%m.start[1,]
}
ft <- t(FF)%*%at
Qt <- t(FF)%*%Rt[,,1]%*%FF + lambda.ping*(1-lambda.ping/nt[1])
At <- as.vector(Rt[,,1]%*%FF)/Qt
et <- Yt[1]-ft
mt[1,] <- at + At*et
Ct[,,1] <- Rt[,,1] - At%*%t(At)*as.vector(Qt)
llh <- llh -0.5*(log(Qt)+et^2/Qt)
for(i in 2:n)
{
G <- Gt(mt[i-1,], Xt[i,])
FF <- Ft(i, Xt[i,], phi)
### at <- G%*%mt[i-1,]
### REPLACED BY
at <- gt(mt[i-1,], Xt[i])
### Rt[,,i] <- G%*%Ct[,,i-1]%*%t(G) + Wt(i,Xt[i,],phi)
### REPLACED BY
Rt[,,i] <- G%*%Ct[,,i-1]%*%t(G)*1.03
if (is.na(m.start[i]))
{
lambda.ping <- FF%*%at
}
else
{
lambda.ping <- FF%*%m.start[i,]
}
if (nt[i]==0)
{
mt[i,] <- at
Ct[,,i] <- Rt[,,i-1]
}
else
{
ft <- t(FF)%*%at
Qt <- t(FF)%*%Rt[,,i]%*%FF + lambda.ping*(1-lambda.ping/nt[i])
At <- as.vector(Rt[,,i]%*%FF)/Qt
et <- Yt[i]-ft
mt[i,] <- at + At*et
Ct[,,i] <- Rt[,,i] - At%*%t(At)*as.vector(Qt)
llh <- llh -0.5*(log(Qt)+et^2/Qt)
}
}
ssm$filtered <- list(mt=mt,Ct=Ct, Rt=Rt, llh=llh)
ssm
}
#####################################################
# MODIFIED KALMAN SMOOTHER FOR BINOMIAL DATA
# IDENTITY LINK, NON-LINEAR SYSTEM EQUATION
#####################################################
modified.smoother <- function(ssm, gt, Gt)
{
# Requires
# ssm: object of class state space model
# requires that the filter has been run!
# gt: state function # NEW
# Gt: Hessian of gt # NEW
# Returns an updated ssm object with
# m.tilde
# C.tilde
if(is.na(ssm$filtered$mt[1])) {
cat("You have not run the filter yet - please do that\n")
break
}
# Extract information from filtered object
# Gt <- ssm$Gt # NOT NEEDED
m.tilde <- ssm$filtered$mt
C.tilde <- ssm$filtered$Ct
Rt <- ssm$filtered$Rt
# number of observations
n <- dim(m.tilde)[1]
# Dimension of latent process
latent.dim <- dim(m.tilde)[2]
# dimension of observation
obs.dim <- 1
# Initialize
# First smoothed is last filtered
### G <- Gt(n, Xt[n, ], phi)
### REPLACED BY
G <- Gt(m.tilde[n-1,], Xt[n,])
for (i in (n-1):2) {
Gtp1 <- G
### G <- Gt(i, Xt[i, ], phi)
### REPLACED BY
G <- Gt(m.tilde[i-1,], Xt[i, ])
### Bt <- C.tilde[, , i] %*% t(Gtp1) %*% solve(Rt[, , i + 1])
### REPLACED BY
Bt <- C.tilde[,,i]%*%t(Gtp1)%*%ginverse(Rt[,,i+1])
### m.tilde[i, ] <- m.tilde[i, ]
### + Bt %*% (m.tilde[i + 1, ] - Gtp1 %*% m.tilde[i, ])
### REPLACED BY
m.tilde[i,] <- m.tilde[i,]
+ Bt%*%(m.tilde[i+1,]-gt(m.tilde[i,], Xt[i+1]))
C.tilde[,,i] <- C.tilde[,,i]
+ Bt%*%(C.tilde[,,i+1]-Rt[,,i+1])%*%t(Bt)
}
### THE FOLLOWING IS NEEDED TO SMOOTH THE FIRST STATE
Gtp1 <- G
G <- Gt(ss$m0, Xt[1, ])
Bt <- C.tilde[,,1]%*%t(Gtp1)%*%ginverse(Rt[,,2])
m.tilde[1,] <- m.tilde[1,]
+ Bt%*%(m.tilde[2,]-gt(m.tilde[1,], Xt[2]))
C.tilde[,,1] <- C.tilde[,,1]
+ Bt%*%(C.tilde[,,2]-Rt[,,2])%*%t(Bt)
ssm$smoothed <- list(m.tilde = m.tilde, C.tilde = C.tilde)
ssm
}
# Apply the modified filter
ss <- modified.filter(ss, Gt=Gt, gt=gt)
# Plot filtered awareness
plot(apply(ss$filtered$mt[,c(1,5)],1,sum),
type="l", ylim=c(0,0.8))
points(ss$Yt/ss$nt)
title("Awareness, first application of modified filter")
# Apply the modified smoother
ss <- modified.smoother(ss,Gt=Gt,gt=gt)
# Plot smoothed awareness
plot(apply(ss$smoothed$m.tilde[,c(1,5)],1,sum),
type="l", ylim=c(0,0.8))
points(ss$Yt/ss$nt)
title("Awareness, first application of modified smoother")
# Iterated extended Kalman filter and smoother
m.start <- ss$smoothed$m.tilde
ss$m.start <- m.start
goon <- T
while (goon)
{
ss <- modified.filter(ss, Gt=Gt, gt=gt)
ss <- modified.smoother(ss,Gt=Gt,gt=gt)
if(max(abs((ss$smoothed$m.tilde - m.start)/m.start))<10^-4)
goon <-F
m.start <- ss$smoothed$m.tilde
}
# Plot final filtered and smoothed awareness
plot(apply(ss$filtered$mt[,c(1,5)],1,sum),
type="l", ylim=c(0,0.8))
points(ss$Yt/ss$nt)
title("Awareness, final application of modified filter")
plot(apply(ss$smoothed$m.tilde[,c(1,5)],1,sum),
type="l", ylim=c(0,0.8))
points(ss$Yt/ss$nt)
title("Awareness, final application of modified smoother")
|