# Fast and Order-invariant Inference in Bayesian VARs with Non-Parametric
Shocks††thanks:  Huber gratefully acknowledges financial support from the
Austrian Science Fund (FWF, grant no. ZK 35).

Florian Huber University of Salzburg International Institute for Applied
Systems Analysis (IIASA) Gary Koop University of Strathclyde

Abstract

The shocks which hit macroeconomic models such as Vector Autoregressions
(VARs) have the potential to be non-Gaussian, exhibiting asymmetries and fat
tails. This consideration motivates the VAR developed in this paper which uses
a Dirichlet process mixture (DPM) to model the shocks. However, we do not
follow the obvious strategy of simply modeling the VAR errors with a DPM since
this would lead to computationally infeasible Bayesian inference in larger
VARs and potentially a sensitivity to the way the variables are ordered in the
VAR. Instead we develop a particular additive error structure inspired by
Bayesian nonparametric treatments of random effects in panel data models. We
show that this leads to a model which allows for computationally fast and
order-invariant inference in large VARs with nonparametric shocks. Our
empirical results with nonparametric VARs of various dimensions shows that
nonparametric treatment of the VAR errors is particularly useful in periods
such as the financial crisis and the pandemic.

JEL Codes: C11, C32, C53

Keywords: Bayesian VARs, infinite mixtures, fast estimation, Markov chain
Monte Carlo.

##  1 Introduction

Bayesian Vector Autoregressions (VARs) are now routinely used with large
numbers of dependent variables. Use of non-conjugate priors or non-Gaussian
error distributions typically requires the use of Markov Chain Monte Carlo
(MCMC) methods, which leads to a large computational burden. This means full
system estimation of the reduced form VAR is difficult or infeasible with
large VARs. This has led many researchers to avoid full system estimation and
instead work with a structural VAR with a diagonal error covariance matrix.
The structural VAR allows for estimating one equation at a time which greatly
reduces the computational burden making Bayesian estimation of large VARs
practical. However, standard specifications for the structural VAR which allow
for equation-by-equation estimation (e.g. Carriero et al. (2019)) suffer from
order dependence (i.e. posterior and predictive densities depend on the manner
in which the variables are ordered in the VAR). The importance of order
dependence, and in particular its impact on predictive variances in larger
VARs, is discussed in papers such as Arias et al. (forthcoming) and Chan et
al. (2021). There have been some order invariant approaches proposed which do
allow for equation-by-equation estimation, including Chan et al. (2021) and Wu
and Koop (2023) but these assume Gaussian errors and the former relies on the
presence of stochastic volatility to identify the model. However, the presence
of large VAR shocks which imply sudden shifts in variances and/or asymmetries
in predictive densities means more flexibility is required. These
considerations motivate the present paper where we develop a VAR with a
nonparametric non-Gaussian error distribution for the shocks. The MCMC
algorithm we derive for this VAR is computationally efficient and order
invariant.

There is a growing VAR literature which wishes to develop flexible models for
the VAR errors. However, our model differs from this literature by its use of
non-parametric methods and its focus on computational efficiency and order
invariance. For instance, several papers, including Chiu et al. (2017) and
Karlsson et al. (2023) work with VARs with parametric, non-Gaussian error
distributions (e.g. Student-t distributions), but do not allow for equation-
by-equation estimation and, thus, are computationally slower than our
approach. A VAR with a non-parametric distribution for the shocks is Braun
(2021) which does allow for equation-by-equation estimation. But the flexible
error distribution in this case is used to identify structural shocks and does
not address the ordering issue. And, in addition, the MCMC algorithm becomes
more involved due to the fact that the contemporaneous relations need to be
sampled using Metropolis-Hastings (MH) updates.

To develop our model, we borrow ideas from the literature on semi- and non-
parametric estimation of random effects in panel data models, see Frühwirth-
Schnatter et al. (2004) or Dunson and Xing (2009). The key insight is that we
can exploit the random effect representation of the covariance matrix of the
system to enable equation-by-equation estimation, see Fox and Dunson (2015).
We decompose the shock vector of the VAR into two components. The first
component is a vector of random effects that feature an unknown multivariate
distribution which exhibits correlation between the errors in different
equations. The second component is a vector of Gaussian random shocks which
are uncorrelated across equations. Conditional on the random effects, the
model becomes a system of uncorrelated regression models. But after
integrating out the random effects, the resulting shock distribution features
cross-sectional dependence. The key implication is that fast estimation is
possible since the VAR coefficients can be drawn one equation at a time
conditional on the random effects. By assuming a Dirichlet Process Mixture
(DPM) for the vector of random effects we achieve great flexibility since the
joint distribution of the shocks can be skewed, feature heavy tails, be
heteroskedastic or be multi-modal. In the multiple equation VAR context, this
flexibility is potentially of great benefit since it allows for the errors in
the different equations to have different properties.

In an exercise using artificial data, we show how our non-parametric VAR can
automatically uncover a variety of departures from Gaussianity. In an
empirical exercise involving a large data set of US macroeconomic variables,
we demonstrate the advantages of our model both for forecasting and for
structural economic analysis.

The remainder of the paper is structured as follows. The next section shows
how a random effects representation of the VAR can be used to facilitate
equation-by-equation estimation and how we treat non-parametric shocks in the
VAR. Section 3 discusses the prior setup, sketches the MCMC algorithm and
discusses some computational details. Sections 4 and 5.4 apply the model to
synthetic and real data, respectively. The final section gives a summary and
concludes the paper.

##  2 VARs with non-parametric shocks

###  2.1 A linear VAR with an additive Gaussian error structure

Before introducing our non-parametric specification for the shocks to the VAR,
it is instructive to begin with a parametric version of our model. This
involves an $M$-dimensional vector of dependent variables,
$\\{\bm{y}_{t}\\}_{t=1}^{T}$, which evolves as:

$\bm{y}_{t}=\bm{A}\bm{X}_{t}+\bm{\epsilon}_{t}+\bm{v}_{t},\quad\bm{\epsilon}_{t}\sim\mathcal{N}(\bm{0}_{M},\bm{\Sigma})\quad\bm{v}_{t}\sim\mathcal{N}(\bm{0}_{M},\bm{\Omega}_{t})$
(1)

where
$\bm{X}_{t}=(\bm{y}^{\prime}_{t-1},\dots,\bm{y}^{\prime}_{t-p})^{\prime}$ is a
$K(=Mp)$ matrix of lagged endogenous variables and $\bm{A}$ denotes an
$M\times K$ matrix of VAR coefficients. This specification differs from a
conventional VAR in that it has two $M$-dimensional errors,
$\bm{\epsilon}_{t}$ and $\bm{v}_{t}$, which are assumed to be independent over
time and of one another at all leads and lags. The only restriction on
$\bm{\Sigma}$ is that it is positive definite whereas $\bm{\Omega}_{t}$ is
restricted to be a diagonal matrix with individual error variances
$\omega_{1t},\dots,\omega_{Mt}$. We assume that the logarithms of
$\omega_{jt}$ evolve according to AR(1) processes leading to a standard SV
specification. The covariance matrix of the VAR errors,
$\bm{\varepsilon}_{t}=\bm{\epsilon}_{t}+\bm{v}_{t}$ is
$\bm{\Xi}_{t}=\bm{\Sigma}+\bm{\Omega}_{t}$. Notice that the SV assumption on
the idiosyncratic shocks implies that the main diagonal elements of
$\bm{\Xi}_{t}$ are given by $\sigma_{ii}^{2}+\omega_{it}$, where
$\sigma_{ii}^{2}$ is the $(i,i)^{th}$ element of $\bm{\Sigma}$.111We note
that, without further restrictions, this specification is not identified. But
this lack of identification poses no problem for Bayesian estimation and
prediction with this model provided a proper prior is used. Furthermore, if
the researcher wishes to impose identification, this can easily be done in
standard cases. For instance, beginning with Cogley and Sargent (2005) many
popular Bayesian VARs have assumed the $\omega_{it}$ to follow stochastic
volatility (SV) processes. Fixing the SV initial conditions to be zero
suffices to identify the model.

The reasons why we adopt this additive error structure relate to efficient
computation and order invariance. To explain these points, we first summarise
some key issues in the Bayesian VAR literature which are particularly acute
with larger VARs. The traditional reduced form VAR is given by (1) but with
$\bm{v}_{t}$ set to zero.222Typically, the assumption of homoskedasticity is
relaxed, but allowing for this does not affect the arguments in this sub-
section. Direct Bayesian estimation of this reduced form VAR is
computationally challenging when $M$ is large due to the need to carry out
matrix manipulations involving the very high dimensional posterior covariance
matrix of $\bm{A}$. Accordingly, it is common to carry out Bayesian estimation
in a structural form involving the use of the Cholesky decomposition of the
reduced form error covariance matrix. That is, decomposing
$\bm{\Sigma}=\bm{B}\bm{D}\bm{B}^{{}^{\prime}}$ where $\bm{B}$ is lower
triangular with ones on the diagonal and $\bm{D}$ is diagonal, the structural
VAR is obtained by multiplying both sides of the VAR by $\bm{B}^{-1}$. The
structural VAR has a diagonal error covariance matrix which means estimation
can be carried out one equation at a time. This leads to huge computational
improvements. For instance, in the specification used in Carriero et al.
(2019) the MCMC algorithm based on the reduced form VAR requires $O(M^{6})$
elementary operations to take one draw of the VAR coefficients but only
$O(M^{4})$ with the structural VAR. Thus, there are enormous computational
benefits from working with VARs specified in such a way as to allow equation-
by-equation estimation.

However, Bayesian results using the structural VAR based on the Cholesky
decomposition are order dependent (i.e. posterior and predictive results
depend on the way the variables are ordered in the VAR). This contrasts with
the reduced form VAR for which standard implementations (e.g. use of an
inverse Wishart prior for $\bm{\Sigma}$) are order invariant. The empirical
importance of ordering issues has been investigated in papers such as Arias et
al. (forthcoming) and Chan et al. (2021) and found to be substantive,
particularly in the case of large VARs and particularly for predictive
variances and higher order predictive moments. Thus, most Bayesian VAR papers
either work in reduced form, and face computational challenges unless the VAR
dimension is very low, or work in structural form and produce empirical
results which depend on the way the variables are ordered in the VAR.

