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.
86,556 characters · 12 sections · 90 citation commands
Inference in Tightly Identified and Large-Scale Sign-Restricted SVARs
Imposing sign restrictions is a central tool for identifying structural vector autoregressions (SVARs). These restrictions are typically implemented by first estimating the posterior distribution of the reduced-form parameters and then mapping posterior draws into the structural parameter space. The most commonly used algorithm faust1998, canova2002,uhlig2005, rubioramirez2010 is an accept-reject sampling procedure that draws independent posterior samples of the reduced-form parameters from a conjugate normal-inverse-Wishart distribution, combines them with independent draws from the uniform distribution over orthogonal matrices, and retains only those that satisfy the imposed sign restrictions. As shown by arias2018, this algorithm generates independent draws of the structural parameters from a well-defined posterior distribution induced by the reduced-form prior under the imposed sign restrictions. It has two appealing properties: in addition to producing independent and identically distributed posterior draws, it can readily accommodate various additional identifying restrictions, such as ranking and elasticity bounds restrictions.
However, the performance of the accept-reject algorithm crucially depends on the size of the admissible set of orthogonal matrices satisfying the restrictions. As the restrictions become tighter, the admissible set shrinks and the probability of accepting a draw declines sharply, rendering the algorithm increasingly inefficient and, in extreme cases, infeasible. This concern is particularly relevant because in recent contributions, strong identifying restrictions that substantially reduce the size of the admissible set are increasingly being imposed inoue2026. In the same vein, accept–reject algorithms face challenges in large-scale SVAR models, where the parameter space is high-dimensional and identification often requires imposing a large number of sign restrictions, frequently combined with additional identifying constraints. For example, chan2025 impose more than one hundred sign restrictions to identify eight shocks in a 35-variable SVAR. Such designs naturally entail a substantial computational burden and highlight the need for inference methods that extend beyond standard accept-reject procedures.
Instead of relying on accept-reject algorithms, we propose imposing inequality restrictions through a reparameterization of the structural model by continuously differentiable mappings. The restrictions are thus imposed by construction, which eliminates the need to discard posterior draws and substantially improves computational efficiency, especially when the admissible set is small. Our approach accommodates a broad range of inequality restrictions commonly used in the SVAR literature. These restrictions, applicable in both small and large models, include shape and ranking restrictions on impulse responses at multiple horizons, as well as bounds on the magnitude of economically relevant elasticities. The framework also allows imposing zero restrictions without requiring importance weighting or other corrections typically needed under accept-reject sampling.
The flexibility of the proposed reparameterization comes at a cost, as the induced posterior distribution of the parameters does not belong to a standard family, which precludes direct sampling. Therefore, posterior inference must rely on Markov chain Monte Carlo (MCMC) methods, and we recommend the No-U-Turn Sampler (NUTS) of hoffman2014, a modern variant of Hamiltonian Monte Carlo (HMC) that has only rarely been employed in the SVAR literature. It exploits gradients of the log-posterior and the smooth posterior geometry induced by the reparameterization to explore the posterior effectively. However, gradient-based sampling is computationally demanding, as each iteration requires repeated evaluations of the log-posterior and its gradient. Ensuring fast computation of the posterior density is therefore crucial, and we show how it can be evaluated efficiently under the proposed reparameterization to facilitate inference even in high-dimensional parameter spaces. For numerical work, we use the software package Stan stanreference, containing an implementation of the NUTS with several desirable features, including adaptive tuning, to facilitate estimation.
Recently, kitagawa2025 have argued that standard MCMC algorithms do not perform well in Bayesian inference on set-identified models because set identification induces flat posterior regions along observationally equivalent parameter manifolds, which can severely impair the mixing of these algorithms. However, while they show that random-walk Metropolis performs poorly in this setting, with mixing deteriorating sharply as the size of the identified set increases, those results do not necessarily carry over to modern, gradient-based samplers with adaptive tuning. Indeed, using the simulation design considered in kitagawa2025, we show that the No-U-Turn Sampler as implemented in Stan is able to effectively explore high-dimensional, set-identified posteriors with flat directions, producing well-mixed chains even when the identified set is large.\footnote{The simulation exercise does not actually involve an SVAR model, but in tightly identified SVARs the identification problem is less severe, which suggests that Hamiltonian Monte Carlo is particularly well suited for posterior exploration in set-identified SVAR models, and this conclusion is also supported by our empirical results.}
We are not the first to depart from standard accept–reject algorithms in sign-identified SVAR inference, as already baumeister2015, baumeister2019 propose to use an MCMC sampler to draw the structural parameters directly. Their approach relies on conjugate prior structures, which produce closed-form conditional posteriors for all parameters except those of the contemporaneous structural relations. The latter are first drawn from their marginal posterior, after which inference for the remaining parameters proceeds using their conditional posteriors. This simplifies computation, but comes at the cost of reduced prior flexibility. Moreover, without substantial computational refinements, their algorithm can become computationally demanding in large VAR systems. More recently hou2024 derives the conditional posteriors for the structural parameters under linear equality and inequality restrictions that can be used to build a Gibbs sampler for an SVAR model. His approach has the limitation that for longer horizons of the IRFs, an accept-reject step is needed, which is not the case with our approach.
bruns2023 develop an importance sampler for sign-restricted Bayesian SVARs that uses proposal draws based on a conjugate reduced-form posterior. Their approach improves sampling efficiency and enables estimation beyond the standard conjugate prior setup. However, importance sampling methods can perform poorly in high-dimensional settings, and, as the authors note, the approach is not intended for large VAR systems. read2025 employ slice sampling for posterior simulation and replace binary sign restrictions with a continuous relaxation. Also their algorithm involves an importance sampling step, which is used to obtain draws from the correct target distribution. Furthermore, although their method yields independent posterior draws, it is formulated under a conditionally uniform prior in the way described in uhlig2017, which is not invariant to the set of imposed restrictions, as shown by arias2025elliptical. Our approach instead uses a uniform prior over orthogonal matrices and therefore guarantees invariance to the restrictions (see arias2025uniform).
The inference methods recently put forward by chan2025 and arias2025elliptical are particularly useful in tightly identified and large-scale SVARs. chan2025 exploits the fact that a uniform distribution over the orthogonal matrices is preserved under reordering the columns or flipping their signs, and thus yields multiple parameter draws at a lower computational cost in an accept-reject algorithm. arias2025elliptical propose embedding elliptical slice sampling within a Gibbs sampler, which allows draws of the orthogonal matrix to remain within the identified set by construction. Both approaches deliver substantial computational gains in tightly identified and large-scale SVARs, while maintaining the conventional reduced-form prior structure. However, because of reliance on an accept–reject framework, the approach of chan2025 is not able to achieve the same level of performance as that of arias2025elliptical.
As the method of arias2025elliptical represents the current state of the art in inference in tightly identified and large-scale SVAR models, we focus on benchmarking our approach against theirs in Section (ref). In particular, we revisit their empirical applications to an oil market SVAR identified using sign and elasticity restrictions, as well as a 35-variable SVAR model in which ten shocks are identified by sign and ranking restrictions. In both settings, our HMC approach exhibits faster convergence, as measured by the $\hat{R}$ statistic of vehtari2021. This is likely to follow from the fact that their Gibbs sampler generates Markov chains with extremely persistent serial correlation, which remains sizeable even at lags in the hundreds, whereas autocorrelations in our sampler decay comparatively quickly.
While our sampler performs better in the two applications, the approach of arias2025elliptical offers complementary advantages. In particular, in contrast to their sampler, our approach does not allow for arbitrary combinations of inequality restrictions; for example, in general, long-run and dynamic sign restrictions cannot both be imposed simultaneously. On the other hand, as shown in Section (ref), we are able to impose exclusion restrictions (even at no additional computational cost) which their approach does not afford. Our method also supports highly flexible prior specifications, which need not yield closed form conditional posteriors, although this may come at a computational cost.
The remainder of the paper is organized as follows. In Section (ref), we present the SVAR model. Section (ref) introduces the identification problem and the restrictions that can be implemented in our framework. In Section (ref), we describe the prior and posterior, our method of imposing identifying restrictions, and an efficient method of computing the posterior density, which is crucial for the implementation of our approach in tightly identified and large-scale SVAR models. Section (ref) illustrates the performance of our approach in empirical applications. Finally, Section (ref) concludes.
We consider the SVAR model
where $y_t$ is an $N$-vector of observables, $c$ is an $N$-vector of constants, $A_i$, $i=1,\dots p$, are $N\times N$ autoregressive coefficient matrices, and the $N\times N$ matrix $A_0$ summarizes the contemporaneous structural relations among the elements of $y_t$. The $N$-vector of structural shocks $\epsilon_{t}$ has mean zero and identity covariance matrix $I_N$, and its elements are independent across time. We assume that $\epsilon_{t}$ is jointly Gaussian, but this assumption can be relaxed, for example, if we wish to make use of non-Gaussianity in identification or to place restrictions on higher moments.
The SVAR model (ref) is in the so-called A-form (see Lutkepohl2005, Chapter 9), but it can equivalently be written in the B-form
where $\tilde{c} = A_0^{-1}c$, $\tilde{A}_i = A_0^{-1}A_i$, and $B = A_0^{-1}$. This representation is useful for estimation because the matrix $B$ gives us the impact responses directly, whereas $A_0$ needs to be inverted to obtain them.
The usual objects of interest in SVAR analysis are the structural impulse response functions (IRF), which trace out the responses of observables at different time horizons $t+i$, $i=0,1,2,\dots$, for a shock occurring at time $t$. The structural IRFs of the SVAR model are obtained recursively as $\Psi_0 = A_0^{-1} = B$, and
The $(l,k)$ element of $\Psi_i$, $\Psi_{i,lk}$, is the response of variable $l$ to shock $k$ at horizon $i$.
This section describes the types of identifying restrictions that can be accommodated within our framework. In Section (ref), we illustrate how these restrictions can be implemented in practice.
In the reduced-form vector autoregression (VAR) corresponding to model (ref),
$u_t = A_0^{-1} \epsilon_{t}$ are the reduced-form errors with mean zero and covariance matrix $\Sigma = A_0^{-1}(A_0^{-1})'$. The structural model suffers from the well-known identification problem that the covariance matrix of the reduced-form residuals is identified, but any $\bar{A}_0$ satisfying $\bar{A}_0^{-1} = A_0^{-1}Q$, for an orthogonal matrix $Q$, yields the same $\Sigma$, that is, the same covariance structure of $u_t$. Therefore, the corresponding models are observationally equivalent.
To identify the structural parameters of an SVAR model, the researcher must use prior information, typically in the form of identifying restrictions. Currently, imposing sign restrictions is probably the most widely used method of identification. In the implementation of uhlig2005, rubioramirez2010, candidate draws of the impact matrix are generated from the product $PQ$, where $P$ is the Cholesky factor of the posterior covariance matrix $\Sigma$ of $u_t$ and $Q$ is an orthogonal rotation matrix drawn from its prior distribution, and only draws satisfying the sign restrictions are kept. However, it is well known that sign restrictions do not point-identify the model, and hence, the prior is not updated within the identified set in Bayesian estimation baumeister2015. This problem can be avoided by sufficiently tight identification of the orthogonal rotation matrix $Q$, as in that case most of the posterior uncertainty comes from uncertainty about the reduced-form parameters inoue2020.
However, tight identification of $Q$ poses computational challenges, as the probability of randomly drawing a candidate $Q$ satisfying all restrictions imposed on the model is inversely proportional to the size of the identified set. Hence, estimation time can be prohibitively long when the restrictions are tight and/or there are a lot of them. To that end, several solutions have recently been proposed to increase the efficiency of sampling the orthogonal matrix $Q$. First, chan2025 exploit the fact that permutations of columns and flips of their signs preserve a uniform distribution over $Q$, which facilitates obtaining multiple candidate draws from a single $Q$. Second, arias2025elliptical use elliptical slice sampling, which is a particular Markov Chain Monte Carlo (MCMC) algorithm, to sample $Q$ such that the draws remain within the identified set. Finally, read2025 approximate sign restrictions by a continuous function, which allows more efficient sampling using the slice sampler. However, their procedure requires an additional importance sampling step to ensure that the draws follow the correct distribution, and it is formulated under a conditionally uniform prior, which is not invariant to the set of imposed restrictions, as discussed in the Introduction.
Instead of sampling $Q$, the method proposed in this paper works with the impact matrix $B$ directly, which facilitates imposing sign restrictions on the parameters of interest using continuously differentiable mappings to transform them into an unconstrained parameter space. Thus, there is no need to perform any rejection sampling. This approach works well even in very tightly identified models, as regardless of the size of the identified set, the sampler always runs on an unrestricted parameter space. Our method allows for sign restrictions on the impulse responses also beyond the impact effect, as long as the horizon to be restricted does not exceed the VAR lag length $p$. We believe that this limitation on the horizon is unlikely to be important in most practical applications, especially in light of recent work by olea2025, which recommends using relatively long lag lengths.
Besides sign restrictions, our approach facilitates imposing other kinds of identifying restrictions. In contrast to arias2018, exclusion restrictions can be imposed directly during sampling without any need for an importance sampling step. Also, the effect of a given shock can be constrained to increase over a specified horizon, or a shock can be required to have the largest impact effect on a given variable. The exact method of placing the latter kinds of restrictions will be presented in Section (ref). Finally, sign restrictions can be combined with other kinds of restrictions, including restrictions on elasticities, as in kilian2012, who impose them to avoid implausibly large estimates of the short-run elasticity of oil supply under sign identification alone. More generally, we can impose ranking restrictions requiring a shock to have a greater or smaller effect on a given variable than on another variable (see amirahmadi2021). The use of these kinds of restrictions is illustrated in Sections (ref) and (ref), respectively.
Our approach is flexible, but there are certain kinds of restrictions and combinations of restrictions that it does not permit, and in these cases, the samplers of arias2025elliptical or chan2025 are a viable choice. In particular, because the sampler operates on IRF parameters, restrictions can only be stated in terms of the IRFs; constraints that jointly involve IRF values and VAR matrices are therefore not feasible. For example, long-run restrictions generally cannot be combined with dynamic sign restrictions on the IRFs. More broadly, only restrictions that admit a diffeomorphism\footnote{A diffeomorphism is a continuously differentiable bijection with a continuously differentiable inverse. This permits a wide range of restrictions, but, for example, restrictions on forecast error variance decompositions (FEVD) are not compatible because the squaring operations required for their calculation are not bijective.} to an unconstrained parameter vector are allowed. It is also not possible to impose restrictions on the likelihood, such as the narrative restrictions of antolindiaz2018.
Following standard practice in the Bayesian VAR literature, we impose a conjugate normal–inverse Wishart (N-IW) prior on the reduced-form parameters $\tilde{A}=\left[\tilde{c}, \tilde{A}_1,\dots,\tilde{A}_p\right]'$, and the covariance matrix $\Sigma$. This choice yields a posterior distribution of the same form and improves computational efficiency. The proposed method, however, does not rely on conjugacy and remains valid under more general prior specifications.
Specifically, we assume $\Sigma \sim IW(\nu_0, S_0)$ and $\operatorname{vec}(\tilde{A}) \mid \Sigma \sim N(\operatorname{vec}(\Phi_0), \Sigma \otimes \Omega_0)$, where $IW(\cdot,\cdot)$ and $N(\cdot,\cdot)$ denote the inverse Wishart and multivariate normal distributions, respectively, and $\operatorname{vec}$ is matrix vectorization. The quantities $\nu_0, S_0, \Phi_0$ and $\Omega_0$ are prior parameters. Under these assumptions, the posterior distribution is again normal–inverse Wishart, with updated parameters given by
The rotation matrix $Q$ is assigned a uniform prior over the space of orthogonal matrices, independent of the reduced-form parameters. Because the likelihood is invariant to orthogonal rotations, the posterior distribution of $Q$ coincides with the prior.
Our aim is to sample from the structural parameters of the model, which requires a transformation of the posterior density of the orthogonal reduced-form parameters into the implied density over the structural parameters of interest. Related transformations from the orthogonal reduced-form parameters to the structural objects have been derived in earlier work. For instance, arias2018 characterize the induced distribution of the structural parameters in the A-model, referring to it as the normal–generalized normal distribution, while inoue2013 derive the distribution of the first $p+1$ structural impulse response matrices. inoue2022 provide a further generalization by allowing the number of IRF matrices to differ from $p+1$.
While we build on inoue2013, inoue2022, we adapt their results to our setting in which the posterior is parameterized in terms of the first $k$ impulse response matrices $\Psi_0,\Psi_1,\dots,\Psi_k$ (on which the identifying restrictions are imposed) and the remaining VAR coefficient matrices $(\tilde{A}_{k+1},\dots,\tilde{A}_p)$. As our interest is in estimating a full SVAR model, we keep the VAR matrices in the density, as opposed to inoue2022 who only consider the IRF parameters. This mixed parameterization also improves the efficiency of posterior simulation because the VAR coefficients tend to exhibit weaker posterior dependence than impulse response parameters. In addition, by exploiting the properties of the $LQ$ decomposition, we obtain a simpler expression for the associated Jacobian determinant than inoue2013. Our expression does not include the orthogonal matrix $Q$, which also saves on calculations as we do not need to evaluate a determinant of another large matrix.
The following proposition characterizes the posterior of the structural parameters collected in $\Pi_k=(\tilde{c},B,\Psi_1,\dots,\Psi_k,\tilde{A}_{k+1},\dots,\tilde{A}_p)$ in terms of the N-IW posterior of the orthogonal reduced-form parameters.
The proof of Proposition (ref) does not rely on the distribution over the reduced-form parameters being N-IW, and thus any other prior or error distribution can be used as well. However, we state this proposition using the N-IW density for compatibility with much of the existing literature, and because it allows fast estimation of large SVARs, as will be demonstrated in Section (ref)
As our goal is to incorporate both inequality and zero restrictions, we also need a result that gives us the density function under these restrictions. The following proposition is a direct result of standard facts, but we state it here for clarity.
In the context of our SVAR model, the density $f$ is the density of the structural parameters in (ref). The set-identifying restrictions (like sign and ranking restrictions) and exclusions restrictions are imposed on $\theta_1$ and $\theta_2$, respectively, while $\theta_3$ collects all free parameters of the model. The key point of Proposition (ref) is that evaluation of the conditional density under the exclusion restrictions requires no normalizing constants. Therefore, in the case of a Metropolis-style algorithm like HMC, it is sufficient to evaluate the unconstrained density on the subspace satisfying the restrictions. It is important to point out that Proposition (ref) is more general than our method affords. In particular, we cannot restrict $\theta_1$ in a totally arbitrary way, as discussed at the end of the previous section. We also confine equality restrictions to zero restrictions to preculde the possibility of disjoint posteriors kitagawa2025.
In this section, we show how to impose identifying restrictions on the parameters using continuously differentiable transformations, which result in an unrestricted parameter vector. Our sampler runs on this vector, which means that rejection sampling is not needed to impose the identifying restrictions. Implementation of exclusion restrictions is straightforward in that no Jacobian adjustment is needed, whereas restrictions based on inequalities do require it as they are based on a change of variables. As a different set of identifying restrictions is considered in each application, there are no general transformation and Jacobian that would cover all possible cases, but we discuss a few building blocks which can be used to construct frequently encountered restrictions.
Exclusion restrictions can be imposed by setting the values of the parameters of interest to zero during the evaluation of the density function, as shown in Proposition (ref). For example, if prior economic knowledge dictates that the response of the first variable to the second shock on impact is zero, we set $B_{1,2} = 0$, where $B_{i,j}$ denotes the $(i,j)$ element of the matrix $B$, and drop the corresponding parameter from the parameter vector. Hence, the posterior density is evaluated only on the subspace where the restriction holds.
Imposing inequality restrictions calls for suitable transformations of the parameters and the associated Jacobian determinant denoted by $|J|$. Combining the Jacobian with the density of the structural parameters in (ref), yields the density over the unrestricted auxiliary vector $\theta$:
This is the density that our algorithm samples from. The exact cost of computing $|J|$ will vary depending on the number and type of restrictions imposed on the model, but in applications involving sign restrictions or elasticity bounds it is negligible in comparison to the cost of evaluating the normal-inverse Wishart density. For example, if we restrict the sign of every parameter in $B$, evaluation of $\log|J|$ collapses to computing the sum over the $N^2$ elements of $\theta$ corresponding to $B$.
Let us consider the implementation of frequently encountered types of inequality restrictions. Although the presentation is not exhaustive, it should give an idea of how to proceed in an empirical application. Suppose first that we have a density $p_\alpha(\alpha)$ over a parameter $\alpha\in\mathbb{R}_+$, which could, for example, be an impact response parameter from the matrix $B$. A Metropolis-type sampler running on $\alpha$ produces proposals $\alpha < 0$, which are always rejected, and this decreases efficiency. Instead, we can reparameterize the sampler in terms of an unrestricted parameter, so that we do not encounter a priori known regions of zero probability. For the sign restricted case, we introduce a parameter $\tilde{\alpha}\in\mathbb{R}$, and set $\alpha = e^{\tilde{\alpha}}$. The Jacobian of this transformation is $e^{\tilde{\alpha}}$, and the corresponding density function is then $p(\tilde{\alpha}) = p_\alpha(e^{\tilde{\alpha}})e^{\tilde{\alpha}}$ by a change of variables. If we want to restrict $\alpha < 0$, we set $\alpha = -e^{\tilde{\alpha}}$, and the corresponding density is $p(\tilde{\alpha}) = p_\alpha(-e^{\tilde{\alpha}})e^{\tilde{\alpha}}$.
Another useful transformation is the logit transformation for a lower and upper bounded scalar. Suppose we want to sample a parameter $\alpha\in(a,b)$, say, an elasticity parameter which needs to be positive and also constrained from above. To that end, we define an $\tilde{\alpha}\in\mathbb{R}$, and obtain $\alpha$ as
where
In this case, the absolute value of the Jacobian is
and the corresponding density function is obtained by a change of variables, as in the sign-restricted case above. Notice that the upper and lower bounds can be parameters themselves, in which case the dependencies need to be accounted for when calculating the Jacobian.
Besides sign and bound restrictions, ordering restrictions related to the relative sizes of the effects of the shocks are often imposed. For instance, we might want to declare a priori that a certain shock has a greater impact effect on a given variable than any other shock. As an example in the context of monetary policy, we could require that the impact effect of a contractionary monetary policy shock on the interest rate be positive and greater in absolute value than that of any other shock.\footnote{The impact effect must be greater in absolute value, not just greater than the impact effect of any other shock because otherwise the estimated effect of the monetary policy shock may change sign and switch columns. In such a case, the column assigned to the monetary policy shock could correspond to a shock with a small estimated effect in magnitude, yet greater than that of the estimated monetary policy shock.}
These ordering restrictions can be implemented by imposing restrictions on a row of the impact matrix $B$. In the case of the monetary policy example, let $(\alpha_1,\dots, \alpha_N)$ denote the row of $B$ corresponding to the interest rate, where $\alpha_N$ denotes the impact effect of the monetary policy shock. Formally, we sample a vector $(\alpha_1,\dots, \alpha_N)\in\mathbb{R}^N$ subject to the restrictions
By defining $\tilde{\alpha}\in\mathbb{R}^N$, the transformation can be written componentwise as
The partial derivatives of $\alpha_N$ needed for the Jacobian are
For the remaining elements $i\in\{1,\dots, N-1\}$, we have
Thus, the Jacobian matrix has non-zero entries only on the diagonal and in the last column, implying that the determinant is just the product of the diagonal elements,
If the positive element is assigned to a position other than $N$, the absolute value of the Jacobian determinant remains unchanged, as a permutation only multiplies the Jacobian by $\pm 1$.
In an application one would decide on the set of the identifying restrictions imposed on the model, and then implement them using a combination of the aforementioned transformations or other transformations that are feasible in that particular setting. The Jacobian of these transformations is then incorporated into the density function in (ref).
To illustrate these methods in practice, let us consider a two-variable model with a single lag and no constant. Suppose that the variables are quantity and price, and that the supply and demand shocks are identified using sign restrictions on the impact responses of the variables. The responses have opposite signs for the supply shock and the same sign for the demand shock. If the first shock is the supply shock, the first column of the impact matrix $B$ contains one positive and one negative element, whereas the second column contains two positive elements.
To sample from the model, we define an unrestricted auxiliary vector $\theta\in\mathbb{R}^8$. The first four entries of this vector correspond to the elements of the impact matrix $B$, and the last four entries correspond to the elements of the VAR matrix $\tilde{A}_1$. As no restrictions are imposed on the elements of $\tilde{A}_1$, they enter $\theta$ as such. Specifically,
To impose the sign restrictions, we use the exponential transformation. Thus, the structural parameter matrices become
and the Jacobian matrix $J$ is diagonal and given by
where $I_4$ denotes the $4\times 4$ identity matrix. Therefore, the absolute Jacobian determinant is
Given the Jacobian and the matrices $B,\tilde{A}_1$, we can now proceed to evaluate the density using formula.(ref).
As the density function (ref) does not facilitate direct sampling, we employ an MCMC method. Specifically, we use the No-U-Turn Sampler (NUTS) of hoffman2014, which is a variant of the Hamiltonian Monte Carlo (HMC) algorithm. An HMC algorithm explores the target distribution by simulating Hamiltonian dynamics, which allows the sampler to construct long trajectories while yielding a high acceptance probability by utilizing the gradients of the target log-density. At a high level, the HMC sampler can be thought of as a Metropolis-Hastings algorithm with a complicated proposal distribution. Its performance depends on tunable parameters $\varepsilon, L$ and $M$. Starting from the current draw, the HMC sampler simulates a trajectory by integrating Hamiltonian dynamics using $L$ steps of size $\varepsilon$ with a mass matrix $M$, whose optimal choice is the inverse of the covariance matrix of the target posterior distribution betancourt2018. The resulting endpoint is then used as a proposal.
The NUTS solves the issue of tuning $L$ by dynamically expanding the path taken within each iteration until it starts to make a U-turn towards the initial state. The implementation of NUTS in Stan further provides routines for adapting the step size and metric matrix during warmup, and calculates the gradient of the target log-density by automatic differentiation. These properties reduce the scope of implementation error and ensure that sampling can be done efficiently in a wide variety of situations. Further details on the HMC sampler and NUTS can be found in Appendix (ref).
Beyond the built-in efficiency gains provided by Stan, sampling performance can be further improved through appropriate parameterization and initialization. As discussed in detail in Appendix (ref), if only impact restrictions are considered, sampling efficiency can be enhanced by adopting a non-centered parameterization, where we sample a standard normal (instead of a general multivariate normal) vector betancourt2013. Moreover, finding reasonable initial values is of utmost importance, and in the applications of Section (ref) we initialize the sampler by drawing from the normal-inverse Wishart posterior and constructing an orthogonal matrix compatible with the imposed restrictions one column at a time, following the method of arias2025elliptical.
It is important to note that set identification poses challenges that persist even under efficient MCMC implementations. In particular, it induces flat posterior regions along observationally equivalent parameter manifolds. Because the likelihood function provides no information in these directions, local MCMC updates receive little guidance on how to move within the identified set, which can hinder effective posterior exploration. kitagawa2025 formally show that this issue can be severe for random-walk Metropolis algorithms, with mixing deteriorating as the size of the identified set increases. However, the mechanism underlying this result does not necessarily extend to Hamiltonian Monte Carlo. This is the case because, in contrast to random-walk Metropolis, HMC algorithms rely on Hamiltonian trajectories rather than diffusive local proposals, and therefore they explore flat posterior regions differently. Moreover, modern HMC implementations (including the one used in this paper) employ automatic step-size and mass-matrix adaptation during warmup, which further helps mitigate these exploration difficulties.
To demonstrate that an HMC sampler can effectively explore the posterior under set identification, we revisit a simulation experiment that kitagawa2025 considered to illustrate the opposite. Although our primary interest lies in VAR models, we adopt this example because it features a known and simple posterior geometry, which allows for a transparent diagnosis of convergence issues. Specifically, we generate data from a normal distribution with mean $\mu=\sum_{i=1}^{k} \mu_i$ and variance unity, where the parameters $\mu_i$, $i=1,\dots,k$, are not identified for $k \geq 2$. We study the performance of Hamiltonian Monte Carlo in this setting for $k = 2$ and $k = 1{,}000$, with the latter case reflecting the high-dimensional parameter spaces encountered in VAR applications. For $k=2$ and $k=1{,}000$, we generate 1{,}000 observations from $N(2,1)$ and $N(10,1)$ distributions, respectively. The larger mean in the latter case reflects the greater number of additive components.
The left panel of Figure (ref) presents a trace plot based on 12{,}000 draws for the case $k=2$, of which 2{,}000 are warmup draws. Following kitagawa2025, the warmup draws are included. The right panel depicts the histograms of the obtained sample (without the warmup iterations) together with the analytical marginal densities, when the support of each $\mu_i$ is bounded between $-10$ and 10.\footnote{In the case $k=2$, the marginal density is proportional to $\Phi(\sqrt{n}(\mu_i+10-\bar{x})) - \Phi(\sqrt{n}(\mu_i-10-\bar{x}))$, where $\Phi$ is the standard normal CDF, $n$ is the sample size, and $\bar{x}$ is the sample mean. For $n=1{,}000$ this density is approximately uniform.} The plots indicate that the NUTS (as implemented in Stan) is able to traverse the entire posterior distribution of $\mu$ without apparent difficulties and produces draws that closely match the analytical marginal densities. This contrasts with the behavior documented by kitagawa2025 for a standard HMC algorithm in the same set-identified setting, where the sampler failed to adequately explore the parameter space associated with $(\mu_1,\mu_2)$.
The NUTS performs well even in a high-dimensional setting, as illustrated in Figure (ref) that presents a scatter plot of the first two components of $\mu$ and the estimated autocorrelation functions for these two parameters when $k = 1{,}000$. The draws are dispersed without visible clustering, and the autocorrelations decay rapidly toward zero. Corresponding plots for other pairs of parameters look essentially identical. These results indicate that, even in this high-dimensional setting, the NUTS produces high-quality samples and efficiently explores the full target distribution.
As bacchiocchi2025 point out, besides flat posterior regions, a related but distinct challenge can arise when the identified set is multimodal because separated regions of posterior mass are difficult for MCMC algorithms (including HMC) to explore effectively. However as noted by kitagawa2025, such local identification arises only in the case of non-homogenous restrictions, that is, when a parameter is set to a non-zero value. As we preclude such restrictions, multimodality does not arise, and hence the sampler does not need to traverse between distinct modes. This absence of multimodality follows from the fact that the likelihood depends on $B$ only through $BB'$ and is therefore invariant under orthogonal rotations, implying a continuum of observationally equivalent factorizations rather than isolated peaks.\footnote{To eliminate any remaining reflection symmetry, we impose the normalization $\det(B)>0$ and restrict attention to a single normalized representation of the structural parameters. This normalization does not restrict the identified set beyond our existing identifying assumptions. To see this, consider the case where at most $N-1$ columns of $B$ are identified, so the last column is left unrestricted. Let $B_1$ satisfy the identifying restrictions and have $\det(B_1)<0$. Define $R=\mathrm{diag}(1,\ldots,1,-1)$ and $B_2 = B_1R$, i.e., $B_2$ equals $B_1$ with the sign of its last column flipped. Then $\det(B_2)=\det(B_1)\det(R)=-\det(B_1)>0$, while the first $N-1$ columns of $B_2$ and $B_1$ are equal. Hence, $B_2$ satisfies the same identifying restrictions and has the same economic interpretation. Moreover, the map $B\mapsto BR$ is bijective with inverse $B\mapsto BR^{-1}=BR$, which guarantees that there always exists exactly one such $B_2$.} As the preceding example indicates, the HMC sampler is well suited for exploring such distributions.
Even though the NUTS has several advantages compared to commonly used MCMC algorithms, it is important to monitor the convergence of the MCMC chain using diagnostic checks. To that end, we recommend the tools of vehtari2021, the split-$\hat{R}$, and the bulk and tail effective sample size (ESS) statistics. We monitor these univariate statistics for all parameters because, as noted by brooks1998, spurious results may be obtained by monitoring only a subset of the parameters.
The $\hat{R}$ measures chain mixing by comparing within- and between-chain variances. For a well-mixed chain, it approaches unity as the number of draws increases. When running a single chain (as is typical in the SVAR literature), $\hat{R}$ is calculated by splitting the chain into two halves. The idea of the bulk and tail ESS is to measure sampling efficiency by adjusting the number of iterations for the autocorrelation in the sample. The bulk ESS estimates the number of effectively independent samples in the MCMC output for estimating the mean or the median of the posterior, while the tail ESS measures sampling efficiency in the tails of the distribution. vehtari2021 recommend accepting an MCMC sample only if $\hat{R} < 1.01$ and the ESS is at least 100 for a single chain. These thresholds ensure sufficient mixing and a sufficiently low Monte Carlo standard error for estimation of the quantities of interest.
The univariate ESS has a multivariate counterpart vats2019 that also accounts for correlations between the parameters when evaluating sample quality. Following arias2025elliptical, we report this statistic to assess the performance of the sampler for the IRF parameters of interest in Section (ref). However, we rely on $\hat{R}$ as the primary convergence diagnostic because it directly assesses chain mixing and can detect non-convergence even when the ESS is large. In contrast, the multivariate ESS is primarily a measure of sampling efficiency and may fail to indicate lack of convergence when the chain has not fully explored the target distribution.
When we use a Hamiltonian Monte Carlo (HMC) algorithm to estimate the parameters, we need to repeatedly evaluate the posterior density of the structural parameters. As is evident from equation (ref), and given the prior we have assumed, this entails repeated evaluation of the normal-inverse Wishart density
We contribute to the literature by presenting a fast evaluation method for the N-IW density. Although the N-IW SVAR model is ubiquitous in the literature, computation of its density function is not routinely performed, as there are methods to generate N-IW random variables using direct sampling. Our method makes extensive use of solvers for triangular systems of equations, and the Cholesky factors of the elements of the covariance matrix of the N-IW distribution. Thus, it has similarities with the work of carriero2016, chan2021, who use the said Cholesky factors to efficiently generate normally distributed random variables with a Kronecker covariance structure.
Evaluation of $|\Sigma|$ is standard and uses the diagonal elements of the Cholesky factor of $\Sigma$. We obtain the Cholesky factor $L_\Sigma$ of $\Sigma$ by a Cholesky decomposition of $\Sigma = BB'$; another option would be to apply the $LQ$ decomposition\footnote{Software packages do not always provide a direct implementation of the $LQ$ decomposition. In that case, we can apply the QR decomposition to $B' = \tilde{Q}R$ and set $L = R'$ and $Q = \tilde{Q}'$. Another aspect to note is that the QR-decomposition for a square invertible matrix is unique only up to the signs of the diagonal of $R$. We assume the convention that these are positive, which yields a unique decomposition.} to the matrix $B$ and take $L$ to be the Cholesky factor of $\Sigma$. In our experience, the former is faster, but the latter provides the matrix $Q$ that might be of interest. Furthermore, as $|\det(B)| = \sqrt{\det(\Sigma)}$, there is no need to evaluate the determinant of $B$ separately to implement the density (ref).
The trace term in the first exponential on the right-hand side of (ref) can be evaluated using the Cholesky decompositions $\tilde{S} = L_{\tilde{S}} L_{\tilde{S}}'$ and $\Sigma = L_\Sigma L_\Sigma'$. The value of $\tilde{S}$ does not change between the iterations of the sampler, so it is pre-computed and cached. Using the cyclic property of trace, we can rearrange the matrices as
where $Z = L_\Sigma^{-1} L_{\tilde{S}} $. Matrix $Z$ can be solved efficiently using forward substitution from the equation $L_\Sigma Z = L_{\tilde{S}}$. Finally, by the properties of the trace of a matrix product,
where $Z_{ij}$ denotes the $(i,j)$ element of $Z$.
The computationally most intensive part is the evaluation of the latter exponential term in (ref). As in the case of $\tilde{S}$, the Cholesky decomposition $\tilde{\Omega} = L_{\tilde{\Omega}} L_{\tilde{\Omega}}'$ is also pre-computed and cached. Using the properties of the Kronecker product, we get
Let us denote $W = L_{\tilde{\Omega}}^{-1} (\tilde{A}-\tilde{\Phi}) L_\Sigma'^{-1}$. We can now solve for $W$ in two steps using forward substitution. First, solve $C$ from $CL_\Sigma' = (\tilde{A}-\tilde{\Phi})$. Next, solve $W$ from $L_{\tilde{\Omega}} W = C$. Finally, we get
where $W_{ij}$ denotes the $(i,j)$ element of $W$. This calculation is substantially faster than explicitly computing the $N(Np+1) \times N(Np+1)$ matrix $\Sigma\otimes\tilde{\Omega}$ and working with it.
This section revisits two empirical applications to showcase the benefits of our approach: In Section (ref), we identify three shocks in the small-scale global oil market model of kilian2014, and in Section (ref), we consider the large-scale U.S. economy model recently analyzed by chan2025. In each case, the empirical findings coincide with prior results based on the same identifying restrictions, but we highlight the advantages of our method over existing alternatives. Also arias2025elliptical consider these applications to illustrate the performance of their Gibbs sampler, which facilitates straightforward comparison. As their approach represents the current state of the art for inference under identifying sign restrictions, it provides a natural benchmark for our analysis. In both applications, our method achieves convergence approximately two to eight times faster than theirs, as measured by the time required for the $\hat{R}$ statistic to fall below the conventional threshold. The obtained samples also exhibit lower autocorrelations, as reflected in greater bulk and tail ESS values in both cases. In the large-scale model, our approach is also significantly faster than the algorithm of chan2025, although they use the asymmetric conjugate prior chan2022 rather than the normal-inverse Wishart prior. In Section (ref), we also show that our approach accommodates exclusion restrictions at no additional computational cost.
To illustrate how our approach handles contemporaneous and dynamic sign restrictions, as well as elasticity restrictions, we estimate the oil market model of kilian2014 on monthly data from 1973:M2 to 2009:M8.\footnote{We obtain the data from the replication package for kilian2014 in the Journal of Applied Econometrics Data Archive.} It contains four variables, namely, the change in global oil production, a measure of global real activity, the real price of oil, and the changes in global oil inventories, as well as an intercept, and monthly dummy variables. Following kilian2014, the lag length is set at 24. The prior is flat and proportional to $|\Sigma|^{-(N+1)/2}$, which results in a conjugate normal-inverse Wishart posterior for the reduced-form parameters. Following the identification strategy of kilian2014, we identify the flow supply, flow demand, and speculative demand shocks. In particular, we impose sign and elasticity restrictions on the impact responses, as well as dynamic sign restrictions, as summarized in Table (ref). These restrictions result in a tightly identified model, and thus random sampling of the orthogonal rotation matrix $ Q $ would result in extremely few accepted posterior draws, as discussed in Section (ref).
Our approach solves the issue of tight identification by transforming the constrained structural parameters into unconstrained auxiliary parameters. Specifically, given that we have dynamic restrictions up to the 12th IRF horizon, we work with the constrained structural parameters
where $ \tilde{D} $ is a $4 \times 11$ matrix of dummy coefficients, with each column corresponding to a dummy variable. These parameters are transformed into an unconstrained vector $\theta$ using the transformations described in section (ref), specifically the exponential and inverse logit transformations on $B$ and $\Psi_{1},\dots, \Psi_{12}$. We run 2{,}000 warmup iterations, followed by 100,000 sampling iterations.
Figure (ref) depicts the posterior IRFs. As expected, the results are similar to those reported by kilian2014, arias2025elliptical. In particular, following a negative flow supply shock, oil production, global economic activity and oil inventories decline, while the real price of oil increases. A positive flow demand shock leads to a persistent increase in global economic activity and the real price of oil. It also induces a positive response in oil production, which peaks after about one year and then returns to its pre-shock level. Finally, following a positive speculative demand shock, the real price of oil and inventories increase, while global economic activity and oil production experience a slight but persistent decline.
We now turn to the properties of our sampler. Figure (ref) presents trace plots of the unnormalized posterior log-density and a parameter ($B_{11}$) from the impact matrix $B$. For clarity, only the first 4{,}000 iterations, including warmup, are shown. As can be seen, both the log-density and the parameter oscillate around a constant value starting from the first iteration onward. This indicates that our method for selecting initial values by direct sampling works well and initializes the chain near the typical set. The traceplots do not exhibit problematic behavior such as nonstationarity or high autocorrelation after warmup. Furthermore, the variance of $B_{11}$ is lower during the early warmup phase than during sampling, which reflects the dynamic adaptation of the step size and inverse metric. This adaptation leads to more efficient sampling after warmup. The remaining structural parameters exhibit similar behavior and their traceplots are therefore not shown.
To compare our approach with the elliptical slice sampler of arias2025elliptical, we implemented their sampler in R and, following their lead, save every tenth draw from the resulting chain. We run this sampler for a total of four million iterations.\footnote{We also implemented the sampler in Matlab, and it turned out faster than the R implementation. All results involving runtimes are based on this faster implementation.} To assess convergence, we use the $\hat{R}$ statistic (see Section (ref)), calculated from a single chain by considering samples of increasing length, with warmup iterations included.\footnote{Because there is no clear way to ensure that the discarded warmup samples are comparable across samplers, we include the warmup iterations in the analysis.} As recommended by vehtari2021, we interpret values of $\hat{R}$ above 1.01 as an indication of potential convergence issues.
Figure (ref) plots the maximum values of the $\hat{R}$ statistic for samples obtained using both samplers. The reported sample size is increased in increments of 1,000 iterations. For the elliptical slice sampler, the maximum $\hat{R}$ is the maximum over the elements of $\Sigma$ and $Q$ , whereas for the NUTS it is the maximum over the elements of $B$.\footnote{In our experience, the unidentified parameters $B$ and $Q$ mix most slowly. To reduce computational cost, we therefore do not compute $\hat{R}$ for all parameters.} As can be seen in the figure, $\hat{R}$ converges in fewer iterations when using the NUTS. For the elliptical slice sampler, $\hat{R}$ first falls below the threshold after approximately 29 thousand iterations (corresponding to roughly 14 minutes), but requires about 100 thousand saved draws (47 minutes) to remain consistently below the threshold of 1.01. Furthermore, it exceeds the threshold again at around 180,000 iterations (84 minutes). In contrast, the NUTS reaches the 1.01 threshold after roughly eight thousand iterations, corresponding to a runtime of 28 minutes.
To investigate the sampling efficiency of the samplers, we consider the bulk and tail ESS measures of vehtari2021. In line with the $\hat{R}$, they are calculated for all elements of $B$ (for the NUTS), and $\Sigma$ and $Q$ (for the slice sampler), again using chains of increasing length and including warmup iterations. Figure (ref) presents the minimum ESS values over these parameters for the NUTS and the elliptical slice sampler. As can be seen, the NUTS is more efficient in terms of ESS per iteration under both measures. This is especially evident in the case of the bulk ESS, indicating that our approach yields more reliable estimates of posterior means for a given number of iterations. Importantly, for both methods, ESS increases approximately linearly with the number of iterations and reaches levels at which the Monte Carlo error is negligible for inference.
The non-monotone behavior of the $\hat{R}$ statistic and the lower ESS values for the elliptical slice sampler are likely due to the high autocorrelation in the sampled chain. Figure (ref) presents the autocorrelation function for $B_{11}$ from the NUTS (top panel), and $Q_{11}$ from the elliptical slice sampler (bottom panel). The latter exhibits substantially higher autocorrelation, which reduces sampling efficiency. Similar patterns are observed for the other elements of the matrices \(B\) and \(Q\) in this application.
We also assess convergence using the minimum multivariate ESS, as suggested by vats2019. To achieve a precision of 0.05 for the 95% confidence region our approach requires 63 minutes of sampling, whereas the elliptical slice sampler requires 29 minutes. Precision of 0.05 implies that the Monte Carlo error is 5% of the variability in target distribution. While this metric suggests that the elliptical slice sampler is more efficient in this setting, it should be interpreted with caution, as the multivariate ESS primarily reflects sampling efficiency and may fail to detect lack of convergence, as discussed in Section (ref). arias2025elliptical report the efficiency in terms of time required to generate 1000 multivariate ESS. The values are 0.100 and 0.080 hours per 1000 ESS for the NUTS (including warmup) and the elliptical slice sampler, respectively. Following arias2025elliptical, the multivariate ESS is calculated jointly for the first three columns of the impulse response parameters, corresponding to the identified shocks in this model.
As an example of a large SVAR model, we consider a 35-variable quarterly model of the U.S. economy with five lags and a constant. This application extends the work of chan2025, which in turn builds on crump2021. The sample period spans from 1975:Q2 to 2020:Q4.\footnote{We thank Christian Matthes for sharing the replication files.} The variables are listed in Table (ref), which also gives the identifying restrictions. Following arias2025elliptical, both sign and ranking restrictions are imposed to identify ten structural shocks: demand, investment, financial, monetary, government spending, technology, labor supply, wage bargaining, oil price, and consumer sentiment shocks.
Following arias2025elliptical, we place a Minnesota prior on the reduced-form parameters.\footnote{We thank Jonas Arias for sharing their prior with us.} As the model contains only contemporaneous restrictions, we estimate it using the non-centered specification, where instead of sampling the VAR coefficient matrices, we sample an auxiliary standard normal vector and transform it into the desired parameters $\tilde{A}$. Further details are provided in the Appendix (ref).
Figure (ref) depicts the impulse response functions of the real GDP, Fed funds rate, unemployment rate, PCE index, nonresidential investment, and real wage to a one standard deviation demand and investment shocks. The impulse responses are visually identical to those obtained using the elliptical slice sampler of arias2025elliptical under the same prior distribution and identifying restrictions (as in the previous section, we implemented their sampler in R, but report times based on a faster Matlab implementation). This suggests that our sampling algorithm correctly targets the same posterior distribution. In particular, the left panel of Figure (ref) shows that the demand shock causes a temporary increase in the real GDP, prices, and the federal funds rate, but reduces unemployment. The demand shock has no effect on nonresidential investment, and a negative impact on real wages on impact. The right panel depicts impulse responses to an investment shock. In this case, nonresidential investment rises for three quarters following the shock, while the other responses are similar to those observed for a demand shock.
Figure (ref) displays the traceplots of the unnormalized log-posterior and the impact response of GDP to a demand shock. For clarity, only the first 4{,}000 iterations of the sampler are shown. The initial value lies within the typical set of the posterior distribution, and, therefore, no iterations are spent moving the chain into the typical set. Furthermore, the chains appear stationary and exhibit low autocorrelation, indicating a good performance of the sampler. Traceplots of other parameters are similar and are therefore not shown.
Figure (ref) plots the maximum values of the $\hat{R}$ statistic for samples obtained using both samplers. The $\hat{R}$ is calculated for the elements of $\Sigma$ and $Q$ for the elliptical slice sampler, and for the elements of $B$ in the case of the NUTS. As in the previous application, we consider a threshold value of 1.01. As can be seen in the figure, $\hat{R}$ converges much faster when sampling using the NUTS. The elliptical slice sampler requires roughly 370 thousand saved draws to reach the threshold of 1.01, while the NUTS only needs 4 thousand total iterations. Hence, the NUTS takes approximately 5 minutes to obtain an acceptable sample, whereas the elliptical slice sampler requires roughly 46 minutes.
The large number of iterations required by the elliptical slice sampler to converge suggests a high degree of autocorrelation also in this empirical application. Indeed, the chain generated by the elliptical slice sampler is highly autocorrelated, as Figure (ref) shows. The autocorrelation function for the $(1,1)$ element of matrix $Q$ crosses zero only after approximately 700 lags, while the chain for the $(1,1)$ element of $B_{11}$ generated by the NUTS is uncorrelated beyond a few dozen lags.
In this large SVAR, the efficiency per iteration of our method in terms of bulk and tail ESS for individual parameters is significantly better than that of the elliptical slice sampler, as indicated by Figure (ref). After 100{,}000 iterations, the lowest bulk ESS for the parameters in the impact matrix $B$ is 25,700 using the NUTS. In contrast, the slice sampler achieves a minimum ESS of 520 for the elements of $\Sigma$ and $Q$ in 400 thousand iterations. In the case of tail ESS, the corresponding values are 35,650 for the NUTS and 1450 for the slice sampler.
As in Subsection (ref), we also assess convergence using the minimum multivariate ESS of vats2019 (the multivariate ESS is calculated for the parameters in columns 1--10 corresponding to the identified shocks of interest). To achieve a precision of 0.05 for the 95% confidence region, the NUTS requires roughly 7,000 iterations, whereas the elliptical slice sampler requires 126,000 iterations. These figures correspond to runtimes of 9 and 15,5 minutes, respectively.
The time required to obtain 1{,}000 effective samples is 0.019 hours for the NUTS, compared to 0.047 hours for our implementation of the elliptical slice sampler. This is, however, slower than the runtime of less than one minute reported by arias2025elliptical. Our approach is faster than the algorithm of chan2025, although their results are not directly comparable to ours, as they use the asymmetric conjugate prior of chan2022, and only identify eight shocks instead of ten. They report a runtime of 140 minutes to obtain 1,000 draws, and as their draws are independent, this corresponds to an ESS of 1,000.
To demonstrate that our approach can simultaneously incorporate both sign and zero restrictions, we augment the identification scheme with the restriction that labor productivity and utilization-adjusted TFP react on impact only to the technology and oil price shocks (assumed to have a positive impact effect). As discussed in Section (ref), introducing these zero restrictions does not require importance sampling to ensure that the sample is from the correct distribution. Instead, the restricted values are simply fixed in the density function and omitted from the parameter vector. Hence, estimation time is largely unchanged by the addition of these restrictions.
In general, the impulse responses depicted in Figure (ref) remain intact. However, the responses to the labor supply shock shown in Figure (ref) change substantially when the zero restrictions are imposed. Under the baseline identification (left panel), labor productivity responds positively on impact and over the following three quarters, while the impact response of TFP is centered around zero but there is a lot of uncertainty about this estimate. In contrast, under the zero restrictions, the 68% posterior credible intervals for both variables include zero at all horizons.
Estimation of the model with additional exclusion restrictions remains efficient: the algorithm reaches a minimum of 19,700 and 31,900 bulk and tail ESS, respectively, while the maximum $\hat{R}$ crosses the 1.01 threshold after 4 thousand iterations, or 5 minutes. In terms of multivariate ESS, the sampling efficiency comes at 0.019 hours for 1{,}000 effective draws (as before, the effective draws are calculated for the impulse response parameters in columns 1--10 corresponding to the identified shocks of interest). These numbers are close to the case with no zero restrictions incorporated in the model.
Tight identifying inequality restrictions may render conventional accept-reject algorithms inefficient and in some cases infeasible for inference in SVAR models. In large-scale models that are becoming increasingly popular in empirical research, the high-dimensional parameter space and the large number of imposed restrictions can give rise to problems, even when individual restrictions are relatively weak, as their cumulative effect may still be substantial. In this paper, we propose an approach based on imposing inequality restrictions through a reparameterization of the structural model by continuously differentiable mappings. This approach accommodates various kinds of inequality restrictions, including shape and ranking restrictions on impulse responses at multiple horizons, as well as bounds on elasticities. Moreover, in contrast to the approaches of arias2025elliptical, chan2025, read2025, our framework allows for a straightforward implementation of exclusion restrictions, with no significant increase in the speed of the sampler. As an additional advantage, the prior can be specified very flexibly, as our approach does not rely on direct sampling procedures, albeit this usually entails an increase in computation time.
While our framework is very flexible, it necessitates relying on MCMC methods because the posterior distribution of the parameters does not belong to a standard family, and we show that the NUTS algorithm, a variant of the HMC sampler, performs extremely effectively in exploring the posterior. However, each iteration of the NUTS requires repeated evaluation of the log-posterior and its gradient, which makes sampling computationally intensive. To address this, we show how the normal-inverse Wishart posterior density can be evaluated efficiently under the proposed reparameterization, facilitating inference even in high-dimensional parameter spaces.
Recently, kitagawa2025 have raised doubts about the effectiveness of MCMC algorithms in Bayesian analysis under set identification, but our simulation experiments indicate that the implementation in the Stan software package stanreference is effective in exploring high-dimensional, set-identified posteriors with flat directions and produces well-mixed chains even when the identified set is large.
Empirical illustrations with small-scale and large-scale SVAR models highlight the good performance of our approach. In particular, it exhibits clearly faster mixing than the Gibbs sampler of arias2025elliptical, while it turns out that the latter can generate Markov chains with extremely persistent serial correlation even at hundreds of lags, resulting in slow mixing and a very small effective sample size.
\printbibliography