State Space Modeling Using Sas

This article provides a brief introduction to the state space modeling capabilities in SAS, a well-known statistical software system. SAS provides state space modeling in a few different settings. SAS/ETS, the econometric and time series analysis module of the SAS system, contains many procedures that use state space models to analyze univariate and multivariate time series data. In addition, SAS/IML, an interactive matrix language in the SAS system, provides Kalman filtering and smoothing routines for stationary and nonstationary state space models. SAS/IML also provides support for linear algebra and nonlinear function optimization, which makes it a convenient environment for general-purpose state space modeling.


Introduction
SAS software (SAS Institute Inc. 2008a) offers a set of solutions for enterprise-wide business users and provides a powerful fourth-generation programming language for performing tasks such as data management, report writing and graphics, statistical and mathematical analysis, and operations research.This article provides a brief introduction to the state space modeling capabilities that are available in SAS/ETS (SAS Institute Inc. 2010), the econometric and time series analysis module of the SAS system, and in SAS/IML (SAS Institute Inc. 2008b), the SAS interactive matrix language.The article is organized as follows.Section 2 uses the SAS/ETS UCM procedure to analyze the Nile data, the common example used in all articles in this volume.Section 3 discusses additional SAS/ETS procedures that are useful for state space modeling.Section 4 illustrates the use of SAS/IML for state space modeling.The final section contains concluding remarks.

Unobserved components modeling
The SAS/ETS UCM procedure can be used to analyze and forecast univariate time series data using an unobserved components model (UCM).A UCM decomposes the response series into components such as trend, seasonals, and the regression effects due to predictors.The next section uses the Nile data to illustrate some of these aspects of unobserved components modeling.

Analysis of Nile data
The Nile data consist of yearly readings of the Nile river water level, from 1871 to 1970, recorded at Aswan, Egypt.The analysis of this series begins by first fitting to it the local level model: That is, the readings y t are noisy observations of the underlying level µ t , which follows a random walk (RW).The observation noise t and the RW disturbance ξ t are assumed to be independent, Gaussian, white noise sequences with variance parameters σ 2 and σ 2 ξ , respectively.The local level model is an example of a UCM that decomposes the response series into two unobserved components: the level component µ t and the irregular component t .It is a natural starting point in this analysis because the yearly water level of a large river, in the absence of major external shocks, can be expected to remain relatively constant for a long time.You can fit this model to the Nile data by using the UCM procedure as follows: proc ucm data = nile; id year interval = year; model waterlevel; irregular plot = smooth; level checkbreak plot = smooth; estimate plot = residual; forecast plot = forecasts lead = 10 alpha = 0.5; run; The PROC UCM statement signifies the start of the UCM procedure and specifies the input data set, nile, which contains the response series.The optional ID statement specifies a date, time, or datetime variable (year in this example) to label the observations.The INTERVAL = YEAR option in the ID statement indicates that the measurements are collected on a yearly basis.The model specification begins with the MODEL statement, where the response series is specified (waterlevel in this case).The components in the model are specified using separate statements that enable you to control their individual properties.The irregular component, t , is specified using the IRREGULAR statement, and the level component, µ t , is specified using the LEVEL statement.The options in the ESTIMATE statement control the model-fitting environment and the diagnostics output-for example, Figure 1 shows the residual plot produced by using the PLOT = RESIDUAL option.The parameters of this model are the variances σ 2 and σ 2 ξ ; Table 1, produced by default, shows their maximum likelihood estimates.The options in the individual component statements serve different purposes-some determine their stochastic properties while others can be used to control their output.For example, Figure 2 shows the plot of the smoothed estimate of level µ t , which is obtained by using the PLOT =   SMOOTH option in the LEVEL statement.The plot seems to indicate a shift in the water level starting at or around 1899.This shift could be attributed to the start of construction of a dam near Aswan in that year.Here it is known that a shift in the water level occurred within the span of the series.However, in many cases such prior information is not available and it is useful to detect such shifts in a data analytic fashion.significant.
An additive outlier, an unusually low water level in 1913, is also detected.The outlier summary reported in Table 2 can also be inferred from the plot of smoothed, standardized auxiliary residuals shown in Figure 3.The series forecasts and the series decomposition can be requested using the FORECAST statement; the forecast plot shown in Figure 4 is produced by the PLOT = FORECASTS option.
The local level model can be corrected for the level shift of 1899 by including a regression variable, say ls1899, that is zero before 1899 and one thereafter.The additive outlier at 1913 can be similarly incorporated in the model by including a variable that is one at 1913 and zero elsewhere.In order to keep the model structurally simple the following program corrects the model for the level shift only.The level shift variable (ls1899) is specified as a regressor variable in the MODEL statement.
proc ucm data = nile; id year interval = year; model waterlevel = ls1899; irregular; level; estimate plot = panel; forecast plot = decomp lead = 10 alpha = 0.1; run; Table 3 shows the parameter estimates of this revised model.Note that the estimate of the noise variance σ 2 has not changed much.However, the estimate of the level disturbance variance σ 2 ξ has become nearly zero, indicating that the level component, after accounting   for the level shift of 1899, is nearly time-invariant.This is also evident from Figure 5, which shows the plot of the sum of smoothed trend and regression effect.In conclusion, you can say that the local level model with the shift of 1899 is a reasonable description for the Nile data.

