source('clarkFunctions2026.r')

A state-space model with covariates

Now suppose that the process is driven by covariates that change over time,

\[ \begin{aligned} x_t &\sim N(x_{t-1} + \mathbf{z}'_{t-1} \boldsymbol{\beta}, \sigma^2) \\ y_t &\sim Poi( e^{x_t} ) \end{aligned} \] Here are simulated data:

nt    <- 300                                     # no. time intervals
time  <- c(1:nt)                                 # sample times
wm    <- 30:60                                   # a sequence of missing values
sig   <- .01                                     # process err variance
x     <- rep(5,nt)
q     <- 3
z     <- matrix( rnorm(q*(nt - 1), .1), nt-1, q) # design matrix
z[,1] <- 1
beta  <- matrix(runif(q, -.01,.01))              # coefficient vector
zb    <- z%*%beta

for(t in 2:nt)x[t] <- x[t-1] + rnorm(1, zb[t-1], sqrt(sig))

y     <- rpois( nt,exp(x) )                      # observations
y[wm] <- NA                                      # missing

plot(time, y, col='grey', ylim = range( y, na.rm=T ) )
lines(time, exp(x) )
*Observed (dots) and latent states (line) for the Poisson model.*

Observed (dots) and latent states (line) for the Poisson model.

Here is model fitting:

cat( "model {

 mu[1] ~ dnorm( log(y[1]),1)
 xg[1] <- mu[1]

 for(t in 2:nt) {
   mu[t] <- xg[t-1] + b1 + b2*z[t-1,2] + b3*z[t-1,3]
   xg[t] ~ dnorm(mu[t],procErr)
   y[t] ~ dpois( exp(xg[t]) )
 }

 procErr ~ dunif(0,10)
 sigma <- 1/procErr
 b1 ~ dunif(-1,1)
 b2 ~ dunif(-1,1)
 b3 ~ dunif(-1,1)
}  ", fill=T, file="stateSpacePoisson.txt" )

Here is the analysis:

library(rjags)
ssData   <- list(y = y, z = z, nt = nt )
parNames <- c('b1','b2','b3', 'sigma')
ssFit    <- jags.model(data=ssData, file="stateSpacePoisson.txt")
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 268
##    Unobserved stochastic nodes: 335
##    Total graph size: 2704
## 
## Initializing model
update(ssFit)

ssFit <- coda.samples(ssFit, variable.names = parNames, n.iter=5000)
par(mar=c(3,3,1,1))
plot( as.mcmc(ssFit) )
*MCMC chains.*

MCMC chains.

Here is a comparison of chains and true values:

chainPlot( as.matrix(ssFit), trueValues= c(beta, sqrt(sig)) )
*Posterior densities with 95% CI (dashed) and true values (red).*

Posterior densities with 95% CI (dashed) and true values (red). 

chain2density( as.matrix(ssFit), trueValues= c(beta, sqrt(sig)) )
*Posterior densities with 95% CI (dashed) and true values (red).*

Posterior densities with 95% CI (dashed) and true values (red).