Large Order-Invariant Bayesian VARs with
Stochastic Volatility^††thanks: We would like to thank Chenghan Hou for many helpful comments.

In this section we extend the stochastic volatility specification in Cogley and Sargent (2005) with the goal of constructing an order-invariant model. The key reason why the ordering issue in Cogley and Sargent (2005) occurs is because of the lower triangular assumption for $\mathbf{B}_{0}$. A simple solution to this problem is to avoid this assumption and assume that $\mathbf{B}_{0}$ is an unrestricted square matrix. To that end, let $\mathbf{y}_{t}=(y_{1,t},\ldots,y_{n,t})^{\prime}$ be an $n\times 1$ vector of variables that is observed over the periods $t=1,\ldots,T.$ To fix ideas, consider the following model with zero conditional mean:

where $\mathbf{D}_{t}=\text{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$ is diagonal. Here we assume that $\mathbf{B}_{0}$ is non-singular, but is otherwise unrestricted. Each of the log-volatilities follows a stationary AR(1) process:

for $t=2,\ldots,T$, where $|\phi_{i}|<1$ and the initial condition is specified as $h_{i,1}\sim\mathcal{N}(0,\omega_{i}^{2}/(1-\phi_{i}^{2}))$. Note that the unconditional mean of the AR(1) process is normalized to be zero for identification. We call this the Order Invariant SV (OI-SV) model.

The assumptions in Proposition 1 are stronger than necessary for identification. In fact, Bertsche and Braun (2020) show that having $n-1$ SV processes is sufficient for identification and derive results relating to partial identification which occurs if there are $n-r$ SV processes with $r>1$. Intuitively, allowing for time-varying volatility generates non-trivial autocovariances in the squared reduced form errors. These additional moments help identify $\mathbf{B}_{0}$. We refer readers to Lanne, Lütkepohl, and Maciejowska (2010) and Lewis (2021) for a more detailed discussion of how identification can be achieved in structural VARs via heteroscedasticity.

Below we outline a few theoretical properties of the stochastic volatility model described in (2). First, we show that the model is order invariant, in the sense that if we permute the order of the dependent variables, the likelihood function implied by this new ordering is the same as that of the original ordering, provided that we permute the parameters accordingly. Stacking $\mathbf{y}=(\mathbf{y}_{1}^{\prime},\ldots,\mathbf{y}_{T}^{\prime})^{\prime}$ and $\mathbf{h}=(\mathbf{h}_{1}^{\prime},\ldots,\mathbf{h}_{T}^{\prime})^{\prime}$, note that the likelihood implied by (2) is

where $|\det\mathbf{C}|$ denotes the absolute value of the determinant of the square matrix $\mathbf{C}$. Now, suppose we permute the order of the dependent variables. More precisely, let $\mathbf{P}$ denote an arbitrary permutation matrix of dimension $n$, and we define $\widetilde{\mathbf{y}}_{t}=\mathbf{P}\mathbf{y}_{t}$. By the standard change of variable result, the likelihood of $\widetilde{\mathbf{y}}=(\widetilde{\mathbf{y}}_{1}^{\prime},\ldots,\widetilde{\mathbf{y}}_{T}^{\prime})^{\prime}$ is

where $\widetilde{\mathbf{B}}_{0}=\mathbf{P}\mathbf{B}_{0}\mathbf{P}^{\prime}$, $\widetilde{\mathbf{D}}_{t}=\mathbf{P}\mathbf{D}_{t}\mathbf{P}^{\prime}=\text{diag}(\text{e}^{\widetilde{\mathbf{h}}_{t}})$ and $\widetilde{\mathbf{h}}=(\widetilde{\mathbf{h}}_{1}^{\prime},\ldots,\widetilde{\mathbf{h}}_{T}^{\prime})$ with $\widetilde{\mathbf{h}}_{t}=\mathbf{P}\mathbf{h}_{t}$. Note that we have used the fact that $\mathbf{P}^{-1}=\mathbf{P}^{\prime}$ and $|\det\mathbf{P}|=1$, as all permutation matrices are orthogonal matrices. In addition, since $\mathbf{P}^{-1}\widetilde{\mathbf{y}}_{t}=\mathbf{y}_{t}$, the first line in the above derivation is equal to $p(\mathbf{y}\,|\,\mathbf{B}_{0},\mathbf{h})$. Hence, we have shown that $p(\mathbf{y}\,|\,\mathbf{B}_{0},\mathbf{h})=p(\widetilde{\mathbf{y}}\,|\,\widetilde{\mathbf{B}}_{0},\widetilde{\mathbf{h}})$. We summarize this result in the following proposition.

Next, we derive the unconditional covariance matrix of the data implied by the stochastic volatility model in (2)-(3). More specifically, given the model parameters $\mathbf{B}_{0}$, $\phi_{i}$ and $\omega^{2}_{i},i=1,\ldots,n,$ suppose we generate the log-volatility processes via (3). It is straightforward to compute the implied unconditional covariance matrix of $\mathbf{y}_{t}$.

Since the unconditional distribution of $h_{i,t}$ implied by (3) is $\mathcal{N}(0,\omega_{i}^{2}/(1-\phi_{i}^{2}))$, the unconditional distribution of $\text{e}^{h_{i,t}}$ is log-normal with mean $\text{e}^{\frac{\omega_{i}^{2}}{2(1-\phi_{i}^{2})}}$. Hence, the unconditional mean of $\mathbf{D}_{t}$ is $\mathbf{V}_{\mathbf{D}}=\text{diag}(\text{e}^{\frac{\omega_{1}^{2}}{2(1-\phi_{1}^{2})}},\ldots,\text{e}^{\frac{\omega_{n}^{2}}{2(1-\phi_{n}^{2})}})$. It then follows that the unconditional covariance matrix of $\mathbf{y}_{t}$ can be written as

where $\overline{\mathbf{B}}_{0}=\mathbf{V}_{\mathbf{D}}^{-\frac{1}{2}}\mathbf{B}_{0}$, i.e., scaling each row of $\mathbf{B}_{0}$ by $\text{e}^{-\frac{\omega_{i}^{2}}{4(1-\phi_{i}^{2})}}$, $i=1,\ldots,n$. From this expression it is also clear that if we permute the order of the dependent variables via $\widetilde{\mathbf{y}}_{t}=\mathbf{P}\mathbf{y}_{t}$, the unconditional covariance matrix of $\widetilde{\mathbf{y}}_{t}$ can be obtained by permuting the rows and columns of $\overline{\mathbf{B}}_{0}$ accordingly. More precisely:

where $\widetilde{\overline{\mathbf{B}}}_{0}=\mathbf{P}\overline{\mathbf{B}}_{0}\mathbf{P}^{\prime}$. Hence, the unconditional variance of any element in $\mathbf{y}_{t}$ does not depend on its position in the $n$-tuple.

This establishes order-invariance in the likelihood function defined by the OI-SV. Of course, the posterior will also be order invariant (subject to permutations and sign switches) if the prior is. This will hold in any reasonable case. All it requires is that the prior for the parameters in equation $i$ under one ordering becomes the prior for equation $j$ in a different ordering if variable $i$ in the first ordering becomes variable $j$ in the second ordering.

We now add the VAR part to our OI-SV model, leading to the OI-VAR-SV which is given by

where $\mathbf{a}$ is an $n\times 1$ vector of intercepts, $\mathbf{A}_{1},\ldots,\mathbf{A}_{p}$ are $n\times n$ matrices of VAR coefficients, $\mathbf{B}_{0}$ is non-singular but otherwise unrestricted $n\times n$ matrix, and $\mathbf{D}_{t}=\text{diag}(\text{e}^{h_{1,t}},\ldots,\text{e}^{h_{n,t}})$ is diagonal. Finally, each of the log-volatility $h_{i,t}$ follows the stationary AR(1) process as specified in (3).

We argued previously that ordering issues are likely to be most important in high dimensional models. With large VARs over-parameterization concerns can be substantial and, for this reason, Bayesian methods involving shrinkage priors are commonly used. In this sub-section, we describe a particular VAR prior with attractive properties. But it is worth noting that any of the standard Bayesian VAR priors (e.g. the Minnesota prior) could have been used and the resulting model would still have been order invariant.

Let $\boldsymbol{\alpha}_{i}$ denote the VAR coefficients in the $i$-th equation, $i=1,\ldots,n$. We consider the Minnesota-type adaptive hierarchical priors proposed in Chan (2021), which have the advantages of both the Minnesota priors (e.g., rich prior beliefs such as cross-variable shrinkage) and modern adaptive hierarchical priors (e.g., heavy-tails and substantial mass around 0). In particular, we construct a Minnesota-type horseshoe prior as follows. For $\alpha_{i,j}$, the $j$-th coefficient in the $i$-th equation, let $\kappa_{i,j}=\kappa_{1}$ if it is a coefficient on an ‘own lag’ and let $\kappa_{i,j}=\kappa_{2}$ if it is a coefficient on an ‘other lag’. Given the constants $C_{i,j}$ defined below, consider the horseshoe prior on the VAR coefficients (excluding the intercepts) $\alpha_{i,j},i=1,\ldots,n,j=2,\ldots,k$ with $k=np+1$:

where $\mathcal{C}^{+}(0,1)$ denotes the standard half-Cauchy distribution. $\kappa_{1}$ and $\kappa_{2}$ are the global variance components that are common to, respectively, coefficients of own and other lags, whereas each $\psi_{i,j}$ is a local variance component specific to the coefficient $\alpha_{i,j}$. Lastly, the constants $C_{i,j}$ are obtained as in the Minnesota prior, i.e.,

where $s_{r}^{2}$ denotes the sample variance of the residuals from an AR(4) model for the variable $r,r=1,\ldots,n$. Furthermore, for data in growth rates $m_{i,j}$ is set to be 0; for level data, $m_{i,j}$ is set to be zero as well except for the coefficient associated with the first own lag, which is set to be one.

It is easy to verify that if all the local variances are fixed, i.e., $\psi_{i,j}\equiv 1$, then the Minnesota-type horseshoe prior given in (5)–(7) reduces to the Minnesota prior. Hence, it can be viewed as an extension of the Minnesota prior by introducing a local variance component such that the marginal prior distribution of $\alpha_{i,j}$ has heavy tails. On the other hand, if $m_{i,j}=0$, $C_{i,j}=1$ and $\kappa_{1}=\kappa_{2}$, then the Minnesota-type horseshoe prior reduces to the standard horseshoe prior where the coefficients have identical distributions. Therefore, it can also be viewed as an extension of the horseshoe prior (Carvalho, Polson, and Scott, 2010) that incorporates richer prior beliefs on the VAR coefficients, such as cross-variable shrinkage, i.e., shrinking coefficients on own lags differently than other lags.

To facilitate estimation, we follow Makalic and Schmidt (2016) and use the following latent variable representations of the half-Cauchy distributions in (6)-(7):

for $i=1,\ldots,n,j=2,\ldots,k$ and $l=1,2$ where $\mathcal{IG}$ denotes the inverse Gamma distribution

Next, we specify the priors on other model parameters. Let $\mathbf{b}_{i}$ denote the $i$-th row of $\mathbf{B}_{0}$ for $i=1,\ldots,n,$ i.e., $\mathbf{B}_{0}^{\prime}=(\mathbf{b}_{1},\ldots,\mathbf{b}_{n})$. We consider independent priors on $\mathbf{b}_{i},i=1,\ldots,n$:

For the parameters in the stochastic volatility equations, we assume the priors for $j=1,\ldots,n$:

In this section we develop a posterior sampler which allows for Bayesian estimation of the order-invariant stochastic volatility model with the Minnesota-type horseshoe prior. For later reference, let $\boldsymbol{\psi}_{i}=(\psi_{i,1},\ldots,\psi_{i,k})^{\prime}$ denote the local variance components associated with $\boldsymbol{\alpha}_{i}$ and define $\boldsymbol{\kappa}=(\kappa_{1},\kappa_{2})^{\prime}$. Furthermore, let $\mathbf{z}_{\boldsymbol{\psi}_{i}}=(z_{\psi_{i,1}},\ldots,z_{\psi_{i,k}})^{\prime}$ denote the latent variables corresponding to $\boldsymbol{\psi}_{i}$ and similarly define $\mathbf{z}_{\boldsymbol{\kappa}}$. Next, stack $\mathbf{y}=(\mathbf{y}_{1}^{\prime},\ldots,\mathbf{y}_{T}^{\prime})^{\prime}$, $\boldsymbol{\alpha}=(\boldsymbol{\alpha}_{1}^{\prime},\ldots,\boldsymbol{\alpha}_{n}^{\prime})^{\prime}$, $\boldsymbol{\psi}=(\boldsymbol{\psi}_{1}^{\prime},\ldots,\boldsymbol{\psi}_{n}^{\prime})^{\prime}$, $\mathbf{z}_{\boldsymbol{\psi}}=(\mathbf{z}_{\boldsymbol{\psi}_{1}}^{\prime},\ldots,\mathbf{z}_{\boldsymbol{\psi}_{n}}^{\prime})^{\prime}$ and $\mathbf{h}=(\mathbf{h}_{1}^{\prime},\ldots,\mathbf{h}_{T}^{\prime})^{\prime}$ with $\mathbf{h}_{t}=(h_{1,t},\ldots,h_{n,t})^{\prime}$. Finally, let $\mathbf{h}_{i,\boldsymbol{\cdot}}=(h_{i,1},\ldots,h_{i,T})^{\prime}$ represent the vector of log-volatility for the $i$-th equation, $i=1,\ldots,n$. Then, posterior draws can be obtained by sampling sequentially from:

1. 1.

   $p(\mathbf{B}_{0}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\psi},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})$ = $p(\mathbf{B}_{0}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{h})$;

