State Space Methods in Ox/SsfPack

The use of state space models and their inference is illustrated using the package SsfPack for Ox. After a rather long introduction that explains the use of SsfPack and many of its functions, four case-studies illustrate the practical implementation of the software to real world problems through short sample programs. 
The first case consists in the analysis of the well-known (at least to time series analysis experts) Nile data with a local level model. The other case-studies deal with ARIMA and RegARIMA models applied to the (also well-known) Airline time series, structural time series models applied to the Italian industrial production index and stochastic volatility models applied to the FTSE100 index. In all applications inference on the model (hyper-) parameters is carried out by maximum likelihood, but in one case (stochastic volatility) also an MCMC-based approach is illustrated. Cubic splines are covered in a very short example as well.


Introduction
SsfPack (Koopman, Shephard, andDoornik 1999, 2008) is a library of routines for state space modelling and inference written in C and linked to Ox, the efficient matrix language developed by Doornik (2007) 1 .SsfPack is written by Siem Jan Koopman, one of the most prolific researchers and developers of state space algorithms, and comes in two different versions.One is free for academic use (the same policy as Ox console) and it is referred to as SsfPack Basic, while the other is commercial and its name is SsfPack Extended.The latter version includes the functionalities of SsfPack Basic plus a set of algorithms that are computationally 1 At the moment of writing this article the current version of SsfPack is 3.0, while that of Ox is 5.1.For SsfPack refer to the web page http://www.ssfpack.com/,while for Ox the relevant site is http://www.doornik.com/ more efficient and allow for exact treatment of the diffuse conditions for the initial state vector.Here only SsfPack Basic will be discussed, but the reader may be interested in knowing that for benefitting from the higher computational efficiency of SsfPack Extended the state space representation of the model must have orthogonal measurement errors (H t diagonal) and the disturbances in the transition equations must be uncorrelated with the measurement errors (Eη t ε t = 0).The book by Koopman et al. (2008) contains the complete documentation for both versions of the package, but the reader only interested in SsfPack Basic may also refer to Koopman et al. (1999).
SsfPack can handle both time-homogeneous and time-varying state space models, and contains various routines for (Kalman) filtering, prediction and smoothing, functions for computing the likelihood, in some cases with analytical scores, and algorithms for generating random sequences from the smoothed distribution of the state vector process.Furthermore, there are functions that facilitate the construction of the system matrices for some classes of frequently used models such as ARIMA, structural time series models (STSM), linear regressions and splines.
The package is usually installed in the folder package\ssfpack in the Ox installation directory, and throughout this article it will be assumed so.In order to use the package in an Ox program the line #include <packages/ssfpack/ssfpack.h> must be present in the program header 2 .
In order to lighten the exposition of the features of SsfPack, we introduce slight modifications to the state space notation used in this Journal of Statistical Software volume: for t = 1, . . ., n, y t = c t + Z t α t + ε t , ε t ∼ NID(0, H t ), where y t , c t and ε t are (p × 1) vectors, α t , d t and η t are (m × 1) vectors, Z t is the (p × m) observation matrix, T t is the (m × m) transition matrix, H t is (p × p), Q t is (m × m) and C t is (m × p).With respect to the notation used in this special issue, we added the two system vectors c t and d t , imposed the restriction3 R = I m and assigned the symbol C t to the covariance matrix of η t with ε t .The system is completed by the initial state specifications α 1 ∼ N(a 1 , P 1 ), α 1 ⊥ ε t , α 1 ⊥ η t , for t = 1, 2, . . ., n.
2 Those not acquainted with Ox, should take notice of the following: • // denotes a line comment; • statements end with a semicolon; • & is used for pointers, that is, for passing a variable through its address; • if pointer is a pointer, the pointed variable can be read or assigned using the syntax pointer[0]; • matrix indexing begins with 0; • <1,2,3;4,5,6> defines a matrix constant: , switches to the next column and ; to the next row; • the operator ~joints two matrices by columns (side by side), while | stacks two matrices; • arithmetic and logic operators on matrices follow linear algebra conventions (e.g., * is matrix product), elementwise operators are preceded by a dot (e.g., .* is Hadamard product).

