R code for forecasting with the Gauss-Hermite Kalman filter

In one of my previous blog posts I showed how to implement and apply the Gauss-Hermite Kalman Filter (GHKF) in R.
In this post I will demonstrate how to predict future system states and observations with the GHKF.

Gauss-Hermite Kalman filter forecasting

Suggestions and/or questions? Please contact Stefan Gelissen (email: info at datall-analyse.nl).
See this page for an overview of all of Stefan’s R code blog posts.

R code

library(fastGHQuad)

#The R functions for the Gauss-Hermite Kalman filter (GHKF)
#can be downloaded from: 
source("http://www.datall-analyse.nl/R/GHKF.R")
#Take a look at the functions GHKF (=the GHKF) and GHKFforecast, and notice
#that at the beginning of the scripts you will find
#a description of the functions' arguments.
GHKF
GHKFforecast



##EXAMPLE 1

#The following example is similar to an example discussed 
#by Arasaratnam, Haykin, and Elliott in their 2007
#paper "Discrete-Time Nonlinear Filtering Algorithms Using
#Gauss-Hermite Quadrature".

#In this example we will estimate the state of the following
#nonlinear discrete-time system:
#x(k+1) = 0.2*x(k) + 0.01*x(k)^2 + 8*cos(1.2*(k+1)) + W
#where x(k+1) is the state at time step k+1, x(k) the state at time
#step k, and W the process noise (variance). We will set W to 10 in
#this example.

#The measurement equation is:
#y(k) = x(k)^2 + V
#where V is the measurement noise (variance). We will set V to 0.01.



##TRUE STATES AND MEASUREMENTS

#sample time
dT <- 1
#you may change the value of dT and see how it influences
#the behavior of the Gauss-Hermite Kalman filter

#observation times
t <- seq(1, 50, dT)

#state at t=0
X0 <- 0

#process noise (standard deviation)
Wk <- sqrt(10)

#measurement noise (standard deviation)
Vk <- sqrt(.01)

#parameters in the dynamic model
phi1 <- 0.2
phi2 <- 0.01

#simulate true state for x at the observation times
Xk <- matrix(NA, nrow=length(t)+1, ncol=1)
Yk <- matrix(NA, nrow=length(t), ncol=1)
Xk[1] <- X0

for (i in 2:(length(t)+1)) {
  Xk[i] <- phi1*Xk[i-1] + phi2*Xk[i-1]^2 + 8*cos(1.2*(i)) + rnorm(1)*Wk
  Yk[i-1] <- Xk[i]^2 + rnorm(1)*Vk}

#plot simulated states
plot(c(0,t), Xk, type="o", xlab="Time, k", ylab="x(k)")

#plot measurements
plot(t, Yk, type="o", xlab="Time, k", ylab="y(k)")

#store measurement data
dataEx1 <- Yk


##FORECASTING WITH THE GAUSS-HERMITE KALMAN FILTER (GHKF)

#Dynamic model:
#specifying 1 state, namely [x].
#Note that the specified initial state estimate (at t=0) below for x
#deviates from the value that was used for generating the data (i.e., 0).
#You may change this initial state estimate too and see how it
#influences the behavior of the Gauss-Hermite Kalman filter.
ex1 <- list(m0=.5, #initial state estimate
            #error covariance of the initial state estimate:
            #this vector reflects the uncertainty in our initial state estimate
            #you may change the value in this vector and see how it influences
            #the behavior of the Kalman filter
            C0=5,
            #measurement noise
            V=Vk^2,
            #process noise
            W=Wk^2)

#Specify the state transition function:
#WARNING: always use arguments x and k when specifying the GGfunction
GGfunction <- function (x, k){
  xk <- x[1]
  phi1*xk + phi2*xk^2 + 8*cos(1.2*(k+1))}

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


##Compute the filtered (a posteriori) state estimates with the GHKF
#If the filter generates an error, then try to decrease the 
#error covariances of the initial state estimates (as specified in C0).
#(For instance, set the value for x from 5 to 1)
ghkf1 <- GHKF(y=dataEx1, mod=ex1, order=3, sqrtMethod="Cholesky",
              GGfunction=GGfunction, FFfunction=FFfunction)


#plot the filtered state
plot(c(0,t), Xk, type="o", col=c(gray(level=.5)), ylim=range(ghkf1$m[,1]),
     ylab="x(k)", xlab="Time, k")
lines(c(0,t), ghkf1$m[,1], lty=2, col="blue", lwd=1)
legend("topright", lty=c(1, 2),
       col=c(gray(level=.5), "blue"),
       legend=c("true state", "GHKF"),
       bty="n", y.intersp=1.2, cex=.7)