UCM procedure: Additional information
The preceding section provided a simple illustration of UC modeling using the UCM procedure.You can use the UCM procedure to fit a wide range of UCMs that can incorporate complex trend, seasonal, and cyclical patterns and can include multiple predictors.For example, you can analyze most univariate models discussed in Harvey (1989)-a classic reference for UC modeling.Moreover, more elaborate models such as the one proposed for modeling hourly electricity demand in Harvey and Koopman (1993) and a model proposed for seasonal adjustment of weekly data in Harvey, Koopman, and Riani (1997) can also be easily handled.Specifically, you can specify UCMs that consist of the following features: a variety of trend models, including random walk and local linear trend, several types of seasonals, including trigonometric seasonals with full control over the included harmonics and spline approximation for long seasonal patterns, multiple stochastic cycles, regression effects with time-invariant and time-varying regression coefficients, nonlinear regression effects that are accomplished using spline specification, an irregular component that can incorporate serial correlation modeled by a two-factor stationary ARMA process (ARMA(p, q)(P, Q) s ), models that include lagged response values -for example, an ARIMA(p, d, q)(P, D, Q) s model.
The UCM procedure properly handles missing values in the response series.In particular, it provides a rigorous treatment of ARIMA models with missing data (see Kohn and Ansley 1986).
The UCMs considered in the UCM procedure are special cases of the (linear) Gaussian state space models (GSSM).The analysis of these GSSMs is based on the Kalman filtering and smoothing (KFS) algorithm.If the initial state of a GSSM is diffuse, as is often the case for UCMs, its treatment requires modification of the traditional KFS.The modification is called the diffuse KFS (DKFS).The DKFS implementation in the UCM procedure is based on de Jong (1991).The outlier detection process implemented in the UCM procedure is based on de Jong and Penzer (1998).For additional information, see the chapter "The UCM Procedure" in SAS Institute Inc. (2010).

State space modeling in SAS/ETS
Many state space models can be formulated as univariate or multivariate autoregressive moving average (ARIMA) models; for examples see Brockwell andDavis (1991) andReinsel (1997).The SAS/ETS ARIMA and VARMAX procedures provide state-of-the-art facilities for univariate and multivariate ARIMA modeling, respectively.
The SAS/ETS STATESPACE procedure can be used for state space modeling of multivariate time series, including automatic model identification based on the canonical correlation analysis.
The canonical correlation analysis implemented in the STATESPACE procedure is based on Aoki (1990).This procedure, however, is not a general purpose state space modeling procedure.In particular, it cannot be used to fit multivariate structural models discussed in Harvey (1989).

