Statistical Software for State Space Methods

In this paper we review the state space approach to time series analysis and establish the notation that is adopted in this special volume of the Journal of Statistical Software . We ﬁrst provide some background on the history of state space methods for the analysis of time series. This is followed by a concise overview of linear Gaussian state space analysis including the modelling framework and appropriate estimation methods. We discuss the important class of unobserved component models which incorporate a trend, a seasonal, a cycle, and ﬁxed explanatory and intervention variables for the univariate and multi-variate analysis of time series. We continue the discussion by presenting methods for the computation of diﬀerent estimates for the unobserved state vector: ﬁltering, prediction, and smoothing. Estimation approaches for the other parameters in the model are also considered. Next, we discuss how the estimation procedures can be used for constructing conﬁdence intervals, detecting outlier observations and structural breaks, and testing model assumptions of residual independence, homoscedasticity, and normality. We then show how ARIMA and ARIMA components models ﬁt in the state space framework to time series analysis. We also provide a basic introduction for non-Gaussian state space models. Finally, we present an overview of the software tools currently available for the analysis of time series with state space methods as they are discussed in the other contributions to this special volume.


Promising trends
In 2001, when considering the possible drawbacks of state space models, Durbin and Koopman wrote: "In our opinion, the only disadvantages are the relative lack in the statistical and econometric communities of information, knowledge, and software regarding these models.",see Durbin and Koopman (2001, p. 52).Ten years later, it is gratifying to see how much progress has been made in the further dissemination of these methods.Not only have state space models been applied in a growing number of scientific fields, but -as is witnessed by this special volume of the Journal of Statistical Software that is completely dedicated to statistical software for state space methods -they have been implemented in STAMP, R, MATLAB, REGCMPNT, SAS, EViews, GAUSS, Stata, RATS, gretl, and SsfPack with links established with S-PLUS and Ox.
State space methods originated in the field of control engineering, starting with the groundbreaking paper of Kalman (1960).They were initially (and still are) deployed for the purpose of accurately tracking the position and velocity of moving objects such as ships, airplanes, missiles, and rockets.A riveting account of this application of state space methods to space travel can be found in the NASA (National Aeronautics and Space Administration) report of McGee and Schmidt (1985), a document that has clearly been written on an old-fashioned type-writer, and can be downloaded as a pdf file from http://ntrs.nasa.gov/.The document explains how state space methods contributed to the successful attempt of the Apollo 11 mission to land Neil Armstrong and Buzz Aldrin on the moon in 1969.
Around the eighties of the last century it was recognized by scientists involved in other fields than control engineering that these ideas could well be applied to time series analysis generally as well.Since then state space methods have been applied in a wide range of subjects, including economics, finance, political science, environmental science, road safety and medicine.
Nowadays state space methods are also used for fitting the autoregressive integrated moving average (ARIMA) models of Box and Jenkins (1976) as these can be put in state space form and then analysed by the Kalman filter.When fitting an ARIMA model to a time series with missing observations in SPSS version 15.0, for example, the SPSS output informs the user that "A Kalman filtering algorithm was used for estimation", confirming the fact that missing values are easily handled in a state space framework while this is much more difficult in the Box-Jenkins approach to time series analysis.
In Section 2 we first provide a general overview of state space models and unobserved component models in particular.State space methods are reviewed in Section 3. Section 4 introduces ARIMA and ARIMA components models, and shows how these can be put in state space form.Generalized or non-Gaussian time series models are discussed in Section 5.In Section 6 we introduce all the software packages currently capable of fitting state space models to time series data.Section 7 concludes.

