PaperlayerView exact revision on arXiv

Source-linked AI summary

KFAS: Exponential Family State Space Models in R

Jouni Helske

arXiv:1612.01907v2stat.COstat.ME

TL;DR

State space modelling covers diverse time-series problems, while non-Gaussian state space modelling in R has been described as somewhat experimental. This paper presents a rewritten KFAS package for exponential-family observations, reports expanded non-Gaussian functionality and a Poisson forecasting example, and notes limits on non-Gaussian filtering efficiency.

  • Problem

    Non-Gaussian modelling in earlier KFAS versions was somewhat experimental, motivating a fully functional package for these observations.

    Source: §1 · p. 2

  • Method

    KFAS provides state space modelling and inference tools for exponential-family observations, with symbolic model descriptions and auxiliary functions for constructing models.

    Source: §1 · p. 1, §9 · p. 33

  • Results

    The rewritten KFAS supports fully functional multivariate non-Gaussian models with different distributions, and its Poisson example estimates a 0.022 slope corresponding to a 2.3% yearly increase in deaths.

    Source: §1 · p. 2, §3.2 · p. 10

  • Takeaways & Limitations

    KFAS offers tools for maximum likelihood estimation, filtering, smoothing, and simulation, while symbolic descriptions can ease construction of models with simple components.

    Source: §9 · p. 33

    Non-Gaussian filtering is relatively slow for large models because importance sampling is performed n times; the paper recommends particle filtering for large-model or online-filtering problems.

    Source: §3 · p. 8

Abstract

from arXiv · show

State space modelling is an efficient and flexible method for statistical inference of a broad class of time series and other data. This paper describes an R package KFAS for state space modelling with the observations from an exponential family, namely Gaussian, Poisson, binomial, negative binomial and gamma distributions. After introducing the basic theory behind Gaussian and non-Gaussian state space models, an illustrative example of Poisson time series forecasting is provided. Finally, a comparison to alternative R packages suitable for non-Gaussian time series modelling is presented.

1. Introduction

State space models provide a unified framework for diverse time-series models, while KFAS extends Gaussian modelling to exponential-family observations. The paper presents a rewritten, more user-friendly package with functional non-Gaussian modelling and additional inference features.

1. Introduction

  • State space models can represent structural time series, ARIMA models, simple regression, generalized linear mixed models, and cubic spline smoothing.Source: §1 · p. 1
  • KFAS supports Gaussian state space models and observations from binomial, Poisson, negative binomial, and Gamma distributions.Source: §1 · p. 1
  • The rewritten package makes model definition more user-friendly through R symbolic formulas.Source: §1 · p. 1
  • Non-Gaussian modelling is now fully functional for multivariate models with different distributions, alongside added residual and prediction methods and improved main-function performance.The rewrite also fixed several bugs.Source: §1 · p. 2
  • The paper introduces the theory, demonstrates KFAS on real data, and compares it with packages for non-Gaussian time series modelling.Source: §1 · p. 2

2. Gaussian state space model

The section introduces Gaussian state space modelling and its filtering and smoothing framework, then describes diffuse initialization and univariate processing. It closes with the diffuse likelihood and an alcohol-deaths example.

2. Gaussian state space model

  • Gaussian state space models specify Gaussian observation and state disturbances and an initial Gaussian state; filtering yields one-step-ahead predictions and smoothing estimates states using the full series.Source: §2 · pp. 2–3
  • Exact diffuse initialization represents diffuse state uncertainty separately and uses special recursions until the diffuse covariance component vanishes.It is less prone to numerical errors than setting P1 = κI, although smoothing can still be numerically unstable with high collinearity.Source: §2 · p. 3
  • Univariate processing addresses non-invertible or computationally costly observation covariance matrices in multivariate filtering, and simplifies diffuse-phase recursions.KFAS uses this approach; its computational advantage is not guaranteed because additional linear-algebra calls and memory management can offset gains.Source: §2 · pp. 3–4

2.1. Log-likelihood of the Gaussian state space model

  • For diffuse initialization, the likelihood includes diffuse-phase terms as well as the usual prediction-error contributions.When Zt or Tt contain unknown parameters in diffuse parts, the standard diffuse likelihood can omit a parameter-dependent term; KFAS supports the recursive correction described by Francke et al. (2010, p.411–412).Source: §2.1 · p. 4