Formulating the state space system
SsfPack represents the state space system through one vector and three matrices: The matrices Φ t and Ω t must always be assigned, while δ t and Σ are optional.The default assignments for the latter two matrices are where k is a very large constant.
In the rest of the article we will refer to the Ox variables containing the system matrices as, respectively, mDelta, mPhi, mOmega and mSigma.
In SsfPack both homogeneous (i.e., time invariant) and inhomogeneous systems can be implemented.In the first case a state space must be defined through either the pair mPhi, mOmega, or the triplet mPhi, mOmega, mSigma, or the quartet mPhi, mOmega, mSigma, mDelta.Since the user is free to chose any one of these forms, we will use the notation {Ssf} to refer to any of the three state space system definitions.
Defining an inhomogeneous system is slightly more complex.For the matrices with timevarying elements additional index matrices of the same dimensions have to be defined.For these matrices, we will use the self-explanatory notation mJ_Phi, mJ_Omega, mJ_Delta.Furthermore, a matrix with the time-varying elements in its rows, say mXt, has to be assigned.The elements of the index matrices are all set to −1 except those elements for which the corresponding elements in Φ t , Ω t , δ t are time varying.For those elements that change with time the row numbers of mXt, in which the needed time-varying elements are stored, must be indicated in the corresponding position of the relative index matrix.
For example, suppose that only the first element of a (3 × 2) Φ matrix is time varying, and the time varying coefficient is stored in the third row of mXt.Then, the index matrices have to be declared as4 mJ_Phi = <2, -1; -1, -1; -1, -1>; mJ_Omega = <>; mJ_Delta = <>; where <> is the empty matrix.As the reader may have noticed from the example, index matrices relative to time-homogeneous system matrices can be set to <>.
The four alternative ways to define a state space model in SsfPack are summarized below: mPhi, mOmega mPhi, mOmega, mSigma mPhi, mOmega, mSigma, mDelta mPhi, mOmega, mSigma, mDelta, mJ_Phi, mJ_Omega, mJ_Delta, mXt The data matrix, that we will denote as mYt, must have the series in rows and the time points in columns.Thus, using the above indexing notation, mYt is a (p × n) matrix.SsfPack has no problem to work with data matrices with missing values.

Generating the system matrices for common classes of models
The system matrices of a state space model can be directly assigned by the user, but for some common classes of models there are functions that simplify their assignment.
ARIMA models One of the most used classes of time series models is the ARMA family.
If a non-seasonal stationary and invertible ARMA(p,q) model is needed, the function to be used is (vAr, vMa, dStDev, &mPhi, &mOmega, &mSigma) where vAr and vMa are vectors containing, respectively, the AR and MA coefficients, dStDev is the standard deviation of the white noise processes and the system matrices are passed through their addresses.This function implements the ARMA model in the form
In addition to the above function, SsfPack Extended provides GetSsfSarima for seasonal integrated ARMA models (SARIMA).The user of SsfPack Basic can design the SARIMA system matrices directly, or take advantage of GetSsfArma by opportunely passing the ARMA parameters for a multiplicative seasonal model and modifying the output matrices for integrating the process.As example, the ubiquitous Airline model for monthly data will be discussed in Section 3.
Structural time series models Another important class of models both for social and natural sciences is represented by the structural time series models (STSM), also referred to as unobserved components (UC) time series models.For a short introduction to these models the reader can refer to the introductory article of this special issue.The books by Harvey (1989), West and Harrison (1997), Durbin and Koopman (2001) and Commandeur and Koopman (2007) are all devoted to these kind of models and the interested reader should refer to them for details.

