KFAS (version 0.6.1)

kf: Kalman Filter With Exact Diffuse Initialisation

Description

Performs Kalman filtering with exact diffuse initialisation using univariate approach (also known as sequential processing). Written in Fortran, uses subroutines from BLAS and LAPACK. See function 'ks' for smoothing.

Usage

kf(yt,  Zt,  Tt,  Rt,  Ht,  Qt,  a1,  P1,  P1inf=0,  optcal=c(TRUE,
TRUE,  TRUE,  TRUE),  tol=1e-7)

Arguments

yt

Matrix or array of observations.

Zt

System matrix or array of observation equation.

Tt

System matrix or array of transition equation.

Rt

System matrix or array of transition equation.

Ht

Variance matrix or array of disturbance terms eps_t of observation equation.

Qt

Variance matrix or array of disturbance terms eta_t.

a1

Initial state vector.

P1

Variance matrix of a1. In diffuse case P1star, the non-diffuse part of P1.

P1inf

Diffuse part of P1. If non-zero, filtering starts with exact diffuse initialisation.

optcal

Vector of length 4. Calculate multivariate vt, Ft, Kt, Lt and their diffuse counterparts Finf, Fstar, Kinf, Kstar, Linf and Lstar. Default is c(TRUE,TRUE,TRUE,TRUE) which calculates all. Note that Kt cannot be calculated without Ft and Lt ca

tol

Tolerance parameter. Smallest covariance/variance value not counted for zero in diffuse phase. Default is 1e-7.

Value

• A list with the following elements:
• yt
  p*n matrix or array of observations.
• ydimt
  array of length n which has dimensions of observation vector at times t=1,...,n.
• tv
  array of length 4 where tv[i]=0 if i is time-invariant and otherwise tv[i]=1, i is dt, Tt, Rt, Qt.
• Zt
  system matrix or array of observation equation.
• Tt
  system matrix or array of transition equation.
• Rt
  system matrix or array of transition equation.
• Ht
  variance matrix or array of disturbance terms eps_t of observation equation.
• Qt
  variance matrix or array of disturbance terms eta_t.
• a1
  initial state vector.
• P1
  variance matrix of a1. In diffuse case P1star, the non-diffuse part of P1 .
• at
  m*(n+1) array of E(alpha_t | y_1, y_2, ... , y_t-1).
• Pt
  m*m*(n+1) array of Var(alpha_t | y_1, y_2, ... , y_t-1).
• vtuni
  p*1*n array of vt of univariate approach.
• Ftuni
  Ft of univariate approach, Var(vtuni).
• Ktuni
  m*p*n array of Kalman gain of univariate approach.
• Pinf, Pstar
  p*p*d+1 arrays of diffuse phase decomposition of Pt.