2.2. Example of Gaussian state space model

  • In the Finland alcohol-deaths example, filtering recursively updates the constant slope estimate, while smoothing uses past and future observations to estimate the time-varying level.The figure contrasts observed deaths with predicted and smoothed estimates.Source: §2.2 · p. 6, Figure 1 · p. 6

3. State space models for the exponential family

KFAS extends state-space models to exponential-family observations and uses Gaussian approximations and importance sampling for non-Gaussian inference.

3. State space models for the exponential family

  • The observation distributions include Gaussian, Poisson, binomial, Gamma, and negative binomial, with the state equation retained from the Gaussian case.Source: §3 · p. 7
  • Non-Gaussian inference iteratively uses a Gaussian model with the same conditional posterior mode, then importance sampling to correct approximation errors and estimate smoothed state functions.The Gaussian approximation’s pseudo-observation variances are based on derivatives of the observation log-density.Source: §3 · pp. 7–8
  • Non-Gaussian filtering repeatedly applies smoothing to data prefixes and is relatively slow for large models; particle filtering is recommended for large or online problems.Source: §3 · p. 8
  • Prediction intervals are formed by importance-weighted predictor draws, sampling observations from the resulting predictors, and taking empirical quantiles.Source: §3 · p. 8

3.1. Log-likelihood of the non-Gaussian state space model

  • The non-Gaussian log-likelihood uses the Gaussian approximating model and an expectation under its conditional state density; often, the mode alone gives a good approximation.Raw importance weights can be numerically unstable when many observations make the discrete likelihood tend to zero.Source: §3.1 · p. 9

3.2. Example of non-Gaussian state space model

  • In the Poisson model, the slope was estimated at 0.022 with standard error 1.4 × 10−4, corresponding to a 2.3% yearly increase in deaths.The η estimates for the Gaussian and Poisson models are not directly comparable because µt has a different interpretation in each model.Source: §3.2 · pp. 9–10

4. Residuals

KFAS provides several residual diagnostics based on Kalman-filter predictions or smoothed states. Multivariate standardization choices and the diffuse phase affect how some residuals should be interpreted.

4. Residuals

  • Recursive residuals use one-step-ahead Kalman-filter predictions, but are defined only after the diffuse phase.A long diffuse phase can make them problematic; a proper, highly non-informative prior for an intervention variable can be used when computing residuals.Source: §4 · pp. 10–11
  • Marginal standardization aims for residual series without autocorrelation; Cholesky standardization aims for the same property in the residual series treated jointly.The Cholesky option uses the prediction-error covariance matrix and can be computed directly from the sequential Kalman filter.Source: §4 · p. 11
  • Multivariate residuals can depend on series ordering, whereas quadratic-form residuals do not.Marginal residuals may be preferred when ordering complicates interpretation of diagnostic deviations.Source: §4 · p. 11
  • Auxiliary residuals provide another option, using smoothed state estimates.The paper refers readers to Harvey and Koopman (1992) and Durbin and Koopman (2012, Chapter 7) for details.Source: §4 · p. 11

5. Functionality of KFAS

KFAS separates model construction, parameter estimation, and state inference, with additional tools for simulation, prediction, and diagnostics. Its functions support both Gaussian and non-Gaussian state space models.

5. Functionality of KFAS

  • SSModel builds a model from a formula and auxiliary components, while fitSSM estimates unknown parameters using an update function.By default, fitSSM estimates NA values in time-invariant H and Q; users can supply their own optimization approach through logLik.Source: §5 · pp. 11–12
  • KFS returns filtered and smoothed state and signal estimates, including inverse-link estimates for non-Gaussian models.For Gaussian models, it also provides disturbance smoothing.Source: §5 · p. 12
  • simulateSSM simulates Gaussian-model states, signals, or disturbances, while importanceSSM returns non-Gaussian states or signals from an approximating Gaussian model with weights.simulateSSM can also simulate missing observations and states from predictive distributions.Source: §5 · p. 12
  • KFAS supplies prediction intervals and residual-inspection tools, alongside methods for extracting model signals, coefficients, fitted values, and residual diagnostics.The extraction operator for SSModel subsets is suggested over the ordinary list extractor to avoid accidentally changing model-matrix dimensions.Source: §5 · p. 12

