SSMMATLAB : A Set of MATLAB Programs for the Statistical Analysis of State Space Models

This article discusses and describes SSMMATLAB , a set of programs written by the author in MATLAB for the statistical analysis of state space models. The state space model considered is very general. It may have univariate or multivariate observations, time-varying system matrices, exogenous inputs, regression eﬀects, incompletely speciﬁed initial conditions, such as those that arise with cointegrated VARMA models, and missing values. There are functions to put frequently used models, such as multiplicative VARMA models, VARMAX models in echelon form, cointegrated VARMA models, and univariate structural or ARIMA model-based unobserved components models, into state space form. There are also functions to implement the Hillmer-Tiao canonical decomposition and the smooth trend and cycle estimation proposed by G´omez (2001). Once the model is in state space form, other functions can be used for likelihood evaluation, model estimation, forecasting and smoothing. A set of examples is presented in the SSMMATLAB manual to illustrate the use of these functions.


Introduction
This article describes SSMMATLAB (Gómez 2014), a set of programs written by the author in MATLAB (The MathWorks Inc. 2014) for the statistical analysis of time series that are assumed to follow state space models. The series can be univariate or multivariate and the state space model can be very general. It may have time-varying system matrices, exogenous inputs, regression effects, incompletely specified initial conditions, such as those that arise with cointegrated VARMA (vector autoregressive moving average) models, and missing values.
The motivation for SSMMATLAB is to provide the time series analyst with a set of programs written in MATLAB that will allow him to work with general state space models. Since many time series models can be put into state space form, special functions have been written for the most usual ones, such as multiplicative VARMA models, VARMAX models in echelon form, cointegrated VARMA models, univariate structural models, like those considered by Harvey (1989, Chapter 4) or Kitagawa and Gersch (1996), and ARIMA model-based (AMB) unobserved components models (Gómez and Maravall 2001b, Chapter 8). But if the user intends to work with more sophisticated state space models that are not available in standard commercial packages for time series analysis or econometrics, he can program his own model in SSMMATLAB and carry out model estimation, interpolation, forecasting and smoothing.
SSMMATLAB provides functions similar to the ones contained in previous packages for linear state space models. In addition, it provides functions for identification, estimation, forecasting and smoothing of VARMAX models, possibly in state space echelon form, and of cointegrated VARMA models. It provides also functions to design digital filters and to estimate smooth trends and cycles in an AMB approach. Moreover, the general functions in SSMMATLAB allow, with careful programming, to do at least all the things that the previous packages can do with linear state space models.
In Section 2, the state space model will be described. In Section 3, the functions to put into state space form multiplicative VARMA models, VARMAX models in echelon form, univariate structural models and AMB unobserved components models will be documented. In Section 4, the identification of VARMAX(p, q, r) models and VARMAX models in echelon form will be considered. Also in Section 4, the estimation of VARX models, the Hannan and Rissanen (1982) method to estimate VARMAX models, as well as the conditional and the exact methods to estimate VARMAX models will be described. The functions for likelihood evaluation, computation of recursive residuals, model estimation, forecasting and smoothing will be described in Section 5. Finally, in Section 6, reference will be made to some examples and case studies using SSMMATLAB.

The state space model
The state space model considered in SSMMATLAB is Y t = X t β + Z t α t + G t t , t = 1, . . . , n, where {Y t } is a multivariate process with Y t ∈ R p , W t , T t , H t , X t , Z t and G t are time-varying deterministic matrices, β ∈ R q is a constant bias vector, α t ∈ R r is the state vector, and { t } is a sequence of uncorrelated stochastic vectors, t ∈ R s , with zero mean and common covariance matrix σ 2 I. The initial state vector α 1 is specified as where c has zero mean and covariance matrix σ 2 Ω, a is a constant vector, A is a constant matrix, and δ has zero mean and covariance matrix kI with k → ∞ (diffuse). It is assumed that the vectors c and δ are mutually orthogonal and that α 1 is orthogonal to the { t } sequence. The vector δ in (3) models uncertainty with respect to the initial conditions. For example, a multivariate random walk model, Y t = Y t−1 + A t , where {A t } is a zero mean normally distributed sequence with common covariance matrix Σ, can be put into state space form as where Σ = LL is the Cholesky decomposition of Σ, L t = A t+1 , σ 2 = 1, and α 1 = δ.
The state space model (1) and (2) is very general. For example, it can be used in macroeconomics for analyzing time-varying parameter VARs as in Primiceri (2005) as well as for forecasting using mixed frequency data as in Aruoba, Diebold, and Scotti (2009).
It is not restrictive that the same term, t , appears in both equations. To see this, suppose the state space model where W t , T t , J t , X t and Z t are time-varying deterministic matrices, δ ts denotes the Kronecker delta, u t ∈ R s , v t ∈ R p , E(u t ) = 0, E(v t ) = 0 and α 1 is as before. To pass from the state space representation (4) and (5) to (1) and (2), let V t be the symmetric covariance matrix Every symmetric matrix, M , satisfies the decomposition M = M 1/2 M 1/2 , where M 1/2 is a square nonunique matrix. For example, let O be an orthogonal matrix such that O M O = D, where D is a diagonal matrix. Then, we can take M 1/2 = OD 1/2 , where D 1/2 is the matrix obtained from D by replacing its nonzero elements with their square roots. This choice of M 1/2 has the advantage of being valid and numerically stable even if M is singular. It follows from this that we can take It could be useful to compare the state space model used in SSMMATLAB with the ones used in some other statistical software packages. For example, the state space model considered in MATLAB corresponds to a VARMAX model. It is of the form where {A t } is an innovations sequence (uncorrelated, zero mean and with common covariance matrix) and {u t } is a sequence of exogenous variables. Since Gu t = vec(Gu t ) = (u t ⊗I)vec (G) and Hu t = vec(Hu t ) = (u t ⊗I)vec(H), if we define β h = vec(H), β g = vec(G), β = (β h , β g ) , W t = (0, u t ⊗ I) and V t = (u t ⊗ I, 0), we see that the previous state space model is as (4) and (5) (4) and (5), but there is no W t β term and the errors u t and v t are uncorrelated. The state space models in STAMP and REGCMPNT are also as in the State Space Models (SSM) toolbox for MATLAB. In the state space model considered in Stata the system matrices are time invariant and the errors are uncorrelated.
To the best of this author's knowledge none of the software packages mentioned in Section 1 handles either VARMAX models in echelon form and Kronecker indices or the use of high pass and band pass filters in a model-based approach.

