Welcome | State space notation | Functions | > kal.fil.mul.pom |
Last update Apr 29, 2003 |
DESCRIPTION:Extended Kalman filter for multinomial data assuming pom link and linear Gaussian latent process. It is a dynamic extension of the proportional odds model of McCullagh (1980). The specific model is
USAGE:kal.fil.mul.pom(ssm)
REQUIRED ARGUMENTS:
OPTIONAL ARGUMENTS:The following attributes of the ssm object might be needed:
VALUE:Returns an object of class ssm with the same attributes as the object in the call, but with attribute filtered updated with:
SIDE EFFECTS:None
DETAILS:The filter works by sequentially approximating the multinomial observation density with a Gaussian density based on a Taylor expansion around the values supplied in the attribute m.start of the ssm object. If m.start is not supplied, the Taylor expansion will be around the one-step forecast mean of the observation.
REFERENCES:Durbin & Koopman, (2002), Time Series Analysis by State Space Models, Oxford Statistical Science Series.Fahrmeir & Tutz, (1994), Multivariate Statistical Modelling Based on Generalized Linear Models, Springer Series in Statistics. Klein, (2003), State Space Models for Exponential Family Data, Ph.D. Thesis, Department of Statistics, University of Southern Denmark. McCullagh, (1980), Regression models for ordinal data (with discussion) Journal of the Royal Statistical Society, Series B, Vol. 42, 109-142. EXAMPLES:# Application of the iterated extended Kalman filter and smoother # to fit a proportional odds model # (example from Chapter 5) # Observations Yt <- matrix(c(0, 0, 1, 7, 8, 8,19, 8, 1, 6, 9,12,11, 7, 6, 1, 0, 0, 1, 1, 6, 8,23, 7, 5, 1, 0, 0, 0, 0, 1, 3, 7,14,16,11), ncol=9, byrow=T) # Number of trials nt <- rep(52,4) # Design matrix Ft <- function(i,xt,phi) { result <- rep(0,11) if (((i-1)%%4+1)==1) result <- t(cbind(diag(8),0,0,0)) if (((i-1)%%4+1)==2) result <- t(cbind(diag(8),-1,0,0)) if (((i-1)%%4+1)==3) result <- t(cbind(diag(8),0,-1,0)) if (((i-1)%%4+1)==4) result <- t(cbind(diag(8),0,0,-1)) result } # Evolution transfer matrix is the identity Gt <- function(i, xt,phi) diag(11) # Zero evolution noise Wt <- function(i,xt,phi) diag(rep(0,11)) # Specify the state space model object cheese.ss <- ssm(Ft=Ft, Gt=Gt, Wt=Wt, m0=c(-5.5,-4.4,-3.3,-2,-1,0,1.5,3,-3.4,-1.7,1.6)/2, C0=diag(rep(1000,11)), Yt=Yt, nt=nt, fam="multinomial", link="pom") # Apply the iterated extended Kalman filter and smoother cheese.ss <- ieks(cheese.ss, eps=10^-9) # Estimates and standard errors cbind(cheese.ss$smoothed$m.tilde[4,], sqrt(diag(cheese.ss$smoothed$C.tilde[,,4]))) ######################################################## # Specify state space model ss <- ssm(Ft = function(i,x,phi) {t(cbind(diag(5),-1))}, Gt = function(i,x,phi) {diag(6)}, Wt = function(i,x,phi) {diag(c(rep(phi[1],5),phi[2]))}, phi = c(0.001,0.1), m0 = c(-2,-1,0,1,2,0), C0 = diag(1,6), nt = floor(runif(50,min=20,max=50)), fam = "multinomial", link = "pom") # Simulate observations ss <- simulate.ssm(ss, n=50) # Plot simulated latent process plot(ss$simulated$theta[,6], type="l", ylim=range(ss$simulated$theta)) for (i in 1:5) lines(ss$simulated$theta[,i], lty=2) title("Simulated latent process") # Plot cumulated observations tmp <- t(apply(ss$Yt,1,cumsum))/ss$nt plot(tmp[,6], ylim=c(0,1), type="l") for(i in 1:5) lines(tmp[,i]) title("Cumulated observations") # Apply the extended Kalman filter ss <- kal.fil.mul.pom(ss) # Plot the filtered states. plot(ss$filtered$mt[,6], xlab="time", ylim=range(ss$filtered$mt), type="l") for(i in 2:6) lines(ss$filtered$mt[,i]) title("Filtered latent process") # Apply the iterated extended Kalman filter and smoother ss <- ieks(ss) # Smoothed trajectory of beta_t upper <- ss$smoothed$m.tilde[,6] + sqrt(ss$smoothed$C.tilde[6,6,]) lower <- ss$smoothed$m.tilde[,6] - sqrt(ss$smoothed$C.tilde[6,6,]) plot(ss$smoothed$m.tilde[,6], type="l", lty=1, ylim=range(upper,lower)) lines(ss$simulated$theta[,6], lty=2) lines(upper, lty=3) lines(lower, lty=3) title("Smoothed trajectory of beta_t") # Plot the smoothed states. plot(ss$smoothed$m.tilde[,6], xlab="time", ylim=range(ss$smoothed$m.tilde), type="l") for(i in 1:5) lines(ss$smoothed$m.tilde[,i], lty=2) title("Smoothed latent process") |