Two Gaussians, Too Many: A bootstrap-based approach to assess identifiability in non-Gaussian structural Vector Autoregressions

Two Gaussians, Too Many: A bootstrap-based approach to assess identifiability in non-Gaussian Structural Vector Autoregressions

\setcounter{page}{1}

\dedicatory{Department of Economics, University of Bologna\\ July 19, 2026}

Introduction

How do we identify the latent structural drivers of the economy? This question has long been a cornerstone of research in empirical macroeconometrics. A prevalent approach with structural vector autoregressions (SVARs) involves linking observed macroeconomic variables to structural shocks. Since the covariance matrix of the reduced-form innovations alone is not sufficient, additional restrictions are necessary to recover the structural shocks from the reduced-form innovations. Traditional approaches to identification rely on economic theory, which include short-run exclusion restrictions blanchard_empirical_2002, cointegration-based constraints blanchard_dynamic_1989, sign-restrictions uhlig_what_2005, or employ external information as instruments for exogenous sources of variations stock_identification_2018. Another strand of literature exploits statistical properties of the structural shocks, while avoiding any economic restrictions, to achieve identification. Derived from the Independent Component Analysis (ICA) literature, see comon_independent_1994, hyvarinen_fast_1999, identification based on non-Gaussianity of the structural shocks focuses on the information beyond the second-order moments of the shocks, such as skewness and excess kurtosis, see lanne_identification_2017, gourieroux_statistical_2017, keweloh_generalized_2021. Moreover, given the nature of statistical identification, we still require economic labeling to gain meaningful interpretation.

If the structural shocks are mutually independent and at most one of them follows a Gaussian distribution, the shocks can be identified up to sign and permutation indeterminacies. However, if two or more independent shocks are Gaussian, rotations within the Gaussian subspace leave the joint distribution of the reduced-form innovations unchanged gourieroux_statistical_2017. This implies that the structural model is not uniquely determined, which makes ensuring the validity of these assumptions essential for valid asymptotic inference. In empirical research, most studies verify these identification conditions indirectly through univariate tests of Gaussianity on the reduced-form innovations, see keweloh_generalized_2021, lanne_identification_2017, among others. We show that this conventional strategy of tests of Gaussianity is insufficient to ensure the validity of these conditions. Since the reduced-form innovations are simply linear combinations of the structural shocks, they remain non-Gaussian also with multiple Gaussian shocks. This misleads the practitioner to presume the validity of the identification conditions when in fact, they are violated.

Instead, we propose a novel bootstrap-based approach to assess the identifiability of non-Gaussian SVARs. The proposed method, under certain regularity conditions, tests the null hypothesis that among the independent shocks, there is at most one Gaussian shock, against the alternative that there are two or more Gaussian structural shocks. The method measures the divergence of the bootstrap distribution of the impact matrix from its asymptotic benchmark: the two coincide under valid identification and diverge systematically under its failure. We build on a growing literature in bootstrap-based diagnostic testing: angelini_bootstrap_2022, angelini_identification_2024 show that bootstrap can be used to detect weak identification and model misspecification in dynamic state-space models, and SVARs identified with external instruments. Moreover, cavaliere_bootstrap_2025 provide a framework for diagnosing general mis-specifications such as non-stationarity and weak instruments, among others, which we specialize and make operational for non-Gaussian identification.

Specifically, we make three contributions: First, we turn the bootstrap-validity principle into an identification diagnostic for non-Gaussian SVARs: under valid identification and certain regularity conditions, the conditional bootstrap distribution of the impact matrix, established within NGML estimation of lanne_identification_2017, is asymptotically standard normal, so the check reduces to a normality test on the bootstrap replications. Second, we prove that this validity extends to the single-Gaussian boundary of the null, where the Gaussian shock's shape parameter lies on the boundary and the full-parameter Fisher information is singular, by showing the impact matrix estimator remains consistent and asymptotically normal; the diagnostic is thus valid across the entire null.\footnote{Beyond the diagnostic, this carries a practical benefit: a practitioner who suspects a Gaussian shock need not determine which shock it is, nor adopt a different specification; estimating all shocks under the NIG family, which nests the Gaussian at its boundary, delivers consistent, asymptotically normal impact-matrix estimates and valid inference regardless of which shock, if any, is Gaussian.} Third, a univariate version of the diagnostic localizes, in large samples, which shocks are responsible for an identification failure, extending the procedure to partial identification maxand_identification_2020. Under the alternative of two or more Gaussian shocks the bootstrap distribution diverges from its limit and the test is consistent, with non-negligible power.

Apart from the incomplete verification through residuals' non-normality, current literature for verifying the identification conditions include a moments-based rank test on reduced-form innovations guay_identification_2021 or a specification test on the estimated shocks amengual_moment_2022,amengual_specification_2024 which may induce severe pre-testing bias, since the finite-sample distribution of a post-selection estimator is not uniformly close to its asymptotic benchmark danilov_harm_2004,roth_pretest_2022,leeb_model_2005. Crucially, we show that this bootstrap-based approach does not induce any pre-testing bias in subsequent inference. This is because, unlike standard bootstrap asymptotic theory where the number of bootstrap replications ($M$) and the sample size ($T$) diverge to infinity sequentially, first $M \rightarrow \infty$, and then $T \rightarrow \infty$, in this framework they diverge jointly i.e., $M,T \rightarrow \infty$ and $M/T \rightarrow 0$ at an appropriate rate. The joint divergence ensures that under the null hypothesis of valid specification, the test statistic, conditional on the data, is asymptotically pivotal and does not distort subsequent inference. This idea forms the motivation for this method: the bootstrap is employed not merely as a tool to gain asymptotic refinements in finite samples over first-order approximations of distributions, but as an instrument capable of detecting identification failures in non-Gaussian SVARs without distorting post-test inference.