State space modeling using SAS/IML
SAS/IML is an interactive matrix language that, among other things, provides support for linear algebra and nonlinear function optimization.It also provides Kalman filtering and smoothing routines for stationary and non-stationary state space models.
Consider the following state space model with a partially diffuse initial state: where z t is the state vector, y t is the observation vector, and for t ≥ 0, The noise vectors (η t , t ) and the (diffuse) random vector δ are assumed to be mutually independent.The system matrices, W t , F t , X t , H t , a, A, b, B, V t , G t , and R t , are assumed to be known.The analysis of state space models of this form requires modification of the traditional Kalman filtering and smoothing (KFS) algorithm, which is called the diffuse KFS.You can analyze state space models of this form by using the SAS/IML subroutine Introductory paper SAS/IML state space model Table 4: Notation for state space models.
KALDFF for diffuse Kalman filtering, and by using the subroutine KALDFS for diffuse Kalman smoothing.The SAS/IML DKFS implementation is based on de Jong (1991).Note that the model described above can be recast into the canonical state space model described in the introductory paper Commandeur, Koopman, and Ooms (2011) of the current volume.This is done by incorporating the regression vector β into the state vector and, if state and observation noise sequences are correlated, by incorporating the observation noise into the state as well.Table 4 provides additional clarification of the difference in notation between the two representations.
The SAS/IML documentation contains detailed information on KALDFF and KALDFS subroutines, including examples of their use.These subroutines are organized such that fairly complex state space models can be easily specified and useful quantities needed in their analysis are obtained as output.A brief description of these routines is given below.
lead is the number of steps to forecast after the end of the data.
int contains information about the regressors X t and W t (t ≥ 1).
coef contains information about the system matrices F t and H t (t ≥ 1).
var contains information about the system matrices R t , V t and G t (t ≥ 1).
intd contains information about the system vectors a and b.
coefd contains information about the system matrices A and B. Output: pred contains the one-step-ahead state forecasts.
vpred contains the covariances of the one-step-ahead state forecasts.
initial contains the estimate of the initial state δ along with its covariance estimate.
Input/Output (these arguments act as input as well as output and are better explained in the documentation): -n0 is an optional integer argument that contains divisor used for computing the estimate of σ 2 .
qt is an optional argument that on output contains quantities needed for the likelihood computation.
at is an optional argument that on output contains certain intermediate quantities needed during the subsequent smoothing computations (KALDFS call).As input it is used to provide certain initialization information.
mt is an optional argument that on output contains certain intermediate quantities needed during the subsequent smoothing computations (KALDFS call) and for the likelihood computation.As input it is used to provide certain initialization information (such as V 0 ).
The KALDFS subroutine computes the smoothed state vector and its mean square error matrix from the one-step-ahead forecast and mean square error matrix computed by KALDFF.It has the following signature: CALL KALDFS( sm,vsm,data,int,coef,var,bvec,bmat,initial,at,mt,s2 ,un,vun ) ; where the arguments have the following meaning: Input: data contains the y-values.
int contains information about the regressors X t and W t .
coef contains information about the system matrices F t and H t .
var contains information about the system matrices R t , V t and G t .
bvec contains information about the system vector b.
bmat contains information about the system matrix B.
quantities computed by the preceding KALDFF call: * initial contains the estimate of the initial state δ along with its covariance estimate.* s2 contains estimate of the scalar σ 2 .* at.* mt.

Output:
sm contains the smoothed state vectors.
vsm contains the covariances of the smoothed state vectors.

Input/Output:
un is an optional argument that on output contains certain quantities used in smoothing recursions.
vun is an optional argument that on output contains certain quantities used in smoothing recursions.
The attached program illustrates how to use SAS/IML to specify and fit a simple multivariate state space model.More general state space models, possibly with time-varying system matrices, can be handled in a similar fashion.The remainder of this section describes the model used in this illustration, a broad outline of the attached program, and the produced output.
Consider the following two-way random effects panel model: for t ≥ 1, where y t are k-dimensional observations, X t are k × p dimensional matrices that contain predictor variables, β is a p-dimensional regression parameter vector, the cross-sectional effects vector γ is a k-dimensional zero-mean Gaussian random vector with covariance σ 2 γ I k , τ t (which denote the time effects that are common across cross-sections) is a sequence of uncorrelated, zero-mean Gaussian random variables with variance σ 2 τ , and t is a sequence of uncorrelated, zero-mean k-dimensional random vectors with covariance σ 2 I k .Here I k denotes a k-dimensional identity matrix, 1 denotes a vector with all entries equal to 1, and the random sequences t and τ t , and the random vector γ are assumed to be mutually independent.
The preceding model can be written as a state space model as y t = X t β + Hz t + t z t+1 = Fz t + η t where the (k + 1)-dimensional state vectors z t consist of τ t as the first element and γ as the remaining elements, the k × (k + 1) matrix H is equal to [1 I k ], the (k + 1) × (k + 1) matrix F is a diagonal matrix with all diagonal elements equal to one except the first (which is zero).For t ≥ 1, the (k + 1)-dimensional disturbance vectors η t in the state equation consist of all zeros except the first element, which is τ t .η 0 is a zero mean Gaussian vector with diagonal covariance consisting of σ 2 τ as its first diagonal element and σ 2 γ as its remaining diagonal elements.The initial state, z 0 , is zero-that is, a and A are zero-and the diffuse vector δ is the same as the regression parameter vector β, which means b is zero and B is identity.The observation and state noise sequences are uncorrelated, that is, G t are zero.The tuning parameters in this model are the disturbance variances σ 2 γ , σ 2 τ , and σ 2 , which can be estimated by optimizing the likelihood of the data.The regression parameter vector β, the cross-sectional effects γ, and the time effects τ t are estimated by using state smoothing.
The attached program uses the airline cost data set from Greene (2000) to illustrate how to use SAS/IML to specify and fit this state space model.In the illustration the panel model studies the dependence of the log transformed cost, lC, on three variables-the log-transformed quantity (lQ), the log-transformed price of fuel (lPF), and the load factor (LF).An overall intercept is also included in the model.The data set contains yearly data on six airlines, which form the panels, over the fifteen-year time span 1970-1984.Note that, for this example, the elements of γ signify airline specific corrections to the overall mean and a time effect τ t represents a similar correction to the overall mean at time t.