If we now return to our VAR with additive errors in (1), it is straightforward
to show that it suffers neither of these drawbacks. Computationally efficient
MCMC algorithms can be developed which exploit the fact that (conditional on
$\bm{\epsilon}_{t}$) the equations are independent of one another. But
marginally (i.e. after integrating out $\bm{\epsilon}_{t}$), the shocks to
$\bm{y}_{t}$ are cross-sectionally correlated. In relation to order-
invariance, it is worth emphasizing that the ordering issue does not relate to
the likelihood (i.e. the structural and reduced form VARs lead to the same
likelihood function) but rather relates to the prior that is placed on the
error variances and covariances. Indeed Carriero et al. (2019) refer to it as
the ”prior ordering issue”. In our additive error setup, order invariance can
be achieved by retaining an inverse Wishart prior for $\bm{\Sigma}$. Since
$\bm{\Omega}_{t}$ is diagonal any conventional set of priors will lead to
order invariance. For instance, in the homoskedastic version of the model,
assuming $\omega_{i}$ for $i=1\ldots M$ to have inverse-Gamma priors which are
independent across $i$ leads to order invariance. In the heteroskedastic case,
assuming $\omega_{it}$ to have independent SV processes leads to order
invariance.

In subsequent sections, we will work with a non-parametric version of this
model and provide full details of the priors we use and our computationally
efficient MCMC algorithm. We stress that the issues discussed in this sub-
section also hold with non-parametric VARs. However, this subsection provides
the basic insights into how these computational benefits are achieved and why
our prior is order invariant. In addition, it may be found useful by Bayesian
VAR researchers who are happy to remain parametric and work with linear VARs
with Gaussian shocks. Specifying the VAR as we have done, with additive
errors, is an attractive way of achieving fast, order invariant inference even
in conventional VARs.

###  2.2 Allowing for shocks of unknown form

The model in the preceding subsection had a linear conditional mean and
Gaussian error structure. In many contexts, linearity and Gaussianity can be
restrictive and this is particularly so in extreme times such as the recent
Covid-19 pandemic. In this paper we focus on relaxing the Gaussianity
assumption relating to the shocks hitting the model. We will maintain the
assumption of a linear conditional mean and assume a standard SV process for
$\bm{\Omega}_{t}$.333Adopting non-parametric approaches for either of these
can easily be done. For instance, Huber and Rossini (2020), Huber et al.
(2020) and Clark et al. (forthcoming) model the conditional mean of a VAR non-
parametrically using regression trees. The last of these papers also uses
regression trees to model the conditional variance. Approaches such as these
could be added to the model of the present paper if extra flexibility is
desired. However, as will be demonstrated below our model is already very
flexible and can model any of the empirical regularities common with
macroeconomic data. In particular our non-parametric VAR is the same as the
one specified in the preceding section except for assumptions relating to
$\bm{\epsilon}_{t}$. The assumption of Gaussianity will be replaced by a
Dirichlet process mixture (DPM) of Gaussians. We will show that by this simple
extension, we will achieve great gains in empirically-relevant flexibility.
But because the error process remains conditionally Gaussian, the benefits
discussed in the preceding sub-section (i.e. equation-by-equation estimation
and order invariance) will be retained.

The ideas underlying our treatment of $\bm{\epsilon}_{t}$ are inspired by
papers such as Frühwirth-Schnatter et al. (2004) and Dunson and Xing (2009)
which develop parametric and non-parametric Bayesian treatments of random
effects in panel data models and, accordingly, we refer to $\bm{\epsilon}_{t}$
as a vector of random effects. We model the random effects by introducing a
base measure $\mathcal{G}$ and defining a parametric family of component
densities $f$ with unknown parameters $\bm{\vartheta}$:

$p(\bm{\epsilon}_{t})=\int
f(\bm{\epsilon}_{t}|\bm{\vartheta})~{}\mathcal{G}(d\bm{\vartheta})=\sum_{j=1}^{\infty}\eta_{j}f(\bm{\epsilon}_{t}|\bm{\vartheta}_{j}),$

where the weights $\sum_{j=1}^{\infty}\eta_{j}$ sum to $1$. We assume that the
component densities $f$ are Gaussian with $M\times 1$ mean vector
$\bm{\mu}_{j}$ and $M\times M$ variance-covariance matrix $\bm{\Sigma}_{j}$
which implies

$p(\bm{\epsilon}_{t})=\sum_{j=1}^{\infty}\eta_{j}\mathcal{N}(\bm{\mu}_{j},\bm{\Sigma}_{j}).$

Alternatively, this model can be written in terms of a discrete latent random
variable $\delta_{t}\in\\{1,2\dots\\}$ that indicates which mixture component
to adopt in time $t$:

$\bm{\epsilon}_{t}=\bm{\mu}_{\delta_{t}}+\bm{Q}_{\delta_{t}}\bm{w}_{t},\quad\bm{w}_{t}\sim\mathcal{N}(\bm{0},\bm{I}_{M})$
(2)

with $\text{Prob}(\delta_{t}=j)=\eta_{j}$ and $\bm{Q}_{\delta_{t}}$ denoting
the lower Cholesky factor of $\bm{\Sigma}_{\delta_{t}}$ with
$\bm{\Sigma}_{\delta_{t}}=\Sigma_{j}$ if $\delta_{t}=j$.444Note that this use
of the Cholesky decomposition does not undermine the order invariance of the
model since we will use an inverse Wishart prior for $\bm{\Sigma}_{j}$ which
is order-invariant. Plugging this expression for $\bm{\epsilon}_{t}$ into (1)
gives

$\bm{y}_{t}=\bm{\mu}_{\delta_{t}}+\bm{A}\bm{X}_{t}+\underbrace{\bm{Q}_{\delta_{t}}\bm{w}_{t}+\bm{v}_{t}}_{\bm{\varepsilon}_{t}}$
(3)

which is VAR with reduced-form shocks
$\bm{\varepsilon}_{t}\sim\mathcal{N}(\bm{0}_{M},\bm{\Xi}_{t})$ and time-
varying error covariance matrix
$\bm{\Xi}_{t}=\bm{\Sigma}_{\delta_{t}}+\bm{\Omega}_{t}$. Another interesting
interpretation of $\bm{\varepsilon}_{t}$ is that it resembles a factor model
with $M$ factors and a particular form of the factor loadings.

The model given in (3) has several properties which make it attractive for use
in empirical macroeconomics. It retains a conditionally (i.e. conditional on
$\delta_{t}$) Gaussian structure which leads to simplicity of computation and
structural economic interpretation. However, it is extremely flexible since
infinite mixtures of Gaussians can approximate any distribution. Notice that
$\bm{\Xi}_{t}$ varies across components in the DPM. Depending on the estimated
values for $\delta_{t}$ this allows for a wide variety of behavior (i.e.
structural breaks, regime switching, outliers, etc.) in the contemporaneous
relationships between the elements in $\bm{y}_{t}$.

By additionally allowing for SV (via our specification for $\bm{\Omega}_{t}$)
we have an error process of great flexibility allowing for impulse responses
and other structural features to differ over time in a way that is estimated
from the data. That is, both $\delta_{t}$ and $\bm{\Omega}_{t}$ allow for
different types of parameter change, the latter only changing the error
variances and being of a smooth nature, with the former additionally relating
to the error covariances and allowing for more abrupt types of regime change
or structural break. To see this feature, notice that a typical main diagonal
element of $\bm{\Xi}_{t}$ is:

$[\bm{\Xi}_{t}]_{ii}=\sigma^{2}_{ii,\delta_{t}}+\omega_{it}.$

In this equation, $\sigma^{2}_{ii,\delta_{t}}$, the $(i,i)^{th}$ element of
$\bm{\Sigma}_{\delta_{t}}$, changes abruptly over time and is thus capable of
handling large outliers whereas $\omega_{it}$ changes smoothly and thus
captures slowly varying trends in the error variances. Using only the latter
implies that in the presence of large shocks, the SV model would only slowly
adapt and would thus imply a higher variance when the large shock has already
faded out. This model resembles the SV with outliers model of Carriero et al.
(2022b).

The conditional representation of the model also gives insights on how the DPM
handles location shifts in the shocks. To see this, note that (3) allows for
the intercept to change over time in a non-parametric manner. Models with
time-varying intercepts are common in macroeconomics (see, e.g., Stock and
Watson (2007) and Antolin-Diaz et al. (2017)). Traditionally, intercepts have
assumed to follow a random walk. However, our non-parametric treatment allows
us to uncover the form of parameter change from the data (see also
Hauzenberger et al., 2022, for a related but parametric treatment of time-
varying parameter regressions).

It is also worth noting that an alternative model for non-parametric shocks
would omit the additive error structure by setting $\bm{v}_{t}=0$ and simply
have one vector of errors which is modelled using a DPM. But this apparently
simpler form would both be more computationally demanding (i.e. since order
invariant equation-by-equation estimation would be difficult to achieve) and
would omit the SV process. That is, a DPM model for the errors on its own
would be very flexible at modelling structural breaks and outliers, but the
assumption that the DPM errors are independent over time means that is less
able to model gradual changes in volatility. In contrast our model combines
the benefits of a very flexible shock distribution with the gradual volatility
change of an SV process.

##  3 Bayesian inference in the VAR with non-parametric shocks

###  3.1 Prior

