Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
143,790 characters · 27 sections · 20 citation commands
Large SVARs
\hypersetup{pageanchor=false}
\hypersetup{pageanchor=true} \onehalfspacing \newgeometry{tmargin=1.25in,bmargin=1.25in,lmargin=1.25in,rmargin=1.25in} \onehalfspacing
The growing availability of large datasets has led to a renewed interest in the use of large-scale time-series models in economics. In univariate settings, \citet*{illusionofsparsity} show that densely parameterized models equipped with appropriate shrinkage priors typically outperform sparse alternatives in terms of predictive accuracy. In multivariate settings, \citet*{banbura2010} and \citet*{koop2013} have demonstrated that Bayesian shrinkage enables the estimation of large vector autoregressions (VARs) without compromising out-of-sample performance, and more recently, \citet*{crump2025large} show how this approach can be used for monetary policy analysis. These insights have direct implications for structural vector autoregressions (SVARs)---one of the major workhorses for studying the propagation of structural shocks in macroeconomics. If prediction tasks are better handled by large models, then inference on structural shocks using SVARs should similarly benefit from broader information sets, in line with early arguments in favor of using large SVARs to understand the macroeconomic effects of monetary policy \citep*[e.g.,][]{leeper1996,bernankefavar}.
Within the SVAR paradigm, sign restrictions have become a particularly popular method for identifying the parameters of interest, typically impulse responses. The conventional Bayesian approach to implementing sign restrictions—pioneered by \citet*{faust1998robustness}, \citet*{canova2002monetary}, and \citet*{uhlig2005effects} and extended by \citet*{RRWZ2010}—relies on sampling from the reduced-form posterior and a uniform prior over the set of orthogonal matrices, combined with an accept-reject approach to impose the restrictions. The method is straightforward to implement and produces independent draws. However, it becomes increasingly infeasible in large systems due to the vanishing probability of sampling admissible orthogonal matrices—especially as the number of sign restrictions needed to identify multiple structural shocks increases, tightening the identified set.\footnote{We use the term admissible to refer to orthogonal matrices that satisfy the sign restrictions given the reduced-form parameters.} Recent work by \citet*{chan2025largestructuralvarsmultiple} introduces algorithmic refinements by exploiting symmetry and permutation invariance in the space of orthogonal matrices, but even these improved methods also face limitations under tight identification. Notably, such computational burdens are not unique to high-dimensional models. These challenges have coincided with advances in identification strategies that also lead to tighter identified sets---even in low-dimensional SVARs---such as ranking restrictions and elasticity bounds \citep*{kilian2012agnostic,AmirAhmadiDrautzburg}, as well as narrative sign restrictions \citep*{antolinrubionarrative,sydneymaserena}. These approaches go beyond traditional sign restrictions while preserving their intuitive appeal and further reducing the volume of the admissible space of orthogonal matrices. Together, these trends---the adoption of larger information sets and the use of tighter identification strategies---and the limitations of accept-reject sampling methods underscore the need for alternative algorithms for Bayesian inference in SVARs identified with sign restrictions.
In this paper, we break with the state-of-the-art accept-reject tradition and show that embedding the elliptical slice sampling method of \citet*{Murray2010} within a Gibbs sampler---yielding an elliptical slice within Gibbs sampler---delivers substantial gains in computational speed, rendering previously infeasible applications tractable. Like in the conventional approach, using the uniform prior over the set of orthogonal matrices, our goal is to draw from the posterior distribution of the orthogonal reduced-form parameters conditional on sign restrictions. However, by eliminating the accept-reject step and directly conditioning on the sign restrictions within the Gibbs sampler, our algorithm overcomes the bottlenecks that arise under tight identification---thus enabling dynamic structural analysis with big data and rich identification schemes. We also show that the algorithm is well-defined by establishing that the stationary distribution of the elliptical slice within Gibbs sampler is the posterior distribution of interest. To illustrate the advantages of our approach, we consider a very simple example similar to the one in \citet*{granziera2018inference} and demonstrate that the efficiency of the {state-of-the-art} accept-reject algorithm hinges critically on the size of the identified set. As the identified set becomes tighter, the accept-reject algorithm slows down dramatically. In contrast, our Gibbs sampler efficiently shrinks the support of the candidate impulse responses toward the identified set, maintaining speed even under stringent restrictions.
For clarity and comparability with the literature, when describing our proposed Gibbs sampler algorithm, we adopt the conjugate normal-inverse-Wishart prior for the reduced-form parameters as our baseline. While this prior is popular due to its analytical convenience, it precludes cross-variable shrinkage. To address this limitation, we extend our algorithm to accommodate alternative priors, including the independent normal-inverse-Wishart and the asymmetric prior of \citet*{Chanassymetric}, both of which support cross-variable shrinkage. We evaluate the performance of our approach using two applications. In the first, we replicate \citet*{kilianmurphy2014}, a model of the world oil market in which the standard accept-reject algorithm fails. To address this infeasibility, \citet*{kilianmurphy2014} adopt an approach similar to that of \citet*{chan2025largestructuralvarsmultiple}, exploiting permutations and sign alternations. Our algorithm handles this application multiple times faster than the accept-reject approach, though the computational times of both approaches are within a range most practitioners would find acceptable. However, once we tighten the identified set by adding a restriction on the price elasticity of oil demand---motivated by \citet*{CALDARA20191}---the difference in performance becomes substantial: The {state-of-the-art} accept-reject algorithm moves from requiring about 20 minutes to produce 1,000 draws to nearly eight hours, whereas the computational time of our Gibbs sampler remains roughly constant, increasing only from about 2 minutes to 5 minutes for the same number of effective draws. In the second application, we revisit the structural analysis in \citet*{chan2025largestructuralvarsmultiple}, who use crump2025large’s crump2025large large SVAR model of the U.S. economy to identify eight structural shocks. While the latter uses the Minnesota prior, the former relies on the asymmetric prior of \citet*{Chanassymetric}. To simplify the comparison, we revert to the Minnesota prior when applying both the accept-reject algorithm and the Gibbs sampler. We show that as the number of shocks under analysis increases, the efficiency of the algorithm in \citet*{chan2025largestructuralvarsmultiple} declines markedly, eventually becoming impractical. With ten shocks, it would take several days to obtain 1,000 draws. In contrast, the computational time of our Gibbs sampler is largely insensitive to the number of identified structural shocks. Even with ten shocks, it would take only a few minutes to obtain 1,000 effective draws. We also show that these striking differences are robust to using the asymmetric prior instead of the Minnesota prior.
We also highlight the contemporaneous contribution of \citet*{readzhu2025}, who propose an algorithm based on slice sampling, but it is limited to the conditionally uniform prior described in \citet*{uhlig2017shocks} and \citet*{AmirAhmadiDrautzburg}. While this prior can deliver substantial speed gains without the need for a Gibbs sampler, it does not satisfy the requirements set out in \citet*{arias2025uniform}. As we explain in Section (ref), this has important consequences: it implicitly alters the prior over the impulse responses in a way that depends on the identification scheme. This entanglement of inference and identification makes it difficult to know whether differences in posterior inference reflect genuine differences in identification or are merely artifacts of unintended changes in the prior distribution. By contrast, the uniform prior over orthogonal matrices satisfies the requirements in \citet*{arias2025uniform} and, therefore, guarantees that inference remains invariant to the set of imposed restrictions and allows researchers to cleanly separate the role of prior beliefs from the role of identification assumptions.
The remainder of the paper is organized as follows. Sections (ref) through (ref) introduce the SVAR model, the sign restrictions, and the baseline conjugate uniform-normal-inverse-Wishart prior. Section (ref) describes the problem in a simple environment. Section (ref) presents our elliptical slice within Gibbs sampler and outlines its theoretical properties. Section (ref) applies the algorithm to two empirical settings: a small SVAR model of the world oil market and a large SVAR model of the U.S. economy. Section (ref) shows the shortcomings of the conditionally uniform prior. Section (ref) concludes. The Appendix adapts the algorithm to two popular priors in SVAR analysis, the independent normal-inverse-Wishart prior and the asymmetric prior proposed by \citet*{Chanassymetric}, and provides a robustness analysis.
Consider the SVAR with the general form,
where $\operatorname{\mathbf{A}}_+^{\prime}=\left[\operatorname{\mathbf{A}}_{1}^{\prime} \;\cdots\;\operatorname{\mathbf{A}}_{p}^{\prime} \;\;{\operatorname{\mathbf{c}}}^{\prime}\right]$ and $\operatorname{\mathbf{x}}_t^{\prime}=\left[\operatorname{\mathbf{y}}_{t-1}^{\prime} \;\cdots\; \operatorname{\mathbf{y}}_{t-p}^{\prime} \;\;1 \right]$ for $1 \leq t \leq T$, and where $\operatorname{\mathbf{y}}_t$ is an $n \times 1$ vector of endogenous variables, $\operatorname{\boldsymbol{\varepsilon}}_t$ is an $n \times 1$ vector of exogenous structural shocks, $\operatorname{\mathbf{A}}_{\ell}$ is an $n \times n$ matrix of parameters for $0 \leq \ell \leq p$ with $\operatorname{\mathbf{A}}_0$ invertible, $\operatorname{\mathbf{c}}$ is a $1 \times n$ vector of parameters, $p$ is the lag length, and $T$ is the sample size. Hence, the dimension of $\operatorname{\mathbf{A}}_+$ is $m \times n$, where $m = np + 1$. The vector $\operatorname{\boldsymbol{\varepsilon}}_t$, conditional on past information and the initial conditions $\operatorname{\mathbf{y}}_0, \ldots, \operatorname{\mathbf{y}}_{1-p}$, is Gaussian with mean zero and covariance matrix $\operatorname{\mathbf{I}}_n$, the $n \times n$ identity matrix.
The reduced-form representation implied by Equation (ref) is
where $\operatorname{\mathbf{B}} =\operatorname{\mathbf{A}}_+\operatorname{\mathbf{A}}_0^{-1}$, $\operatorname{\mathbf{u}}_t^{\prime}=\operatorname{\boldsymbol{\varepsilon}}_t^{\prime}\operatorname{\mathbf{A}}_0^{-1}$, and $\mathbb{E}\left[\operatorname{\mathbf{u}}_t\operatorname{\mathbf{u}}_t^{\prime}\right]=\bfS=\left(\operatorname{\mathbf{A}}_0\operatorname{\mathbf{A}}_0^{\prime}\right)^{-1}$. The matrices $\operatorname{\mathbf{B}}$ and $\bfS$ are the reduced-form parameters, while $\operatorname{\mathbf{A}}_0$ and $\operatorname{\mathbf{A}}_+$ are the structural parameters. While $\operatorname{\mathbf{B}}$ is an $m \times n$ matrix, $\bfS$ belongs to the set $\mathcal{S}(n)$, which is the set of $n \times n$ positive definite matrices. It will be useful to partition $\operatorname{\mathbf{B}}$ as follows: $\operatorname{\mathbf{B}}=\left[\operatorname{\mathbf{B}}_{1}^{\prime} \;\cdots\; \operatorname{\mathbf{B}}_{p}^{\prime} \;\;{\bf{d}}^{\prime}\right]^{\prime}$ where $\operatorname{\mathbf{B}}_{\ell}$ is an $n \times n$ matrix of parameters for $1 \le \ell \le p$, and ${\bf{d}}$ is a $1 \times n$ vector of parameters.
It is well known that for linear Gaussian models of the type studied in this paper, $(\operatorname{\mathbf{A}}_0, \operatorname{\mathbf{A}}_+)$ and $(\tilde\operatorname{\mathbf{A}}_0, \tilde\operatorname{\mathbf{A}}_+)$ are observationally equivalent if and only if they have the same reduced-form representation. This implies that the structural parameters $(\operatorname{\mathbf{A}}_0, \operatorname{\mathbf{A}}_+)$ and $(\tilde\operatorname{\mathbf{A}}_0, \tilde\operatorname{\mathbf{A}}_+)$ are observationally equivalent if and only if $\operatorname{\mathbf{A}}_0 = \tilde\operatorname{\mathbf{A}}_0 \operatorname{\mathbf{Q}}$ and $\operatorname{\mathbf{A}}_+ = \tilde\operatorname{\mathbf{A}}_+ \operatorname{\mathbf{Q}}$ for some $\operatorname{\mathbf{Q}} \in \mathcal{O}(n)$, where $\mathcal{O}(n)$ is the set of all $n \times n$ orthogonal matrices. To solve the identification problem, one often imposes sign restrictions on either the structural parameters or some function of the structural parameters, such as the impulse responses. To simplify the notation, we summarize the sign restrictions by $\mathbf{S}_S(\operatorname{\mathbf{A}}_0, \operatorname{\mathbf{A}}_+) > \bf0$, where $\bf0$ is a column vector of the appropriate dimension, and let $\left[\mathbf{S}_S(\operatorname{\mathbf{A}}_0, \operatorname{\mathbf{A}}_+) > \bf0\right]$ be an indicator function that equals one if the sign restrictions are satisfied and zero otherwise. We will consider continuous functions $\mathbf{S}_S$.
Equation (ref) represents the SVAR in terms of the structural parameterization, which is characterized by $(\operatorname{\mathbf{A}}_0,\operatorname{\mathbf{A}}_+)$. The SVAR can alternatively be written in what we call the orthogonal reduced-form parameterization; see \citet*{arrw2018}. This parameterization is characterized by the reduced-form parameters $(\operatorname{\mathbf{B}},\bfS)$ together with an orthogonal matrix $\operatorname{\mathbf{Q}}$, and is given by the following equation:
where the $n \times n$ matrix $h(\bfS)$ is any decomposition of the covariance matrix $\bfS$ satisfying $h(\bfS)^{\prime}h(\bfS) = \bfS$. We take $h$ to be the Cholesky decomposition, though any differentiable decomposition would suffice.
Given Equations (ref) and (ref), we can define a mapping between $(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}})$ and $(\operatorname{\mathbf{A}}_0, \operatorname{\mathbf{A}}_+)$ by
This mapping makes clear how the structural parameters depend on the reduced-form parameters and orthogonal matrices. Given the reduced-form parameters, each value of $\operatorname{\mathbf{Q}} \in \mathcal{O}(n)$ can be viewed as a particular choice among observationally equivalent structural parameters. Thus, we can always write the sign restrictions in terms of the orthogonal reduced-form parameterization. Hence, let $\left[\mathbf{S}_R(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}}) > {\bf0}\right]$ be an indicator function in terms of the orthogonal reduced-form parameterization that equals one if the sign restrictions are satisfied and zero otherwise, where $\mathbf{S}_R(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}})=\mathbf{S}_S(f(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}}))$. Because both $\mathbf{S}_S$ and $f$ are continuous, $\mathbf{S}_R$ is also continuous.
We can also define the impulse responses. Let $\operatorname{\mathbf{u}}_t = \operatorname{\mathbf{L}}_0 \operatorname{\boldsymbol{\varepsilon}}_t$ for $1 \le t \le T$, where $\operatorname{\mathbf{L}}_0$ is an $n \times n$ invertible matrix that represents the impulse responses at horizon zero. Given $\operatorname{\mathbf{L}}_0$ and $\operatorname{\mathbf{B}}$, it is possible to obtain the impulse responses beyond horizon zero recursively as
We combine the impulse responses from horizons one through $p$ and the constant term ${\operatorname{\mathbf{c}}}$ into a single matrix, $\operatorname{\mathbf{L}}_+ = \left[ \operatorname{\mathbf{L}}_1^{\prime} \;\cdots\; \operatorname{\mathbf{L}}_p^{\prime} \;\; \operatorname{\mathbf{c}}^{\prime} \right]^{\prime}$, where the maximum horizon of the impulse response in ${\bf L}_+$ matches the lag length in Equation (ref). The impulse response parameterization is characterized by $(\operatorname{\mathbf{L}}_0,\operatorname{\mathbf{L}}_+)$. Given the function $f$ and Equation (ref), we can also define a mapping from $(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}})$ to $(\operatorname{\mathbf{L}}_0, \operatorname{\mathbf{L}}_+)$ by
where ${\bf L}_\ell(\operatorname{\mathbf{B}}, \bfS, \operatorname{\mathbf{Q}})$ for $1 \le \ell \le p$ is implicitly defined by Equation (ref), and $\operatorname{\mathbf{d}} = \operatorname{\mathbf{c}}\operatorname{\mathbf{L}}_0^{\prime}$. The functions $f$ and $\phi$ are invertible, and $f$, $\phi$, and their inverses are differentiable.
For the reduced-form representation in Equation (ref), the normal-inverse-Wishart family of distributions is conjugate. A conjugate normal-inverse-Wishart distribution over the reduced-form parameters is characterized by four parameters: a scalar $\nu\geq n$, an $n\times n$ symmetric and positive definite matrix $\bfPhi$, an $m\times n$ matrix $\bfPsi$, and an $m\times m$ symmetric and positive definite matrix $\bfOme$. We denote this distribution by $NIW(\nu,\bfPhi,\bfPsi,\bfOme)$ and its density by $NIW_{(\nu,\bfPhi,\bfPsi,\bfOme)}(\operatorname{\mathbf{B}},\bfS)$. Furthermore,
where
If the prior {distribution} {over the reduced-form parameters} is $NIW(\bar{\nu},\bar\bfPhi,\bar\bfPsi,\bar\bfOme)$, then the posterior {distribution} over the reduced-form parameters is $NIW(\tilde{\nu},\tilde\bfPhi,\tilde\bfPsi,\tilde\bfOme)$, where
for $\operatorname{\mathbf{Y}}=[\operatorname{\mathbf{y}}_1\;\;\cdots\;\;\operatorname{\mathbf{y}}_T]^\prime$ and $\operatorname{\mathbf{X}}=[\operatorname{\mathbf{x}}_1\;\;\cdots\;\;\operatorname{\mathbf{x}}_T]^\prime$.
The conjugate normal-inverse-Wishart prior is widely used in Bayesian VARs due to its computational convenience and desirable properties \citep*[see][]{uhlig1994macroeconomists,faust1998robustness,uhlig2005effects,sims1998bayesian,RRWZ2010,kilian2012agnostic,kilianmurphy2014}. When combined with the conventional accept-reject approach, it produces independent draws from the posterior, which makes it especially attractive. However, it also imposes a Kronecker structure on the prior distribution of $\operatorname{\mathbf{B}}$, thereby constraining its covariance matrix, and it rules out {cross-variable} shrinkage. Consequently, researchers {often} consider: (i) the independent normal-inverse-Wishart prior, which avoids the Kronecker covariance structure and allows for greater flexibility; and (ii) the asymmetric priors proposed by \citet*{Chanassymetric}, {which are becoming popular because they accommodate} cross-variable shrinkage. {Although} we present the methodology using the conjugate normal-inverse-Wishart prior, due to its prevalence in the literature, in Online Appendix (ref) we also adapt the algorithm to these alternative priors.
Given the results in \citet*{arias2025uniform}, we will combine the conjugate prior with the following uniform density over the set of orthogonal matrices:
where $\int_{\mathcal{O}(n)} \kappa d\operatorname{\mathbf{Q}} = 1$. This choice can be motivated by the fact that it assigns equal prior weight to both observationally equivalent vectors of impulse responses and observationally equivalent structural parameters \citep*[see][]{arias2025uniform}. We call this combination the conjugate uniform-normal-inverse-Wishart distribution over the orthogonal reduced-form parameterization; denote it by $UNIW(\nu,\bfPhi,\bfPsi,\bfOme)$, and denote its density over the orthogonal reduced-form parameterization by $UNIW_{(\nu,\bfPhi,\bfPsi,\bfOme)}(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})$. It is the case that
Our objective will be to draw from the posterior of the orthogonal reduced-form parameters conditional on the sign restrictions,
where
and then use $f$ and $\phi$ to transform the draws to the desired impulse responses. The traditional approach to obtain draws from Equation (ref) uses the following accept-reject algorithm:
As mentioned above, while this algorithm has been widely adopted, it is well known that there are cases in which the identified set is narrow, limiting the efficiency of the algorithm \citep*[{see, e.g.,}][]{kilianmurphy2014,Baumeister_Hamilton_2024,chan2025largestructuralvarsmultiple,readzhu2025}. In the next section, we use a simple example to show its shortcomings. We will also show how a carefully designed elliptical slice within Gibbs sampler is not subject to this limitation and delivers dramatic speed gains. Importantly, \citet*{chan2025largestructuralvarsmultiple} present a numerically efficient version of Algorithm (ref). Therefore, when comparing our algorithm to the traditional accept-reject approach, we will use this efficient version as the benchmark.
For the purposes of demonstrating the limitations of the accept-reject approach, it suffices to work with a simple example similar to the one explored by \citet*{granziera2018inference}. Thus, consider the following SVAR, with $n = 2$, $p = 0$, and ${\operatorname{\mathbf{c}}}={\bf 0}$ written under the orthogonal reduced-form parameterization:
where we let $\bfS_{tr} = h(\bfS)^{\prime}$. Initially, we assume $\bfS_{tr}$ is known, but we will later relax this assumption. Let $\sigma_{tr, ij}$ denote the $i$-th row and $j$-th column entry of $\bfS_{tr}$. For simplicity, we set $\sigma_{tr, 11} = \sigma_{tr, 22} = 1$ and $\sigma_{tr, 21} = -0.9$. Note that the contemporaneous impact matrix ${\bf L}_{0}$ is defined as ${\bf L}_{0} = \bfS_{tr} \operatorname{\mathbf{Q}}$. Henceforth, we focus on the impulse responses to the first shock—it is straightforward to extend our analysis to the second shock.
Given the above, it is easy to see that the impact of the first shock on $y_{t,1}$ and $y_{t,2}$ can be written as $\ell_{11} = q_{11}$ and $\ell_{21} = -0.9 q_{11} + q_{21}$, where $\ell_{i1}$ and $q_{i1}$ are the $(i,1)$ entries of ${\bf L}_{0}$ and $\operatorname{\mathbf{Q}}$, respectively. We now impose sign restrictions requiring that $\ell_{11}$ and $\ell_{21}$ are nonnegative. These sign restrictions imply $q_{11} \geq 0$ and $q_{21} \geq 0.9 q_{11}$. Figure (ref) illustrates this setup graphically. The green circle represents the domain of $\operatorname{\mathbf{Q}}_1 = \left(q_{11}, q_{21}\right)^{\prime}$, while the red arc highlights the identified set that satisfies the imposed sign restrictions.
When using the popular accept-reject sampling approach described in Algorithm (ref), obtaining a draw from the posterior distribution of impulse responses satisfying the sign restrictions involves drawing a $2 \times 1$ vector $\operatorname{\mathbf{x}}_1$ from a $N({\bf 0}, \operatorname{\mathbf{I}}_{2})$ distribution and converting it into a unit vector $\operatorname{\mathbf{q}}_1$ via the normalization $\operatorname{\mathbf{q}}_1 = \operatorname{\mathbf{x}}_1 / ||\operatorname{\mathbf{x}}_1||$. The draw is accepted only if $\operatorname{\mathbf{q}}_1$ satisfies the sign restrictions. Unrestricted draws $\left(q_{11}, q_{21}\right)^{\prime}$ lie uniformly on the entire unit circle (depicted in green), whereas the accepted draws are uniformly distributed only over the subset of the unit circle that meets the sign restrictions (the red arc).
The efficiency of the posterior simulator based on this type of accept-reject algorithm depends heavily on the size of the identified set. As the identified set becomes tighter, we naturally expect to discard a larger number of draws. Indeed, the expected number of draws required to satisfy the sign restrictions is inversely proportional to the probability of meeting those restrictions. Figure (ref) illustrates this relationship analytically: the green line plots the expected number of draws needed to satisfy the sign restrictions as a function of the size of the identified set (i.e., the length of the red arc). More specifically, we generate smaller identified sets by gradually moving the left endpoint of the red arc toward its right endpoint. As shown in the figure, the expected number of draws required increases hyperbolically as the identified set shrinks. In realistic scenarios, as illustrated later in our empirical applications, the number of draws required can become quite large, rendering the algorithm inefficient.
In this paper, we propose an elliptical slice within Gibbs sampler, which draws from the identified set more efficiently. This method can be viewed as an {advanced} Metropolis-Hastings algorithm that transitions from the previous draw $\operatorname{\mathbf{x}}_1^{(0)}$ to a new draw using the following elliptical slice sampling proposal:
where $\boldsymbol{\nu}$ is a $2 \times 1$ vector drawn from $N({\bf 0}, \operatorname{\mathbf{I}}_{2})$. The scalar parameter $\theta$ controls the step size of the proposed move. For instance, when $\theta$ is close to $0$, the proposal is closer to the previous draw $\operatorname{\mathbf{x}}_1^{(0)}$, whereas when $\theta$ approaches $\pi/2$, the proposal is closer to the newly drawn random vector $\boldsymbol{\nu}$. Unlike a conventional Metropolis-Hastings algorithm, the elliptical slice sampling {automatically} searches for a suitable step size to guarantee acceptance of the proposed draw at every iteration. Intuitively, given that the previous draw lies within the identified set, the elliptical slice sampling ensures that the new proposal $\operatorname{\mathbf{x}}_1^{(\star)}$ also remains within the identified set by uniformly drawing $\theta$ from a candidate set that shrinks exponentially. Under appropriate regularity conditions, this procedure ensures the validity and convergence of the algorithm as long as the random variable of interest, in this case $\operatorname{\mathbf{q}}_1$, can be written as a transformation of a normally distributed random variable \citep*[see][for details]{Murray2010,natarovskii2021,hasenpflug2025reversibility}.
The fact that the candidate set for $\theta$ shrinks exponentially is an appealing feature, as it significantly reduces the number of candidate draws of $\operatorname{\mathbf{q}}_1$ needed to satisfy the restrictions. This efficiency gain becomes particularly important as the dimension of the model increases, since generating new draws of $\operatorname{\mathbf{q}}_1$ is computationally costly. Figure (ref) (red line) displays the average number of trials required by the elliptical slice sampling to generate an accepted draw of $\operatorname{\mathbf{q}}_1$ within the identified set as a function of the length of the identified set. The number of required trials for the conventional accept-reject sampler grows hyperbolically, whereas that for elliptical slice sampling increases at a much slower rate.\footnote{To make a fair comparison, one should account for the autocorrelation introduced by {the elliptical slice within Gibbs sampler}, since the accept-reject algorithm generates independent draws. As demonstrated in Section (ref), we compute the effective sample size and find that in this example the number of draws required to obtain one effective draw ranges from 1.04 to 1.35. Therefore, this adjustment does not alter the main conclusion illustrated in the figure.}
In the following section, we extend this simple example into a more realistic and useful setting by: (1) identifying multiple shocks simultaneously rather than just a single shock; (2) allowing sign restrictions to take a general form; and (3) developing an elliptical slice within Gibbs sampler that targets the posterior of the orthogonal reduced-form parameters conditional on the sign restrictions.
In this section, we develop an elliptical slice within Gibbs sampler to obtain draws that approximate the posterior of the orthogonal reduced-form parameters conditional on the sign restrictions defined in Equation (ref). As we will demonstrate in Section (ref), this algorithm can achieve order-of-magnitude gains when compared with the state-of-the-art accept-reject algorithm. For ease of exposition, we organize the section into four parts. Section (ref) establishes a mapping between the posterior distribution of interest and its implied posterior over a transformed random vector of matrices (defined below) that forms the basis of our algorithm. Section (ref) describes the design of our sampler. Section (ref) describes the algorithm. Finally, Section (ref) provides an argument for the well-definedness of our elliptical slice within Gibbs sampler.
Let $\operatorname{\mathbf{Z}}=(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}_{\bfS}},\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})\in \operatorname{\mathcal{Z}}$, with $\operatorname{\mathcal{Z}} = \mathbb{R}^{m\times n} \times \mathbb{R}^{n\times \tilde{\nu}} \times \mathbb{R}^{n \times n}$, and let $p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}}\mid ({\mathbf{y}}_t)_{t=1}^{T}, \mathbf{S}_R(T(\operatorname{\mathbf{Z}})) > {\bf0})$ denote the following conditional Lebesgue density:
where
and $T(\operatorname{\mathbf{Z}})= (\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}},\zeta(\operatorname{\mathbf{Z}_{\bfS}}),\gamma(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}))$ is a transformed random vector in which $\gamma(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$ denotes the orthogonal matrix obtained from the QR decomposition of $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}$ normalized so that the diagonal of the associated upper triangular matrix is positive, and $\zeta(\operatorname{\mathbf{Z}_{\bfS}}) = (\operatorname{\mathbf{Z}_{\bfS}} \operatorname{\mathbf{Z}_{\bfS}}^{\prime} )^{-1}$ denotes a positive definite matrix.
Considering the transformed random matrices $\zeta(\operatorname{\mathbf{Z}_{\bfS}})$ and $\gamma(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$, the measurable sets $\operatorname{{\mathcal{A}}}\subseteq \mathbb{R}^{m\times n}$, $\operatorname{{\mathcal{B}}}\subseteq \mathcal{S}(n)$, and $\operatorname{{\mathcal{C}}}\subseteq \mathcal{O}(n)$, and defining $(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})= T(\operatorname{\mathbf{Z}})$, we have
where
denote the push-forward measures associated with the mappings $\zeta \colon \mathbb{R}^{n\times \tilde{\nu}} \to \mathcal{S}(n)$ and $\gamma\colon \mathbb{R}^{n\times n} \to\mathcal{O}(n)$.\footnote{The push-forward characterization is introduced only for notational consistency. The explicit forms of these measures follow immediately from Theorem 2 in arrw2018: $\zeta$ yields an inverse-Wishart distribution and $\gamma$ yields the Haar measure on $\mathcal{O}(n)$.} As a consequence, the density of the random vector $T(\operatorname{\mathbf{Z}})$ associated with Equation (ref) coincides with the density in Equation (ref), that is, the posterior distribution of interest. This implies that we can draw from Equation (ref) and transform the draws using the mapping $T$ to obtain approximate draws from Equation (ref).
The sampling strategy has at its core an elliptical slice within Gibbs sampler that generates approximate draws from the target density $p_{\operatorname{\mathbf{Z}}}$. More specifically, the sampler produces draws of a Markov chain $(\operatorname{\mathbf{Z}}^{i})_{i\in\mathbb{N}}$ whose {single random variable distribution} $Law(\operatorname{\mathbf{Z}}^{i})$ is close to the target law $Law(\operatorname{\mathbf{Z}})$ induced by $p_{\operatorname{\mathbf{Z}}}$. Our ultimate object of interest, however, is the posterior distribution of $(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})$ with density $p$ in Equation (ref). Hence, to accomplish our goal we obtain draws in this parameterization of interest via the transformation $T(\operatorname{\mathbf{Z}}^{i})$. For now, we {justify this approach by assessing} the approximation error using the total variation distance between the target $Law(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})$ and the induced $Law(T(\operatorname{\mathbf{Z}}^{i}))$. Since $(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})=T(\operatorname{\mathbf{Z}})$, $Law(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})$ is the push-forward of $Law(\operatorname{\mathbf{Z}})$ under $T$, and similarly for $Law(T(\operatorname{\mathbf{Z}}^{i}))$ and $Law(\operatorname{\mathbf{Z}}^{i})$. A standard contraction property of total variation under measurable maps implies
Therefore, whenever the right-hand side is small, the transformed draws $T(\operatorname{\mathbf{Z}}^{i})$ approximate the target posterior $Law(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})$ in total variation distance.
This section describes the aforementioned elliptical slice within Gibbs sampler. We {begin by describing} the three building blocks of the approach: the conditional posterior $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}$ given $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}_{\bfS}})$, the conditional posterior for $\operatorname{\mathbf{Z}_{\bfS}}$ given $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$, and the conditional posterior for $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}$ given $(\operatorname{\mathbf{Z}_{\bfS}},\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$. Then we combine these posteriors into a single algorithm.
\paragraph{Conditional Posterior for $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}$.} We first derive the posterior for $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}$ conditional on $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\bfS}}$, and the sign restrictions. Using Bayes’ rule and independence of $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\bfS}})$ and $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}$, Equation (ref) implies:
The right-hand side determines a well-defined conditional Lebesgue density. To see this, define the set
and notice that by the continuity of $\zeta$, $\gamma$ and $\mathbf{S}_R$ the set $\mathcal{P}$ is open. Consequently, for any fixed $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\bfS}})$, the set
is open in $\mathbb{R}^{n \times n}$.
\paragraph{Conditional Posterior for $\operatorname{\mathbf{Z}_{\bfS}}$.} We next derive the posterior of $\operatorname{\mathbf{Z}_{\bfS}}$ conditional on $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$ and the sign restrictions:
As before, the proportionality argument (well-definedness of the denominator) follows from the fact that, for any fixed $(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$, the set
is open in $\mathbb{R}^{n\times \tilde{\nu}}$.
\paragraph{Conditional Posterior for $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}$.} Finally, the posterior of $\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}$ conditional on $(\operatorname{\mathbf{Z}_{\bfS}}, \operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$ and the sign restrictions is:
Finally, the proportionality argument in this case follows from the fact that, for any fixed $(\operatorname{\mathbf{Z}_{\bfS}}, \operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$, the set
is open in $\mathbb{R}^{m \times n}$.
\paragraph{Algorithm.} Crucially, each of the conditional densities defined above shares the feature of being proportional to the product of a multivariate Gaussian and an arbitrary function that we can evaluate. This is exploited in our elliptical slice sampling-within-Gibbs approach. We derive a Markov chain $(\operatorname{\mathbf{Z}}^i)_{i\in\mathbb{N}}=(\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}^{i}_{\bfS}},\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}})_{i\in\mathbb{N}}$ and the transformed chain $(\operatorname{\mathbf{B}}^i,\bfS^i,\operatorname{\mathbf{Q}}^i)_{i\in\mathbb{N}}=(T(\operatorname{\mathbf{Z}}^i))_{i\in\mathbb{N}}$ for targeting the distribution of interest determined by Equation (ref).
As mentioned above, we have used the conjugate normal-inverse-Wishart prior over the reduced-form parameters to describe the algorithm. Appendix (ref) shows that the approach can be easily adapted to two alternative priors: the independent normal-inverse-Wishart prior and the asymmetric conjugate priors of Chanassymetric.
As in other MCMC algorithms, our algorithm must start from parameter values that satisfy the sign restrictions. In particular, the elliptical slice sampling in Step 2 must be initialized at values $\operatorname{\mathbf{Z}}^{0} \in \mathcal{P}$. Depending on the sign restrictions imposed, finding such parameter values to initialize the algorithm can be computationally costly. Therefore, we develop a novel and robust method that combines the i.i.d.\ sampler with our proposed Gibbs sampler. Appendix (ref) describes the procedure in detail.
In this section, we justify the well-definedness of Algorithm (ref) by establishing stationarity of the elliptical slice sampling-within-Gibbs transition mechanism with respect to the desired distribution $p_{\operatorname{\mathbf{Z}}}$ given in Equation (ref). To motivate the argument, suppose that $\operatorname{\mathbf{Z}}^{i-1}$ follows the desired posterior distribution. We then want to construct the transition to $\operatorname{\mathbf{Z}}^{i}$ so that this distribution is preserved exactly. If so, $\operatorname{\mathbf{Z}}^{i}$ also follows the posterior of interest. Lemma (ref) formalizes the intuition provided above.
We now illustrate the performance of our algorithm using two empirical applications. The first is a small-scale SVAR of the global oil market, based on the model in \citet*{kilianmurphy2014}, which identifies flow supply, flow demand, and speculative demand shocks using a combination of sign and elasticity bounds. The tight identifying assumptions in this model render traditional accept-reject methods computationally intensive, whereas our algorithm improves efficiency while replicating the main results. The second application revisits the large-scale SVAR model of the U.S. economy developed by \citet*{crump2025large} and analyzed structurally by \citet*{chan2025largestructuralvarsmultiple}, which includes 35 macroeconomic and financial variables and identifies up to eight structural shocks. We show that our algorithm remains computationally stable as the number of restrictions increases, in contrast to the exponential rise in computation time exhibited by the accept-reject method. Both applications highlight the scalability of our approach in distinct empirical settings. For each application, we first demonstrate that our approach replicates the main results reported in the original papers and then analyze the computational timing to show that our method can be more efficient than the traditional accept-reject algorithm.
In our first application, we replicate the results of kilianmurphy2014, who extend the framework of kilian2012agnostic by incorporating oil inventories to identify speculative demand shocks. The identification strategy in kilianmurphy2014 relies on tight sign and elasticity bound restrictions, which result in a small identified set and may render standard accept-reject algorithms slow. Notably, kilianmurphy2014 adopt an approach similar to that of \citet*{chan2025largestructuralvarsmultiple}, relying on permutations and sign alternations. Therefore, when using the same set of sign and elasticity bound restrictions, the computation times we report for the accept-reject algorithm are comparable to those in the original study and similar to those obtained with our {Gibbs-based} approach. However, when an additional elasticity bound is introduced, the accept-reject approach becomes computationally unfeasible, while the Gibbs sampler remains as fast, highlighting the speed advantages of our algorithm in such settings.
We begin by describing the model specification in \citet*{kilianmurphy2014}. They model the global market for crude oil using a four-variable SVAR featuring the percent change in global crude oil production, a measure of global real activity, the real price of crude oil, and the change in global above-ground crude oil inventories. The SVAR is specified at a monthly frequency, with an estimation sample covering 1973M2--2009M8. The model includes 24 lags, a constant, and seasonal dummies to remove seasonal variation. \citet*{kilianmurphy2014} adopt a weak conjugate normal-inverse-Wishart prior distribution \citep*[see, e.g.,][]{uhlig2005effects} for the reduced-form parameters.
Turning to identification, the goal of \citet*{kilianmurphy2014} is to identify three structural shocks using a combination of sign restrictions on impact impulse responses, sign restrictions at horizons 1 through 12, and elasticity bounds. Table (ref) summarizes the identifying assumptions. The structural shocks are labeled flow supply shock, flow demand shock, and speculative demand shock.
We adopt the exact same specification and use Algorithm (ref). We obtain one million draws, saving one every 10; hence, the figures are produced using one hundred thousand draws. Figure (ref) presents the impulse responses to the three shocks implied by the Gibbs sampler and the conventional accept-reject approach in RRWZ2010 with the efficiency gains proposed by chan2025largestructuralvarsmultiple. {As can be seen,} the results are nearly identical and they are in line with those in \citet*{kilianmurphy2014}. In particular, a negative flow supply shock causes a persistent decline in global economic activity and oil inventories, and a persistent increase in the real price of oil. The response of oil production is persistently negative. A positive flow demand shock is associated with a persistent increase in global economic activity, a persistent increase in the real price of oil, and a positive response of oil production. Oil production increases sluggishly, given the imposed elasticity bounds, and peaks at about one year after the shock before declining to pre-shock levels. Finally, a positive speculative demand shock causes a persistent increase in the real price of oil and a large increase in inventories. Global real activity and oil production decline persistently in response to this shock, although the effects are modest.
We next compare the computational time of our Gibbs sampler to that of the accept-reject algorithm. Table (ref) reports the time (in hours) per 1,000 effective draws---defined as total computation time divided by the effective sample size and scaled by 1,000---using Algorithm (ref) and the accept-reject method.\footnote{Since the accept-reject approach produces independent draws, the number of effective draws equals the total number of draws in that case.} As highlighted above, we implement the efficient variant of the accept-reject algorithm proposed by chan2025largestructuralvarsmultiple. We approximate the effective sample size using the multivariate effective sample size metric of \citet*{vats2019multivariate}, which is well suited for SVAR analysis where inference often targets high-dimensional objects such as vectors of impulse responses. In particular, we estimate effective sample size under the impulse response parameterization using the multivariate batch means approach described in that work, with a batch size of $N^{1/4}$, where $N$ denotes the number of stored draws.\footnote{We consider the first three columns of the impulse response parameters, which correspond to the shocks of interest in this application. Results are robust to using impulse responses for all shocks.} Because the accept-reject algorithm produces independent draws, effective draws and sample size are the same.
The column “Benchmark Model” in Table (ref) compares the time (in hours) required to obtain 1,000 effective draws under the specification of \citet*{kilianmurphy2014}. As shown in the table, the Gibbs sampler requires less than 2 minutes to produce 1,000 effective draws. In contrast, the accept-reject algorithm takes approximately 20 minutes to achieve the same number of draws.\footnote{All computations were performed in MATLAB on an Intel Xeon Platinum 8488C processor with 16 active cores running at $2.4$ GHz on an x86_64 architecture.} To further illustrate the gains of our proposed Gibbs sampler, we consider a scenario in which a researcher imposes an additional restriction on the price elasticity of oil demand in response to a flow supply shock. As emphasized by \citet*{CALDARA20191}, such a restriction is empirically important in SVAR models of the oil market. Following their work, we constrain the price elasticity of crude oil demand to lie within a narrow interval around the point estimate of $-0.08$ reported by \citet*{CALDARA20191}; specifically, we impose the restriction that the elasticity must lie in the interval $(-0.09,-0.07)$. The column “Benchmark Model + Additional Restriction” in Table (ref) reports the results under this added constraint. The time under the Gibbs sampler is under 6 minutes, while the performance of the accept-reject algorithm deteriorates sharply, requiring nearly 8 hours to obtain 1,000 draws.
These results indicate that, for this model, the accept-reject algorithm is already near its computational limit under the benchmark specification, and that introducing even a single additional restriction dramatically reduces its efficiency. In contrast, the Gibbs sampler maintains its performance even as additional identifying restrictions are introduced.
In our second application, we replicate and extend the analysis of \citet*{chan2025largestructuralvarsmultiple}, who build on the large-scale SVAR framework of \citet*{crump2025large} to study the structural dynamics of the U.S. economy. Their model incorporates 35 macroeconomic and financial variables commonly monitored by the Federal Reserve and identifies eight structural shocks using an extensive set of sign and ranking restrictions. \citet*{chan2025largestructuralvarsmultiple} employ an accept-reject algorithm, which becomes computationally intensive as the number of identifying restrictions increases.\footnote{We thank Christian Matthes for sharing their replication files and data with us.} As we show below, our algorithm is more efficient.
To demonstrate this, we extend the baseline model by identifying two additional shocks---an oil price shock and a consumer sentiment shock---bringing the total number of sign restrictions from 105 to 129. This provides a stringent test of our algorithm's performance relative to the accept-reject method. Importantly, \citet*{chan2025largestructuralvarsmultiple} use the asymmetric priors defined in \citet*{Chanassymetric} for the reduced-form parameters, instead of the Minnesota prior used in \citet*{crump2025large}. To simplify the comparison, we revert to the Minnesota prior when applying both the accept-reject algorithm and the Gibbs sampler. Online Appendix (ref) compares both approaches under the asymmetric prior.
The SVAR used in this section is specified at a quarterly frequency, includes a constant and five lags, and uses an estimation sample that spans from 1977Q4 to 2019Q4.\footnote{The starting date corresponds to the first date with data availability for all variables after taking into account the initial observations required by the specification.} As mentioned, we assume a Minnesota prior for the reduced-form parameters and set the hyperparameters following \citet*{giannone2015prior}. Turning to the identification, Table (ref) summarizes the variables and the sign and ranking restrictions imposed on the contemporaneous impulse responses.
A table entry equal to $+1$ indicates a positive sign restriction, a $-1$ indicates a negative one, and $*$ indicates that the corresponding contemporaneous impulse response has been left unrestricted. \citet*{chan2025largestructuralvarsmultiple} consider only the first eight shocks (demand, investment, financial, monetary, government spending, technology, labor supply, and wage bargaining). In total, 105 sign restrictions are imposed in their baseline specification. We have added two additional shocks (labeled oil price and consumer sentiment) to assess the performance of our algorithm. With the inclusion of these two shocks, the total number of sign restrictions increases to 129. When using the Gibbs sampler, we obtain one million draws and retain one every 10.
Figures (ref) through (ref) report the impulse responses for each of the eight structural shocks identified by \citet*{chan2025largestructuralvarsmultiple} obtained using the Gibbs sampler described in Algorithm (ref). The solid red lines and the orange-shaded areas depict the point-wise posterior median and 68 percent probability bands, respectively. For comparison, we also report the point-wise posterior median (gray dashed lines) and point-wise 68% probability bands (gray-shaded area) obtained with chan2025largestructuralvarsmultiple's chan2025largestructuralvarsmultiple accept-reject algorithm.\footnote{When using their algorithm we use 1,000 draws to report the posterior impulse responses.} The results are virtually identical as expected given that both algorithms are sampling from the same posterior distribution. Let us begin by describing the selected impulse responses to a unit standard deviation expansionary demand shock, shown in Figure (ref). The signs of the impact responses of real GDP, the PCE price index, the federal funds rate, and the unemployment rate are restricted. The remaining horizons, as well as the responses of non-residential investment and the real wage, are unrestricted. As can be seen, the demand shock causes a transient increase in output and prices, and a decrease in the unemployment rate. The federal funds rate increases in response to the shock. The restrictive stance of monetary policy eventually lowers economic activity, as seen, for example, in the decline of non-residential investment. The real wage decreases in the short run in response to the shock, as nominal wage increases are not sufficient to offset higher prices---possibly due to sluggish nominal wage adjustment.
The investment shock, shown in Figure (ref), appears similar to the demand shock in terms of economic consequences for real GDP, the federal funds rate, the price level, and the unemployment rate. However, the impulse response of non-residential investment is substantially different. In particular, the investment shock causes a short-run boom in non-residential investment. This finding partly reflects the ranking restriction requiring that the impact response of non-residential investment be larger than the impact response of real GDP. As with the demand shock, the investment shock causes a persistently negative response of the real wage.
Turning to the financial shock, shown in Figure (ref), it is worth highlighting that this shock is identified using the same sign restrictions as the investment shock, except for the impact response of the S&P 500, which is assumed to be positive instead of negative. Overall, the impulse responses are similar, except that the decline in non-residential investment after five quarters is slightly less pronounced under the financial shock, consistent with the positive response of asset prices.
The impulse responses to a unit standard deviation contractionary monetary policy shock are depicted in Figure (ref). This shock causes the federal funds rate to remain above zero for more than two years, reflecting inertia in the conduct of monetary policy. Real GDP and prices decline persistently, and the unemployment rate jumps upon impact before slowly returning to baseline. Non-residential investment drops on impact and recovers after about one year, in line with a less restrictive monetary policy stance. The real wage increases, driven by a decrease in the price level. A notable aspect of these responses is that they suggest monetary policy can operate with shorter lags than traditionally assumed under the “long and variable lags” view.
The government spending shock is shown in Figure (ref). An expansionary one unit standard deviation government spending shock leads to an increase in real GDP for about two quarters and to a long-lasting increase in the price level. To conclude, we discuss the impulse responses to the supply-related structural shocks, that is, the technology, labor supply, and wage bargaining shocks. A unit standard deviation positive technology shock leads to a protracted increase in real GDP and non-residential investment (see Figure (ref)). The higher level of output is accompanied by a sustained decline in the unemployment rate and a sustained increase in the real wage. The federal funds rate rises marginally, indicating that monetary policy remains roughly neutral in response to technology shocks.
The responses to a unit standard deviation positive labor supply shock are shown in Figure (ref). This shock induces a hump-shaped response of real GDP and leads to persistently lower prices. The responses to a unit standard deviation negative wage bargaining shock are shown in Figure (ref). The identifying assumptions for this shock are identical to those of an expansionary labor supply shock, except that the unemployment rate is assumed to decrease upon impact. When a negative wage bargaining shock occurs, workers experience a decline in their nominal wage alongside a decrease in the unemployment rate. The real wage remains unaffected on impact, as the lower wages are offset by the assumed decrease in the price level. Subsequently, the price level remains below zero, inducing an increase in the real wage.
We begin by comparing the efficiency of the Gibbs sampler algorithm relative to the accept-reject algorithm when replicating the identification scheme in \citet*{chan2025largestructuralvarsmultiple}. Figure (ref) reports the time (in minutes) per 1,000 effective draws using Algorithm (ref) as a function of the number of identified shocks.\footnote{When computing the multivariate effective sample size, we only consider the columns of the impulse response parameters corresponding to the shocks of interest in this application.}
To assess the computational time as a function of the size of the identified set, we proceed incrementally: we first obtain draws by identifying only the demand shock, then add the investment shock, the financial shock, and so on, until all eight shocks in Table (ref) are included. As shown, the time per 1,000 effective draws remains computationally feasible even as the number of sign restrictions increases.
Figure (ref) replicates the same figure but using the efficient accept-reject version of Algorithm (ref) proposed by \citet*{chan2025largestructuralvarsmultiple}. In this case, the computation time increases dramatically, as shown in the figure. Although the runtime will vary depending on the hardware architecture and the number of variables, the main conclusion from comparing Figures (ref)–(ref) remains unchanged: the performance of the accept-reject algorithm can deteriorate sharply as the identified set narrows. Figure (ref) combines the timings to facilitate visual comparison.
To further emphasize this point, we now consider additional shocks to illustrate that the accept-reject approach can eventually become impractical. Specifically, we extend the number of shocks identified in \citet*{chan2025largestructuralvarsmultiple} by adding the oil price shock and the consumer sentiment shock described in Table (ref). Figures (ref)–(ref) replicate the exercise shown in Figures (ref)–(ref) for the cases of nine and ten shocks.\footnote{The runtime of the accept-reject algorithm for the nine- and ten-shock cases is extrapolated based on ten draws.}
As the reader can see, computation time does not increase exponentially when we use the elliptical slice sampling-within-Gibbs approach. In contrast, when we consider nine shocks under the accept-reject approach, the times are now measured in hours, and when we consider ten shocks, the times are measured in days---Figures (ref) and (ref) provide the detailed timings. To facilitate the comparison, Figure (ref) overlays both sets of timings to make clear that our algorithm can handle settings (in terms of the number of variables and shocks) that the traditional accept-reject approach cannot.
The main computational cost of our approach stems from running the Gibbs sampler, particularly in large models, since it produces {autocorrelated} draws. As a result, one might be tempted to bypass this cost and instead address the bottleneck issues inherent in the accept-reject approach by employing the conditionally uniform prior approach described in \citet*{uhlig2017shocks}, \citet*{AmirAhmadiDrautzburg}, and \citet*{readzhu2025}, among others. The main appeal of this simpler approach is that, like the accept-reject algorithm, it generally yields independent draws. While such a simplification is indeed attractive due to its lower computational burden, it is essential for the researcher to be aware of two critical drawbacks, which we will explain in this section. Before turning to these pitfalls, we first describe the conditionally uniform prior and outline the main algorithmic steps commonly used in its implementation.
Let $\mathbb{Q}_n(\operatorname{\mathbf{B}},\bfS)=\{\operatorname{\mathbf{Q}} \in \mathcal{O}(n): \mathbf{S}_R(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})> {\bf0}\}$ and define $\kappa(\operatorname{\mathbf{B}},\bfS)$ implicitly by $\int_{\mathbb{Q}_n(\operatorname{\mathbf{B}},\bfS)} \kappa(\operatorname{\mathbf{B}},\bfS) d\operatorname{\mathbf{Q}} = 1$. Unlike the case of the uniform prior, this constant of integration depends on $(\operatorname{\mathbf{B}},\bfS)$. The prior underlying the conditionally uniform approach is then given by:
Crucially, $\kappa(\operatorname{\mathbf{B}},\bfS)$ depends on $(\operatorname{\mathbf{B}},\bfS)$, while $\kappa$ in Equation (ref) does not. This is due to the fact that the conditionally uniform approach combines the conjugate normal-inverse-Wishart prior over the reduced-form parameters with $\pi(\operatorname{\mathbf{Q}} \mid \operatorname{\mathbf{B}}, \bfS)$ of the form:
This conditional uniform prior has the property that it overweights reduced-form parameters with smaller identified sets, measured by $\kappa^{-1}(\operatorname{\mathbf{B}}, \bfS)$ \citep*[see][]{uhlig2017shocks}. Unlike the prior in Equation (ref), the prior in Equation (ref) cannot be justified using the results of arias2025uniform, who aim to construct priors that, among other things, separate inference from identification, as it is undesirable for the prior to change when the restrictions are modified, since this makes it impossible to determine whether differences in results stem from changes in the prior or from changes in the identification restrictions.
Under the conditional uniform prior, the objective is to draw from the following posterior of the orthogonal reduced-form parameters conditional on the sign restrictions:
where
and then use $f$ and $\phi$ to transform the draws to the desired impulse responses. It is straightforward to adapt the traditional approach described in Algorithm (ref) to obtain draws from Equation (ref) as follows:
While Algorithm (ref) does not sample from the posterior distribution defined in Equation (ref), it draws from the posterior distribution defined in Equation (ref) and it can be justified under a prior different from that in Equation (ref).
The first issue with Algorithm (ref) is that Step 2 may {never} succeed for some draws of $(\operatorname{\mathbf{B}},\bfS)$. Recall $\mathbb{Q}_n(\operatorname{\mathbf{B}},\bfS)=\{\operatorname{\mathbf{Q}}\in\mathcal{O}(n):\mathbf{S}_R(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})> {\bf0}\}$. If, for a given $(\operatorname{\mathbf{B}},\bfS)$, $\mathbb{Q}_n(\operatorname{\mathbf{B}},\bfS)=\varnothing$, then drawing $\operatorname{\mathbf{Q}}$ uniformly over $\mathcal{O}(n)$ until $\mathbf{S}_R(\operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}})> {\bf0}$ results in non-termination. This situation is more pervasive than one might expect. For example, requiring the impulse response of a given variable to a particular shock to maintain the same sign for more than one period rules out sign-flipping impulse-response patterns (i.e., responses that change sign every period). Consequently, any $\operatorname{\mathbf{B}}$ that implies such every-period oscillations in the impulse responses (e.g., via complex eigenvalues of the companion-form VAR coefficient matrix) will be associated with an empty identified set.
In contrast, our algorithm conditions on the open set $\mathcal{P}$, defined in Section (ref). As a consequence, by construction it cannot produce draws outside {that set}.
The second issue with the conditionally uniform prior approach is that inference and identification become intertwined. Consider an SVAR with $n = 3$ and $m = 0$, that is, without lags or constant terms. Suppose that Researcher A aims to identify three structural shocks using identification scheme A, as defined in Table (ref), while Researcher B employs identification scheme B, also defined in the same table. Clearly, any set of impulse responses satisfying scheme B will also satisfy scheme A.
Define $\mathbb{Q}_n^j(\bfS)=\{\operatorname{\mathbf{Q}} \in \mathcal{O}(n): \mathbf{S}^j_R(\bfS,\operatorname{\mathbf{Q}})> {\bf0}\}$ as {the set of orthogonal matrices} consistent with the identification restrictions of researcher $j \in \{A,B\}$, with associated $\kappa^j(\bfS)$ implicitly defined by $\int_{\mathbb{Q}_n^j(\bfS)} \kappa^j(\bfS) d\operatorname{\mathbf{Q}} = 1$. Even if both researchers specify the same reduced-form prior, changing the identification scheme alters the implied prior over the impulse responses. To illustrate, set $\nu=100$ and $\bfPhi=\operatorname{\mathbf{I}}_n$, and consider ten values of $\{\bfS^i\}_{i=1}^{10}$ such that when evaluated at the prior distribution $IW_{(\nu,\bfPhi)}$ they all have identical prior density. Let $\{\operatorname{\mathbf{L}}_0^i\}_{i=1}^{10}$ denote ten corresponding impact impulse response matrices {that happen to satisfy identification B}. Because the volume element from $(\bfS,\operatorname{\mathbf{Q}})$ to $\operatorname{\mathbf{L}}_0$ only depends on the determinant of $\bfS$ and we have also restricted $\{\bfS^i\}_{i=1}^{10}$ to have the same determinant, all $\{\operatorname{\mathbf{L}}_0^i\}_{i=1}^{10}$ are equally favored by the prior before any sign restrictions are imposed. After the sign restrictions are introduced, we have:
where $\pi^j(\operatorname{\mathbf{L}}_0)$ is the prior density over impulse responses under researcher $j \in \{A,B\}$. This expression implies that any variation in this ratio is solely attributable to the additional negative sign restriction on shock 2 to variable 3 associated with identification scheme B. As shown in Table (ref), the difference in the identification schemes leads to differences in the implied priors. On the one hand, Researcher A’s prior favors $\operatorname{\mathbf{L}}^5_0$ $1.45$ times as much as $\operatorname{\mathbf{L}}^{1}_{0}$, while Researcher B’s prior favors $\operatorname{\mathbf{L}}^1_0$ about twice as much as $\operatorname{\mathbf{L}}^4_0$. On the other hand, Researcher B’s prior favors $\operatorname{\mathbf{L}}^3_0$ twice as much as $\operatorname{\mathbf{L}}^{1}_{0}$, while Researcher A’s prior favors $\operatorname{\mathbf{L}}^1_0$ $1.13$ times as much as $\operatorname{\mathbf{L}}^3_0$. These findings demonstrate that, under the conditionally uniform approach, changing identification schemes alters the implied prior over parameters of interest (such as impulse responses), thereby entangling estimation and identification.
The reason for this unfortunate result is that $\kappa(\bfS)$ varies across identification schemes in a way that disproportionately favors $\operatorname{\mathbf{L}}_{0}$ values associated with $\bfS$ that induce smaller identified sets. Since different sign restrictions affect the size of the identified set differently for each $\bfS$, different identification restrictions will imply different priors over impulse response functions. This problem does not arise under the uniform prior described in Equation (ref), since in that case $\kappa$ does not depend on the reduced-form parameters. More broadly, this highlights the cost of not adopting a uniform prior over the orthogonal matrices as described in \citet*{arias2025uniform}. As shown in that paper, specifying a uniform prior over the set of orthogonal matrices {ensures that inference is disentangled from identification}.
This paper proposes an elliptical slice within Gibbs sampler for inference based on structural vector autoregressions identified with sign restrictions. We show that the algorithm effectively overcomes the computational bottlenecks associated with conventional accept-reject methods, especially as the number of identifying restrictions increases or as the identified set becomes tight. Our empirical applications illustrate how the proposed algorithm can extend existing analyses in the literature, including SVARs with a {large number} of macroeconomic and financial variables. Overall, the paper provides contributions to the implementation of sign-restricted SVARs, offering tools that are applicable across a wide range of empirical models.
\addappheadtotoc \setcounter{section}{0} \setcounter{equation}{0}
{ \bf{Appendix}}
Let us begin by noting that given $\operatorname{\mathbf{Z}}^{i-1}$, we obtain $\operatorname{\mathbf{Z}}^{i}$ by crucially drawing from the elliptical slice sampling in Steps (ref)-(ref) of Algorithm (ref). Recall that the transition kernels associated with each of the elliptical slice samplers are $P_{\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}^{i-1}_{\bfS}}}( \operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{Q}}}},d\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}})$, $P_{\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}}}(\operatorname{\mathbf{Z}^{i-1}_{\bfS}},d\operatorname{\mathbf{Z}^{i}_{\bfS}})$, and $P_{\operatorname{\mathbf{Z}^{i}_{\bfS}}, \operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}}}(\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}},d\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{B}}}})$, respectively. Our proof of stationarity builds on reversibility and on three block-conditional representations of Equation (ref). We describe these properties in Sections (ref)-(ref), and then we use them to establish stationarity in Section (ref). For notational convenience in this proof, sometimes we use measure notation and other times we use density notation.
Importantly, these {kernels} satisfy by construction a reversibility property, cf. hasenpflug2025reversibility\footnote{In hasenpflug2025reversibility Theorem 3.1 is formulated for densities that are strictly positive. The issue of applying the result to densities that may attain the value zero can be resolved by modifying the objects and arguments of hasenpflug2025reversibility as in Durmus2026geodesic.}. For the latter result to be applicable we exploit the facts, provided in Section (ref), that the sets $\mathcal{P}(\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{B}}},\operatorname{\mathbf{Z}}_{\bfS})$, $\mathcal{P}(\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{B}}},\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{Q}}})$ and $\mathcal{P}(\operatorname{\mathbf{Z}}_{\bfS},\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{Q}}})$ are open and therefore yield that the conditional posterior densities for $\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{Q}}}$, $\operatorname{\mathbf{Z}}_{\bfS}$ and $\operatorname{\mathbf{Z}}_{\operatorname{\mathbf{B}}}$ are lower semi-continuous. This is essential to guarantee that the shrinkage procedure used within elliptical slice sampling terminates almost surely, cf. hasenpflug2025reversibility. The aforementioned reversibility property implies the following crucial stationarity features with respect to the appearing conditionals:
where we shorten the notation using a placeholder $\Delta\in \{ \operatorname{\mathbf{B}},\bfS,\operatorname{\mathbf{Q}}\}$ to write $H^{\Delta}_{\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\bfS}},\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}}$ for the “conditioning expression”
$H_{\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}_{\bfS}},\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}}}$ for $({\mathbf{y}}_t)_{t=1}^{T}, \mathbf{S}_R(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}, \zeta(\operatorname{\mathbf{Z}_{\bfS}}), \gamma(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})) > {\bf0}$, and
In addition to Equations (ref)-(ref), we also need the following relations regarding the different conditional densities. The next displayed expression rewrites the same joint density in three equivalent block-conditional forms, one for each Gibbs step.
where
Given that $\operatorname{\mathbf{Z}}^{i-1}$ is distributed according to $p_{\operatorname{\mathbf{Z}}}$, then the distribution of $\operatorname{\mathbf{Z}}^{i}$ is determined by the probability of $(\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}^{i}_{\bfS}},\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}}) \in \mathcal{A} \times \mathcal{B} \times \mathcal{C}$, for arbitrary measurable sets $\mathcal{A}\subseteq \mathbb{R}^{m \times n}$, $\mathcal{B} \subseteq \mathbb{R}^{n \times \tilde{\nu}}$, $\mathcal{C} \subseteq \mathbb{R}^{n \times n}$, given by
with
Combining the two equations above, we obtain:
Subsequently, we will exploit the equations described in Sections (ref)-(ref) along with multiple changes of the order of integration to demonstrate stationarity. Let us begin by substituting Equation (ref) for $p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}}^{i-1} \mid H_{\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}}, \operatorname{\mathbf{Z}^{i-1}_{\bfS}},\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{Q}}}}})$ in order to obtain:
Then, by Equation (ref) we obtain:
From Equation (ref), it follows that:
By Equation (ref), we can substitute out $\int_{\mathbb{R}^{n \times \tilde{\nu}}} P_{\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}},\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}}} (\operatorname{\mathbf{Z}^{i-1}_{\bfS}}, d \operatorname{\mathbf{Z}^{i}_{\bfS}})P_{\operatorname{\mathbf{Z}^{i}_{\bfS}},\operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{Q}}}}} (\operatorname{\mathbf{Z}^{i-1}_{\operatorname{\mathbf{B}}}}, d \operatorname{\mathbf{Z}^{i}_{\operatorname{\mathbf{B}}}})$:
Using Equation (ref), we get:
Applying Equation (ref), it follows that:
Finally, Equations (ref), (ref) and (ref) together imply that we can write:
which is equal to:
Thus, we have shown that:
From this we can deduce that the distribution after a single transition from $\operatorname{\mathbf{Z}}^{i-1}$ to $\operatorname{\mathbf{Z}}^{i}$ does not change; hence, the chain is stationary. This proves Lemma (ref).
As discussed in Section (ref), our algorithm must start from parameter values that satisfy the sign restrictions. Depending on the sign restrictions imposed, finding such parameter values to initialize the algorithm can be computationally costly. For example, in our applications in Section (ref), it can take up to 2 hours to find valid initial values using the state-of-the-art accept-reject sampler developed by chan2025largestructuralvarsmultiple. Hence, it is desirable to have a robust and more efficient method to initialize the algorithm.
To achieve this goal, inspired by adaptive Sequential Monte Carlo (SMC), we begin with draws (particles) from the unrestricted posterior and move them toward the posterior conditional on the sign restrictions.
More specifically, consider a sequence of posterior distributions indexed by $n \in \left\lbrace 1,\dots,N \right\rbrace$:
where $p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}} \mid ({\mathbf{y}}_t)_{t=1}^{T})$ is the unrestricted posterior density of $\operatorname{\mathbf{Z}}$ in the transformed space and $\{\operatorname{\mathbf{c}}_{n}\}^N_{n=1}$ is a sequence of threshold vectors with ${\bf{-\infty}} = {\bf{c}}_{1} < {\bf{c}}_{2} < \ldots < {\bf{c}}_{N} \leq {\bf{0}}$. For instance, under the conjugate normal-inverse-Wishart prior, $p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}} \mid ({\mathbf{y}}_t)_{t=1}^{T}) = N_{(\tilde\bfPsi, \tilde\bfOme, \zeta(\operatorname{\mathbf{Z}_{\bfS}}))}(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{B}}}}) \, N_{(\mathbf{0}_{n \times \tilde{\nu}}, \tilde{\bfPhi}^{-1}, \operatorname{\mathbf{I}}_{\tilde{\nu}})}(\operatorname{\mathbf{Z}_{\bfS}}) \, N_{(\mathbf{0}_{n \times n}, \operatorname{\mathbf{I}}_n, \operatorname{\mathbf{I}}_n)}(\operatorname{\mathbf{Z}_{\operatorname{\mathbf{Q}}}})$ as defined in Section (ref), Equation (ref). Under alternative priors, this density can be adjusted accordingly; see Appendix (ref). When ${\bf{c}}_{1} = {\bf{-\infty}}$, the indicator is always one, so that $\pi_{1}(\operatorname{\mathbf{Z}})=p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}} \mid ({\mathbf{y}}_t)_{t=1}^{T})$. As ${\bf{c}}_{n}$ increases, the restrictions tighten and $\pi_{n}$ converges to the targeted posterior that imposes the exact sign restrictions at ${\bf{c}}_{N} = {\bf{0}}$. That is, we have $\pi_{N}(\operatorname{\mathbf{Z}})=p_{\operatorname{\mathbf{Z}}}(\operatorname{\mathbf{Z}} \mid ({\mathbf{y}}_t)_{t=1}^{T}, \mathbf{S}_R(T(\operatorname{\mathbf{Z}})) > {\bf0})$. Rather than specifying the sequence of threshold vectors in advance, we devise an adaptive rule that $\operatorname{\mathbf{c}}_{n}$ is determined based on particles from $\pi_{n-1}$, which we will discuss after we introduce the algorithm.
Our algorithm for the initialization starts from drawing particles from the unrestricted posterior. Then, we move through the ladder of intermediate posteriors using an importance-resampling-mutation scheme. We terminate as soon as we find a particle $\operatorname{\mathbf{Z}}_{n}^{i^{*}}$ that satisfies the original sign restrictions. Conceptually, this scheme nests the traditional accept-reject method as a special case with $N=2$ and $\operatorname{\mathbf{c}}_{2} = {\bf{0}}$: if a feasible draw appears at the first stage, the outcome is identical to accept-reject; otherwise, our procedure proceeds through intermediate posteriors with progressively tighter thresholds until feasibility is reached. Several remarks on Algorithm (ref) are described in Online Appendix (ref).
\setcounter{section}{0}