Putting some common models into state space form
Given that the state space model considered by SSMMATLAB is very general, it is advisable to have some functions that allow to put some of the most commonly used models in practice into state space form. In this section, we will document some functions that can be used for this purpose. The user can of course modify these functions or write his own functions in order to suit his needs, but in many cases these functions will be sufficient.

Theoretical introduction
Suppose a vector ARMA (VARMA) model given by that can be written more compactly as The model (6) is stationary if the roots of det[Φ(z)] are all outside the unit circle and the model is invertible when the roots of det[Θ(z)] are all outside the unit circle.
To obtain initial conditions for the Kalman filter, the mean and the covariance matrix of the initial state vector are needed. If the series is stationary, the mean is obviously zero. As for the covariance matrix, letting VAR(α 1 ) = V , the matrix V satisfies the Lyapunov equation where T and H are given by (7). In SSMMATLAB, this equation is solved in a numerically stable manner.
The VARMA models considered in SSMMATLAB can be multiplicative, i.e., they can be of the form (8) where s is the number of observations per year.
There also exists the possibility to incorporate regression variables into the model. More specifically, models of the form where U t follows a VARMA model (8) and β is a vector of regression coefficients, can be handled in SSMMATLAB.

SSMMATLAB implementation
In SSMMATLAB, the matrix polynomials in (8)  Model estimation Once the model has been defined and, therefore, the structure str exists, it can be estimated. Before estimation, the user has to decide whether there are fixed parameters in the model or not. How to fix some parameters has been explained in the previous section. The parameters to estimate are in the array str.xv, and the fixed parameters are in str.xf. It is assumed that the values entered by the user for the parameters to be estimated are reasonable initial values. In any case, the estimation function checks at the beginning whether the model is stationary and invertible and issues a warning message if the model is nonstationary or noninvertible.
One method that usually provides good initial estimates for VARMA models is the Hannan and Rissanen (1982) method for univariate series or its generalization to multivariable series Kavalieris 1984, 1986). This method has been recently implemented in SSMMATLAB and will be described later.
It should be emphasized that in SSMMATLAB, the (1, 1) parameter in the covariance matrix of the innovations is always concentrated out of the likelihood.
During the estimation process, each time the log-likelihood is evaluated SSMMATLAB checks whether the model is stationary and invertible. In case any of these conditions is not satisfied, the variable in the corresponding matrix polynomial is multiplied by a small number so that all its roots are outside the unit circle. This guarantees that the solution will always be stationary and invertible.
The following function can be used for parameter estimation.
function result = varmapqPQestim(y, str, Y) After model estimation, the function pr2varmapqPQ can be used to set up the estimated model in VARMA form. For example, the following commands achieve this.
Note that n ip specifies the number of free coefficients in the polynomial φ ip (z) for i = p. The numbers {n i : i = 1, . . . , s} are called Kronecker indices and l = max{n i : i = 1, . . . , s}.
The state space echelon form corresponding to the previous VARMA echelon form is It is described in more detail in the SSMMATLAB manual. Note that the A t in the state space form (13) and (14) are the model innovations.
In a similar way, VARMAX models in echelon form can be handled in SSMMATLAB, as described in the manual.

SSMMATLAB implementation
As in the case of VARMA models, in SSMMATLAB the matrix polynomials of a VARMA or VARMAX model in echelon form are given as three dimensional arrays in MATLAB.
Once the Kronecker indices for model (9) or for a VARMAX model have been specified, we can use the following function to put this model into state space form, using NaN to represent the parameters that have to be estimated.

function [str, ferror] = matechelon(kro, s, m)
Fixing of parameters to zero The user can fix some parameters to zero in a VARMAX model after the structure str has been created using function matechelon or any other method. To this end, he can set the corresponding parameters of the appropriate matrix polynomials to zero and subtract the number of fixed parameters from str.nparm. For example, in the following lines of MATLAB code first some parameters are fixed to zero after the model has been estimated using the Hannan-Rissanen method. Then, in the next step, the model is re-estimated.