2. 2.

   $p(\boldsymbol{\alpha}\,|\,\mathbf{y},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\psi},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})$ = $p(\boldsymbol{\alpha}\,|\,\mathbf{y},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\psi},\boldsymbol{\kappa})$;

3. 3.

   $p(\boldsymbol{\psi}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})$ = $\prod_{i=1}^{n}\prod_{j=2}^{k}p(\psi_{i,j}\,|\,\alpha_{i,j},\boldsymbol{\kappa},z_{\psi_{i,j}})$;

4. 4.

   $p(\boldsymbol{\kappa}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\psi},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})$ = $\prod_{l=1}^{2}p(\kappa_{l}\,|\,\boldsymbol{\alpha},\boldsymbol{\psi},\mathbf{z}_{\kappa_{l}})$;

5. 5.

   $p(\mathbf{z}_{\boldsymbol{\psi}}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\kappa},\boldsymbol{\psi},\mathbf{z}_{\boldsymbol{\kappa}})$ = $\prod_{i=1}^{n}\prod_{j=2}^{k}p(z_{\psi_{i,j}}\,|\,\psi_{i,j})$;

6. 6.

   $p(\mathbf{z}_{\boldsymbol{\kappa}}\,|\,\mathbf{y},\boldsymbol{\alpha},\mathbf{B}_{0},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\kappa},\boldsymbol{\psi},\mathbf{z}_{\boldsymbol{\psi}})$ = $\prod_{l=1}^{2}p(z_{\kappa_{l}}\,|\,\kappa_{l})$.