In this sub-section we describe our prior. We emphasize that the innovations
in this paper relate to the parameters in the random effects. For the
remaining parameters, any standard Bayesian priors can be used. In this paper,
we use the Normal-Gamma prior of Brown and Griffin (2010) for the VAR
coefficients (see, e.g., Huber and Feldkircher, 2019) although any other
common Bayesian VAR prior could be used (e.g. the Minnesota prior or a global
local shrinkage prior such as the Horseshoe). We also adopt standard priors
for the parameters of the log-volatilities $\log\omega_{it}$. These are
Gaussian priors on the unconditional mean $\mu_{\omega,j}$, (transformed) Beta
priors on the persistence parameter $\rho_{\omega,j}$ and Gamma priors on the
innovation variances $\sigma^{2}_{\omega,i}$, respectively.555More precisely,
we set $\mu_{\omega,j}\sim\mathcal{N}(0,10)$,
$\frac{\rho_{\omega,j}+1}{2}\sim\mathcal{B}(25,5)$ and
$\sigma^{2}_{\omega,i}\sim\mathcal{G}(1/2,1/2)$. We also consider a
homoskedastic version of our model where $\bm{\Omega}_{t}=\bm{\Omega}$ for all
$t$ and in this case assume priors
$\omega_{i}\sim\mathcal{G}^{-1}(a_{\omega},b_{\omega})$ for $i=1\ldots M$ with
$a_{\omega}=b_{\omega}=10^{-3}$.

Our prior for the mean of the random effects is similar to one developed in
Malsiner-Walli et al. (2017). It is a Gaussian prior on $\bm{\mu}_{j}$ that
shrinks the different elements of $\bm{\mu}_{j}$ towards a common location:

$\bm{\mu}_{j}\sim\mathcal{N}(\bm{\mu}_{0},\bm{B}_{0})\text{ for
}j=1,\dots,\infty,$

with $\bm{B}_{0}=\text{diag}(b_{1},\dots,b_{M})$ being a diagonal prior
variance matrix with

$b_{j}\sim\mathcal{G}(c_{b},d_{b}).$

The hyperparameters $c_{b},d_{b}$ are greater than zero. This is the Normal-
Gamma prior proposed in Brown and Griffin (2010). If $c_{b}=1$ we obtain the
LASSO (Park and Casella, 2008). However, the LASSO is known to overshrink
significant signals and undershrink irrelevant ones. Hence, we set
$c_{b}=d_{b}=0.6$. This leads to a model that implies more shrinkage and
flexible tail behavior. Since our prior is fully hierarchical, we also require
another prior on $\bm{\mu}_{0}$. This is assumed to be
$\mathcal{N}(\bm{0},c^{-1}\bm{I}_{M})$ with $c\to 0$, yielding a non-
informative prior. In all our empirical work, we set $c=10^{-3}$ to render the
prior relatively non-informative but proper.

The combination of a flexible shrinkage prior that forces the component-
specific means towards a common location has implications on the clustering
behavior of the mixture model. To illustrate this, let $\mu_{ji}$ denote the
$j^{th}$ element of $\bm{\mu}_{i}$ for $i=k$ or $\tilde{k}$ (i.e. these are
the two intercepts in the $j^{th}$ equation for two different components of
the infinite mixture $k\neq\tilde{k}$). The prior above implies the following
in terms of the distance between $\mu_{jk}$ and $\mu_{j\tilde{k}}$ (see Yau
and Holmes, 2011):

$\frac{(\mu_{jk}-\mu_{j\tilde{k}})}{\sqrt{2}}\sim\mathcal{N}(0,b_{j}).$

Thus, our prior is centred over intercept homogeneity and $b_{j}$ controls the
strength of this belief. If $b_{j}$ is close to zero, the intercepts collapse
to a common value (which is $\mu_{0j}$, the $j^{th}$ element of
$\bm{\mu}_{0}$). For larger $b_{j}$, we allow for more heterogeneity in the
intercepts. This feature is crucial since the presence of the non-zero
location parameter allows us to capture skewness in the shocks. If we use a
non-informative prior on the component means we would risk overfitting the
data. Our shrinkage prior effectively enables us to investigate how much
asymmetries are in the data in a fully automatic manner.

For the covariance matrices for each component in the DPM, we use a conjugate
Wishart prior on $\bm{\Sigma}^{-1}$ as this leads to order invariance within
each component which implies order invariance in the VAR as a whole. Thus we
assume

$\bm{\Sigma}^{-1}_{k}\sim\mathcal{W}(c_{0},\bm{\Sigma}_{0}^{-1}).$

Note that we parameterize the Wishart such that the prior mean equals
$c_{0}\bm{\Sigma}_{0}^{-1}$ with $c_{0}$ being its degrees of freedom. The
prior hyperparameters can be chosen in any way. In our empirical work, we use
a relatively noninformative data-based prior inspired by the Minnesota prior.
In particular, we set the prior degrees of freedom as $c_{0}=2(2.5+(M-1)/2)$.
The prior scaling matrix is estimated from the data and set as
$\bm{\Sigma}_{0}=\text{diag}(\hat{\sigma}^{2}_{1},\dots,\hat{\sigma}^{2}_{M})$
where $\hat{\sigma}^{2}_{j}$ are the OLS error variances obtained by running
an AR($p$) model for $y_{jt}$.

In terms of the weights in the DPM, we use a stick breaking process (SBP)
prior on $\eta_{j}$. The SBP prior introduces additional auxiliary random
variables $\nu_{j}$ (called sticks) such that the weights are obtained
sequentially:

$\eta_{1}=\nu_{1},\quad\eta_{j}=\nu_{j}\prod_{i=1}^{j-1}(1-\nu_{j}),\quad\nu_{j}\sim\mathcal{B}(1,\alpha).$

The parameter $\alpha$ determines the clustering behavior of the mixture
model. To see this, notice that the prior probability of forming a new cluster
when assigning $\bm{\varepsilon}_{t}$ conditional on all
$\bm{\varepsilon}_{\tau}(\tau\neq t)$ is (Lau and Green, 2007; Frühwirth-
Schnatter and Malsiner-Walli, 2019):

$\frac{\alpha}{T-1+\alpha},$ (4)

and thus decreases in $T$. This implies that the DPM has the potential to
create few larger clusters and then has a low probability of opening up new
clusters that are populated by relatively few observations. In macroeconomic
data, this behavior might be necessary to single out events that are are
different from those produced by the DGP in normal times such as the Covid-19
pandemic. Since $\alpha$ crucially impacts this behavior, we treat it as an
unknown parameter and estimate it from the data. We assume that $\alpha$
arises from a Gamma distribution a priori, i.e., $\alpha\sim\mathcal{G}(2,4)$
which implies a prior mean $0.5$ and a prior variance $0.125$. This choice was
originally suggested by Escobar and West (1995) and encourages clustering
behavior of the mixture model.

###  3.2 Posterior simulation of the VAR coefficients and random effects

In this sub-section we describe how we simulate the VAR coefficients and
random effects in more detail. The other steps are relatively standard and we
sketch them in the next sub-section. A key theme of this paper is
computational efficiency and, to achieve this end, we need to draw the VAR
coefficients one equation at a time. This requires knowledge of the random
effects which serve to establish correlations across the shocks. Accordingly,
we describe how we draw the VAR coefficients and the random effects in more
detail.

Conditional on the random effects $\\{\bm{w}_{t}\\}_{t=1}^{T}$ and the
estimates of $\\{\bm{Q}_{\delta_{t}}\\}_{t=1}^{T}$ we draw the equation-
specific VAR coefficients $\bm{A}_{i}$ from

$\bm{A}_{i}|\bullet\sim\mathcal{N}(\bm{m}_{i},\bm{V}_{a,i}),\quad
i=1,\dots,M,$

where $\bullet$ denotes all arguments necessary to define the full conditional
posterior distribution and:

$\displaystyle\bm{V}_{a,i}$
$\displaystyle=(\tilde{\bm{X}}^{\prime}_{i}\tilde{\bm{X}}_{i}+\underline{\bm{V}}^{-1}_{a,i})^{-1},$
$\displaystyle\bm{m}_{i}$
$\displaystyle=\bm{V}_{a,i}\tilde{\bm{X}}^{\prime}_{i}\tilde{\bm{Y}}_{i},$

denote the posterior moments with $\tilde{\bm{X}}_{i}$ being a $T\times K$
matrix with typical $t^{th}$ row $\bm{X}_{t}/\sqrt{\omega}_{it}$,
$\underline{\bm{V}}_{a,i}$ denotes the prior variance matrix,
$\tilde{\bm{Y}}_{i}$ denotes a $T\times 1$ vector with typical element given
by $(y_{it}-\bm{q}^{\prime}_{i\delta_{t}}\bm{w}_{t})/\sqrt{\omega}_{it}$ and
$\bm{q}_{i\bullet,\delta_{t}}$ denotes the $i^{th}$ row of
$\bm{Q}_{\delta_{t}}$. For large models, the inversion of the posterior
covariance (and computation of its Cholesky factor) becomes computationally
cumbersome. However, in these cases (characterized by $T\ll K$) we can use the
algorithm outlined in Bhattacharya et al. (2016) that speeds up computation
enormously.

Next we describe the sampling steps involved in simulating the random effects
in more detail. As stated in the previous sub-section, the model in (3)
resembles a factor model with the particular structure on the factor loadings
given in (2). We can use this representation to back out the conditional
posterior distribution of $\bm{w}_{t}$, $p(\bm{w}_{t}|\bullet)$. The key point
to notice is that the random effects are conditionally independent over time
and hence
$p(\bm{w}_{1},\dots,\bm{w}_{T}|\bullet)=\prod_{t=1}^{T}p(\bm{w}_{t}|\bullet)$
with time $t$ posteriors given by:

$p(\bm{w}_{t}|\bullet)=\mathcal{N}(\overline{\bm{w}}_{t},\overline{\bm{V}}_{w,t}),$
(5)

and moments given by

$\displaystyle\overline{\bm{V}}_{w,t}$
$\displaystyle=(\tilde{\bm{Q}}_{\delta_{t}}\tilde{\bm{Q}}^{\prime}_{\delta_{t}}+\bm{I}_{M})^{-1},$
$\displaystyle\overline{\bm{w}}_{t}$
$\displaystyle=\overline{\bm{V}}_{w,t}(\tilde{\bm{Q}^{\prime}}_{\delta_{t}}\tilde{\bm{\varepsilon}_{t}}),$

