Thursday, 28 December 2017

Gaussian State Space Model in R using Kalman Filtering and Smoothing

State space model (SSM) refers to a class of probabilistic graphical model (Koller and Friedman, 2009) that describes the probabilistic dependence between the latent state variable and the observed measurement. The state or the measurement can be either continuous or discrete. The term “state space” originated in 1960s in the area of control engineering (Kalman, 1960). SSM provides a general framework for analyzing deterministic and stochastic dynamical systems that are measured or observed through a stochastic process. Now what is state space and dynamic system?

State space is the set of all possible states of a dynamical system; each state of the system corresponds to a unique point in the state space. For example, the state of an idealized pendulum is uniquely defined by its angle and angular velocity, so the state space is the set of all possible pairs (angle, velocity). A dynamical system is a rule for time evolution on a state space.
The most well studied SSM is the Kalman filter, which defines an optimal algorithm for inferring linear Gaussian systems.

Like hidden markov models, in SSM we have ,

  1. Observation Equation
  2. State Equation.
where yt is observations at time t. Zt, Tt and Rt are system matrices.  Alpha is latent state variable. et is normally distributed with 0 mean and Ht standard deviation (sd) and nt is also normally distributed with 0 mean and Qt sd.

Example of Gaussian State space model ( using KFAS package in R)-

Providing Initial values of parameters-

deaths <- window(alcohol[,2], end =2007)
population <- window(alcohol[,6], end = 2007)
## defining all system matrices
Zt <- matrix(c(1,0),1,2) #matrix(data = NA, nrow = 1, ncol = 1)
Ht<- matrix(NA)
Tt <- matrix(c(1,0,1,1),2,2)
Rt <- matrix(c(1,0), 2,1)
Qt <- matrix(NA)
a1 <- matrix(c(1,0),2,1)
P1 <- matrix(0,2,2)
P1inf <- diag(2)

Create a State Space Model Object of Class SSModel-

model_gaussian <- SSModel(deaths / population ~ -1 + SSMcustom(Z = Zt, T = Tt, R = Rt, Q = Qt, a1 = a1, P1 = P1,P1inf = P1inf),H = Ht)

y; univariate time-series is deaths / population. Function SSModel (SSMarima, SSMcustom , SSMcycle, SSMregression, SSMseasonal¸ SSMtrend)  creates a state space object of class SSModel which can be used as an input object for various functions like fitSSM of KFAS package.

Fitting the model-

fit_gaussian <- fitSSM(model_gaussian, inits = c(0, 0), method = "BFGS")

Function fitSSM finds the maximum likelihood estimates for unknown parameters of an arbitary state space model, given the user-defined model updating function. One of the arguments of fitSSM is ‘updatefn’ which is used to update parameters. This is wrapper on optim function which is a General-purpose optimization based on Nelder–Mead, quasi-Newton and conjugate-gradient algorithms.

out_gaussian <- KFS(fit_gaussian$model)

This will provide state filtering and state smoothing. KFS will also do next step forecasting. Forecasted value can be taken from ts_data_forecasted <- out_gaussian$muhat.The plot of fitted and actual values can be drawn by below code-

ts_data_forecasted <- ts_data_forecasted[1: length(ts_data_forecasted)-1]
x<- 1:length(ts_data_forecasted)
alcohol_plots <- data.frame(X=x,actuals=ts_data,forecasted=ts_data_forecasted)

ggplot(alcohol_plots, aes(x=X)) +
  geom_line(aes(y = forecasted, colour = "red")) + geom_line(aes(y = actuals, colour = "blue"))

red-forecasted & black-actual values
More about state based Models ( Hidden Markov) -

HMM- 1 ( Hidden Markov Model 1)

HMM- 2 ( Hidden Markov Models 2)