7. 7.

   $p(\mathbf{h}\,|\,\mathbf{y},\mathbf{B}_{0},\boldsymbol{\alpha},\boldsymbol{\phi},\boldsymbol{\omega}^{2},\boldsymbol{\psi},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})=\prod_{i=1}^{n}p(\mathbf{h}_{i,\boldsymbol{\cdot}}\,|\,\mathbf{y},\mathbf{B}_{0},\boldsymbol{\alpha},\boldsymbol{\phi},\boldsymbol{\omega}^{2})$;

8. 8.

   $p(\boldsymbol{\omega}^{2}\,|\,\mathbf{y},\mathbf{B}_{0},\boldsymbol{\alpha},\mathbf{h},\boldsymbol{\phi},\boldsymbol{\psi},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})=\prod_{i=1}^{n}p(\omega_{i}^{2}\,|\,\mathbf{h}_{i,\boldsymbol{\cdot}},\phi_{i})$;

9. 9.

   $p(\boldsymbol{\phi}\,|\,\mathbf{y},\mathbf{B}_{0},\boldsymbol{\alpha},\mathbf{h},\boldsymbol{\omega}^{2},\boldsymbol{\psi},\boldsymbol{\kappa},\mathbf{z}_{\boldsymbol{\psi}},\mathbf{z}_{\boldsymbol{\kappa}})=\prod_{i=1}^{n}p(\phi_{i}\,|\,\mathbf{h}_{i,\boldsymbol{\cdot}},\sigma^{2}_{i})$;

