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: Wed, 28 Oct 2026 AIC -696.078
Time: 15:24:38 BIC -665.951
Sample: 04-01-1960 HQIC -684.049
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.04, 10.09 Jarque-Bera (JB): 11.24, 2.44
Prob(Q): 0.83, 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.2408 0.093 -2.592 0.010 -0.423 -0.059
L1.dln_inc 0.2884 0.449 0.642 0.521 -0.592 1.169
L2.dln_inv -0.1642 0.155 -1.059 0.290 -0.468 0.140
L2.dln_inc 0.0810 0.422 0.192 0.848 -0.746 0.908
beta.dln_consump 0.9523 0.639 1.490 0.136 -0.301 2.205
Results for equation dln_inc
====================================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------------
L1.dln_inv 0.0631 0.036 1.761 0.078 -0.007 0.133
L1.dln_inc 0.0816 0.107 0.761 0.447 -0.129 0.292
L2.dln_inv 0.0106 0.033 0.321 0.748 -0.054 0.075
L2.dln_inc 0.0342 0.134 0.255 0.799 -0.229 0.297
beta.dln_consump 0.7734 0.112 6.886 0.000 0.553 0.994
Error covariance matrix
============================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv 0.0434 0.004 12.294 0.000 0.036 0.050
sqrt.cov.dln_inv.dln_inc 5.494e-05 0.002 0.027 0.978 -0.004 0.004
sqrt.var.dln_inc 0.0109 0.001 11.215 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.887
+ intercept AIC -683.774
Date: Wed, 28 Oct 2026 BIC -655.964
Time: 15:26:37 HQIC -672.670
Sample: 04-01-1960
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.00, 0.07 Jarque-Bera (JB): 12.73, 13.42
Prob(Q): 0.95, 0.79 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.01, 4.83
Results for equation dln_inv
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0182 0.005 3.794 0.000 0.009 0.028
L1.e(dln_inv) -0.2583 0.106 -2.441 0.015 -0.466 -0.051
L1.e(dln_inc) 0.4999 0.631 0.792 0.428 -0.737 1.737
L2.e(dln_inv) 0.0313 0.149 0.210 0.834 -0.261 0.324
L2.e(dln_inc) 0.1927 0.476 0.405 0.686 -0.741 1.126
Results for equation dln_inc
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0207 0.002 13.138 0.000 0.018 0.024
L1.e(dln_inv) 0.0476 0.042 1.145 0.252 -0.034 0.129
L1.e(dln_inc) -0.0777 0.140 -0.557 0.578 -0.351 0.196
L2.e(dln_inv) 0.0184 0.042 0.433 0.665 -0.065 0.102
L2.e(dln_inc) 0.1229 0.153 0.804 0.421 -0.177 0.422
Error covariance matrix
==================================================================================
coef std err z P>|z| [0.025 0.975]
----------------------------------------------------------------------------------
sigma2.dln_inv 0.0020 0.000 7.354 0.000 0.001 0.003
sigma2.dln_inc 0.0001 2.33e-05 5.829 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_4461/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.294
+ intercept AIC -682.588
Date: Wed, 28 Oct 2026 BIC -652.461
Time: 15:33:27 HQIC -670.559
Sample: 04-01-1960
- 10-01-1978
Covariance Type: opg
===================================================================================
Ljung-Box (L1) (Q): 0.00, 0.07 Jarque-Bera (JB): 11.16, 14.06
Prob(Q): 0.98, 0.79 Prob(JB): 0.00, 0.00
Heteroskedasticity (H): 0.43, 0.91 Skew: 0.02, -0.45
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.0095 0.060 0.157 0.875 -0.109 0.128
L1.dln_inv -0.0088 0.646 -0.014 0.989 -1.275 1.258
L1.dln_inc 0.4316 2.584 0.167 0.867 -4.633 5.497
L1.e(dln_inv) -0.2439 0.656 -0.372 0.710 -1.529 1.041
L1.e(dln_inc) 0.1143 2.842 0.040 0.968 -5.456 5.685
Results for equation dln_inc
=================================================================================
coef std err z P>|z| [0.025 0.975]
---------------------------------------------------------------------------------
intercept 0.0163 0.025 0.653 0.514 -0.033 0.065
L1.dln_inv -0.0299 0.258 -0.116 0.908 -0.536 0.476
L1.dln_inc 0.2424 1.017 0.238 0.812 -1.750 2.235
L1.e(dln_inv) 0.0865 0.265 0.327 0.744 -0.432 0.605
L1.e(dln_inc) -0.2406 1.055 -0.228 0.820 -2.307 1.826
Error covariance matrix
============================================================================================
coef std err z P>|z| [0.025 0.975]
--------------------------------------------------------------------------------------------
sqrt.var.dln_inv 0.0448 0.003 14.537 0.000 0.039 0.051
sqrt.cov.dln_inv.dln_inc 0.0017 0.003 0.660 0.510 -0.003 0.007
sqrt.var.dln_inc 0.0116 0.001 11.821 0.000 0.010 0.013
============================================================================================
Warnings:
[1] Covariance matrix calculated using the outer product of gradients (complex-step).