where $\tilde{\bm{Q}}_{\delta_{t}}=\bm{Q}_{\delta_{t}}\bm{\Omega}_{t}^{-1/2}$
and $\tilde{\bm{\varepsilon}_{t}}=\bm{\Omega}_{t}^{-1/2}\bm{\varepsilon}_{t}$.

###  3.3 Full conditional MCMC sampling

Our MCMC sampler is relatively straightforward and contains several steps
which involve standard full conditional distributions. Here, we summarize the
different updating steps necessary to sample from the joint distribution of
the latent states and coefficients of the model.

Our sampler iterates between the following steps:

  1. 1.

Sample $\bm{A}_{i}|\bullet$ for each equation from
$p(\bm{A}_{i}|\bullet)=\mathcal{N}(\bm{m}_{i},\bm{V}_{a,i})$ as described in
the previous sub-section.

  2. 2.

Sample $\bm{\Sigma}^{-1}_{k}$ from
$p(\bm{\Sigma}_{k}|\bullet)=\mathcal{W}(\overline{c}_{k},\overline{\bm{\Sigma}}_{k})$.
The posterior degrees of freedom are $\overline{c}_{k}=c_{0}+T_{k}/2$ where
$T_{k}=\sum_{t=1}^{T}\mathbb{I}(\delta_{t}=k)$ denotes the number of
observations allocated to cluster $k$. The posterior scaling matrix is given
by

$\overline{\bm{\Sigma}}_{k}=\frac{1}{2}{\sum_{t:\delta_{t}=k}(\tilde{\bm{y}}_{t}-\bm{\mu}_{\delta_{t}})(\tilde{\bm{y}}_{t}-\bm{\mu}_{\delta_{t}})^{\prime}}+\bm{\Sigma}^{-1}_{0},$

where $\tilde{\bm{y}}_{t}=\bm{y}_{t}-\bm{A}^{\prime}\bm{x}_{t}-\bm{v}_{t}$.

  3. 3.

Sample $\bm{\mu}_{k}$ from
$p(\bm{\mu}_{k}|\bullet)=\mathcal{N}(\overline{\bm{\mu}}_{k},\overline{\bm{V}}_{\mu,k})$.
The posterior moments of the random intercept terms are given by

$\displaystyle\overline{\bm{V}}_{\mu,k}$
$\displaystyle=\left(\bm{\Sigma}_{k}^{-1}T_{k}+\bm{B}_{0}^{-1}\right)^{-1},$
$\displaystyle\overline{\bm{\mu}}_{k}$ $\displaystyle=\overline{\bm{V}}_{\mu
k}\left(\bm{\Sigma}_{k}^{-1}\sum_{t:\delta_{t}=k}\tilde{\bm{y}}_{t}+\bm{B}_{0}^{-1}\bm{\mu}_{0}\right).$

  4. 4.

Sample $\bm{\mu}_{0}$ from
$p(\bm{\mu}_{0}|\bullet)=\mathcal{N}(\overline{\bm{\mu}}_{0},\overline{\bm{V}}_{\mu,0})$
with

$\displaystyle\overline{\bm{V}}_{\mu,0}=\frac{1}{J}\bm{B}_{0},\quad\overline{\bm{\mu}}_{0}=\frac{1}{J}\sum_{j=1}^{J}\bm{\mu}_{j}.$

  5. 5.

Sample $\bm{w}_{t}$ from
$p(\bm{w}_{1},\dots,\bm{w}_{T}|\bullet)=\prod_{t=1}^{T}\mathcal{N}(\overline{\bm{w}}_{t},\overline{V}_{wt})$
with moments given below Eq. (missing)5.

  6. 6.

Sample $b_{i}$ (for $i=1,\dots,M$) from
$p(b_{j}|\bullet)=\text{GIG}(p_{0},p_{K},z_{j})$ with $p_{0}=2d_{b}$,
$p_{K}=c_{b}-J/2$ and $z_{j}=\sum_{j=1}^{J}(\mu_{j,i}-\mu_{0,i})^{2}$ and GIG
denoting the generalized inverse Gaussian distribution.666The density of the
GIG is given by $f(x)=x^{\lambda-1}e^{-\frac{1}{2}(\chi/x+\psi x)}$.

  7. 7.

Sample $\omega_{1t},\dots,\omega_{Mt}$ from
$p(\omega_{1t},\dots,\omega_{Mt}|\bullet)$. For the SV case,
$\omega_{j1},\dots,\omega_{jT}$ (for $j=1,\dots,M$) and the parameters of the
state equation are drawn using the algorithm outlined in Kastner and
Frühwirth-Schnatter (2014) and implemented in the R package stochvol. In case
we assume homoskedasticity, we simply draw $\omega_{j}$ from an inverse Gamma
posterior.

  8. 8.

Sample $\nu_{1},\dots,\nu_{J-1}$ from
$p(\nu_{1},\dots,\nu_{J-1}|\bullet)=\prod_{j=1}^{J}\mathcal{B}\left(1+T_{j},\alpha+\sum_{l=j+1}^{J}T_{l}\right)$.

  9. 9.

Sample $\delta_{1},\dots,\delta_{T}$ from
$p(\delta_{1},\dots,\delta_{T}|\bullet)=\prod_{t=1}^{T}p(\delta_{t}|\bullet)$
using the Slice sampler (Kalli et al., 2011) in two steps. First, let
$u_{t}|\delta_{t}\sim\mathcal{U}(0,\zeta_{\delta_{t}})$ denote a set of
auxiliary random variables with $\zeta_{k}=(1-w)w^{k-1}$ and $w=0.8$. Then,
conditional on $u_{t}$, we simulate $\delta_{t}$ from its discrete
distribution as follows:

$\text{Prob}(\delta_{t}=k|\bullet)\propto\frac{\mathbb{I}(u_{t}<\xi_{k})}{\xi_{k}}~{}\eta_{k}~{}f_{\mathcal{N}}(\tilde{\bm{y}}_{t}|\bm{\mu}_{k},\bm{\Sigma}_{k}).$

  10. 10.

Sample $\alpha$ using a Metropolis-Hastings step.

Notice that $J$ is a truncation parameter that determines the effective number
of regimes which is obtained by solving
$1-\sum_{j=1}^{J}\zeta_{j}<\text{min}(u_{1},\dots,u_{T})$. This implies that
our infinite mixture model becomes effectively finite dimensional and thus
computationally tractable.

In most our applications we repeat this algorithm 20,000 times and discard the
first 10,000 draws as burn-in. For the actual data application and in
simulations we did not encounter any mixing issues if we consider functions of
the parameters such as impulse responses or forecast densities.

Before applying our model to the data, it is worth stressing that the mixture
model is not identified with respect to relabeling the latent discrete
indicators. This does not cause any issues if interest is exclusively on
either unconditional impulse responses or forecast distributions. In our
structural application, we discuss impulse responses conditional on a given
cluster. To avoid label switching in this case, we introduce restrictions that
ensure a unique ordering ex-post. More details are provided below.

###  3.4 Computational aspects of our algorithm

The algorithm described in the previous sub-section is efficient and has the
same complexity as the original (but incorrect) version of the algorithm
proposed in Carriero et al. (2019) and its corrected variant proposed in
Carriero et al. (2022a).777Notice that Carriero et al. (2019) do not use the
efficient sampling algorithm of Bhattacharya et al. (2016) to simulate from
the equation-specific coefficient posteriors if $T<K$.  This is because we
sample from the equation-specific posteriors. This implies that the posterior
covariance matrices are $K\times K$ and inversion of such matrices has
computational complexity $O(M^{4})$. A similar computational advantage can be
obtained by applying the algorithm outlined in Kastner and Huber (2021) which
uses a factor model to render equation-by-equation estimation possible. In any
case, all these algorithms are much more efficient than the one based on
treating the VAR as a full system of equations in large dimensions.

Figure 1 illustrates the magnitude of efficiency improvements of equation-by-
equation estimation relative to full system estimation and the increase in
computation required to add the DPM. It compares our model (BVAR-DPM) to a
restricted version of it which is Gaussian and obtained by setting $G=1$
(BVAR-G1) and a BVAR estimated under a non-conjugate independent inverse
Wishart prior (BVAR-NIW). All models feature a single lag. Figure 1 compares
the computation times necessary to generate $10$ draws from the posterior.
Note that the figure is in log-scales.

Figure 1: Comparison of computation times between the BVAR-DPM, BVAR-G1 and
the BVAR with an independent Normal-Wishart prior (BVAR-NIW) for
$M\in\\{5,\dots,40\\}$. The y-axis is in log scale.

The figure shows that BVAR-DPM and BVAR-G1 are much faster than full-system
estimation when $M$ exceeds $8$ (in the case of BVAR-G1) and when $M$ exceeds
$12$ (in the case of the BVAR-DPM). In larger panels, the BVAR-NIW becomes
much slower and the computational burden increases substantially. Differences
between the blue and orange line reflect the additional computational burden
from adding the DPM piece to the model. In most empirical work, the effective
number of regimes is small (i.e. the infinite mixture reduces to a four-
components mixture of Gaussians). This implies that if the true number of
regimes is small, the computational burden would also decrease under the DPM
specification.

The algorithms of Carriero et al. (2019) and Carriero et al. (2022a) have the
same computational complexity as our approach that sets $G=1$. In this case,
computation times would be very similar (and the shape of the curves would, in
fact, be identical) and our algorithm thus scales as well as theirs. However,
it is worth stressing that ours is order invariant, can have non-parametric
errors, and is applicable to any conditional mean function.

##  4 Artificial data exercise

We illustrate the merits of our approach by simulating data from a set of
different DGPs. These DGPs differ both in terms of model size and the error
distributions. With respect to model size, we consider three sizes that
capture typical situations in applied macroeconomic work. The smallest DGP has
$M=5$ endogenous variables, the medium-sized DGP features $M=10$ and the
largest includes $M=20$ endogenous variables.

With respect to the error distributions we consider four different shock
assumptions. Three of these four feature substantial departures from
homoskedasticity and normality while one assumes Gaussian and homoskedastic
shocks. All of these DGPs assume that the conditional mean of the process is
given by:

