The State Space Models Toolbox for MATLAB

State Space Models ( SSM ) is a MATLAB toolbox for time series analysis by state space methods. The software features fully interactive construction and combination of models, with support for univariate and multivariate models, complex time-varying (dynamic) models, non-Gaussian models, and various standard models such as ARIMA and structural time-series models. The software includes standard functions for Kalman ﬁl-tering and smoothing, simulation smoothing, likelihood evaluation, parameter estimation, signal extraction and forecasting, with incorporation of exact initialization for ﬁlters and smoothers, and support for missing observations and multiple time series input with common analysis structure. The software also includes implementations of TRAMO model selection and Hillmer-Tiao decomposition for ARIMA models. The software will provide a general toolbox for time series analysis on the MATLAB platform, allowing users to take advantage of its readily available graph plotting and general matrix computation capabilities.


Introduction
State Space Models (SSM) is a MATLAB (The MathWorks, Inc. 2007) toolbox for time series analysis using general state space models and the Kalman filter (Durbin and Koopman 2001). The goal of this software package is to provide users with an intuitive, convenient and efficient way to do general time series modeling within the state space framework. Specifically, it seeks to provide users with easy construction and combination of arbitrary models without having to explicitly define every component of the model, and to maximize transparency in their data analysis usage so no special consideration is needed for any individual model. This is achieved through the unification of all state space models and their extension to non-Gaussian and nonlinear special cases (those which can be linearized). The user creation of custom mod-els is also implemented to be as general, flexible and efficient as possible. Thus, there are often multiple ways of defining a single model and choices as to the parametrization versus initialization and to how the model update functions are implemented. Stock model components are also provided to ease user extension to existing predefined models. Functions that implement standard algorithms such as the Kalman filter and state smoother, log likelihood calculation and parameter estimation will work across all models, including any user defined custom models. The state space model manipulation procedures are implemented through object-oriented programming primitives provided by MATLAB and classes in the toolbox are defined to conform to MATLAB conventions whenever possible. The standard Kalman filter-based state space algorithms are implemented in C. Thus the toolbox combines efficient state space algorithms (by allowing to run C functions in MATLAB) with full integration into the MATLAB computing environment. Model building is intuitive and model objects behave like standard MATLAB objects, whereas application of state space algorithms to these objects is as fast as comparable programs in C code. Data exploration and custom model building using SSM is also greatly facilitated by MATLAB's interactive environment.
The result is an integrated toolbox with support for general state space models and standard state space algorithms, complemented by the built-in matrix computation and graphic plotting capabilities of MATLAB.

The state space model equations
This section presents a summary of the basic definition of models supported by SSM. Currently SSM implements the Kalman filter and related algorithms for model and state estimation, hence non-Gaussian or nonlinear models need to be approximated by linear Gaussian models prior to or during estimation. However the approximation is done automatically and seamlessly by the respective routines, even for user-defined non-Gaussian or nonlinear models.
The following notation for various sequences will be used throughout the paper: y t p × 1 Observation sequence ε t p × 1 Observation disturbance (unobserved) α t m × 1 State sequence (unobserved) η t r × 1 State disturbance (unobserved) Linear Gaussian models SSM supports linear Gaussian state space models in the form y t = Z t α t + ε t , ε t ∼ N (0, H t ), α t+1 = c t + T t α t + R t η t , η t ∼ N (0, Q t ), α 1 ∼ N (a 1 , P 1 ), t = 1, . . . , n. (1) Thus the matrices Z t , c t , T t , R t , H t , Q t , a 1 , P 1 are required to define a linear Gaussian state space model. For details of these matrices refer to Commandeur et al. (2011). The matrices and their dimensions are listed here for reference: Non-Gaussian models SSM supports non-Gaussian state space models in the form α 1 ∼ N (a 1 , P 1 ), t = 1, . . . , n. (2) The sequence θ t is the signal and Q t is a non-Gaussian distribution (e.g., heavy-tailed distribution). The non-Gaussian observation disturbance can take two forms: an exponential family distribution or a non-Gaussian additive noise With model combination it is also possible for H t and Q t to be a combination of Gaussian distributions (represented by variance matrices) and various non-Gaussian distributions.
(5) Z t and T t are functions that map m × 1 vectors to p × 1 and m × 1 vectors respectively, and both functions should possess first derivatives. With model combination it is also possible for Z t and T t to be a combination of linear functions (matrices) and nonlinear functions.

