# Chapter 4: Fitting state-space models

## Kalman Filter, (Negative log) Likelihood, and MLE code

for NDLM with 2 component state and observation vectors (See Table 4.1 and Section 4.4.1)

```# ---------- Part 1:  Simulating the NDLM -----------------------
#   n[t] = A x n[t-1] + "Q"
#   y[t] = B x n[t] + "R"
set.seed(22)
library(mvtnorm)

# State process matrix
A <- matrix(data=c(0.0, 1.5, 0.3, 0.55),nrow=2,ncol=2,byrow=TRUE)
# State variance
sigma.p <- 0.5

# Obs'n model matrix
B <- diag(x=c(0.8,1.1),nrow=2)
# Obs'n variance
sigma.o <- 2.0

# Initial state values
n0 <- cbind(c(30,50))
num.times <- 10
n.vec <- y.vec <- matrix(data=0,nrow=2,ncol=num.times+1,
n.vec[,1] <- n0

#-- Simulate states and observations
for(i in 1:num.times) {
n.vec[,i+1] <- A %*% n.vec[,i] + rnorm(n=2,mean=0,sd=sigma.p)
y.vec[,i+1] <- B %*% n.vec[,i+1] + rnorm(n=2,mean=0,sd=sigma.o)
}
print(n.vec)
print(y.vec)

# ---------- Part 2:  Kalman Filter Function -----------------------
kf <- function(obs,A,B,Q,R,n0) {
#-- Kalman filter estimates of state vectors

# Output:
#  Kalman filter estimates of state vector for static SSM
#  Also returns matrix of 1 step ahead predicted states and corresponding
#  array of 1 step ahead predicted covariance matrices

# Input:
# obs is a p x [T+1] matrix
# A is the state transition matrix
# Q is the state covariance matrix, n[t] = A x n[t-1] + "Q"
# B is the observation "linkage" matrix
# R is the obs'n covariance matrix, y[t] = B x n[t] + "R"
# n0 is the initial state vector

dim.state <- dim(Q)[1]
filter.n <- matrix(data=0,nrow=dim.state,ncol=num.times+1)
pred.n <- update.n <- matrix(0,nrow=dim.state,ncol=1)
cov.pred.n <- cov.update.n <- diag(x=rep(0,dim.state),nrow=dim.state)

#-- the following are used for likelihood evaluation
pred.n.matrix <- matrix(data=0,nrow=dim.state,ncol=num.times+1)
cov.pred.array <- array(data=0,dim=c(dim.state,dim.state,num.times+1))
pred.n.matrix[,1] <- n0

# the iterations
update.n <- n0
filter.n[,1] <- update.n
identity.mat <- diag(x=1,nrow=dim.state)
for(i in 1:num.times) {
pred.n <- A %*% update.n
cov.pred.n <- A %*% cov.update.n %*% t(A) + Q

Kalman.gain <- cov.pred.n %*% t(B) %*% solve(B %*% cov.pred.n %*% t(B) + R)

update.n <- pred.n + Kalman.gain %*% (y.vec[,i+1,drop=FALSE]-B%*%pred.n)
cov.update.n <- (diag(x=1,nrow=2)-Kalman.gain %*% B) %*% cov.pred.n

filter.n[,i+1] <- update.n

pred.n.matrix[,i+1] <- pred.n
cov.pred.array[,,i+1] <- cov.pred.n
}
out <- list(filter.n=filter.n,pred.n.matrix=pred.n.matrix,
cov.pred.array=cov.pred.array)
return(out)
}

# ---------- Part 3: Applying KF with true parameter values  ------#
Q <- diag(sigma.p^2,nrow=2)
R <- diag(sigma.o^2,nrow=2)
out <- kf(obs=y.vec,A=A,B=B,Q=Q,R=R,n0=n0)
print(out\$filter.n, digits=3)

# Juvenile Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[1,],n.vec[1,],out\$filter.n[1,]))
dimnames(temp) <- list(c("y-juv","n-juv","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

# adult Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[2,],n.vec[2,],out\$filter.n[2,]))
print(temp)

# ---------- Part 3: Negative Log-Likelihood Function  ------
#  For state matrix A, with independent state components.
#  Need to have mvtnorm package installed
library(package=mvtnorm)
like.fn.2state.SSM <- function(theta,obs,B,Q,R,n0) {
num.times  <- dim(obs)[2]-1  # time starts at 0
num.states <- dim(Q)[1]
A <- matrix(data=theta,nrow=num.states,ncol=num.states,byrow=FALSE)
out <- kf(obs=obs,A=A,B=B,Q=Q,R=R,n0=n0)
pred.n.matrix  <- out\$pred.n.matrix
cov.pred.array <- out\$cov.pred.array

log.l <- 0
for(i in 1:num.times) {
cond.mean.vector <- B %*% pred.n.matrix[,i]
cond.cov.matrix  <- B %*% cov.pred.array[,,i] %*% t(B) + R
log.l.comp <- dmvnorm(x=obs[,i+1], mean=cond.mean.vector,
sigma=cond.cov.matrix, log=TRUE)
log.l      <- log.l + log.l.comp
}
out <- -log.l
return(out)
}

#-- testing log likelihood function
Q <- diag(sigma.p^2,nrow=2)
R <- diag(sigma.o^2,nrow=2)
like.fn.2state.SSM(theta=as.vector(A),obs=y.vec,B=B,Q=Q,R=R,n0=n0)

A.alt <- jitter(A)
like.fn.2state.SSM(theta=as.vector(A.alt),obs=y.vec,B=B,Q=Q,R=R,n0=n0)

# ---------- Part 4: Calculating MLE of State Transition Matrix  ------
A.alt <- jitter(A)
mle <- optim(par=as.vector(A.alt),fn=like.fn.2state.SSM,
obs=y.vec,B=B,Q=Q,R=R,n0=n0)

A.est <- matrix(data=mle\$par,nrow=dim(Q)[1],ncol=dim(Q)[1],byrow=FALSE)
cat("MLE of State Transition Matrix \n"); print(A.est)
cat("True value of S.T.M. \n"); print(A)

# checks on optimization
like.fn.2state.SSM(theta=as.vector(A.est),obs=y.vec,B=B,Q=Q,R=R,n0=n0)
like.fn.2state.SSM(theta=as.vector(A),obs=y.vec,B=B,Q=Q,R=R,n0=n0)

# Applying Kalman Filter with mle for A
out <- kf(obs=y.vec,A=A.est,B=B,Q=Q,R=R,n0=n0)

# Juvenile Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[1,],n.vec[1,],out\$filter.n[1,]))
dimnames(temp) <- list(c("y-juv","n-juv","n.kf"),0:(dim(y.vec)[2]-1))
print(temp)

# adult Obs'ns, True Values, and Filtered Values
temp <- round(rbind(y.vec[2,],n.vec[2,],out\$filter.n[2,]))
print(temp)

```