Linear Gaussian state space models
A time series is a set of observations which are sequentially ordered over time.In a state space analysis the time series observations are assumed to depend linearly on a state vector that is unobserved and is generated by a stochastically time-varying process (a dynamic system).The observations are further assumed to be subject to measurement error that is independent of the state vector.The state vector can be estimated or identified once a sufficient set of observations becomes available.In this section we concentrate on the state space model and its special cases.In Section 3 we discuss methods for estimation, residual analysis and forecasting on the basis of state space models.The expositions rely mostly on the introductory textbook by Commandeur and Koopman (2007) and on the more advanced textbooks by Harvey (1989) and Durbin and Koopman (2001).
The general linear Gaussian state space model for the n-dimensional observation sequence y 1 , . . ., y n will be given in each contribution of this special volume by where α t is the state vector, ε t and η t are disturbance vectors and the system matrices Z t , T t , R t , H t and Q t are fixed and known but a selection of elements may depend on an unknown parameter vector.Equation ( 1) is called the observation or measurement equation, while ( 2) is called the state or transition equation.The p × 1 observation vector y t contains the p observations at time t and the m × 1 state vector α t is unobserved.The p × 1 irregular vector ε t has zero mean and p × p variance matrix H t .
The p × m matrix Z t links the observation vector y t with the unobservable state vector α t and may consist of regression variables.The m × m transition matrix T t in (2) determines the dynamic evolution of the state vector.The r × 1 disturbance vector η t for the state vector update has zero mean and r × r variance matrix Q t .The observation and state disturbances ε t and η t are assumed to be serially independent and independent of each other at all time points.In many standard cases, r = m and matrix R t is the identity matrix I m .In other cases, matrix R t is an m×r selection matrix with r < m.Although matrix R t can be specified freely, it is often composed of a selection from the first r columns of the identity matrix I m .
The initial state vector α 1 is assumed to be generated as independently of the observation and state disturbances ε t and η t .Mean a 1 and variance P 1 can be treated as given and known in almost all stationary processes for the state vector.For nonstationary processes and regression effects in the state vector, the associated elements in the initial mean a 1 can be treated as unknown and need to be estimated.For an extensive discussion of initialization in state space analysis, we refer to Durbin and Koopman (2001, Chapter 5)

Local level model and other unobserved component models
By appropriate choices of the vectors α t , ε t and η t , and of the matrices Z t , T t , H t , R t and Q t , a wide range of different time series models can be derived from ( 1) and ( 2).Here we discuss the class of unobserved components time series models.A number of special cases will be discussed in some detail.Special attention is given to the univariate local level model.
ξ , (all of order 1 × 1) for t = 1, . . ., n, model (1)-(2) reduces to the local level model as given by for t = 1, . . ., n.The level component µ t can be conceived of as the equivalent of the intercept in the classical linear regression model y t = µ + ε t which is obtained by setting all the level disturbances ξ t in (3) equal to zero and with µ = µ 1 .The key difference is that the intercept µ in a regression model is fixed whereas the level component µ t in (3) is allowed to change from time point to time point.
Since the second equation in (3) defines a random walk, the local level model is also referred to as the random walk plus noise model (where the noise refers to the irregular component ε t ).
It can be shown that the dynamic process for x t = y t+1 − y t = η t + ε t+1 − ε t , for t = 1, . . ., n, reduces to the moving average process x t = ε t + θε t−1 where θ relates to the signal-tonoise ratio q = σ 2 ξ / σ 2 ε via a quadratic function.Furthermore, the forecasting function of observations generated by the local level model is equivalent to the exponentially weighted moving average scheme or exponential smoothing.