Getting started
The easiest and most frequent way to start using SSM is by constructing predefined models, as opposed to creating a model from scratch. This section presents some examples of simple time series analysis using predefined models; the complete list of available predefined models can be found in Appendix A.

Model construction
To create an instance of a predefined model, call the functions ssm_* where the wildcard is the short name for the model, with arguments as necessary. For example: . . .
can be combined by horizontal concatenation More details on the class SSMODEL can be found in Peng and Aston (2007).

Model and state estimation
With the model created, estimation can be performed. SSM expects the data y to be a matrix of dimensions p × n, where n is the data size (or time duration). The model parameters are estimated by maximum likelihood, the SSMODEL class method estimate performs the estimation. For example: model1 = estimate(y, model0) estimates the model and stores the result in model1, where the parameter values of model0 is used as initial value.
After the model is estimated, state estimation can be performed, this is done by the SS-MODEL class method kalman and statesmo, which is the Kalman filter and state smoother respectively.
[a P] = kalman(y, model) applies the Kalman filter on y and returns the one-stepahead state prediction and variance.
[alphahat V] = statesmo(y, model) applies the state smoother on y and returns the expected state mean and variance.
The filtered and smoothed state sequences a and alphahat are m × n+1 and m × n matrices respectively, and the filtered and smoothed state variance sequences P and V are m × m × n+1 and m × m × n matrices respectively, except if m = 1, in which case they are squeezed and transposed. The complete list of data analysis functions can be found in Appendix B.

Toolbox structure and implementation
SSM is composed of two parts, the core C library of state space algorithms (built on top of LAPACK Anderson et al. 1999) and the MATLAB model object library. This section provides a brief overview of the toolbox, for complete details as well as usage particulars refer to the State Space Models manual (Peng and Aston 2007).
C state space algorithms library The core C library has functions for Kalman filtering and smoothing, with support for exact diffuse initialization, missing values and dynamic observation vector length (p t varies through time), and unconditional sampling from state space models. A set of MATLAB C mex functions are built on top of the C library implementing the Kalman filter, state smoother, disturbance smoother, fast state smoother, simulation smoother, loglikelihood and gradient calculation, unconditional sampling, signal extraction and filtering and smoothing weights calculation (support for dynamic p t is suppressed). The C source files are stored in the csrc subdirectory, and can all be automatically compiled by running the script mexall.m.