##Compute the filtered (a posteriori) state estimates
#GHKF for the first 40 observations
ghkf2 <- GHKF(y=dataEx1[1:40], mod=ex1, order=3, sqrtMethod="Cholesky",
              GGfunction=GGfunction, FFfunction=FFfunction)

#Compute forecasts for 10 time steps ahead
#(If this forecasting function generates an error, then decrease nAhead
#from 10 to, for instance, 5)
fc <- GHKFforecast(ghkf2, nAhead=10)


#95% confidence intervals for forecasts

#Future system states
seFoX <- sqrt(unlist(lapply(fc$R, function(x) x[1,1])))
ciFoX <- fc$a[,1] + qnorm(.05/2)*seFoX%o%c(1, -1)

#Plot future system states including confidence intervals
plot(c(0,t), Xk, type="o", lty=2, col="red",
     ylim=range(ciFoX),
     xlab="Time, k", ylab="x(k)")
lines(t[1:40], ghkf2$m[-1,1], col=gray(level=.5), lwd=1)
lines(t[41:50], fc$a[,1], col="blue", lwd=2)
lines(t[41:50], ciFoX[,1], col="blue", lty=2, lwd=1)
lines(t[41:50], ciFoX[,2], col="blue", lty=2, lwd=1)
legend("bottomleft", lty=c(2, 1), 
       col=c("red", "blue"),
       legend=c("true state", "future state"),
       bty="n", y.intersp=1.2, cex=.7)


#Future observations
seFoY <- sqrt(unlist(lapply(fc$Q, function(x) x[1,1])))
ciFoY <- fc$f[,1] + qnorm(.05/2)*seFoY%o%c(1, -1)

#Plot future observations including prediction intervals
plot(c(0,t), Xk^2, type="o", lty=2, col="red", ylim=range(ciFoY),
     xlab="Time, k", ylab="y(k)")
lines(t, Yk, type="l", col=gray(level=.7)) #measurements
lines(t[1:40], ghkf2$f[,1], col=gray(level=.1), lwd=1)
lines(t[41:50], fc$f[,1], col="blue", lwd=2)
lines(t[41:50], ciFoY[,1], col="blue", lty=2, lwd=1)
lines(t[41:50], ciFoY[,2], col="blue", lty=2, lwd=1)
legend("topleft", lty=c(2, 1, 1), 
       col=c("red", col=gray(level=.7), "blue"),
       legend=c("true state", "measurements", "future observations"),
       bty="n", y.intersp=1.2, cex=.7)





##EXAMPLE 2

#The following example can be found in Crassidis and Junkins's
#2012 book "Optimal estimation of dynamic systems".

#In this example we will estimate the state of the following
#nonlinear discrete-time system:
#x(k+1) = x(k)/2 + 25*x(k)/(1+x(k)^2) + 8*cos(1.2*k) + W
#where x(k+1) is the state at time step k+1, x(k) the state at time
#step k, and W the process noise (variance). We will set W to 10^2 in
#this example.

#Note, that the nonlinearity of this dynamic system is more severe
#than that of the dynamic system in Example 1.

#The measurement equation is:
#y(k) = x(k)^2/20 + V
#where V is the measurement noise (variance). We will set V to 1^2.



##TRUE STATES AND MEASUREMENTS

#sample time
dT <- 1
#you may change the value of dT and see how it influences
#the behavior of the Gauss-Hermite Kalman filter

#observation times
t <- seq(1, 75, dT)

#state at t=0
X0 <- 0

#process noise (standard deviation)
Wk <- 10

#measurement noise (standard deviation)
Vk <- 1

#simulate true state for x at the observation times
Xk <- matrix(NA, nrow=length(t)+1, ncol=1)
Yk <- matrix(NA, nrow=length(t), ncol=1)
Xk[1] <- X0

for (i in 2:(length(t)+1)) {
  Xk[i] <- Xk[i-1]/2 + (25*Xk[i-1])/(1+Xk[i-1]^2) + 8*cos(1.2*(i-1)) + rnorm(1)*Wk
  Yk[i-1] <- (Xk[i]^2)/20 + rnorm(1)*Vk}

#plot simulated states
plot(c(0,t), Xk, type="o", xlab="Time, k", ylab="x(k)")

#plot measurements
plot(t, Yk, type="o", xlab="Time, k", ylab="y(k)")

#store measurement data
dataEx2 <- Yk