$\bm{y}_{t}=\bm{A}\bm{y}_{t-1}+\bm{\epsilon}_{t},\quad\bm{y}_{0}=\bm{0},\quad
t=1,\dots,250,$

where $\bm{A}$ has diagonal elements $A_{ii}=0.75$ and off-diagonal elements
sampled from a Gaussian distribution $A_{ij}\sim\mathcal{N}(0,0.1^{2})(i\neq
j)$. To ensure stationarity we reject draws of $\bm{A}$ which imply an
unstable model.

The errors $\bm{\epsilon}_{t}$ differ across the DGPs in the following four
ways.

  1. 1.

The first DGP assumes that the shocks follow a multivariate skew Gaussian
distribution: $\bm{\epsilon}_{t}\sim
s\mathcal{N}(\bm{\mu}_{t},\bm{W},\bm{\kappa})$ where $\bm{\mu}_{t}$ is a
location parameter that equals $\bm{0}$ with probability $98.5$ percent and
$-3\times z$, $z\sim\mathcal{U}(0,1)$ with probability $1.5$ percent of the
time. This increases the complexity of this error distribution further by
capturing a situation where we have a few, large and negative shocks. The
matrix $\bm{W}$ is created as follows. We specify a lower uni-triangular
matrix $\bm{U}$ with $u_{ij}\sim\mathcal{N}(0,0.1^{2})$ for
$i=2,\dots,M;j=1,\dots,M-1$ and then set $\bm{W}=\bm{U}\bm{U}^{\prime}$. The
parameter $\bm{\kappa}$ is obtained by simulating from a zero-mean Gaussian
distribution with variance $9$ and then rounded to the nearest integer.

  2. 2.

For the second DGP we assume that $\bm{\varepsilon}_{t}$ is multivariate
Student t with three degrees of freedom and covariance matrix $\bm{W}$.

  3. 3.

The third DGP has a common stochastic volatility specification (Carriero et
al., 2016). We assume that the shocks are Gaussian distributed with zero mean
and time-varying covariance matrix
$\bm{W}_{t}=e^{\mathfrak{s}_{t}}\times\bm{W}$. The log-volatility process
$\mathfrak{s}_{t}$ evolves according to a random walk with innovation variance
$0.25^{2}$

  4. 4.

Finally, we also consider a homoskedastic DGP that assumes that
$\bm{\varepsilon}_{t}\sim\mathcal{N}(\bm{0},\bm{W})$.

For comparison, we consider the BVAR-DPM and the BVAR-DPM with $G=1$ and set
the number of lags equal to 5. All models are estimated with homoskedastic and
SV measurement errors. This combination allows us to analyze whether the BVAR-
DPM is capable of recovering non-Gaussian features in the DGP without
overfitting. In particular, the final DGP can be used to focus on the question
whether adding the non-parametric component to the model leads to a
deterioration in estimation accuracy or whether the DPM can recover the simple
Gaussian case without overfitting.

Table 1: Simulation results | Skew | t | SV | Homosk.
---|---|---|---|---
$M\downarrow$ | DPM | $G=1$ | DPM | $G=1$ | DPM | $G=1$ | DPM | $G=1$
Homosk.
5 | 1.49 | 1.54 | 1.19 | 1.93 | 4.96 | 5.71 | 1.58 | 1.58
10 | 1.52 | 1.56 | 1.22 | 1.71 | 4.83 | 5.16 | 1.65 | 1.65
20 | 1.45 | 1.48 | 1.84 | 1.92 | 3.81 | 3.90 | 1.56 | 1.55
SV
5 | 1.86 | 1.99 | 1.36 | 1.57 | 1.93 | 1.80 | 1.94 | 1.92
10 | 1.78 | 1.80 | 1.23 | 1.35 | 1.71 | 1.68 | 1.82 | 1.81
20 | 1.59 | 1.60 | 1.21 | 1.22 | 1.57 | 1.56 | 1.67 | 1.67
Notes: This table shows mean absolute errors (MAEs) between the posterior
median of the VAR coefficients and the true VAR coefficients. DPM is the BVAR
with a DPM of Gaussians and $G=1$ is a BVAR with a set of random coefficients
(G=1).

To rank models in terms of estimation accuracy we focus on the mean absolute
error (MAE) between the true set of coefficients and the posterior median of
the estimated VAR coefficients. Table 1 shows the accuracy of the BVAR-DPM
versus the BVAR with Gaussian errors in terms of their estimation of the VAR
coefficients. All results are means across MAEs from 50 replications from each
of the DGPs.

In general, our results indicate that if the Gaussian-errored model is mis-
specified, the DPM model is consistently more accurate. But even for cases
where the BVAR $G=1$ is correctly specified (i.e. with the homoskedastic model
and homoskedastic DGP and with the SV model with the SV DGP) our BVAR-DPM is
only slightly less accurate. This illustrates how the BVAR-DPM can
successfully model features such as skewness and fat tails, but when the DGP
does not have such features it can successfully uncover the underlying
Gaussian model with very little over-fitting. When we consider differences
across model sizes, we find that accuracy gains from setting $G=\infty$ vis-á-
vis $G=1$ decline with larger information sets. This finding points towards
the fact that large models can soak up non-Gaussian features in the data.

Table 2: Effective number of clusters | Skew | t | SV | Homosk.
---|---|---|---|---
Homosk.
5 | 2.03 | 3.98 | 2.22 | 1.00
10 | 2.05 | 2.51 | 1.92 | 1.00
20 | 2.09 | 1.00 | 1.42 | 1.01
SV
5 | 2.10 | 4.22 | 2.65 | 1.00
10 | 2.08 | 2.71 | 1.54 | 1.00
20 | 2.10 | 1.15 | 1.64 | 1.03
Notes: This table shows the effective number of clusters. Numbers represent
means over the different replications from the DGP.

When we use the DPM we can infer the effective number of clusters. This is
achieved as follows. For the $i^{th}$ run of our MCMC algorithm, we compute

$G^{(i)}=\sum_{j}^{J}\mathbb{I}\left(T^{(i)}_{j}>0\right).$

We can then compute the posterior median of these runs to obtain an estimate
of the number of clusters. Table 2 shows the mean over these posterior median
estimates across the different realizations from the DGP. Note that, when the
DGP is Gaussian and the DPM extension is unnecessary, our algorithm is
correctly selecting $G=1$. When the DGP is skewed, the DPM is modelling this
with a mixture of $G=2$ distributions, regardless of VAR dimension or whether
the DGP has SV or not. When the DGP involves a Student t distribution (and to
a lesser extent when the DGP is SV) we are finding an interesting pattern
where the number of clusters is inversely related to the VAR dimension. This
is consistent with our conjecture that, as the VAR dimension increases and
more explanatory variables appear on the right hand side of each equation, the
extra variables can fit some of the fat tailed behavior of the DGP. This
tradeoff between VAR dimension and the need for a non-parametric distribution
is something we find in our forecasting exercise and we will say more about
this issue then.

##  5 Empirical application using US data

###  5.1 Data and specifications

We use US quarterly macroeconomic data from 1960Q1 to 2022Q1 taken from the
FRED database, see McCracken and Ng (2020). A full list of variables along the
transformations we use is given in the Data Appendix. In our forecasting
exercise, we evaluate forecast performance beginning in $1977Q1$ and rely on
the iterated method of forecasting and consider forecast horizons of one-
quarter and one-year.

We work with small ($M=4$), medium ($M=7$) and large ($M=28$) dimensional VARs
with the variables included in each being given in the Data Appendix. All
variables are transformed to be approximately stationary (with detailed
information provided in the Data Appendix). In the forecasting exercise we
evaluate the performance of the models by focusing on the variable-specific
performance for GDP growth, the unemployment rate and inflation. These three
variables form our set of focus variables. We choose a long lag length, $p=5$,
and trust our shrinkage prior to prevent overfitting.

In Sub-section 5.4 we carry out a forecasting exercise to assess whether
adding the DPM to the VAR improves predictive accuracy. To this end, we
compare the performance of the model with non-parametric shocks (BVAR-DPM) to
models with Gaussian errors which we obtain by taking our BVAR-DPM and setting
$G=1$. This is the model described in Sub-section 2.1. All other specification
details, including prior choice, are the same in all of our models. In this
way, we can focus on the specific issue of what the use of the DPM adds to the
model. Note that we are not also including a BVAR with full error covariance
since, as established previously, it is computationally much more burdensome
in larger models and is expected to give results very similar to the BVAR with
$G=1$. We consider two versions of every model, one with SV and one
homoskedastic.

###  5.2 Full sample analysis

In order to illustrate the properties of our models, we begin by carrying out
full sample estimation using a single model: the BVAR-DPM without SV using the
medium data set. This model works well empirically and allows for a
straightforward comparison between impulse responses in the next sub-section.

An advantage of the DPM is that it can be used to estimate the effective
number of components in the Gaussian mixture and our simulation results show
that it does so accurately.

Table 3: Probability of a given effective number of clusters $G=$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
---|---|---|---|---|---|---|---|---|---
| 0.00 | 0.00 | 0.00 | 0.95 | 0.05 | 0.00 | 0.00 | 0.00 | 0.00

Table 3 presents evidence relating to this. It provides quantitative
information on the posterior distribution of the effective number of clusters.
The table suggests that the BVAR-DPM-SV is allocating almost all of the
probability to $G=4$ with the remaining probability allocated to $G=5$. Thus,
we are finding no evidence in favour of the conventional Gaussian VAR with SV
which has $G=1$.

Figure 2: Probability of being in a given cluster over time

The characteristics of the four clusters estimated by the DPM can be seen in
Figure 2. This figure shows the posterior probability that $\delta_{t}=j$ for
$j=1,\dots,4$. To introduce persistence and simplify the discussion, we report
yearly rolling averages of these probabilities.

