Welcome | State space notation | Functions | > Advance example |
Last update Apr 29, 2003 |
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") |