By defining
, and R t = 1 0 0 1 , the scalar notation of (1) and (2) leads to for t = 1, . . ., n, and we obtain the local linear trend model.
The local linear trend model requires a 2 × 1 state vector α t : one element for the level component µ t and one element for the slope component ν t .The slope component can be conceived of as the equivalent of the regression coefficient in the classical regression model where the observed time series y t is regressed on the independent variable time t: y t = µ + νt + ε t with µ = µ 1 and ν = ν 1 .Again, the important difference is that the regression coefficient or weight ν is fixed in classical linear regression, whereas the slope ν t in the local linear trend model is allowed to change over time.
In the situation that the observed time series consists of quarterly or monthly data, for example, the local level and the local linear trend model can be extended with a stochastic seasonal dummy component denoted here by γ t .In the case of a quarterly time series (the seasonal length is 4), by defining and expanding (1) and (2) in scalar notation, we obtain for t = 1, . . ., n, which is a local level and dummy seasonal model for a quarterly time series where the seasonal component is allowed to change over time.The seasonal dummy model is not the only approach to incorporate time-varying seasonal effects in unobserved components time series model; see Proietti (2000) for a review of different seasonal specifications and their properties.
When a slope component is included in (5) as well, Harvey (1989) refers to this model as the basic structural time series model.A typical application of this model is for the seasonal adjustment of time series.A seasonally adjusted time series is defined in this context simply by the estimate of y t − γ t = µ t + ε t for t = 1, . . ., n.
Another extension is to include one or more cycles to any of the special models within the class of unobserved components time series models.By defining in (1) and ( 2), we obtain the following local level plus cycle model as given by for t = 1, . . ., n, where 0 < ρ ≤ 1 is the damping factor and λ c is the frequency of the cycle measured in radians so that 2π / λ c is the period of the cycle.In case ρ = 1, the cycle reduces to a fixed sine-cosine wave but the component is still stochastic since the initial values c 1 and c * 1 are stochastic variables with mean zero and variance σ 2 c .A typical application of this model is for the signal extraction of business cycles from macro-economic time series.

Regression and intervention effects
Another extension of the local level and local linear trend models concerns the incorporation of fixed explanatory and intervention variables.In the case of one regression variable x t and one intervention variable w t , for example, we have y t = µ t + βx t + λw t + ε t for the local level model and a state vector of three elements is required: one for the level component µ t , one for the regression coefficient β, and one for the regression coefficient λ.The substitution of in ( 1) and ( 2) results in where β = β 1 = β t and λ = λ 1 = λ t for t = 1, . . ., n.This is the local level model with one continuous explanatory variable x t and one intervention variable w t .By adding disturbance terms to the state equation for β t in (7), this regression weight is effectively subjected to a random walk, thus allowing for the estimation of time-varying regression effects (see Harrison and Stevens (1976) for an early application of such time-varying regression effects).
Letting τ denote the time point at which an intervention effect occurred, variable w t can either be coded as a pulse intervention: (to model an outlier observation), or as a level intervention: (to model a structural break in the level of the series), or as a slope intervention: (to model a structural break in the slope of the series).Other types of intervention effects can be modelled as well, see Box and Tiao (1975).

Structural time series analysis
From the previous sections, a structural approach to time series analysis emerges which is facilitated by the state space framework.In this approach, different unobserved components or building blocks responsible for the dynamics of the series such as trend, seasonal, cycle, and the effects of explanatory and intervention variables are identified separately before being put together in a state space model.It is the responsibility of the researcher to decide what components are required in a specific situation, and then to consider whether they apply to the time series under consideration.This explains why state space models are often referred to as structural time series models.The monograph of Harvey (1989) has been instrumental in the dissemination of state space models outside the field of control engineering.

Multivariate models
The structural approach to time series analysis is easily extended to multivariate time series.For example, when we denote y t as a p × 1 vector of observations, a multivariate local level model can be adopted for modelling the p time series simultaneously: for t = 1, . . ., n, where µ t , ε t , and ξ t are p×1 vectors and Σ ε and Σ ξ are p×p variance matrices.
In what is known as the seemingly unrelated time series equations model (8), the series are modelled as in the univariate situation, but the disturbances driving the level components are allowed to be instantaneously correlated across the p series.When slope, seasonal, or cycle components are also included in the model, each of these three components has an associated p × p variance matrix for their disturbance vectors, allowing for correlation between the series.
If the rank r of Σ ξ in ( 8) is assumed smaller than p, then the model implies that the p series have r common levels.Such common factors may not only have a nice and interesting interpretation, but may also result in more efficient inferences and forecasts since the number of parameters reduces.

