Add Hierarchical BVAR example #456
Conversation
…example Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
Check out this pull request on See visual diffs & provide feedback on Jupyter Notebooks. Powered by ReviewNB |
|
View / edit / reply to this conversation on ReviewNB jessegrabowski commented on 2022-11-02T14:15:29Z I don't think you need to credit me, I'd rather just help on the commit.
Also there's a typo in "fake fata". NathanielF commented on 2022-11-02T14:19:57Z Sure thing. Thanks |
|
View / edit / reply to this conversation on ReviewNB jessegrabowski commented on 2022-11-02T14:15:31Z Is this data real or fake? If it's fake, I suggest to not use actual country names, because it's confusing. If it's real, why have the fake data part above? NathanielF commented on 2022-11-02T14:21:44Z This is real data retrieved using the R package WDI:
NathanielF commented on 2022-11-02T14:39:16Z The reason for the inclusion of the fake data was to just demonstrate that the pattern really works to pull back the actual parameters governing the data generating process. jessegrabowski commented on 2022-11-02T14:45:02Z Yeah good point. I think that could just be more clear, perhaps by segmenting the notebook into two parts. The first part would go over generating the fake data and fitting the model on it. When you make the posterior plots, you could then add some fat red lines showing the true parameter values to really emphasize that it recovered them.
The second part could largely repeat the exercise (highlighting that different values of n_eqs and n_lags have been chosen), and the comparison could be to the statsmodels baseline instead of the known true parameters.
It might also be nice to set some zeros into the coefficient matrix, to show how a researcher could easily implement a restricted VAR (since those are important). |
|
View / edit / reply to this conversation on ReviewNB jessegrabowski commented on 2022-11-02T14:15:32Z Can we just call GFCF Investment? It might also be more interesting to look at a financial variable (like the interest rate) as the 3rd variable, since Y, C, and I are deterministically connected via accounting identities.
NathanielF commented on 2022-11-02T14:23:12Z I'm open to that. The dataset was just pulled using the pattern here: https://rpubs.com/jimsavage/hierarchical_var
But if you think there is a more interesting set of variables to look at, we can do that too. |
|
View / edit / reply to this conversation on ReviewNB jessegrabowski commented on 2022-11-02T14:15:34Z In econometric literature, "beta" always refers to the parameters. Here it seems it refers to the conditional mean of the VAR model -- I think the "make_betas" function should be called something like "compute_mean" or "autoregressive_mean".
NathanielF commented on 2022-11-02T14:24:59Z Sure. That's an easy enough change. |
|
View / edit / reply to this conversation on ReviewNB jessegrabowski commented on 2022-11-02T14:15:35Z It's not very clear which data is being fit here. I guess it's the fake data? I would just pull out one country and fit that instead, eliminating the fake data entirely.
NathanielF commented on 2022-11-02T14:25:44Z Yeah, it was the fake data. Just wanted to demonstrate how the functions could handle model specifications for different inputs of n_lags and n_equations. Will tidy up to make clearer. |
|
Sure thing. Thanks View entire conversation on ReviewNB |
|
This is real data retrieved using the R package WDI:
View entire conversation on ReviewNB |
|
I'm open to that. The dataset was just pulled using the pattern here: https://rpubs.com/jimsavage/hierarchical_var
But if you think there is a more interesting set of variables to look at, we can do that too. View entire conversation on ReviewNB |
|
Sure. That's an easy enough change. View entire conversation on ReviewNB |
|
Yeah, it was the fake data. Just wanted to demonstrate how the functions could handle model specifications for different inputs of n_lags and n_equations View entire conversation on ReviewNB |
|
The reason for the inclusion of the fake data was to just demonstrate that the pattern really works to pull back the actual parameters governing the data generating process. View entire conversation on ReviewNB |
|
Yeah good point. I think that could just be more clear, perhaps by segmenting the notebook into two parts. The first part would go over generating the fake data and fitting the model on it. When you make the posterior plots, you could then add some fat red lines showing the true parameter values to really emphasize that it recovered them.
The second part could largely repeat the exercise (highlighting that different values of n_eqs and n_lags have been chosen), and the comparison could be to the statsmodels baseline instead of the known true parameters.
View entire conversation on ReviewNB |
…nd fake data and irish example compared to statsmodels. Included example of small hierarchical case. Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
Just a note here. I've not yet incorporated all your suggestions. I have restructured it to have a fake data section and an Ireland section where I compared the output to statsmodels. I've also included a third example of a hierarchical VAR which converged pretty well. I think we should stick to the plan of 3 posts, so I will remove this just wanted to show that the code is broadly sound. Note I used one global prior for the LKJ prior, rather than using the rho pattern to modulate it at the country level. I still need to do more of an introductory write up on the mathematical nature of a VAR model and maybe a discussion of their typical use cases. On the comparison to statsmodels I think I need to tweak my priors a bit more. I left them at the same priors used with the fake data and it still broadly converged to the "same" store as statsmodels, but it is too uncertain. Other than that... I think it is taking shape. Curious about your thoughts on the hierarchical example? |
…imated hierarchical model using blackjax nuts sampler Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
@jessegrabowski I think the above notebook is basically complete in so far as it is has all the pieces i want in it: (1) fake data demo, (2) ireland gdp and cons VAR and (3) simple hierarchical VAR with 10 countries. The hierarchical VAR fits quickly now ~30mins and convergence diagnostics look good on all models. Posterior predictive checks seem sound too... One question i had was the comparison with statsmodels: the beta coefficients are markedly different in that e.g. the BVAR fit deems the Gdp <- Gdp_lag_1 ~0, whereas statsmodels has it to be non negative. Not sure if it makes economic sense to have the prior gdp be so weakly related to current GDP, but.... it could just be because of Ireland's weird corporate tax windfalls.....? In any case, if you can have a look and let me know your thoughts that'd be good. I'm fairly happy with it. I can't add the inference data object because it is too big, but i'm hoping i won't need to since the model converges now in a reasonable time. It should be straightforward for you to recreate and explore the Impulse response functions. |
…ies to hierarchical model fit Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
… fit. Fixed some typos too. Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
… and improved discussion. Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
…image for full model Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
change filepath to full_model.png
Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
This pull request is intended as part of a series with @jessegrabowski where we demonstrate the particularities and capabilities of Bayesian VAR models. In this (the first of the series) i've demonstrated how to build and estimate a simple Bayesian VAR and a hierarchical Bayesian VAR. Drawing on the Pymc labs post I've shown: (i) How to build a simple Bayesian VAR with 2 lags on fake data and recover the same data. I think all the ingredients are there to start a discussion in the next post of the series about how to use the implied posterior distribution over the correlation parameters to (i) model the impulse response function and (ii) build forecasts from VAR models. I'm marking this pull request as ready for review mainly because there is allot here already. I'm open to suggestions of re-writes or a different focus, but i'd like to avoid displaying more functionality in this particular post because we already have covered allot of ground.... @lucianopaz, @ricardoV94 can i request a review from whomever is most appropriate? |
|
View / edit / reply to this conversation on ReviewNB ricardoV94 commented on 2022-12-09T08:38:00Z Make the PyMC labs a reference, similar to how it shows in this other example: https://pymcio--456.org.readthedocs.build/projects/examples/en/456/causal_inference/excess_deaths.html#references NathanielF commented on 2022-12-10T20:13:59Z This is resolved now: |
|
View / edit / reply to this conversation on ReviewNB ricardoV94 commented on 2022-12-09T08:38:01Z You can add dims to the OriolAbril commented on 2022-12-09T21:56:40Z You can already define dimentions in deterministics. You can't yet for cholesky but you can still define their dims via idata kwargs. https://www.pymc.io/projects/examples/en/latest/case_studies/multilevel_modeling.html does it (search for idata_kwargs to find the model that does this). NathanielF commented on 2022-12-10T20:15:09Z Added in dims for the |
|
View / edit / reply to this conversation on ReviewNB ricardoV94 commented on 2022-12-09T08:38:02Z Why the different colormap here, compared to the ones below? OriolAbril commented on 2022-12-09T22:22:58Z Also, you can plot dots with borders in a different color without needing to call plot two times, you can use markerdegecolor and markeredgewidth (or shortened aliases like below):
plt.plot(..., "o", mfc="black", mec="white", mew=2)
Ref: https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.plot.html (you need to scroll down to the "other parameters" section to see al options available) NathanielF commented on 2022-12-10T20:17:12Z reverted this plot to the same color scheme as the rest, and used your tip about white edge on markers in matplotlib. Reserved this colour scheme just for Ireland in the final plot to highlight it against other countries. |
Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
You can already define dimentions in deterministics. You can't yet for cholesky but you can still define their dims via idata kwargs. https://www.pymc.io/projects/examples/en/latest/case_studies/multilevel_modeling.html does it (search for idata_kwargs to find the model that does this). View entire conversation on ReviewNB |
|
Also, you can plot dots with borders in a different color without needing to call plot two times, you can use markerdegecolor and markeredgewidth (or shortened aliases like below):
plt.plot(..., "o", mfc="black", mec="white", mew=2)
Ref: https://matplotlib.org/stable/api/_as_gen/matplotlib.axes.Axes.plot.html (you need to scroll down to the "other parameters" section to see al options available) View entire conversation on ReviewNB |
|
View / edit / reply to this conversation on ReviewNB OriolAbril commented on 2022-12-09T22:26:41Z I would use NathanielF commented on 2022-12-10T20:20:31Z Have added this. |
|
View / edit / reply to this conversation on ReviewNB OriolAbril commented on 2022-12-09T22:26:42Z I think all of this can be reduced to
ax = az.plot_posterior(
idata,
var_names="noise_chol_corr",
hdi_prob="hide",
point_estimate=None,
grid=(2,2),
labeller=az.labels.NoVarLabeller(),
kind="hist",
ec="black", # to keep similar look
)
Plus potentially looping over elements in
I believe you could loop over the patches with:
for ax in ax.ravel():
left_bin_edges = [p.get_x() for p in ax.patches]
cs = left_bin_edges + (left_bin_edges[1] - left_bin_edges[0]) / 2
...
NathanielF commented on 2022-12-10T20:21:09Z Thank you for this, it's much neater. Also dropped the colormap pattern. |
|
Thanks @OriolAbril, will try and address these points tonight. |
…improved prior of model for statsmodels comparison Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
This is resolved now: View entire conversation on ReviewNB |
|
Added in dims for the View entire conversation on ReviewNB |
|
reverted this plot to the same color scheme as the rest, and used your tip about white edge on markers in matplotlib. Reserved this colour scheme just for Ireland in the final plot to highlight it against other countries. View entire conversation on ReviewNB |
|
I've dropped mvnorm argument here, and removed the alternative likelihood. View entire conversation on ReviewNB |
|
Moved this plot. View entire conversation on ReviewNB |
|
Good idea. I've added that now in as a red line. View entire conversation on ReviewNB |
|
That's fair. Have changed that now and use View entire conversation on ReviewNB |
|
Have adopted this suggestion. I think it looks better now too. View entire conversation on ReviewNB |
|
Have added this. View entire conversation on ReviewNB |
|
Thank you for this, it's much neater. Also dropped the colormap pattern. View entire conversation on ReviewNB |
|
I think it will be ready to merge after fixing the tag. Thanks so much for the notebook! |
Signed-off-by: Nathaniel <NathanielF@users.noreply.github.com>
|
Thanks so much for your time and work on this request @ricardoV94 and @OriolAbril. |
|
Thank you @OriolAbril! |
OriolAbril
changed the title
Signed-off-by: Nathaniel NathanielF@users.noreply.github.com
Example of Hierarchical Bayesian VAR model
An example adapting the PYMC labs by @ricardoV94 post on the same topic and which includes a generalisation of the functionality seen there, plus an application of the hierarchical modelling strategy to the VAR model. I will also add a discussion section on the benefits of the Bayesian approach for these kinds of model, but wanted to share this for the purposes of discussion and refinement.
@jessegrabowski has already done some amazing work on explaining VAR modeling and showing how they are a type of convolution. I see this example as building on the discussion in his work and trying to exhibit one of the key advantages of the Bayesian approach to this kind of modelling.
One thing to note is that the while my hierarchical model here samples. It is quite slow, and it'd be good to improve that somehow.