MATLAB object library
The MATLAB object library consists of five classes SSMAT, SSDIST, SSFUNC, SSPARAM and SSMODEL.
The class SSMAT represents a state space matrix, with elements marked as variable (dependent on model parameters) and/or dynamic (dependent on time). SSMAT objects are designed to mimic MATLAB matrices and, in particular, concatenation and addition are supported.
The classes SSDIST and SSFUNC are derived from SSMAT and represent non-Gaussian distributions and nonlinear functions, respectively. Because SSM currently handles both non-Gaussian and nonlinear model elements by linear approximation of the loglikelihood function, derivation of these classes from state space matrices is natural. User definition of custom SSDIST objects can be a bit involved, as functions that calculate the approximating (heteroscedastic) Gaussian variances from observation data and the smoothed state, and the (non-Gaussian) loglikelihood from estimation innovation, are required. But once these are defined, the model can then be used in all SSM functions without further coding. User definition of SSFUNC requires the nonlinear function itself and its first derivative. The class methods of SSMODEL provides the glue that joins the state space algorithms (ssa) mex functions with the object library. Overhead is minimized by only passing the resolved state space matrices H t , Z t , etc. to the mex functions (in other words, the mex functions have no knowledge of model parameters). Parameter estimation is achieved through a set of update functions that transform parameters to model matrices, which are all stored in the SSMODEL object level for efficiency and flexibility, since often multiple model elements are dependent on parameter values through the same mechanism, storing the update functions in individual SSMATs leads to the need for either messy global variables or duplicated computation. The disadvantage for this solution is increased complexity in defining custom models. Reducing this complexity will be the objective of future work. The parameter estimation class method itself is implemented in MATLAB and makes use of existing optimization routines (e.g., fminsearch), and calls the mex function kalman_int as part of its operation.
A set of functions for the construction of predefined objects are provided (see Appendix A for a list of predefined models). The functions ssm_* returns predefined SSMODEL objects for various models (see the appendix for the complete list of predefined models). Auxiliary functions mat_*, fun_* and dist_* constructs predefined individual MATLAB matrices, update functions and SSDIST objects respectively, and can be variously combined for more convenient and flexible construction of custom models.

Paper overview
The rest of the article proceeds as follows. In Section 2, as in all the papers in this special volume, the local level model will be examined and used to analyze the Nile river data. In the remaining sections, many data analysis examples are conducted, based on the data in the book by Durbin and Koopman (2001), to facilitate easy comparison and understanding of the methods used here. In the interests of brevity, most of the model descriptions are either omitted or only cursorily introduced, as further details can be found either in Durbin and Koopman (2001) or in the introductory paper of this volume (Commandeur et al. 2011).
In addition, all the examples (plus others as well) are included as demos in the software package. Specifically, Section 3 contains a univariate analysis of the road accident data from (Harvey and Durbin 1986), while Section 4 contains a bivariate analysis. Section 5 contains an analysis using ARMA models, while Section 6 contains an analysis using cubic spline methods. Section 7 reviews an application of the software to seasonal adjustment using ARIMA models. Section 8 shows how the software extends to non-Gaussian models.

Local level models -The Nile river data
In this example the Nile data (as seen in the other papers of this volume) will be analyzed using the local level model as described by Harvey (1989). Briefly, in line with the description given in Commandeur et al. (2011), the equations for the local level model are Code for the model construction and estimation for the data y are as follows:

Univariate analysis
In this and the following example, data on road accidents in Great Britain (Harvey and Durbin 1986) is analyzed using structural time series models following Durbin and Koopman  (2001, Section 9.2). The purpose of the analysis is to assess the effect of seat belt laws on road accident casualties, with individual monthly figures for drivers, front seat passengers and rear seat passengers. The monthly price of petrol and average number of kilometers traveled will be used as regression variables. The data is from January 1969 to December 1984.
The drivers series will be analyzed with a univariate structural time series model, which consists of local level component µ t , trigonometric seasonal component γ t , regression component (on price of petrol) and intervention component (introduction of seat belt law) βx t . The model equation is where ε t is the observation noise. The following code example constructs this model: lvl = ssm_llm; seas = ssm_seasonal('trig1', 12); With the model constructed, estimation can proceed with the code: [bstsm logL] = estimate(y, bstsm); [alphahat V] = statesmo(y, bstsm); irr = disturbsmo(y, bstsm); ycom = signal(alphahat, bstsm); ylvl = sum(ycom([1 3 4], :), 1); yseas = ycom(2, :); The estimated model parameters are [0.0037862 0.00026768 1.162e-006], which can be obtained by displaying bstsm.param, the loglikelihood is 175.7790. State and disturbance smoothing is performed with the estimated model, and the smoothed state is transformed into signal components by signal. Because p = 1, the output ycom is M × n, where M is the number of signal components. The level, intervention and regression are summed as the total data level, separating seasonal influence. Using MATLAB graphic functions, the individual signal components and data are shown in Figure 4. The coefficients for intervention and regression are defined as part of the state vector in this model, so they can be obtained from the last two elements of the smoothed state vector (due to the order of component concatenation) at the last time point. The coefficient for the intervention is alphahat(13, n) = -0.23773, and the coefficient for the regression (price of petrol) is alphahat(14, n) = -0.2914. In this way diagnostics of these coefficient estimates can also be obtained by the smoothed state variance V.

