KFAS (version 1.0.3)

KFAS: KFAS: Functions for Gaussian and Non-Gaussian State Space Models

Description

Package KFAS contains functions for Kalman filtering, smoothing and simulation of linear state space models with exact diffuse initialization.

Details

The linear gaussian state space model is given by $$y_t = Z_t \alpha_t + \epsilon_t,$$ $$\alpha_{t+1} = T_t \alpha_t + R_t \eta_t,$$ where $\epsilon_t \sim N(0,H_t)$, $\eta_t \sim N(0,Q_t)$ and $\alpha_1 \sim N(a_1,P_1)$ independently of each other. All system and covariance matrices Z, H, T, R and Q can be time-varying, and partially or totally missing observations $y_t$ are allowed. Covariance matrices H and Q has to be positive semidefinite (although this is not checked). Dimensions of system matrices are rl{ Z $p \times m \times 1$ or $p \times m \times n$ in time varying case H $p \times p \times 1$ or $p \times p \times n$ in time varying case (Omitted in non-gaussian models) T $m \times m \times 1$ or $m \times m \times n$ in time varying case R $m \times k \times 1$ or $m \times k \times n$ in time varying case Q $k \times k \times 1$ or $k \times k \times n$ in time varying case u $n \times p$ (Omitted in gaussian models) } In case of any of the series in model is defined as non-gaussian, the observation equation is of form $$\prod_i^p p_i(y_{i,t}|\theta_t)$$ with $\theta_{i,t}=Z_{i,t}\alpha_t$ being one of the following: If observations $y_{i,1},\ldots,y_{i,n}$ are distributed as $N(\mu_t,u_t)$, then $\theta_t=\mu_t$. Note that now variances are defined using u, not H. If correlation between gaussian observation equations is needed, one can use $u_t=0$ and add correlating disturbances into state equation (although care is needed when making inferences as then $y_t=\theta_t$) If observations are distributed as $Poisson(u_t\lambda_t)$, where $u_t$ is offset term, then $\theta_t = log(u_t\lambda_t)$. If observations are distributed as $binomial(u_t,\pi_t)$, then $\theta_t = log[\pi_t/(1-\pi_t)]$, where $\pi_t$ is the probability of success at time $t$. If observations are distributed as $gamma(u_t,\mu_t)$, then $\theta_t = log(\mu_t)$, where $\mu[t]$ is the mean parameter and $u$ is the shape parameter. If observations are distributed as $negative binomial(u_t,\mu_t)$ (with expected value $\mu_t$ and variance $\mu_t+ \mu_t^2/u_t$, see dbinom), then $\theta_t = log[\mu_t]$. For exponential family models $u_t=1$ as a default. For completely gaussian models, parameter is omitted. For the unknown elements of initial state vector $a_1$, KFS uses exact diffuse initialization by Koopman and Durbin (2000, 2001, 2003), where the unknown initial states are set to have a zero mean and infinite variance, so $$P_1 = P_{\ast,1} + \kappa P_{\infty,1},$$ with $\kappa$ going to infinity and $P_{\infty,1}$ being diagonal matrix with ones on diagonal elements corresponding to unknown initial states. Diffuse phase is continued until rank of $P_{\infty,t}$ becomes zero. Rank of $P_{\infty}$ decreases by 1, if $F_\infty>tol>0$. Usually the number of diffuse time points equals the number unknown elements of initial state vector, but missing observations or time-varying Z can affect this. See Koopman and Durbin (2000, 2001, 2003) for details for exact diffuse and non-diffuse filtering. To lessen the notation and storage space, KFAS uses letters P, F and K for non-diffuse part of the corresponding matrices, omitting the asterisk in diffuse phase. All functions of KFAS use the univariate approach (also known as sequential processing, see Anderson and Moore (1979)) which is from Koopman and Durbin (2000, 2001). In univariate approach the observations are introduced one element at the time. Therefore the prediction error variance matrices F and Finf does not need to be non-singular, as there is no matrix inversions in univariate approach algorithm. This provides more stable and possibly more faster filtering and smoothing than normal multivariate Kalman filter algorithm. If covariance matrix H is not diagonal, it is possible to transform the model by either using LDL decomposition on H, or augmenting the state vector with $\epsilon$ disturbances. See transformSSM for more details.

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

# Example of local level model for Nile series

modelNile<-SSModel(Nile~SSMtrend(1,Q=list(matrix(NA))),H=matrix(NA))
modelNile
modelNile<-fitSSM(inits=c(log(var(Nile)),log(var(Nile))),model=modelNile,
                  method='BFGS',control=list(REPORT=1,trace=1))$model
# Filtering and state smoothing
out<-KFS(modelNile,filtering='state',smoothing='state')
out

# Confidence and prediction intervals for the expected value and the observations.
# Note that predict uses original model object, not the output from KFS.
conf<-predict(modelNile,interval='confidence')
pred<-predict(modelNile,interval='prediction')

