# R code for implementing the extended Kalman filter

The code below implements the discrete-time extended Kalman filter (EKF) in R.

For numerical stability and precision the implemented EKF uses a Singular Value Decomposition (SVD) based square root filter. For a description of this SVD-based square root filter see Appendix B of Petris and colleagues’ 2009 book Dynamic linear models with R.

See this blog post for implementing the Unscented Kalman Filter (UKF) in R, this post for implementing the Gauss-Hermite Kalman Filter (GHKF), and this post for implementing the Ensemble Kalman Filter (EnKF).

See this page for an overview of all of Stefanâ€™s R code blog posts.

R code (last update: 2016/02/15)

```library(dlm)
library(numDeriv)

#The R functions for the extended Kalman filter (EKF) and its smoother
source("http://www.datall-analyse.nl/EKF.R")
#Take a look at the functions dlmExtFilter (=the EKF) and dlmExtSmooth (=the
#RTS smoother), and notice that at the beginning of the scripts you will find
#a description of the functions' arguments.
dlmExtFilter
dlmExtSmooth
#You may save these functions to your computer using the save function in R.

##EXAMPLE 1

#In this example we will apply the Extended Kalman Filter (EKF) to the
#Lotka-Volterra system of nonlinear differential equations. The Lotka-Volterra
#equations describe how prey and predators interact.
#en.wikipedia.org/wiki/Lotka-Volterra_equations

##TRUE STATES

#Generate the amount of prey and predators based on the Lotka-Volterra equations.
#Data are generated using the forward Euler method.
dt <- 1/5000 #time step for Euler integration
tt <- 7 #upper bound of time window
st <- seq(0, tt, by=dt) #lower time bounds of the integration intervals
ns <- length(st) #number of Euler integrations
x <- matrix(0, ncol=2, nrow=ns) #specify matrix for states
x[1,] <- c(400, 200) #respective amounts of prey and predators at t=0
colnames(x) <- c("Prey", "Predators")

#parameters in the Lotka-Volterra model
alpha <- 1
beta <- 1/300
delta <- 1/200
gamma <- 1

#simulate true states
for (i in 2:ns) {
#prey population
x[i,1] <- x[i-1,1] + (alpha*x[i-1,1] - beta*x[i-1,1]*x[i-1,2])*dt
#predator population
x[i,2] <- x[i-1,2] + (delta*x[i-1,1]*x[i-1,2] - gamma*x[i-1,2])*dt
}

#phase-space plot of simulated predator-prey system
plot(x[,1], x[,2], type="l", xlab="Prey", ylab="Predators")

##MEASUREMENTS

#Take measurements with a sample time of .01
dT <- .01 #sample time for measurements
#you may change the value of dT and see how it influences
#the behavior of the extended Kalman filter
nm <- tt/dT #total number of measurements
mt <- seq(dT, tt, dT) #measurement times

#standard deviations for the measurement noise
sigmaPrey <- 7 #prey
sigmaPred <- 10 #predators

#measurements at specified measurement times
yprey <- sapply(1:nm, function(i) x[i*((ns-1)/nm) + 1, 1] + rnorm(1, 0, sigmaPrey))
ypred <- sapply(1:nm, function(i) x[i*((ns-1)/nm) + 1, 2] + rnorm(1, 0, sigmaPred))

#store measurement data
dataEx1 <- cbind(yprey, ypred)

#plot the generated measurements
plot(x[,1], x[,2], type="l", xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred))
points(yprey, ypred, cex=.5, col="red")

#plot of time against measurements
par(mfrow=c(2, 1))
plot(st, x[,1], type="l", xlab="Time", ylab="Prey",
ylim=range(yprey))
lines(mt, yprey, col="darkgreen")

plot(st, x[,2], type="l", xlab="Time", ylab="Predators",
ylim=range(ypred))
lines(mt, ypred, col="brown")
par(mfrow=c(1, 1))

##EXTENDED KALMAN FILTER (EKF)

#Dynamic model:
#specifying 2 states, namely [amount of prey, amount of predators].
#We will first use the dlm function from the dlm package for specifying the
#dynamic model. The dlm function checks for us if the dimensions
#of our specified matrices are correct, and if the specified covariance matrices
#are positive semidefinite.
ex1 <- dlm(m0=c(400, 200), #initial state estimates for prey and predators
#error covariances of the initial state estimates:
#this matrix reflects the uncertainty in our initial state estimates
C0=diag(rep(100,2)),
#observation matrix:
#we will not use this FF matrix in the extended Kalman filter,
#so we set the values in this matrix to zero
FF=matrix(0, nrow=2, ncol=2),
#measurement noise
V=diag(c(sigmaPrey^2, sigmaPred^2)),
#state transition matrix:
#we will not use this GG matrix in the extended Kalman filter,
#so we also set the values in this matrix to zero
GG=matrix(0, nrow=2, ncol=2),
#process noise
W=diag(rep(0,2)))

#For the EKF we will use a list-object (instead of the dlm-object above),
#and remove the FF and GG matrix (which we do not need anyway for the EKF).
#Note that the specified initial state estimates (at t=0) below for prey and
#predators deviate from the values that were used for generating the data
#(i.e., 400 for prey and 200 for predators).
#You may change these initial state estimates too and see how they
#influence the behavior of the extended Kalman filter.
ex1 <- list(m0=c(390, 220), #initial state estimates
#error covariances of the initial state estimates:
#this matrix reflects the uncertainty in our initial state estimates
#you may change the values in this matrix and see how they influence
#the behavior of the Kalman filter
C0=diag(rep(100,2)),
#measurement noise
V=diag(c(sigmaPrey^2, sigmaPred^2)),
#process noise
W=diag(rep(0,2)))
#Note that all covariances in the process noise matrix (W) are set to zero.
#This makes sense since the change in the amount of prey and predators at
#each time step is fully explained by our Lotka-Volterra model describing the
#prey-predator interaction.

#Specify the state transition function:
#note that we will use as state functions the difference equations
#given by Euler's forward method. These difference equations will yield valid
#estimates for the amounts of prey and predators at each time step as long
#as the specified value for dT above is relatively small.
#WARNING: always use arguments x and k when specifying the GGfunction
GGfunction <- function (x, k){
prey <- x[1]; pred <- x[2]
c(prey + (alpha*prey - beta*prey*pred)*dT,
pred + (delta*prey*pred - gamma*pred)*dT)}

#Specify the observation/measurement function
#WARNING: always use arguments x and k when specifying the FFfunction
FFfunction <- function (x, k){
prey <- x[1]; pred <- x[2]
c(prey, pred)}

##Compute the filtered (a posteriori) state estimates
ext1 <- dlmExtFilter(y=dataEx1, mod=ex1,
GGfunction=GGfunction, FFfunction=FFfunction)

#Instead of relying on a numerical method for approximating the Jacobians
#in the EKF, it is also possible to calculate the Jacobians by hand and
#subsequently use these in the EKF.
#WARNING: always use arguments x and k when specifying the GGjacobian
GGjacobian <- function (x, k){
prey <- x[1]; pred <- x[2]
c(1 + alpha*dT - beta*pred*dT, -beta*prey*dT,
delta*pred*dT, 1 + delta*prey*dT - gamma*dT)}

#WARNING: always use arguments x and k when specifying the FFjacobian
FFjacobian <- function (x, k){
prey <- x[1]; pred <- x[2]
c(1, 0,
0, 1)}

#Use these latter Jacobians in the EKF
ext1 <- dlmExtFilter(y=dataEx1, mod=ex1,
GGfunction=GGfunction, FFfunction=FFfunction,
GGjacobian=GGjacobian, FFjacobian=FFjacobian)

#plot the filtered state estimates
plot(x[,1], x[,2], type="l", lwd=2, xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred), col=gray(level=.5))
points(yprey, ypred, col="red", cex=.5)
lines(ext1\$m[,1], ext1\$m[,2], lty=2, col="blue", lwd=2)
legend("topright", pch=c(NA, 1, NA), lty=c(1, NA, 2), lwd=c(2, NA, 2),
col=c(gray(level=.5), "red", "blue"),
legend=c("true state", "measured", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#compute the confidence intervals for the filtered state estimates
varcovFilteredState <- dlmSvd2var(ext1\$U.C, ext1\$D.C)
#95% ci for prey
seFY <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[1,1])))
ciFY <- ext1\$m[,1] + qnorm(.05/2)*seFY%o%c(1, -1)
#95% ci for predators
seFD <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[2,2])))
ciFD <- ext1\$m[,2] + qnorm(.05/2)*seFD%o%c(1, -1)

#prey population
plot(st, x[,1], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFY[,1], lty=2, col="blue")
lines(c(0,mt), ciFY[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#predator population
plot(st, x[,2], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFD[,1], lty=2, col="blue")
lines(c(0,mt), ciFD[,2], lty=2, col="blue")
legend("topright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

##Diagnostics

#Diagnostics below are computed for the prey estimates only
#(diagnostics for the predator estimates are calculated in a similar way)

#raw residuals (=innovations)
plot(mt, yprey-ext1\$f[,1], type="l", col="gray", xlab="time", ylab="Raw residuals")
abline(h=0, lty=2, col="blue")

#standardized residuals
varcovPredState <- dlmSvd2var(ext1\$U.R, ext1\$D.R)
Qt <- sapply(1:nm, function (i)
ex1\$V + ext1\$dFF.dx[[i]]%*%varcovPredState[[i]]%*%t(ext1\$dFF.dx[[i]]),
simplify=FALSE)
sqrtQTx <- sqrt(unlist(lapply(Qt, function(x) x[1,1])))
serrx <- (yprey-ext1\$f[,1])/sqrtQTx
plot(mt, serrx, type="l", col="gray", xlab="time", ylab="Standardized residuals")
abline(h=0, lty=2, col="blue")

#qq-plot for standardized residuals
qqnorm(serrx, cex=.5)
abline(a=0, b=1, col="darkgray", lty=2, lwd=2)

#acf plot for standardized residuals
acf(serrx, lag.max=10, main="Standardized residuals")

##Compute smoothed state estimates

#apply fixed-interval smoothing
smoothState <- dlmExtSmooth(ext1)
plot(x[,1], x[,2], type="l", lwd=2, xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred), col=gray(level=.5))
lines(ext1\$m[,1], ext1\$m[,2], type="l", lty=2, col="blue")
lines(smoothState\$s[,1], smoothState\$s[,2], type="l", lty=2, col="darkgreen")
legend("topright", lty=c(1, 2, 2), lwd=c(1, 1, 1),
col=c("black", "blue", "darkgreen"),
legend=c("true state", "filtered state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#confidence intervals for smoothed state estimates
varcovSmoothState <- dlmSvd2var(smoothState\$U.S, smoothState\$D.S)
#95% ci for prey population
seSY <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[1,1])))
ciSY <- smoothState\$s[,1] + qnorm(.05/2)*seSY%o%c(1, -1)
#95% ci for predator population
seSD <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[2,2])))
ciSD <- smoothState\$s[,2] + qnorm(.05/2)*seSY%o%c(1, -1)

#plot prey
plot(st, x[,1], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciSY[,1], lty=2, col="blue")
lines(c(0,mt), ciSY[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#plot predators
plot(st, x[,2], type="l", xlab="time", ylab="Predator",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciSD[,1], lty=2, col="blue")
lines(c(0,mt), ciSD[,2], lty=2, col="blue")
legend("topright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

##EXAMPLE 2

#Similar to Example 1, but in this example we assume that we only
#measured the number of prey at each time step.

#Data for Example 2
dataEx2 <- yprey

#Dynamic model:
#specifying 2 states, namely [number of prey, number of predators].
#Note that the specified initial state estimates (at t=0) for prey and predators
#deviate from the values that were used for generating the data (i.e., 400
#for prey and 200 for predators).
ex2 <- list(m0=c(390, 220), #initial state estimates
#error covariances of the initial state estimates:
#this matrix reflects the uncertainty in our initial state estimates
#you may change the values in this matrix and see how they influence
#the behavior of the Kalman filter
C0=diag(rep(100,2)),
#measurement noise
V=sigmaPrey^2,
#process noise
W=diag(rep(0,2)))

#Specify the state transition function
#REMEMBER: always use arguments x and k when specifying the GGfunction
GGfunction <- function (x, k){
prey <- x[1]; pred <- x[2]
c(prey + (alpha*prey - beta*prey*pred)*dT,
pred + (delta*prey*pred - gamma*pred)*dT)}

#Specify the observation function
#REMEMBER: always use arguments x and k when specifying the FFfunction
FFfunction <- function (x, k){
prey <- x[1]
c(prey)}

##Compute the filtered (a posteriori) state estimates
ext2 <- dlmExtFilter(y=dataEx2, mod=ex2,
GGfunction=GGfunction, FFfunction=FFfunction)

#plot the filtered state estimates
plot(x[,1], x[,2], type="l", lwd=2, xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred), col=gray(level=.5))
lines(ext2\$m[,1], ext2\$m[,2], lty=2, col="blue", lwd=2)
legend("topright", lty=c(1, 2), lwd=c(2, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#compute the confidence intervals for the filtered state estimates
varcovFilteredState <- dlmSvd2var(ext2\$U.C, ext2\$D.C)
#95% ci for prey
seFY <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[1,1])))
ciFY <- ext2\$m[,1] + qnorm(.05/2)*seFY%o%c(1, -1)
#95% ci for predators
seFD <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[2,2])))
ciFD <- ext2\$m[,2] + qnorm(.05/2)*seFD%o%c(1, -1)

#prey population
plot(st, x[,1], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFY[,1], lty=2, col="blue")
lines(c(0,mt), ciFY[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#predator population
plot(st, x[,2], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFD[,1], lty=2, col="blue")
lines(c(0,mt), ciFD[,2], lty=2, col="blue")
legend("topright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

##Compute smoothed state estimates

#apply fixed-interval smoothing
smoothState <- dlmExtSmooth(ext2)
plot(x[,1], x[,2], type="l", lwd=2, xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred), col=gray(level=.5))
lines(ext2\$m[,1], ext2\$m[,2], type="l", lty=2, col="blue")
lines(smoothState\$s[,1], smoothState\$s[,2], type="l", lty=2, col="darkgreen")
legend("topright", lty=c(1, 2), lwd=c(1, 1),
col=c("black", "blue", "darkgreen"),
legend=c("true state", "filtered state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#confidence intervals for smoothed state estimates
varcovSmoothState <- dlmSvd2var(smoothState\$U.S, smoothState\$D.S)
#95% ci for prey population
seSY <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[1,1])))
ciSY <- smoothState\$s[,1] + qnorm(.05/2)*seSY%o%c(1, -1)
#95% ci for predator population
seSD <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[2,2])))
ciSD <- smoothState\$s[,2] + qnorm(.05/2)*seSY%o%c(1, -1)

#plot prey
plot(st, x[,1], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciSY[,1], lty=2, col="blue")
lines(c(0,mt), ciSY[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#plot predators
plot(st, x[,2], type="l", xlab="time", ylab="Predator",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciSD[,1], lty=2, col="blue")
lines(c(0,mt), ciSD[,2], lty=2, col="blue")
legend("topright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#EXAMPLE 3

#In the above examples we assumed that the values for the parameters
#beta and gamma were somehow known to us. Let us now include these
#two parameters as states in our model and have the EKF estimate their values.

#Data
dataEx3 <- cbind(yprey, ypred)

#Dynamic model:
#specifying 4 states, namely [number of prey, number of predators, beta, delta].
#Note that the specified initial state estimates (at t=0) deviate from the
#values that were used for generating the data.
#Remember that the values for beta and delta used in generating the data
#were 1/300 and 1/200, respectively.
ex3 <- list(m0=c(390, 220, 1/400, 1/100), #initial state estimates
#error covariances of the initial state estimates:
#this matrix reflects the uncertainty in our initial state estimates
#you may change the values in this matrix and see how they influence
#the behavior of the Kalman filter
C0=diag(rep(1000,4)),
#measurement noise
V=diag(c(sigmaPrey^2, sigmaPred^2)),
#process noise
W=diag(rep(0,4)))

#Specify the state transition function
#REMEMBER: always use arguments x and k when specifying the GGfunction
GGfunction <- function (x, k){
prey <- x[1]; pred <- x[2]; beta <- x[3]; delta <- x[4]
c(prey + (alpha*prey - beta*prey*pred)*dT,
pred + (delta*prey*pred - gamma*pred)*dT,
beta,
delta)}

#Specify the observation function
#REMEMBER: always use arguments x and k when specifying the FFfunction
FFfunction <- function (x, k){
prey <- x[1]; pred <- x[2]
c(prey, pred)}

##Compute the filtered (a posteriori) state estimates
ext3 <- dlmExtFilter(y=dataEx3, mod=ex3,
GGfunction=GGfunction, FFfunction=FFfunction)

#plot the filtered state estimates
plot(x[,1], x[,2], type="l", lwd=2, xlab="Prey", ylab="Predators",
xlim=range(yprey), ylim=range(ypred), col=gray(level=.5))
lines(ext3\$m[,1], ext3\$m[,2], lty=2, col="blue", lwd=2)
legend("topright", lty=c(1, 2), lwd=c(2, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#compute the confidence intervals for the filtered state estimates
varcovFilteredState <- dlmSvd2var(ext3\$U.C, ext3\$D.C)
#95% ci for prey
seFY <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[1,1])))
ciFY <- ext3\$m[,1] + qnorm(.05/2)*seFY%o%c(1, -1)
#95% ci for predators
seFD <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[2,2])))
ciFD <- ext3\$m[,2] + qnorm(.05/2)*seFD%o%c(1, -1)
#95% ci for beta
seFB <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[3,3])))
ciFB <- ext3\$m[,3] + qnorm(.05/2)*seFB%o%c(1, -1)
#95% ci for delta
seFE <- sqrt(unlist(lapply(varcovFilteredState, function(x) x[4,4])))
ciFE <- ext3\$m[,4] + qnorm(.05/2)*seFE%o%c(1, -1)

#prey population
plot(st, x[,1], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFY[,1], lty=2, col="blue")
lines(c(0,mt), ciFY[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#predator population
plot(st, x[,2], type="l", xlab="time", ylab="Prey",
col=gray(level=.7), lwd=1.5)
lines(c(0,mt), ciFD[,1], lty=2, col="blue")
lines(c(0,mt), ciFD[,2], lty=2, col="blue")
legend("topright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#beta
plot(x=NA, y=NA, type="n", xlab="time", ylab="Beta",
xlim=c(0,7), ylim=c(0, .005))
abline(h=beta, lty=2, col=gray(level=.5))
lines(c(0,mt), ciFB[,1], lty=2, col="blue")
lines(c(0,mt), ciFB[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

#delta
plot(x=NA, y=NA, type="n", xlab="time", ylab="Gamma",
xlim=c(0,7), ylim=c(.004, .006))
abline(h=delta, lty=2, col=gray(level=.5))
lines(c(0,mt), ciFE[,1], lty=2, col="blue")
lines(c(0,mt), ciFE[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "filtered state"),
bty="n", y.intersp=1.2, cex=.7)

##Compute smoothed state estimates

#apply fixed-interval smoothing
smoothState <- dlmExtSmooth(ext3)

#confidence intervals for smoothed state estimates
varcovSmoothState <- dlmSvd2var(smoothState\$U.S, smoothState\$D.S)
#95% ci for beta
seSD <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[3,3])))
ciSD <- smoothState\$s[,3] + qnorm(.05/2)*seSD%o%c(1, -1)
#95% ci for delta
seSE <- sqrt(unlist(lapply(varcovSmoothState, function(x) x[4,4])))
ciSE <- smoothState\$s[,4] + qnorm(.05/2)*seSE%o%c(1, -1)

#beta
plot(x=NA, y=NA, type="n", xlab="time", ylab="Beta",
xlim=c(0,7), ylim=c(0, .005))
abline(h=beta, lty=2, col=gray(level=.5))
lines(c(0,mt), ciSD[,1], lty=2, col="blue")
lines(c(0,mt), ciSD[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)

#delta
plot(x=NA, y=NA, type="n", xlab="time", ylab="Gamma",
xlim=c(0,7), ylim=c(.004, .006))
abline(h=delta, lty=2, col=gray(level=.5))
lines(c(0,mt), ciSE[,1], lty=2, col="blue")
lines(c(0,mt), ciSE[,2], lty=2, col="blue")
legend("bottomright", lty=c(1, 2),
col=c(gray(level=.5), "blue"),
legend=c("true state", "smoothed state"),
bty="n", y.intersp=1.2, cex=.7)
```

In my next two blog posts I will show 1) how to estimate unknown parameters of an EKF-model using likelihood maximization, and 2) how to predict future states and observations with the EKF.