Bivariate analysis
The front seat passenger and rear seat passenger series (y2) will be analyzed together using bivariate structural time series model following Durbin and Koopman (2001, Section 9.3), with components as before. Specifically, separate level and seasonal components are defined for both series, but the disturbances are assumed to be correlated. To reduce the number of parameters estimated the seasonal component is assumed to be fixed, so that the total number of parameters is six. We also include regression on the price of petrol (petrol) and kilometers traveled (km), and intervention for only the first series, since the seat belt law only effects the front seat passengers. The following is the model construction code: bilvl = ssm_mvllm (2) The model is then estimated with carefully chosen initial values, and state smoothing and signal extraction proceeds as before: When p > 1 the output from signal is of dimension p × n × M, where M is the number of components. The level, regression and intervention are treated as one component of data level as before, separated from the seasonal component. The components estimated for the two series is shown in Figure 5. The intervention coefficient for the front passenger series is obtained by alphahat(25, n) = -0.30025, the next four elements are the coefficients of the regression of the two series on the price of petrol and kilometers traveled.

ARMA models
In this example the difference of the number of users logged on an internet server (Makridakis et al. 1998) is analyzed by ARMA models (Commandeur et al. 2011, Section 4) following Durbin and Koopman (2001, Section 9.4), and model selection via BIC and missing data analysis are demonstrated. To select an appropriate ARMA(p, q) model for the data various values for p and q are tried, and the BIC of the estimated model for each is recorded, the model with the lowest BIC value is chosen.  Table 1: BIC values of ARMA(p, q) models for users logged on an internet server.
Next missing data is simulated by setting some time points to NaN, model and state estimation can still proceed normally with missing data present. Forecasting is equivalent to treating future data as missing, thus the data set y is appended with as many NaN values as the steps ahead to forecast. Using the previous estimated ARMA(1, 1) model the Kalman filter will then effectively predict future data points.

Cubic spline smoothing
The general state space formulation can also be used to do cubic spline smoothing (Durbin and Koopman 2001, Sections 3.11 and 9.5), by putting the cubic spline into an equivalent state space form, and accounting for the continuous nature of such smoothing procedures. The discrete state space representation of the cubic spline model is given by Here the continuous acceleration data of a simulated motorcycle accident (Silverman 1985) is smoothed by the cubic spline model, where delta is the distance between the observations y, not necessarily regularly spaced. The smoothed data and standardized irregular are plotted and shown in Figure 7.
It is seen from Figure 7 that the irregular may be heteroscedastic, an easy ad hoc solution is to model the changing variance of the irregular by a continuous version of the local level model. Assume the irregular variance is σ 2 ε h 2 t at time point t and h 1 = 1, then we model the absolute value of the smoothed irregular abs(eps) with h t as its level. As defined in Harvey and Koopman (2000), the continuous local level model with level h t needs to be constructed from scratch. The smoothed data and standardized irregular with heteroscedastic assumption is shown in Figure 8, where it is seen that the confidence interval shrank, especially at the start and end of the series, and the irregular is slightly more uniform.
In summary the motorcycle series is first modeled by a cubic spline model with Gaussian irregular assumption, then the smoothed irregular magnitude itself is modelled with a local level model. Using the irregular level estimated at each time point as the relative scale of irregular variance, a new heteroscedastic irregular continuous time model is constructed with the estimated level built-in, and plugged into the cubic spline model to obtain new estimates for the motorcycle series.

