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 \(\mathrm{k_endog} \times 1\) vector.
[1]:
%matplotlib inline
[2]:
import pandas as pd
import statsmodels.api as sm
[3]:
import shutil
import requests
def download_file(url):
local_filename = url.split("/")[-1]
with requests.get(url, stream=True, timeout=30) as r:
with open(local_filename, "wb") as f:
shutil.copyfileobj(r.raw, f)
return local_filename
filename = download_file("https://www.stata-press.com/data/r12/lutkepohl2.dta")
dta = pd.read_stata(filename)
dta.index = dta.qtr
dta.index.freq = "QS"
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.
[4]:
exog = endog["dln_consump"]
mod = sm.tsa.VARMAX(endog[["dln_inv", "dln_inc"]], order=(2, 0), trend="n", exog=exog)
res = mod.fit(maxiter=1000, disp=False)
print(res.summary())
Statespace Model Results
==================================================================================
Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75
Model: VARX(2) Log Likelihood 361.039
Date: Sun, 13 Sep 2026 AIC -696.078
Time: 00:03:20 BIC -665.950
Sample: 04-01-1960 HQIC -684.048
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.05, 10.18 Jarque-Bera (JB): 11.22, 2.41
Prob(Q): 0.82, 0.00 Prob(JB): 0.00, 0.30
Heteroskedasticity (H): 0.45, 0.40 Skew: 0.16, -0.38
Prob(H) (two-sided): 0.05, 0.03 Kurtosis: 4.87, 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.2905 0.449 0.646 0.518 -0.590 1.171
L2.dln_inv -0.1660 0.155 -1.070 0.285 -0.470 0.138
L2.dln_inc 0.0685 0.421 0.163 0.871 -0.757 0.894
beta.dln_consump 0.9622 0.639 1.506 0.132 -0.290 2.214
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.0820 0.107 0.764 0.445 -0.128 0.292
L2.dln_inv 0.0102 0.033 0.308 0.758 -0.054 0.075
L2.dln_inc 0.0310 0.134 0.231 0.817 -0.232 0.294
beta.dln_consump 0.7768 0.112 6.921 0.000 0.557 0.997
Error covariance matrix
============================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv 0.0434 0.004 12.300 0.000 0.036 0.050
sqrt.cov.dln_inv.dln_inc 5.183e-05 0.002 0.026 0.979 -0.004 0.004
sqrt.var.dln_inc 0.0109 0.001 11.214 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.
[5]:
ax = res.impulse_responses(10, orthogonalized=True, impulse=[1, 0]).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.
[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())
Statespace Model Results
==================================================================================
Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75
Model: VMA(2) Log Likelihood 353.886
+ intercept AIC -683.773
Date: Sun, 13 Sep 2026 BIC -655.963
Time: 00:09:36 HQIC -672.669
Sample: 04-01-1960
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.00, 0.06 Jarque-Bera (JB): 12.84, 13.32
Prob(Q): 0.97, 0.80 Prob(JB): 0.00, 0.00
Heteroskedasticity (H): 0.44, 0.81 Skew: 0.06, -0.48
Prob(H) (two-sided): 0.04, 0.60 Kurtosis: 5.02, 4.83
Results for equation dln_inv
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0182 0.005 3.810 0.000 0.009 0.028
L1.e(dln_inv) -0.2543 0.106 -2.401 0.016 -0.462 -0.047
L1.e(dln_inc) 0.5239 0.633 0.828 0.408 -0.716 1.764
L2.e(dln_inv) 0.0258 0.150 0.172 0.863 -0.268 0.319
L2.e(dln_inc) 0.1684 0.474 0.355 0.723 -0.762 1.098
Results for equation dln_inc
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0208 0.002 13.084 0.000 0.018 0.024
L1.e(dln_inv) 0.0478 0.042 1.152 0.249 -0.034 0.129
L1.e(dln_inc) -0.0780 0.139 -0.562 0.574 -0.350 0.194
L2.e(dln_inv) 0.0183 0.042 0.431 0.666 -0.065 0.101
L2.e(dln_inc) 0.1284 0.152 0.842 0.400 -0.170 0.427
Error covariance matrix
==================================================================================
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------
sigma2.dln_inv 0.0020 0.000 7.341 0.000 0.001 0.003
sigma2.dln_inc 0.0001 2.33e-05 5.826 0.000 9e-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.
[7]:
mod = sm.tsa.VARMAX(endog[["dln_inv", "dln_inc"]], order=(1, 1))
res = mod.fit(maxiter=1000, disp=False)
print(res.summary())
/tmp/ipykernel_4750/242493548.py:1: EstimationWarning: Estimation of VARMA(p,q) models is not generically robust, due especially to identification issues.
mod = sm.tsa.VARMAX(endog[["dln_inv", "dln_inc"]], order=(1, 1))
Statespace Model Results
==================================================================================
Dep. Variable: ['dln_inv', 'dln_inc'] No. Observations: 75
Model: VARMA(1,1) Log Likelihood 354.288
+ intercept AIC -682.575
Date: Sun, 13 Sep 2026 BIC -652.448
Time: 00:14:16 HQIC -670.546
Sample: 04-01-1960
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.00, 0.06 Jarque-Bera (JB): 11.11, 14.14
Prob(Q): 0.95, 0.81 Prob(JB): 0.00, 0.00
Heteroskedasticity (H): 0.43, 0.91 Skew: 0.01, -0.46
Prob(H) (two-sided): 0.04, 0.81 Kurtosis: 4.89, 4.92
Results for equation dln_inv
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0105 0.065 0.162 0.871 -0.117 0.138
L1.dln_inv -0.0063 0.697 -0.009 0.993 -1.372 1.359
L1.dln_inc 0.3795 2.739 0.139 0.890 -4.988 5.747
L1.e(dln_inv) -0.2475 0.707 -0.350 0.726 -1.634 1.139
L1.e(dln_inc) 0.1251 2.988 0.042 0.967 -5.732 5.982
Results for equation dln_inc
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0165 0.027 0.607 0.544 -0.037 0.070
L1.dln_inv -0.0333 0.278 -0.120 0.905 -0.579 0.512
L1.dln_inc 0.2347 1.104 0.212 0.832 -1.930 2.399
L1.e(dln_inv) 0.0890 0.285 0.312 0.755 -0.469 0.647
L1.e(dln_inc) -0.2365 1.139 -0.208 0.836 -2.469 1.996
Error covariance matrix
============================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv 0.0449 0.003 14.523 0.000 0.039 0.051
sqrt.cov.dln_inv.dln_inc 0.0017 0.003 0.651 0.515 -0.003 0.007
sqrt.var.dln_inc 0.0116 0.001 11.724 0.000 0.010 0.013
============================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).