### Section 4.5.2.1 (WinBUGS code)

``````
#(i) Model Statement
model {
#Prior distribution for parameters
beta1  ~ dnorm(0,0.01)
beta0  ~ dnorm(0,0.01)
sigma  ~ dunif(0.01,10)
tau <- 1/(sigma*sigma) #the precision
#Likelihood defined
for(i in 1:n) {
mu[i] <- beta0+beta1*x[i]
y[i]  ~ dnorm(mu[i],tau)
}
}

#(ii) Data Statement
list(n=15,
y=c( 29.4, 33.8, 24.8, 26.1, 31.6, 29.4, 14.6, 12.6, 29.3, 27.5,
28.9, 28.3 ,24.7, 28.5, 14.9),
x=c(8, 10, 7, 7, 9, 9, 3, 3, 8, 7, 8, 8, 7, 8, 3))
#(iii) Initial Values Statement
list(beta0=0.5, beta1=1.0, sigma=3.0)
``````

### Section 4.5.3 salmon SSM in WinBUGS

Save this block of WinBUGS code locally as C:/PopDynBook/BUGSCode/Ricker-Poisson-logN.txt)

``````
model {
#WinBUGS model code for Ricker model Poisson-LogN
# assuming known obs'n noise CV

# Priors
alpha    ~ dunif(1,2.5)
beta     ~ dunif(0.00001,0.001)
n.init   ~ dunif(50,500)
sigma.sq <- log(CV.obs*CV.obs+1)
tau.obs  <- 1/sigma.sq

mu[1]    <- alpha*n.init*exp(-beta*n.init)
state[1] ~ dpois(mu[1])
temp[1]  <- log(state[1])-sigma.sq/2
obs[1]   ~ dlnorm(temp[1],tau.obs)

for(i in 2:T) {
mu[i]    <- alpha*state[i-1]*exp(-beta*state[i-1])
state[i] ~  dpois(mu[i])
temp[i]  <- log(state[i])-sigma.sq/2
obs[i]   ~  dlnorm(temp[i],tau.obs)
}
}
}
``````