6. Constructing common state space models with KFAS

KFAS uses auxiliary functions or direct system-matrix definitions to construct common state space models, including structural, ARIMA, regression, mixed, and custom components.

6.1. Structural time series

  • Structural models decompose a series into trend, seasonal, and cycle components, built with SSMtrend, SSMseasonal, and SSMcycle.Seasonality can be dummy-variable or trigonometric; non-Gaussian structural models can also include Gaussian noise through SSMcustom.Source: §6.1 · pp. 13–14

6.2. ARIMA models

  • SSMarima defines ARIMA models through autoregressive and moving-average coefficients, differencing, and a stationarity option; the example's drift and log-likelihood match the structural model.The example estimates θ1 = −0.5 and σ = 17.Source: §6.2 · pp. 15–16

6.3. Linear and generalized linear models

  • State space GLMs extend classical regression to multivariate and dynamic settings, including random-walk coefficients and time-varying dispersion.For an ordinary GLM, the final coefficient estimates and standard errors coincide with those from ordinary GLM fitting.Source: §6.3 · pp. 17–18

6.4. Generalized linear mixed models

  • KFAS mixed-model estimates of likelihood, variance and covariance parameters, and fixed and random effects are practically identical to lme4 estimates.The random-effect covariance standard errors differ because KFAS accounts for fixed-effect estimation uncertainty, while lme4 conditions on the fixed effects.Source: §6.4 · p. 20

6.5. Arbitrary state space models

  • The Gaussian and Poisson structural models with added white noise yield nearly indistinguishable smoothed trend estimates.The paper attributes this similarity to the Poisson process's relatively high intensity.Source: Figure 3 · p. 21, §6.5 · p. 21

7. Illustration

The illustration fits a multivariate Poisson state-space model to Finnish alcohol-related deaths across four age groups, then examines smoothed states, residuals, and forecasts. Forecasts for 2008–2013 differed from subsequent observed trends across age groups.

7. Illustration

  • The example models four age groups using data through 2007 and forecasts deaths for 2008–2013.The modeled series covers ages 30–39, 40–49, 50–59, and 60–69.Source: §7 · pp. 21–22, Figure 4 · p. 22
  • The multivariate Poisson model combines a random-walk level with drift and an additional white-noise component for extra variation.The level and noise covariance structures are left unrestricted.Source: §7 · p. 23
  • Estimates from the initial fit and the subsequent importance-sampling step were practically identical.The procedure first fits without simulation, then uses those estimates to initialize the importance-sampling fit.Source: §7 · p. 23
  • Parameter estimation is sensitive to initial values, and iterative non-Gaussian likelihood noise can affect BFGS gradient calculations.The paper recommends trying multiple initial configurations; BFGS is faster than Nelder–Mead and is preferred for preliminary analysis.Source: §7 · p. 24
  • Smoothed estimates show clear trends in three age groups, while the slope is not significant for ages 30–39.The additional variation is highly correlated across all four groups.Source: §7 · p. 25, Figure 5 · p. 26
  • Residual checks found occasional lagged cross-correlation, and the noise component’s large 1973 drop may reflect an outlier.The authors nevertheless describe the model as relatively satisfactory overall.Source: §7 · p. 25, Figure 6 · p. 27
  • Comparing 2008–2013 forecasts with true observations, alcohol-related deaths slightly increased among ages 60–69, while deaths at another age decreased substantially during the forecasting period.This was partly explained by an almost monotonic decline in total alcohol consumption, which might have been caused by alcohol tax increases in 2008, 2009 and 2012.Source: §7 · p. 27

8. Other packages for non-Gaussian time series modelling

Non-Gaussian time-series packages differ in their modelling approaches and supported structures, so package choice depends on the problem. In a Salmonella example, KFAS and INLA produced practically the same results despite their different frameworks.