GetSsfStsm
Notice that CMP_LEVEL, CMP_SLOPE, CMP_SEAS_DUMMY etc. are predefined constants, the notation CMP_SEAS_DUMMY implies a stochastic dummy specification of the seasonality and is alternative to CMP_SEAS_TRIG, which defines a (stochastic) trigonometric seasonal component.Up to ten stochastic cycles may be used in the same model.The ordering of the components specification may be different, but the system matrices are built according to the ordering given above (trend, slope, seasonality, cycles, irregular).
The second column of the matrix contains the standard deviations of the components' disturbance with one exception: the standard deviations of the cycle components refer directly to the square root of the cycle variance σ 2 ψ = σ 2 κ /(1 − ρ 2 ), where σ 2 κ is the variance of the cycle disturbance.
The third column of the matrix is used only for the seasonal and cyclical components and, in the former case, s is the seasonal periodicity (e.g., 12 for monthly data and 4 for quarterly series), while in the latter case λ c is the modal frequency of the cycle (e.g., if the modal cycle length is 60 months the corresponding frequency is given by λ c = 2π/60).
The fourth column is used only in the specification of cycle components and contains the damping factors ρ ∈ (0, 1).The fifth column is not used in the SsfPack Basic implementation, while in SsfPack Extended it is used to specify higher order cycles of the type proposed by Harvey and Trimbur (2003).
The use of GetSsfStsm is illustrated with an example in Section 4.
Regression models Regression models can be easily represented in inhomogeneous state space form as where x t is a vector of regressors at time t and α t is a vector of regression coefficients.A regression with time-varying coefficients may be obtained by changing the state equations to α t+1 = α t + η t .Unless better information is available, the state vector at time t = 1 is given a diffuse distribution.
The SsfPack function for building the relevant system matrices is (mXt, &mPhi, &mOmega, &mSigma, &mJ_Phi) where mXt is the (k × n) data matrix.Notice that the function uses only the row dimension of mXt and not the content.

GetSsfReg
Adding regressors to time series models Suppose the matrices mPhi, mOmega and mSigma have been properly assigned for a time-invariant state space model.If regressors have to be added to the model, the function AddSsfReg(mXt, &mPhi, &mOmega, &mSigma, &mJ_Phi) can be used.Again, the function uses only the row dimension, say k, of mXt and not its content.
The function call produces where the system matrices T , Z, Q, H, P 1 and a 1 are extracted from the passed mPhi, mOmega and mSigma.The returned index matrix mJ_Phi is where −1 denotes a matrix of −1's with the indicated dimensions and i = (0, 1, . . ., k − 1).
In Section 3 it will be shown how to build a regression model with ARMA errors (generally referred to as RegARMA or ARMAX ) using the functions GetSsfArma and AddSsfReg jointly.
Nonparametric cubic splines A cubic spline is a function used in interpolation problems and for the approximation of smooth functions observed with noise (a form of nonparametric regression).The interested reader should refer to Green and Silverman (1994) for a general treatment on cubic splines and to Wecker and Ansely (1983) for their relation with signal extraction.
The SsfPack function for computing cubic splines is (dq, vxt, &mPhi, &mOmega, &mSigma, &mJ_Phi, &mJ_Omega, &mXt) where dq is the signal to noise ratio, vxt is a (row) vector of abscissa points for which we have noisy observations of the smooth function:

GetSsfSpline
The approximation to the function can be recovered by passing the system matrices generated by GetSsfSpline and the ordinate data y i to the smoothing functions discussed below.A case-study for this function is not covered in the text of this article, but the interested reader can take a look at the code in Spline.ox, in which a (random) continuous curve is generated and sampled at a finite number of points, then some noise is added and a cubic spline is fitted.