Hillmer-Tiao decomposition
As an example of extensions to the SSM object library, TRAMO model selection (Gómez and Maravall 2001a) and Hillmer-Tiao decomposition for ARIMA type models as used in SEATS (Gómez and Maravall 2001b) are implemented. These ideas are widely used as part of the TRAMO/SEATS seasonal adjustment framework. The functions implemented here either take SSMODEL objects as arguments, or return SSMODEL objects.
In this example seasonal adjustment is performed by the Hillmer-Tiao decomposition (Hillmer and Tiao 1982) of airline models, a particular type of regComponent model (see Bell (2011) in this volume for a more in depth discussion of RegComponent models). Briefly, The ARIMA Model Based (AMB) decomposition for the airline model with seasonal period s, can be expressed in the following way using the decomposition of Hillmer and Tiao (1982) and ω t , η t , t are independent white noise processes and Here By t = y t−1 , with S t being the seasonal component, T t the trend and I t the irregular component, and θ S (B) and θ T B are seasonal and trend MA components, respectively. In order to define a unique solution (which may or may not exist) for θ S (B) and θ T (B), the ARIMA component parameters, in terms of the parameters of the original airline model, with MA parameter θ, seasonal MA parameter Θ, and variance σ 2 , the restriction is taken that the pseudo-spectral densities of the seasonal and trend components have a minimum of zero (in line with the admissible decompositions of Hillmer and Tiao 1982). When this condition cannot be met without resulting in a negative variance for I t , the decomposition is said to be inadmissible.
The data (manufacturing and reproducing magnetic and optical media, US Census Bureau) is fitted with the airline model, then the estimated model is Hillmer-Tiao decomposed into an ARIMA components model with trend and seasonal components. The same is done for the Frequency Specific ARIMA model , also known as the generalized airline model, and the seasonal adjustment results are compared. The following are the code to perform seasonal adjustment with the airline model: air = estimate(y, ssm_airline, 0.1); aircom = ssmhtd(air); ycom = signal(statesmo(y, aircom), aircom); airseas = ycom(2, :); aircom is the decomposed ARIMA components model corresponding to the estimated airline model, and airseas is the seasonal component, which will be subtracted out of the data y to obtain the seasonal adjusted series. ssmhtd automatically decompose ARIMA type models into trend, seasonal and irregular components, plus any extra MA components as permissible, which typically result from sARIMA models other than the airline model.
The same seasonal adjustment procedure is then done with the generalized airline model, using parameter estimates from the airline model as initial parameters: The code creates a generalized airline (3-5-1) model, Hillmer-Tiao decomposition produces the same components as for the airline model since the two models have the same order. From the code it can be seen that the various functions work transparently across different ARIMA type models. Figure 9 shows the comparison between the two seasonal adjustment results.

Non-Gaussian linear models 8.1. Poisson distribution error models
The road accident casualties and seat belt law data analyzed in Sections 3 and 4 also contains a monthly van driver casualties series. Due to the smaller numbers of van driver casualties the Gaussian assumption is not justified in this case, previous methods cannot be applied. Here a poisson distribution is assumed for the data (Durbin and Koopman 2001, Section 14.2), the mean is exp(θ t ) and the log density of the observation y t is log p(y t |θ t ) = θ t y t − exp(θ t ) − log y t ! Model estimation using the function estimate will automatically calculate the Gaussian approximation to the poisson model. Since this is an exponential family distribution the data y t also need to be transformed toỹ t , a time-varying transformation based on an appropriate approximating model (see Durbin and Koopman 2001, Table 11.1), which is stored in the output argument output, and used in place of y t for all functions implementing linear Gaussian (Kalman filter related) algorithms. The following is the code for model and state estimation: [pbstsm logL output] = estimate(y, pbstsm, 0.0006, [], 'fmin', 'bfgs', 'disp', 'iter'); [alphahat V] = statesmo(output.ytilde, pbstsm); thetacom = signal(alphahat, pbstsm); Note that the original data y is input to the model estimation routine, which also calculates the transform ytilde. The model estimated then has its Gaussian approximation built-in, and will be treated by the state smoother as a linear Gaussian model, hence the transformed data ytilde needs to be used as input. The loglikelihood logL here is based on importance sampling using Gaussian approximation of the model, with 1000 samples (default setting). The state estimates alphahat are calculated from both the approximation and original models. The signal components thetacom obtained from the smoothed state alphahat is the components of θ t , the mean of y can then be estimated by exp (sum(thetacom, 1)).
The exponentiated level exp(θ t ) with the seasonal component eliminated is compared to the original data in Figure 10. The effect of the seat belt law can be clearly seen.