Cluster 1 is predominant and holds with high probability in most periods.
Cluster 2 is associated with times of high volatility such as periods in the
early 1970s and around 1980 as well as the financial crisis and, to a lesser
extent, the pandemic. Clusters 3 and 4 are completely associated with the
pandemic. Another interesting point to note is that, much of the time Cluster
1 applies with probability near one. But in more volatile times like the
financial crisis and the pandemic, no single cluster holds with probability
near one. For instance, early in the Covid-19 pandemic each of the four
clusters receives roughly equal probability. Thus, in normal times a single
Gaussian distribution suffices to model the error distribution, but in less
stable times a mixture of two or more Gaussians is required.

A deeper understanding of the properties of the clusters can be obtained by
looking at the posteriors of $\bm{\Sigma}_{k}$ for $k=1\ldots 5$. Figure 3
contains box and whisker plots of the posterior of the log of the determinant
of each $\bm{\Sigma}_{k}$. It can be seen that Cluster 1 is the low volatility
cluster whereas Cluster 2 has much higher volatility (as evidenced by much
larger log-determinants). Clusters 2 through 5 also have much more uncertainty
(e.g. wider credible intervals) than Cluster 1.

Figure 3: Log-Determinant of $\bm{\Sigma}_{k}$

###  5.3 Structural impulse responses to a monetary policy shock

Next we investigate the dynamic effects of monetary policy shocks using a
standard identification scheme with a Cholesky ordering where the Federal
Funds rate is ordered above financial market variables and below real
quantities. This ordering implies that real quantities belong to the ’slow-
moving’ block whereas financial markets are ’fast-moving’ (Bernanke et al.,
2005). To economize on space, we focus on the IRFs of the three focus
variables (output growth, unemployment and inflation).

A feature of our model is that, conditional on the mixture indicators
$\delta_{t}$, our model can be interpreted as a constant parameter VAR with a
fully time-varying covariance matrix. This implies that the structural form of
the model features time-varying parameters. To see this, multiply (3) by
$\bm{\Psi}_{t}^{-1}$, the inverse of the lower Cholesky factor of
$\bm{\Xi}_{t}$, from the left. This yields:

$\bm{\Psi}_{t}^{-1}\bm{y}_{t}=\tilde{\bm{A}}_{t}\bm{X}_{t}+\tilde{\bm{\varepsilon}}_{t}.$

Here we let $\tilde{\bm{A}}_{t}=\bm{\Psi}_{t}^{-1}\bm{A}$ and
$\tilde{\bm{\varepsilon}}_{t}\sim\mathcal{N}(\bm{0}_{M},\bm{I}_{M})$. The key
implication is that the structural coefficients of the model are time-varying
and all parameters in the structural form of the model change if the regime
shifts.

(a) GDPC1

(b) UNRATE

(c) CPIAUCSL

Notes: The black lines show the posterior median of the impulse responses for
the different clusters. Figure 4: Median impulse responses to a monetary
policy shock for the different clusters.

We start our analysis by considering median impulse responses over time. Since
within each regime we have a standard VAR with constant coefficients, we
compute the IRFs per regime. This gives a set of $G$ dynamic responses to a
monetary policy shock. To analyze differences in shapes and magnitudes we
focus on the posterior median in the main body of the paper. Results that
include posterior credible intervals and the reactions of the other variables
in $\bm{y}_{t}$ are provided in the Empirical Appendix.

Figure 4 presents posterior medians of the impulse responses implied by each
cluster. In general, the model produces impulse responses that are consistent
with our economic intuition. In response to unexpected increases in the policy
rate, output growth declines and the unemployment rate increases. Inflation,
unexpectedly, increases, pointing towards a price puzzle. It is worth
stressing that this price puzzle is most pronounced in the second cluster,
which is the high volatility cluster.

When we compare differences across clusters, it can be seen that Cluster 2
stands out as implying very different impulse responses but in ways that are
different for the different variables. For instance, all of the clusters are
very similar for long-run impulse responses (e.g. greater than 10 quarters),
Cluster 2 differs greatly from other clusters at short horizons (e.g. less
than one year) for most of the variables. In particular, short-run responses
in Cluster 2 appear to be much more pronounced. This is driven by the fact
that the variances of the structural shocks are much larger and the monetary
shock implies a stronger impact reaction of the Federal funds rate (see panel
(d) of Figure B.1 in the Empirical Appendix). For medium-run responses, the
second cluster also yields responses of GDP growth and unemployment rates
which differ from the remaining clusters. For GDP growth, our results indicate
an overshoot in real activity after around two years whereas the unemployment
reaction appears to be much more persistent. This points towards differences
in the transmission of monetary policy to real activity in turbulent periods
(i.e. periods that are allocated to Cluster 2).

(a) GDPC1

(b) CPIAUCSL

(c) UNRATE

Notes: The black lines show the posterior median of the impulse responses for
the BVAR-DPM. The shaded region is a credible interval. The orange lines
denote credible intervals for the Gaussian BVAR. Both intervals cover
$16^{th}$ to $84^{th}$ percentiles. Figure 5: Median impulse responses to a
monetary policy shock, average across clusters.

Next we turn to the question about the overall effects of monetary shocks on
the three focus variables. The overall impulse responses (i.e. averaged over
the different clusters) for our core variables (i.e. the ones common to all
VARs) are given in Figure 5. These are calculated by taking the impulse
responses for cluster $i$ and weighting it by the posterior mean of
$\eta_{i}$. This figure also contains credible intervals for both our BVAR-DPM
and the BVAR with $G=1$ so that the reader can gauge whether the differences
in impulse responses between our model and the Gaussian-errored equivalent are
substantial in a statistical sense.

The figure suggests that differences between the model with $G=\infty$ and
$G=1$ are not substantial in the sense that the credible intervals of the
BVAR-DPM include the ones of the BVAR in almost all cases. The main exception
is the reaction of inflation. Here, we observe a much stronger immediate
increase (i.e. more evidence for a price puzzle) but also a slightly more
pronounced one-year-ahead decline in inflation. It is also worth stressing
that short-run unemployment reactions of the BVAR with $G=1$ suggest
substantial posterior mass of the IRFs are located below zero, suggesting a
decline in the unemployment rate to a monetary tightening. The BVAR-DPM
allocates appreciably less posterior evidence of declines in unemployment
rates.

###  5.4 Forecasting performance

In this sub-section we present the results of our forecasting exercise. Our
measure of point forecast performance is the mean squared forecast error (MSE)
and our measure of density forecast performance is the average of log
predictive likelihoods (LPL). Tail forecast performance is measured using
absolute quantile scores (QSs). Results are reported relative to a benchmark
model which is the large BVAR with SV and Gaussian errors.

(a) One-quarter-ahead

(b) One-year-ahead

Notes: This figure shows cumulative log predictive likelihoods of the BVAR-DPM
(across model sizes) relative to the large BVAR-G1 model. All models have SV.
Figure 6: Density forecasting performance over time

Our discussion starts by considering the overall density forecasting
performance for the three focus variables over the hold-out period (1977:Q1 to
2022:Q1). Figure 6 presents an overall summary of how our BVAR-DPM with SV
models are performing over time. It contains the cumulative differences in the
joint predictive likelihood for the three main variables being forecast
(inflation, GDP growth and the unemployment rate) for the BVAR-DPM with SV and
the BVAR-SV. Numbers greater than zero indicate that the DPM component of the
model yields accuracy gains whereas negative numbers suggest that it hurts
predictive performance.

When we consider the entire hold-out period, all of our BVAR-DPM with SV
models beat the benchmark for both forecast horizons. In terms of VAR
dimension, it is the medium VAR which forecasts best at both horizons. But
this overall performance masks interesting patterns over time. It can be seen
that the excellent performance of the DPM models is driven largely by the
pandemic and the financial crisis. Clearly, adding the non-parametric
treatment of shocks is particularly useful in such periods. This is because
the likelihood of observing outliers under the predictive density increases if
the training sample includes similar outliers. A model that does not
appropriately capture outliers is not capable of producing abrupt changes in
the predictive variance and thus underestimates the probability of tail
events.

In tranquil periods, the only DPM model which beats the benchmark is the small
BVAR-DPM. Hence, we have a story where, in normal times, smaller, more
parsimonious models suffice but in unusual times larger, non-Gaussian models
forecast better. An advantage of the BVAR-DPM approach is that it can decide
automatically when to introduce non-Gaussian errors and when to simply choose
$G=1$ and use Gaussian errors.

Next we zoom into variable-specific forecasting performance. We do so by
focusing on MSE ratios and LPL differences between a given model and the large
VAR with SV. Table 4 present the relative MSEs and LPLs (in parentheses) for
the three variables being forecast for both homoskedastic and heteroskedastic
models, allowing for a variable-specific examination of forecast performance.

Table 4: Forecasting performance: 1977:Q1 to 2022:Q1.

|  | GDPC1 | UNRATE | CPIAUCSL
---|---|---|---|---
|  | DPM | $G=1$ | DPM | $G=1$ | DPM | $G=1$
One-quarter-ahead
Homosk. | S | 0.951 | 1.219 | 1.007 | 1.535 | 0.912 | 0.951
|  | (-0.019) | (-0.055) | ( 0.550) | (-0.664) | ( 0.065) | ( 0.020)
| M | 0.872 | 1.003 | 0.959 | 1.291 | 0.836 | 0.863
|  | (-0.023) | (-0.014) | ( 0.965) | (-0.596) | ( 0.133) | ( 0.060)
| L | 0.955 | 1.022 | 1.080 | 1.055 | 0.804 | 0.824
|  | (-0.090) | (-0.165) | ( 0.073) | (-0.037) | ( 0.086) | ( 0.070)
SV | S | 0.898 | 0.981 | 0.940 | 1.135 | 0.948 | 0.919
|  | ( 0.182) | ( 0.219) | ( 0.339) | (-0.125) | ( 0.117) | ( 0.163)
| M | 0.802 | 0.919 | 1.023 | 1.082 | 0.867 | 0.865
|  | ( 0.039) | ( 0.278) | ( 0.979) | (-0.265) | ( 0.181) | ( 0.140)
| L | 0.959 | 1.374 | 1.193 | 1.884 | 1.050 | 1.286
|  | (-0.188) | (-1.489) | ( 0.058) | (-2.660) | ( 0.070) | (-1.613)
One-year-ahead
Homosk. | S | 0.950 | 1.007 | 0.978 | 1.243 | 1.045 | 1.039
|  | ( 0.317) | ( 0.267) | ( 0.496) | (-0.155) | ( 0.107) | ( 0.071)
| M | 0.988 | 0.984 | 0.979 | 1.067 | 1.020 | 1.036
|  | ( 0.312) | ( 0.227) | ( 0.516) | (-0.147) | ( 0.159) | ( 0.048)
| L | 0.969 | 0.973 | 0.997 | 0.990 | 0.953 | 0.954
|  | ( 0.163) | (0.133) | ( 0.339) | ( 0.096) | ( 0.072) | ( 0.070)
SV | S | 0.952 | 0.948 | 0.979 | 1.023 | 1.047 | 1.049
|  | ( 0.076) | ( 0.171) | ( 0.118) | ( 0.105) | ( 0.125) | ( 0.106)
| M | 0.951 | 0.950 | 0.999 | 0.995 | 1.005 | 1.026
|  | ( 0.215) | (-0.036) | ( 0.473) | (-0.046) | ( 0.181) | ( 0.120)
| L | 0.966 | 1.149 | 1.024 | 1.351 | 1.018 | 1.292
|  | (-0.013) | (-1.924) | ( 0.363) | (-2.198) | ( 0.120) | (-1.760)