Filtering, smoothing, predicting and forecasting
SsfPack is extremely rich in routines for carrying out predictions, filtering and smoothing.
Here, we will just cover those functions useful to the typical "end-user" of state space models.
Recall that in a Gaussian state space system, conditionally on the observations, the state variables are normal and the Kalman filter and smoother are algorithms to compute the conditional moment pairs with Y s = {y 1 , . . ., y s }, which are named predictor, filter or smoother when, respectively, s < t, s = t, s > t.
The main function for estimating the first two conditional moments of the state variables is where iSel is a selection variable whose possible values are ST_SMO, ST_FIL, ST_PRED, DS_SMO, for moment smoothing, filtering, prediction and disturbance smoothing, respectively.The matrix mOutput, passed to the function by its address, will contain the main output of the function.Naming the first two conditional moments of the signal s t = c t + Z t α t as the matrix mOutput will have the following structure: where the subscript t|• must be substituted either with t|t − 1, t|t or t|n, and the diag operator returns the main diagonal of a square matrix as column vector.
For forecasting (multi-step out-of-sample predictions), the analyst can augment the data matrix with columns of missing values at the end of the observed sample, and then use the function SsfMomentEst for either smoothing, filtering or prediction.For example, for forecasting q steps ahead, if the matrix mYt contains the observed data, the following line of code could be used5 : mYt ~= constant(.NaN, p, q); where the Ox function constant generates a (p × q) matrix with all elements equal to its first argument, that in this case is the missing value (.NaN, not a number).

Likelihood evaluation
SsfPack offers the following three functions for evaluating the log-likelihood of a state space model: bSuccess = SsfLik(&dLogLik, &dVar, mYt, {Ssf}) bSuccess = SsfLikSco(&dLogLik, &dVar, &mSco, mYt, {Ssf}) bSuccess = SsfLikConc(&dLogLikConc, &dVar, mYt, {Ssf}) All functions return 1 if successful and 0 otherwise.The first two functions evaluate the log-likelihood, , and write its value to the variable dLogLik passed through its address.SsfLikSco computes also the matrix of scores with respect to the elements of the covariance matrix Ω, and assign it to mSco through its address.
Suppose that ψ i is an unknown parameter associated only with elements of the matrix Ω, then Koopman and Shephard (1992, equation 3.2) show that the score with respect to ψ i is given by By all functions the variable dVar is assigned the value σ2 = (np − d) −1 n t=1 ν t F −1 t ν t , where d is the number of state variables with diffuse initial conditions, ν t is the innovation vector and F t its covariance matrix.
The function SsfLikConc computes the profile log-likelihood where a scale parameter is concentrated out, reducing the dimensionality of the estimation problem by one parameter (see Harvey 1989, pp. 126-127).Notice that, while for SsfLik and SsfLikSco the variable dVar will be approximately equal to 1 after maximum likelihood (ML) estimation, for SsfLikSco dVar contains the scale parameter concentrated out of the likelihood.In order to make the examples in the next sections easily readable by readers who are not acquainted with Ox, a few lines are now dedicated to the maximization function MaxBFGS, upon which we rely for ML estimation of models in state space form.The optimization package is to be imported using #import <maximize> in the header of the program.The BFGS maximization function in Ox is iCode = MaxBFGS(func, &vP, &dFunc, 0, bNumDer) where func is the name of the objective function, vP is the parameter (column) vector (passed through its address) that contains the initial values at input and the value of the parameters at maximum at output, dFunc is a variable that will contain the maximum of the objective function at the end of the maximization process, and bNumDer is a boolean variable that should be set to 0 if func returns analytical derivatives and to 1 otherwise (numerical derivatives will then be computed).The fourth argument of the function is not used in this paper and will always be set to 0. MaxBFGS returns a number that can be passed to MaxConvergenceMsg(iCode) which returns a human-readable string for assessing the degree of success of the convergence process.
The function func is expected by MaxBFGS to have the following syntax: bSuccess = func (vP,&dFunc,&vScore,0) where vP is a variable through which the parameter vector is passed, in dFunc the function func must write the value of the objective function at vP, and if the third argument is different from 0, then func is expected to write the value of the gradient at vP in the variable vScore.The function must return 1 in case of success or 0 in case the evaluation failed.