t distribution models
In this example another kind of non-Gaussian models, t distribution models is used to analyze Since t distribution is not an exponential family distribution, data transformation is not necessary, and y is used throughout. A plot of the irregular in Figure 11 shows that potential outliers with respect to the Gaussian assumption have been detected by the use of heavy-tailed distribution.

Concluding remarks
SSM is an open source MATLAB toolbox for data analysis using state space methods. It combines the usability and flexibility of interactive model construction with the efficiency of C-implemented state space algorithms. In addition to the large number of predefined models available for immediate use, any models expressible in state space form, with appropriate approximations if they are nonlinear or non-Gaussian, can be defined and used within SSM, without having to write additional code for optimization, estimation or combination with other models. While extensive support for MCMC or other sampling based approximation techniques is not currently implemented beyond that of simulation smoothing, these could be combined using other MATLAB implementations of such methods.
By adapting open source licenses, it is hoped that the software package can attain wide-spread use (as of September 2010 well over 2000 downloads of the software have been made since SSM went online on October 2007), with user contribution of new models and even new state space algorithms. Possible future features under consideration include incorporation of particle filters to analyze more non-Gaussian or heteroscedastic model types, and incorporation of Markov switching for general regime or change-point analysis. On the software side, usability will continue to be improved, while it is also important to consider porting to a completely open source environment such as R or Python.

Software availability
The software is freely available from http://sourceforge.net/projects/ssmodels/ along with detailed documentation and numerous demo files, including full code and instructions for performing the analysis of the examples given in the paper. The toolbox was written for MATLAB 7.0 R14 and later, and any C compiler compatible with MATLAB can be used.

A. Predefined model reference
The predefined models can be organized into two categories, observation disturbance models and normal models. The former contains only specification of the observation disturbance, and is used primarily for model combination; the latter contains all other models, and can in turn be partitioned into structural time series models, ARIMA type models, and other models.
The following is a list of each model and function name.
Observation disturbance models -Gaussian noise: ssm_gaussian ([p, cov]) or ssm_normal ([p, cov]) p is the number of variables, default 1. cov is true if they are correlated, default true.
-Poisson error: ssm_poisson -Binary error: ssm_binary -Binomial error: ssm_binomial(k) k is the number of trials, can be a scalar or row vector.
-Negative binomial error: ssm_negbinomial(k) k is the number of trials, can be a scalar or row vector.
-Exponential error: ssm_exp -Multinomial error: ssm_multinomial(h, k) h is the number of cells. k is the number of trials, can be a scalar or row vector.
t-distribution noise: ssm_t([nu]) nu is the degree of freedom, will be estimated as model parameter if not specified.
- x is a m × n matrix, m is the number of regression variables. varname is the name of the variables.
-Intervention components: ssm_intv(n, type, tau) n is the total time duration. type can be 'step', 'pulse', 'slope' or 'null'. tau is the onset time.
-Length-of-month variables: ssm_lom(y, m, N) -Leap-year variables: ssm_ly(y, m, N) -Easter effect variables: ssm_ee(y, m, N, d) y is the starting year. m is the starting month. N is the total number of months. d is the number of days before Easter.
-Structural time series models: ssm_stsm(lvl, seas, s[, cycle, x]) lvl is 'level' or 'trend'. seas is the seasonal type (see seasonal components). s is the seasonal period. cycle is true if there's a cycle component in the model, default false.
x is explanatory (regression) variables (see regression components).
-Common levels models: ssm_commonlvls(p, A_ast, a_ast) p is the number of variables (length of y t ).
A_ast is A * , a (p − r) × r matrix. a_ast is a * , a (p − r) × 1 vector. s is the period, default 12. -Frequency specific SARIMA models: ssm_freqspec(p, d, q, P, D, nparam, nfreq, freq) -Frequency specific airline models: ssm_genair (nparam, nfreq, freq) nparam is the number of parameters, 3 or 4. nfreq is the size of the largest subset of frequencies sharing the same parameter. freq is an array containing the members of the smallest subset. -ARIMA component models: ssm_arimacom (d, D, s, phi, theta, ksivar) The arguments match those of the function htd, see its description for details.
Other models -Cubic spline smoothing (continuous time): ssm_spline(delta) delta is the time duration of each data point.
-1/f noise models (approximated by AR): ssm_oneoverf(m) m is the order of the approximating AR process.