Theoretical introduction
Cointegrated VARMA models can be handled in SSMMATLAB. The VARMA models can be ordinary, multiplicative, or in echelon form. The following discussion is valid for all these types of models. Let the k-dimensional VARMA model be given by where B is the backshift operator, BY t = Y t−1 , Φ(z) = Φ 0 + Φ 1 z + · · · + Φ l z l , Θ(z) = Θ 0 + Θ 1 z + · · · + Θ l z l , Θ 0 = Φ 0 , Φ 0 is lower triangular with ones in the main diagonal, and det[Φ(z)] = 0 implies |z| > 1 or z = 1. We assume that the matrix Π, defined by has rank r such that 0 < r < k and that there are exactly k − r roots in the model equal to one. Then, Π can be expressed (non uniquely) as and define the matrix P as Then, it is not difficult to verify that Q = P −1 is given by and that if we further define U 1 = P 1 Q 1 and U 2 = P 2 Q 2 , the following relations hold Thus, we can write The error correction form corresponding to model (15) is where It follows from (18) that β Y t−1 is stationary because all the terms in this equation different from ΠY t−1 = αβ Y t−1 are stationary. Therefore, there are r cointegrated relations in the model given by β Y t .
Considering (16) and (17), the following relation between the autoregressive polynomials in (15) and (18) holds and the model (15) as Since both U 1 = β ⊥ (β ⊥ β ⊥ ) −1 β ⊥ and U 2 = β(β β) −1 β are idempotent and symmetric matrices of rank k − r and r, respectively, the eigenvalues of these two matrices are all equal to one or zero. In particular, (19) is a "differencing" matrix polynomial because it contains all the unit roots in the model. This implies in turn that the matrix polynomial Φ * (z) in (19) has all its roots outside the unit circle and the series D(B)Y t in (20) is stationary. Thus, the matrix polynomial Φ * (z) can be inverted so that the following relation holds Premultiplying the previous expression by β ⊥ , we can see that there are k − r linear combinations of Y t that are I(1) given by β ⊥ Y t . In a similar way, premultiplying by β , it follows as before that there are r linear combinations of Y t that are I(0) given by β Y t .
The series D(B)Y t can be considered as the "differenced series", and the notable feature of (20) is that the model followed by D(B)Y t is stationary. Therefore, we can specify and estimate a stationary VARMA model if we know the series D(B)Y t .
It is shown in the SSMMATLAB manual how the matrix U 1 can be parameterized.

SSMMATLAB implementation
In SSMMATLAB there are two ways to handle cointegrated VARMA models. The first one parameterizes model (18) in terms of the matrix polynomials Γ(z) and Θ(z) and the matrices α and β ⊥ . The second one parameterizes model (20) in terms of the matrix polynomials Φ * (z) and Θ(z) and the matrix β ⊥ . The advantage of the latter parametrization is that we can specify a stationary VARMA model in echelon form for the "differenced" series by directly specifying Φ * (z) and Θ(z). There is no need for a reverse echelon form considered by some authors (Lütkepohl 2007).
Once the cointegration rank has been identified, the matrix β ⊥ and the differencing matrix polynomial D(z) can be estimated in SSMMATLAB using the following function, that also gives the "differenced" series.
function [D, DA, yd, ferror] = mdfestim1r(y, x, prt, nr) If a model such as (20) has been identified for the "differenced" series, D(B)Y t , the following function can be used in SSMMATLAB to obtain the matrix polynomial Γ(z) and the matrices α and β ⊥ corresponding to the error correction model If a model in error correction form (18) has been identified, the following function can be used in SSMMATLAB to obtain the matrix polynomial Φ * (z), the matrix β ⊥ and the differencing matrix polynomial, D(z), corresponding to the model for the "differenced" series.
function [phi, D, DA, ferror] = mecf2mid(Lambda, alpha, betap) Estimating the number of unit roots in the model The number of unit roots in the model can be obtained in SSMMATLAB using a generalization to multivariate series of the criterion based on different rates of convergence proposed by Gómez (2013) for univariate series. The following function can be used in SSMMATLAB for that purpose.
function [D, nr, yd, DA, ferror] = mcrcregr(y, x) Model estimation Before estimating a model parameterized in terms of the matrix polynomials Γ(z) and Θ(z) and the matrices α and β ⊥ corresponding to the error correction model (18) If the model is parameterized in terms of the matrix polynomials Φ * (z) and Θ(z) and the matrix β ⊥ corresponding to the model for the "differenced" series (20), we can put the model into state space form in SSMMATLAB using first the functions suvarmapqPQ, described earlier, and then the following function. In the SSMMATLAB manual, it is described how to set up the estimated model and how to forecast with it if desired.