We validate the theoretical results with a Monte Carlo study across three different sets of distributional specifications for the structural shocks: (i) ICA 0, where all three shocks are non-Gaussian\footnote{Here, non-Gaussian refers to the NIG distribution, nonetheless this diagnostic procedure is also valid for other heavy-tailed distributions, such as Student-$t$ distribution.}, (ii) ICA 1, where one shock is Gaussian and the other two are non-Gaussian, and (iii) No ICA, where two of the shocks are Gaussian and one is non-Gaussian. The estimates of the impact matrix and their probability coverages confirm consistency and asymptotic normality under ICA 0 and ICA 1, while under the identification failure of No ICA, the estimates and their bootstrap analogs diverge sharply from Gaussianity. The bootstrap diagnostic test of normality, under the appropriate rate of joint divergence of $M$ and $T$, achieves correct empirical size under the null and increasing power in $T$ against the alternative of invalid identification; $M$ governs the size-power trade off in finite samples. In contrast, conventional residual-based normality tests cannot discriminate between valid and invalid identification, due to \textit{extra} Gaussian shocks. The bootstrap-based approach thus provides a reliable and implementable tool for evaluating identification conditions for non-Gaussian SVARs in empirical applications. Furthermore, under weak identification where one structural shock is \textit{nearly} Gaussian in the presence of another Gaussian shock, the probability coverages (conditional on the rejection of the residuals' normality test) of the estimates are incoherent i.e., either severely under covered or artificially inflated. In contrast, the probability coverages conditioned on the validity of the bootstrap diagnostic remain coherent and consistently track their unconditional analogs.

We further demonstrate the potential of our diagnostic in an empirical illustration, where we investigate whether uncertainty is an exogenous driver of business cycle fluctuations or an endogenous response to the macroeconomic drivers. We consider a trivariate SVAR model with measures of macro uncertainty, real activity and financial uncertainty from the U.S. macroeconomic data. The bootstrap diagnostic does not suggest rejection of the null hypothesis of non-Gaussianity of the independent shocks, across different significance levels and choices of $M$. This supports the evidence of ludvigson_uncertainty_2021 where they document significant excess kurtosis in their identified shocks. Furthermore, we find that a negative shock to real activity causes a significant and persistent increase in macro uncertainty, whereas the impact on financial uncertainty is negligible. We also find that there is a persistent negative effect on real activity (albeit with positive immediate impact) from an increase in macro uncertainty, whereas the effect of financial uncertainty on real activity is insignificant on impact and becomes negative only with a delay. This is partly in contrast with the findings in ludvigson_uncertainty_2021, where financial uncertainty is found as an exogenous driver that causes significant impact on real activity. However, our findings are in line with the zero restrictions of angelini_uncertainty_2019 where the contemporaneous impact of an exogenous impulse in financial uncertainty on real activity, and vice versa, are set to zero.

The remainder of the paper is organized as follows. Section (ref) elaborates on identification in non-Gaussian SVARs, outlining the NGML estimation procedure and its asymptotic properties, and provides the asymptotic distribution of the bootstrap estimator under the null of valid identification, Section (ref) introduces the bootstrap diagnostic for assessing identification validity, Section (ref) presents Monte Carlo results, Section (ref) discusses the case of weak identification where apart from one existing Gaussian shock, there is another nearly Gaussian shock, Section (ref) details an empirical illustration and Section (ref) concludes. Appendix contains notational preliminaries, proofs of the main results and additional Monte Carlo evidence.

Identification in non-Gaussian SVARs

This section elaborates identification in non-Gaussian SVARs through the NGML estimator of lanne_identification_2017. The structural shocks follow a Normal-inverse Gaussian distribution, which nests the Gaussian distribution as a limiting case. We extend the consistency and asymptotic normality of the structural impact matrix to the boundary case in which one shock lies exactly at the Gaussian limit (Proposition (ref)). Furthermore, under the null of valid identification, we also provide the limiting distribution of the bootstrap estimator of the impact matrix (Proposition (ref)).\\ Consider a VAR($p$) model in its companion form:

where,

and $u_t = B\varepsilon_t$ with,

Here $y_t$ is an $n \times 1$ vector of endogenous variables; $\Pi_i, \ i = 1, 2,\cdots, p$ are $n \times n$ autoregressive parameter matrices; $u_t$ is a vector of reduced-form innovations, $B$ is a non-singular, $n \times n$ impact/ \enquote{mixing} matrix. It links the structural shocks, $\varepsilon_t$, an $n \times 1$ vector, to the reduced-form VAR innovations. With $u_t = B\varepsilon_t$, we note that the covariance matrix of the reduced-form innovations, $\mathbb{E}(u_t u_t') = B \Sigma_\varepsilon B'$. The VAR model (ref) is assumed to be stable i.e., the companion matrix $\mathbf{\Pi}$ satisfies:

Beyond a stable VAR, we assume that the structural shocks are mutually independent and non-Gaussian, and exploit the information in addition to the unconditional variances of the reduced-form innovations i.e., skewness and/or excess kurtosis of the structural shocks. This yields identification of the impact matrix up to sign and column-permutation indeterminacies, without the need to impose any additional identifying restrictions. This idea follows from the darmois_analyse_1953, skitovich_linear_1953 theorem, which states that if $X = (X_1, X_2, \cdots, X_n)$ are independent random variables and $\alpha'X$ and $\beta'X$ are independent for $\alpha_i, \beta_i \in \mathbb{R}/ \{0\}$, then all $X_i, \ i = 1,2, \cdots, n$ are Gaussian. Extending this result and formalizing in Independent Component Analysis (ICA) literature, comon_independent_1994 shows that we can extract a unique decomposition of $u_t = B\varepsilon_t$ up to sign and column permutation, if the $\varepsilon_t$ has at most one Gaussian component and all components are independent\footnote{It should be noted that the assumption of independence of structural shocks is more restrictive than uncorrelatedness. This is because it is not guaranteed to recover independent structural shocks from non-Gaussian reduced-form innovations solely through linear transformation kilian_structural_2017}. However, recent literature shows that it is not necessary to assume complete independence of the shocks to achieve identification. jarocinski_estimating_2024 relaxes the assumption of complete independence and assumes that the shocks are drawn from partially dependent multivariate $t$-distributions. Other moment-based estimators, instead of independence, assume restrictions on higher order moments and co-moments of the shocks, see guay_identification_2021, lanne_identifying_2023, and/or cumulant tensors mesters_non-independent_2024.

NGML estimation:

Apart from non-parametric algorithms, see hyvarinen_independent_2000, hyvarinen_independent_2013 for a detailed review, other likelihood-based approaches can be used for estimation. With NGML we assume that the structural shocks, $\varepsilon_{it}$, are i.i.d.\ sequences\footnote{lanne_identification_2017 relax the temporal independence assumption to temporal uncorrelatedness, allowing for conditionally heteroscedasticity in shocks.}, with at most one sequence following Gaussian distribution.

The densities of the structural shocks may depend on their individual set of parameters, $\theta_i$. The maximum likelihood estimator is defined in terms of the density $f_{i,\sigma_i}(x;\theta_i)$ such that $f_{i,\sigma_i}(x;\theta_i) \coloneqq \sigma_i^{-1}f_i(\sigma_i^{-1}x;\theta_i)$. As noted before, ICA ensures identification of the impact matrix up to sign and column permutations. To alleviate this indeterminacy, we follow a series of transformations, see ilmonen_semiparametrically_2011, where $B$ is column-wise normalized (with the corresponding standard deviation of the structural shock, $\sigma_i$) to have a unitary diagonal, thereby removing the sign and column-order indeterminacies. Nonetheless, we still require economic labeling to provide interpretation to the statistically-identified shocks.

In our framework, the structural shocks are assumed to be from the Normal-inverse Gaussian (NIG) distribution, a member of the generalized hyperbolic distribution family. It is suitable for modeling processes where the probability of obtaining extreme values is higher, vis-à-vis the normal distribution schrodinger1915theorie, and jarocinski_estimating_2024, andrade_higher-order_2025,ludvigson_uncertainty_2021 note non-trivial excess kurtosis in the distribution of macroeconomic shocks, and their proxies respectively\footnote{jarocinski_estimating_2024, lanne_identification_2017 and other studies employ independent Student-$t$ distributions, another member of generalized hyperbolic distribution family, to capture the excess kurtosis (and skewness), however NIG distribution allows for more flexible heavy-tailed processes lahcene_extended_2019.}. The NIG distribution is defined by four parameters: $\alpha$ (tail heaviness), $\gamma$ (asymmetry), $\delta$ (scale) and $\mu$ (location), and is denoted as $NIG(\alpha, \gamma, \delta, \mu)$. Its density function is given by:

where, $K_1$ is the modified Bessel function of the second kind with index 1. The first two central moments of NIG distribution are, $\mu + \frac{\delta\gamma} {\sqrt{\alpha^2 - \gamma^2}}$ and $\frac{\delta\alpha^2}{(\alpha^2 - \gamma^2)^{3/2}}$, respectively. Moreover, given its flexibility, it nests the standard normal distribution when $\gamma = 0, \ \delta = \sigma^2\alpha$ and $\alpha \rightarrow \infty$.

Standard asymptotic-normality results for these estimators, whether the full-parametric NGML of lanne_identification_2017 or the pseudo-ML of gourieroux_statistical_2017, are derived under a non-singular information matrix. This condition is violated at the single-Gaussian boundary, where the estimated shape parameter of the Gaussian shock lies on the edge of the parameter space and its information block vanishes, even though the impact matrix itself remains identified. We show that the impact matrix still remains consistent and asymptotically normal, which eventually allows the bootstrap diagnostic to be consistent across the null hypothesis of at most one Gaussian shock.

Before setting up the estimation procedure, we summarize the parameters of the model: with the unit diagonal normalization of $B$, the off-diagonal parameters of $B$ are collected in an $n \times (n-1)$ vector, $\beta \coloneqq \operatorname{vecd}(B), \ \beta \in \mathbb{R}^{n(n-1)}$, where $\operatorname{vecd}(\cdot)$ is the operator that stacks the columns of an $n \times n$ matrix, excluding the diagonal elements from its $\operatorname{vec}(\cdot)$ form. We can collect the parameters of the individual shocks' densities $(\alpha_i, \gamma_i, \delta_i, \mu_i)$ in a vector $\theta_i$ i.e., $\theta_i \coloneqq (\alpha_i, \gamma_i, \delta_i, \mu_i), \ i = 1,2, \cdots,n$. Furthermore, for brevity, all the model parameters (including the autoregressive coefficients) are collected in a vector $\lambda \coloneqq (\Pi,\beta,\sigma,\theta)$, where $\Pi = \operatorname{vec}(\Pi_1, \Pi_2, \cdots, \Pi_p)$ are the vectorized-autoregressive coefficients, $\theta = (\theta_1', \theta_2', \cdots, \theta_n')'$, where each $\theta_i$ is a $4 \times 1$ vector, are the distribution parameters and $\sigma = (\sigma_1, \cdots, \sigma_n)$ are the standard deviations of the structural shocks, $i = 1, 2, \cdots, n$. We denote $\lambda_0$ as the true parameter vector.

NGML is a two-step estimation procedure\footnote{lanne_identification_2017 also consider a three-step estimation method, for large $n$ and small $T$, provided the distribution of the underlying $\varepsilon_{i,t}$ is symmetric. This estimator is asymptotically efficient.}. First, the least squares estimates of the reduced-form VAR innovations, $\hat{u}_t$, are obtained. With the estimated $\hat{u}_t$, the log-likelihood:

where,

is maximized with respect to the parameter vector $(\beta',\sigma,\theta)'$. Here, $\imath_j$ is the $j^{th}$ unit vector i.e., $\imath_j = (0,\cdots,1,\cdots,0)'$, where $1$ is in the $j^{th}$ position. The estimator is given by:

The estimated structural shocks are obtained by multiplying a transposed $j^{th}$ unit vector, $\imath_j'$,

Here $B(\hat{\beta}_T)$ is the estimated (column-normalized) matrix $B$ as a function of its off-diagonal elements $\hat{\beta}_T$. Once the structural shocks are identified and estimated, we can obtain the estimates of IRFs, given by:

where, $R \coloneqq (I_n , \mathbf{0}_{n \times n(p-1)})$ is a selection matrix, $\hat{\mathbf{\Pi}}$ is the companion form of the estimated autoregressive coefficients, $\hat{\Pi}_1, \hat{\Pi}_2, \cdots, \hat{\Pi}_p$.

The instantaneous responses of one standard deviation shocks are given by $B(\hat{\beta}_T)\operatorname{diag}(\hat{\sigma}_T)$. Henceforth, for brevity, we denote the impact matrix to one standard deviation shocks as:

We summarize the set of regularity conditions on the non-Gaussian SVAR model which allow for standard asymptotic and bootstrap inference for NGML estimates in the Appendix (ref). These regularity conditions include the stability conditions for the SVAR model, the true parameter values lie in a compact, permissible parameter space, conventional differentiability assumptions on the density functions, suitable integrability conditions ensuring the score function has zero mean and finite covariance matrix when evaluated at the true parameter value, and the covariance matrix of the limiting distribution of the ML estimator is positive definite.

Consistency and Asymptotic Normality:

Consider the model representation in equation (ref), satisfying the regularity conditions ((ref)) with Assumptions ((ref)) - ((ref)), then the NGML estimator, $\hat{\lambda}_T = (\hat{\Pi}_T, \hat{\beta}_T, \hat{\sigma}_T, \hat{\theta}_T)$, obtained by maximising the log-likelihood function (ref) and (ref), satisfies\footnote{For detailed proofs, see lanne_identification_2017, Theorem 1. When one structural shock is Gaussian, (ref) holds for the structural sub-vector $\psi=(\operatorname{vec}(\Pi),\beta,\sigma)$ with the robust covariance of Proposition (ref) in place of $\mathcal{I}(\lambda_0)^{-1}$; the impact-matrix estimator $\hat{B}_T$ remains asymptotically normal, which is the object of interest.}:

where $\xrightarrow{d}$ denotes convergence in distribution, $\mathcal{I}(\lambda_0)$ is the Fisher information matrix, defined as, $\mathcal{I}(\lambda_0) \coloneq -\mathbb{E}[\nabla^2 \ell_{\lambda\lambda,t}(\lambda_0)]$, where $\nabla^2 \ell_{\lambda\lambda,t}(\lambda_0) = {\partial^2 \ell_t(\lambda_0)}/{\partial \lambda \partial \lambda'}$. Furthermore, a consistent estimator of the asymptotic covariance matrix can be obtained by the ML estimator, $\hat{\lambda}_T$ and the Hessian matrix of the log-likelihood function i.e.,

We extend these results of consistency and asymptotic normality of the estimator to the case of single, exact Gaussian shock. The identification conditions permit at most one Gaussian shock, and the boundary of this null is precisely where the standard asymptotics require attention. When one shock's fitted NIG density reaches its Gaussian limit ($\alpha_k\to\infty$, equivalently $\rho_k=1/\alpha_k=0$), that shock's shape parameters are no longer identified: the corresponding block of the Fisher information collapses, so the full-parameter information matrix $\mathcal{I}(\lambda_0)$ is singular. This boundary is not an isolated knife-edge; the full-parameter information is already ill-conditioned whenever a single shock is merely {near}-Gaussian. The exact-Gaussian point is simply the limit of this near-singular region. The case is therefore a genuine part of the null hypothesis over which the diagnostic must remain valid.

Proposition (ref) establishes that this ill-conditioning is confined to the Gaussian shock's own shape parameters and does {not} propagate to the impact matrix. Although the full information is singular at the boundary, the autoregressive-impact parameters $\psi=(\operatorname{vec}(\Pi),\beta,\sigma)$ retain a non-singular information $\mathcal{I}_\psi$. Since $\mathcal{I}_\psi$ is non-singular at the boundary, the impact-matrix estimator $\hat{B}_T$ is uniformly consistent and asymptotically normal as a shock ranges from strongly non-Gaussian down to and including the Gaussian limit. Impact-matrix inference is thus stable across the entire null, not merely at its interior. The estimator loses identification only when a second shock also approaches Gaussianity: the rotational near-indeterminacy of the weak-identification regime studied in Section (ref) which lies in the alternative hypothesis.

Proposition (ref) has an implication for practitioners that is independent of the diagnostic. A practitioner who believes one structural shock may be Gaussian need not decide which one, nor switch to a partly Gaussian or otherwise restricted specification. Estimating all $n$ shocks under the NIG family, which nests the Gaussian at the boundary, leaves the impact-matrix estimator $\hat{B}_T$ consistent and asymptotically normal, with valid asymptotic standard errors from $\mathcal{V}_B$ in (ref), even when a shock lies at the boundary of the parameter space.

By delta method, see kilian_structural_2017, bruggemann_inference_2016, the covariance matrix of $\hat{B}_T = B(\hat{\beta}_T)\operatorname{diag}(\hat{\sigma}_T)$, denoted as $\hat{\Sigma}_{{B}_T}$, can be obtained as: for $i = 1,2,\cdots,n$ and $j = 1,2,\cdots,n$,

where $B(\hat{\beta}_T)_{i,j}$ is the $(i,j)^{th}$ element of the matrix $B(\hat{\beta}_T)$. Throughout the diagnostic, $\hat{\Sigma}_{B_T}$ denotes the robust estimator of $\mathcal{V}_B$ (Proposition (ref)); when all shocks are non-Gaussian, it coincides with the efficient covariance (ref).

Bootstrap Inference:

Now we discuss the bootstrap estimator, $\hat{B}_T^*$ and provide its limiting distribution under the null of valid identification conditions. We strictly follow the residual-based MBB algorithm of bruggemann_inference_2016, see Appendix section (ref), where the residual-based MBB yields asymptotically valid inference for the reduced-form VAR under conditional heteroscedasticity of unknown form (the NGML score conditions being verified under the i.i.d.\ design; see Appendix (ref)). In the Monte Carlo simulations, we consider structural shocks are drawn i.i.d.\ from a non-Gaussian distribution which allows for parametric i.i.d.\ bootstrap algorithms as well. However in empirical applications with small samples, the general setup of residual-based MBB and drawing sample from the empirical distribution of the residuals is expected to be more robust. In the Monte Carlo simulations, with valid specifications, the results are similar with both bootstrap algorithms.

Let $\hat{B}_T^{\ast}=B(\hat{\beta}_T^{\ast})\operatorname{diag}(\hat{\sigma}_T^{\ast})$ denote the bootstrap analog of the impact-matrix estimator $\hat{B}_T$, obtained from Algorithm (ref), and define the studentized bootstrap statistic

where $\hat{\Sigma}_{B_T}$ is a consistent estimator of the covariance $\mathcal{V}_B$ of $T^{1/2}(\hat{B}_T-B_0)$, and $B_0$ is the true impact matrix. The distribution of $Q_T^{\ast}$ conditional on the observed data is denoted by $\hat{G}_T^{\ast}(\cdot)$. It is used to approximate the sampling distribution $G_T(\cdot)$ of the corresponding sample statistic

which, under valid identification, is asymptotically standard normal, see (ref). $\hat{G}_T^{\ast}$ can be computed with arbitrary precision from the bootstrap replicates $\{Q_{T,b}^{\ast}\}_{b=1}^{N}$ of Algorithm (ref) by the empirical distribution function

where $\mathbf{1}(\cdot)$ is the indicator function by the Glivenko-Cantelli theorem, $\sup_x\big|\hat{G}_{T,N}^{\ast}(x)-\hat{G}_T^{\ast}(x)\big| \xrightarrow{\text{a.s.}}0$ conditionally on the data as $N\to\infty$. The following proposition and corollary establish the asymptotic validity of this bootstrap approximation; under valid identification $\hat{G}_T^{\ast}$ converges to the standard normal, the property exploited by the diagnostic of Section (ref).

Proposition (ref) is stated for the full vector $\lambda$ under the regularity conditions; when one structural shock is Gaussian these hold for the structural sub-vector $\psi$ rather than for $\lambda$ (Proposition (ref)), and Corollary (ref) is the specialization relevant to the diagnostic. Here, $\xrightarrow{p^*}_p$ and $\xrightarrow{d^*}_p$ denote convergence in probability and in distribution, respectively, conditional on the observed data. For a more detailed discussion on the notation, we refer the reader to Appendix (ref).

Bootstrap Diagnostic

This section develops the bootstrap diagnostic for the identification conditions of a non-Gaussian SVAR. Proposition (ref) and Corollary (ref) established that, under valid identification, the studentized statistic $Q_T^{\ast}$ is asymptotically standard normal; equivalently, its bootstrap distribution satisfies

where $\Phi$ is the standard normal distribution function. The diagnostic exploits the converse: when identification fails, $\hat{G}_T^{\ast}$ departs from normality, so a test of the normality of the bootstrap replications detects the failure.

The mechanism underlying the diagnostic is as follows. Under valid identification and the maintained regularity conditions, the log-likelihood admits a unique maximum and the NGML estimator concentrates at rate $T^{1/2}$; the standardized bootstrap replications are then asymptotically standard normal. When two or more shocks are Gaussian, the log-likelihood is flat along the rotations of the Gaussian subspace, the affected columns of $\hat{B}_T$ are not identified, and the estimator fails to concentrate (Remark (ref)). The replications then follow the non-degenerate distribution induced by the unidentified rotation, so they do not converge to normality (Lemma (ref)). Conditional on the maintained assumptions of stable VAR and shocks' independence and regularity conditions, a rejection is therefore attributed to identification failure, that is, the presence of two or more Gaussian shocks.

To measure the deviation of $\hat{G}^{\ast}_T$ from $\Phi$, we use the studentized statistic

where $\hat{\Omega}_T(x)$ is a consistent estimator of $\Omega_T(x) \coloneqq\hat{G}_T^{\ast}(x)\bigl(1-\hat{G}_T^{\ast}(x)\bigr)$.

At finite $T$, the bootstrap distribution $\hat{G}_T^{\ast}$ is only approximately standard normal: under the null it differs from $\Phi$ by an Edgeworth term of order $T^{-\rho}$ (Condition (ref)), and this discrepancy is a function of the observed data. The statistic decomposes into a resampling term, $M^{1/2}\hat{\Omega}_T^{-1/2}(\hat{G}_{T,M}^{\ast}-\hat{G}_T^{\ast})$, and a centering term, $M^{1/2}\hat{\Omega}_T^{-1/2}(\hat{G}_T^{\ast}-\Phi)$ in (ref), which makes the problem transparent. If $M\to\infty$ for fixed $T$, the empirical distribution resolves $\hat{G}_T^{\ast}$ exactly, so the test effectively checks whether $\hat{G}_T^{\ast}=\Phi$ exactly, which may hold only for a finite but very large $T$. The centering term is then of order $M^{1/2}T^{-\rho}$ and diverges, so the statistic reflects the $O(T^{-\rho})$ approximation error rather than an identification failure, and the test is too conservative under the null. Moreover, since $\hat{G}_T^{\ast}$ is a function of the data, conditioning subsequent inference on such a test induces pre-testing bias roth_pretest_2022, which can be especially difficult to control because the impulse responses of interest are non-linear functions of the estimator.

Following cavaliere_bootstrap_2025, we therefore let $M$ and $T$ diverge jointly, with $MT^{-2\rho}=o_p(1)$ for some $\rho>0$.\footnote{Under the sequential regime in which $M\to\infty$ after $T\to\infty$ the pre-testing bias likewise vanishes asymptotically; we adopt the joint regime because the practitioner controls $M$ more readily than the sample size.} This requirement is exactly $M^{1/2}T^{-\rho}\to0$, which forces the centering term to vanish. The limiting randomness of $d^{\ast}_{T,M}$ then derives solely from the resampling term, which is free of the data, so the statistic is pivotal and the test attains correct asymptotic size (Remark (ref)).

Monte Carlo Study

SVAR Design:

In this section, we demonstrate the performance of the bootstrap diagnostic test in the detection of invalidity of identification in non-Gaussian SVARs due to the presence of two or more Gaussian shocks. Consider the SVAR model in (ref), with $n=3$ and $p=1$ i.e., a trivariate VAR(1) model:

where $Y_t = (Y_{1,t}, Y_{2,t}, Y_{3,t})'$, $\Pi_1$ is a $3 \times 3$ matrix of autoregressive coefficients, $u_t = (u_{1,t}, u_{2,t}, u_{3,t})'$ is the vector of reduced-form innovations, $B$ is the $3 \times 3$ impact matrix, and $\varepsilon_t = (\varepsilon_{1,t}, \varepsilon_{2,t}, \varepsilon_{3,t})'$ is the vector of structural shocks. And,

As noted in equation (ref), the impact matrix $B$ is assumed to be full rank and column-wise normalized to have a unitary diagonal, $\beta_0$ is a vector containing the off-diagonal elements of $B$, of dimension $n(n-1) = 6$. Hence, $B = B(\beta_0)$ is parameterized as a function of $\beta_0$. We will focus on the estimates of $\hat{B}_T$: the instantaneous impact matrix of a one standard deviation shocks, given by:

where $\hat{\sigma}_T$ is the estimate of the standard deviations of structural shocks.

Distributional Specifications:

We consider three different specifications for the structural shocks' distributions. Specifically\footnote{We also consider a set of specifications where the non-Gaussian structural shocks are drawn from the Student-$t$ distribution, the results are qualitatively similar.}:

Though the NIG distribution is defined by four parameters, we fix the location and asymmetry parameters respectively, $\mu = \gamma = 0$, to center the respective structural shocks around zero\footnote{The choice to restrict the asymmetry parameter $\gamma$ to zero is not consequential, but it allows direct comparison with results from Student-$t$ distributed shocks. The results are robust when $\gamma$ is allowed to be non-zero.}. Hence, the above-mentioned parameterization is of the form: $NIG(\alpha, \delta)$, see (ref). In each specification, the structural shocks have different levels of standard deviations and kurtosis\footnote{The excess kurtosis of the shocks is: $\textbf{ICA 0} \sim (20, 6, 2.85)$, $\textbf{ICA 1} \sim (20, 6, 0)$ and $\textbf{No ICA} \sim (20, 0, 0)$.}. The NIG distribution has finite moments of all orders for $|\gamma_i|<\alpha_i$ barndorff-nielsen_normal_1997, so the maintained assumption $\mathbb{E}\|\nabla_\lambda\ell_t\|^{2+\delta}<\infty$ holds. Since our object of interest is the impact matrix, computed from the estimates of $\hat{\beta}_T$ and $\hat{\sigma}_T$, the standard deviations $\sigma_i$ of the structural shocks $\varepsilon_t$ are:

From the Assumptions (ref) and (ref), we can infer that the first two specifications are valid. However, in the third specification, two of the independent structural shocks are Gaussian, thereby violating the identifying conditions.

Insufficient verification with reduced-form innovations:

Given the relation $u_t = B \varepsilon_t$, it can be tempting to verify the assumptions of non-Gaussianity of structural shocks by testing the estimated\footnote{Unlike the moment-based specification tests of amengual_moment_2022,amengual_specification_2024, which test the fitted parametric shock distribution and can induce pre-testing bias, the proposed diagnostic tests the bootstrap distribution of the impact matrix estimator, and is asymptotically ancillary (Remark (ref)), so conditioning on it does not distort subsequent inference.} reduced-form innovations $\hat{u}_t$ for non-Gaussianity. Indeed, many studies, including lanne_identification_2017, andrade_higher-order_2025, jarocinski_deconstructing_2020, among others, have verified the non-Gaussianity of reduced-form innovations as a pre-test for valid identification. Under the null of non-Gaussian structural shocks, their linear combination would naturally allow the reduced-form innovations to inherit the non-Gaussianity. However as noted earlier, given that valid identification allows at most one Gaussian shock, the presence of non-Gaussianity of reduced-form innovations does not preclude the presence of more than one Gaussian shocks. The following table illustrates this point.

Tables (ref) and (ref) report the empirical rejection frequencies of the multivariate and univariate doornik_omnibus_2008 and jarque_test_1987 tests for normality of $\hat{u}_t$. Table (ref) shows the insufficiency of verifying the assumption of at most one Gaussian shock through testing the reduced-form innovations. Both specifications, ICA 1 and No ICA, predominantly reject both univariate tests of normality across all three residual series. Moreover, with increasing sample sizes, the increase in rejection frequencies of the null hypothesis of Gaussianity exacerbates the insufficiency to verify the assumptions.

Bootstrap Probability Coverages:

Tables (ref) and (ref) report the estimates of the true parameter on-impact matrix $B_{0}$, from the two distributional specifications ICA 1 and No ICA respectively, with sample size $T = 100, 500$, Monte Carlo simulations $N_S = 500$ and residual-based MB bootstrap replications $N = 999$. We also report the empirical coverage probabilities of the nominal 90% confidence intervals (CIs) for each specification\footnote{Since we are concerned with invalidity of identification due to extra Gaussian shock, in this section we show results pertaining only to specifications ICA 1 and No ICA. Results for the specification ICA 0, where all three structural shocks are non-Gaussian can be found in Appendix (ref).}.

Table (ref) shows that for the specification with one Gaussian shock, the estimates $\hat{B}_T$ are consistent, and the empirical coverage probabilities of the nominal 90% CIs are close to the nominal level. The bootstrap standard errors also closely match those from NGML estimation, confirming the theoretical bootstrap validity established in Proposition (ref) under the regularity conditions ensuring asymptotic consistency of the NGML estimator.

In contrast, Table (ref) shows that when two Gaussian shocks are present, the estimates of $\hat{B}_T$ become inconsistent, and the empirical CI coverages fall well below the nominal level. Moreover, for large samples ($T = 500$), the estimates corresponding to the instantaneous impact of the single non-Gaussian shock in the specification No ICA ($\varepsilon_{1,t}$), remain consistent and comparable to those obtained from valid specifications (ICA 0 and ICA 1). The studentized-bootstrap CI coverages also align with their asymptotic analogs, indicating asymptotic validity of these individual estimates.

These findings are consistent with the literature on partial identification of non-Gaussian structural shocks in the presence of multiple Gaussian shocks maxand_identification_2020. However, inference from such partially identified models requires prior knowledge of which shocks are Gaussian, which is generally infeasible in empirical settings. Nonetheless, the bootstrap diagnostic can detect such misspecification, in large samples, allowing practitioners to identify the Gaussian shocks in the No ICA specification through univariate diagnostics since their corresponding estimates' distribution deviate from their asymptotic distributions.

Bootstrap Diagnostic:

We consider the multivariate and univariate (over each element of the estimated on-impact matrix) normality test of Doornik-Hansen for the sequences, ${\{\hat{B}_{T,1}^*, \hat{B}_{T,2}^*, \cdots, \hat{B}_{T,M}^*\}}_{s}, \ s = 1, 2, \cdots, S, \ b = 1,2,\cdots,M$, where $M$ is the number of bootstrap replications considered for the number of tests $S$. As noted earlier, we choose, (i.) $M = (1/3)T^{3/5}$ and (ii.) $M = (1/2) T^{3/5}$, to maintain a balance between power and size control in finite samples, while ensuring that $M/T\to 0$. We consider the empirical rejection frequencies of the tests across the whole range of significance level i.e., $\alpha \in [0,1]$, and the empirical distribution function of the $p$-values, $p^*_{M,T,s}$, $s = 1, 2, \cdots, S$. Under the null hypothesis, for continuous test statistics, it is well known that the $p$-values are uniformly distributed i.e., under valid specification, $p^*_{M,T,s} \xrightarrow{d^*}_p \mathbb{U}(0,1)$, as $S \to \infty$ and the average rejection frequencies can be computed as $\pi^*_{M,T,S}(x) = \frac{1}{S} \sum_{s=1}^S \mathbf{1}(p^*_{M,T,s} \leq x)$. Hence, under valid specification of at most one Gaussian shock, $\pi^*_{M,T,S}(x) \xrightarrow{p} x$, as $S \to \infty$. However, under the alternative hypothesis, the $p$-values are stochastically smaller than $\mathbb{U}(0,1)$ and hence, $\pi^*_{M,T,S}(x) > x$ for some $x \in (0,1)$. (hung_behavior_1997, tang_general_2021).

Multivariate Diagnostic:

Figures (ref) show plots of the empirical distribution function of the $p$-values $p^*_{M,T,s}$ of the multivariate normality test\footnote{Empirical rejection frequencies for nominal significance level 5% are also tabulated in Appendix, Table (ref)} for the three specifications (ICA 0, ICA 1 and No ICA), with sample size\footnote{The empirical distribution plots for sample size $T=300$ can be found in Appendix (ref) (Figure (ref))} $T = 500$, the number of bootstrap replications in a sequence $M = (1/2)T^{3/5}$, and the number of sequences $S = 1000$, across Monte Carlo simulations $N_S = 500$. The dashed yellow line indicates the $45^{\circ}$ line, with the $p$-values on $x$-axis and their corresponding empirical rejection frequencies on the $y$-axis. The black (red) line is the median (mean) of the empirical rejection frequencies, while the dark and light shaded areas indicate the $50$ and $90$-percentile intervals, respectively, across Monte Carlo simulations. The plots of the empirical distribution function of the $p$-values reinforce the previous assertion that in specifications with 0 and 1 Gaussian shocks, (ICA 0 and ICA 1), the $p$-values of the multivariate normality test on bootstrap replications are uniformly distributed, whereas under the invalid specification i.e., with extra Gaussian shocks, the $p$-values diverge from their asymptotic distribution\footnote{Under the alternative hypothesis i.e., when there are two Gaussian shocks, the choice of $M$ becomes less relevant as the normality test soundly rejects the null hypothesis, even for unusually large sample sizes of $T = 10,000$, and $M \in [10,1000]$. So the choice of $M$ is made to control size of the test.}. Hence, a multivariate test of normality of the bootstrap replications of the estimated on-impact matrix can detect the invalidity of the identifying assumptions i.e., the presence of excess Gaussian shocks.

Univariate Diagnostic:

The bootstrap diagnostic can also be implemented in a univariate manner, by applying the univariate Doornik-Hansen test for normality\footnote{Results with the univariate Jarque-Bera test for normality can be found in Appendix (ref).} to each element of the estimated impact matrix $\hat{B}_T$. Figures (ref) and (ref) show the empirical distribution function of $p^*_{M,T,s}(x)$, of the tests over bootstrap sequences ${\{\hat{B}_{i,j,T,1}^*, \hat{B}_{i,j,T,2}^*, \cdots, \hat{B}_{i,j,T,M}^*\}}, \ i = 1,2,3, \ j = 1,2,3$ for the specifications ICA 1 and No ICA, with sample size $T = 500$, the number of bootstrap replications in a sequence $M = (1/2)T^{3/5}$, and the number of tests $S = 1000$ across Monte Carlo simulations $N_S = 500$.

We can observe that, under the null of valid specification with 1 Gaussian shock, the empirical distributions of $p^*_{M,T,s}(x)$ are centered around the $45^{\circ}$ line, indicating that the empirical rejections maintain the nominal level of significance and the $p$-values are uniformly distributed. However, under the alternative of invalid specification of No ICA i.e., with 2 Gaussian shocks, the empirical distributions of $p^*_{M,T,s}(x)$ corresponding only to the elements of $\hat{B}_T$ associated with the two Gaussian shocks ($\varepsilon_{2,t}, \varepsilon_{3,t}$) are significantly above the $45^{\circ}$ line, indicating that the test has power against the invalidity that increases with $T$ and $M$. This is consistent with the results of partial identification from Table (ref), where the probability coverages (asymptotic and bootstrap) of estimates corresponding to the one non-Gaussian ($\varepsilon_{1,t}$) closely align with the nominal 90% levels, whereas the estimates corresponding to the two Gaussian shocks were inconsistent, thereby deviating from their asymptotic distributions.

The Case of Weak Identification

This section examines the performance of the bootstrap diagnostic under the conditions of weak identification. We define weak identification as a sequence of data-generating processes under ICA 1 for a sample size, $T$, in which one of the non-Gaussian structural shocks, say $\varepsilon_{2,t}$, follows a NIG distribution with excess kurtosis $\kappa_T \rightarrow 0$ as $T \rightarrow \infty$, while remaining structurally non-Gaussian for any finite $T$. Specifically, we parameterize the approach to Gaussianity through the NIG parameters $\alpha$ and $\delta$: as $\alpha \rightarrow \infty$ and $\delta/\alpha = \sigma^2$, the NIG distribution converges to $N(0, \sigma^2)$. In our simulation, we calibrate $\kappa_{2,T} \sim 0.75$ at $T = 300$, which represents a regime where the information content of the fourth-order cumulant for identification is substantially reduced relative to the strongly identified case ($\kappa_{2,T} \sim 2$). The analysis is simulation-based, consistent with the approach of moneta_identification_2022.

Under the conventional strategy of verification through univariate normality tests, subsequent inference regarding the estimates of the structural impact matrix, $B_0$ is conditional upon the rejection of these normality tests. In finite samples, this pre-testing bias is particularly acute under weak identification, where estimated residuals more closely approximate Gaussian distributions. This implies an increase in the probability of failing to reject the null hypothesis of normality, as documented in Table (ref).

To assess the performance of the bootstrap diagnostic, we compute and compare conditional and unconditional probability coverages for estimates of the impact matrix for a realistic finite sample size of $T = 300$ across $N_S = 500$ Monte Carlo simulations. The conditional coverage based on residual non-normality is calculated by conditioning on the rejection of the univariate Jarque-Bera normality test at the 1% significance level for all estimated reduced-form innovations' series\footnote{Nearly four percent of the total Monte Carlo replications, $N_S = 500$, reject the univariate Jarque-Bera test of normality, at 1% nominal significance level, for all three VAR residuals' series, $\hat{u}_{i,t}, \, i = 1,2,3$.}, $\hat{u}_{i,t}, \ i = 1,2,3$. The conditioning on bootstrap diagnostic procedure operates as follows. For each Monte Carlo replication, we randomly choose $S = 1000$ bootstrap sequences of the on-impact matrix estimates ${\{\hat{B}_{T,1}^*, \hat{B}_{T,2}^*, \cdots, \hat{B}_{T,M}^*\}}_{s}, s = 1, 2, \cdots, S$ of length $M = (1/2)T^{3/5}$. We apply the multivariate Doornik-Hansen normality test at nominal 1% significance level to each bootstrap sequence and compute the rejection rate across the $S$ bootstrap sequences for each simulation. Conditional coverage (bootstrap diagnostic) is then calculated over only those simulations in which the bootstrap-based rejection rate is equal or lower than the overall rejection rate, computed across all $N_S$ Monte Carlo simulations. Coverage probabilities are reported for: nominal percentile confidence intervals (constructed from the empirical distribution of the point estimates across Monte Carlo simulations), asymptotic confidence intervals (based on the asymptotic normal approximation), studentized bootstrap confidence intervals, and percentile bootstrap confidence intervals. All confidence intervals employ a nominal 90% confidence level. Results are presented in Table (ref) for the strongly identified scenario (excess kurtosis approximately 20 and 2) and Table (ref) for the weakly identified scenario (excess kurtosis approximately 3 and 0.75).

Finite-Sample Performance Under Strong Identification

Under the strongly identified specification (Table (ref)), the unconditional coverage probabilities for all four inference methods are close to their nominal levels. Point estimates exhibit negligible bias relative to true parameter values (column a), with discrepancies not exceeding 0.03 in absolute value. We note a minimal impact of residual-based pre-test conditioning. The conditional coverage probabilities on rejection of univariate normality tests, reported in columns (g)-(j), differ only marginally from their unconditional counterparts in columns (c)-(f). Analogously, the bootstrap diagnostic conditioning criterion yields coverage probabilities in columns (k)-(n) that are virtually indistinguishable from unconditional measures.

Inferential Distortion Under Weak Identification:

We observe substantial undercoverage across parameters which is consistent with the breakdown of non-Gaussian identification as more than one structural shock approaches Gaussianity. The heterogeneity in coverage across parameters reflects the differential identification strength: elements of $B_0$ corresponding to the more strongly non-Gaussian shock (excess kurtosis approximately 3) retain moderate identification and exhibit relatively modest undercoverage, whereas elements associated with the nearly Gaussian shock (excess kurtosis approximately 0.75) suffer a severe lack of identification strength.

The third element provides a stark illustration: unconditional asymptotic coverage stands at 0.84, yet conditioning on residuals' normality rejection reduces this to 0.69. This adverse selection extends to the bootstrap inference as well. The third element displays studentized bootstrap probability coverage declining from 0.83 to 0.75 upon conditioning, while the fourth element exhibits a reduction from 0.88 to 0.75. Paradoxically, the same conditioning criterion that induces severe undercoverage for some parameters generates artificial inflation in probability coverages for other parameters, and hence an illusion of precision. The ninth element displays conditional asymptotic coverage of 0.94 compared to unconditional coverage of 0.84. This spurious elevation of coverage probabilities for poorly identified parameters represents an incoherent form of diagnostic failure: researchers examining these conditional probabilities might erroneously infer that the parameter is well-identified and precisely estimated, when in fact the conditioning mechanism has merely selected those \enquote{fortunate} samples in which estimation happened to succeed despite the fundamental lack of identification strength.

The heterogeneity of coverage distortions across parameters is equally problematic from a practical standpoint. Partial identification, see maxand_identification_2020, guay_identification_2021, involves focusing on a subset of estimates corresponding to a target structural shock, and deriving implications from the inferred impulse response functions. Unlike specifications exhibiting uniform coverage inflation, where the artificial precision might alert careful researchers to selection bias, the mixed pattern signals a common source: the unreliability of residual-based pre-testing under weak identification.

The Bootstrap Diagnostic as a Stable Alternative:

In contrast, conditioning on the bootstrap diagnostic exhibits robustness across the parameter space, inference methods, and identification regimes. The conditional coverage probabilities reported in columns (k)-(n) track their unconditional counterparts in columns (c)-(f) without exhibiting any systematic pattern of inflation or deflation. This coherence suggests that the bootstrap diagnostic conditioning event does not interact arbitrarily with an inference method but preserves the relative strength of identification, or lack thereof, under unconditional inference.

Overall, the bootstrap-diagnostic conditioning evaluates the stability of the identified structure rather than the distributional features of the reduced-form innovations. Therefore, it avoids both the adverse selection that generates undercoverage and the favorable selection that generates spurious precision under residual-based conditioning, and yields probability coverage that reflects the inferential uncertainty inherent in the weakly identified scenarios.

Empirical Illustration

In this section, we consider an empirical illustration to demonstrate the potential of our bootstrap-based approach to verify the identification conditions for non-Gaussian SVARs, and investigate whether uncertainty is an exogenous source of business cycle fluctuations or an endogenous response to the dynamics of macroeconomic drivers.

We consider a trivariate SVAR which includes (i) a measure of macroeconomic uncertainty ($U_{Mt}$), (ii) a measure of real economic activity ($Y_{t}$), and (iii) a measure of financial uncertainty ($U_{Ft}$), adopted from ludvigson_uncertainty_2021 and angelini_uncertainty_2019. To achieve identification, the former imposes a series of inequality restrictions to allow simultaneous feedback between uncertainty and real activity (shock-based restrictions), which involve identified shocks to have certain statistical properties during historically influential events, such as the Black Monday Crash of 1987 and the Global Financial Crisis of 2008. It also includes external variables and their correlations with uncertainty shocks, such as stock market returns and gold prices, to generate additional inequality constraints. On the other hand, the latter extends the approach of heteroscedasticity-based identification by exploiting breaks in the unconditional volatility and merges it with non-recursive zero restrictions.

Instead, we adopt identification based on non-Gaussianity and hence, do not require economic identification restrictions. The monthly data spans from 1960:07 to 2025:06 ($T = 780$). For a detailed exposition on the construction of uncertainty measures, we refer the reader to jurado_measuring_2015 and ludvigson_uncertainty_2021. For the measure of real activity, we consider the log of real industrial production (FRED database). Based on AIC and BIC criteria, we use $p = 4$ lags in the VAR specification. Also, the multivariate Ljung-Box test for serial correlation (up to 12 lags) in the estimated residuals does not show evidence against the null of no serial correlation. For the non-Gaussian identification, we use the NGML estimator with the structural shocks assumed to follow NIG distribution\footnote{This choice is motivated by the empirical evidence suggesting that the identified shocks in ludvigson_uncertainty_2021 display significant skewness and excess kurtosis. Moreover, using the Student-$t$ distribution for the structural shocks we obtain similar results.}. Hence, the estimated reduced-form VAR is specified as:

where $X_{t} = (U_{Mt}, Y_{t}, U_{Ft})'$ is the vector of measures of uncertainty and real activity, $B$ is the impact matrix and $u_{t} = (u_{1,t}, u_{2,t}, u_{3,t})'$ are the reduced-form innovations. We are interested in the dynamic impulse responses of $X_{t+h}$, at horizon $h = 0,1,2,\cdots$, to structural shocks, $(\varepsilon_{Mt}, \varepsilon_{Yt},\varepsilon_{Ft})'$ corresponding to macro uncertainty, real activity and financial uncertainty shocks, respectively. Once identified, we can compute the IRFs to one-standard deviation structural shocks as:

where, $R = [I_3 , 0_{3 \times 9}]$ is a selection matrix such that $R'R = I_3$, $\hat{\mathbf{\Pi}}$ is the companion form of the estimated VAR coefficient matrices, and $\hat{B}_T$ is the estimated mixing matrix from NGML. The instantaneous response i.e., response of $X_{t+h}$, $h = 0$, to one-standard deviation shock in $\varepsilon_t$ is captured by $B$ ($\hat{\Psi}_0 = \hat{B}_T$).

Table (ref) reports the estimates of the impact matrix $B$ (to one-standard deviation in structural shocks) and the excess kurtosis $\kappa$ for the structural shocks. The estimates show significant excess kurtosis in all three structural shocks, supporting the assumption of non-Gaussianity in the structural shocks.

Before interpreting the dynamic responses of the variables to the estimated structural shocks, we assess the validity of the identification conditions through our bootstrap-based approach. Table (ref) reports the relative rejection frequencies of the multivariate Doornik-Hansen test for normality of $S = 1000$ bootstrapped estimate sequences, ${\{\hat{B}_{T,1}^*, \hat{B}_{T,2}^*, \cdots, \hat{B}_{T,M}^*\}}_{s}, s = 1, 2, \cdots, S$, for $T = 780$ and for different values of $M = 10,12,15$. The choice of $M$ is deliberately small relative to $T$. The results show that for all nominal significance levels, the rejection frequencies mostly follow the nominal levels, suggesting no evidence against the null hypothesis of valid identification conditions. This is consistent with the simulation results where under the null of valid identification conditions, the conditional (on the data) distribution of the $p$-values of the bootstrap tests converges to the Uniform distribution in probability, and the relative rejection frequencies are close to their corresponding significance levels\footnote{The rejection frequencies increase with $M$. For $T = 780$, the three choices $M \in \{10, 12, 15\}$ yield $M/T \in \{0.013, 0.015, 0.019\}$, all satisfying the joint divergence requirement. Though the rejection frequencies at $M = 15$ and the 5% and 10% significance levels (0.10 and 0.19, respectively) exceed their nominal levels, it does not imply test degeneracy since they have small-sample validity but suggest that even within the valid $M/T$ range, practitioners should prefer smaller values of $M$ when $T$ is moderate. This is consistent with our simulation evidence, where $M = (1/3)T^{3/5}$ provides better size control than $M = (1/2)T^{3/5}$ across both $T = 300$ and $T = 500$.}. This is consistent with the findings of ludvigson_uncertainty_2021 who document significant heavy tails (excess kurtosis) and skewness in the identified shocks.

IRFs: Although exact identification is achieved through NGML, the statistically identified shocks lack direct economic interpretation. Consequently, we use sign restrictions to label the identified shocks. Consistent with the literature in which macroeconomic and financial uncertainty co-move positively jurado_measuring_2015, carriero_measuring_2018, we assume that a positive shock to either uncertainty index elicits a positive response in both indices. To economically distinguish the two uncertainty shocks, we impose that macroeconomic uncertainty is countercyclical i.e., its response to a positive real activity shock is negative bloom_impact_2009, while leaving the response of financial uncertainty unrestricted. Figure (ref) reports the estimated IRFs of $X_{t}$ to one-standard deviation positive shocks in $\varepsilon_{Mt}$, $\varepsilon_{Yt}$ and $\varepsilon_{Ft}$ respectively, along with the 68% bootstrap CIs (dashed red lines). The bootstrap-based CIs are constructed using residual-based MB bootstrap, with $N = 10,000$ replications and block length\footnote{We choose the block length, $l$ as the largest integer smaller than $5.03T^{1/4}$, in line with the current literature see jentsch_proxy_2016, mertens_dynamic_2019, angelini_identification_2024. The results are robust to varying block lengths between $10$ and $40$.} = $26$.

First, the response of real activity to a financial uncertainty shock is statistically insignificant on impact and turns significantly negative only with a lag. This delayed effect is consistent with financial uncertainty acting as an endogenous response rather than an exogenous impulse, in contrast to ludvigson_uncertainty_2021, who find that an increase in financial uncertainty is an exogenous impulse causing an immediate and persistent decline in real activity. Notably, the immediate response of financial uncertainty to a real activity shock is statistically insignificant too, so the procedure independently recovers the contemporaneous zero restrictions $B_{FY}=B_{YF}=0$ of angelini_uncertainty_2019 without imposing them. Second, a positive real activity shock lowers macroeconomic uncertainty significantly and persistently, consistent with the documented counter-cyclicality of macro uncertainty. Third, real activity declines significantly and persistently following a positive macroeconomic uncertainty shock, although its immediate response is positive. This short-run increase is consistent with ludvigson_uncertainty_2021 and with growth-options theories, in which a mean-preserving spread in risk raises expected profits and can induce firms to invest and hire kraft_growth_2018, segal_good_2015; the subsequent persistent decline reflects the conventional contractionary effect of heightened uncertainty.

Overall, non-Gaussian identification suggests that macroeconomic uncertainty acts as a driver of business-cycle fluctuations whereas financial uncertainty behaves largely as an endogenous response, affecting real activity only indirectly through its positive co-movement with macroeconomic uncertainty.

Conclusion

This paper develops a bootstrap-based approach to assess the identification conditions of non-Gaussian SVARs, in which identification through independent shocks requires that at most one structural shock be Gaussian. Since the reduced-form innovations are linear mixtures of the structural shocks, their non-Gaussianity is uninformative about how many of the underlying shocks are Gaussian. This implies that the conventional approach of pre-testing the reduced-form innovations for normality cannot detect the failure of identification due to extra Gaussian shocks. We instead propose a bootstrap-based approach where under valid identification and maintained regularity conditions, the conditional bootstrap distribution of the NGML impact-matrix estimator is asymptotically normal, and it departs from normality when identification fails. The verification thus reduces to a test of the normality of the bootstrap replications of the impact matrix.

We provide three theoretical results which underpin the diagnostic. First, the diagnostic is valid across the entire null: even in the single-Gaussian configuration that lies on the boundary of the null, where the full-parameter information matrix is singular but the impact-matrix estimator retains a non-singular profile information and remains consistent and asymptotically normal (Proposition (ref)). Second, by letting the number of bootstrap replications and the sample size diverge jointly (at an appropriate rate) rather than sequentially, the test statistic is asymptotically pivotal and ancillary, so that conditioning on the diagnostic does not induce any pre-testing bias. Third, a univariate version of the diagnostic localizes, in large samples, which shocks are responsible for an identification failure, extending its use to partially identified non-Gaussian SVARs.

The Monte Carlo evidence, based on Normal-inverse Gaussian shocks, shows that the diagnostic attains near-nominal size under valid identification and power that increases in both the sample size and the number of bootstrap replications, whereas conventional residual-normality tests cannot discriminate between valid and invalid specifications. Under weak identification with a near-Gaussian shock, conditioning on the diagnostic preserves the probability coverage of the estimates, while conditioning on residual pre-tests substantially distorts it. Finally, in an empirical illustration investigating the relationship between business cycle fluctuations and macro-financial uncertainty measures, the diagnostic finds no evidence against the validity of non-Gaussian identification. The implied impulse responses suggest that macroeconomic uncertainty acts as an exogenous driver of real activity while financial uncertainty behaves largely as an endogenous response.

Several extensions remain open. The finite-sample power of the diagnostic is governed by the size-power trade-off intrinsic to the joint regime of $M$ and $T$, and its behavior against local, near-Gaussian alternatives requires further research. The framework also allows extension beyond the fully independent, likelihood-based setting considered here, for instance, to moment-based estimators that relax complete independence, for which analogous bootstrap diagnostic could be developed to verify their identifying assumptions.