Simulation
Since the advent of fast and cheap personal computers, simulation has become a standard tool in statistics.SsfPack is extremely rich in functions for generating random quantities implied by models in state space form.
For simulating artificial observations from the joint distribution implied by a state space model, the SsfPack function is mD = SsfRecursion(mR, {Ssf}) is a sequence of random vectors with covariance matrices Ω t to be generated by the user.For example, in the Gaussian case, i.e., for (η t ε t ) ∼ NID(0, Ω), the sequence for t = 1, . . ., n can be generated using the following line of Ox code: choleski(mOmega)*rann(m+p, n); Notice that, while ε 0 is irrelevant, since no y 0 vector will be generated, η 0 must have zero mean and identity covariance matrix, since α 1 is generated according to α 1 = a + Lη 0 , and the mean vector a and the covariance matrix P of α 1 must be specified in the initial moments matrix mSigma as follows: where L satisfies LL = P and can be computed with the function mL = choleski(mP).
The output mD contains a pseudo-random realization of the state variables process and of the observable variables process: mD = α 1 α 2 . . .α n+1 0 y 1 . . .y n .
Carrying There are two ways to generate samples from the conditional random processes above: one is more efficient, but more involved, the other one is straightforward to implement, but somewhat slower.Here, we cover just the latter method, and refer to the SsfPack manual (Koopman et al. 2008) for the former.
The function for easily drawing from the smoothed distributions is mD = SsfCondDens(iSel, mYt, {Ssf}) iSel can be set equal to ST_SIM or DS_SIM and the function returns a sample from, respectively,

Case 1: The local level model applied to the Nile data
In this section it will be shown how the local level model (LLM) can be fit to the volume of the Nile river at Aswan from 1871 to 1970, using maximum likelihood (ML) estimates for the two unknown variances in the model.At the end of the section the model will be extended to allow for a break in the trend.
The SsfPack functions treated are SsfLik, SsfLikSco and SsfMomentEst.The complete code can be found in the file LLM_Nile.ox,that makes use of the data in Nile.dat.
The code is organized in three functions plus the main().Furthermore, since the optimization functions do not allow to pass variables to the objective function other than the parameters with respect to which the optimization is carried out, a global (static) variable is needed: static decl s_vYt; The first function takes the two variances of the LLM (state disturbance first) as arguments and uses them to write the Φ and Ω system matrices to variables passed by address.
The second function is the objective function to be passed to the maximizer and, as noted above, must have a specific structure (see Section 1.4).Before proceeding with the discussion of the function, two things need to be emphasized.Firstly, since the MaxBFGS optimizer does not allow constraints on the parameter space, we pass the two variances in their (unconstrained) logarithms, and take the exponents of these log-variances to ensure their nonnegativity.This implies that the scores for the log-variances have to be obtained from the scores of the variances.So, if (σ 2 ) is the log-likelihood and σ 2 = exp(p) is the variance, than the score with respect to p is given by where ∂ (σ 2 )/∂σ 2 is computed by SsfLikSco.Secondly, for numerical reasons we prefer to return (through the pointer adLogLik) the mean log-likelihood.Of course, this choice does not affect the estimates.
Notice that in the function below the variable/pointer amHessian is not used, and avScore will contain the (mean) score if an address is passed, while if equal to 0 only the (mean) log-likelihood will be computed.The code is rather self-explanatory and illustrates the use of both SsfLik and SsfLikSco.
The third function, acov(const vLogP) (not reported in the text) computes estimates of the asymptotic covariance matrix of the ML estimators.This function will be used in many other short programs.
The last function is the main(), which loads the data, calls the maximization routine, prints the estimates and plots the graphs.

