The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
119,407 characters
Two Gaussians, Too Many: A bootstrap-based approach to assess identifiability in non-Gaussian Structural Vector Autoregressions
\setcounter{page}{1}
\title[\tiny Two Gaussians, Too Many]{Two Gaussians, Too Many: A bootstrap-based approach to assess identifiability in non-Gaussian Structural Vector Autoregressions}
\author[]{Paritosh Shankarrao Junare}
\thanks{E-mail: [email removed]. I am immensely grateful to Luca Fanelli, Giuseppe Cavaliere and Giovanni Angelini for their guidance and support throughout this work. I am also thankful to Markku Lanne, Jani Luoto, Mika Meitz and Savi Virolainen for their insightful comments. Useful feedback from Alessandra Luati and other participants of 2026 RCEA Time Series Workshop, Madrid; Katerina Petrova and other participants of the IAAE Annual Conference 2026, Lisbon is gratefully acknowledged.}
\dedicatory{Department of Economics, University of Bologna\\
July 19, 2026}
\begin{abstract}
\noindent Standard pre-tests of normality on reduced-form innovations are insufficient to detect two or more Gaussian shocks and hence, the failure of identification in non-Gaussian SVARs. We instead propose a bootstrap-based approach to evaluate the asymptotic validity of this condition by measuring the divergence between the conditional bootstrap distribution of a maximum likelihood estimator and its limiting distribution under valid identification. We show that, under valid identification and certain regularity conditions, the conditional bootstrap distribution of the impact matrix is asymptotically normal, so the diagnostic reduces to a test of normality of the bootstrap replications. The diagnostic remains valid in the single-Gaussian case, where the shape parameter of the Gaussian shock lies on the boundary, and the full-parameter information is singular; this establishes its validity across the entire null. Under the null of valid identification, the diagnostic induces no pre-testing bias as bootstrap replications and sample size diverge jointly at an appropriate rate. The joint divergence ensures that the test statistic, conditional on the data, is asymptotically pivotal, so conditioning on the diagnostic does not distort subsequent inference. Monte Carlo simulations with Normal-Inverse Gaussian shocks show that the diagnostic attains near-exact nominal size under valid identification and detects the failure due to multiple Gaussian shocks with power increasing in the sample size. Under weak identification with a near-Gaussian shock, conditioning on the bootstrap diagnostic, unlike on residual-based normality pre-tests, preserves the probability coverage of the estimates. Based on estimates of a SVAR model in the macroeconomic and financial uncertainty literature, we demonstrate its potential as a practical, robust tool for validating non-Gaussian identification without pre-testing bias.
\end{abstract} \maketitle
\section{Introduction}{\label{sec:intro}}
\noindent 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 \citep{blanchard_empirical_2002}, cointegration-based constraints
\citep{blanchard_dynamic_1989}, sign-restrictions \citep{uhlig_what_2005}, or employ external information as instruments for exogenous sources of variations \citep{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 \citet{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 \citet{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 \citep{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 \citet{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 \textit{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: \citet{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, \citet{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: \emph{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 \citet{lanne_identification_2017}, is asymptotically standard
normal, so the check reduces to a normality
test on the bootstrap replications. \emph{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.} \emph{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
\citep{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
\citep{guay_identification_2021} or a specification test on the estimated shocks \citep{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
\citep{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) \textbf{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) \textbf{ICA 1}, where one shock is Gaussian and the other two are non-Gaussian, and (iii) \textbf{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 \textbf{ICA 0} and \textbf{ICA 1}, while under the identification failure of \textbf{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 \citet{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 \citet{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 \citet{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{sec:ICA} 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{sec:BSdiag} introduces the bootstrap diagnostic for
assessing identification validity, Section \ref{sec:MC} presents Monte Carlo results, Section \ref{sec:WKId} discusses the case of weak identification where apart from one existing Gaussian shock, there is another \textit{nearly Gaussian} shock, Section \ref{sec:EMP} details an empirical illustration and Section \ref{sec:concl} concludes. Appendix contains notational preliminaries, proofs of the main results and additional Monte Carlo evidence.
\section{Identification in non-Gaussian SVARs}\label{sec:ICA}
\noindent This section elaborates identification in non-Gaussian SVARs through the NGML estimator of \citet{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{prop:onegaussian}). Furthermore, under the null of valid identification, we also provide the limiting distribution of the bootstrap estimator of the impact matrix (Proposition \ref{prop:BSvalidity}).\\
Consider a VAR($p$) model in its companion form:
\begin{equation}
\mathbf{Y}_t = \mathbf{\Pi} \mathbf{Y}_{t-1} + \mathbf{U}_t
\label{eq:companionVAR}
\end{equation}
where,
\begin{equation}
\mathbf{Y}_t = \underset{np \times 1}{\begin{pmatrix}
y_t \\
y_{t-1} \\
\vdots \\
y_{t-p+1}
\end{pmatrix}}, \quad
\mathbf{\Pi} = \underset{np \times np}{\begin{pmatrix}
\Pi_1 & \Pi_2 & \cdots & \Pi_p \\
I_n & 0 & \cdots & 0 \\
0 & I_n & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & I_n
\end{pmatrix}}, \quad
\mathbf{U}_t = \underset{np \times 1}{\begin{pmatrix}
u_t \\
0 \\
0 \\
\vdots \\
0
\end{pmatrix}}
\label{eq:setup_companion_matrices}
\end{equation}
and $u_t = B\varepsilon_t$ with,
\begin{equation}
\underset{n \times n}{B} = \begin{pmatrix}
\beta_{11} & \beta_{12} & \cdots & \beta_{1n} \\
\beta_{21} & \beta_{22} & \cdots & \beta_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
\beta_{n1} & \beta_{n2} & \cdots & \beta_{nn}
\end{pmatrix}
\label{eq:setup_Bmatrix}
\end{equation}
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 \eqref{eq:companionVAR} is assumed to be stable i.e., the companion matrix $\mathbf{\Pi}$ satisfies:
\begin{assumption}$\det(I_{np} - \mathbf{\Pi}z) \neq 0, \ \forall z \in \mathbb{C}, |z| \leq 1$.
\label{as:stability}
\end{assumption}
Beyond a stable VAR, we assume that the structural shocks are mutually
\textit{independent} and \textit{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 \citet{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, \citet{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 \citep{kilian_structural_2017}}. However, recent literature shows that it is not necessary to assume complete independence of the shocks
to achieve identification. \citet{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 \citet{guay_identification_2021}, \citet{lanne_identifying_2023}, and/or cumulant tensors
\citet{mesters_non-independent_2024}.
\subsection{NGML estimation:}\label{subsec:NGMLest}
Apart from non-parametric algorithms, see \citet{hyvarinen_independent_2000}, \citet{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{\citet{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.
\begin{assumption}
\textit{The structural shocks, $\varepsilon_{i,t}, i = 1,2,\cdots n$, are mean-zero, mutually independent and identically distributed (i.i.d.) processes,
with diagonal covariance matrix i.e., $\mathbb{E}(\varepsilon_t\varepsilon_t') = \operatorname{diag}(\sigma_i^2), \ \sigma_i^2 \in \mathbb{R}^{+}/\{0\}$.}
\label{as:indep}
\end{assumption}
\begin{assumption}
\textit{At most one of the components of $\varepsilon_t$ is Gaussian.}
\label{as:gau}
\end{assumption}
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 \citet{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 \citep{schrodinger1915theorie}, and \citet{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{\citet{jarocinski_estimating_2024},
\citet{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 \citep{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:
\begin{equation}
f(x; \alpha, \gamma, \delta, \mu) = \frac{\alpha\delta K_1(\alpha\sqrt{\delta^2 + (x - \mu)^2})}{\pi\sqrt{\delta^2 + (x - \mu)^2}} \exp(\delta\sqrt{\alpha^2 - \gamma^2} + \gamma(x - \mu))
\label{eq:NIGpdf}
\end{equation}
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 \citet{lanne_identification_2017} or the pseudo-ML of \citet{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 \emph{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{\citet{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:
\begin{equation}
\mathcal{L}_T(\lambda) = T^{-1}\sum_{t=1}^{T}\ell_t(\lambda),
\label{eq:loglik}
\end{equation}
where,
\begin{equation}
\ell_t(\lambda) = \sum_{j=1}^{n} \log f_j(\sigma_j^{-1}\imath_j'B(\beta)^{-1}\hat{u}_t; \theta_j) - \log \det(B(\beta)) -
\sum_{j=1}^{n} \log \sigma_j
\label{eq:lik}
\end{equation}
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:
\begin{equation}
\hat{\lambda}_T = (\hat{\Pi}_T, \hat{\beta}_T, \hat{\sigma}_T, \hat{\theta}_T) = \underset{\beta, \sigma, \theta}{\arg\max} \
\mathcal{L}_T(\lambda)
\label{eq:NGMLestimator}
\end{equation}
The estimated structural shocks are obtained by multiplying a transposed $j^{th}$ unit vector, $\imath_j'$,
\begin{equation}
\hat{\varepsilon}_{j,t} = \hat{\sigma}_j^{-1}\imath_j'B(\hat{\beta}_T)^{-1}\hat{u}_t, \ j = 1,2,\cdots,n.
\label{eq:est_struc_shock}
\end{equation}
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:
\begin{equation}
\underset{n \times n}{\hat{\Psi}_i} = R \ \hat{\mathbf{\Pi}}^i \ R' \ B(\hat{\beta}_T)\operatorname{diag}(\hat{\sigma}_T), \quad i = 0,1,2,\cdots,h,
\label{eq:IRF}
\end{equation}
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:
\begin{equation}
\hat{B}_T \coloneqq B(\hat{\beta}_T)\underset{n \times n}{\operatorname{diag}}(\hat{\sigma}_T)
\label{eq:BT}
\end{equation}
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{subsec:regcond}. 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.
\begin{remark}[Identification Failure and the Fisher Information Matrix]\label{rem:identification_failure}
Assumptions \ref{as:indep} and \ref{as:gau} are jointly necessary for the positive definiteness of the Fisher information matrix $\mathcal{I}(\lambda_0)$ stated in regularity condition \ref{subsec:regcond} (7). To see this, suppose Assumption \ref{as:gau} fails and two shocks, say $\varepsilon_{i,t}$ and $\varepsilon_{j,t}$ are Gaussian. Then any rotation of the structural space within the Gaussian subspace $\{i,j\}$ leaves the joint distribution of $u_t$ unchanged \citep{gourieroux_statistical_2017}. This implies that the log-likelihood is invariant along the rotation of the Gaussian subspace. Hence, $\mathcal{I}(\lambda_0)$ is singular in that
rotational direction and this causes the NGML estimator to be inconsistent and its bootstrap distribution to deviate from normality: the mechanism which our diagnostic exploits.
\end{remark}
\subsection{Consistency and Asymptotic Normality:} Consider the model representation in equation \eqref{eq:companionVAR},
satisfying the regularity conditions (\ref{subsec:regcond}) with Assumptions (\ref{as:stability}) - (\ref{as:gau}), 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
\eqref{eq:loglik} and \eqref{eq:lik}, satisfies\footnote{For detailed proofs, see \citet{lanne_identification_2017}, Theorem 1. When one structural shock is Gaussian, \eqref{eq:asymp} holds
for the structural sub-vector $\psi=(\operatorname{vec}(\Pi),\beta,\sigma)$
with the robust covariance of Proposition \ref{prop:onegaussian} in place of
$\mathcal{I}(\lambda_0)^{-1}$; the impact-matrix estimator $\hat{B}_T$
remains asymptotically normal, which is the object of interest.}:
\begin{equation}
T^{1/2}(\hat{\lambda}_T - \lambda_0) \xrightarrow{d} \mathcal{N}(0, \mathcal{I}(\lambda_0)^{-1}), \quad \text{as} \ T \rightarrow \infty,
\label{eq:asymp}
\end{equation}
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.,
\begin{equation}
- \mathcal{L}_{\lambda\lambda,T}^{-1}(\hat{\lambda}_T) \coloneqq -\left(T^{-1} \sum_{t=1}^{T}\nabla^2 \ell_{\lambda\lambda,t}(\hat{\lambda}_T)\right)^{-1} \rightarrow \quad \mathcal{I}(\lambda_0)^{-1}, \quad \text{as} \ T \rightarrow \infty,
\label{eq:fisher}
\end{equation}
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{prop:onegaussian} 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 \emph{second} shock also
approaches Gaussianity: the rotational near-indeterminacy of the weak-identification
regime studied in Section~\ref{sec:WKId} which lies in the alternative hypothesis.
\begin{lemma}[Information at the single-Gaussian boundary]\label{lem:profileinfo}
Let the regularity conditions of Appendix~\ref{subsec:regcond} and Assumptions
\ref{as:stability}--\ref{as:indep} hold, and suppose exactly one structural shock,
say shock $k$, is Gaussian. Partition $\lambda=(\psi',\theta')'$ with
$\psi=(\operatorname{vec}(\Pi)',\beta',\sigma')'$ the autoregressive--impact
parameters and $\theta=(\theta_1',\dots,\theta_n')'$ the shape parameters, and let
$\mathcal{I}(\lambda_0)=-\mathbb{E}[\nabla^2_{\lambda\lambda}\ell_t(\lambda_0)]$ and
$\mathcal{J}(\lambda_0)=\mathbb{E}[\nabla_\lambda\ell_t(\lambda_0)\nabla_\lambda\ell_t(\lambda_0)']$
be partitioned conformably. Then:
\begin{enumerate}
\item[(i)] the shape block $\mathcal{I}_{\theta\theta}(\lambda_0)$ is singular, with
null space $\mathcal{N}=\ker\mathcal{I}_{\theta\theta}$ equal to the shape directions
along which shock $k$'s Gaussian-limit density is invariant;
\item[(ii)] the cross-information blocks annihilate $\mathcal{N}$:
$\mathcal{I}_{\psi\theta}(\lambda_0)\,\mathcal{N}=\{0\}$ and
$\mathcal{J}_{\psi\theta}(\lambda_0)\,\mathcal{N}=\{0\}$;
\item[(iii)] the Schur-complement information
$\mathcal{I}_\psi=\mathcal{I}_{\psi\psi}-\mathcal{I}_{\psi\theta}\mathcal{I}_{\theta\theta}^{+}\mathcal{I}_{\theta\psi}$
is independent of the choice of generalized inverse $\mathcal{I}_{\theta\theta}^{+}$
and positive definite, with $\mathcal{J}_\psi$ the corresponding score
covariance.
\end{enumerate}
\end{lemma}
\begin{proposition}[]\label{prop:onegaussian}
Let the regularity conditions of Appendix~\ref{subsec:regcond} and Assumptions
\ref{as:stability}--\ref{as:indep} hold, and suppose exactly one structural shock is
Gaussian. Then the NGML estimator of the autoregressive--impact parameters
$\psi=(\operatorname{vec}(\Pi)',\beta',\sigma')'$, the impact-matrix estimator
$\hat{B}_T=B(\hat\beta_T)\operatorname{diag}(\hat\sigma_T)$ is consistent
and asymptotically normal,
\begin{equation}
T^{1/2}\bigl(\hat{B}_T-B_0\bigr)\;\xrightarrow{d}\;\mathcal{N}\bigl(0,\;\mathcal{V}_B\bigr),
\label{eq:Bblock_normal}
\end{equation}
where $\mathcal{V}_B$ is obtained from the robust estimator
$\mathcal{I}_\psi^{-1}\mathcal{J}_\psi\mathcal{I}_\psi^{-1}$, and $\mathcal{I}_\psi$ is positive-definite from Lemma~\ref{lem:profileinfo}. When all
structural shocks are non-Gaussian the information identity
$\mathcal{I}_\psi=\mathcal{J}_\psi$ holds and $\mathcal{V}_B$ reduces to the efficient
covariance \eqref{eq:SigmaB}.
\end{proposition}
\begin{proof}
See Appendix \ref{proof:onegaussian}.
\end{proof}
\begin{remark}\label{rem:onegaussian_bridge}
Because the diagnostic of Section~\ref{sec:BSdiag} enters only through $\hat{B}_T$
and requires only the asymptotic normality \eqref{eq:Bblock_normal} of its bootstrap
distribution, Proposition~\ref{prop:onegaussian} is exactly what makes the procedure
valid under the full null of at most one Gaussian shock. Consistency of the
residual-based moving block bootstrap (MBB) estimator of $\mathcal{V}_B$ is established in
Appendix~\ref{subsec:proofBSvalidity}.
\end{remark}
Proposition~\ref{prop:onegaussian} has an implication for practitioners that is
independent of the diagnostic. A practitioner who believes one structural shock may be
Gaussian need not decide \emph{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 \eqref{eq:Bblock_normal},
even when a shock lies at the boundary of the parameter space.
By delta method, see \citet{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$,
\begin{align}
\hat{\Sigma}_{{B}_T}(i,j) &= \hat{\Sigma}_{{\sigma}_T}(i,j), \quad \text{for} \ i = j, \\
\hat{\Sigma}_{{B}_T}(i,j) &= \operatorname{Var}(B(\hat{\beta}_T)_{i,j})(\hat{\sigma}_T(j))^2 + \operatorname{Var}(\hat{\sigma}_T(j))(B(\hat{\beta}_T)_{i,j})^2 \nonumber\\
&+2\operatorname{Cov}(B(\hat{\beta}_T)_{i,j}, \hat{\sigma}_T(j))B(\hat{\beta}_T)_{i,j}\hat{\sigma}_T(j),
\quad \text{for} \ i \neq j,
\label{eq:SigmaB}
\end{align}
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{prop:onegaussian}); when all shocks are non-Gaussian, it coincides with the efficient covariance \eqref{eq:SigmaB}.
\subsection{Bootstrap Inference:}\label{sec:BSINF}
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 \citet{bruggemann_inference_2016}, see Appendix section \ref{subsec:RMBB}, 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{subsec:proofBSvalidity}). 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{subsec:RMBB}, and define the studentized bootstrap statistic
\begin{equation}
Q_T^{\ast} \coloneqq \hat{\Sigma}_{B_T}^{-1/2}\,T^{1/2}\big(\hat{B}_T^{\ast}-\hat{B}_T\big),
\end{equation}
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
\begin{equation}
Q_T \coloneqq \hat{\Sigma}_{B_T}^{-1/2}\,T^{1/2}\big(\hat{B}_T-B_0\big),
\end{equation}
which, under valid identification, is asymptotically standard normal, see \eqref{eq:asymp}. $\hat{G}_T^{\ast}$ can be computed with arbitrary precision from the bootstrap replicates
$\{Q_{T,b}^{\ast}\}_{b=1}^{N}$ of Algorithm~\ref{subsec:RMBB} by the empirical
distribution function
\begin{equation}\label{eq:edf}
\hat{G}_{T,N}^{\ast}(x)
\coloneqq \frac{1}{N}\sum_{b=1}^{N}
\mathbf{1}\!\big(\operatorname{vec}(Q_{T,b}^{\ast})\le x\big),
\qquad x\in\mathbb{R}^{n^{2}},
\end{equation}
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{sec:BSdiag}.
\begin{proposition}\label{prop:BSvalidity}
Consider the NGML estimator, $\hat{\lambda}_T = (\hat{\Pi}_T, \hat{\beta}_T, \hat{\sigma}_T, \hat{\theta}_T)$, defined in \eqref{eq:NGMLestimator} where, $\hat{B}_T = B(\hat{\beta}_T)\operatorname{diag}(\hat{\sigma}_T)$, and its bootstrap analog
$\hat{\lambda}_T^{\ast}$ (and $\hat{B}_T^{\ast}$), defined in Algorithm (\ref{subsec:RMBB}). Under the Assumptions \ref{as:stability} - \ref{as:gau}, consistent estimator \eqref{eq:asymp} and
\eqref{eq:fisher} with the regularity conditions (\ref{subsec:regcond}), as $T \rightarrow \infty$:
\begin{equation*}
\hat{\Sigma}^{-1/2}_{\lambda_T}T^{1/2}(\hat{\lambda}_T^* - \hat{\lambda}_T)
\xrightarrow{d^*}_p\ \mathcal{N}(\mathbf{0}_{\dim (\lambda)}, I_{\dim (\lambda)}), \ \text{in probability}
\end{equation*}
where, $\hat{\Sigma}_{\lambda_T}$ is a consistent estimator of the asymptotic covariance
matrix of $\hat{\lambda}_T$.
\end{proposition}
\begin{proof}
See Appendix \ref{subsec:proofBSvalidity}
\end{proof}
\begin{corollary}\label{cor:Bblock}
Under Assumptions \ref{as:stability}--\ref{as:gau} (in particular, at most one
Gaussian structural shock) and the regularity conditions of \ref{subsec:regcond}
holding for the structural sub-vector $\psi=(\operatorname{vec}(\Pi),\beta,\sigma)$ (Proposition \ref{prop:onegaussian}), the bootstrap impact-matrix estimator satisfies, as
$T\to\infty$,
\begin{equation*}
\hat{\Sigma}^{-1/2}_{B_T}\,T^{1/2}\bigl(\hat{B}_T^{\ast}-\hat{B}_T\bigr)
\;\xrightarrow{d^*}_p\;\mathcal{N}\bigl(\mathbf{0}_{n^2},\,I_{n^2}\bigr),
\qquad\text{in probability},
\end{equation*}
where $\hat{\Sigma}_{B_T}$ is a consistent estimator of the robust covariance $\mathcal{V}_B$.
\end{corollary}
Proposition \ref{prop:BSvalidity} 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{prop:onegaussian}), and Corollary \ref{cor:Bblock} 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{appendix:a}.
\section{Bootstrap Diagnostic}\label{sec:BSdiag}
\noindent This section develops the bootstrap diagnostic for the identification conditions of a
non-Gaussian SVAR. Proposition \ref{prop:BSvalidity} and Corollary~\ref{cor:Bblock} established that, under valid
identification, the studentized statistic $Q_T^{\ast}$ is asymptotically standard
normal; equivalently, its bootstrap distribution satisfies
\begin{equation}\label{eq:BSconsist}
\sup_{x\in\mathbb{R}^{n^2}}\bigl|\hat{G}_T^{\ast}(x)-\Phi(x)\bigr|
\xrightarrow{p^{\ast}}_{p}0,\qquad T\to\infty,
\end{equation}
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{rem:identification_failure}). The replications then follow the non-degenerate
distribution induced by the unidentified rotation, so they do not converge to normality
(Lemma~\ref{lemma:altlaw}). 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
\begin{equation}\label{eq:teststat}
d^{*}_{T,M}(x)\coloneqq M^{1/2}\,\hat{\Omega}_T^{-1/2}(x)
\bigl(\hat{G}^{*}_{T,M}(x)-\Phi(x)\bigr),
\end{equation}
\begin{align}
d^{*}_{T,M}(x)
&= M^{1/2}\,\hat{\Omega}^{-1/2}_{T}(x)
\Bigl(\hat{G}^{*}_{T,M}(x) - \hat{G}^{*}_{T}(x)\Bigr)
\notag \\
&+ M^{1/2}\,\hat{\Omega}^{-1/2}_{T}(x)
\Bigl(\hat{G}^{*}_{T}(x) - \Phi(x)\Bigr).
\label{eq:decomp2}
\end{align}
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{cond:BErate}), and this discrepancy is a
function of the observed data. The statistic decomposes into a \emph{resampling}
term, $M^{1/2}\hat{\Omega}_T^{-1/2}(\hat{G}_{T,M}^{\ast}-\hat{G}_T^{\ast})$, and a
\emph{centering} term, $M^{1/2}\hat{\Omega}_T^{-1/2}(\hat{G}_T^{\ast}-\Phi)$ in \eqref{eq:decomp2}, 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$ \emph{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 \citep{roth_pretest_2022}, which can be especially difficult to control
because the impulse responses of interest are non-linear functions of the
estimator.
Following \citet{cavaliere_bootstrap_2025}, we therefore let $M$ and $T$ diverge
\emph{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{rem:pivotal}).
\begin{proposition}\label{prop:bootstraptest}
Let Assumptions \ref{as:stability}--\ref{as:gau} and the regularity conditions of
\ref{subsec:regcond} hold, and let $d^{*}_{T,M}(x)$ be given by \eqref{eq:teststat}.
As $T,M\to\infty$ jointly with $MT^{-2\rho}=o_p(1)$ for some $\rho>0$, for each
$x\in\mathbb{R}^{n^2}$:
\begin{enumerate}
\item[\textbf{(i)}] Under the null $\mathcal{H}_0$ of valid identification i.e., at most one structural
shock is Gaussian,
\begin{equation}\label{eq:testdist}
d^{*}_{T,M}(x)\xrightarrow{d^{*}}_{p}\mathcal{N}(0,1);
\end{equation}
\item[\textbf{(ii)}] Under $\mathcal{H}_1$ where there are two or more
Gaussian shocks,
\begin{equation}
\label{eq:altdist}
\bigl|d^{*}_{T,M}(x)\bigr| \;\xrightarrow{p}\; \infty,
\end{equation}
\end{enumerate}
\end{proposition}
\begin{proof}
See Appendix \ref{subsec:proofbootstraptest}.
\end{proof}
\begin{corollary}[Valid size-$\alpha$ diagnostic]\label{cor:sizetest}
Let $\mathcal{T}_{T,M}$ be any consistent test statistic for the normality of the
bootstrap replications, computed on the raw $\{\hat{B}^{\ast}_{T,b}\}_{b=1}^{M}$ or
equivalently the standardized $\{Q^{\ast}_{T,b}\}_{b=1}^{M}$ and let $c_{1-\alpha}$
be the $(1-\alpha)$-quantile of its standard ($\chi^{2}$) null distribution. Under
the conditions of Proposition~\ref{prop:bootstraptest}(i), the rule
$\phi_{T,M}\coloneqq\mathbf{1}(\mathcal{T}_{T,M}>c_{1-\alpha})$ has asymptotic size
$\alpha$ and, by Proposition~\ref{prop:bootstraptest}(ii), power tending to one
against all fixed alternatives in $\mathcal{H}_1$.
\end{corollary}
\begin{proof}
See Appendix \ref{subsec:proofcorroboottest}.
\end{proof}
\begin{remark}[Asymptotic pivotality and ancillarity]\label{rem:pivotal}
The limit $\mathcal{N}(0,1)$ obtained above is free of the data-generating
parameters, so $d^{\ast}_{T,M}$ is {asymptotically pivotal}, and the test has
asymptotically correct size. The surviving bootstrap-resampling term,
$M^{1/2}\hat{\Omega}_T^{-1/2}(\hat{G}^{\ast}_{T,M}-\hat{G}^{\ast}_T)$, is,
conditional on the data, independent of the original-sample estimator, so the
diagnostic is asymptotically ancillary for the inference target. This implies that, unlike a
conventional residual-based normality pre-test, this diagnostic does not distort subsequent inference \citep[Theorem~3.1]{cavaliere_bootstrap_2025}. Section~\ref{sec:WKId}
confirms by simulation evidence that conditioning on the diagnostic leaves the probability coverage of
the estimates undistorted.
\end{remark}
\begin{remark}[The rate $\rho$ and the choice of $M$]\label{rem:Mchoice}
The choice of $M$ relative to $T$ is governed by the rate $\rho$ in the
joint requirement $MT^{-2\rho}=o_p(1)$. When the bootstrap admits an Edgeworth
expansion, $\rho=\tfrac12$ (Condition~\ref{cond:BErate}) and the requirement reduces
to $M/T\to0$. If $M$ is too large relative to $T^{2\rho}$, the centering term is
not negligible anymore, and $d^{\ast}_{T,M}$ fails to converge to $\mathcal{N}(0,1)$ even
when the bootstrap is consistent, inflating false rejections
\citep{angelini_bootstrap_2022}. $M/T$ should therefore be kept small in finite
samples. In the Monte Carlo study, we choose $M=(\frac{1}{3}) T^{3/5}$ and
$M=(\frac{1}{2}) T^{3/5}$, which balance power against size control across sample sizes.
\end{remark}
\begin{remark}[Implementation: a normality test on the replications]\label{rem:DH}
Since, by Corollary~\ref{cor:Bblock} and Proposition~\ref{prop:bootstraptest}, the
standardized replications $\{Q^{\ast}_{T,b}\}_{b=1}^{M}$ are, conditional on the
data, asymptotically i.i.d.\ standard normal under $\mathcal{H}_0$, any consistent normality test, multivariate,
or univariate applied element-wise, inherits the conclusions of the proposition:
its standard ($\chi^{2}$) null distribution controls asymptotic size, while the
statistic diverges under $\mathcal{H}_1$. Because the Doornik--Hansen and
Jarque--Bera statistics are invariant to non-singular affine maps (respectively to
location and scale), they take the same value on the raw replications
$\{\hat{B}^{\ast}_{T,b}\}$ and on the standardized $\{Q^{\ast}_{T,b}\}$. We apply
them to the raw $\{\hat{B}^{\ast}_{T,b}\}$, the standardized version being equivalent
and the size guarantee of Corollary~\ref{cor:sizetest} transferring accordingly. We
use the multivariate omnibus\footnote{See
\citet{fresoli_bootstrap_2022,angelini_bootstrap_2022} for the Doornik--Hansen test
in this context; the univariate test is applied to each sequence
$\{\hat{B}^{\ast}_{i,j,T,b}\}_{b=1}^{M}$, $i,j=1,\dots,n$.} test of \citet{doornik_omnibus_2008} and the univariate
test of \citet{jarque_test_1987}, both are moment-based
normality tests with $\chi^{2}$ null distributions and small-sample corrections.
\end{remark}
\section{Monte Carlo Study}\label{sec:MC}
\noindent \subsection{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 \eqref{eq:companionVAR}, with $n=3$ and $p=1$ i.e., a trivariate VAR(1) model:
\begin{equation}
Y_t = \Pi_1 Y_{t-1} + u_t, \quad {t = 1,2,\cdots,T}, \quad u_t = B \varepsilon_t,
\end{equation}
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,
\begin{equation}
\Pi_1 = \begin{pmatrix}
0.30 & 0.60 & -0.40 \\
-0.50 & 0.20 & 0.10 \\
0.70 & 0.10 & -0.30
\end{pmatrix}, \quad
B = \begin{pmatrix}
1.00 & 0.20 & 0.00 \\
0.63 & 1.00 & -0.60 \\
-0.51 & 0.55 & 1.00
\end{pmatrix}, \quad
\label{eq:model_param}
\end{equation}
As noted in equation \eqref{eq:setup_Bmatrix}, 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:
\begin{equation}
\hat{B}_T = B(\hat{\beta}_T) \operatorname{diag}(\hat{\sigma}_T),
\end{equation}
where $\hat{\sigma}_T$ is the estimate of the standard deviations of structural shocks.
\subsection{Distributional Specifications:}\label{sec:dist_spec}
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.}:
\begin{itemize}
\item \textbf{ICA 0:} $\varepsilon_{1,t} \sim NIG(0.3,0.5), \varepsilon_{2,t} \sim NIG(0.5,1)$, and
$\varepsilon_{3,t} \sim NIG(0.7,1.5)$.
\item \textbf{ICA 1:} $\varepsilon_{1,t} \sim NIG(0.3,0.5), \varepsilon_{2,t} \sim NIG(0.5,1)$, and
$\varepsilon_{3,t} \sim \mathcal{N}(0,1)$.
\item \textbf{No ICA:} $\varepsilon_{1,t} \sim NIG(0.3,0.5), \varepsilon_{2,t} \sim \mathcal{N}(0,1)$, and $\varepsilon_{3,t} \sim \mathcal{N}(0,1)$.
\label{dis spec}
\end{itemize}
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 \eqref{eq:NIGpdf}.
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$ \citep{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:
\begin{itemize}
\item \textbf{ICA 0:} $\sigma_1 = 1.29, \sigma_2 = 1.41, \sigma_3 = 1.46$.
\item \textbf{ICA 1:} $\sigma_1 = 1.29, \sigma_2 = 1.41, \sigma_3 = 1.00$.
\item \textbf{No ICA:} $\sigma_1 = 1.29, \sigma_2 = 1.00, \sigma_3 = 1.00$.
\end{itemize}
From the Assumptions \ref{as:indep} and \ref{as:gau}, 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.
\subsection{Insufficient verification with reduced-form innovations:} \label{subsec:insuff_test} 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
\citet{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{rem:pivotal}), so conditioning on it does not
distort subsequent inference.} reduced-form innovations $\hat{u}_t$ for non-Gaussianity.
Indeed, many studies, including \citet{lanne_identification_2017}, \citet{andrade_higher-order_2025}, \citet{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 \textit{at most} one Gaussian shock, the presence of non-Gaussianity of reduced-form innovations
does not preclude the presence of \textit{more than} one Gaussian shocks. The following table illustrates this point.
\begin{table}[H]
\centering
\scalebox{0.9}{
\begin{tabular}{|c|ccc|ccc|}
\hline
& \multicolumn{6}{c|}{\textbf{Specification}} \\
& \multicolumn{3}{c|}{\textbf{ICA 1}} & \multicolumn{3}{c|}{\textbf{No ICA}} \\
\textbf{Sample Size} & \textbf{T = 100} & \textbf{T = 300} & \textbf{T = 500} & \textbf{T = 100} & \textbf{T = 300} & \textbf{T = 500}\\
\hline
\textit{DH} & & & & & & \\
\rule{0pt}{2.5ex}
$\hat{u}_{1,t}$ & 0.99 & 1.00 & 1.00 & 0.99 & 1.00 & 1.00\\
$\hat{u}_{2,t}$ & 0.76 & 0.99 & 1.00 & 0.38 & 0.68 & 0.88\\
$\hat{u}_{3,t}$ & 0.41 & 0.74 & 0.84 & 0.26 & 0.55 & 0.70\\
\hline
\textit{JB} & & & & & & \\
\rule{0pt}{2.5ex}
$\hat{u}_{1,t}$ & 0.99 & 1.00 & 1.00 & 0.99 & 1.00 & 1.00\\
$\hat{u}_{2,t}$ & 0.75 & 0.99 & 1.00 & 0.35 & 0.69 & 0.88\\
$\hat{u}_{3,t}$ & 0.41 & 0.73 & 0.85 & 0.26 & 0.57 & 0.70\\
\hline
\end{tabular}
}
\caption{\scriptsize Empirical rejection frequencies of the univariate Doornik-Hansen and Jarque-Bera tests for normality of the reduced-form innovations, $\hat{u}_t$, for different sample sizes, ${T} = 100, 300, 500$, across ${N_S} = 500$ Monte Carlo simulations,
at the nominal 5\% level, for the two distributional specifications, \textbf{ICA 1} and \textbf{No ICA}, described in Section \ref{sec:dist_spec}.}
\label{tab:unireducedformtest}
\end{table}
\begin{table}[H]
\centering
\scalebox{1}{
\begin{tabular}{|c|ccc|}
\hline
& \multicolumn{3}{c|}{\textbf{Specification}} \\
\textbf{Sample Size} & \textbf{ICA 0} & \textbf{ICA 1} & \textbf{No ICA}\\
\hline
\rule{0pt}{2.5ex}
\textbf{100} & 1.00 & 0.99 & 0.96\\
\textbf{300} & 1.00 & 1.00 & 1.00\\
\textbf{500} & 1.00 & 1.00 & 1.00\\
\hline
\end{tabular}
}
\caption{\scriptsize Empirical rejection frequencies of the multivariate Doornik-Hansen
tests for normality of the reduced-form innovations, $\hat{u}_t$, for different sample sizes, \\
${T} = 100, 300, 500$, across ${N_S} = 500$ Monte Carlo simulations,
at the nominal 5\% level, for the three different distributional specifications described in Section \ref{sec:dist_spec}.}
\label{tab:reducedformtest}
\end{table}
Tables \ref{tab:unireducedformtest} and \ref{tab:reducedformtest} report the empirical rejection frequencies of the multivariate and univariate \citet{doornik_omnibus_2008} and \citet{jarque_test_1987} tests for normality
of $\hat{u}_t$.
Table \ref{tab:unireducedformtest} shows the insufficiency of verifying the assumption of at most one Gaussian shock through testing the reduced-form innovations. Both specifications, \textbf{ICA 1} and \textbf{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.
\subsection{Bootstrap Probability Coverages:}
Tables \ref{tab:EstandCovICA1} and \ref{tab:EstandCovnoICA} report the estimates of the true parameter on-impact matrix $B_{0}$, from the two distributional specifications \textbf{ICA 1} and \textbf{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 \textit{extra} Gaussian shock, in this section we show results pertaining only to specifications \textbf{ICA 1} and \textbf{No ICA}. Results for the specification \textbf{ICA 0}, where all three structural shocks are non-Gaussian can be found in Appendix \ref{subsec:resultsICA0}.}.
Table \ref{tab:EstandCovICA1} 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{prop:BSvalidity} under the regularity conditions ensuring asymptotic consistency of the NGML estimator.
In contrast, Table \ref{tab:EstandCovnoICA} 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 \textbf{No ICA} ($\varepsilon_{1,t}$), remain consistent and comparable to those obtained from valid specifications (\textbf{ICA 0} and \textbf{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 \citep{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 \textbf{No ICA} specification through univariate diagnostics since their corresponding estimates' distribution deviate from their asymptotic distributions.
\begin{table}[H]
\centering
\scalebox{0.8}{
\begin{tabular}{|c|c|cc|c|ccc|}
\hline
& & & & & & \multicolumn{2}{c|}{\textbf{Bootstrap Coverage}}\\
& $\mathbf{B_{0}}$ & $\mathbf{\hat{B}_T}$ & $\mathbf{\hat{B}^*_T}$ & $\mathbf{90\%}$\textbf{CI} & \textbf{Asymp. Cov}. & \textbf{Studentized} & \textbf{Percentile} \\
& (a) & (b) & (c) & (d) & (e) & (f) & (g) \\
\hline
\rule{0pt}{1ex}
& 1.29 & 1.25 & 1.22 & [0.86, 1.72] & 0.80 & 0.66 & 0.70 \\
&& (0.26) & (0.26) & &[0.82,1.61] & [0.55,1.42] & [0.83,1.48]\\
& 0.81 & 0.80 & 0.79 & [0.43, 1.24] & 0.89 & 0.85 & 0.89 \\
&& (0.23) & (0.24) & &[0.43,1.14] & [0.29,1.04] & [0.39,1.10]\\
& -0.66 & -0.64 & -0.62 & [-0.99, -0.29] & 0.91 & 0.84 & 0.85 \\
&& (0.21) & (0.22) & &[-0.95,-0.32] & [-0.87,-0.20] & [-0.92,-0.26]\\
& 0.28 & 0.27 & 0.25 & [0.09, 0.44] & 0.85 & 0.89 & 0.90 \\
&& (0.10) & (0.11) & &[0.13,0.42] & [0.07,0.41] & [0.05,0.44]\\
$\mathbf{T = 100}$ & 1.41 & 1.36 & 1.29 & [0.96, 1.78] & 0.89 & 0.77 & 0.76 \\
&& (0.25) & (0.25) & &[1.01,1.72] & [0.88,1.59] & [0.87,1.61]\\
& 0.78 & 0.72 & 0.66 & [0.25, 1.09] & 0.88 & 0.87 & 0.87 \\
&& (0.24) & (0.26) & &[0.43,1.10] & [0.31,1.04] & [0.26,1.06]\\
& 0.00 & 0.01 & 0.01 & [-0.16, 0.19] & 0.88 & 0.92 & 0.95 \\
&& (0.11) & (0.12) & &[-0.15,0.16] & [-0.18,0.19] & [-0.19,0.22]\\
& -0.60 & -0.54 & -0.48 & [-0.88, -0.11] & 0.87 & 0.90 & 0.93 \\
&& (0.25) & (0.26) & &[-0.89,-0.27] & [-0.92,-0.12] & [-0.91,-0.05]\\
& 1.00 & 0.95 & 0.89 & [0.72, 1.17] & 0.86 & 0.86 & 0.85 \\
&& (0.15) & (0.15) & &[0.75,1.18] & [0.67,1.15] & [0.61,1.13]\\
\hline
& 1.29 & 1.28 & 1.27 & [1.09, 1.49] & 0.87 & 0.81 & 0.84 \\
&& (0.12) & (0.12) & &[1.08,1.47] & [1.01,1.42] & [1.09,1.45]\\
& 0.81 & 0.81 & 0.80 & [0.64, 0.97] & 0.90 & 0.86 & 0.87 \\
&& (0.10) & (0.10) & &[0.64,0.96] & [0.61,0.93] & [0.64,0.95]\\
& -0.66 & -0.66 & -0.65 & [-0.79, -0.53] & 0.92 & 0.89 & 0.89 \\
&& (0.09) & (0.09) & &[-0.80,-0.51] & [-0.77,-0.49] & [-0.78,-0.51]\\
& 0.28 & 0.28 & 0.28 & [0.23, 0.34] & 0.91 & 0.92 & 0.91 \\
&& (0.04) & (0.04) & &[0.22,0.34] & [0.21,0.34] & [0.21,0.34]\\
$\mathbf{T = 500}$ & 1.41 & 1.41 & 1.40 & [1.26, 1.57] & 0.90 & 0.85 & 0.87 \\
&& (0.10) & (0.10) & &[1.25,1.57] & [1.23,1.53] & [1.24,1.55]\\
& 0.78 & 0.77 & 0.77 & [0.63, 0.92] & 0.90 & 0.89 & 0.90 \\
&& (0.09) & (0.09) & &[0.63,0.91] & [0.63,0.90] & [0.63,0.91]\\
& 0.00 & 0.00 & 0.00 & [-0.06, 0.07] & 0.87 & 0.88 & 0.89 \\
&& (0.04) & (0.04) & &[-0.06,0.06] & [-0.06,0.06] & [-0.06,0.07]\\
& -0.60 & -0.59 & -0.59 & [-0.72, -0.46] & 0.91 & 0.88 & 0.90 \\
&& (0.08) & (0.08) & &[-0.72,-0.46] & [-0.72,-0.46] & [-0.72,-0.45]\\
& 1.00 & 0.99 & 0.99 & [0.90, 1.08] & 0.89 & 0.86 & 0.86 \\
&& (0.06) & (0.06) & &[0.90,1.09] & [0.90,1.08] & [0.89,1.08]\\
\hline
\hline
\end{tabular}
}
\caption{\scriptsize Estimates of true parameter on-impact matrix, $B_{0}$, from the process: \\
\textbf{ICA 1:} $\epsilon_{1,t} \sim NIG(0.3,0.5), \epsilon_{2,t} \sim NIG(0.5,1)$, and \\
$\epsilon_{3,t} \sim \mathcal{N}(0,1)$
(see Section \ref{sec:dist_spec}),
with ${T} = 100, 500$, Monte Carlo Simulations, ${N_S} = 500$ and residual-based MB bootstrap replications, ${N} = 999$.\\
(a) True parameter values.\\
(b) NGML estimates, averaged across Monte Carlo simulations.\\
(c) Bootstrap estimates, averaged across Monte Carlo simulations.\\
The parentheses contain their respective standard errors, averaged across Monte Carlo simulations.\\
(d) Empirical 90\% percentile-confidence intervals (CIs). The square brackets contain the median of the upper and lower limits of CIs, averaged across Monte Carlo simulations.\\
(e) Empirical probability coverage, i.e. the frequencies that the asymptotic $90\%$ CIs contain the true parameter values.\\
(f) Bootstrap (studentized) probability coverage of nominal $90\%$ CIs.\\
(g) Bootstrap (percentile) probability coverage of nominal $90\%$ CIs.\\
The square brackets contain the median of the respective CIs, across Monte Carlo simulations.
}
\label{tab:EstandCovICA1}
\end{table}
\begin{table}[H]
\centering
\scalebox{0.8}{
\begin{tabular}{|c|c|cc|c|ccc|}
\hline
& & & & & & \multicolumn{2}{c|}{\textbf{Bootstrap Coverage}}\\
& $\mathbf{B_{0}}$ & $\mathbf{\hat{B}_T}$ & $\mathbf{\hat{B}^*_T}$ & $\mathbf{90\%}$\textbf{CI} & \textbf{Asymp. Cov}. & \textbf{Studentized} & \textbf{Percentile} \\
& (a) & (b) & (c) & (d) & (e) & (f) & (g) \\
\hline
\rule{0pt}{1ex}
& 1.29 & 1.25 & 1.30 & [0.85, 1.73] & 0.79 & 0.61 & 0.65 \\
&& (0.26) & (0.22) & &[0.81,1.59] & [0.52,1.39] & [0.77,1.43]\\
& 0.81 & 0.79 & 0.83 & [0.44, 1.18] & 0.88 & 0.76 & 0.77 \\
&& (0.23) & (0.19) & &[0.42,1.11] & [0.26,1.00] & [0.31,1.01]\\
& -0.66 & -0.65 & -0.66 & [-0.98, -0.29] & 0.91 & 0.80 & 0.79 \\
&& (0.20) & (0.18) & &[-0.96,-0.32] & [-0.85,-0.18] & [-0.87,-0.19]\\
& 0.20 & 0.10 & 0.13 & [-0.13, 0.29] & 0.80 & 0.83 & 0.85 \\
&& (0.18) & (0.15) & &[-0.08,0.33] & [-0.10,0.35] & [-0.07,0.32]\\
$\mathbf{T = 100}$ & 1.00 & 0.71 & 0.92 & [0.07, 1.20] & 0.75 & 0.80 & 0.85 \\
&& (0.35) & (0.27) & &[0.33,1.32] & [0.14,1.24] & [0.11,1.14]\\
& 0.55 & -0.03 & -0.03 & [-1.11, 1.13] & 0.58 & 0.71 & 0.61 \\
&& (1.06) & (0.71) & &[-1.04,0.97] & [-1.08,0.99] & [-0.80,0.72]\\
& 0.00 & 0.07 & 0.07 & [-0.20, 0.30] & 0.77 & 0.85 & 0.95 \\
&& (0.21) & (0.17) & &[-0.15,0.31] & [-0.16,0.30] & [-0.15,0.27]\\
& -0.60 & 0.02 & -0.00 & [-1.15, 1.15] & 0.59 & 0.71 & 0.61 \\
&& (1.09) & (0.69) & &[-1.03,0.99] & [-1.08,1.08] & [-0.81,0.82]\\
& 1.00 & 0.69 & 0.90 & [0.07, 1.16] & 0.74 & 0.82 & 0.85 \\
&& (0.34) & (0.25) & &[0.38,1.25] & [0.17,1.22] & [0.11,1.12]\\
\hline
& 1.29 & 1.28 & 1.36 & [1.10, 1.49] & 0.88 & 0.80 & 0.75 \\
&& (0.12) & (0.10) & &[1.08,1.47] & [1.00,1.41] & [1.05,1.43]\\
& 0.81 & 0.81 & 0.86 & [0.65, 0.96] & 0.91 & 0.84 & 0.80 \\
&& (0.10) & (0.08) & &[0.65,0.97] & [0.61,0.93] & [0.62,0.93]\\
& -0.66 & -0.66 & -0.70 & [-0.80, -0.53] & 0.91 & 0.84 & 0.80 \\
&& (0.09) & (0.08) & &[-0.80,-0.52] & [-0.77,-0.49] & [-0.77,-0.49]\\
& 0.20 & 0.09 & 0.13 & [-0.09, 0.23] & 0.64 & 0.77 & 0.64 \\
&& (0.10) & (0.09) & &[-0.01,0.23] & [-0.01,0.26] & [0.00,0.22]\\
$\mathbf{T = 500}$ & 1.00 & 0.66 & 0.94 & [0.08, 1.17] & 0.62 & 0.76 & 0.87 \\
&& (0.30) & (0.27) & &[0.31,1.25] & [0.16,1.28] & [0.14,1.15]\\
& 0.55 & -0.10 & -0.08 & [-1.14, 1.11] & 0.60 & 0.71 & 0.58 \\
&& (0.85) & (0.82) & &[-1.15,0.94] & [-1.13,1.00] & [-0.89,0.74]\\
& 0.00 & 0.07 & 0.09 & [-0.17, 0.23] & 0.62 & 0.70 & 0.83 \\
&& (0.14) & (0.13) & &[-0.10,0.24] & [-0.10,0.24] & [-0.06,0.20]\\
& -0.60 & 0.08 & 0.04 & [-1.18, 1.16] & 0.62 & 0.74 & 0.59 \\
&& (0.88) & (0.83) & &[-1.02,1.08] & [-1.11,1.11] & [-0.79,0.88]\\
& 1.00 & 0.65 & 0.92 & [0.07, 1.15] & 0.64 & 0.76 & 0.87 \\
&& (0.29) & (0.27) & &[0.30,1.23] & [0.11,1.27] & [0.13,1.13]\\
\hline
\hline
\end{tabular}
}
\caption{\scriptsize Estimates of true parameter on-impact matrix, $B_{0}$, from the process: \\
\textbf{No ICA:} $\epsilon_{1,t} \sim NIG(0.3,0.5), \epsilon_{2,t} \sim \mathcal{N}(0,1)$, and
$\epsilon_{3,t} \sim \mathcal{N}(0,1)$
(see Section \ref{sec:dist_spec}),
with ${T} = 100, 500$, Monte Carlo Simulations, ${N_S} = 500$ and residual-based MB bootstrap replications, ${N} = 999$.\\
(a) True parameter values.\\
(b) NGML estimates, averaged across Monte Carlo simulations.\\
(c) Bootstrap estimates, averaged across Monte Carlo simulations.\\
The parentheses contain their respective standard errors, averaged across Monte Carlo simulations.\\
(d) Empirical 90\% percentile-confidence intervals (CIs). The square brackets contain the median of the upper and lower limits of CIs, averaged across Monte Carlo simulations.\\
(e) Empirical probability coverage, i.e. the frequencies that the asymptotic $90\%$ CIs contain the true parameter values.\\
(f) Bootstrap (studentized) probability coverage of nominal $90\%$ CIs.\\
(g) Bootstrap (percentile) probability coverage of nominal $90\%$ CIs.\\
The square brackets contain the median of the respective CIs, across Monte Carlo simulations.
}
\label{tab:EstandCovnoICA}
\end{table}
\subsection{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 \emph{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)$.
(\citet{hung_behavior_1997}, \citet{tang_general_2021}).
\subsubsection{Multivariate Diagnostic:}
\begin{figure}[h]
\centering
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width = \textwidth]{Images/fan_NIGT500.png}
\subcaption{\textbf{ICA 0}}
\label{fig:500ICA0}
\end{subfigure}
\quad
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width = \textwidth]{Images/fan_NIGOGT500.png}
\subcaption{\textbf{ICA 1}}
\label{fig:500ICA1}
\end{subfigure}
\quad
\begin{subfigure}{0.45\textwidth}
\centering
\includegraphics[width = \textwidth]{Images/fan_NIGWGT500.png}
\subcaption{\textbf{No ICA}}
\label{fig:500NoICA}
\end{subfigure}
\caption{Empirical distribution functions of the $p$-values, $p^*_{M,T,s}(x)$, $s = 1,2,\cdots, \mathbf{S}$ of the multivariate Doornik-Hansen normality test for the estimates of $B_T$ for the three specifications - \textbf{ICA 0}, \textbf{ICA 1} and \textbf{No ICA} - with sample size $T = 500$ and $M = (1/2)T^{3/5}$, for the number of bootstrap sequences, $S = 1000$, across Monte Carlo simulations, $N_S = 500$ The darker and lighter shades indicate the 50 and 90- percentile bands.}
\label{fig:fanchart500}
\end{figure}
Figures \ref{fig:fanchart500} 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{tab:mvnormtest}} for the three specifications (\textbf{ICA 0}, \textbf{ICA 1} and \textbf{No ICA}), with sample size\footnote{The empirical distribution plots for sample size $T=300$ can be found in Appendix \ref{subsec:resultsICA0} (Figure \ref{fig:fanchart300})} $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, (\textbf{ICA 0} and \textbf{ICA 1}), the $p$-values of the multivariate normality test on bootstrap replications are uniformly distributed, whereas under the invalid specification i.e., with \textit{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.
\subsubsection{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{subsec:JBuni}.} to each element of the estimated impact matrix $\hat{B}_T$. Figures \ref{fig:uni_NIGOGT500} and \ref{fig:uni_NIGWGT500} 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 \textbf{ICA 1} and \textbf{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 \textbf{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{tab:EstandCovnoICA}, 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.
\begin{figure}[h]
\centering
\includegraphics[scale = 0.45]{Images/uni_NIGOGT500.png}
\caption{Univariate \textbf{ICA 1}: Empirical distribution function of $p^*_{M,T,S}(x)$ from the univariate Doornik-Hansen test for normality of each element of (bootstrapped) estimated
on-impact matrix $\hat{B}_{T}^*$ for the (valid) specification \textbf{ICA 1} (see Section \ref{sec:dist_spec}), with sample size $\mathbf{T} = 500$, the number of bootstrap replicates in a sequence, $\mathbf{M} = (1/2)T^{3/5}$,
and the number of sequences/tests, ${S} = 1000$, across Monte Carlo Simulations ${N_S} = 500$. The dashed red line indicates the $45^{\circ}$ line. The black line is the median of the empirical distribution of the $p$-values, $\pi^*_{M,T,S}(x)$, while the dark and light shaded areas indicating the $50$ and $90$-percentile intervals, respectively, with the $p$-values on $x$-axis and their corresponding empirical frequencies on the $y$-axis.}
\label{fig:uni_NIGOGT500}
\end{figure}
\begin{figure}[H]
\centering
\includegraphics[scale = 0.45]{Images/uni_NIGWGT500.png}
\caption{Univariate \textbf{No ICA}: Empirical distribution function of $p^*_{M,T,S}(x)$ from the univariate Doornik-Hansen test for normality of each element of (bootstrapped) estimated
on-impact matrix $\hat{B}_{T}^*$ for the (invalid) specification \textbf{No ICA} (see Section \ref{sec:dist_spec}), with sample size $\mathbf{T} = 500$, the number of bootstrap replicates in a sequence, $\mathbf{M} = (1/2)T^{3/5}$,
and the number of sequences/tests, ${S} = 1000$, across Monte Carlo Simulations ${N_S} = 500$. The dashed red line indicates the $45^{\circ}$ line. The black line is the median of the empirical distribution of the $p$-values, $\pi^*_{M,T,S}(x)$, while the dark and light shaded areas indicating the $50$ and $90$-percentile intervals, respectively, with the $p$-values on $x$-axis and their corresponding empirical frequencies on the $y$-axis.}
\label{fig:uni_NIGWGT500}
\end{figure}
\section{The Case of Weak Identification}\label{sec:WKId}
\noindent
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 \textbf{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 \citet{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{tab:reducedformtestcompare}.
\begin{table}[H]
\centering
\scalebox{0.8}{
\begin{tabular}{|c|c|c|}
\hline
& \multicolumn{2}{c|}{\textbf{Specification ICA 1}} \\
\rule{0pt}{2.5ex}
& \textbf{Strongly Identified} & \textbf{Weakly Identified} \\
\hline
\textit{DH} & & \\
\rule{0pt}{2.5ex}
$\hat{u}_{1,t}$ & 1 & 0.88\\
$\hat{u}_{2,t}$ & 0.85 & 0.17\\
$\hat{u}_{3,t}$ & 0.50 & 0.09\\
\hline
\textit{JB} & & \\
\rule{0pt}{2.5ex}
$\hat{u}_{1,t}$ & 1 & 0.87 \\
$\hat{u}_{2,t}$ & 0.84 & 0.16\\
$\hat{u}_{3,t}$ & 0.49 & 0.08\\
\hline
\end{tabular}
}
\caption{\scriptsize Relative rejection frequencies of the univariate Doornik-Hansen and Jarque-Bera tests for normality of the reduced-form innovations, $\hat{u}_t$, for sample size, ${T} = 300$, across ${N_S} = 500$ Monte Carlo simulations,
at the nominal 1\% significance level, for the distributional specification, \textbf{ICA 1} (two non-Gaussian shocks and one Gaussian shock), under the strongly identified and weakly identified scenarios. The two non-Gaussian shocks are independent, NIG distributed with \textbf{excess kurtosis ~ 20 and 2} for the strongly identified scenario, and \textbf{excess kurtosis ~ 3 and 0.75} for the weakly identified scenario.}
\label{tab:reducedformtestcompare}
\end{table}
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{tab:strongidencondcov} for the strongly identified scenario (excess kurtosis approximately 20 and 2) and Table \ref{tab:weakidencondcov} for the weakly identified scenario (excess kurtosis approximately 3 and 0.75).
\subsection{Finite-Sample Performance Under Strong Identification}
Under the strongly identified specification (Table \ref{tab:strongidencondcov}), 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.
\subsection{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 \citet{maxand_identification_2020}, \citet{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.
\subsection{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.
\begin{landscape}
\begin{table}[p]
\vspace*{\fill}
\centering
\scalebox{0.68}{
\begin{tabular}{|cc|cccc|cccc|ccc@{\hspace{2pt}}c|}
\hline
& & \multicolumn{4}{c|}{\textbf{Unconditional Coverage}} & \multicolumn{4}{c|}{\textbf{Conditional Coverage (Residuals' Non-Normality)}} & \multicolumn{4}{c|}{\textbf{Conditional Coverage (Bootstrap Diagnostic)}} \\
$\mathbf{B_{0}}$ & $\mathbf{\hat{B}_T}$ & $\mathbf{90\%}$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} & \textbf{Perc. BootS} & $\mathbf{90\% }$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} &\textbf{Perc. BootS} & $\mathbf{90\% }$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} & \textbf{Perc. BootS} \\
(a) & (b) & (c) & (d) & (e) & (f) & (g) & (h) & (i) & (j) & (k) & (l) & (m) & (n) \\
\hline
1.30 & 1.28 & [1.03, 1.59] & 0.85 & 0.77 & 0.80 & [1.09, 1.60] & 0.92 & 0.86 & 0.90 & [1.04, 1.56] & 0.86 & 0.77 & 0.81 \\
& (0.16) & & (0.49) & (0.52) & (0.45) & & (0.54) & (0.58) & (0.51) & & (0.49) & (0.49) & (0.44) \\
0.82 & 0.81 & [0.62, 1.02] & 0.88 & 0.82 & 0.83 & [0.66, 1.05] & 0.93 & 0.90 & 0.88 & [0.62, 1.01] & 0.89 & 0.82 & 0.83 \\
& (0.12) & & (0.37) & (0.38) & (0.37) & & (0.40) & (0.41) & (0.40) & & (0.37) & (0.37) & (0.36) \\
-0.66 & -0.66 & [-0.85, -0.48] & 0.89 & 0.84 & 0.84 & [-0.87, -0.53] & 0.93 & 0.94 & 0.90 & [-0.84, -0.48] & 0.90 & 0.84 & 0.85 \\
& (0.11) & & (0.36) & (0.34) & (0.33) & & (0.36) & (0.37) & (0.37) & & (0.36) & (0.33) & (0.33) \\
0.16 & 0.16 & [0.08, 0.23] & 0.89 & 0.89 & 0.91 & [0.08, 0.24] & 0.88 & 0.89 & 0.91 & [0.08, 0.23] & 0.86 & 0.88 & 0.90 \\
& (0.05) & & (0.15) & (0.16) & (0.18) & & (0.15) & (0.16) & (0.19) & & (0.15) & (0.17) & (0.19) \\
0.82 & 0.80 & [0.61, 0.95] & 0.91 & 0.92 & 0.92 & [0.60, 0.95] & 0.89 & 0.90 & 0.90 & [0.58, 0.95] & 0.90 & 0.91 & 0.93 \\
& (0.11) & & (0.34) & (0.33) & (0.35) & & (0.35) & (0.34) & (0.36) & & (0.34) & (0.33) & (0.36) \\
0.45 & 0.43 & [0.16, 0.67] & 0.89 & 0.91 & 0.92 & [0.15, 0.68] & 0.87 & 0.90 & 0.91 & [0.15, 0.70] & 0.86 & 0.89 & 0.92 \\
& (0.16) & & (0.50) & (0.60) & (0.66) & & (0.50) & (0.63) & (0.68) & & (0.50) & (0.62) & (0.72) \\
0.00 & 0.00 & [-0.10, 0.09] & 0.88 & 0.90 & 0.93 & [-0.08, 0.09] & 0.88 & 0.89 & 0.93 & [-0.10, 0.10] & 0.85 & 0.90 & 0.93 \\
& (0.05) & & (0.17) & (0.20) & (0.21) & & (0.18) & (0.20) & (0.21) & & (0.18) & (0.20) & (0.22) \\
-0.60 & -0.57 & [-0.78, -0.36] & 0.86 & 0.90 & 0.90 & [-0.78, -0.36] & 0.85 & 0.89 & 0.90 & [-0.79, -0.36] & 0.83 & 0.89 & 0.90 \\
& (0.12) & & (0.37) & (0.47) & (0.53) & & (0.38) & (0.48) & (0.55) & & (0.37) & (0.47) & (0.60) \\
1.00 & 0.97 & [0.82, 1.10] & 0.89 & 0.89 & 0.87 & [0.80, 1.11] & 0.88 & 0.89 & 0.86 & [0.79, 1.10] & 0.85 & 0.87 & 0.86 \\
& (0.08) & & (0.25) & (0.26) & (0.30) & & (0.25) & (0.26) & (0.31) & & (0.25) & (0.25) & (0.33) \\
\hline
\hline
\end{tabular}
}
\caption{\scriptsize \textbf{Strong Identification}: A comparison of asymptotic and bootstrap coverage probabilities (unconditional and conditional) of the NGML estimates of the on-impact matrix, $\mathbf{B_0}$, under the identification scenario \textbf{ICA 1} (two non-Gaussian shocks and one Gaussian shock) in the trivariate SVAR model \ref{eq:model_param}, with sample size $T = 300$, based on ${N = 999}$ residual-based MBB (Section \ref{subsec:RMBB}) replications and ${N_S = 500}$ Monte Carlo simulations. The two non-Gaussian shocks are independent, NIG distributed with \textbf{excess kurtosis ~ 20 and 2}, respectively.\\
The true parameter values, $\mathbf{B_{0}}$ and the NGML estimate, $\mathbf{\hat{B}_T}$ are reported in columns \textit{(a)} and \textit{(b)}. The unconditional coverage columns \textit{(c)-(f)} report the coverage probabilities based on all the DGPs of the simulation (unconditional on any pre-test or specification test).\\
The conditional coverage (residuals' normality) columns \textit{(g)-(j)} report the coverage probabilities conditional on the \textit{rejection} of the univariate Jarque-Bera normality test at 1\% significance level for \textit{all estimated VAR residual series} - $\hat{u}_{i,t}, \ i = 1,2,3$.\\
The bootstrap diagnostic is conducted by multivariate Doornik-Hansen test for normality at 1\% significance level over the number of bootstrap sequences, $S = 1000$, of length $M = (1/2)T^{3/5}$. The conditional coverage (bootstrap diagnostic) columns \textit{(k)-(n)} report the coverage probabilities conditional on the relative rejection rates (across the bootstrap sequences, ${\{\hat{B}_{T,1}^*, \hat{B}_{T,2}^*, \cdots, \hat{B}_{T,M}^*\}}_{s}, s = 1, 2, \cdots, S$) being \textit{equal to or lower than} the overall rejection rates (across the Monte Carlo replications). The rejection rate is computed at 1\% nominal significance level for the multivariate normality test.\\
The coverage probabilities are reported for Percentile CIs (across the replications), Asymptotic CIs, Studentized Bootstrap CIs and Percentile Bootstrap CIs, respectively. All coverage probabilities are reported for nominal 90\% confidence levels.
}
\label{tab:strongidencondcov}
\vspace*{\fill}
\end{table}
\end{landscape}
\newpage
\begin{landscape}
\begin{table}[p]
\vspace*{\fill}
\centering
\scalebox{0.68}{
\begin{tabular}{|cc|cccc|cccc|ccc@{\hspace{2pt}}c|}
\hline
& & \multicolumn{4}{c|}{\textbf{Unconditional Coverage}} & \multicolumn{4}{c|}{\textbf{Conditional Coverage (Residuals' Non-Normality)}} & \multicolumn{4}{c|}{\textbf{Conditional Coverage (Bootstrap Diagnostic)}} \\
$\mathbf{B_{0}}$ & $\mathbf{\hat{B}_T}$ & $\mathbf{90\%}$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} & \textbf{Perc. BootS} & $\mathbf{90\% }$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} &\textbf{Perc. BootS} & $\mathbf{90\% }$\textbf{CI} & \textbf{Asymp.} & \textbf{Stu. BootS} & \textbf{Perc. BootS} \\
(a) & (b) & (c) & (d) & (e) & (f) & (g) & (h) & (i) & (j) & (k) & (l) & (m) & (n) \\
\hline
1.47 & 1.43 & [1.25, 1.60] & 0.84 & 0.74 & 0.75 & [1.24, 1.72] & 0.81 & 0.81 & 0.81 & [1.25, 1.60] & 0.86 & 0.75 & 0.76 \\
& (0.10) & & (0.32) & (0.32) & (0.36) & & (0.38) & (0.45) & (0.43) & & (0.32) & (0.32) & (0.36) \\
0.92 & 0.89 & [0.57, 1.15] & 0.84 & 0.86 & 0.87 & [0.50, 1.27] & 0.75 & 0.88 & 0.94 & [0.57, 1.14] & 0.86 & 0.88 & 0.87 \\
& (0.15) & & (0.48) & (0.55) & (0.57) & & (0.40) & (0.46) & (0.49) & & (0.48) & (0.55) & (0.57) \\
-0.75 & -0.73 & [-1.00, -0.42] & 0.84 & 0.83 & 0.87 & [-1.05, -0.69] & 0.69 & 0.75 & 0.94 & [-1.00, -0.42] & 0.86 & 0.85 & 0.88 \\
& (0.15) & & (0.47) & (0.52) & (0.54) & & (0.39) & (0.50) & (0.53) & & (0.47) & (0.52) & (0.54) \\
0.20 & 0.18 & [-0.08, 0.47] & 0.83 & 0.88 & 0.93 & [-0.05, 0.59] & 0.69 & 0.75 & 0.81 & [-0.08, 0.48] & 0.84 & 0.90 & 0.94 \\
& (0.18) & & (0.51) & (0.61) & (0.64) & & (0.43) & (0.62) & (0.64) & & (0.51) & (0.61) & (0.65) \\
1.00 & 0.93 & [0.40, 1.24] & 0.84 & 0.89 & 0.94 & [0.47, 1.35] & 0.88 & 0.88 & 0.75 & [0.48, 1.25] & 0.86 & 0.91 & 0.94 \\
& (0.22) & & (0.59) & (0.69) & (0.79) & & (0.47) & (0.61) & (0.94) & & (0.59) & (0.69) & (0.78) \\
0.55 & 0.40 & [-0.75, 0.98] & 0.79 & 0.87 & 0.89 & [-0.26, 0.94] & 0.88 & 0.88 & 0.88 & [-0.72, 0.95] & 0.80 & 0.89 & 0.91 \\
& (0.41) & & (0.86) & (1.24) & (1.41) & & (0.80) & (1.03) & (1.58) & & (0.86) & (1.24) & (1.42) \\
0.00 & 0.01 & [-0.31, 0.35] & 0.85 & 0.89 & 0.95 & [-0.22, 0.54] & 0.69 & 0.81 & 0.81 & [-0.31, 0.35] & 0.86 & 0.91 & 0.95 \\
& (0.20) & & (0.56) & (0.62) & (0.64) & & (0.45) & (0.49) & (0.57) & & (0.56) & (0.62) & (0.64) \\
-0.60 & -0.44 & [-1.04, 0.72] & 0.78 & 0.88 & 0.89 & [-1.07, 0.22] & 0.69 & 0.69 & 0.81 & [-1.02, 0.69] & 0.79 & 0.90 & 0.90 \\
& (0.39) & & (0.87) & (1.27) & (1.45) & & (0.84) & (1.28) & (1.53) & & (0.87) & (1.26) & (1.46) \\
1.00 & 0.92 & [0.42, 1.20] & 0.84 & 0.89 & 0.92 & [0.31, 1.10] & 0.94 & 0.94 & 0.94 & [0.51, 1.20] & 0.86 & 0.91 & 0.92 \\
& (0.20) & & (0.53) & (0.62) & (0.71) & & (0.47) & (0.62) & (0.74) & & (0.53) & (0.61) & (0.71) \\
\hline
\hline
\end{tabular}
}
\caption{\scriptsize \textbf{Weak Identification}: A comparison of asymptotic and bootstrap coverage probabilities (unconditional and conditional) of the NGML estimates of the on-impact matrix, $\mathbf{B_0}$, under the identification scenario \textbf{ICA 1} (two non-Gaussian shocks and one Gaussian shock) in the trivariate SVAR model \ref{eq:model_param}, with sample size $T = 300$, based on ${N = 999}$ residual-based MBB (Section \ref{subsec:RMBB}) replications and ${N_S = 500}$ Monte Carlo simulations. The two non-Gaussian shocks are independent, NIG distributed with \textbf{excess kurtosis ~ 3 and 0.75}, respectively.\\
The true parameter values, $\mathbf{B_{0}}$ and the NGML estimate, $\mathbf{\hat{B}_T}$ are reported in columns \textit{(a)} and \textit{(b)}. The unconditional coverage columns \textit{(c)-(f)} report the coverage probabilities based on all the DGPs of the simulation (unconditional on any pre-test or specification test).\\
The conditional coverage (residuals' normality) columns \textit{(g)-(j)} report the coverage probabilities conditional on the \textit{rejection} of the univariate Jarque-Bera normality test at 1\% significance level for \textit{all estimated VAR residual series} - $\hat{u}_{i,t}, \ i = 1,2,3$.\\
The bootstrap diagnostic is conducted by multivariate Doornik-Hansen test for normality at 1\% significance level over the number of bootstrap sequences, $S = 1000$, of length $M = (1/2)T^{3/5}$. The conditional coverage (bootstrap diagnostic) columns \textit{(k)-(n)} report the coverage probabilities conditional on the relative rejection rates (across the bootstrap sequences, ${\{\hat{B}_{T,1}^*, \hat{B}_{T,2}^*, \cdots, \hat{B}_{T,M}^*\}}_{s}, s = 1, 2, \cdots, S$) being \textit{equal to or lower than} the overall rejection rates (across the Monte Carlo replications). The rejection rate is computed at 1\% nominal significance level for the multivariate normality test.\\
The coverage probabilities are reported for Percentile CIs (across the replications), Asymptotic CIs, Studentized Bootstrap CIs and Percentile Bootstrap CIs, respectively. All coverage probabilities are reported for nominal 90\% confidence levels.
}
\label{tab:weakidencondcov}
\vspace*{\fill}
\end{table}
\end{landscape}
\section{Empirical Illustration}\label{sec:EMP}
\noindent 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 \citet{ludvigson_uncertainty_2021} and \citet{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 \citet{jurado_measuring_2015} and \citet{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 \citet{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:
\begin{equation}
X_{t} = \Pi_0 + \Pi_{1}X_{t-1} + \Pi_{2}X_{t-2} + ... + \Pi_{4}X_{t-4} + u_{t}
\label{eq:ludvig_VAR}
\end{equation}
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:
\begin{equation}
\underset{3 \times 3}{\hat{\Psi}_i} = (R \ \hat{\mathbf{\Pi}}^i \ R') \ \hat{B}_T, \quad i = 0,1,2,\cdots,h,
\label{eq:IRF_emp}
\end{equation}
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{tab:ludvig_est} 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.
\begin{table}[H]
\centering
\begin{tabular}{ccc @{\hspace{1.5cm}} lc}
\hline\hline
\multicolumn{3}{c}{$\hat{B}$} & \multicolumn{2}{c}{$\hat{\kappa}$} \\
\hline
& & & Equations & \\
0.0100 & -0.0060 & 0.0027 & $U_{Mt}$ & 5.36 \\
(0.0008) & (0.0010) & (0.0005) & & \\
0.0030 & 0.0073 & -0.0004 & $Y_{t}$ & 5.04 \\
(0.0007) & (0.0005) & (0.0003) & & \\
0.0017 & 0.0008 & 0.0256 & $U_{Ft}$ & 9.27 \\
(0.0009) & (0.0008) & (0.0015) && \\
\hline
\hline
\end{tabular}
\caption{Estimates of the contemporaneous on-impact matrix $B$ (to one-standard deviation in structural shocks) and excess kurtosis $\hat{\kappa}$ for the estimated structural shocks in the trivariate SVAR system $X_{t} = (U_{Mt}, Y_{t}, U_{Ft})'$ identified through ICA with NGML estimation with NIG distributed structural shocks along with standard errors in parentheses. The sample period is from 1960:07 to 2025:06 ($T = 780$) and the VAR lag length is set to $p = 4$.}
\label{tab:ludvig_est}
\end{table}
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{tab:emp_test} 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 \citet{ludvigson_uncertainty_2021} who document significant heavy tails (excess kurtosis) and skewness in the identified shocks.
\begin{table}[H]
\centering
\scalebox{1}{
\begin{tabular}{|c|ccc|}
\hline
\textbf{Significance Levels}& \multicolumn{3}{c|}{\textbf{M}} \\
& $\mathbf{10}$ & $\mathbf{12}$ & $\mathbf{15}$ \\
\hline
0.01 & 0.021 & 0.024 & 0.030 \\
0.02 & 0.034 & 0.036 & 0.055 \\
0.03 & 0.042 & 0.051 & 0.080 \\
0.05 & 0.055 & 0.072 & 0.103 \\
0.10 & 0.110 & 0.125 & 0.192 \\
\hline
\end{tabular}
}
\caption{Relative rejection frequencies of multivariate Doornik-Hansen tests 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 consciously small relative to $T$. The results are for the trivariate SVAR system $X_{t} = (U_{Mt}, Y_{t}, U_{Ft})'$ identified through ICA with NGML estimation with NIG-distributed structural shocks, using monthly data from 1960:07 to 2025:06. The bootstrap estimates are computed using the residual-based moving block bootstrap, with number of replications, $N = 10,000$, and block length = $26$. The nominal significance levels are reported in the first column, with the corresponding relative rejection frequencies reported in the subsequent columns for different choices of $M$.}
\label{tab:emp_test}
\end{table}
\textbf{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 \citep{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 \emph{negative} \citep{bloom_impact_2009}, while leaving the response of financial uncertainty unrestricted. Figure \ref{fig:IRF_NGMLBS} 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 \citet{jentsch_proxy_2016}, \citet{mertens_dynamic_2019}, \citet{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 \emph{endogenous} response rather than an exogenous impulse, in contrast
to \citet{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
\citet{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 \citet{ludvigson_uncertainty_2021} and with
\emph{growth-options} theories, in which a mean-preserving spread in risk raises expected profits and can induce firms
to invest and hire \citep{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.
\begin{figure}[H]
\centering
\includegraphics[scale = 0.60]{Images/IRF_NIGMLBS_extended.png}
\caption{Estimated IRFs of the trivariate SVAR system $X_{t} = (U_{Mt}, Y_{t}, U_{Ft})'$ identified through ICA with NGML estimation with NIG-distributed structural shocks, along with the bootstrap-based confidence intervals (CIs) at $68\%$ level (dashed red lines). The solid lines report the impulse responses of $X_t$ to one-standard deviation positive shock. The bootstrap-based CIs are constructed using residual-based MB bootstrap, with $N = 10,000$ replications. The $x$-axis represents the time horizon (in months) and the $y$-axis represents the magnitude of the impulse responses. The sample period is from 1960:07 to 2025:06 ($T = 780$) and the VAR lag length is set to $p = 4$.}
\label{fig:IRF_NGMLBS}
\end{figure}
\section{Conclusion}\label{sec:concl}
\noindent 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{prop:onegaussian}). 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.
\pagebreak
\bibliographystyle{apalike}
\bibliography{Ref1}
\pagebreak