[pystatsmodels] MC Error Bands for VAR IRF
Bart Baker
I've come to a major checkpoint in integrating Monte Carlo error bands
for impulse response functions.
Here is some quick code to get VAR IRF MC standard errors:
import scikits.statsmodels.api as sm
import 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 = mymodel.fit(maxlags=3,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)
plt.show()
I added functions to the VARResults, IRAnalysis, and also added to
some of the pre-existing functions in these classes and also to
plotting.py. Because Monte Carlo standard errors are in general not
symmetric, I had to alter the plot_with_error function from the
tsa.vector_ar.plotting.py 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.
Functions added to:
BaseIRAnalysis.plot
Added specification of error type stderr_type='asym' or 'mc' (see
example). Also repl (replications) and seed options added if 'mc'
errors are specified.
in tsa.vector_ar.plotting.py
plot_with_error()
irf_grid_plot()
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 series 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).
It doesn't look like SAS or Stata include MC standard errors as an
option for impulse response functions (assuming that what are
generated are asymptotic). Eviews does include it, but it takes way
less time to perform the same task.
Any suggestions?
I'm going to add a blog post with some nice graphics soon (hopefully today).
Bart
Skipper Seabold
Method names should be lower case. See below at the end.
I think you should give the user the option to set disc. I think we're
also trying to call this parameter burn, burnin, nburn, or something.
This is somewhat of a convention for this option across packages.
> 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.
>
Some comments on the IRF code in general. Maybe this question for Wes.
I don't understand why there are two separate classes here with
identical signatures. Also, BaseIRAnalysis is calling cov and cov_mc
from the child class, which suggests there shouldn't be two classes...
What am I missing?
> Functions added to:
>
> BaseIRAnalysis.plot
> Added specification of error type stderr_type='asym' or 'mc' (see
> example). Also repl (replications) and seed options added if 'mc'
> errors are specified.
>
I'm not quite sure yet how we want to best handle seeds. It doesn't
look like it's used yet though. What did scikit-learn settle on?
> in tsa.vector_ar.plotting.py
> plot_with_error()
> irf_grid_plot()
>
> 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 series 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).
>
Let's get everything right, then we can worry about speed. But we
might want to move varsim and the MC stuff into Cython eventually.
I'm not quite sure about the algorithm though. They seem to suggest
that you can just draw from the joint distribution of B and Sigma to
sample from the posterior IRF (I think this is why they say their
method uses the likelihood and the fit) rather than generating new
data and fitting a new VAR. Unless I've misunderstood, which is
entirely possible though. I need to spend more time with the paper. If
we can just draw from some distributions then we can probably avoid a
lot of calls to Python. More on this later.
> It doesn't look like SAS or Stata include MC standard errors as an
> option for impulse response functions (assuming that what are
> generated are asymptotic). Eviews does include it, but it takes way
> less time to perform the same task.
The R code is almost instantaneous. It's written in C++ according to the docs.
http://cran.r-project.org/web/packages/MSBVAR/index.html
I've attached a script to get the IRFs in R, if you want more
comparisons. There are 3 different types described in the docs. They
might describe the algorithm more as well, I didn't look. They do seem
to be a bit more informative in their working paper than the
Econometrica one.
>
> Any suggestions?
>
> I'm going to add a blog post with some nice graphics soon (hopefully today).
>
> Bart
>
Nitpicky stuff.
The TypeError in line 111 here, should be a ValueError. You use these
for bad arguments. Something like
raise ValueError("Give an informative error message")
https://github.com/bartbkr/statsmodels/blob/cum_irf_MC/scikits/statsmodels/tsa/vector_ar/irf.py#L111
http://docs.python.org/library/exceptions.html#exceptions.ValueError
The method stderr_MC_irf should probably be named something like
irf_stderr_mc. Method names should be lower case.
The call signature for this is also over 80 characters. You want 80
characters max for line lengths.
See the style guide (PEP-8): http://www.python.org/dev/peps/pep-0008/
You can run this over your code to remind you of the style
conventions, but it will complain about some things that we don't deem
to be too bad. You can override these to make it not be so picky.
http://pypi.python.org/pypi/pep8
Another general comment that's not important right this second, but
make sure you get in the habit of making docstrings that explain the
API. Don't follow Wes's example ;) Most seem to be okay though. I also
don't have the Lutkepohl book, so page numbers don't help me. I'd
almost prefer that the equation numbers, etc. were comments in the
code rather than in the docs. We can put the actual math in the docs.
The reference is at the top-level and that should be good enough.
Skipper
Wes McKinney
Maybe there don't need to be. I wanted models such as VECM to inherit
from the VAR IRF class where possible. So maybe just need a single
class for now to do VAR IRFs and we can refactor / extend with
additional models.
Ouch, hurting my feelings =P I intend to do something about the lack
of comprehensive docs. Despite my best efforts I like to write code a
lot more than I like to write docs. A project for this summer as I'll
be spending a lot of time on statsmodels and pandas.
> Skipper
>
Skipper Seabold
<snip>
>> Some comments on the IRF code in general. Maybe this question for Wes.
>> I don't understand why there are two separate classes here with
>> identical signatures. Also, BaseIRAnalysis is calling cov and cov_mc
>> from the child class, which suggests there shouldn't be two classes...
>> What am I missing?
>
> Maybe there don't need to be. I wanted models such as VECM to inherit
> from the VAR IRF class where possible. So maybe just need a single
> class for now to do VAR IRFs and we can refactor / extend with
> additional models.
>
I'm just trying to understand the abstraction and how best to fit in
these other ways to get the standard errors. So IRAnalysis is intended
to be VARIRAnalysis(BaseIRAnalysis) or whatever and then we would have
a VECMIRAnalysis(BaseIRAnalysis) as well that overwrites what it needs
to? So really BaseIR is just a template to be inherited and is
intended to be more general than just for VARs? That makes sense.
<snip>
>> Another general comment that's not important right this second, but
>> make sure you get in the habit of making docstrings that explain the
>> API. Don't follow Wes's example ;) Most seem to be okay though. I also
>> don't have the Lutkepohl book, so page numbers don't help me. I'd
>> almost prefer that the equation numbers, etc. were comments in the
>> code rather than in the docs. We can put the actual math in the docs.
>> The reference is at the top-level and that should be good enough.
>
> Ouch, hurting my feelings =P I intend to do something about the lack
> of comprehensive docs. Despite my best efforts I like to write code a
> lot more than I like to write docs. A project for this summer as I'll
> be spending a lot of time on statsmodels and pandas.
>
Ha, nothing personal. I guess most of the stuff I'm talking about is
intended to be private anyway, so I should RTFS. And, yes, writing
documentation isn't fun, but I am usually complaining for
(hypothetical) users. If I pick up something that doesn't tell me what
it's doing (especially wrt assumptions) I tend not to use it
(eViews...). This is something I'm more sensitive too now that I've
torn my hair out trying to replicate results from other packages. The
user interface code is actually fairly well documented though.
Skipper
Bart Baker
Ok, good to know. Just changed it.
>
> I think you should give the user the option to set disc. I think we're
> also trying to call this parameter burn, burnin, nburn, or something.
> This is somewhat of a convention for this option across packages.
Ok, I'll set it to burn for now.
The error bands that I am producing here are the pure MC error bands
from Lutkepohl (2005) 3.7.4 (also see Appendix D). I was planning on
moving onto the Sims and Zha error bands next. The R package that you
list below includes the MC error bands produced here as an option and
the three Sims and Zha types.
>
>> It doesn't look like SAS or Stata include MC standard errors as an
>> option for impulse response functions (assuming that what are
>> generated are asymptotic). Eviews does include it, but it takes way
>> less time to perform the same task.
>
> The R code is almost instantaneous. It's written in C++ according to the docs.
>
> http://cran.r-project.org/web/packages/MSBVAR/index.html
>
> I've attached a script to get the IRFs in R, if you want more
> comparisons. There are 3 different types described in the docs. They
> might describe the algorithm more as well, I didn't look. They do seem
> to be a bit more informative in their working paper than the
> Econometrica one.
Thanks. This will come in handy once I get going on these SIms and Zha
specific error bands.
>
>>
>> Any suggestions?
>>
>> I'm going to add a blog post with some nice graphics soon (hopefully today).
>>
>> Bart
>>
>
> Nitpicky stuff.
>
> The TypeError in line 111 here, should be a ValueError. You use these
> for bad arguments. Something like
>
> raise ValueError("Give an informative error message")
Ok, fixed this. Right, ValueError makes a lot more sense.
>
> https://github.com/bartbkr/statsmodels/blob/cum_irf_MC/scikits/statsmodels/tsa/vector_ar/irf.py#L111
> http://docs.python.org/library/exceptions.html#exceptions.ValueError
>
> The method stderr_MC_irf should probably be named something like
> irf_stderr_mc. Method names should be lower case.
>
> The call signature for this is also over 80 characters. You want 80
> characters max for line lengths.
Ok, fixed.
>
> See the style guide (PEP-8): http://www.python.org/dev/peps/pep-0008/
>
> You can run this over your code to remind you of the style
> conventions, but it will complain about some things that we don't deem
> to be too bad. You can override these to make it not be so picky.
>
> http://pypi.python.org/pypi/pep8
>
> Another general comment that's not important right this second, but
> make sure you get in the habit of making docstrings that explain the
> API. Don't follow Wes's example ;) Most seem to be okay though. I also
> don't have the Lutkepohl book, so page numbers don't help me. I'd
> almost prefer that the equation numbers, etc. were comments in the
> code rather than in the docs. We can put the actual math in the docs.
> The reference is at the top-level and that should be good enough.
Added more comprehensive docstring for irf_stderr_mc.
>
> Skipper
>
Thanks for all your comments. Just pushed.
Bart
Skipper Seabold
Ah, ok. I thought I was going crazy. I don't see these error bands
(via parameter uncertainty?) in R? Why do you say you don't think they
look good versus eViews? They look close to me. Are you looking at the
standard errors for a suitably high number of replications? I don't
see the actual error bands output anywhere. If the std errors, how are
you calculating them (with a dof correction)? Also, now that I think
about it irf_stderr_mc returns the error bands not the standard errors
right?
I'm also not sure that the percentile calculations are entirely
correct with the rounding. Have a look at numpy.percentile or
scipy.stats.scoreatpercentile to handle weighting for the cases that
aren't integers for rounding (or just use scoreatpercentile, we can't
use numpy.percentile because of versions I think, though it might be
nicer here).
It's hard to tell what eViews is doing, and I don't see it mentioned
anywhere in their docs. Unless there is somewhere else to look for
more details? Assuming that the Monte Carlo is the one mentioned in
Lutkepohl and it's not doing something like MCMC on the posterior
distribution, does it say the default significance level? Once we're
sure about the details, I think we should just document what we're
doing.
Changes look good.
Skipper
PS. I just realized I can make source code comments on github, so I'll
be doing that mostly on there from now on unless it warrants
discussion. I thought I had tried this last night but apparently not.
josef...@gmail.com
I don't think this is necessary with a 1000 or a few thousand
replications. For a fixed significance level, e.g. alpha = 0.05, it's
also possible to target that the threshold falls on an integer.
For percentiles far out in the tail some smoothing would be useful but
I don't think it's necessary in this case.
My impression was Bart was so far doing just standard parametric
bootstrap, but I haven't looked carefully.
Josef
(almost done with another round of orthonormal polynomials :)
Bart Baker
Yeah, you are right. They are error bands and not standard errors. I
changed the function names and docstrings to reflect this.
So, in regard to the Eviews output, I think they are in fact
generating MC standard errors and not error bands, meaning the
standard errors are centered somehow (no documentation as far as I
know).
In the MSBVAR R package, when specifying the error band type in the
mc.irf function, I believe that the "Percentile" method is the
standard MC error band method that I produced here.
>
> I'm also not sure that the percentile calculations are entirely
> correct with the rounding. Have a look at numpy.percentile or
> scipy.stats.scoreatpercentile to handle weighting for the cases that
> aren't integers for rounding (or just use scoreatpercentile, we can't
> use numpy.percentile because of versions I think, though it might be
> nicer here).
I looked at scipy.stats.scoreatpercentile, but this generates a value
that could be halfway between two ints. I need this function to
generate an int, so that the index picks out the corresponding matrix
with the given percentiles at each lag length. This would still seem
some sort of rounding procedure to guarantee an int for the index.
>
> It's hard to tell what eViews is doing, and I don't see it mentioned
> anywhere in their docs. Unless there is somewhere else to look for
> more details? Assuming that the Monte Carlo is the one mentioned in
> Lutkepohl and it's not doing something like MCMC on the posterior
> distribution, does it say the default significance level? Once we're
> sure about the details, I think we should just document what we're
> doing.
Good thinking. The Eviews standard errors, are +- 2 standard errors,
so it would make sense that the Monte Carlo bands generated here would
not match up completely. My fault.
Bart