Closed NathanielF closed 1 year ago
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.
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
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
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?
@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.
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. (ii) I've applied the same technique to assessing the relationship between GDP and Consumption for Ireland. Ireland is noteworthy in this respect because the corporate tax receipts in recent years have inflated Irish GDPs. The modelling shows that there is a lesser than expected correlation between GDP and consumption in Irelands case. I've compared this result to the one achieved by Statsmodels on the same data (iii) I've fit a non-centrered hierarchical model Bayesian VAR model over 8 countries including Ireland and shown that the estimated relationship between GDP and Consumption is much stronger than we saw in the simple Irish case. Again i show that the model seems to recover the data well in the posterior predictive distribution for all countries.
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 betaX
eterministic andobs
variable in the VAR model so that they show in the graphviz (I think you can't do this for the CholeskyCov and deterministics, maybe we should open an issue for that)
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 obs
and betaX
params now.
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.
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:03Z ----------------------------------------------------------------
Shouldn't this be above the plot?
Also, I would suggest using az.plot_posterior(..., ref_val=true_vals)
to more clearly show that the parameters are being well recovered. It's difficult to compare the az.summary
table with the real values because they are located far from each other in the notebook (also comparing numbers is harder than looking at plots).
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:03Z ----------------------------------------------------------------
When you mention similar magnitudes maybe add some numbers in parenthesis to help guide the user (pymc ~ mean_posterior, statsmodels = x)?
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:04Z ----------------------------------------------------------------
Is the non mv_norm
doing what you want here. It doesn't seem to do the same type of pooling of the noise term? Shouldn't the omega_global
and rho
only be created when mv_norm=True
. I think it's also fine to just not include the mv_norm=False
option here.
NathanielF commented on 2022-12-10T20:17:54Z ----------------------------------------------------------------
I've dropped mvnorm argument here, and removed the alternative likelihood.
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:05Z ----------------------------------------------------------------
Perhaps mention the exact variables / expressions where this is happening to help guiding readers.
by allowing the us to shift the country specific estimates away from the hierarchical mean. This is done by `rho * omega_global + (1 - rho) * noise_col
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:06Z ----------------------------------------------------------------
This one looks terrible when rendered in the docs, perhaps don't include it? https://pymcio--456.org.readthedocs.build/projects/examples/en/456/time_series/bayesian_var_model.html#adding-a-bayesian-twist-hierarchical-vars
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:07Z ----------------------------------------------------------------
You mention (and motivate) posterior predictive plots below, so maybe move this plot over there?
NathanielF commented on 2022-12-10T20:18:11Z ----------------------------------------------------------------
Moved this plot.
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:08Z ----------------------------------------------------------------
I suggest adding a vertical line at 0, to help "situating" the ridgeplots that are closer to the top
NathanielF commented on 2022-12-10T20:18:44Z ----------------------------------------------------------------
Good idea. I've added that now in as a red line.
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:09Z ----------------------------------------------------------------
Would it be too much trouble to plot the posterior of the correlation in a similar fashion as you did above? Would help make it more consistent (also distributions are always nicer than tables) :)
NathanielF commented on 2022-12-10T20:19:30Z ----------------------------------------------------------------
That's fair. Have changed that now and use az.plot_posterior(...)
instead of my clunkier plot.
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:10Z ----------------------------------------------------------------
I would emphasize exactly how the differ. E.g., Correlation between GDF and CONS is now much higher/lower...
View / edit / reply to this conversation on ReviewNB
ricardoV94 commented on 2022-12-09T08:38:11Z ----------------------------------------------------------------
How bad would this plots look with a white border around the black circles? Trying to make them a bit more visible when they fall farther way from the prediction mass (black against dark gray)
NathanielF commented on 2022-12-10T20:20:04Z ----------------------------------------------------------------
Have adopted this suggestion. I think it looks better now too.
This is amazing @NathanielF! I left a bunch of comments in the nbviewer. Not asking for any fundamental changes, just tweaks to the plots / explanations.
Thanks so much @ricardoV94, will work through over the weekend!
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 vector autoregressive model
as tag instead to match better the pattern. sphinx sorts first uppercase then lowercase so combining lower and uppercase looks weird, and not everyone will be privy to all acronyms.
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 ax
to recolor the histogram, I would probably advise against coloring the histogram with a colormap though, there is no colorbar and there seems to be no information encoded in the colormap other than the x.
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 ...
Thank you for this, it's much neater. Also dropped the colormap pattern.
Thanks @OriolAbril, will try and address these points tonight.
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
That's fair. Have changed that now and use az.plot_posterior(...)
instead of my clunkier plot.
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!
Thanks so much for your time and work on this request @ricardoV94 and @OriolAbril.
Thank you @OriolAbril!
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.