Load the following code into R

``````
library(R2WinBUGS)
input.data <- list(T=num.times, obs=estimates,CV.obs=CV.obs)
init.value.generator <- function(num.starts=1,n,alpha.params=c(1,2),
beta.params=c(0.00001,0.001),n.init.params=c(10,2000)) {
init.values <- NULL
for(i in 1:num.starts) {
alpha<- runif(1,alpha.params[1],alpha.params[2])
beta <- runif(1,beta.params[1],beta.params[2])
n.init <- runif(1,n.init.params[1],n.init.params[2])
state <- numeric(n)
state[1] <- rpois(1,alpha*n.init*exp(-beta*n.init))
for(j in 2:n)
state[j] <-rpois(1,alpha*state[j-1]*exp(-beta*state[j-1]))
init.values[[i]] <-
list(alpha=alpha, beta =beta,n.init=n.init,state=state)
}
return(init.values)
}

init.values <- init.value.generator(num.starts=3,n=num.times,
n.init.params=c(0.5*estimates[1],1.2*estimates[1]))

params <- c("alpha","beta","n.init","state")

out.ssm <- bugs(data=input.data, inits=init.values,
parameters.to.save=params, model.file=
"C:/PopDynBook/BUGSCode/Ricker-Poisson-logN.txt",
n.chains=3, n.iter=50000, n.burnin=10000,n.thin=10,debug=TRUE)
``````

### Section 4.5.4 WinBUGS model statement for the BRS SSM

``````
model {
#WinBUGS model code for BRS model (Bin-Bin-Bin)-LogN assuming known obs'n noise CV

# Priors
#Survival parameters
phi1     ~ dunif(0.05,0.95)
phi2     ~ dunif(0.05,0.95)

#Growth
Pi       ~ dunif(0.05,0.95)
not.Pi   <- 1-Pi

#Birth
lambda   ~ dunif(0.05,0.95)

n1.init  ~ dunif(10,500)
n2.init  ~ dunif(10,500)

sigma.sq <- log(CV.obs*CV.obs+1)
tau.obs  <- 1/sigma.sq

#Year 1 calculations
u1s[1]   ~ dbin(phi1,n1.init)
u2s[1]   ~ dbin(phi2,n2.init)

u1r[1]   ~ dbin(not.Pi,u1s[1])
n2[1]    <- u2s[1] + u1s[1]-u1r[1]

births[1] ~ dbin(lambda,n2[1])
n1[1]    <- u1r[1] + births[1]

temp1[1]    <- log(n1[1]) - sigma.sq/2
obs.imm[1]   ~ dlnorm(temp1[1],tau.obs)
temp2[1]    <- log(n2[1]) - sigma.sq/2
obs.mat[1]  ~ dlnorm(temp2[1],tau.obs)

#Iterate over remaining years
for(i in 2:T) {
u1s[i] ~ dbin(phi1,n1[i-1])
u2s[i] ~ dbin(phi2,n2[i-1])

u1r[i] ~ dbin(not.Pi,u1s[i])
n2[i] <- u2s[i] + u1s[i]-u1r[i]

births[i] ~ dbin(lambda,n2[i])
n1[i] <- u1r[i] + births[i]

temp1[i]    <- log(n1[i]) - sigma.sq/2
obs.imm[i]   ~ dlnorm(temp1[i],tau.obs)
temp2[i]    <- log(n2[i]) - sigma.sq/2
obs.mat[i]  ~ dlnorm(temp2[i],tau.obs)
}
}
``````