Theoretical introduction
Univariate structural models are models in which the observed univariate process, {y t }, is assumed to be the sum of several unobserved components. In its general form, the model is where p t is the trend, s t is the seasonal, u t is the cyclical, v t is the autoregressive, and e t is the irregular component. Each of these components follows an ARIMA model. All the models described later in this section can be handled in SSMMATLAB.
The trend component is usually specified as where {c t } and {d t } are two mutually and serially uncorrelated sequences of random variables with zero mean and variances σ 2 c and σ 2 d . The idea behind the previous model is to make the slope and the intercept stochastic in a linear equation, p t = p + b(t − 1), and to let them vary according to a random walk. In fact, if σ 2 c = 0 and σ 2 d = 0 we get the deterministic linear There are basically two specifications for the seasonal component. The first one is called "stochastic dummy seasonality" and, according to it, s t follows the model where S(B) = 1 + B + · · · + B f −1 , B is the backshift operator, By t = y t−1 , f is the number of observations per year and {r t } is an uncorrelated sequence of random variables with zero mean and variance σ 2 r . The idea behind this model is that the seasonal component is periodic and its sum should be approximately zero in one year. The other representation is called "trigonometric seasonality" and in this case s t follows the model where [x] denotes the greatest integer less than or equal x, f is, as before, the number of observations per year, s it follows the model ω i = 2πi/f , and {j i,t } and {j * i,t } are two mutually and serially uncorrelated sequences of random variables with zero mean and common variance σ 2 i . If f is even, ω f /2 = 2π[f /2]/f = π and the model followed by the component s f /2,t in (21), corresponding to the frequency ω f /2 , collapses to s f /2,t+1 = −s f /2,t + j f /2,t . In SSMMATLAB, it is assumed that all seasonal components have a common variance, This representation of the seasonal component has its origin in the observation that, from the theory of difference equations, we know that the solution of the equation S(B)s t = 0 is the sum of [f /2] deterministic harmonics, each one corresponding to a seasonal frequency ω i .
If the cyclical component, u t , is present, it can be modeled in two different ways. The first one corresponds to that proposed by Harvey (1993), namely where 0 < ρ < 1, θ ∈ [0, π] is the cyclical frequency, and {k t } and {k * t } are two mutually and serially uncorrelated sequences of random variables with zero mean and common variance σ 2 k . The ρ factor ensures that the cycle is stationary.
It can be shown that the initial conditions for the cycle (22) satisfy where I 2 is the unit matrix of order two and (1 − ρ 2 )σ 2 u = σ 2 k . Here, the notation ∼ refers to the first two moments of the distribution of u 1 and u * 1 , meaning that these two variables have zero mean, are uncorrelated and have common variance σ 2 u . The second way to model the cycle has its origin in the model-based interpretation of a bandpass filter derived from a Butterworth filter based on the sine function. See Gómez (2001) for details. The model for the cycle is in this case where ρ and θ are as described earlier and {k t } is an uncorrelated sequence of random variables with zero mean and variance σ 2 k . The previous model can be put into state space form as To obtain the initial conditions in this case, we can consider that, clearly, u 1 and u 1|0 have zero mean and their covariance matrix, V , satisfies the Lyapunov equation The matrix V is obtained in SSMMATLAB by solving the previous Lyapunov equation in a numerically safe manner.
It is to be noticed that the variables s * i,t , u * i,t and u t|t−1 in (21), (22) and (23) are auxiliary variables used to define the state space forms for the components of interest. In addition, it can be shown that the cycle specified as in (23) can be obtained from (22) if we let {k * t } be deterministic and equal to zero while maintaining {k t } stochastic and without change.
The autoregressive component, v t , is assumed to follow an autoregressive model, i.e., where the polynomial φ(B) = 1 + φ 1 B + · · · + φ p B p has all its roots outside the unit circle and {w t } is an uncorrelated sequence of random variables with zero mean and variance σ 2 w . In SSMMATLAB only one cycle at a time can be specified in the structural model. The reason for this is that cycles are usually difficult to specify and to estimate. Thus, if one believes that there are several cycles in the model, one can specify one cycle and let the autoregressive component take account of the other cycles by specifying a sufficiently high autoregressive order.
There exists the possibility to incorporate regression variables into structural models. More specifically, models of the form where w t follows a structural model and β is a vector of regression coefficients, can be handled in SSMMATLAB. It is also possible to incorporate interventions that affect some component. For example, an impulse to accommodate a sudden change in the slope of the series that takes place at one observation only. This type of intervention can be modeled by defining a proper W t matrix in Equation 1. This procedure will be illustrated in Case Study 3.

SSMMATLAB implementation
The following function can be used in SSMMATLAB to put a univariate structural model into state space form. Model estimation Once the model has been defined and, therefore, the structure str exists, it can be estimated. Before estimation, the user has to decide whether to fix some parameters in the model or not. How to fix some parameters has been explained in the previous paragraph. The parameters to be estimated are in the array str.xv, and the fixed parameters are in str.xf. It is assumed that the values entered by the user for the parameters to be estimated are reasonable initial values.
Initial values that can be used for the standard deviations are all equal to .1, except the slope standard deviation that is usually smaller and can be set to 0.005. Initial values that can be used for the autoregressive parameters are all equal to .1. Initial values for the cycle ρ and frequency parameters can be .9 and a frequency that can be considered reasonable by the user.
If the user has not selected a variance to be concentrated out using the field comp.conout, the program will select the biggest variance to that effect. After calling function suusm, the index for the parameter to be concentrated out is in the field str.conc.
The following function can be used for parameter estimation.

function [result, str] = usmestim(y, str)
It is to be noticed that, even if the user has selected a variance to be concentrated out using the field comp.conout, the program will always check whether the selected variance is the biggest one. To this end, a preliminary estimation is performed in usmestim. After it, if the biggest estimated variance does not correspond to the initially selected parameter to be concentrated out, the program will change this parameter and will make the necessary adjustments in structure str. Therefore, structure str can change after calling usmestim. The actual estimation is performed after the previous check.
After model estimation, function pr2usm can be used to set up the estimated structural model. For example, the following commands achieve this. xvf = result.xvf; xf = result.xf; [X, Z, G, W, T, H, ins, ii, ferror] = pr2usm(xvf, xf, str); Recursive residuals As explained in Section 5.2, recursive residuals can be computed using functions scakff or scakfff. For example, the following command can be used to compute recursive residuals after estimation of a structural model, assuming that X, Z, G, W, T and H are the estimated matrices and ins and ii contain the initial conditions.
[Xt, Pt, g, M, initf, recrs] = scakff(y, X, Z, G, W, T, H, ins, ii); Forecasting As described in Section 5.4, forecasts can be obtained using function ssmpred. For example, the following commands can be used to obtain ten forecasts after estimating a structural model, assuming that X, Z, G, W, T, H, ins and ii are as in the previous paragraph. Smoothing As described in Section 5.5, smoothing can be performed using function scakfs.
For example, assuming that ten forecasts have been previously obtained after estimating a structural model and that the series forecasts are in array pry and the state vector forecasts are in array alpr, the following commands can be used to estimate the trend using smoothing, extend it with the forecasts, and display both the extended original and trend series. It is further assumed that X, Z, G, W, T and H are the estimated matrices and ins and ii contain the initial conditions. The following function can be used after smoothing to select a desired smoothed component. This function works with structural as well as with AMB unobserved components models. These last models will be introduced in Section 3.5.
function Cc = dispcomp (KKP, str, comp, varargin) For example, the following lines can be used to select the smoothed cycle and to plot it after the smoothed components have been obtained using function scakfs.