• Finfuni, Fstaruni
  p*p*d (p*d arrays of diffuse phase decomposition of Ftuni.
• Kinfuni, Kstaruni
  m*p*d arrays of diffuse phase decomposition of Ktuni.
• d
  the last index of diffuse phase, ie. the non-diffuse phase began from time d+1.
• j
  the index of last y_t,i of diffuse phase.
• p
  the dimension of observation vector.
• m
  the dimension of state vector.
• r
  the dimension of variance matrix Qt.
• n
  the number of observations.
• lik
  Value of the log-likelihood function. If NaN, Ftuni_i,t was zero at some t=d+1,...,n or Finfuni_i,t and Fstaruni_i,t was zero at some t=1,...,d.
• optcal
• info
  if info[1]=1, could not diagonalize Ht. If info[i]=1, i=2,3,4, Finf, Fstar or Ft was singular.
• vt
  p*1*n array of vt = yt - Zt * at.
• Ft
  p*p*n array of Ft = Var(vt) of Kalman filter.
• Kt
  m*p*n array of Kalman gain: Kt = Tt * Pt * Zt' * Ft^-1.
• Lt
  the m*m*n array, Lt = Tt - Kt * Zt.
• Finf, Fstar
  p*p*d arrays of diffuse phase decomposition of Ft.
• Kinf, Kstar
  m*p*d arrays of diffuse phase decomposition of Kt.
• Linf, Lstar
  m*m*d arrays of diffuse phase decomposition of Lt.
• tol
  Tolerance parameter.

Details

Function kf performs Kalman filtering of gaussian multivariate state space model using the univariate approach from Koopman and Durbin (2000, 2001). Univariate approach is also known as sequential processing, see Anderson and Moore (1979). In case where the distributions of some or all of the elements of initial state vector are not known, kf uses exact diffuse initialisation using univariate approach by Koopman and Durbin (2000, 2003). Note that in univariate approach the prediction error variance matrices Ft, Finf and Fstar does not need to be non-singular, as there is no matrix inversions in univariate approach algorithm. This provides faster and more general filtering than normal multivariate Kalman filter algorithm.

Filter can deal partially or totally missing observation vectors. If y_t,i is NA, it is interpreted as missing value, and the dimensions of vtuni and Ftuni (or Fstaruni and Finfuni) are decreased and the corresponding elements of vt are marked as NA.

The state space model is given by

y_t = Z_t * alpha_t + eps_t (observation equation)

alpha_t+1 = T_t * alpha_t + R_t * eta_t(transition equation)

where eps_t ~ N(0,H_t) and eta_t ~ N(0,Q_t). Note that error terms eps_t and eta_t are assumed to be uncorrelated.

When P1inf is non-zero, exact diffuse initialisation is used. Matrices Pt, Kt, Ft and Lt are decomposed to Pinf, Pstar, Kinf, Kstar, Finf, Fstar, Linf and Lstar. Diffuse phase is continued until Pinf becomes zero-matrix. See Koopman and Durbin (2000, 2001, 2003) for details for exact diffuse and non-diffuse filtering.

Notice that vtuni, Ftuni, Fstaruni, Finfuni, Ktuni, Kinfuni and Kstaruni are usually not the same as those calculated in usual multivariate Kalman filter. Also the Lt, Linf and Lstar are not calculated explicitly. If usual vt, Ft, Finf, Fstar, Kt, Kinf, Kstar, Lt, Linf and Lstar are needed, use optcal=c(TRUE, TRUE, TRUE, TRUE). When estimating parameters, it is suggested to use optcal=c(FALSE,FALSE,FALSE,FALSE) for maximum speed.

Dimensions of variables are: 'yt' p*n 'Zt' p*m or p*m*n 'Tt' m*m or m*m*n 'Rt' m*r or m*r*n 'Ht' p*p or p*p*n 'Qt' r*r or r*r*n 'a1' m*1 'P1' m*m 'P1inf' m*m where p is dimension of observation vector, m is dimension of state vector and n is number of observations.

References

Koopman, S.J. and Durbin J. (2000). Fast filtering and smoothing for non-stationary time series models, Journal of American Statistical Assosiation, 92, 1630-38. Koopman, S.J. and Durbin J. (2001). Time Series Analysis by State Space Methods. Oxford: Oxford University Press. Koopman, S.J. and Durbin J. (2003). Filtering and smoothing of state vector for diffuse state space models, Journal of Time Series Analysis, Vol. 24, No. 1. Shumway, Robert H. and Stoffer, David S. (2006). Time Series Analysis and Its Applications: With R examples.

Examples

library(KFAS)

#Example of local level model
#Using the Nile observations
data(Nile)

yt<-t(data.matrix(Nile))
s2_eps<-15099
s2_eta<-1469.1

f.out<-kf(yt = yt,  Zt = 1,  Tt=1,  Rt=1,  Ht= s2_eps,  Qt=s2_eta,  a1 =
0,  P1=1e7)

#a1 and P1 are not really estimated, 
#should actually use exact diffuse initialisation:

fd.out<-kf(yt = yt,  Zt = 1,  Tt=1,  Rt=1,  Ht=s2_eps,  Qt=s2_eta,  a1 =
0, P1=0,  P1inf=1)

#No stationary elements, P1=0, P1inf=1

#Plotting observations, non-diffuse and diffuse at,  not plotting the a1=0:

ts.plot(Nile,  ts(f.out$at[1,2:length(Nile)], start=1872, end=1970),
ts(fd.out$at[1,2:length(Nile)], start=1872, end=1970), col=c(1,2,3))

#Looks identical. Actually start of series differs little bit:

f.out$at[1,1:20]
fd.out$at[1,1:20]



#Example of multivariate local level model
#Two series of average global temperature deviations for years 1880-1987
#See Shumway and Stoffer (2006), p. 327 for details

data(GlobalTemp)
yt<-array(GlobalTemp,c(2,108))


#Estimating the variance parameters

likfn<-function(par, yt, a1, P1, P1inf) #Function to optim
{
L<-matrix(c(par[1],par[2],0,par[3]),ncol=2)
H<-L%*%t(L)
q11<-exp(par[4])
lik<-kf(yt = yt, Zt = matrix(1,nrow=2), Tt=1, Rt=1, Ht=H, Qt=q11, a1 =
a1, P1=P1, P1inf=P1inf, optcal=c(FALSE,FALSE,FALSE,FALSE))
lik$lik
}

#Diffuse initialisation
#Notice that diffuse initialisation without univariate approach does not
#work here because Finf is non-singular and non-zero

est<-optim(par=c(.1,0,.1,.1), likfn, method="BFGS",
control=list(fnscale=-1), hessian=TRUE, yt=yt, a1=0, P1=0, P1inf=1) 

pars<-est$par
L<-matrix(c(pars[1],pars[2],0,pars[3]),ncol=2)
H<-L%*%t(L)
q11<-exp(pars[4])

kfd<-kf(yt = yt, Zt = matrix(1,nrow=2), Tt=1, Rt=1, 
	Ht=H, Qt=q11, a1 = 0, P1=0, P1inf=1)

#Same as non-diffuse, initial values from Shumway and Stoffer (2006):

est<-optim(par=c(.1,0,.1,.1), likfn, method="BFGS",
control=list(fnscale=-1), hessian=TRUE, yt=yt, a1=-0.35, P1=0.01, P1inf=0)

pars<-est$par
L<-matrix(c(pars[1],pars[2],0,pars[3]),ncol=2)
H<-L%*%t(L)
q11<-exp(pars[4])

kfnd<-kf(yt = yt, Zt = matrix(1,nrow=2), Tt=1, Rt=1, Ht=H, Qt=q11, 
a1 = -0.35, P1=0.01, P1inf=0)

kfd$Qt
kfnd$Qt
kfd$Ht
kfnd$Ht
#Estimated Qt and Ht differs

ts.plot(ts(yt[1,],start=1880),ts(yt[2,],start=1880),ts(kfd$at[1,],start=1880),ts(kfnd$at[1,],start=1880),col=c(1,2,3,4))
#differs at start


#Example of stationary ARMA(1,1)

n<-1000
ar1<-0.8 
ma1<-0.3
sigma<-0.5
yt<-arima.sim(model=list(ar=ar1,ma=ma1), n=n, sd=sigma)

dt <- matrix(0, nrow = 2)
Zt<-matrix(c(1,0),ncol=2)
Tt<-matrix(c(ar1,0,1,0),ncol=2)
Rt<-matrix(c(1,ma1),ncol=1)
a1<-matrix(c(0,0),ncol=1)
P1<-matrix(0,ncol=2,nrow=2)
P1[1,1]<-1/(1-ar1^2)*(1+ma1^2+2*ar1*ma1)
P1[1,2]<-ma1
P1[2,1]<-ma1
P1[2,2]<-ma1^2
f.out<-kf(yt = array(yt,dim=c(1,n)), Zt = Zt, Tt=Tt, Rt=Rt, Ht= 0,
Qt=sigma^2, a1 = a1, P1=P1)