R code for generating initial values to fit BRS SSM with WinBUGS.

``````
init.value.generator <- function(num.starts=1,n,
phi1.params=c(0.1,0.9),  phi2.params=c(0.1,0.9),
Pi.params = c(0.1,0.9), lambda.params=c(0.1,0.9),
n1.init.params=c(10,500), n2.init.params=c(10,500)) {
init.values <- NULL
for(i in 1:num.starts) {
phi1   <- runif(1,phi1.params[1],phi1.params[2])
phi2   <- runif(1,phi2.params[1],phi1.params[2])
Pi     <- runif(1,Pi.params[1],Pi.params[2])
lambda <- runif(1,lambda.params[1],lambda.params[2])

n1.init <- round(runif(1,n1.init.params[1],n1.init.params[2]))
n2.init <- round(runif(1,n2.init.params[1],n2.init.params[2]))
n1 <- n2 <- u1s <- u2s <- u1r <- births <- numeric(n)

u1s[1]    <- rbinom(1,size= n1.init,prob=phi1)
u2s[1]    <- rbinom(1,size= n2.init,prob=phi2)

u1r[1]    <-rbinom(1,size=u1s[1],prob=1-Pi)
n2[1]     <- u2s[1] + u1s[1]-u1r[1]
births[1] <- rbinom(1,size=n2[1],prob=lambda)
n1[1]     <- u1r[1]+ births[1]

for(j in 2:n) {
u1s[j] <- rbinom(1,n1[j-1],phi1)
u2s[j] <- rbinom(1,n2[j-1],phi2)

u1r[j] <- rbinom(1,u1s[j],1-Pi)
n2[j]  <- u2s[j] + (u1s[j]-u1r[j])

births[j] <- rbinom(1,n2[j],lambda)
n1[j] <- u1r[j] + births[j]
}

init.values[[i]] <-
list(phi1=phi1, phi2 =phi2,Pi=Pi,lambda=lambda,
n1.init=n1.init,n2.init=n2.init,
u1s=u1s,u2s=u2s,u1r=u1r,births=births)
}
return(init.values)
}

init.values <- init.value.generator(num.starts=3,n=num.yrs,
phi1.params=c(0.45,0.9),  phi2.params=c(0.55,0.9),
Pi.params = c(0.5,0.9), lambda.params=c(0.7,0.85),
n1.init.params=c(40,100), n2.init.params=c(40,100))
``````

### Section 4.5.5.2 Sequential importance sampling

R code for univariate coho salmon SSM (revised January 2016)

``````
# Chapter SIS example complete
# Section 4.5.5.2 Sequential importance sampling
# R code for using SIS to fit univariate coho salmon SSM

ricker.sis <- function(N=1000,params,obs,CV.obs=0.3,smth.param=0.95) {
r.a <- params[,1]; r.b <- params[,2]; n0 <- params[,3]
obs.sd <- sqrt(log(CV.obs^2+1))
num.times <- length(obs)
states <- matrix(NA,nrow=N,ncol=num.times)
uni.tally <- 1:N

states[,1] <- rpois(n=N,lambda=r.a*n0*exp(-r.b*n0))
wts <- dlnorm(x=obs,meanlog=log(states[,1])-obs.sd^2/2,obs.sd)
for(i in 2:num.times) {
states[,i] <- rpois(n=N,lambda=r.a*states[,i-1]*exp(-r.b*states[,i-1]))
wts <- wts*dlnorm(x=obs[i],meanlog=log(states[,i])-obs.sd^2/2,obs.sd)
new.set <- sample(1:N,N,replace=TRUE,prob=wts)
uni.tally <- uni.tally[new.set]
wts <- 1 #re-initialize wts
states <- states[new.set,]
params <- params[new.set,]

#kernel smoothing of params
mean.param <- apply(params,2,mean)
cov.param  <- cov(params)
smth.params <- smth.param*params + (1-smth.param)*mvrnorm(n=N,mu=cbind(mean.param),Sigma=cov.param)
params <- smth.params
r.a <- params[,1]; r.b <- params[,2]; n0 <- params[,3]
}
out <- list(states=states,params=params,tally=uni.tally)
return(out)
}