Theoretical introduction
The ARIMA model-based (AMB) method to decompose a given time series that follows an ARIMA model into several unobserved components that also follow ARIMA models is described in, for example, Gómez and Maravall (2001b, Chapter 8).
This approach was originally proposed by Hillmer and Tiao (1982). The idea is based on a partial fraction expansion of the pseudospectrum of an ARIMA model specified for the series at hand, {y t }. According to this decomposition, terms with denominators originating peaks at the low frequencies should be assigned to the trend component, terms with denominators originating peaks at the seasonal frequencies should be assigned to the seasonal component, and the other terms should be grouped into a so-called "stationary component". This last component can in turn be decomposed into an irregular (white noise) plus some other, usually moving average, component. For example, consider the model where ∇ = 1 − B, B is the backshift operator, By t = y t−1 , and {a t } is a white noise sequence with zero mean and VAR(a t ) = σ 2 . Given that (1 − z)(1 − z 4 ) = (1 − z) 2 (1 + z + z 2 + z 3 ), the pseudoespectrum is where A(x) and B(x) are polynomial functions of cos(x) to be determined. To see this, consider that, setting y = e −ix + e ix as the new variable, any pseudospectrum can be written as a quotient of polynomials in y = 2 cos(x).
In the previous decomposition of f (x), the first term on the right hand side becomes infinite at the zero frequency and should be assigned to the trend, whereas the second term becomes infinite at the seasonal frequencies, π and π/2, and should, therefore, be assigned to the seasonal component. However, both the seasonal and the trend components are not identified because it is possible that one may subtract some positive quantity from each of the terms on the right hand side and at the same time add it as a new term in the decomposition of f (x), so that we would obtain where A(x) and B(x) are new polynomial functions in cos(x) and k is a positive constant. This positive constant gives rise to a new white noise component.
To identify the components, the so-called canonical decomposition is performed. According to this decomposition, a positive constant, as big as possible, is subtracted from each term on the right hand side. In this way, the components are made as smooth as possible and become identified. The resulting components are called canonical components. The canonical decomposition does not always exist and this constitutes a flaw in the procedure. However, there are simple solutions to this problem.
It can be shown that the trend and seasonal components, p t and s t , corresponding to the previous example are of the form