State space analysis
For given values of all system matrices -and for given initial conditions a 1 and P 1 -the state vector can be estimated in three different ways, yielding what are known as the filtered, the predicted, and the smoothed estimates of the state vector.Depending on the types of state estimates required in the analysis, the estimates of the state vector can be obtained by performing one or two passes through the observed time series: 1. a forward pass, from t = 1, . . ., n, using a recursive algorithm known as the Kalman filter enables the computation of predicted states, based on y 1 , . . ., y t−1 , filtered states, based on y 1 , . . ., y t , and observation prediction errors; 2. a backward pass, from t = n, . . ., 1, using output of the Kalman filter and using recursive algorithms known as state and disturbance smoothers enables the computation of smoothed estimates of states and disturbances that are based on all observations y 1 , . . ., y n .

Kalman filter for prediction, filtering and forecasting
The forward pass through the data with the well-known Kalman (1960) filter provides all estimates that are relevant for the filtered and the predicted state.The main purpose of the

Statistical Software for State Space Methods
Kalman filter is to obtain optimal estimates of the state at time point t, only considering the observations {y 1 , y 2 , . . ., y t−1 }.A key property of the predicted state and its related estimates is therefore that they are only based on past values of the observed time series.The recursive formulas for the Kalman filter are for t = 1, . . ., n.The values of a t in (9) represent the predicted state, while the values of P t quantify the estimation error variance matrix of the predicted state a t .Under the assumption of normality, the latter variances are useful for the construction of confidence intervals for the predicted state, which -assuming that we are interested in their 90% confidence limits for example -can be calculated as a t ± 1.64 P t , for t = 1, . . ., n.A modification of the Kalman filter also allows the computation of the filtered estimate of the state vectors, that is where a t|t is the optimal estimate of the state at time point t by considering the observations {y 1 , y 2 , . . ., y t } while P t|t is the state filtered estimation error variance matrix.The values of v t in ( 9) are called the one-step ahead prediction or forecast errors, since they quantify the lack of accuracy of a t in predicting the observed value of y t at time point t; the values of F t are the variances of these one-step ahead prediction errors v t .
One of the convenient features of state space methods is the ease with which they deal with two important aspects of time series analysis -forecasting and missing observations: by treating them in exactly the same way.Missing observations are handled by setting K t and v t in (9) equal to 0. Forecasts for y n+1 , . . ., y n+k given y 1 , . . ., y n are simply obtained by applying the Kalman filter for t = 1, . . ., n, n + 1, . . ., n + k and by treating y n+1 , . . ., y n+k as missing observations.

State and disturbance smoothing
The backward pass through the data is only required for smoothing that leads to estimates such as the smoothed states and smoothed disturbances.The main purpose of state and disturbance smoothing is to obtain estimated values of the state and disturbance vectors at time point t, considering all available observations {y 1 , y 2 , . . ., y n }.The recursive formulas for state smoothing are for t = n, . . ., 1.The recursive formulas for smoothing (10) are initialized with r n = 0 and N n = 0.The state smoothing equations (11) yield the smoothed state estimate αt and is defined as the optimal estimate of α t using the full set of observations {y 1 , y 2 , . . ., y n }; the state smoothing equations also yield the corresponding smoothed state estimation error variance matrix V t .
Analogous to the predicted state, under the assumption of normality the smoothed state estimation error variance matrix V t is useful for the construction of confidence intervals for the smoothed state components, which -should we happen to be interested in their 90% confidence limits for example -can be calculated as αt ± 1.64 V t , for t = 1, . . ., n.
The recursions for r t−1 and N t−1 in (10) also enable the computation of the smoothed estimates of the disturbances ε t and η t in the following way, for t = n, . . ., 1.The equations ( 12) and ( 13) compute the smoothed observation disturbances εt , the smoothed state disturbances ηt , and their corresponding smoothed estimation error variance matrices Var(ε t ) and Var(η t ).

