Last update Jul 10, 2004

DESCRIPTION:

USAGE:

      kal.fil.mul.pom(ssm)

REQUIRED ARGUMENTS:

ssm:

Object of class ssm (state space model).
The following attributes must be specified in ssm:

Yt:

n × k matrix containing the multinomial observations arranged row-wise.

nt:

1 × n vector containing the number of trials at each time-point.
Should be nt <- apply(Yt,1,sum).

Ft:

Function returning a p × (k-1) design matrix.

Gt:

Function returning a p × p evolution transfer matrix.

Wt:

Function returning a p × p evolution variance matrix.

m0:

p × 1 vector (prior mean).

C0:

p × p matrix (prior variance).

OPTIONAL ARGUMENTS:

Xt:

Matrix containing the covariates needed in Ft, Gt and/or Wt. This matrix must have n rows, such that row t contains the covariates needed for generating the system matrices at time t.

psi:

Parameter vector must be supplied if it is needed in the calculation of the system matrices.

fam:

Tacitly assumed to have value multinomial.

link:

Tacitly assumed to have value pom.

m.start:

Used to specify the Taylor expansion points. These are not meant to be set by the user. Applying the Kalman smoother on a ssm object that has been filtered will set the m.start attribute to the smoothed states. This is needed in the iterated extended Kalman filter and smoother.

VALUE:

mt:

Posterior means m_t given the observations up to and including time t. These are arranged row-wise in a matrix.

Ct:

Posterior variances C_t given the observations up to and including time t. These are arranged in a 3-dimensional array such that C_t is placed in filtered$Ct[,,t].

llh:

Log-likelihood of of the approximating Gaussian state space model.

Rt:

Prior variances for the states arranged in 3-dimensional array such that R_t is placed in filtered$Rt[,,t]. The reason for returning these is that they are needed in the Kalman smoother.

SIDE EFFECTS:

DETAILS:

REFERENCES:

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")