and
( where {b t } and {c t } are two uncorrelated white noises and the polynomial 1+β 1 z +β 2 z 2 +β 3 z 3 has at least one root in the unit circle. In addition, the equality y t = p t + c t + i t holds, where {i t } is white noise. If logs of the series, y t , are taken, then the procedure is applied to the transformed series. Thus, in order to obtain the multiplicative components one has to exponentiate the components obtained from the decomposition of log(y t ). This may cause problems with the estimated trend because usually the annual trend sums are lower than the annual sums of the original series, a phenomenon due to geometric means being smaller than arithmetic means. For this reason, some kind of "bias" correction is usually applied to the estimated trend. This problem is also present in structural models.

SSMMATLAB implementation
Before performing the canonical decomposition, it is necessary to select the roots in the autoregressive polynomial that should be assigned to the trend and the seasonal components.
The following function can be used in SSMMATLAB for that purpose.
function [phir, phis, thr, ths, phirst] = arima2rspol(phi, Phi, th, Th, ... freq, dr, ds) Once the model has been decomposed into its canonical components, one can put the unobserved components model into state space form and perform forecasting and smoothing in the same way as that previously described for structural models.
To put the model into state space form, the following SSMMATLAB function can be used.
function [X, Z, G, W, T, H, ins, ii, strc, ferror] = sucdm(comp, y, Y, ... stra, npr) The following lines of MATLAB code illustrate how to first decompose an ARIMA model into its canonical components and then how to smooth these components. Finally, the original series as well as the trend-cycle are displayed. The model is given by the polynomials phi, Phi, th and Th. The regular and seasonal differences are one and one, respectively. The standard deviation of the residuals is sconp.

Estimation of smooth trends and cycles
In the AMB approach, it is not usually possible to directly estimate cycles. This is due to the fact that the majority of the ARIMA models fitted in practice do not have autoregressive components with complex roots that may give rise to cyclical components. Trend components given by the AMB approach are for this reason also called "trend-cycle" components. For similar reasons, it is also usually not possible to estimate smooth trends using only the unobserved components given by the canonical decomposition.
To estimate smooth trends and cycles within the AMB approach, one possibility is to incorporate fixed filters into the approach in the manner proposed by Gómez (2001).
The filters considered in SSMMATLAB for smoothing trends are two-sided versions of Butterworth filters. Butterworth filters are low-pass filters and they are of two types. The first one is based on the sine function (BFS), whereas the second one is based on the tangent function (BFT). See, for example, Otnes and Enochson (1978).
The squared gain of a BFS is given by where x denotes angular frequency and x c is such that |G(x c )| 2 = 1/2. These filters depend on two parameters, d and x c . If x c is fixed, the effect of increasing d is to make the fall of the squared gain sharper. BFSs are autoregressive filters of the form H(B) = 1/θ(B), where B is the backshift operator, By t = y t−1 , θ(B) = θ 0 + θ 1 B + · · · + θ d B d and |G(x)| 2 = H(e −ix )H(e ix ). Thus, if {y t } is the input series, the output series, {z t }, is given by the recursion To start the recursion at t = 1 say, some initial values, z 1−d , . . . , z 0 , are needed. It can be shown that H s (B, F ) can be given a model-based interpretation. It is the Wiener-Kolmogorov filter to estimate the signal in the signal plus noise model

The BFSs used in
under the assumption that the signal s t follows the model ∇ d s t = b t , where {b t } is a white noise sequence with zero mean and unit variance and {b t } is independent of the white noise sequence {n t }. The estimator of s t is given bŷ The weights ν k in (27) can be obtained from the signal extraction formula where λ = VAR(n t ). The frequency response function,Ĥ s (x), of the filter H s (B, F ) is obtained from (28) by replacing B and F with e −ix and e ix , respectively. After some manipulation, it is obtained thatĤ where λ = [2 sin(x c /2)] −2d . Thus, the gain, |Ĥ s (x)|, of H s (B, F ) coincides with the squared gain of a BFS. See Gómez (2001) for details.
For BFT, the squared gain function is given by (25), but replacing the sine function by the tangent function. The filter is of the form H( To design a BFS in SSMMATLAB, the following function can be used. The following function can be applied in SSMMATLAB to design a BFT.
function ggsintanbut (D, Thetap, Thetas, d, thc) The following lines of MATLAB code can be used to first specify a BFS that coincides with the Hodrick-Prescott filter and then to plot the gain function of the filter. The filter is specified giving the parameters Lambda and Di. To estimate cycles in SSMMATLAB, one can use band-pass filters derived from BFT. These are two-sided filters that can be obtained by estimating signals which follow the model (1 − 2 cos xB + B 2 ) d s t = (1 − B 2 ) d b t in the signal plus noise model (26). Details regarding the design of these band-pass filters and their model-based interpretation can be found in Gómez (2001).
To design a band-pass filter in SSMMATLAB, the following function can be used. To plot the gain function of a band-pass filter in SSMMATLAB, the following function can be used.
function ggbptanbut (D, omp1, omp2, oms2, d, alph, lambda) The following MATLAB code lines illustrate how to first design a band-pass filter to be applied to quarterly data to estimate a cycle with frequencies in the business cycle frequency band (periods between a year and a half and eight years). Then, it is shown how the gain of the designed filter can be plotted. All of the previously described filters are fixed filters. However, they can be incorporated into the AMB approach as described in Gómez (2001). See the SSMMATLAB manual in Gómez (2014) for more details.
The following function can be used in SSMMATLAB to set up a state space model for an unobserved components model, where the components are obtained in the manner previously described given information from both the canonical decomposition of an ARIMA model and a designed BFS or BFT.
function [X, Z, G, W, T, H, ins, ii, strc, ferror] = sucdmpbst (comp, ... compf, y, Y, stra, npr) If instead of a low-pass filter (BFS or BFT), as in the previous function, a band-pass filter based on BFT is applied, the following function can be used to set up the appropriate state space model. Once we have a state space model in which the trend-cycle given by the AMB approach has been further decomposed into a smooth trend and a cycle by means of a fixed filter of the type BFS, BFT or band-pass filter based on BFT, we can use the Kalman filter to smooth the components. This can be done in SSMMATLAB by using the function scakfs. For example, the following lines of MATLAB code can be used to first design a low-pass filter of the BFS type, the Hodrick-Prescott filter, and then to estimate the unobserved components. Finally, the smooth trend and the cycle, estimated as the difference between the trend-cycle and the smooth trend, are plotted.

Identification and estimation of VAR(MA)X models
In this section, we will describe a number of tools available in SSMMATLAB for the identification and estimation of VARX and VARMAX models.

VARX identification and estimation
The VARX models considered in SSMMATLAB are of the form These models are important because every VARMAX model can be approximated to any degree of accuracy by a VARX model with a sufficiently high order.
Although a VARX model can be put into state space form, VARX models are estimated in SSMMATLAB using OLS. The reason why these models are included in SSMMATLAB is that they usually constitute a good starting point when analyzing multivariate models, like VARMAX models. To estimate a VARX model in SSMMATLAB, the following function can be used.
function res = varx_est(y, nlag, x, test, xx) When estimating a VARX model, sometimes only the residuals are desired. In this case, the following function can be used in SSMMATLAB.

function resid = varx_res(y, nlag, x)
To identify the order of a possibly nonstationary VARX model, the likelihood ratio criterion can be used. The following function can be applied in SSMMATLAB for this purpose.

function [lagsopt, initres] = lratiocrx(y, maxlag, minlag, prt, x)
When there are no exogenous variables, that is, when the model is a VAR model, the following function can be used in SSMMATLAB for model estimation.
function res = var_est(y, nlag, test, x) If only the residuals are desired when estimating a VAR model, the following function can be used in SSMMATLAB.

function resid = var_res(y, nlag, x)
To determine the optimal lag length of a possibly nonstationary VAR model using the likelihood ratio criterion, the following function can be used in SSMMATLAB.

function [lagsopt, initres] = lratiocr(y, maxlag, minlag, prt, x)
In the case of a possibly nonstationary VAR model, the following function can be used in SSMMATLAB to determine the optimal lag length using the AIC or BIC criterion.

Multivariate residual diagnostics
To estimate the covariances and the autocorrelations, as well as the portmanteau statistics of a multivariate time series, the following function can be used in SSMMATLAB.

Identification of VARMAX(p, q, r) models
The following function can be used in SSMMATLAB to identify VARMAX(p, q, r) models. It applies a sequence of likelihood ratio tests to obtain the orders.

Identification of VARMAX models in echelon form
The following two functions can be used in SSMMATLAB to identify the Kronecker indices for VARMAX models in echelon form. The first one identifies and estimates in a preliminary step a VARMAX(p, q, r) model, whereas the second one starts by identifying and estimating a VARMAX(p, p, p) model.

The Hannan-Rissanen method to estimate VARMAX models
Although state space models can be directly estimated using regression techniques, like subspace methods, these methods involve the estimation of a large number of parameters as soon as the dimension of the state vector increases. For this reason, the approach adopted in SSMMATLAB is to use the Hannan-Rissanen method, that applies regression techniques only and is based on the VARMAX specification of the model. Even though it does not use state space models, it usually gives very good starting values when estimating a VARMAX model in state space echelon form by maximum likelihood.

Theoretical introduction
Suppose that the process {Y t } follows the VARMAX model in echelon form where Φ 0 = Θ 0 is a lower triangular matrix with ones in the main diagonal. Equation 31 can be rewritten as

SSMMATLAB: State Space Models in MATLAB
Applying the vec operator to (32), it is obtained that where The parameter restrictions given by the echelon form (31) can be incorporated into Equation 33 by defining a selection matrix, R, containing zeros and ones such that where β is the vector of parameters that are not restricted in the matrices Φ i , Ω i or Θ i , i = 0, 1, . . . , r. Using (34), Equation 33 can be rewritten as where X t = W t R. Notice that, as mentioned earlier, X t is uncorrelated with A t in (35) and that if we knew X t , we could estimate β by OLS. The idea behind the Hannan-Rissanen method is to estimate β in (35) after we have replaced the unknown innovations in X t with those estimated using a VARX model.
The Hannan-Rissanen method is described in more detail in the SSMMATLAB manual.

SSMMATLAB implementation
As mentioned earlier, the Hannan-Rissanen method usually gives very good starting values when estimating a VARMAX model in state space echelon form by maximum likelihood. The following function can be used in SSMMATLAB to estimate a VARMAX model using this method.
function [str, ferror] = estvarmaxkro (y, x, seas, kro, hr3, finv2, ... mstainv, nsig, tsig) For example, the following lines of MATLAB code can be used to estimate a transfer function model with one input series. Both, the input and the output series, are differenced prior to estimation. All polynomials have degree two. yd = diferm(y, 1); xd = diferm(x, 1); kro = 2; hr3 = 0; finv2 = 1; strv = estvarmaxkro(yd, xd, seas, kro, hr3, finv2); As described earlier, if there are parameters that are not significant after estimation, it is possible to fix them to zero and estimate the model again. The following lines can be used in the previous example to fix some parameters to zero and re-estimate the model. The estimation is performed using function mhanris, that will be described later in this section.
To re-estimate a VARMAX model in SSMMATLAB after having fixed some parameters to zero, the following function can be used.
function str = mhanris (y, x, seas, str, hr3, finv2, mstainv, nsig, tsig) Sometimes a VARMAX model is given as a multiplicative VARMA model with exogenous inputs. In this case, the following function can be used in SSMMATLAB to estimate the models in this form.
The following function can be used in SSMMATLAB to estimate a VARMAX model using the conditional method.
function [xvf, str, ferror] = mconestim(y, x, str) Conditional residuals After estimating a VARMAX model using the conditional method, the conditional residuals are in the field residcon of the structure str given as output by function mconestim. For example, the following lines of MATLAB code can be used to plot the conditional residuals of a bivariate model and their simple and partial correlograms after estimation.

The exact ML method to estimate VARMAX models
After a VARMAX model has been estimated using the Hannan-Rissanen or the conditional method, the user may be interested in estimating the model using the exact maximum likelihood (ML) method.
The following function can be used in SSMMATLAB to estimate a VARMAX model using the exact ML method.
function [xvf, str, ferror] = mexactestimc(y, x, str, Y) Recursive residuals After estimating a VARMAX model using the exact ML method, the following function can be used to obtain the recursive residuals.
function [ff, beta, e, f, str, stx, recrs] = exactmedfvc(beta, y, x, ... str, Y, chb) For example, the following lines of MATLAB code illustrate the estimation of a bivariate VARMAX model using the exact ML method. Then, some diagnostic statistics based on the recursive residuals are computed. Finally, the recursive residuals and their simple and partial correlograms are plotted. Forecasting To obtain some forecasts after estimating a VARMAX model using the conditional or the exact method, the observed series, y t , that is assumed to follow the state space model in echelon form where V t is the exogenous part, that depends on the inputs x t and their initial condition only, and U t is the endogenous part, that depends on the innovations a t and their initial condition only. Then, the forecasts can be obtained separately by forecasting V t and U t , that are uncorrelated. To this end, one can use functions ssmpredexg and ssmpred, respectively. The latter is described in Section 5.4.2, whereas the former is as follows.
function [ypr, mypr, alpr, malpr] = ssmpredexg(n, x, stx, sts) It is to be noted that if the inputs are stochastic, a model for them must be provided by the user. This model will be used in function ssmpredexg to obtain the input forecasts. If the inputs are not stochastic, the user must provide the forecasts.

Theoretical introduction
As described in the SSMMATLAB manual, the Kalman filter can be used for likelihood evaluation. Assuming β = 0 in (1) and (2) and δ = 0, a = 0 in (3), the Kalman filter is given for t = 1, . . . , n by the recursions initialized withα 1|0 = a and P 1 = Ω. More details are given in the SSMMATLAB manual.

SSMMATLAB implementation
The MATLAB function used in SSMMATLAB for likelihood evaluation is function [e, f, hb, Mb, A, P, qyy, R] = scakfle2(y, X, Z, G, W, T, H, ... ins, i, chb) For parameter estimation, we first use function scakfle2 to obtain the residual vector e and the constant f . Then, we multiply e by f to get F = ef , the vector of nonlinear functions that has to be minimized. Using the notation of the Kalman filter equations (36), e t = Σ −1/2 t E t , e = (e 1 , . . . , e n ) and f = n t=1 |Σ t | 1/(2np) . More details are given in the SSMMATLAB manual.

Missing values
If the series has missing values, these should be replaced with the symbol NaN in MATLAB. The algorithms in SSMMATLAB are designed to take account of the missing values. For example, for a univariate series that follows an ARIMA model, each time the Kalman filter encounters a missing value, it skips this observation, sets K t = 0 and continues filtering.

Theoretical introduction
Recursive residuals can be of two types, depending on whether one considers the estimated regression parameters fixed, together with the other parameters of the model, or not. More details are given in the SSMMATLAB manual.

SSMMATLAB implementation
When the estimated regression parameters are not considered fixed, in SSMMATLAB a square root information filter is applied to obtain the recursive residuals. The following function can be used for that purpose. It is to be noted that this function also provides the filtered state estimates, that is, the estimate of the state α t based on the observations Y 1 , . . . , Y t , as well as their MSE. When the estimated regression parameters are considered fixed, together with the other parameters of the model, the following function can be used in SSMMATLAB to obtain the recursive residuals.

Theoretical introduction
Once the state space model has been defined and assuming that reasonable initial parameter values are available, the model can be estimated. It is to be emphasized that in SSMMATLAB we always concentrate one parameter out of the likelihood in the covariance matrix of the errors of the state space model. As shown in Section 5.1, this allows for the transformation of the log-likelihood maximization problem into a minimization of a nonlinear sum of squares function. In SSMMATLAB, the optimization method used is that of Levenberg-Marquardt (Levenberg 1944;Marquardt 1963). This method has been proved in practice to be a reliable method for minimizing a nonlinear sum of squares function.

SSMMATLAB implementation
The following function can be used in SSMMATLAB for parameter estimation.

Theoretical introduction
Let the forecast or, equivalently, the orthogonal projection of α n+h onto the sample Y = (Y 1 , . . . , Y n ) beα n+h|n , where h ≥ 1. Then, the h-period-ahead forecasts and their mean squared error, P n+h , can be recursively obtained bŷ whereγ n+1 and Π n+1 are the GLS estimator of γ based on Y and its MSE and for h > 1

SSMMATLAB implementation
In SSMMATLAB, the following function can be used for forecasting.

Theoretical introduction
For smoothing, the Bryson-Frazier recursions are used, as described in the SSMMATLAB manual.

SSMMATLAB implementation
The following function can be used in SSMMATLAB for smoothing of the state vector.
function [KKP, PT, hd, Md] = scakfs(y, X, Z, G, W, T, H, ins, i) If it is of interest to smooth a general vector of the form Y t = U t β + C t α t +D t t , the following function can be used.

Examples and case studies
In the SSMMATLAB manual in Gómez (2014) some examples are given on how to use SSM-MATLAB in practice. These include estimation, computation of recursive residuals, forecasting and smoothing using univariate ARMA and ARMAX models, VARMA and VARMAX models, AMB unobserved components models and univariate structural models. In addition, five case studies illustrate the use of SSMMATLAB to analyze sophisticated state space models that cannot be dealt with standard commercial packages.
As an illustration, we will present in this section one example and one case study. The example is Example 5 in the SSMMATLAB manual. In this example, the series of ozone levels used by Box and Tiao (1975) to introduce Intervention Analysis is considered. Unlike these authors, a structural model instead of an ARIMA model is specified. The model considered has a deterministic level, a seasonal component that is modeled as trigonometric seasonality and an autoregressive component of order one. In addition, there are three intervention variables corresponding to the interventions in Box and Tiao (1975). The MATLAB script file usm2_d.m contains the instructions for putting the model into state space form, and for model estimation, computation of recursive residuals, forecasting and smoothing of the trend. In Figure 1, one can see the original series together with the seasonally differenced series and its sample autocorrelations and partial autocorrelations. In Figure 2, the original and the trend series are displayed together with twelve forecasts of both series.
The case study is Case Study 4 in the SSMMATLAB manual. The purpose of the analysis is to first estimate the business cycle of the US Industrial Production Index for the periods 1946.Q1 through 2011.Q3. Then, to obtain bootstrap samples of the estimated cycle. The estimated cycle can be used as business cycle indicator while studying the cyclical comovements of different series. The cycle is estimated using two different methods. The first one consists of fitting a structural model that includes a cycle. The second one applies the AMB methodology described in Section 6.4 of the SSMMATLAB manual.
Two script files are used. In the first one, USIPIstscl_d.m, a structural model that includes a cycle is first fitted to the data and the cycle is estimated. Then, bootstrap samples of this estimated cycle are obtained using the algorithm proposed by Stoffer and Wall (2004).
In the second file, USIPIcdstcl_d.m, the procedure proposed by Gómez (2001) is applied. Prior to the analysis, an ARIMA model was identified using the program TRAMO of Gómez and Maravall (2001a). In the script file, first the identified model is estimated by exact maximum likelihood and the models for the unobserved components, trend-cycle, seasonal and irregular, are obtained by means of the canonical decomposition. Then, based on the trend-cycle model and the model corresponding to the band-pass filter, models for the two subcomponents of the trend-cycle, the cycle and the smooth trend, are obtained as explained in Gómez (2001). After putting the new model into state space form, the Kalman filter and smoother are applied to estimate the smooth trend, the cycle, the seasonal and the irregular. Finally, bootstrap samples of the estimated cycle are obtained using the same method as in the case of the structural model.
The cycle estimated with the structural model is displayed in Figure 3. In Figure 4, the cycle estimated with a band-pass filter within the AMB approach is shown. It is seen that this cycle is smoother than the cycle estimated with the structural model. Finally, in Figure 5, two cycles are displayed. They correspond to two filters applied within the AMB approach. The first filter is the band-pass filter mentioned previously and the second one is a low-pass filter well known to economists, the Hodrick-Prescott filter corresponding to the parameter λ = 1600. Again, the cycle obtained with the band-pass filter is smoother than the one obtained with the low-pass filter. This is due to the fact that the band-pass filter extracts the components corresponding to the business cycle frequencies better than the low-pass filter.