Diagnostic checking
All significance tests in linear Gaussian state space models -and the construction of confidence intervals -are based on three assumptions concerning the residuals of the analysis.
The residuals should satisfy independence, homoscedasticity, and normality, in this order of importance.Whether the residuals satisfy these three assumptions can be established by diagnosing what are known as the standardized prediction errors.They are defined as for t = 1, . . ., n.For the computations of the one step-ahead prediction errors v t and their variances F t in ( 14), we refer to the recursive formulas for the Kalman filter given in (9).In case y t is a vector of time series, the prediction error v * t is defined by t .The assumptions of independence and normality of the residuals can be diagnosed using the Box-Ljung test statistic and the Bowman and Shenton test statistic, respectively.The assumption of homoscedasticity can be checked by testing whether the variance of the standardized prediction errors in the first third part of the series is equal to the variance of the errors corresponding to the last third part of the series.For further details concerning these diagnostic tests, we refer to Harvey (1989), Durbin and Koopman (2001) and Commandeur and Koopman (2007).
A second diagnostic tool for determining the appropriateness of a model is provided by inspection of what are known as the auxiliary residuals.As already mentioned above, the disturbance smoothing filters applied in the backward pass through the data yield, amongst others, estimates of the smoothed observation and state disturbances, and of their variances.
The auxiliary residuals are obtained by dividing the smoothed observation and state disturbances with the square root of their corresponding variances, as follows: for t = 1, . . ., n, resulting in standardized smoothed disturbances.Visual inspection of the standardized smoothed observation disturbances e * t allows for the detection of possible outlier observations in a time series, while inspection of the standardized smoothed state disturbances r * t enables the detection of structural breaks in the underlying development of a time series.Each auxiliary residuals can be considered as a t test for the null hypothesis that there was no outlier observation (when inspecting the auxiliary residuals at the left of ( 15)) or as a t test for the null hypothesis that there was no structural break in the corresponding unobserved component of the observed time series (when inspecting the auxiliary residuals at the right of ( 15)).Applying the usual 95% confidence limits of ±1.96 corresponding to a two-tailed t test, possible outlier observations or structural breaks in the unobserved components making up the state vector are thus easily detected.A more detailed discussion on outlier and break detection in a state space analysis is provided by Harvey and Koopman (1992) and de Jong and Penzer (1998).

Parameter estimation
So far, we have presented all of the results that can be obtained with state space methods as if the disturbance variances, the fixed regression effects, the parameters ρ and λ c associated with cycles, etcetera, are given and known.In practice, of course, these parameters are unknown, and have to be estimated.
It can be shown that the Kalman filter presented in ( 9) also provides the necessary ingredients required for evaluating the log-likelihood function, which is given by log where the v t are the one-step ahead prediction errors, the F t are their variances for t = 1, . . ., n defined in (9), and ψ denotes the vector of unknown parameters.The log-likelihood ( 16) is maximized with respect to ψ numerically using the score vector or the EM algorithm, see Durbin and Koopman (2001) and Commandeur and Koopman (2007) for details.In this way we obtain the maximum likelihood estimate of ψ.