Step 1. To implement Step 1, we adapt the sampling approach in Waggoner and Zha (2003) and Villani (2009) to the setting with stochastic volatility. More specifically, we aim to draw each row of $\mathbf{B}_{0}$ given all other rows. To that end, we first rewrite (4) as:

where $\mathbf{Y}$ is the $T\times n$ matrix of dependent variables, $\mathbf{X}$ is the $T\times k$ matrix of lagged dependent variables with $k=1+np$, $\mathbf{A}=(\mathbf{a},\mathbf{A}_{1},\ldots,\mathbf{A}_{n})^{\prime}$ is the $k\times n$ matrix of VAR coefficients and $\mathbf{E}$ is the $T\times n$ matrix of errors. Then, for $i=1,\ldots,n$, we have

where $\mathbf{E}_{i}$ is the $i$-th column of $\mathbf{E}$ and $\boldsymbol{\Omega}_{\mathbf{h}_{i,\boldsymbol{\cdot}}}=\text{diag}(\text{e}^{h_{i,1}},\ldots,\text{e}^{h_{i,T}}).$ Hence, the full conditional distribution of $\mathbf{b}_{i}$ is given by

where

The above full conditional posterior distribution is non-standard and direct sampling is infeasible. Fortunately, Waggoner and Zha (2003) develop an efficient algorithm to sample $\mathbf{b}_{i}$ in the case when $\widehat{\mathbf{b}}_{i}=\mathbf{0}.$ Villani (2009) further generalizes this sampling approach for non-zero $\widehat{\mathbf{b}}_{i}$.

