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.

[1]:
%matplotlib inline
[2]:
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
[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']]
---------------------------------------------------------------------------
ConnectionRefusedError                    Traceback (most recent call last)
/usr/lib/python3.8/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1349             try:
-> 1350                 h.request(req.get_method(), req.selector, req.data, headers,
   1351                           encode_chunked=req.has_header('Transfer-encoding'))

/usr/lib/python3.8/http/client.py in request(self, method, url, body, headers, encode_chunked)
   1254         """Send a complete request to the server."""
-> 1255         self._send_request(method, url, body, headers, encode_chunked)
   1256

/usr/lib/python3.8/http/client.py in _send_request(self, method, url, body, headers, encode_chunked)
   1300             body = _encode(body, 'body')
-> 1301         self.endheaders(body, encode_chunked=encode_chunked)
   1302

/usr/lib/python3.8/http/client.py in endheaders(self, message_body, encode_chunked)
   1249             raise CannotSendHeader()
-> 1250         self._send_output(message_body, encode_chunked=encode_chunked)
   1251

/usr/lib/python3.8/http/client.py in _send_output(self, message_body, encode_chunked)
   1009         del self._buffer[:]
-> 1010         self.send(msg)
   1011

/usr/lib/python3.8/http/client.py in send(self, data)
    949             if self.auto_open:
--> 950                 self.connect()
    951             else:

/usr/lib/python3.8/http/client.py in connect(self)
   1416
-> 1417             super().connect()
   1418

/usr/lib/python3.8/http/client.py in connect(self)
    920         """Connect to the host and port specified in __init__."""
--> 921         self.sock = self._create_connection(
    922             (self.host,self.port), self.timeout, self.source_address)

/usr/lib/python3.8/socket.py in create_connection(address, timeout, source_address)
    807         try:
--> 808             raise err
    809         finally:

/usr/lib/python3.8/socket.py in create_connection(address, timeout, source_address)
    795                 sock.bind(source_address)
--> 796             sock.connect(sa)
    797             # Break explicitly a reference cycle

ConnectionRefusedError: [Errno 111] Connection refused

During handling of the above exception, another exception occurred:

URLError                                  Traceback (most recent call last)
<ipython-input-3-a26cf3e09ec4> in <module>
----> 1 dta = sm.datasets.webuse('lutkepohl2', 'https://www.stata-press.com/data/r12/')
      2 dta.index = dta.qtr
      3 endog = dta.loc['1960-04-01':'1978-10-01', ['dln_inv', 'dln_inc', 'dln_consump']]

/usr/lib/python3/dist-packages/statsmodels/datasets/utils.py in webuse(data, baseurl, as_df)
     41     """
     42     url = urljoin(baseurl, data+'.dta')
---> 43     return read_stata(url)
     44
     45

/usr/lib/python3/dist-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

/usr/lib/python3/dist-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

/usr/lib/python3/dist-packages/pandas/io/stata.py in read_stata(filepath_or_buffer, convert_dates, convert_categoricals, encoding, index_col, convert_missing, preserve_dtypes, columns, order_categoricals, chunksize, iterator)
    219 ):
    220
--> 221     reader = StataReader(
    222         filepath_or_buffer,
    223         convert_dates=convert_dates,

/usr/lib/python3/dist-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

/usr/lib/python3/dist-packages/pandas/util/_decorators.py in wrapper(*args, **kwargs)
    206                 else:
    207                     kwargs[new_arg_name] = new_arg_value
--> 208             return func(*args, **kwargs)
    209
    210         return wrapper

/usr/lib/python3/dist-packages/pandas/io/stata.py in __init__(self, path_or_buf, convert_dates, convert_categoricals, index_col, convert_missing, preserve_dtypes, columns, order_categoricals, encoding, chunksize)
   1093         path_or_buf = _stringify_path(path_or_buf)
   1094         if isinstance(path_or_buf, str):
-> 1095             path_or_buf, encoding, _, should_close = get_filepath_or_buffer(path_or_buf)
   1096
   1097         if isinstance(path_or_buf, (str, bytes)):

/usr/lib/python3/dist-packages/pandas/io/common.py in get_filepath_or_buffer(filepath_or_buffer, encoding, compression, mode)
    194
    195     if _is_url(filepath_or_buffer):
--> 196         req = urlopen(filepath_or_buffer)
    197         content_encoding = req.headers.get("Content-Encoding", None)
    198         if content_encoding == "gzip":

/usr/lib/python3.8/urllib/request.py in urlopen(url, data, timeout, cafile, capath, cadefault, context)
    220     else:
    221         opener = _opener
--> 222     return opener.open(url, data, timeout)
    223
    224 def install_opener(opener):

/usr/lib/python3.8/urllib/request.py in open(self, fullurl, data, timeout)
    523
    524         sys.audit('urllib.Request', req.full_url, req.data, req.headers, req.get_method())
--> 525         response = self._open(req, data)
    526
    527         # post-process response

/usr/lib/python3.8/urllib/request.py in _open(self, req, data)
    540
    541         protocol = req.type
--> 542         result = self._call_chain(self.handle_open, protocol, protocol +
    543                                   '_open', req)
    544         if result:

/usr/lib/python3.8/urllib/request.py in _call_chain(self, chain, kind, meth_name, *args)
    500         for handler in handlers:
    501             func = getattr(handler, meth_name)
--> 502             result = func(*args)
    503             if result is not None:
    504                 return result

/usr/lib/python3.8/urllib/request.py in https_open(self, req)
   1391
   1392         def https_open(self, req):
-> 1393             return self.do_open(http.client.HTTPSConnection, req,
   1394                 context=self._context, check_hostname=self._check_hostname)
   1395

/usr/lib/python3.8/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
   1351                           encode_chunked=req.has_header('Transfer-encoding'))
   1352             except OSError as err: # timeout error
-> 1353                 raise URLError(err)
   1354             r = h.getresponse()
   1355         except:

URLError: <urlopen error [Errno 111] Connection refused>

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())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-4-75f089c68993> in <module>
----> 1 exog = endog['dln_consump']
      2 mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(2,0), trend='n', exog=exog)
      3 res = mod.fit(maxiter=1000, disp=False)
      4 print(res.summary())

NameError: name 'endog' is not defined

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

[5]:
ax = res.impulse_responses(10, orthogonalized=True).plot(figsize=(13,3))
ax.set(xlabel='t', title='Responses to a shock to `dln_inv`');
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-5-5a59fc5576cf> in <module>
----> 1 ax = res.impulse_responses(10, orthogonalized=True).plot(figsize=(13,3))
      2 ax.set(xlabel='t', title='Responses to a shock to `dln_inv`');

NameError: name 'res' is not defined

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())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-6-fc41e9fb664e> in <module>
----> 1 mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(0,2), error_cov_type='diagonal')
      2 res = mod.fit(maxiter=1000, disp=False)
      3 print(res.summary())

NameError: name 'endog' is not defined

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())
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-7-38abad3e3464> in <module>
----> 1 mod = sm.tsa.VARMAX(endog[['dln_inv', 'dln_inc']], order=(1,1))
      2 res = mod.fit(maxiter=1000, disp=False)
      3 print(res.summary())

NameError: name 'endog' is not defined