ARIMA and ARIMA components models
In this section we consider more traditional approaches to the analysis of time series proposed in Box and Jenkins (1976).We will show how the autoregressive integrated moving average (ARIMA) model can be put into state space form.An essential feature of the Box-Jenkins approach to time series analysis is that the observed times series is assumed to be stationary, meaning that its means and covariances should be invariant when the series is shifted or translated through time.If an observed time series satisfies stationarity, then the corresponding series can be analysed with an autoregressive moving average (ARMA) model.However, if the observed time series also contains non-stationary features like a trend and a seasonal, for example, then stationarity has to be enforced before the analysis can start by differencing the observed time series: key features of the series such as a trend and a seasonal will be removed from the series y t by transforming it into ∆y t = y t − y t−1 , to remove the trend in the series, and to remove a seasonal with periodicity s from the series.In some cases a combined removal of trend and seasonal is necessary and is achieved by Should the differenced time series still not be stationary, then the differencing procedure can be continued by taking second differences such as or even combinations with third differences, etcetera.
We denote the stationary time series by y * t which can be equal to the series y t but can also be the series y t after it is appropriately differenced.The Box-Jenkins ARMA(p, q) model for the analysis of univariate time series is given by where p and q are non-negative integers, and ζ t is a series of independent disturbances (also referred to as white noise) that is normally distributed here.In ( 17) the values of p and q denote the number of autoregressive and moving average terms required to model the observed series, respectively.If the observed series already satisfies stationarity, then y * t = y t in (17), and the corresponding model is called an ARMA model.Should this not be the case, then y * t is obtained after some differencing operations and we are dealing with an ARIMA model.
To illustrate how an ARMA model can be formulated in state space form, we consider the ARMA(2, 1) model for t = 1, . . ., n.By considering the state space model ( 1)-( 2) and by defining for t = 1, . . ., n, we have represented the ARMA(2, 1) model in the state space framework.For a general overview of how any of the general ARMA (p, q) models can be represented in the state space framework, we refer to Durbin and Koopman (2001, Section 3.3).
A more general approach to time series analysis can be based on the ARIMA components model as given by for t = 1, . . ., n, where the µ (j) t 's follow mutually independent ARIMA processes with different levels of differencing and different values for p and q for j = 1, . . ., m.The dynamic properties of y t can be derived from the different ARIMA processes while the time series y t can be expressed as an ARIMA process itself.Unobserved components time series models can be regarded as special cases of the ARIMA components model.For example, the local level model (3) can be expressed as an ARIMA components model with m = 2, µ (1) t = µ t and µ (2) t = ε t .We take µ (1) t as an ARMA(0, 0) process after first differencing (with p = 0 and q = 0) and we take µ (2) t as an ARMA(0, 0) process without differencing.A more general and more detailed discussion of ARIMA components models can be found in Bell (2004).