The program can be broadly described by the following steps:
A subroutine, rantwo, is defined that does filtering, smoothing, and likelihood computation for the two-way random effects panel model.This is done by creating appropriate system matrices corresponding to this model and calling the KALDFF and KALDFS subroutines to do the actual computations.More precisely, the signature of rantwo is as follows: CALL RANTWO( y, X, va, vb, ve, smooth, pred, vpred, initial, dff_logl, sm, vsm); where the arguments have the following meaning: -Input * y contains the y-values, for example, lC values for the six airlines in the case of this example * X contains the regressors, for example, the intercept column and the lQ, lPF, and LF values for the six airlines in the case of this example * va, vb, and ve contain the variances σ 2 γ , σ 2 τ , and σ 2 * smooth, a zero or one flag for deciding whether the smoothing step is to be performed or not -Output * pred and vpred contain the one-step-ahead state predictions and their covariances * sm and vsm contain the smoothed states and their covariances * initial contains the estimate of β and its covariance * dff_logl contains the diffuse log-likelihood A very simple function, rantwo_2ll, is defined, which computes minus two times the log-likelihood of the data for given values of σ 2 γ , σ 2 τ , and σ 2 , by calling rantwo.It assumes that the data, y and X, remain fixed (global variables).
Obtain the maximum likelihood estimates of the variance parameters by minimizing rantwo_2ll with respect to σ 2 γ , σ 2 τ , and σ 2 .This is accomplished by using a SAS/IML non-linear optimization subroutine NLPQN.SAS/IML contains a set of optimization subroutines for minimizing or maximizing a continuous nonlinear function.The parameters can be subject to boundary constraints and linear or nonlinear equality and inequality constraints.
Obtain the estimates of the regression parameter vector β, the cross-sectional effects γ, and the time effects τ t by state smoothing.This is done by calling rantwo at the optimized values of σ 2 γ , σ 2 τ , and σ 2 .
Finally, a portion of the output from this program is shown next.This output has been slightly modified for readability purposes.Table 5 shows the maximum likelihood estimates of σ 2 γ , σ 2 τ , and σ 2 .Similarly, Table 6 shows the estimates of regression effects and their standard errors, and Table 7 shows the smoothed estimate of the panel effect vector γ.

Conclusion
This article and the attached SAS program provide a brief overview of the state space modeling functionality available in SAS/IML and SAS/ETS.Generally speaking, if a state space model belongs to one of the standard forms such as ARIMA, VARMAX, or UCM, it is easier to use the SAS/ETS procedures that are specially designed for their handling.Otherwise, you can use SAS/IML software for general state space modeling.For additional information you can consult SAS online documentation available at http://support.sas.com/documentation/.There you can find SAS Institute Inc. (2010) for SAS/ETS documentation, and SAS Institute Inc. (2008b) for SAS/IML documentation.

Figure 3 :
Figure 3: Local level model: Smoothed and standardized t and ξ t .

Figure 4 :
Figure 4: Local level model: Nile water level forecasts.

Finally, Figure 6 ,
obtained by using the PLOT = PANEL option in the ESTIMATE statement, shows a panel of residual diagnostic plots that are useful for checking the normality and the lack of autocorrelation in the residuals.Figure6does not show any major violations of the model assumptions.

Table 1 :
Parameter estimates for the local level model.

Table 2 :
Table2shows a summary of the most likely outliers in the series based on the fitted model.Additive outliers are reported by default; in addition, level shifts are also reported if the CHECKBREAK option is used in the LEVEL statement.The outlier summary shows that the level shift in the year 1899 is statistically Outlier summary.

Table 3 :
Parameter estimates for the local level plus LS1899 model.

Table 5 :
ML estimates of the variance parameters.

Table 7 :
Smoothed estimate of the panel effects vector γ.