# simulate coho data
set.seed(2345)
num.times <- 20
Ricker.a <- 1.5
Ricker.b <- 0.0003
equilibrium <-  log(Ricker.a)/Ricker.b
cat("equilibrium=", round(equilibrium),"\n")

initial.n <- round(0.1*equilibrium)
juveniles <- estimates <- numeric(num.times)
juveniles[1] <- rpois(1,lambda = Ricker.a*initial.n*exp(-Ricker.b*initial.n))
CV.obs <- .30
obs.sd <- sqrt(log(CV.obs^2+1))
estimates[1] <- round(rlnorm(1,meanlog=log(juveniles[1]),sdlog=obs.sd))
for(i in 2:num.times) {
juveniles[i] <- rpois(1,lambda=Ricker.a*juveniles[i-1]*exp(-Ricker.b*juveniles[i-1]))
estimates[i] <-  round(rlnorm(1,meanlog=log(juveniles[i])-obs.sd^2/2,sdlog=obs.sd))
}

# Look at state and obs'n time series
my.ylim <- range(c(juveniles,estimates) )
plot(as.ts(juveniles),xlab="Time",ylab="",pch="s",col="blue",type="b",
main="Simulated Coho Salmon Population and Estimates",ylim=my.ylim)
lines(as.ts(estimates),pch="o",col="red",type="b")

# set-up for SIS
num.particles <- 200000
set.seed(814)
init.params <-  matrix(NA,nrow=num.particles,ncol=3,dimnames=list(NULL,c("alpha","beta","n0")))
init.params[,1] <- runif(num.particles,min=1,max=2.5)
init.params[,2] <- runif(num.particles,min= .00001,max=0.001 )
init.params[,3] <- runif(num.particles,min=50,max=500)
dimnames(init.params) <- list(NULL,c("alpha","beta","n0"))

par(mfrow=c(2,2))
for(i in 1:3) hist(init.params[,i],main=dimnames(init.params)[[2]][i])
par(mfrow=c(1,1))

# run SIS
library(MASS)

temp <- ricker.sis(N=num.particles,params=init.params,obs=estimates,CV.obs=CV.obs,smth.param=0.95)

# Examine results
cat("particle depletion=",1-length(unique(temp\$tally))/num.particles,"\n")

post.states <- apply(temp\$states,2,mean)

par(mfrow=c(2,2),oma=c(0,0,3,0))
true.par <- c(Ricker.a,Ricker.b,initial.n)
for(i in 1:3) {
my.xlim <- range(c(init.params[,i],temp\$params[,i]))
my.ylim <- range(c(density(init.params[,i])\$y, density(temp\$params[,i])\$y))
plot(density(init.params[,i]),xlim=my.xlim,ylim=my.ylim,
main=dimnames(init.params)[[2]][i],xlab="")
lines(density(temp\$params[,i]),col=2,lty=2)
text(c(true.par[i],mean(temp\$params[,i])),c(0,0),"X",col=1:2)
}
my.ylim <- range(c(juveniles,estimates,post.states))
plot(1:num.times,juveniles,xlab="Year",ylab="",ylim=my.ylim,type="l")
lines(1:num.times,estimates,lty=2,col=2)
lines(1:num.times,post.states,lty=3,col=3,lwd=2)
legend(7,0.45*max(my.ylim),legend=c("State","Observation","Smooth"),lty=1:3,col=1:3)
par(mfrow=c(1,1))
print(rbind(true.par,apply(temp\$params,2,mean)))``````