.] }
Notice the different features of the function SsfMomentEst used for smoothing both the state variables and their disturbances, getting one-step-ahead predictions and forecasts.
Since the optimization has been carried out with respect to log-variances, standard errors for the variances have been obtained using the delta method6 .By running LLM_Nile.oxthe following estimates are obtained, while the graphs are shown in Figure 1.

3139.1
As can be noticed by looking at the auxiliary residuals for the trend and relative CI in Figure 1, a break occurred at the end of the XIX century.
Atkinson, Koopman, and Shephard (1997, Section 6.3) suggest a break between 1897 and 1900 as well.The file LLM_Nile_Break.oxcontains a version of the LLM seen above, where a further trend-break parameter is to be estimated for a date fixed by the user.
We discuss the main modifications made to the previous code.Now there are two global variables: static decl s_vYt; // data static decl s_iIndex; // index of the break The function that builds the (now inhomogeneous) system matrices is   The other functions (not reported) are not too different, even though in the loglik function we show how to compute analytical scores with respect to the variance and numerical scores with respect to the break impact.The textual output for a break in year 1897 is:

2292.5
Thus, when a break is allowed in the LLM, the trend becomes deterministic since its disturbance variance is virtually zero.The same result was found by Atkinson et al. (1997).
Figure 2 illustrates this finding with a line plot.

Case 2: The airline model applied to the airline time series
State space models are often used for carrying out exact ML estimation of ARIMA models.
The following function computes the log-likelihood at the parameters' value in vP.Notice that only the two MA parameters are to be passed since the innovation variance is concentrated out by the function SsfLikConc and assigned to the global variable s_dVar.Indeed, the function airline is called with arbitrary (unitary) innovation standard deviation.loglik(const vP, const adLogLik, const avScore, const amHessian) { decl cT = columns(s_vYt); // Number of observations decl iRetCode, dLogLik, dVar, mPhi, mOmega, mSigma; airline(vP|1, s_cPeriod, &mPhi, &mOmega, &mSigma); iRetCode = SsfLikConc(adLogLik, &s_dVar, s_vYt, mPhi, mOmega, mSigma); adLogLik[0] /= cT; // Mean loglik return iRetCode; } The next function produces forecasts given the model's parameters, the forecast horizon and the desired confidence level in (0, 1).Moreover, if the airline model is built for data in log, the function produces forecasts (under log-normality) and confidence bands for the original data.We just mention that the function regressors(cPeriod, cT, &mXt) in the code fills the matrix pointed by &mXt (by row) with a constant, a time-counter and cPeriod-1 seasonal sinusoids.In cPeriod the seasonal periodicity integer must be passed, while the integer in cT represents the number of time points to be generated.
The function that computes the log-likelihood value for the ARMA parameters in the vector vP is reported below.The function uses the global variables mentioned above to retrieve the orders of the AR and MA parts.Since the regression coefficients are specified as (constant) state variables, they are automatically concentrated out of the likelihood.In our implementation, the sum of the AR and MA orders must be at least one (otherwise no ARMA is specified).

Airline
Regression 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957 1958 1959  We do not report the main function, but we suggest the reader to go through the code therein to see how the inference on the regression coefficients is obtained through the SsfMomentEst function.The code produces Figure 4 and the estimates reported below.

Case 3: A STSM for the Italian industrial production
In this section we show how to fit a complete structural time series model to real datathe Italian industrial production index (IPI) -using the SsfPack function GetSsfStsm.By complete STSM it is meant that the data generating process is the sum of a local linear trend, estimates and smoothed components seem quite reasonable and interpretable.For example, the modal frequency of the cycle (0.1678)corresponds to a period of 37 months, which is not surprising for the Italian industrial production.

