Modelling Population Dynamics: Chapter 9

Integrated Methods

Section 9.3.1 Data, models and integrated modelling

Kalman filter with simulated data

y <- seq(10,150,10)
n <- length(y)
Z <- matrix(c(0,1),1,2)
H <- 1
a <- matrix(c(25,50),2,1)
P <- diag(c(100,100))
p <- 1
phi1 <- 0.6
phiA <- 0.8

logL <- 0 
for (t in 1:(n-1)){ 
  v <- y[t]-Z %*% a
  F <- Z %*% P%*% t(Z) + H
  M <- P %*% t(Z) 
  print(c(v,F)) 
  logL <- logL - log(det(F)) - t(v) %*% solve(F) %*% v
  att <- a + M %*% solve(F) %*% v  
  Ptt <- P - M %*% solve(F) %*%t (M)  
  Ptt <- (Ptt+t(Ptt))/2
  T <- matrix(c(0, p*phi1, phiA, phiA), 2, 2, byrow=TRUE)
  Q <- diag(c((T * (1-T)) %*% att)) 
  Q[1,1] <- T[1,] %*% att 
  a <- T %*% att 
  P <- T %*% Ptt %*% t(T) + Q 
  P <- (P+t(P))/2 
  }
 v <- y[n] - Z %*% a
 F <- Z %*% P %*% t(Z) + H 
 logL <- logL - log(det(F)) - t(v) %*% solve(F) %*% v 
 c(-0.5*logL) # minus log-likelihood