8. Other packages for non-Gaussian time series modelling

  • KFAS and pomp use parameter-driven models with latent states evolving as stochastic processes, while observation-driven models can yield closed-form likelihood evaluation and computational gains for non-Gaussian models.The comparison notes that pomp applies to a broader model class, but has a relatively high learning curve because users write their own simulation and likelihood functions.Source: §8 · p. 29
  • acp, glarma, and tscount implement observation-driven count or time-series regression approaches, and all three assume univariate responses.acp models Poisson counts; glarma supports Poisson, binomial, and negative binomial distributions, while tscount supports Poisson and negative binomial distributions.Source: §8 · p. 29
  • Package choice depends on whether time-dependence is a nuisance or the latent time-varying processes are of interest; KFAS and pomp support flexible latent structures such as stochastic trends and seasonal and cyclic components.The authors note that glarma and tscount offer more limited options because of their general model specification.Source: §8 · p. 30
  • INLA, KFAS, and pomp can incorporate time-varying regression coefficients and random effects and handle missing response observations straightforwardly.The paper contrasts these packages with others that may omit missing time points without accounting for uneven spacing during estimation.Source: §8 · p. 30

8.1. Comparison to INLA

  • For the Salmonella random intercept, KFAS offers two equivalent constructions, but the 18-state SSMregression option is less efficient on large data than a single-state SSMcustom construction.The example uses SSMregression; the two constructions give identical results.Source: §8.1 · pp. 30–31
  • KFAS and INLA gave practically the same results on the Salmonella data, although KFAS estimated σ2 = 0.066 and INLA estimated σ2 = 0.048.Substituting INLA’s σ2 estimate into KFAS produced a slightly lower log-likelihood (-73.50 versus -73.66).Source: §8.1 · p. 32
  • INLA’s Bayesian framework accounts for parameter uncertainty in inference, but numerical integration over hyperparameters can become infeasible as their number grows.The paper notes that defining good priors for those hyperparameters is also non-trivial in practice.Source: §8.1 · p. 33

9. Discussion

KFAS provides tools for inference in linear state space models with exponential-family observations. Its auxiliary functions and formula-based symbolic descriptions let users construct simple components without explicitly defining system matrices.

9. Discussion

  • KFAS provides relatively fast tools for maximum likelihood estimation, filtering, smoothing, and simulation in linear state space models with exponential-family observations.The package is intended to support multiple inference problems across this general modelling framework.Source: §9 · p. 33
  • For simple model components, KFAS users can construct state space models with auxiliary functions and formula-based symbolic descriptions instead of explicitly defining system matrices.The paper says this should greatly ease package use compared with early KFAS versions.Source: §9 · p. 33

A. Appendix: Filtering and smoothing recursions

The appendix collects the Kalman filtering and smoothing formulas used in KFAS for diffuse and sequential cases, drawing on Durbin and Koopman (2012) and related articles.

A. Appendix: Filtering and smoothing recursions

  • The appendix gathers KFAS’s diffuse and sequential Kalman filtering and smoothing formulas, based on Durbin and Koopman (2012) and related articles.The author collected formulas that were scattered across references and used slightly different notation.Source: A. Appendix: Filtering and smoothing… · p. 34

A.1. Filtering

The filtering appendix presents recursions for the general Gaussian model, including a univariate approach and diffuse filtering.

A.1. Filtering

  • The appendix presents Kalman filter recursions for the general Gaussian model and separately gives the univariate filtering equations.It specifies that the univariate quantities are indexed by time t and observation i, with scalar prediction errors and variances.Source: §A.1 · p. 34
  • The diffuse univariate filtering equations include transition equations from time t to t + 1.When F∞,t,i = 0, the passage states that the transition equations apply.Source: §A.1 · pp. 34–35

A.2. Smoothing

KFAS smoothing is built from backward recursions that produce state and disturbance smoothing recursions. A separate set of recursions handles the diffuse phase.

A.2. Smoothing

  • The smoothing algorithms are based on recursions that yield state and disturbance smoothing recursions.Source: §A.2 · p. 35
  • The diffuse phase has its own recursions for state and disturbance smoothing.Source: §A.2 · pp. 35–36
  • During diffuse-phase transition from time t to t − 1, the recursion uses N and Tᵀ with the specified terminal diffuse-phase conditions.Source: §A.2 · p. 36
Loading 1612.01907v2…