The key idea is to show that $\mathbf{b}_{i}$ has the same distribution as a linear transformation of an $n$-vector that consists of one absolute normal and $n-1$ normal random variables. A random variable $Z$ follows the absolute normal distribution $\mathcal{AN}(\mu,\rho)$ if it has the density function

where $c$ is a normalizing constant.

To formally state the results, let $\mathbf{C}_{i}$ denote the Cholesky factor of $T^{-1}\mathbf{K}_{\mathbf{b}_{i}}$ such that $\mathbf{K}_{\mathbf{b}_{i}}=T\mathbf{C}_{i}\mathbf{C}_{i}^{\prime}$. Let $\mathbf{F}_{-i}$ represent the matrix $\mathbf{F}$ with the $i$-th column deleted and $\mathbf{F}^{\perp}$ be the orthogonal complement of $\mathbf{F}$. Furthermore, define $\mathbf{v}_{1}=\mathbf{C}_{i}^{-1}(\mathbf{B}_{0,-i}^{\prime})^{\perp}/||\mathbf{C}_{i}^{-1}(\mathbf{B}_{0,-i}^{\prime})^{\perp}||$, where $(\mathbf{B}_{0,-i}^{\prime})^{\perp}\equiv((\mathbf{B}_{0}^{\prime})_{-i})^{\perp}$, and let $(\mathbf{v}_{2},\ldots,\mathbf{v}_{n})=\mathbf{v}_{1}^{\perp}$. Then, we claim that

where $\xi_{1}\sim\mathcal{AN}(\widehat{\xi}_{1},T^{-1})$ and $\xi_{j}\sim\mathcal{N}(\widehat{\xi}_{j},T^{-1}),j=2,\ldots,n$ with $\widehat{\xi}_{j}=\widehat{\mathbf{b}}_{i}^{\prime}\mathbf{C}_{i}\mathbf{v}_{j}$.

To prove the claim, let $\mathbf{b}_{i}=(\mathbf{C}_{i}^{\prime})^{-1}\sum_{j=1}^{n}\xi_{j}\mathbf{v}_{j}$. It suffices to show that if we substitute $\mathbf{b}_{i}$ into (12), $\xi_{1},\ldots,\xi_{n}$ are independent with $\xi_{1}\sim\mathcal{AN}(\widehat{\xi}_{1},T^{-1})$ and $\xi_{j}\sim\mathcal{N}(\widehat{\xi}_{j},T^{-1}),j=2,\ldots,n$. To that end, note that by construction, $\mathbf{v}_{1},\ldots,\mathbf{v}_{n}$ forms an orthonormal basis of $\mathbb{R}^{n}$, particularly $\sum_{j=1}^{n}\mathbf{v}_{j}\mathbf{v}_{j}^{\prime}=\mathbf{I}_{n}$. We first write the quadratic form in (12) as:

Next, note that by construction, $(\mathbf{v}_{2},\ldots,\mathbf{v}_{n})$ spans the same space as $\mathbf{B}_{0,-i}^{\prime}$. Hence, it follows that

Finally, substituting the quadratic form and the determinant into (12), we obtain