ts.plot(cbind(Nile,pred,conf[,-1]),col=c(1:2,3,3,4,4),
        ylab='Predicted Annual flow', main='River Nile')


# Missing observations, using same parameter estimates

y<-Nile
y[c(21:40,61:80)]<-NA
modelNile<-SSModel(y~SSMtrend(1,Q=list(modelNile$Q)),H=modelNile$H)

out<-KFS(modelNile,filtering='mean',smoothing='mean')

# Filtered and smoothed states
plot.ts(cbind(y,fitted(out,filtered=TRUE),fitted(out)), plot.type='single',
        col=1:3, ylab='Predicted Annual flow', main='River Nile')


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

data(GlobalTemp)

model<-SSModel(GlobalTemp~SSMtrend(1,Q=NA,type='common'),H=matrix(NA,2,2))

# Estimating the variance parameters
inits<-chol(cov(GlobalTemp))[c(1,4,3)]
inits[1:2]<-log(inits[1:2])
fit<-fitSSM(inits=c(0.5*log(.1),inits),model=model,method='BFGS')

out<-KFS(fit$model)

ts.plot(cbind(model$y,coef(out)),col=1:3)
legend('bottomright',legend=c(colnames(GlobalTemp), 'Smoothed signal'), col=1:3, lty=1)



# Seatbelts data
model<-SSModel(log(drivers)~SSMtrend(1,Q=list(NA))+
               SSMseasonal(period=12,sea.type='trigonometric',Q=NA)+
               log(PetrolPrice)+law,data=Seatbelts,H=NA)

# As trigonometric seasonal contains several disturbances which are all
# identically distributed, default behaviour of fitSSM is not enough,
# as we have constrained Q. We can either provide our own
# model updating function with fitSSM, or just use optim directly:

# option 1:
ownupdatefn<-function(pars,model,...){
  model$H[]<-exp(pars[1])
  diag(model$Q[,,1])<-exp(c(pars[2],rep(pars[3],11)))
  model #for option 2, replace this with -logLik(model) and call optim directly
}

fit<-fitSSM(inits=log(c(var(log(Seatbelts[,'drivers'])),0.001,0.0001)),
            model=model,updatefn=ownupdatefn,method='BFGS')

out<-KFS(fit$model,smoothing=c('state','mean'))
out
ts.plot(cbind(out$model$y,fitted(out)),lty=1:2,col=1:2,
main='Observations and smoothed signal with and without seasonal component')
lines(signal(out,states=c("regression","trend"))$signal,col=4,lty=1)
legend('bottomleft',
legend=c('Observations', 'Smoothed signal','Smoothed level'),
col=c(1,2,4), lty=c(1,2,1))


# Multivariate model with constant seasonal pattern,
# using the the seat belt law dummy only for the front seat passangers,
# and restricting the rank of the level component by using custom component

# note the small inconvinience in regression component,
# you must remove the intercept from the additional regression parts manually

model<-SSModel(log(cbind(front,rear))~ -1 + log(PetrolPrice) + log(kms)
               + SSMregression(~-1+law,data=Seatbelts,index=1)
               + SSMcustom(Z=diag(2),T=diag(2),R=matrix(1,2,1),
                           Q=matrix(1),P1inf=diag(2))
               + SSMseasonal(period=12,sea.type='trigonometric'),
                 data=Seatbelts,H=matrix(NA,2,2))

likfn<-function(pars,model,estimate=TRUE){
  model$H[,,1]<-exp(0.5*pars[1:2])
  model$H[1,2,1]<-model$H[2,1,1]<-tanh(pars[3])*prod(sqrt(exp(0.5*pars[1:2])))
  model$R[28:29]<-exp(pars[4:5])
  if(estimate) return(-logLik(model))
  model
}
fit<-optim(f=likfn,p=c(-7,-7,1,-1,-3),method='BFGS',model=model)
model<-likfn(fit$p,model,estimate=FALSE)
model$R[28:29,,1]%*%t(model$R[28:29,,1])
model$H

out<-KFS(model)
out
ts.plot(cbind(signal(out,states=c('custom','regression'))$signal,model$y),col=1:4)

# For confidence or prediction intervals, use predict on the original model
pred <- predict(model,states=c('custom','regression'),interval='prediction')
ts.plot(pred$front,pred$rear,model$y,col=c(1,2,2,3,4,4,5,6),lty=c(1,2,2,1,2,2,1,1))

 # Poisson model
model<-SSModel(VanKilled~law+SSMtrend(1,Q=list(matrix(NA)))+
               SSMseasonal(period=12,sea.type='dummy',Q=NA),
               data=Seatbelts, distribution='poisson')

# Estimate variance parameters
fit<-fitSSM(inits=c(-4,-7,2), model=model,method='BFGS')

model<-fit$model

# use approximating model, gives posterior mode of the signal and the linear predictor
out_nosim<-KFS(model,smoothing=c('signal','mean'),nsim=0)
# State smoothing via importance sampling
out_sim<-KFS(model,smoothing=c('signal','mean'),nsim=1000)