Non-Gaussian state space models
The state space model ( 1)-( 2 The generalized state space model ( 19)-( 2)-( 20) is more general than the exponential family only.We can consider continuous density functions for p(•) such as those related to the heavytailed distributions Student's t, mixture of normals and general error.In these cases, we can often return to the linear observation equation ( 19) but written as where p(ε t ) = p(y t |θ t ).When we take p(•) as the Gaussian density with mean zero and variance H t , the generalized model reduces to the Gaussian state space model ( 1)-(2).
A particular example of a non-Gaussian model of interest for financial markets is the stochastic volatility model where the observations are financial returns with a constant mean and a stochastically time-varying variance.The univariate stochastic volatility model is given by where p(ε t ) can be Gaussian, Student's t or mixture of normals.The log-variance θ t can be generally specified as θ t = Z t α t but it is usually considered as a stationary autoregressive process.For more detailed discussions concerning the stochastic volatility model ( 21), see Harvey, Ruiz, and Shephard (1994) and Shephard (2005).
The statistical analysis based on generalized state space models cannot rely on linear methods such as the Kalman filter and smoothing recursions of Section 3. A number of statistical approaches can be adopted to treat non-Gaussian and nonlinear features of the model adequately.Durbin and Koopman (2000) have developed an estimation methodology based on importance sampling and simulated maximum likelihood methods.Bayesian estimation methods for generalized state space models rely typically on Markov chain Monte Carlo (MCMC) methods and are developed by Frühwirth-Schnatter (1994), Carter and Kohn (1994) and Shephard and Pitt (1997).
More recent developments on the analysis of generalized or nonlinear state space models are reported in the engineering literature.These developments are based on methods collectively known as particle filtering and originate from the work by Gordon, Salmond, and Smith (1993), Kitagawa (1996) and with a noteworthy contribution by Pitt and Shephard (1999).Particle filtering and related Monte Carlo methods are still in development as can be learned from the recent review article of Creal (2011).

Software packages
An overview of all the software packages covered in this special volume is provided in

Structure of the papers
When inviting the authors to contribute to this special volume, in order to obtain some degree of uniformity we asked them all to adhere to the following structure in their papers: Title: State space methods in <name of software package> Author: <author>

Conclusions
We have reviewed state space methods for linear Gaussian state space models and discussed their generalizations for models with nonlinear non-Gaussian features.In this special volume of the Journal of Statistical Software twelve contributions are given that discuss the currently available software implementations of state space methods and their extensions.The most important aim of this volume is to highlight the presence of such implementations for different statistical software platforms.We hope it will encourage more applied researchers to use these software tools in order to improve the quality of their time series analyses.
(including possibilities and limitations of the software at hand) For each case study in their contribution, we asked them to present and discuss the (important parts of the) required code and the analysis results.Moreover, for case study 1 we asked the contributors to show how to apply the local level model to the Nile data (a series of readings of the annual flow volume at Aswan from 1871 to 1970 also used as an illustration inDurbin and Koopman 2001) with the software at hand, and when discussing this analysis in their paper we also asked them to include the following results: the ML estimates of the two variances including asymptotic standard errors the plot with data and smoothed state vector with 90% CI the plot with standardized prediction errors the two plots with auxiliary residuals and 95% CI a plot with forecasts until and including 1980 with 50% CI After discussing the code and results of the analysis of the Nile data with the local level model in Section 2 of their paper, we finally asked them to add one or two examples of their own design.Of these added examples we mention the following highlights, including references to the corresponding papers in this volume: illustrations on ocean and climate time series, specifically El Niño and a space-time data set on sea surface temperature (with STAMP in Mendelssohn 2011); a basic structural model with cycle for the Italian industrial production and a stochastic volatility model for the FTSE100 data (with Ox/SsfPack in Pelagatti 2011); a Bayesian analysis of the Nile data and a multivariate dynamic capital asset pricing model for four assets (with R in Petris and Petrone 2011); an affine model for the term structure of interest rates (with S-PLUS in Zivot 2011) a cubic spline analysis for the motorcycle acceleration data of Silverman (1985) and some non-Gaussian illustrations (with MATLAB in Peng and Aston 2011); a time series model with a sampling error component (with REGCMPNT in Bell 2011); a bivariate latent risk model for the simultaneous analysis of mobility and road traffic fatalities (with EViews in Van den Bossche 2011); an illustration of the dynamic stochastic general equilibrium (DSGE, with RATS in Doan 2011); a vector autoregressive moving average (VARMA) model, and a dynamic factor model (with Stata in Drukker and Gates 2011); a multivariate structural time series model for the labour market applied to the number of employees and the number of hours worked (with gretl in Lucchetti 2011); the use of diffuse priors in the state vector with an application for a two-way random effects panel model (with SAS in Selukar 2011); a Bayesian analysis of time varying volatility for Standard & Poors 500 returns (with Ox in Bos 2011).
) is for Gaussian observations y t .In case the observations y t are discrete or have other distributional properties than Gaussian, we can consider the generalized state space modely t ∼ p(y t |Z t α t ),(19)for t = 1, ..., n, where p(•) is a particular non-Gaussian density function, possibly for discrete observations, and with the state vector α t and system matrix Z t as defined in the previous section.Conditional on the state α t , we assume that the observations are serially independent.This assumption is similar to the linear Gaussian model.The state vector α t evolves stochastically over time according to the linear Gaussian state equation (2).Time series observations which come from the exponential family distribution require density function p(yt |θ t ) = exp y t θ t − b t (θ t ) + c t (y t ) ,(20)whereθ t = Z t α t , b t (•)is a twice differentiable function and c t (•) is a function of y t only.The modelling framework (19)-(2)-(20) makes it possible to analyse time series observations coming from, among others, Poisson, binary, binomial and multinomial distributions.For example, in case a Poisson density is considered, we have b t (θ t ) = exp θ t and c t (y t ) = log(y t !).

Table 1 ,
including their current version number, and indicating whether they currently can handle univariate unobserved component models, multivariate unobserved component models, the univariate treatment of multivariate time series presented inDurbin and Koopman (2001,  Chapter 6), the exact initialization procedures described inDurbin and Koopman (2001,  Chapter 5), and non-Gaussian state space models.

Table 1 :
Overview of state space software packages and their options.The options are univariate models (UM), multivariate models (MM), exact initialization (EI), univariate treatment of multivariate time series (UTMTS), and non-linear and non-Gaussian models (NLNGM).