VARMAX models

This is a brief introduction notebook to VARMAX models in Statsmodels. The VARMAX model is generically specified as: $$ y_t = \nu + A_1 y_{t-1} + \dots + A_p y_{t-p} + B x_t + \epsilon_t + M_1 \epsilon_{t-1} + \dots M_q \epsilon_{t-q} $$

where $y_t$ is a $\text{k_endog} \times 1$ vector.

In [1]:
%matplotlib inline
In [2]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
In [3]:
dta = sm.datasets.webuse('lutkepohl2', 'https://www.stata-press.com/data/r12/')
dta.index = dta.qtr
endog = dta.loc['1960-04-01':'1978-10-01', ['dln_inv', 'dln_inc', 'dln_consump']]

Model specification

The VARMAX class in Statsmodels allows estimation of VAR, VMA, and VARMA models (through the order argument), optionally with a constant term (via the trend argument). Exogenous regressors may also be included (as usual in Statsmodels, by the exog argument), and in this way a time trend may be added. Finally, the class allows measurement error (via the measurement_error argument) and allows specifying either a diagonal or unstructured innovation covariance matrix (via the error_cov_type argument).

Example 1: VAR

Below is a simple VARX(2) model in two endogenous variables and an exogenous series, but no constant term. Notice that we needed to allow for more iterations than the default (which is maxiter=50) in order for the likelihood estimation to converge. This is not unusual in VAR models which have to estimate a large number of parameters, often on a relatively small number of time series: this model, for example, estimates 27 parameters off of 75 observations of 3 variables.

In [4]:
exog = endog['dln_consump']
mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(2,0), trend='nc', exog=exog)
res = mod.fit(maxiter=1000, disp=False)
print(res.summary())
/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  % freq, ValueWarning)
                             Statespace Model Results                             
==================================================================================
Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75
Model:                            VARX(2)   Log Likelihood                 361.037
Date:                    Sun, 04 Nov 2018   AIC                           -696.075
Time:                            21:37:09   BIC                           -665.947
Sample:                        04-01-1960   HQIC                          -684.045
                             - 10-01-1978                                         
Covariance Type:                      opg                                         
===================================================================================
Ljung-Box (Q):                61.23, 39.38   Jarque-Bera (JB):          11.01, 2.45
Prob(Q):                        0.02, 0.50   Prob(JB):                   0.00, 0.29
Heteroskedasticity (H):         0.45, 0.40   Skew:                      0.16, -0.38
Prob(H) (two-sided):            0.05, 0.03   Kurtosis:                   4.85, 3.44
                            Results for equation dln_inv                            
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
L1.dln_inv          -0.2391      0.093     -2.571      0.010      -0.421      -0.057
L1.dln_inc           0.2883      0.449      0.643      0.520      -0.591       1.167
L2.dln_inv          -0.1676      0.155     -1.081      0.280      -0.471       0.136
L2.dln_inc           0.0547      0.420      0.130      0.896      -0.768       0.878
beta.dln_consump     0.9853      0.636      1.549      0.121      -0.262       2.232
                            Results for equation dln_inc                            
====================================================================================
                       coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------------
L1.dln_inv           0.0632      0.036      1.766      0.077      -0.007       0.133
L1.dln_inc           0.0817      0.107      0.762      0.446      -0.128       0.292
L2.dln_inv           0.0109      0.033      0.330      0.741      -0.054       0.076
L2.dln_inc           0.0315      0.134      0.235      0.814      -0.232       0.295
beta.dln_consump     0.7754      0.112      6.894      0.000       0.555       0.996
                                  Error covariance matrix                                   
============================================================================================
                               coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv             0.0433      0.004     12.308      0.000       0.036       0.050
sqrt.cov.dln_inv.dln_inc  5.653e-05      0.002      0.028      0.978      -0.004       0.004
sqrt.var.dln_inc             0.0109      0.001     11.212      0.000       0.009       0.013
============================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

From the estimated VAR model, we can plot the impulse response functions of the endogenous variables.

In [5]:
ax = res.impulse_responses(10, orthogonalized=True).plot(figsize=(13,3))
ax.set(xlabel='t', title='Responses to a shock to `dln_inv`');

Example 2: VMA

A vector moving average model can also be formulated. Below we show a VMA(2) on the same data, but where the innovations to the process are uncorrelated. In this example we leave out the exogenous regressor but now include the constant term.

In [6]:
mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(0,2), error_cov_type='diagonal')
res = mod.fit(maxiter=1000, disp=False)
print(res.summary())
/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  % freq, ValueWarning)
                             Statespace Model Results                             
==================================================================================
Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75
Model:                             VMA(2)   Log Likelihood                 353.888
                              + intercept   AIC                           -683.775
Date:                    Sun, 04 Nov 2018   BIC                           -655.965
Time:                            21:37:15   HQIC                          -672.671
Sample:                        04-01-1960                                         
                             - 10-01-1978                                         