##FORECASTING WITH THE GAUSS-HERMITE KALMAN FILTER (GHKF)

#Dynamic model:
#specifying 1 state, namely [x].
#Note that the specified initial state estimates (at t=0) below for x
#deviates from the value that was used for generating the data (i.e., 0).
#You may change this initial state estimate too and see how it
#influences the behavior of the Gauss-Hermite Kalman filter.
ex2 <- list(m0=2, #initial state estimate
            #error covariance of the initial state estimate:
            #this vector reflects the uncertainty in our initial state estimate
            #you may change the value in this vector and see how it influences
            #the behavior of the Kalman filter
            C0=5^2,
            #measurement noise
            V=Vk^2,
            #process noise
            W=Wk^2)

#Specify the state transition function:
#WARNING: always use arguments x and k when specifying the GGfunction
GGfunction <- function (x, k){
  x/2 + (25*x)/(1+x^2) + 8*cos(1.2*(k))}

#Specify the observation/measurement function:
#WARNING: always use arguments x and k when specifying the FFfunction
FFfunction <- function (x, k){
  (x^2)/20}



##Compute the filtered (a posteriori) state estimates with the GHKF
#If the filter generates an error, then try to decrease the 
#error covariance of the initial state estimate (as specified in C0).
#(For instance, set the value for x from 5^2 to 10)
ghkf3 <- GHKF(y=dataEx2, mod=ex2, order=20, sqrtMethod="Cholesky",
              GGfunction=GGfunction, FFfunction=FFfunction)


#plot the filtered state
plot(c(0,t), Xk, type="o", col=c(gray(level=.5)),
     ylim=range(cbind(Xk, ghkf3$m[,1])),
     ylab="x(k)", xlab="Time, k")
lines(c(0,t), ghkf3$m[,1], lty=2, col="blue", lwd=1)
legend("topright", lty=c(1, 2),
       col=c(gray(level=.5), "blue"),
       legend=c("true state", "GHKF"),
       bty="n", y.intersp=1.2, cex=.7)



##Compute the filtered (a posteriori) state estimates
#GHKF for the first 40 observations
ghkf4 <- GHKF(y=dataEx2[1:65], mod=ex2, order=20, sqrtMethod="Cholesky",
              GGfunction=GGfunction, FFfunction=FFfunction)

#Compute forecasts for 10 time steps ahead
#(If this forecasting function generates an error, then decrease nAhead
#from 10 to, for instance, 5)
fc <- GHKFforecast(ghkf4, nAhead=10)


#95% confidence intervals for forecasts

#Future system states
seFoX <- sqrt(unlist(lapply(fc$R, function(x) x[1,1])))
ciFoX <- fc$a[,1] + qnorm(.05/2)*seFoX%o%c(1, -1)

#Plot future system states including confidence intervals
plot(c(0,t), Xk, type="o", lty=2, col="red",
     ylim=range(ciFoX),
     xlab="Time, k", ylab="x(k)")
lines(t[1:65], ghkf4$m[-1,1], col=gray(level=.5), lwd=1)
lines(t[66:75], fc$a[,1], col="blue", lwd=2)
lines(t[66:75], ciFoX[,1], col="blue", lty=2, lwd=1)
lines(t[66:75], ciFoX[,2], col="blue", lty=2, lwd=1)
legend("bottomleft", lty=c(2, 1), 
       col=c("red", "blue"),
       legend=c("true state", "future state"),
       bty="n", y.intersp=1.2, cex=.7)


#Future observations
seFoY <- sqrt(unlist(lapply(fc$Q, function(x) x[1,1])))
ciFoY <- fc$f[,1] + qnorm(.05/2)*seFoY%o%c(1, -1)

#Plot future observations including prediction intervals
plot(c(0,t), (Xk^2)/20, type="o", lty=2, col="red", ylim=range(ciFoY),
     xlab="Time, k", ylab="y(k)")
lines(t, Yk, type="l", col=gray(level=.7)) #measurements
lines(t[1:65], ghkf4$f[,1], col=gray(level=.1), lwd=1)
lines(t[66:75], fc$f[,1], col="blue", lwd=2)
lines(t[66:75], ciFoY[,1], col="blue", lty=2, lwd=1)
lines(t[66:75], ciFoY[,2], col="blue", lty=2, lwd=1)
legend("topleft", lty=c(2, 1, 1), 
       col=c("red", col=gray(level=.7), "blue"),
       legend=c("true state", "measurements", "future observations"),
       bty="n", y.intersp=1.2, cex=.7)