B. Data analysis functions reference
Most functions in this section accepts analysis settings options, specified as option name and option value pairs (e.g., ('disp', 'off')). These groups of arguments are specified at the end of each function that accepts them, and are represented by opt in this section.  [, param0, alpha0, opt]) estimates the parameters of model starting from the initial parameter value param0. y is the data of dimension p × n, and model is a SSMODEL. param0 can be empty if the current parameter values of model is used as initial value, and a scalar param0 sets all parameters to the same value. Alternatively param0 can be a logical row vector specifying which parameters to estimate, or a 2 × w matrix with the first row as initial value and second row as estimated parameter mask. The initial state sequence estimate alpha0 is needed only when model is non-Gaussian or nonlinear. output is a structure that contains optimization routine information, approximated observation sequenceỹ if non-Gaussian or nonlinear, and the AIC and BIC of the output model. y is the data of dimension p × n, and model is a SSMODEL. The output is respectively the filtered state (m × n+1), filtered state variance (m × m × n+1, or 1 × n+1 if m = 1), one-step prediction error (innovation) (p × n), one-step prediction variance (p × p × n, or 1 × n if p = 1).
linear [model ytilde] = LINEAR(y, model [, alpha0, opt]) calculates the linear approximation. y is the data of dimension p × n, and model is a SSMODEL. alpha0 is the initial state sequence estimate and can be empty or omitted.
loglik LOGLIK(y, model[, ytilde, opt]) returns the log likelihood of model given y. y is the data of dimension p × n, and model is a SSMODEL. ytilde is the approximating observationỹ and is needed for some non-Gaussian or nonlinear models. generates observation samples from model conditional on data y. y is the data of dimension p × n, and model is a SSMODEL. antithetic should be set to 1 if antithetic variables are used. The output is respectively the sampled state sequence (m × n × N), sampled observation disturbance (p × n × N), and sampled state disturbance (r × n × N).
statesmo [alphahat V] = STATESMO(y, model[, opt]) performs state smoothing. y is the data of dimension p×n, and model is a SSMODEL. The output is respectively the smoothed state (m × n), smoothed state variance (m × m × n, or 1 × n if p = 1). If only the first output argument is specified, fast state smoothing is automatically performed instead.

oosforecast
[yf err SS] = OOSFORECAST(y, model, n1, h) performs out-of-sample forecast. y is the data of dimension p × n, model is a SSMODEL, n1 is the number of time points to exclude at the end, and h is the number of steps ahead to forecast, which can be an array. The output yf is the forecast obtained, err is the forecast error, and SS is the forecast error cumulative sum of squares.