Covariance Type:                      opg                                         
===================================================================================
Ljung-Box (Q):                68.61, 39.16   Jarque-Bera (JB):         12.55, 13.54
Prob(Q):                        0.00, 0.51   Prob(JB):                   0.00, 0.00
Heteroskedasticity (H):         0.44, 0.81   Skew:                      0.05, -0.48
Prob(H) (two-sided):            0.04, 0.60   Kurtosis:                   5.00, 4.84
                           Results for equation dln_inv                          
=================================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------
intercept         0.0182      0.005      3.789      0.000       0.009       0.028
L1.e(dln_inv)    -0.2594      0.106     -2.454      0.014      -0.467      -0.052
L1.e(dln_inc)     0.5207      0.632      0.824      0.410      -0.718       1.759
L2.e(dln_inv)     0.0331      0.149      0.223      0.824      -0.258       0.324
L2.e(dln_inc)     0.1723      0.475      0.363      0.716      -0.758       1.102
                           Results for equation dln_inc                          
=================================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------
intercept         0.0207      0.002     13.019      0.000       0.018       0.024
L1.e(dln_inv)     0.0483      0.042      1.160      0.246      -0.033       0.130
L1.e(dln_inc)    -0.0739      0.140     -0.529      0.597      -0.347       0.200
L2.e(dln_inv)     0.0182      0.042      0.428      0.669      -0.065       0.101
L2.e(dln_inc)     0.1291      0.153      0.843      0.399      -0.171       0.429
                             Error covariance matrix                              
==================================================================================
                     coef    std err          z      P>|z|      [0.025      0.975]
----------------------------------------------------------------------------------
sigma2.dln_inv     0.0020      0.000      7.360      0.000       0.001       0.003
sigma2.dln_inc     0.0001   2.33e-05      5.811      0.000    8.99e-05       0.000
==================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).

Caution: VARMA(p,q) specifications

Although the model allows estimating VARMA(p,q) specifications, these models are not identified without additional restrictions on the representation matrices, which are not built-in. For this reason, it is recommended that the user proceed with error (and indeed a warning is issued when these models are specified). Nonetheless, they may in some circumstances provide useful information.

In [7]:
mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(1,1))
res = mod.fit(maxiter=1000, disp=False)
print(res.summary())
/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/statespace/varmax.py:159: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.
  EstimationWarning)
/home/travis/build/statsmodels/statsmodels/statsmodels/tsa/base/tsa_model.py:165: ValueWarning: No frequency information was provided, so inferred frequency QS-OCT will be used.
  % freq, ValueWarning)
                             Statespace Model Results                             
==================================================================================
Dep. Variable:     ['dln_inv', 'dln_inc']   No. Observations:                   75
Model:                         VARMA(1,1)   Log Likelihood                 354.281
                              + intercept   AIC                           -682.562
Date:                    Sun, 04 Nov 2018   BIC                           -652.435
Time:                            21:37:17   HQIC                          -670.533
Sample:                        04-01-1960                                         
                             - 10-01-1978                                         
Covariance Type:                      opg                                         
===================================================================================
Ljung-Box (Q):                68.80, 38.96   Jarque-Bera (JB):         10.76, 14.04
Prob(Q):                        0.00, 0.52   Prob(JB):                   0.00, 0.00
Heteroskedasticity (H):         0.43, 0.91   Skew:                     -0.00, -0.46
Prob(H) (two-sided):            0.04, 0.81   Kurtosis:                   4.86, 4.91
                           Results for equation dln_inv                          
=================================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------
intercept         0.0112      0.069      0.161      0.872      -0.124       0.147
L1.dln_inv       -0.0101      0.735     -0.014      0.989      -1.450       1.430
L1.dln_inc        0.3508      2.889      0.121      0.903      -5.311       6.012
L1.e(dln_inv)    -0.2499      0.747     -0.335      0.738      -1.714       1.214
L1.e(dln_inc)     0.1262      3.129      0.040      0.968      -6.006       6.258
                           Results for equation dln_inc                          
=================================================================================
                    coef    std err          z      P>|z|      [0.025      0.975]
---------------------------------------------------------------------------------
intercept         0.0166      0.029      0.570      0.569      -0.040       0.074
L1.dln_inv       -0.0329      0.292     -0.113      0.910      -0.605       0.539
L1.dln_inc        0.2294      1.180      0.194      0.846      -2.084       2.542
L1.e(dln_inv)     0.0884      0.299      0.296      0.767      -0.497       0.674
L1.e(dln_inc)    -0.2342      1.211     -0.193      0.847      -2.609       2.140
                                  Error covariance matrix                                   
============================================================================================
                               coef    std err          z      P>|z|      [0.025      0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv             0.0450      0.003     14.522      0.000       0.039       0.051
sqrt.cov.dln_inv.dln_inc     0.0017      0.003      0.646      0.518      -0.003       0.007
sqrt.var.dln_inc             0.0115      0.001     11.648      0.000       0.010       0.013
============================================================================================

Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).