Case 4: The stochastic volatility model
The aim of this last case-study is to show how simple it is for an SsfPack user to implement MCMC methods in a non-Gaussian state space framework.The function we focus on is SsfCondDens that generates from the state vector distribution conditionally on the observations.
In financial econometrics there are two popular classes of models that deal with time series with time-varying variances: the class of GARCH-type models, where the variance h t of the process y t is measurable with respect to the information set F t−1 generated by the past of y t , and the stochastic volatility (SV, also stochastic variance) models, in which the variance follows a process of its own, which is not measurable with respect to F s for any s.
The simplest specification of a SV model is where the two disturbances are independent.This model is a nonlinear state space form, but can be easily linearized by taking the log of the square of the observable time series: Now the state space is linear, but not Gaussian, since the observation error is a log χ 2 1 distribution.Harvey, Ruiz, and Shephard (1994) propose to approximate the observation error distribution with a normal with the same first two moments of log χ 2 1 : log y 2 t = −1.27+ h t + ε t , ε t ∼ NID(0, π 2 /2).The Gaussian likelihood obtained through prediction error decomposition for this model becomes a pseudo-likelihood (also quasi-likelihood) and the Kalman filter and smoother are no more conditional expectations, but only linear projections.Thus, the relative inference provides consistent although inefficient estimators.
The code in SV_Harvey_ftse.ox,that makes use of the data in ftse100.xls,implements the simple SV model for the daily returns of the FTSE100 index.In order to reduce the dimension of the parameter vector with respect to which the likelihood is maximized, the model has been implemented in the equivalent state space form log y 2 t = −1.27 where the constant µ = µ t plays the role of the marginal expectation of h t and h t = h * t + µ.The Ox/SsfPack code for this analysis is very simple and so similar to the previous examples that we do not discuss it in the text.The returned QML estimates are µ = −0.023,φ = 0.925 and σ 2 = 0.289.Figure 6 plots the FTSE100 absolute percentage returns together with their smoothed standard deviations.Kim, Shephard, and Chib (1998) consider the same model, but propose a different estimation technique based on Gibbs sampling.Their approach is excellent for demonstrating the simulation facilities of SsfPack.The code implementing a slight modification of their Gibbs sampler that we are about to discuss is in the file SV_Kim_ftse.ox,and uses the data in ftse100.xls.In Section 3 of their paper, Kim et al. (1998) provide a seven-component Gaussian mixture that approximates the log χ 2 1 distribution very well.In our code the mixture is defined as the (static) global variable s_mMix, a (7 × 3) matrix, in which the first column contains the component probabilities, and the second and the third columns contain, respectively, mean and variance of each component.
The Gibbs sampling recursion for estimating the parameters and the unobservable variance of the SV model consists of the following three steps: (this is achieved using the SsfPack function SsfCondDens);