Notes. The table reports mean squared forecast errors relative to the large
BVAR with SV and average log predictive likelihood differences (in
parentheses) between a given model and the large BVAR with SV. The BVAR-SV is
estimated by setting $G=1$. The blue shaded cells are absolute MSEs and LPLs.
Bold numbers indicate the best performing model for a given variable and
horizon. Table 5: Tail Forecasting performance: 1977:Q1 to 2022:Q1

|  | GDPC1 | UNRATE | CPIAUCSL
---|---|---|---|---
|  | DPM | $G=1$ | DPM | $G=1$ | DPM | $G=1$
One-quarter-ahead
Homosk. | S | 0.882 | 0.890 | 1.295 | 1.705 | 0.919 | 0.931
|  | (0.996) | (1.092) | (0.957) | (1.037) | (1.029) | (1.059)
| M | 0.883 | 0.897 | 1.410 | 1.453 | 0.913 | 0.941
|  | (1.038) | (1.016) | (1.085) | (0.974) | (0.943) | (0.942)
| L | 0.944 | 0.948 | 1.436 | 1.261 | 0.941 | 0.944
|  | (1.034) | (1.065) | (0.959) | (0.953) | (0.919) | (0.947)
SV | S | 0.915 | 0.928 | 0.900 | 0.916 | 0.941 | 0.943
|  | (0.832) | (0.881) | (0.966) | (1.035) | (1.004) | (0.945)
| M | 0.941 | 0.914 | 1.462 | 0.940 | 0.942 | 0.973
|  | (1.021) | (0.864) | (1.148) | (1.002) | (0.865) | (0.881)
| L | 1.000 | 0.183 | 1.573 | 0.109 | 1.003 | 0.211
|  | (1.139) | (0.163) | (1.053) | (0.148) | (1.061) | (0.175)
One-year-ahead
Homosk. | S | 0.817 | 0.836 | 0.919 | 0.963 | 1.028 | 1.038
|  | (0.777) | (0.830) | (0.953) | (0.947) | (0.982) | (0.972)
| M | 0.871 | 0.834 | 1.139 | 0.866 | 1.031 | 1.040
|  | (0.865) | (0.805) | (1.143) | (0.940) | (0.951) | (0.953)
| L | 0.887 | 0.890 | 0.897 | 0.896 | 1.006 | 1.012
|  | (0.891) | (0.892) | (0.929) | (0.929) | (0.934) | (0.939)
SV | S | 0.846 | 0.852 | 0.905 | 0.924 | 0.997 | 1.022
|  | (0.817) | (0.844) | (0.955) | (1.036) | (0.957) | (0.991)
| M | 0.888 | 0.861 | 1.211 | 0.944 | 1.021 | 0.997
|  | (0.944) | (0.843) | (1.208) | (1.017) | (0.985) | (0.967)
| L | 1.019 | 0.216 | 1.179 | 0.162 | 1.031 | 0.204
|  | (1.072) | (0.199) | (1.112) | (0.191) | (1.017) | (0.208)

Notes. The table reports relative quantile scores (10% and 90 %, in
parentheses) to the large BVAR with SV. The BVAR-SV is estimated by setting
$G=1$. The blue shaded cells are absolute quantile scores. Bold numbers
indicate the best performing model for a given variable and horizon.

With some exceptions, it can be seen that DPM models perform well for all
variables. The gains are particularly pronounced for point forecasts, where
GDP growth and inflation forecasts are sometimes close to 20 percent more
precise than the forecast produced by the benchmark model. The only variable
where gains are muted (and sometimes even negative) is the unemployment rate.
In this case, we sometimes find slight improvements in point forecasting
accuracy but these never exceed 10 percent. When we consider density forecasts
we often find positive numbers in the parentheses. This indicates that the
gains in terms of point forecasts also translate into more accurate density
predictions.

To analyze whether the DPM component improves forecasts we can compare each of
the BVAR-DPM models with the equivalent model with $G=1$. For all three
variables, two forecast horizons and for both forecast metrics, we almost
always find the DPM version of the model to forecast better than the Gaussian
one. In addition, versions of all models with SV tend to forecast better than
homoskedastic versions of each model. This corroborate the finding in the
forecasting literature that points towards the necessity for using flexible
assumptions on the shocks to improve forecasts (Clark, 2011; Huber and
Feldkircher, 2019; Carriero et al., 2022b).

In terms of VAR dimension, it is interesting to note that the strongest
performance for the DPM models arises with small and medium VARs. With large
VARs, the evidence is more mixed both with regards to the need for DPM and
with regards to the need for SV. We conjecture that in the large VAR, the
explanatory power of the right hand side variables can mop up some (but not
all) of the need to allow for non-Gaussianity or volatility change.

Table 5 presents evidence on tail forecasting performance. The table shows 10%
and 90% QSs relative to the large VAR with SV so that numbers smaller than one
indicate that a given model produces more accurate tail forecasts than the
benchmark.

We would expect the DPM to be particularly good in capturing tail behavior in
a way that the Gaussian model cannot. And, with some exceptions, Table 5
confirms this expectation. The pattern of results is similar to those found in
Table 4, but are slightly stronger in favor of BVAR-DPM models. More
specifically, we again find strong gains of the different models relative to
the large BVAR with SV. Depending on the model size, we also find accuracy
gains from using the DPM specification relative to setting $G=1$. These gains
are more pronounced for one-year-ahead tail forecasts of output growth and the
unemployment rate and smaller-sized models.

##  6 Conclusions

In this paper we propose a new specification for the errors in a VAR which
takes a particular additive form involving two components. The first is a
homoskedastic error and a conventional inverse Wishart prior can be used for
its covariance matrix. The second is a diagonal error covariance matrix with
diagonal elements following SV processes. We show that, by adopting this
additive form, we gain two major advantages. First, computation is much faster
since it allows for equation-by-equation estimation. Second, posterior and
predictive inference does not depend on the way the variables are ordered in
the VAR. We then extend this model to allow the first error to follow a DPM.
We discuss, both theoretically and empirically, the great flexibility that is
obtained by doing so. In addition, we develop a computationally fast MCMC
algorithm which allows for posterior and predictive inference in high
dimensional non-parametric VARs. Our empirical results, using artificial and
real data, show the flexibility of our approach. In our forecasting exercise,
we find that the added flexibility of our BVAR-DPM in particularly useful in
extreme times such as the financial crisis and the Covid-19 pandemic.