In other words, $\xi_{1}\sim\mathcal{AN}(\widehat{\xi}_{1},T^{-1})$ and $\xi_{j}\sim\mathcal{N}(\widehat{\xi}_{j},T^{-1}),j=2,\ldots,n$, and we have proved the claim.

In order to use the above result, we need an efficient way to sample from $\mathcal{AN}(\widehat{\xi}_{1},T^{-1})$. This can be done by using the 2-component normal mixture approximation considered in Appendix C of Villani (2009). In addition, the orthogonal complement of $\mathbf{v}_{1}$, namely, $\mathbf{v}_{1}^{\perp}$, can be obtained using the singular value decomposition.²²2In MATLAB, the orthogonal complement of $\mathbf{v}_{1}$ can be obtained using null$(\mathbf{v}_{1}^{\prime})$. In the sampler we fix the sign of the $i$-th element of $\mathbf{b}_{i}$ to be positive (so that the diagonal elements of $\mathbf{B}_{0}$ are positive). This is done by multiplying the draw $\mathbf{b}_{i}$ by $-1$ if its $i$-th element is negative.

Step 2. We sample the reduced-form VAR coefficients row by row along the lines of Carriero, Clark, and Marcellino (2019). In particular, we extend the triangular algorithm in Carriero, Chan, Clark, and Marcellino (2021) that is designed for a lower triangular impact matrix $\mathbf{B}_{0}$ to the case where $\mathbf{B}_{0}$ is a full matrix. To that end, define $\mathbf{A}_{i=0}$ to be a $k\times n$ matrix that has exactly the same elements as $\mathbf{A}$ except for the $i$-th column, which is set to be zero, i.e., $\mathbf{A}_{i=0}=(\boldsymbol{\alpha}_{1},\ldots,\boldsymbol{\alpha}_{i-1},\mathbf{0},\boldsymbol{\alpha}_{i+1},\ldots,\boldsymbol{\alpha}_{n}).$ Then, we can rewrite (4) as

where $\mathbf{x}_{t}=(1,\mathbf{y}_{t-1}^{\prime},\ldots,\mathbf{y}_{t-p}^{\prime})^{\prime}$ and $\mathbf{B}_{0,i}$ is the $i$-th column of $\mathbf{B}_{0}$. Let $\mathbf{z}^{i}_{t}=\mathbf{B}_{0}(\mathbf{y}_{t}-\mathbf{A}_{i=0}^{\prime}\mathbf{x}_{t})$ and stack $\mathbf{z}^{i}=(\mathbf{z}^{i^{\prime}}_{1},\ldots,\mathbf{z}^{i^{\prime}}_{T})^{\prime}$, we have

where $\mathbf{W}^{i}=\mathbf{X}\otimes\mathbf{B}_{0,i}$ and $\mathbf{D}=\text{diag}(\mathbf{D}_{1},\ldots,\mathbf{D}_{T})$.

Next, the Minnesota-type horseshoe prior on $\boldsymbol{\alpha}_{i}$ is conditionally Gaussian given the local variance component $\boldsymbol{\psi}_{i}$. More specifically, we can rewrite the conditional prior on $\boldsymbol{\alpha}_{i}$ in (5) as:

where $\mathbf{m}_{i}=(m_{i,1},\ldots,m_{i,k})^{\prime}$ and $\mathbf{V}_{\boldsymbol{\alpha}_{i}}=\text{diag}(C_{i,1},\kappa_{i,2}\psi_{i,2}C_{i,2}\ldots,\kappa_{i,k}\psi_{i,k}C_{i,k})$ (the prior on the intercept is Gaussian). Then, by standard linear regression results, we have

where $\boldsymbol{\alpha}_{-i}=(\boldsymbol{\alpha}_{1}^{\prime},\ldots,\boldsymbol{\alpha}_{i-1}^{\prime},\boldsymbol{\alpha}_{i+1}^{\prime},\ldots,\boldsymbol{\alpha}_{n}^{\prime})^{\prime}$,

The computational complexity of this step is the same, i.e., $\mathcal{O}(n^{4})$, as in the triangular case considered in Carriero, Chan, Clark, and Marcellino (2021).

Step 3. First, note that the elements of $\boldsymbol{\psi}_{i}$ are conditionally independent and we can sample them one by one without loss of efficiency. Next, combining (5) and (9), we obtain