given the previously sampled {h
and generate first σ 2 η and than β using their posterior distributions under noninformative priors: where IG(a, b) denotes the Inverse Gamma distribution with parameters a and b.
Before going through the SV_Kim_ftse.oxroutines that carry out the above steps, it is useful to notice that the following auxiliary functions are therein coded and used: phi(mX, mMu, mSigma) returns a matrix containing the values of the Gaussian densities of the points in mX computed for mean and standard deviations, respectively, in mMu and mSigma; randiscrete(mProbs) returns a (1×n) vector of random numbers chosen from 0, 1, . . ., k−1, according to the probabilities specified in the (k ×n) matrix mProbs (every column sums to 1); ranmvn(vM, mS, c) returns a matrix (k × c) of c random vectors generated from a normal distribution with mean (k × 1) vector vM and covariance (k × k) matrix mS.
The function that accomplishes step 1 of the Gibbs sampler is w_gen(const vH, const vLogY2) that generates the indices of the mixture given the variances (vH) and the transfomed observations (vLogY2).Step 3 is carried out by arpars_gen(const vH) that generates the autoregression parameters given the log-variances in vH.
The Gibbs sampler iterations are implemented in the function sv_sim(const vLogY2, const cN, const amH, const amArPars, const amW), where the transformed data are passed through vLogY2, the number of simulations is specified in cN and the other three parameters are pointers to variables that will host the simulated log-variance time series, the AR parameters and the mixture weights time series, respectively.
In principle, the initial values for the Gibbs sampler could be arbitrary, since, eventually, the ergodic Markov chain will start exploring the sample space with the right marginal probabilities, but initial values in a probability-dense neighbor of the sample space will speed up this convergence.In the function sv_sim we chose to fix the initial {h t } through some kind of exponential smoothing of log y 2 t and then generate the AR parameters and weights using the steps 3 and 1 of the Gibbs sampler.Generally a burn-in sample of simulations is discarded in order to make the Markov chain "forget" the initial values.
The main function loads the data, calls the Gibbs sampler (sv_sim) and, after 21000 iterations, the first 1000 of which are discarded, produces the plot in Figure 7 and the following textual output.The smoothed variance estimates in Figure 7 are obtained as averages of the 20000 Gibbs samples for exp{h t /2}.By taking sample quantiles of the same quantities, confidence intervals (or credible intervals if we take a Bayesian approach) can be computed.The main difference between these results and our pseudo-ML estimates is the higher persistence of the log-variance process h t .This can be seen both from the larger value of the AR(1) coefficient φ and in the smoothed variance graphs.

Mean
In order to assess if the Gibbs sample is large enough to sufficiently explore the parameter space the code in SV_Kim_ftse.oxproduces also autocorrelation functions and kernel density estimates of the MCMC samples (not reported here).

Conclusions
In this tour trough the Ox package SsfPack, we visited just the higher-level functions which are those that the typical user will utilize.Advanced users will find a host of other functions for simulation and filtering/smoothing purposes that can make these operations even more efficient and accurate (see Durbin and Koopman 2001;Koopman et al. 2008).
In particular, the commercial version of the software, SsfPack Extended, implements exact algorithms to deal with diffuse initial conditions for the state variables and efficient routines to carry out computations in systems with large matrices, in some cases by exploiting the fact that these are typically sparse (i.e., contain many zeros).
Furthermore, SsfPack Extended works under Windows 32-bit, Windows 64-bit, OS X, Linux 32-bit and Linux 64-bit, while SsfPack Basic is only for Windows 32-bit.
If compared to other software for state space modelling we are acquainted with, SsfPack is faster and wider-ranging.On the other hand, something that is missing in SsfPack, but present in software written with engineering in mind, such as MATLAB, is some implementation of the Extended Kalman Filter (EKF), that is, the possibility of letting the system matrices depend on predicted state variables.
Any researcher that makes extensive use of models in state space form can only benefit from trying SsfPack, and after some use he/she will probably find it irreplaceable for the speed and stability of its algorithms.
Information, updates and other example programs can be found at the official SsfPack Internet site http://www.ssfpack.com/,which is maintained by S. J. Koopman, the author of SsfPack.

Figure 1 :
Figure 1: Nile with smoothed trend, standardized prediction errors, auxiliary residuals for the observation equation (for outliers), auxiliary residuals for the trend (for breaks), forecasts.

Figure 4 :
Figure 4: Regression fit and smoothed ARMA errors in the RegARMA model fitted to the airline data.

Figure 5 :
Figure 5: Smoothed components of the Italian industrial production index.

Figure 6 :
Figure 6: FTSE100 absolute percentage returns and smoothed standard deviation based on pseudo-ML estimation.

Figure 7 :
Figure 7: FTSE100 absolute percentage returns and smoothed standard deviation based on MCMC estimation.
out Bayesian inference using Markov Chain Monte Carlo methods often requires generating sample values from the joint distribution of the state variables conditional on the observed data.If we set Y s = {y 1 , . . ., y s }, SsfPack has functions to generate from α t |Y n and η t |Y n , as well as from s t |Y n and ε t |Y n .
The code we are about to illustrate is contained in the file Airline.ox and uses the data in airline.xls.The global variables needed in the rest of the code are declared as The SsfPack functions treated are GetSsfArma, SsfLikConc, SsfMomentEst and AddSsfReg.