Saturday, June 11, 2011

Monte Carlo Standard Errors for Impulse Response Functions

I've come to a major checkpoint in integrating Monte Carlo error bands
for impulse response functions (this is only non-orthogonal right now).

Here is some quick code to get VAR IRF MC standard errors:

import numpy as np import scikits.statsmodels.api as sm from scikits.statsmodels.tsa.api import VAR import matplotlib.pyplot as plt mdata = sm.datasets.macrodata.load().data mdata = mdata[['realgdp','realcons','realinv']] #mdata = mdata[['realgdp','realcons','realinv', 'pop', 'realdpi', 'tbilrate']] names = mdata.dtype.names data = mdata.view((float,3)) data = np.diff(np.log(data), axis=0) mymodel = VAR(data,names=names) res =,ic=None) #first generate asymptotic standard errors (to compare) res.irf(periods=30).plot(orth=False, stderr_type='asym') #then generate monte carlo standard errors res.irf(periods=30).plot(orth=False, stderr_type='mc', seed=29, repl=1000)

Produces the following plots of a shock to (delta) realgdp.


Monte Carlo:

I added functions to the VARResults, IRAnalysis, and also added to
some of the pre-existing functions in these classes and also to Because Monte Carlo standard errors are in general not
symmetric, I had to alter the plot_with_error function from the file and number of other functions.

Functions added:

VARResults.stderr_MC_irf(self, orth=False, repl=1000, T=25,
signif=0.05, seed=None)
This function generates a tuple that holds the lower and upper error
bands generated by Monte Carlo simulations.

IRAnalysis.cov_mc(self, orth=False, repl=1000, signif=0.05, seed=None)
This just calls the stderr_MC_irf function from the original model
using the number of periods specified when the irf() class is called
from the VAR class.

Modified functions:

Added specification of error type stderr_type='asym' or 'mc' (see
example). Also repl (replications) and seed options added if 'mc'
errors are specified.


These functions now take the error type specified in the plot()
function above and treats the errors accordingly. Because of how all
of the VAR work (esp. plotting) assumed asymptotic standard errors,
all of these irf plot functions assumed that errors were passed in as
a single matrix with each standard error depending on the MA lag
length and shock variable. Now, depending on which error is specified,
the function will take a tuple of arrays as the standard error if
(stderr_type='mc') rather than a single array.

A serious issue right now is speed. While the asymptotic standard
errors take about a half second on my home laptop to run, the Monte
Carlo standard errors with 3 variables and 1000 replications takes
about 13 seconds to run. Each replication discards the first 100
observations. The most taxing aspect of generating the errors is
re-simulating the data assuming normally distributed standard errors
1000 times (using the util.varsim() function).



  1. Hello, I tried running you code but stumpled on some errors. First, you don't import numpy which you need for np.log(data). Second, your second line, import scikits.statsmodels.tsa.api import VAR; shouldn't the second import be an "as"?

    After correcting these errors, I try to run the code, but get an error on the line, mymodel = VAR(data, names = names) saying that 'module' object is not callable.

    Any ideas? Thanks.

  2. Hi, are you using Bart's fork on github? This might not work with trunk yet.

    The import should be

    from scikits.statsmodels.tsa.api import VAR

    Thanks for the feedback.

  3. Hi,

    Thanks for the feedback. My mistakes. I added import numpy as np and it should be "from scikits.statsmodels.api import VAR".

    If you still have issues, let me know.