VARMAX models¶
This is a brief introduction notebook to VARMAX models in statsmodels. The VARMAX model is generically specified as:
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)
1318 try:
-> 1319 h.request(req.get_method(), req.selector, req.data, headers,
1320 encode_chunked=req.has_header('Transfer-encoding'))
/usr/lib/python3.8/http/client.py in request(self, method, url, body, headers, encode_chunked)
1229 """Send a complete request to the server."""
-> 1230 self._send_request(method, url, body, headers, encode_chunked)
1231
/usr/lib/python3.8/http/client.py in _send_request(self, method, url, body, headers, encode_chunked)
1275 body = _encode(body, 'body')
-> 1276 self.endheaders(body, encode_chunked=encode_chunked)
1277
/usr/lib/python3.8/http/client.py in endheaders(self, message_body, encode_chunked)
1224 raise CannotSendHeader()
-> 1225 self._send_output(message_body, encode_chunked=encode_chunked)
1226
/usr/lib/python3.8/http/client.py in _send_output(self, message_body, encode_chunked)
1003 del self._buffer[:]
-> 1004 self.send(msg)
1005
/usr/lib/python3.8/http/client.py in send(self, data)
943 if self.auto_open:
--> 944 self.connect()
945 else:
/usr/lib/python3.8/http/client.py in connect(self)
1391
-> 1392 super().connect()
1393
/usr/lib/python3.8/http/client.py in connect(self)
914 """Connect to the host and port specified in __init__."""
--> 915 self.sock = self._create_connection(
916 (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)
216 ):
217
--> 218 reader = StataReader(
219 filepath_or_buffer,
220 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)
1088 path_or_buf = _stringify_path(path_or_buf)
1089 if isinstance(path_or_buf, str):
-> 1090 path_or_buf, encoding, _, should_close = get_filepath_or_buffer(path_or_buf)
1091
1092 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)
1360
1361 def https_open(self, req):
-> 1362 return self.do_open(http.client.HTTPSConnection, req,
1363 context=self._context, check_hostname=self._check_hostname)
1364
/usr/lib/python3.8/urllib/request.py in do_open(self, http_class, req, **http_conn_args)
1320 encode_chunked=req.has_header('Transfer-encoding'))
1321 except OSError as err: # timeout error
-> 1322 raise URLError(err)
1323 r = h.getresponse()
1324 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