## References

  * (1)
  * Antolin-Diaz et al. (2017) Antolin-Diaz, J., T. Drechsel, and I. Petrella (2017): “Tracking the Slowdown in Long-Run GDP Growth,” _The Review of Economics and Statistics_ , 99(2), 343–356.
  * Arias et al. (forthcoming) Arias, J. E., J. F. Rubio-Ramirez, and M. Shin (forthcoming): “Macroeconomic forecasting and variable ordering in multivariate stochastic volatility models,” _Journal of Econometrics_.
  * Bernanke et al. (2005) Bernanke, B. S., J. Boivin, and P. Eliasz (2005): “Measuring the effects of monetary policy: a factor-augmented vector autoregressive (FAVAR) approach,” _The Quarterly Journal of Economics_ , 120(1), 387–422.
  * Bhattacharya et al. (2016) Bhattacharya, A., A. Chakraborty, and B. K. Mallick (2016): “Fast sampling with Gaussian scale mixture priors in high-dimensional regression,” _Biometrika_ , 103, 985–991.
  * Braun (2021) Braun, R. (2021): “The importance of supply and demand for oil prices: evidence from non-Gaussianity,” Bank of England working papers 957, Bank of England.
  * Brown and Griffin (2010) Brown, P. J., and J. E. Griffin (2010): “Inference with normal-gamma prior distributions in regression problems,” _Bayesian Analysis_ , 5(1), 171 – 188.
  * Carriero et al. (2022a) Carriero, A., J. Chan, T. E. Clark, and M. Marcellino (2022a): “Corrigendum to “Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors”[Journal of Econometrics 212 (1)(2019) 137–154],” _Journal of Econometrics_ , 227(2), 506–512.
  * Carriero et al. (2016) Carriero, A., T. Clark, and M. Marcellino (2016): “Common drifting volatility in large Bayesian VARs,” _Journal of Business & Economic Statistics_, 34(3), 375–390.
  * Carriero et al. (2019) (2019): “Large Bayesian vector autoregressions with stochastic volatility and non-conjugate priors,” _Journal of Econometrics_ , 212(1), 137–154.
  * Carriero et al. (2022b) Carriero, A., T. E. Clark, M. Marcellino, and E. Mertens (2022b): “Addressing COVID-19 outliers in BVARs with stochastic volatility,” _Review of Economics and Statistics_ , pp. 1–38.
  * Chan et al. (2021) Chan, J. C., G. Koop, and X. Yu (2021): “Large Order-Invariant Bayesian VARs with Stochastic Volatility,” _arXiv preprint arXiv:2111.07225_.
  * Chiu et al. (2017) Chiu, C.-W. J., H. Mumtaz, and G. Pintér (2017): “Forecasting with VAR models: Fat tails and stochastic volatility,” _International Journal of Forecasting_ , 33(4), 1124–1143.
  * Clark (2011) Clark, T. (2011): “Real-time density forecasts from Bayesian vector autoregressions with stochastic volatility,” _Journal of Business & Economic Statistics_, 29(3), 327–341.
  * Clark et al. (forthcoming) Clark, T. E., F. Huber, G. Koop, M. Marcellino, and M. Pfarrhofer (forthcoming): “Tail Forecasting with Multivariate Bayesian Additive Regression Trees,” _International Economic Review_.
  * Cogley and Sargent (2005) Cogley, T., and T. J. Sargent (2005): “Drifts and volatilities: monetary policies and outcomes in the post WWII US,” _Review of Economic Dynamics_ , 8(2), 262–302.
  * Dunson and Xing (2009) Dunson, D. B., and C. Xing (2009): “Nonparametric Bayes Modeling of Multivariate Categorical Data,” _Journal of the American Statistical Association_ , 104(487), 1042–1051.
  * Escobar and West (1995) Escobar, M. D., and M. West (1995): “Bayesian density estimation and inference using mixtures,” _Journal of the american statistical association_ , 90(430), 577–588.
  * Fox and Dunson (2015) Fox, E. B., and D. B. Dunson (2015): “Bayesian Nonparametric Covariance Regression,” _Journal of Machine Learning Research_ , 16(77), 2501–2542.
  * Frühwirth-Schnatter and Malsiner-Walli (2019) Frühwirth-Schnatter, S., and G. Malsiner-Walli (2019): “From here to infinity: sparse finite versus Dirichlet process mixtures in model-based clustering,” _Advances in data analysis and classification_ , 13, 33–64.
  * Frühwirth-Schnatter et al. (2004) Frühwirth-Schnatter, S., R. Tüchler, and T. Otter (2004): “Bayesian Analysis of the Heterogeneity Model,” _Journal of Business & Economic Statistics_, 22(1), 2–15.
  * Hauzenberger et al. (2022) Hauzenberger, N., F. Huber, G. Koop, and L. Onorante (2022): “Fast and flexible Bayesian inference in time-varying parameter regression models,” _Journal of Business & Economic Statistics_, 40(4), 1904–1918.
  * Huber and Feldkircher (2019) Huber, F., and M. Feldkircher (2019): “Adaptive shrinkage in Bayesian vector autoregressive models,” _Journal of Business & Economic Statistics_, 37(1), 27–39.
  * Huber et al. (2020) Huber, F., G. Koop, L. Onorante, M. Pfarrhofer, and J. Schreiner (2020): “Nowcasting in a pandemic using non-parametric mixed frequency VARs,” _Journal of Econometrics_ , in-press.
  * Huber and Rossini (2020) Huber, F., and L. Rossini (2020): “Inference in Bayesian Additive Vector Autoregressive Tree models,” _arXiv_ , 2006.16333.
  * Kalli et al. (2011) Kalli, M., J. E. Griffin, and S. G. Walker (2011): “Slice sampling mixture models,” _Statistics and computing_ , 21, 93–105.
  * Karlsson et al. (2023) Karlsson, S., S. Mazur, and H. Nguyen (2023): “Vector autoregression models with skewness and heavy tails,” _Journal of Economic Dynamics and Control_ , 146, 104580.
  * Kastner and Frühwirth-Schnatter (2014) Kastner, G., and S. Frühwirth-Schnatter (2014): “Ancillarity-sufficiency interweaving strategy (ASIS) for boosting MCMC estimation of stochastic volatility models,” _Computational Statistics & Data Analysis_, 76, 408–423.
  * Kastner and Huber (2021) Kastner, G., and F. Huber (2021): “Sparse Bayesian Vector Autoregressions in huge dimensions,” _Journal of Forecasting_ , 39, 1142–1165.
  * Lau and Green (2007) Lau, J. W., and P. J. Green (2007): “Bayesian model-based clustering procedures,” _Journal of Computational and Graphical Statistics_ , 16(3), 526–558.
  * Malsiner-Walli et al. (2017) Malsiner-Walli, G., S. Frühwirth-Schnatter, and B. Grün (2017): “Identifying Mixtures of Mixtures Using Bayesian Estimation,” _Journal of Computational and Graphical Statistics_ , 26(2), 285–295, PMID: 28626349.
  * McCracken and Ng (2020) McCracken, M., and S. Ng (2020): “FRED-QD: A quarterly database for macroeconomic research,” Discussion paper, National Bureau of Economic Research.
  * Park and Casella (2008) Park, T., and G. Casella (2008): “The bayesian lasso,” _Journal of the American Statistical Association_ , 103(482), 681–686.
  * Stock and Watson (2007) Stock, J. H., and M. W. Watson (2007): “Why Has U.S. Inflation Become Harder to Forecast?,” _Journal of Money, Credit and Banking_ , 39(S1), 3–33.
  * Wu and Koop (2023) Wu, P., and G. Koop (2023): “Fast, order-invariant Bayesian inference in VARs using the eigendecomposition of the error covariance matrix,” _manuscript_.
  * Yau and Holmes (2011) Yau, C., and C. Holmes (2011): “Hierarchical Bayesian nonparametric mixture models for clustering with variable relevance determination,” _Bayesian Analysis_ , 6(2), 329 – 351.

##  Appendix A Data Appendix

Table A1: Description of the Dataset

FRED Code | Description | Transformation Codes | S | M | L |
---|---|---|---|---|---|---
GDPC1 | Real Gross Domestic Product | 5 | X | X | X |
PCECC96 | Real Personal Consumption Expenditures | 5 |  |  | X |
FPIx | Real private fixed investment | 5 |  |  | X |
GCEC1 | Real Government Consumption Expenditures and Gross Investment | 5 |  |  | X |
INDPRO | IP:Total index Industrial Production Index (Index 2012=100) | 5 |  |  | X |
CUMFNS | Capacity Utilization: Manufacturing (SIC) (Percent of Capacity) | 1 |  |  | X |
PAYEMS | Emp:Nonfarm All Employees: Total nonfarm (Thousands of Persons) | 5 |  |  | X |
CE16OV | Civilian Employment (Thousands of Persons) | 5 |  |  | X |
UNRATE | Civilian Unemployment Rate (Percent) | 2 | X | X | X |
AWHMAN | Average Weekly Hours of Production and Nonsupervisory Employees: Manufacturing (Hours) | 1 |  |  | X |
CES0600000007 | Average Weekly Hours of Production and Nonsupervisory Employees: Goods-Producing | 2 |  |  | X |
CLAIMSx | Initial Claims | 5 |  |  | X |
GDPCTPI | Gross Domestic Product: Chain-type Price Index | 6 |  |  | X |
CPIAUCSL | Consumer Price Index for All Urban Consumers: All Items | 6 | X | X | X |
PPIACO | Producer Price Index for All Commodities | 6 |  |  | X |
WPSID61 | Producer Price Index by Commodity Intermediate Materials: Supplies & Components | 6 |  |  | X |
WPSID62 | Producer Price Index: Crude Materials for Further Processing | 6 |  |  | X |
COMPRNFB | Nonfarm Business Sector: Real Compensation Per Hour (Index 2012=100) | 5 |  |  | X |
ULCNFB | Nonfarm Business Sector: Unit Labor Costs for All Workers | 5 |  |  | X |
CES0600000008 | Average Hourly Earnings of Production and Nonsupervisory Employees: | 6 |  | X | X |
FEDFUNDS | Effective Federal Funds Rate (Percent) | 2 | X | X | X |
BAA10YM | Moody’s Seasoned Baa Corporate Bond Yield Relative to Yield on 10-Year Treasury | 1 |  | X | X |
GS10TB3Mx | 10-Year Treasury Constant Maturity Minus 3-Month Treasury Bill, secondary market | 1 |  | X | X |
CPF3MTB3Mx | 3-Month Commercial Paper Minus 3-Month Treasury Bill, secondary market | 1 |  |  | X |
M2REAL | Real M2 Money Stock | 5 |  |  | X |
BUSLOANSx | Real Commercial and Industrial Loans, All Commercial Banks | 5 |  |  | X |
CONSUMERx | Real Consumer Loans at All Commercial Banks | 5 |  |  | X |
S.P.500 | S&P’s Common Stock Price Index: Composite | 5 |  | X | X |

Notes: This table provides an overview of the dataset employed. The
transformation codes are applied to each time series $\bm{Y}_{j}$ and
described in McCracken and Ng (2020): (1) no transformation; (2) $\Delta
y_{jt}$; (3) $\Delta^{2}y_{jt}$; (4) $\log(y_{jt})$; (5) $\Delta\log(y_{jt})$;
(6) $\Delta^{2}\log(y_{jt})$; (7) $\Delta(y_{jt}/y_{jt-1}-1)$. ’X’ marks the
inclusion of a variable into one of the datasets.

##  Appendix B Additional empirical results

###  B.1 Impulse responses across clusters

(a) GDPC1

(b) CPIAUCSL

(c) UNRATE

(d) FEDFUNDS

Notes: The black solid line refers to the posterior mean IRF from the BVAR-DPM
with SV across different clusters. The gray shaded area are the 16th/84th
credible intervals of the BVAR-DPM while the orange lines refer to the 16th
and 84th of the BVAR with $G=1$.

(a) GS10

(b) S&P 500

Notes: The black solid line refers to the posterior mean IRF from the BVAR-DPM
with SV across different clusters. The gray shaded area are the 16th/84th
credible intervals of the BVAR-DPM while the orange lines refer to the 16th
and 84th of the BVAR with $G=1$. Figure B.2: Impulse responses to a monetary
policy shock for the different clusters.