Adventures with State Space Models 4: Multiple Dependent Variables

[#1]

Welcome to the fourth installment of our State Space model (SSM) post series. Previous posts in this series (links at the end) introduced SSMs as a collection of independent, additive components and detailed differences between the two component types: dynamic and static. Post #3 built on the first two, and the focus was on fitting a SSM with a dynamic input variable component. Variances or hyperparameters in the SSM enable the model’s dynamic components. This post will focus more on the specification of model’s hyperparameters than previous ones and presents a multivariate or multiple dependent variable example. The purpose of this post is to show that the ability to specify model variance parameters in a customizable and flexible way is another distinct advantage of SSMs.

Two model variances we’re introduced in the basic, single dependent variable SSM, shown below.

Select any image to see a larger version.
Mobile users: To view the images, select the "Full" version at the bottom of the page.

Sigma squared is the variance on EPSILON. This is the white noise error term in the observation equation. Delta squared, the variance on the state, regulates how the dynamic part of the model evolves over time. In previous posts, you may have noticed that these hyperparameters are a main feature of estimation results. This focus on variances is somewhat different than other model specifications for timeseries, but it is another feature of SSMs that make them particularly useful. To illustrate this, we’ll move from a model with a single dependent variable to a system of dependent variables and show that by specifying the variances in the model, the SSM easily generalizes to a multivariate framework.

Multivariate systems of equations have more than one dependent variable, and these variables move together or are correlated. A feedback loop is the most straight-forward way to describe a multivariate relationship; current or contemporaneous variation in one dependent variable impacts current variation in other dependent variables in the model and vice versa. Dynamics of mink and muskrat populations provide a textbook example. Minks eat muskrats, and if the mink population grows above a steady state rate, the muskrat population declines. This depletes the mink’s main food source. The mink population declines as a result which allows the muskrat population to grow more rapidly, and so on.

A plot of mink and muskrat populations is shown below. Both variables have a year interval and have been log transformed. This is a subset of the historical fur-return data collected by the Hudson’s Bay Company.

A model for the evolution of the mink and muskrat populations can be specified as below. This specification is known as a vector auto-regression, lag 1 or VAR1.

The model has two dependent or endogenous variables, Minks and Muskrats at time t. They are specified as a function of dependent variable lags. Lags of dependent variables are considered known or exogenous, so, as written, the multivariate component is not obvious. That is, there’s nothing to indicate that the two equations are not separate and independent.

Epsilons 1 & 2 are white noise error terms, but they are not necessarily independent. The connection or link between the equations comes from the variance-covariance matrix shown below.

The feedback or contemporaneous relationship between the mink and muskrat populations is accommodated in the off-diagonal terms in the matrix. If sigma(1,2) = sigma(2,1) = 0, then the mink and muskrat equations are separate, independent equations. That is, they don’t systematically co-vary contemporaneously. If the cross-covariance terms are not equal to zero, then the equations are a dependent system.

Next, we’ll specify and fit two examples of a VAR1 in the SSM procedure. The first specification reproduces the VAR1 as listed above and shows how the variance-covariance matrix can be specified. The second example adds flexibility by allowing multiple sources of co-variation between the dependent variables.

Demonstration 1

Here’s the syntax for the VAR1 as written above.

Notes on the syntax.

Specifying the bi-variate white noise component in the VAR1.

• The PARMS statements create parameters to be estimated, and that go into the variance-covariance matrix.
• The 2x2 array, wnQ, specifies the variance-covariance matrix. The (1,1) element (first row, first column) represents the variance of the mink observation equation. The (1, 2) and (2,1) elements represent the covariance between the mink and muskrat observation equations. RHO is included to ensure that the estimated matrix will be positive definite. The (2,2) element represents the variance of the muskrat observation equation.
• The two-dimensional state named ERROR is a white noise (WN) type. It has a general type (g) covariance matrix that is equal to wnQ.
• The two COMPONENT statements specify the bi-variate white noise error terms that go into the observation equations. We’ve named the first component MINKWN, and it’s equal to the first element of the two-dimensional state, ERROR.

Specifying the other parts of the VAR1.

• The DEPLAG statement is another convenience provided by the SAS developers. It creates model components that contain lags of dependent variables in the observation equations. Here, the component named LAGSFORMINK is anchored to the LMINK (log of MINK) dependent variable. It creates a linear combination of lag 1 values for LMINK and LMUSK.
• The two MODEL statements specify the two observation equations in the VAR1.

Results

The Estimates of Named Parameters show that all of the elements of the variance-covariance matrix are significantly different from zero. The scaling factor, RHO, is not close to a boundary value.

The Model Parameter Estimates table shows the estimated relationships associated with the dependent lags in the VAR1. The estimates for autoregressive lags in the mink and muskrat equations show significant, positive auto-correlation. This indicates that, for example, when the mink population jumps above average in a given year, it tends to stay above average with decay back to the status-quo in following years. The lag 1 correlation between mink and muskrat is estimated to be 0.270. This makes sense, because when the mink food source increases this year, the mink population will tend to increase in the following year. The lag 1 correlation between muskrat and mink is estimated to be -0.283. This indicates that when the mink population increases this year, the muskrat population declines in the following year.

The Information Criteria provide baseline fit statistics.

The estimated variance-covariance and correlation matrices are shown below. The off-diagonal terms suggest that a multivariate relationship exists between the mink and muskrat populations.

Demonstration 2

Now, we’ll add some flexibility to the VAR1 specification and also generalize the model to better accommodate some other aspects of the data. You may have noticed in the plot of the data, shown above, that both mink and muskrat populations tend to drift upwards over time. Many textbook treatments include a step to first difference both dependent variables before fitting the VAR1 to the data to account for this. Here, we’ll accommodate the trend by adding a random walk with drift trend component to each equation in the model. We’ll also include another source of covariation between the two dependent variables by trying the evolution of two trend components together using a variance-covariance matrix similar to the one used to create bi-variate white noise component above.

Here’s the syntax for the modified VAR1.

Notes on the new syntax.

•  New defined parameters in the PARMS statement are ESD1 and ESD2. These represent the standard deviations to be estimated for the dynamic trend terms. The new scaling parameter, RHO2, will be associated with block of the variance-covariance matrix that corresponds to the bi-variate trend terms.
• Syntax associated with the ARRAY, rwQ creates the 2x2 block of the variance-covariance matrix corresponding to the trend terms in the model.
• The two dimensional STATE element ALPHA is a random walk (rw) type. COV specifies the block of the corresponding variance-covariance matrix as equal to rwQ. In general, the W(.) term specifies regression effects that enter the model through the state equation. Here, we’re specifying it as a 2x2 identify matrix, and it effectively adds an intercept term to each of the random walk trends. This converts the random walks into random walks with drift trends which is consistent with the patterns seen in the plot of the data.
• The two new COMPONENT statements specify the bi-variate trend components that go into the observation equations. We’ve named the first component MINKLEVEL, and it’s equal to the first element of the two-dimensional state, ALPHA.

Note, the syntax for this model was provided by Rajesh Selukar.

Results

The Estimates of the State Equation Regression Vector indicate that the drift terms in the bi-variate trend are not significantly different from zero.

The Estimates of Named Parameters indicates that MSD1, the estimated standard deviation of the trend term in the mink equation, is not significantly different from zero in the presence of other parameters estimated in the model. RHO1, the scaling factor in the variance-covariance matrix for the bi-variate white noise components, is also at a boundary value. In spite of these potential issues, the Likelihood optimization algorithm converged. All other terms associated with the model’s variance-covariance matrix blocks are significant.

The Model Parameter Estimates associated with the dependent lags are similar to the ones described above. The negative correlation between the mink population this year and the muskrat population next year is stronger than in the previous model.

The Information Criteria suggest that adding flexibility to the model via the bi-variate trend improves the penalized, overall fit relative to the previous model.

The block of the estimated variance-covariance matrix corresponding to the bi-variate white noise component is similar to the previous model. The estimated variance associated with muskrats in this block (2,2) has diminished in the presence of other parameters in the model.

The estimated values in the block of the variance-covariance matrix corresponding to the bi-variate trend are relatively small. However, the multivariate or feedback effects in the system are a function of off diagonal elements of both variance-covariance blocks, and the feedback effects are estimated to be stronger in this model than in the baseline model.

It’s interesting that the estimated covariance effects in both models are positive, indicating that the mink and muskrat populations move together contemporaneously. This could be the result of milder or harsher winters leading to better or worse breeding conditions for both species. Both species were also harvested for their fur, and the positive covariance could also be a result of tastes and preferences of fur wearing consumers at the time.

Smoothed estimates of LMINK and LMUSKRAT provide forecasts, shown below. Smoothed and filtered estimates for the dependent variables and other model components are contained in the output table generated in the syntax, MKNMSK2. We’ll describe the differences between smoothed and filtered estimates as well as other SSM details in the next post in this series.

In the next post in this series, we’ll take a step back from applications of SSMs and spend some time exploring how the Kalman-filter-smoother algorithm works. Stay tuned for more SSM action!

SSM Post 1, https://communities.sas.com/t5/SAS-Communities-Library/Adventures-with-State-Space-Models-Introducti...

SSM Post 2, https://communities.sas.com/t5/SAS-Communities-Library/Adventures-with-State-Space-Models-2-More-Dyn...

SSM Post 3, https://communities.sas.com/t5/SAS-Communities-Library/Adventures-with-State-Space-models-3-Dynamic-...

Find more articles from SAS Global Enablement and Learning here.