accessing/editing Cholesky decomposition matrix

[laurencehendry]

Hi @mboeck11 @hdarjus @edoardochiarotti @eddelbuettel,

Per the get_shockinfo("chol") method in Recursive Identification and GIRFs, please could we ask your guidance on how to access and edit the Cholesky decomposition matrix in the BGVAR R package?

Background to this is a project at the Reserve Bank of Australia (RBA), where we are interested in a method for defining relationships via the lower triangle of our matrix of variance-covariance estimates (a.k.a. the Cholesky method of decomposition). The ordering of the variables in our matrix is important, since the variable in the first variable/row will never be sensitive to a contemporaneous shock of any other variable, and the last row/variable in the system will be sensitive to shocks of all the other variables.

Thank you very much for contributing this great package, also in advance for any consideration. Kind regards, Laurence

[mboeck11]

Hi!

Thanks for your intereset in the BGVAR. Can you elaborate a bit, what do you mean with editing the Cholesky decomposition?

Currently, in the package if you run the irf command, e.g.,

shockinfo_chol <- get_shockinfo("chol") irf <- irf(run, n.ahead = 24, shockinfo = shockinfo_chol)

you can access with irf$struc.obj some structural objects. Specifically, the struc.obj is a list with A, Fmat, Smat and Ginv. Here, Fmat refers to the $F1$ matrix in Eq. (1) in the paper, while Smat($\Sigma$) and Ginv ($G^{-1}$) refer to the matrices building $\Sigma{e} = G^{-1} \Sigma G^{-1'}$. Note that this only contains the posterior median of these matrices since storing the full posterior distribution of these matrices would come with substantial computational cost (actually, the computation of BGVAR would fail on most standard computers because the required memory would be too big).

If you're interested in computing the Cholesky decomposition on your own, you can access the posterior draws directly from the BGVAR run. Suppose, you run a model e.g.,

run <- bgvar(Data = testdata, W = W.test, plag = plag)

you can access the same matrices via run$stacked.results, which is a list containing S_large($\Sigma$), F_large ($F$), and Ginv_large($G^{-1}$). These object are arrays and the third dimension corresponds to the posterior draws.

As a last point, note also how we perform structural inference. We describe this in Section 2.3 of the JSS paper (which I attach here). We identify shocks only locally with recursive ordering within a country, leaving the other elements of the covariance matrix (of other countries) untouched. Then, we premultiply with $G^{-1}$ to get the global solution and back out contemporaneous relationships on a global level.

I hope this helps!

Boeck, Feldkircher, and Huber (2022) JSS.pdf

[laurencehendry]

Hi @mboeck11,

Thank you very much, especially for the guidance on accessing the model’s ex post structural objects (we can access A / Fmat / Ginv / Smat now).

Please could I ask where in the code/package it refers to the ordering of the local variables for Choleski decomposition?

Many thanks, Laurence

[mboeck11]

Hi Laurence,

Happy to help. Sure; the ordering in the Cholesky decomposition is directly determined by the order of the variables in the dataset. So, for instance, in the example dataset (as described in the vignette here of the eerData dataset, which you can access as follows

data(eerData colnames(data$US)

we have the following order of variables imposed: y, Dp, rer, stir, ltir, tb, poil. Here, y is real GDP, Dp is inflation, rer the real effective exchange rate, stir short-term interest rates, ltir long-term interest rates, tb the trade balance, and finally poil the oil price. If you would use the Cholesky decomposition without changing the dataset, then this is also the ordering of the variables. Suppose, we want to identify an monetary policy shock as an exogenous innovation to the short-term interest rate; then y, Dp, and rer are allowed to react contemporaneously to this shock, while the remaining variables do not react contemporaneously. Since we identify the global locally, you only have to change the ordering of the variables of the respective country where the shock originates from. In case you would impose a global shock, the country-ordering also plays a role. We do not have implemented this but given that all the estimated matrices are available, it is rather straightfoward to do.

Hope this helps!

[laurencehendry]

GVAR-introduction-and-demonstration.zip Hi @mboeck11, Thank you (very much!) for your continued support. Per section 1.6 in the attached Markdown (zipped), please could I ask your advice where and/or how to configure the lower triangle of our matrix of variance-covariance estimates? For background, although acknowledging the possibility to configure the ordering of the variables, I remain unsure how to configure the contemporaneous correlations between variables so that they don't correlate. Thank you, Laurence

[mboeck11]

Hi Laurence, I think I don't fully understand your question. If you have a reduced-form VAR process, e.g., $yt = A y{t-1} + u_t$ with $u_t \sim N(0,\Sigma)$, then you're error are correlated. Your goal ist to transform the model in its structural form, such that $B_0 yt = B y{t-1} + e_t$ with $e_t \sim N(0,I)$. Identificaton of this process means that you find a matrix $B_0$, such that $B_0^{-1} e_t = u_t$ and $\Sigma = B_0^{-1} B_0^{-1'}$ holds. To find $B_0^{-1}$ you have to impose identification restrictions. More generally, all matrices $B_0^{-1} = LQ$ with $L$ being the Cholesky decomposition and $Q$ being an orthogonal random matrix drawn from a uniform distribution (or Haar measure), and $QQ'=I$ holds. If you assume that $Q=I$, then this corresponds to the Cholesky decomposition translating into restrictions of the upper diagonal. Economically, this means imposing timing restrictions. This directly corresponds to the global solution of the GVAR, as in Eq (1) of our JSS paper. Again, the BGVAR package provides local identification as opposed to global identification, as dicussed in one of my previous answers.

So, if you want that local shocks do not correlate in the BGVAR, you have to identifiy the structural form. The package provides short-run (timing) restrictions and sign-restrictions. You could of course also work with proxy identification, as recently been used quite heavily in macroeconomics. If you compute impulse response functions with the irf function, and provide the Cholesky scheme, then this decomposition will be used and the shocks will not correlate.

I hope this clarifies, how the procedure works in the BGVAR package. May I also use this opportunity to say that we would appreciate if you cite our JSS paper when using the BGVAR package. Thanks in advance.

[laurencehendry]

Hi @mboeck11,

Thank you for the clarification – if I understood correctly, there are then two theoretical avenues for using a Cholesky decomposition in the BGVAR R package:

(1) Correlated approach (2) Proxy identification by supplying the Cholesky scheme (i.e. with short-run timing restrictions)

Practically, using the BGVAR R package, the short-run timing restriction for a Cholesky decomposition can be introduced through the irf function, per your vignette and the expanded example below.

Please, two follow-up questions

1. When imposing an orthogonal shock using BGVAR::irf() and defining a Cholesky scheme, does a method exist for precisely configuring the timing of these short-run timing restrictions, for when an ordered list of multiple shock variables is defined?
2. Please does any guidance exist for better understanding of the matrices (/objects) stored within the BGVAR bgvar.irf object in R more broadly? (I.e. model.obj and IRF_store; the descriptions provided in your vignette for posterior, also for struc.obj in this thread, have been very helpful!)

On JSS citation – for reference, we are currently three internal staff exploring/testing appropriate R packages at this stage only. We will for sure cite the JSS paper wherever we use the BGVAR package in future!

Many thanks, Laurence

--

OIRF

\ Source: Below code adapted from BGVAR: Bayesian Global Vector Autoregressions with Shrinkage Priors in R (2022) paper by Maximilliian Boeck, Martin Feldkircher, Florian Huber. JSS and CRAN vignette. \ Orthogonalized impulse response functions using a Cholesky decomposition of the variance covariance matrix. \ The direction of this orthogonal shock is:

1. US short-term interest rates (US.stir) first (meaning the orthogonal shock applied commences with this shock);
2. US real GDP (US.y);
3. US inflation (US.Dp).

Note the ordering of the variables here is very important, since the variable in the first row will never be sensitive to a contemporaneous shock of any other variable, and the last row in the list will be sensitive to shocks of all the other variables.

# US economy shocks
shockinfo_chol<-get_shockinfo(ident = "chol",
                              nr_rows = 3) # Quarters

shockinfo_chol$shock<-c("US.stir",
                        "US.y",
                        "US.Dp")

# Indicate scale of shock
shockinfo_chol$scale <-c(100,
                         200,
                         300)

# We are not programming a global shock here, so we default to FALSE for our three shock variables
# shockinfo_chol$global <-c(TRUE,
#                           TRUE,
#                           TRUE)
shockinfo_chol |> attributes()

Examine IRF shock with defined Cholesky decomposition scheme

shockinfo_chol

Run the shock through the model:

irf.chol.us.mp<-BGVAR::irf(model.ssvs.1, 
                    n.ahead=24, # Project out 24 quarters
                    shockinfo=shockinfo_chol, 
                    expert=list(save.store=TRUE)) # Preserve the metadata for our Cholesky shock

Examine shock matrices

names(irf.chol.us.mp)

It is possible to retrieve the matrices used in the global representation of the model, for instance $G^{-1}$ matrix that reveals the $k$ × $k$ dimensional global coefficient matrix used to generate the example GVAR model here.

DT::datatable(irf.chol.us.mp$struc.obj$Ginv)

Alternatively, it is possible to retrieve the matrix of the contemporaneous shocks used if defining a OIRF Cholesky Decomposition shock:

DT::datatable(irf.chol.us.mp$IRF_store[, , 2, 2])

[mboeck11]

Hi Laurence,

No, I think there is a misunderstanding. The Cholesky decomposition is just a decomposition of the variance-covariance matrix. By using this decomposition, economically you introduce timing assumptions. This means, that the first shock in the system moves all variables in the system contemporaneously. The second shock in the system has a zero resetriction for the first variable in the system, i.e., no contemporaneous reaction of the first variable in the system. Timing refers thus to the ordering of the variables: Depending on the ordering of the variables, different variables are allowed to contemporaneously react (or not). The order determines the timing. For instance, for monetary policy shocks, CEE use the following ordering: slow-moving variables (e.g., GDP, inflation), interest rate, fast-moving variables (stock prices). This then translates into: A monetary policy shock does not move slow-moving variables on impact / contemporaneously but fast-moving variables are allowed to react instantenously. I am not sure what you mean with correlated approach, because the structural errors are then by definition uncorrelated.

Regarding your second point of proxy identification: If you have a proxy (instrument) for your shock, then you can order it as the first variable in the system and achieve identification. This has been shown by this paper.. But this has then nothing to do with timing restrictions, but with the mere fact that the proxy is exogenous by definition and all the variables are allowed to react contemporaneously (since it's exogenous).

Regarding your two questions. Again, the ordering of the variables defines the timing assumptions imposed by the Cholesky factorization. I am not sure what you mean with ordered list of multiple shock variables. If you want to re-order the variables, you have to re-estimate the model with a different ordering. Regarding the R objects, they always contain those matrices. If they are labeld equally, then you have the same matrices at hand. But we do not document this in any detail since most users do not have to access those matrices. Otherwise, I would relegate to my earlier answer describing the structure of the model coefficients.

Thanks regarding the citation. I am just pointing this out, because people tend to forget citing packages/articles of packages but a lot of work goes into it.

I hope this clarifies some issues and all the best with your project.

[laurencehendry]

Hi @mboeck11 , thank you very much. I have sent you a message via email.