Add example on using BVAR model for economists

[twiecki]

Notebook proposal

Title: BVAR model for economists

Why should this notebook be added to pymc-examples?

Economists seem to like BVAR and are looking for this example as a starting point. They are also doing a lot of Bayesian modeling but rarely using PPLs.

[drbenvincent]

@ricardoV94 has done this on the PyMC Labs website https://www.pymc-labs.io/blog-posts/bayesian-vector-autoregression/. Is this issue now irrelevant?

[OriolAbril]

If this is something that we believe is useful and important for users to be part of the docs then we should adapt the notebook and include it to the collection here so it is discoverable by users (who will use the search bar, tags... and generally won't know about the pymc labs blog) and is kept updated as new versions of pymc are released (I assume the blog post won't be updated recurrently but left frozen instead). The adapted notebook should obviously link to the original notebook and give proper attribution which should also drive some traffic to pymc labs blog too I think.

If this is a one-off more obscure functionality that you don't think is needed as part of the pymc examples collection then it's fine as is and there is no need to reopen the issue

[drbenvincent]

tagging @ricardoV94

[twiecki]

Oh I got confused with the repos, should definitely add it here too.

[NathanielF]

Not an economist, but I think this is a cool kind of model is under served in the python infrastructure. I've used the statsmodels version to interrogate questions about the lagged relationship between app performance and customer NPS. It'd be great to see an authoritative example of how to conduct such an analysis in pymc. I think a good write could help popularise this kind of model. Even the statsmodels docs are sparse , how to do X rather than why we might want to do X.

I'm happy to help adapting the pymc labs blog post if that's something you want to pursue?

[drbenvincent]

Ha, the pre-commit checks haven't scared you away!

Go for it 👍🏻

[NathanielF]

Will do a bit of research later this week, and try to pick it up this weekend.

[NathanielF]

So I looked into this a little bit yday and I have a few questions re: (a) the implementation in the pymc labs post and (b) bayesian vars in the literature.

So on (a) I noticed that pymc lab used a normal likelihood. Most packages and the literature seem interested in the covariance matrix to do downstream analysis of impulse response functions. Any reason we didn't use MvNormal? The covariance matrix and moving average representation of AR models should be at least mentioned.

On (b) I found a nice write up by Jim Savage of how Bayesian VARs are helpful for forecasting due to the regularisation of hierarchical priors over short time series e.g. the Minnesota prior.

I think for the write up it would be good to demonstrate this. Am working on trying to convert a stan model, but struggling a bit with shape handling at the moment.

Current idea is to (1) show classical example with statsmodels, (2) show pymc lab model can estimate the same. (3) discuss moving average representation of VARs as beyond scope of post but crucial to know (4) demo hierarchical bvar as argument for why the bayesian part is important.

Any thoughts or questions?

I'll try to share an example or the Bayesian hierarchical var in stan later.

[NathanielF]

Looking at this Stan implementation of a VAR(1) hierarchicall bayesian model:

https://rpubs.com/jimsavage/hierarchical_var

We can see a panoply of priors and hyper priors of various dimensions. I can fit the model to data and retrieve similar parameter estimates to those reported. Broadly the model seems sensible, but i don't know how to translate the auto-regressive step where they seem to just apply an element wise multiplication of the lags with a matrix of beta parameters.

Estimating the model I recover the dimensions for the priors and hyperpriors in stan:

and i can specify a similar set of priors in pymc with the same dimensions the key dimension here is country in the economic data. There are 8 in our data set and three variables.

Which gives a promising looking graph:

but the autoregressive component where they multiply the lags by the beta components will fail in pymc because we can't braodcast:

alpha_8_3 + beta_8_3_3*y_lagged_370_3

and i'm unsure how to best replicate the "double" indexing into country and observation that they are able to do in how they specify the stan likelihood. I feel like i'm missing something crucial, so any tips appreciated.

[ricardoV94]

CC @jessegrabowski

[jessegrabowski]

I have something ready to go on this, I would just need to clean it up a bit and re-focus it away from some of the tangential issues. Thoughts?

https://github.com/jessegrabowski/VAR-Blog/blob/main/Big%20VAR%20Blog%20Post.ipynb

[NathanielF]

That's amazing work @jessegrabowski , you're 10 steps ahead of me on this! Those unit-circle plots in particular are amazing!

I'm happy to bow out if you want to do the write up for this example? But i'd love to see the argument made for why Bayesian VARs in particular are useful as per the example of a hierarchical case discussed in the post above.