out_nosim
out_sim

# Example of generalized linear modelling with KFS

# Same example as in ?glm
counts <- c(18,17,15,20,10,20,25,13,12)
outcome <- gl(3,1,9)
treatment <- gl(3,3)
print(d.AD <- data.frame(treatment, outcome, counts))
glm.D93 <- glm(counts ~ outcome + treatment, family = poisson())


model<-SSModel(counts ~ outcome + treatment, data=d.AD,
               distribution = 'poisson')

out<-KFS(model)
coef(out,start=1,end=1)
coef(glm.D93)

summary(glm.D93)$cov.s
out$V[,,1]

outnosim<-KFS(model,smoothing=c('state','signal','mean'))
set.seed(1)
outsim<-KFS(model,smoothing=c('state','signal','mean'),nsim=1000)


## linear
# GLM
glm.D93$linear.predictor
# approximate model, this is the posterior mode of p(theta|y)
c(outnosim$thetahat)
# importance sampling on theta,  gives E(theta|y)
c(outsim$thetahat)



## predictions on response scale
# GLM
fitted(glm.D93)
# approximate model with backtransform, equals GLM
c(fitted(outnosim))
# importance sampling on exp(theta)
fitted(outsim)

# prediction variances on link scale
# GLM
as.numeric(predict(glm.D93,type='link',se.fit=TRUE)$se.fit^2)
# approx, equals to GLM results
c(outnosim$V_theta)
# importance sampling on theta
c(outsim$V_theta)


# prediction variances on response scale
# GLM
as.numeric(predict(glm.D93,type='response',se.fit=TRUE)$se.fit^2)
# approx, equals to GLM results
c(outnosim$V_mu)
# importance sampling on theta
c(outsim$V_mu)

data(sexratio)
model<-SSModel(Male~SSMtrend(1,Q=list(NA)),u=sexratio[,'Total'],data=sexratio,
               distribution='binomial')
fit<-fitSSM(model,inits=-15,method='BFGS',control=list(trace=1,REPORT=1))
fit$model$Q #1.107652e-06

# Computing confidence intervals in response scale
# Uses importance sampling on response scale (4000 samples with antithetics)

pred<-predict(fit$model,type='response',interval='conf',nsim=1000)

ts.plot(cbind(model$y/model$u,pred),col=c(1,2,3,3),lty=c(1,1,2,2))

# Now with sex ratio instead of the probabilities:
imp<-importanceSSM(fit$model,nsim=1000,antithetics=TRUE)
sexratio.smooth<-numeric(length(model$y))
sexratio.ci<-matrix(0,length(model$y),2)
w<-imp$w/sum(imp$w)
for(i in 1:length(model$y)){
 sexr<-exp(imp$sample[i,1,])
 sexratio.smooth[i]<-sum(sexr*w)
 oo<-order(sexr)
 sexratio.ci[i,]<-c(sexr[oo][which.min(abs(cumsum(w[oo]) - 0.05))],
                      + sexr[oo][which.min(abs(cumsum(w[oo]) - 0.95))])
}

# Same by direct transformation:
out<-KFS(fit$model,smoothing='signal',nsim=1000)
sexratio.smooth2 <- exp(out$thetahat)
sexratio.ci2<-exp(c(out$thetahat)
                  + qnorm(0.025) * sqrt(drop(out$V_theta))%o%c(1, -1))

ts.plot(cbind(sexratio.smooth,sexratio.ci,sexratio.smooth2,sexratio.ci2),
        col=c(1,1,1,2,2,2),lty=c(1,2,2,1,2,2))
# Example of Cubic spline smoothing
require(MASS)
data(mcycle)

model<-SSModel(accel~-1+SSMcustom(Z=matrix(c(1,0),1,2),
                                 T=array(diag(2),c(2,2,nrow(mcycle))),
                                 Q=array(0,c(2,2,nrow(mcycle))),
                                 P1inf=diag(2),P1=diag(0,2)),data=mcycle)

model$T[1,2,]<-c(diff(mcycle$times),1)
model$Q[1,1,]<-c(diff(mcycle$times),1)^3/3
model$Q[1,2,]<-model$Q[2,1,]<-c(diff(mcycle$times),1)^2/2
model$Q[2,2,]<-c(diff(mcycle$times),1)


updatefn<-function(pars,model,...){
  model$H[]<-exp(pars[1])
  model$Q[]<-model$Q[]*exp(pars[2])
  model
}

fit<-fitSSM(model,inits=c(4,4),updatefn=updatefn,method="BFGS")

pred<-predict(fit$model,interval="conf",level=0.95)
plot(x=mcycle$times,y=mcycle$accel,pch=19)
lines(x=mcycle$times,y=pred[,1])
lines(x=mcycle$times,y=pred[,2],lty=2)
lines(x=mcycle$times,y=pred[,3],lty=2)