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.
75,140 characters
Bandwidth-Free Inference for Recursive Nonlinear Impulse Response Functions
\maketitle
\begin{abstract}
Recursive nonlinear impulse responses require an estimated innovation law whenever the impact shock is normalized by innovation ranks and future innovations are integrated out. The closest semiparametric recursive construction in the literature estimates the relevant innovation quantile functions smoothly and discusses a direct empirical-residual implementation without developing its complete first-order inference theory. We tackle this gap in a finite-dimensional nonlinear structural autoregression with unrestricted continuous marginal innovation distributions and a fixed normal-rank shock. Our estimator replaces each innovation quantile function with the empirical quantile of generated structural residuals and iterates the same structural transition. For any fixed collection of responses, we establish a joint \(\sqrt{T}\) asymptotic linear representation with four components: direct transition estimation, the effect of transition estimation on residual order statistics, ordinary innovation-quantile estimation, and the shifted impact quantile. After projection through the recursion, the quantile terms admit a residual-rank-and-spacing representation, yielding feasible inference without innovation-density estimation or quantile smoothing. We then characterize the propagated bias from smoothing, establish validity of a full recursive residual bootstrap, and derive the additional covariance contribution from a finite number of simulated paths, providing bandwidth-free inference for the empirical-residual version of the same normal-rank response used in the smooth recursive construction.
\end{abstract}
\noindent\textbf{Keywords:} recursive nonlinear impulse responses; bandwidth-free inference; empirical quantiles; generated residuals; residual bootstrap; simulation error.
\medskip
\noindent\textbf{JEL classification:} C14; C22; C32.
\section{Introduction}
\label{sec:introduction}
\subsection{The unresolved empirical-residual problem}
\label{sec:introduction-problem}
Impulse response functions summarize how a structural disturbance changes the future path of an economic system. In a nonlinear model, the response is commonly computed by iterating an estimated structural transition along shocked and unshocked paths and averaging their difference over the innovations that arrive after impact. When the structural innovation distributions are unknown, those distributions are part of the response calculation and their quantile functions determine both the impact shock in structural units and the future innovations through which the response propagates.
The closest semiparametric recursive analysis represents structural innovations by marginal quantile functions evaluated at Gaussian ranks and computes the response by iterating the nonlinear transition \citep{GourierouxLee2025NonlinearIRF}. Its direct implementation estimates the innovation quantile functions smoothly and even though it discusses drawing innovations from estimated residuals, it does not develop the corresponding first-order or bootstrap theory. We address this and propose a simple method: replace each unknown innovation quantile function with the empirical quantile of estimated structural residuals and use those order statistics in the recursive response.
The empirical-residual implementation is not covered by a standard parametric delta method. The residuals depend on estimated transition parameters; their order statistics are nonsmooth functions of those parameters; the impact shock evaluates one quantile function at a transformed rank; and each perturbation is carried through a nonlinear recursion. Treating the estimated residuals as observed innovations would omit a first-order channel. Treating the impact shock as fixed in structural units would instead change the population response.
We study a finite-dimensional nonlinear structural autoregression with unrestricted continuous marginal innovation distributions and an identified residual map. A selected innovation percentile is mapped to a standard-normal rank, shifted by a fixed amount, and mapped back through the same innovation quantile function. The shocked and unshocked paths start from the same state and receive the same future innovation ranks. The initial state, response variable, shocked component, shock size, future innovation law, and path coupling are fixed before any estimator is introduced. Accordingly, the empirical-quantile and smoothed-quantile procedures estimate one population response and differ only in how they estimate the unknown quantile functions.
The principal estimator first estimates the transition parameter and recovers the structural residuals. It then supplies the componentwise residual order statistics to the recursive response map. The central technical observation is that a fixed-horizon response does not require an unrestricted first-order approximation of the entire quantile process. After the recursive map is differentiated, quantile perturbations enter through a finite collection of weighted integrated-quantile functionals. These weights measure how innovations at different dates and ranks affect the terminal response. This projection permits an observation-level expansion and, after a change of variables, removes the innovation densities from feasible inference.
\subsection{Contributions}
\label{sec:introduction-contributions}
Our main contribution concerns the empirical-residual implementation of a fixed normal-rank recursive response.
First, we establish a joint \(\sqrt{T}\) asymptotic linear representation for the empirical-quantile estimator over any fixed collection of horizons, initial states, response variables, shocked components, and shock sizes. The representation separates direct transition estimation, the effect of transition estimation on generated residual order statistics, estimation of the ordinary impact and future innovation quantiles, and estimation of the shifted impact quantile. Primitive sufficient conditions are given for a scalar nonlinear location--scale autoregression, while the general structural result is stated under the corresponding componentwise generated-quantile expansion. After the quantile effects are projected through the recursive derivative, their influence contributions can be written using residual ranks and adjacent spacings. Feasible covariance estimation therefore requires neither innovation-density estimation nor a quantile-smoothing bandwidth. We are not aware of a prior result that combines these four first-order channels for this empirical-residual recursive estimator.
Next, we derive a same-target comparison with a smoothed-quantile estimator. The leading difference is the smoothing error weighted by the recursive sensitivity of the reported response to innovation quantiles at different dates and ranks. If the propagated bias is \(o\left(T^{-1/2}\right)\), the empirical and smoothed estimators have the same first-order distribution. A bias of order \(T^{-1/2}\) shifts the limiting distribution, while a larger bias dominates sampling uncertainty. This result isolates the effect of smoothing without changing the shock definition, the future innovation law, or the response being estimated.
Finally, we establish conditional validity of a full recursive residual bootstrap. Each replication regenerates a complete time series, re-estimates the transition, reconstructs the residual quantile functions, and recomputes the response. This re-estimation reproduces all four components of the influence representation. We separately characterize numerical integration error. When the number of simulated paths is proportional to the sample size, it contributes an additional covariance term; when it grows faster, the term is asymptotically negligible. Path-only resampling is shown to estimate numerical integration uncertainty conditional on the fitted model rather than sampling uncertainty in the recursive estimator.
\subsection{Relation to the literature}
\label{sec:introduction-literature}
Our object of study belongs to the literature on nonlinear structural impulse responses defined by comparing recursively simulated shocked and unshocked paths \citep{KoopPesaranPotter1996,GourierouxLee2025NonlinearIRF,Ballarin2025}. That literature establishes the role of the initial state, shock size, future innovations, and nonlinear transition. We take structural identification as maintained and focus on inference when the innovation distributions used in the recursion are estimated from generated residuals. Relative to the closest normal-rank construction, the main change is the empirical-residual estimator and its first-order, smoothing, bootstrap, and finite-simulation theory.
Nonparametric and semiparametric local projections estimate nonlinear responses through horizon-specific conditional-mean relations rather than by iterating a structural transition \citep{Jorda2005,PlagborgMollerWolf2021,JordaTaylor2025,GoncalvesHerreraKilianPesavento2024,GoncalvesHerreraKilianPesaventoHolban2026}. These methods address related questions through a different estimation route. The present results neither require nor imply a general ranking between recursive and local-projection procedures; they provide inference for the innovation-distribution component of a recursively defined structural response.
\subsection{Organization}
\label{sec:introduction-organization}
Section~2 defines the population response and the empirical-quantile estimator, with the smoothed estimator introduced as a paired comparison. Section~3 develops the generated-residual expansion, the recursive influence representation, and density-free feasible inference. Section~4 characterizes the effect of quantile smoothing under the same target. Section~5 studies the full recursive residual bootstrap and numerical integration error. Section~6 concludes. The appendices provide primitive conditions, derivations, proofs, and implementation conventions.
\section{Population response and empirical-quantile estimation}
\label{sec:model-estimators}
This section fixes the population impulse-response function before estimation. A scalar model first displays the two roles of the innovation quantile function. The vector formulation then defines the normal-rank structural response. The empirical-quantile estimator is introduced as the principal procedure, followed by a smoothed estimator that changes only the quantile estimate supplied to the same recursive map.
\subsection{A scalar model that displays the estimation problem}
\label{sec:scalar-example}
Consider the nonlinear location-scale autoregression
\begin{equation}
Y_t
=
\mu\left(Y_{t-1};\boldsymbol{\beta}_0\right)
+
\sigma\left(Y_{t-1};\boldsymbol{\beta}_0\right)U_t,
\qquad
\sigma\left(y;\boldsymbol{\beta}_0\right)>0,
\label{eq:scalar-model}
\end{equation}
where \(\boldsymbol{\beta}_0\) is finite-dimensional. The innovations are independent and identically distributed with continuous distribution function \(F_0\) and quantile function \(Q_0=F_0^{-1}\). The positivity of \(\sigma\) yields the residual map
\[
G\left(y_t,y_{t-1};\boldsymbol{\beta}\right)
=
\frac{y_t-\mu\left(y_{t-1};\boldsymbol{\beta}\right)}{\sigma\left(y_{t-1};\boldsymbol{\beta}\right)},
\]
so that \(U_t=G\left(Y_t,Y_{t-1};\boldsymbol{\beta}_0\right)\).
Let \(P_1=F_0\left(U_1\right)\), and define a shock of size \(\delta\in\mathbb{R}\) through
\begin{equation}
\tau_{\delta}\left(p\right)
=
\Phi\left(\Phi^{-1}\left(p\right)+\delta\right),
\qquad
p\in\left(0,1\right).
\label{eq:rank-shift}
\end{equation}
The impact innovations are \(Q_0\left(P_1\right)\) and \(Q_0\left(\tau_{\delta}\left(P_1\right)\right)\). Thus, \(\delta\) is measured on the standard-normal rank scale, while the change in structural innovation units depends on \(Q_0\) and the realized rank.
Fix an initial value \(y\) and let \(P_1,\ldots,P_h\) be independent uniform random variables. For a generic parameter \(\boldsymbol{\beta}\) and quantile function \(Q\), construct the unshocked and shocked paths from the same future ranks. Their difference satisfies
\begin{equation}
\begin{aligned}
D_1
&=
\sigma\left(y;\boldsymbol{\beta}\right)
\left[
Q\left(\tau_{\delta}\left(P_1\right)\right)-Q\left(P_1\right)
\right],\\
D_j
&=
\mu\left(Y_{j-1}^{\delta};\boldsymbol{\beta}\right)
-
\mu\left(Y_{j-1}^{0};\boldsymbol{\beta}\right)
+
\left[
\sigma\left(Y_{j-1}^{\delta};\boldsymbol{\beta}\right)
-
\sigma\left(Y_{j-1}^{0};\boldsymbol{\beta}\right)
\right]Q\left(P_j\right),
\qquad
j=2,\ldots,h.
\end{aligned}
\label{eq:scalar-response-recursion}
\end{equation}
Equation~\eqref{eq:scalar-response-recursion} displays the two uses of the unknown quantile function. It determines the impact difference and the future innovations through which that difference propagates. The scalar impulse-response function is \(\psi_h\left(y,\delta\right)=\mathbb{E}\left[D_h\right]\).
\subsection{Structural response in the general model}
\label{sec:general-response}
Let \(\boldsymbol{Y}_t\in\mathbb{R}^{d}\) satisfy
\begin{equation}
\boldsymbol{Y}_t
=
g\left(\boldsymbol{Y}_{t-1},\boldsymbol{U}_t;\boldsymbol{\beta}_0\right),
\label{eq:structural-transition}
\end{equation}
where \(\boldsymbol{U}_t\in\mathbb{R}^{n}\) is the vector of structural innovations. Higher-order dynamics are represented by augmenting the state. We assume that the structural representation is identified and that a known residual map \(G\) satisfies
\[
\boldsymbol{U}_t
=
G\left(\boldsymbol{Y}_t,\boldsymbol{Y}_{t-1};\boldsymbol{\beta}_0\right).
\]
The innovation vectors are independent and identically distributed over time. Their components are mutually independent, with continuous marginal distribution functions \(F_{j0}\) and quantile functions \(Q_{j0}\), for \(j=1,\ldots,n\). Define
\[
\boldsymbol{Q}_0\left(\boldsymbol{p}\right)
=
\left(
Q_{10}\left(p_1\right),
\ldots,
Q_{n0}\left(p_n\right)
\right)^\prime.
\]
Then \(\boldsymbol{U}_t=\boldsymbol{Q}_0\left(\boldsymbol{P}_t\right)\), where the components of \(\boldsymbol{P}_t\) are independent uniform random variables. For a selected component \(k\), let \(\boldsymbol{\tau}_{k,\delta}\left(\boldsymbol{p}\right)\) apply the map in equation~\eqref{eq:rank-shift} to \(p_k\) and leave the remaining coordinates unchanged.
Fix an initial state \(\boldsymbol{y}\), a response vector \(\boldsymbol{a}\), and a horizon \(h\). For \(r\in\left\{0,\delta\right\}\), set \(\boldsymbol{P}_1^0=\boldsymbol{P}_1\), \(\boldsymbol{P}_1^{\delta}=\boldsymbol{\tau}_{k,\delta}\left(\boldsymbol{P}_1\right)\), and \(\boldsymbol{P}_j^0=\boldsymbol{P}_j^{\delta}=\boldsymbol{P}_j\) for \(j\geqslant2\). For a generic \(\boldsymbol{\beta}\) and \(\boldsymbol{Q}\), define
\begin{equation}
\boldsymbol{Y}_0^r=\boldsymbol{y},
\qquad
\boldsymbol{Y}_j^r
=
g\left(\boldsymbol{Y}_{j-1}^r,\boldsymbol{Q}\left(\boldsymbol{P}_j^r\right);\boldsymbol{\beta}\right),
\qquad
j=1,\ldots,h.
\label{eq:paired-paths}
\end{equation}
The path response is
\[
D_h\left(\boldsymbol{P}_{1:h};\boldsymbol{\beta},\boldsymbol{Q},\boldsymbol{y},\boldsymbol{a},k,\delta\right)
=
\boldsymbol{a}^\prime\left(\boldsymbol{Y}_h^{\delta}-\boldsymbol{Y}_h^0\right),
\]
and the population impulse-response function is
\begin{equation}
\psi_h\left(\boldsymbol{y},\boldsymbol{a},k,\delta\right)
=
\mathbb{E}\left[
D_h\left(\boldsymbol{P}_{1:h};\boldsymbol{\beta}_0,\boldsymbol{Q}_0,\boldsymbol{y},\boldsymbol{a},k,\delta\right)
\right].
\label{eq:population-irf}
\end{equation}
The initial state, response variable, shocked component, shock map, and innovation law in equation~\eqref{eq:population-irf} remain fixed throughout the estimator comparison.
\subsection{Empirical-quantile estimator and smoothing comparator}
\label{sec:estimators}
Let \(\boldsymbol{\hat{\beta}}\) estimate \(\boldsymbol{\beta}_0\), and recover the structural residuals by
\[
\boldsymbol{\hat{U}}_t
=
G\left(\boldsymbol{Y}_t,\boldsymbol{Y}_{t-1};\boldsymbol{\hat{\beta}}\right),
\qquad
t=1,\ldots,T.
\]
For component \(j\), let \(\hat{U}_{j,\left(1\right)}\leqslant\cdots\leqslant\hat{U}_{j,\left(T\right)}\) denote the residual order statistics.
\subsubsection{Empirical-quantile estimator}
The empirical innovation quantile is
\[
\hat{Q}_{j}^{\mathrm{E}}\left(p\right)
=
\hat{U}_{j,\left(\lceil Tp\rceil\right)},
\qquad
p\in\left(0,1\right).
\]
This estimate uses the generated residual ranks directly and requires no smoothing parameter. Collect the componentwise estimates in \(\boldsymbol{\hat{Q}}^{\mathrm{E}}\).
\subsubsection{Paired smoothed estimator}
For comparison, apply a Gaussian rank smoother \(\mathcal{S}_{b_T}\) to the same residual order statistics,
\[
\hat{Q}_{j,b_T}^{\mathrm{S}}
=
\mathcal{S}_{b_T}\hat{Q}_{j}^{\mathrm{E}},
\]
where \(b_T\) is the bandwidth. Section~\ref{sec:smoothing} states the conditions on this operator. No smoothing choice enters the empirical estimator.
Let \(\boldsymbol{P}_{s,1:h}\), for \(s=1,\ldots,S\), be simulation draws independent of the sample. For \(m\in\left\{\mathrm{E},\mathrm{S}\right\}\), collect the componentwise quantile estimates in \(\boldsymbol{\hat{Q}}^{m}\) and define
\begin{equation}
\hat{\psi}_{h,S}^{m}\left(\boldsymbol{y},\boldsymbol{a},k,\delta\right)
=
\frac{1}{S}
\sum_{s=1}^{S}
D_h\left(
\boldsymbol{P}_{s,1:h};
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{m},
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right).
\label{eq:unified-estimator}
\end{equation}
Both estimators in equation~\eqref{eq:unified-estimator} use the same transition estimate, response specification, and simulation ranks. Their difference is therefore attributable to the quantile estimate supplied to the recursive response map. Common simulation ranks reduce numerical variation in this paired comparison; they do not replace sampling inference.
\section{First-order inference with empirical residual quantiles}
\label{sec:main-theory}
For a fixed response specification, write
\[
\Psi_h\left(\boldsymbol{\beta},\boldsymbol{Q}\right)
=
\mathbb{E}\left[
D_h\left(
\boldsymbol{P}_{1:h};
\boldsymbol{\beta},
\boldsymbol{Q},
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right)
\right]
\]
and let \(\Psi_{h,S}\left(\boldsymbol{\beta},\boldsymbol{Q}\right)\) denote the corresponding average over the \(S\) simulated rank paths. Thus, \(\psi_h=\Psi_h\left(\boldsymbol{\beta}_0,\boldsymbol{Q}_0\right)\) and \(\hat{\psi}_{h,S}^{\mathrm{E}}=\Psi_{h,S}\left(\boldsymbol{\hat{\beta}},\boldsymbol{\hat{Q}}^{\mathrm{E}}\right)\). The arguments \(\left(\boldsymbol{y},\boldsymbol{a},k,\delta\right)\) are suppressed in this section.
\begin{assumption}[First-order regularity]
\label{ass:first-order}
The following conditions hold for every response specification considered below.
\begin{enumerate}[label=(\roman*),leftmargin=2.2em]
\item The transition in equation~\eqref{eq:structural-transition} has a unique strictly stationary and ergodic solution. Its dependence and moments are sufficient for laws of large numbers and a joint central limit theorem for the influence sequences defined below. The maps \(g\) and \(G\) are continuously differentiable in the state, innovation, and parameter arguments used in the analysis. For any fixed maximum horizon, the derivatives of the paired paths are bounded by a random variable with a finite \(2+\eta\) moment for some \(\eta>0\).
\item The transition estimator is consistent and asymptotically linear:
\begin{equation}
\sqrt{T}\left(\boldsymbol{\hat{\beta}}-\boldsymbol{\beta}_0\right)
=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\boldsymbol{L}_t
+
o_{\mathbb{P}}\left(1\right),
\qquad
\mathbb{E}\left[\boldsymbol{L}_t\right]
=
\boldsymbol{0},
\label{eq:beta-linearization}
\end{equation}
where \(\left\{\boldsymbol{L}_t\right\}\) is stationary, has a finite \(2+\eta\) moment, and satisfies the joint central limit theorem in part~\textup{(i)}.
\item Each marginal innovation distribution \(F_{j0}\) has a positive, continuously differentiable density \(f_{j0}\) on the interior of its support. The tails of \(Q_{j0}\), the reciprocal densities, and the path derivatives are controlled so that the derivative and influence terms defined below are square integrable. These conditions also make \(\Psi_h\) continuous and first-order differentiable in a weighted uniform norm that controls the ranks generated by \(\tau_{\delta}\).
\item For each innovation component, the empirical quantile function computed from generated residuals satisfies
\begin{equation}
\left\|
\sqrt{T}\left(\hat{Q}_{j}^{\mathrm{E}}-Q_{j0}\right)
-
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\left[
\Xi_{j,t}
+
\boldsymbol{r}_j\left(\cdot\right)^{\prime}\boldsymbol{L}_t
\right]
\right\|_j
=
o_{\mathbb{P}}\left(1\right),
\label{eq:generated-quantile-expansion}
\end{equation}
where \(\|\cdot\|_j\) is the norm in part~\textup{(iii)}, \(\boldsymbol{r}_j\) is a deterministic vector-valued function, and
\[
\Xi_{j,t}\left(p\right)
=
\frac{
p-
\mathbf{1}\left\{
U_{jt}\leqslant Q_{j0}\left(p\right)
\right\}
}{
f_{j0}\left(Q_{j0}\left(p\right)\right)
},
\qquad
p\in\left(0,1\right).
\]
For the smooth residual maps considered here, the generated-residual correction can be written as
\[
\boldsymbol{r}_j\left(p\right)
=
-
\frac{
\boldsymbol{\Gamma}_j\left(Q_{j0}\left(p\right)\right)
}{
f_{j0}\left(Q_{j0}\left(p\right)\right)
},
\qquad
\boldsymbol{\Gamma}_j\left(u\right)
=
\left.
\frac{\partial}{\partial\boldsymbol{\beta}}
\mathbb{E}\left[
\mathbf{1}\left\{
G_j\left(
\boldsymbol{Y}_t,
\boldsymbol{Y}_{t-1};
\boldsymbol{\beta}
\right)
\leqslant u
\right\}
\right]
\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}_0}.
\]
Thus, \(\boldsymbol{r}_j\left(p\right)^{\prime}\boldsymbol{L}_t\) is the first-order effect of estimating \(\boldsymbol{\beta}_0\) before forming residual order statistics.
\item The number of response specifications and their horizons are fixed as \(T\) increases. The simulation draws are independent of the sample. Consistency requires \(S\rightarrow\infty\); the first-order results below impose \(S/T\rightarrow\infty\).
\end{enumerate}
\end{assumption}
Assumption~\ref{ass:first-order} separates the requirements imposed by the dynamic model, the transition estimator, the innovation distributions, and the generated residuals. The vector theorem below uses equation~\eqref{eq:generated-quantile-expansion} as a componentwise high-level condition. The next proposition verifies that condition under the primitive scalar assumptions and records the corresponding vector extension.
\begin{proposition}[Generated-residual empirical quantiles]
\label{prop:generated-residual-quantile-expansion}
Under Assumptions~\ref{ass:scalar-dynamics}--\ref{ass:scalar-estimator} and~\ref{ass:response-tail-continuity}, equation~\eqref{eq:generated-quantile-expansion} holds for the scalar model. In the vector model, the same conclusion holds jointly across components under the residual-process conditions in Lemma~\ref{lem:generated-residual-empirical-process}, the componentwise density conditions in Assumption~\ref{ass:first-order}\textup{(iii)}, and Assumption~\ref{ass:response-tail-continuity}.
\end{proposition}
Proposition~\ref{prop:generated-residual-quantile-expansion} is the bridge between residual empirical-process theory and the recursive response. Its proof is given in Appendix~\ref{app:residual-quantiles}; Appendix~\ref{app:assumptions} states the primitive scalar conditions in full.
We next introduce the derivatives needed for the first-order result. Let
\[
\Lambda_{h,s,j}^{r}
=
\left.
\frac{
\partial\left(
\boldsymbol{a}^{\prime}\boldsymbol{Y}_h^{r}
\right)
}{
\partial u_{j,s}^{r}
}
\right|_{\left(\boldsymbol{\beta}_0,\boldsymbol{Q}_0\right)},
\qquad
r\in\left\{0,\delta\right\},
\]
be the sensitivity of the date-\(h\) outcome along path \(r\) to structural innovation \(j\) at date \(s\). These sensitivities are obtained by differentiating the finite recursion in equation~\eqref{eq:paired-paths}. Also define the direct transition derivative
\[
\boldsymbol{A}_h
=
\mathbb{E}\left[
\left.
\frac{\partial}{\partial\boldsymbol{\beta}}
D_h\left(
\boldsymbol{P}_{1:h};
\boldsymbol{\beta},
\boldsymbol{Q}_0,
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right)
\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}_0}
\right],
\]
which holds the innovation quantile functions fixed.
For a collection of quantile perturbations \(\boldsymbol{q}=\left(q_1,\ldots,q_n\right)^{\prime}\), split the derivative with respect to \(\boldsymbol{Q}_0\) into
\begin{equation}
\begin{aligned}
\dot{\Psi}_{h}^{\mathrm{dist}}\left[\boldsymbol{q}\right]
&=
\mathbb{E}\left[
\sum_{j=1}^{n}
\left(
\mathbf{1}\left\{j\neq k\right\}
\Lambda_{h,1,j}^{\delta}
-
\Lambda_{h,1,j}^{0}
\right)
q_j\left(P_{j1}\right)
+
\sum_{s=2}^{h}
\sum_{j=1}^{n}
\left(
\Lambda_{h,s,j}^{\delta}
-
\Lambda_{h,s,j}^{0}
\right)
q_j\left(P_{js}\right)
\right],\\
\dot{\Psi}_{h}^{\mathrm{imp}}\left[\boldsymbol{q}\right]
&=
\mathbb{E}\left[
\Lambda_{h,1,k}^{\delta}
q_k\left(
\tau_{\delta}\left(P_{k1}\right)
\right)
\right].
\end{aligned}
\label{eq:response-quantile-derivatives}
\end{equation}
The first map covers the unshifted impact innovations and the common future innovations. The second isolates the shifted innovation in the shocked path at impact.
For later use, collect the empirical-quantile influence functions in
\[
\boldsymbol{\Xi}_t
=
\left(
\Xi_{1,t},
\ldots,
\Xi_{n,t}
\right)^{\prime},
\]
and define the generated-residual perturbation
\[
\boldsymbol{R}_t
=
\left(
p\mapsto
\boldsymbol{r}_1\left(p\right)^{\prime}\boldsymbol{L}_t,
\ldots,
p\mapsto
\boldsymbol{r}_n\left(p\right)^{\prime}\boldsymbol{L}_t
\right)^{\prime}.
\]
Before turning to the first-order result, consistency follows from the decomposition
\[
\begin{aligned}
\hat{\psi}_{h,S}^{\mathrm{E}}-\psi_h
&=
\left[
\Psi_{h,S}\left(
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
-
\Psi_h\left(
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
\right]\\
&\quad+
\left[
\Psi_h\left(
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
-
\Psi_h\left(
\boldsymbol{\beta}_0,
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
\right]\\
&\quad+
\left[
\Psi_h\left(
\boldsymbol{\beta}_0,
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
-
\Psi_h\left(
\boldsymbol{\beta}_0,
\boldsymbol{Q}_0
\right)
\right].
\end{aligned}
\]
The three terms are, respectively, Monte Carlo integration error, transition-estimation error, and residual-quantile estimation error.
\begin{proposition}[Consistency]
\label{prop:consistency}
Under Assumption~\ref{ass:first-order}\textup{(i)--(iv)}, if \(S\rightarrow\infty\), then for any fixed finite collection of response specifications,
\[
\max_{1\leqslant m\leqslant M}
\left|
\hat{\psi}_{m,S}^{\mathrm{E}}
-
\psi_m
\right|
\xrightarrow{\mathbb{P}}
0.
\]
If the smoother satisfies
\[
\left\|
\hat{Q}_{j,b_T}^{\mathrm{S}}
-
Q_{j0}
\right\|_j
\xrightarrow{\mathbb{P}}
0
\]
for every \(j\), the same conclusion holds for \(\hat{\psi}_{m,S}^{\mathrm{S}}\).
\end{proposition}
\begin{theorem}[Asymptotic linear representation]
\label{thm:main-linearization}
Suppose Assumption~\ref{ass:first-order} holds and \(S/T\rightarrow\infty\). For each fixed response specification,
\begin{equation}
\sqrt{T}\left(
\hat{\psi}_{h,S}^{\mathrm{E}}
-
\psi_h
\right)
=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
Z_{h,t}
+
o_{\mathbb{P}}\left(1\right),
\label{eq:empirical-irf-linearization}
\end{equation}
where the observation-level influence contribution is
\begin{equation}
\begin{aligned}
Z_{h,t}
&=
Z_{h,t}^{\mathrm{tr}}
+
Z_{h,t}^{\mathrm{res}}
+
Z_{h,t}^{\mathrm{dist}}
+
Z_{h,t}^{\mathrm{imp}},\\
Z_{h,t}^{\mathrm{tr}}
&=
\boldsymbol{A}_h^{\prime}\boldsymbol{L}_t,\\
Z_{h,t}^{\mathrm{res}}
&=
\dot{\Psi}_{h}^{\mathrm{dist}}\left[
\boldsymbol{R}_t
\right]
+
\dot{\Psi}_{h}^{\mathrm{imp}}\left[
\boldsymbol{R}_t
\right],\\
Z_{h,t}^{\mathrm{dist}}
&=
\dot{\Psi}_{h}^{\mathrm{dist}}\left[
\boldsymbol{\Xi}_t
\right],\\
Z_{h,t}^{\mathrm{imp}}
&=
\dot{\Psi}_{h}^{\mathrm{imp}}\left[
\boldsymbol{\Xi}_t
\right].
\end{aligned}
\label{eq:influence-decomposition}
\end{equation}
More generally, let \(m=1,\ldots,M\) index any fixed collection of horizons, initial states, response vectors, shocked components, and shock sizes. Stack the corresponding estimators, population responses, and influence contributions in \(\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}\), \(\boldsymbol{\psi}\), and \(\boldsymbol{Z}_t\). Then
\begin{equation}
\sqrt{T}\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
-
\boldsymbol{\psi}
\right)
=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\boldsymbol{Z}_t
+
o_{\mathbb{P}}\left(1\right).
\label{eq:joint-empirical-irf-linearization}
\end{equation}
\end{theorem}
Theorem~\ref{thm:main-linearization} preserves the covariance among all four sources of uncertainty because they are assembled at the observation level before the long-run covariance is computed. In particular, the direct and generated-residual terms share the transition influence \(\boldsymbol{L}_t\), while the distribution and impact terms are formed from the same residual observation.
\begin{corollary}[Gaussian limit]
\label{cor:main-gaussian-limit}
Under the conditions of Theorem~\ref{thm:main-linearization},
\begin{equation}
\sqrt{T}\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
-
\boldsymbol{\psi}
\right)
\xrightarrow{\mathrm{d}}
\mathcal{N}\left(
\boldsymbol{0},
\boldsymbol{\Omega}
\right),
\qquad
\boldsymbol{\Omega}
=
\sum_{\ell=-\infty}^{\infty}
\operatorname{Cov}\left(
\boldsymbol{Z}_0,
\boldsymbol{Z}_{\ell}
\right).
\label{eq:empirical-irf-limit}
\end{equation}
For one response, the asymptotic variance is the corresponding diagonal element of \(\boldsymbol{\Omega}\). If \(\left\{\boldsymbol{Z}_t\right\}\) is a martingale difference sequence, the sum reduces to
\[
\mathbb{E}\left[
\boldsymbol{Z}_t
\boldsymbol{Z}_t^{\prime}
\right].
\]
\end{corollary}
\subsection{Interpretation of the influence decomposition}
\label{sec:main-theory-interpretation}
Equation~\eqref{eq:influence-decomposition} keeps four uses of estimated objects separate while preserving their covariance at the observation level. Table~\ref{tab:influence-components} summarizes the source and scope of each component.
\begin{table}[t]
\centering
\caption{Sources of first-order uncertainty}
\label{tab:influence-components}
\begin{tabular}{p{2.4cm}p{7.0cm}p{5.0cm}}
\toprule
Component & Source & Vanishes when \\
\midrule
\(Z_{h,t}^{\mathrm{tr}}\) & The transition parameter changes the shocked and unshocked recursive paths directly. & \(\boldsymbol{\beta}_0\) is known. \\
\addlinespace
\(Z_{h,t}^{\mathrm{res}}\) & The transition estimate changes the generated residual values and their order statistics before the response is simulated. & \(\boldsymbol{\beta}_0\) is known, or \(\boldsymbol{Q}_0\) is treated as known. \\
\addlinespace
\(Z_{h,t}^{\mathrm{dist}}\) & The unshifted impact innovations and the common future innovations use estimated quantiles. & \(\boldsymbol{Q}_0\) is known. \\
\addlinespace
\(Z_{h,t}^{\mathrm{imp}}\) & The shocked impact innovation uses the estimated quantile evaluated at the shifted rank. & \(\boldsymbol{Q}_0\) is known, or the impact intervention is fixed in structural units. \\
\bottomrule
\end{tabular}
\end{table}
The direct and generated-residual components share the transition influence \(\boldsymbol{L}_t\), while the distribution and shifted-impact components are formed from the same residual observation. Computing their covariance after separate aggregation would therefore lose first-order cross terms. The decomposition also clarifies the nested cases: known transition parameters remove the first two components; known innovation distributions leave only direct transition uncertainty; and a fixed additive structural shock removes the separate shifted-impact component without generally removing distribution uncertainty from future innovations.
\subsection{Density-free feasible inference}
\label{sec:main-theory-feasible}
Theorem~\ref{thm:main-linearization} yields feasible inference once its observation-level contributions are estimated. A direct use of the empirical-quantile influence function in equation~\eqref{eq:generated-quantile-expansion} appears to require the innovation densities. Here those densities cancel after the quantile effects are integrated over the ranks entering the impulse-response function.
Let \(\rho_{\delta}=\tau_{\delta}^{-1}=\tau_{-\delta}\). If \(\phi\) denotes the standard normal density, then
\[
\rho_{\delta}\left(p\right)
=
\Phi\left(\Phi^{-1}\left(p\right)-\delta\right),
\qquad
\rho_{\delta}'\left(p\right)
=
\frac{\phi\left(\Phi^{-1}\left(p\right)-\delta\right)}
{\phi\left(\Phi^{-1}\left(p\right)\right)}.
\]
For each innovation component, define
\begin{equation}
\begin{aligned}
\omega_{h,j}^{\mathrm{dist}}\left(p\right)
&=
\mathbb{E}\left[
\left(
\mathbf{1}\left\{j\neq k\right\}\Lambda_{h,1,j}^{\delta}
-
\Lambda_{h,1,j}^{0}
\right)
\mathrel{\big|}
P_{j1}=p
\right]
+
\sum_{s=2}^{h}
\mathbb{E}\left[
\Lambda_{h,s,j}^{\delta}-\Lambda_{h,s,j}^{0}
\mathrel{\big|}
P_{js}=p
\right],\\
\omega_{h,k}^{\mathrm{imp}}\left(p\right)
&=
\mathbb{E}\left[
\Lambda_{h,1,k}^{\delta}
\mathrel{\big|}
P_{k1}=p
\right].
\end{aligned}
\label{eq:propagation-weights}
\end{equation}
The first weight collects the effects of the unshifted impact innovations and the common future innovations. The second collects the effect of the shifted impact innovation. Consequently,
\[
\dot{\Psi}_{h}^{\mathrm{dist}}\left[\boldsymbol{q}\right]
=
\sum_{j=1}^{n}
\int_{0}^{1}
\omega_{h,j}^{\mathrm{dist}}\left(p\right)
q_j\left(p\right)
\,\mathrm{d}p,
\qquad
\dot{\Psi}_{h}^{\mathrm{imp}}\left[\boldsymbol{q}\right]
=
\int_{0}^{1}
\omega_{h,k}^{\mathrm{imp}}\left(p\right)
q_k\left(\tau_{\delta}\left(p\right)\right)
\,\mathrm{d}p.
\]
Let
\[
\boldsymbol{J}_{j,t}
=
\left.
\frac{\partial}{\partial\boldsymbol{\beta}}
G_j\left(
\boldsymbol{Y}_t,
\boldsymbol{Y}_{t-1};
\boldsymbol{\beta}
\right)
\right|_{\boldsymbol{\beta}=\boldsymbol{\beta}_0}
\]
denote the derivative of structural residual \(j\) with respect to the transition parameter.
\begin{proposition}[Influence contributions without density estimation]
\label{prop:density-free-influence}
Suppose Assumption~\ref{ass:first-order} holds and the generated-residual correction is induced by the differentiable residual map \(G\). Then the residual, innovation-distribution, and shifted-impact terms in equation~\eqref{eq:influence-decomposition} satisfy
\begin{equation}
\begin{aligned}
\boldsymbol{B}_{h}^{\mathrm{res}}
&=
\sum_{j=1}^{n}
\mathbb{E}\left[
\omega_{h,j}^{\mathrm{dist}}\left(P_{jt}\right)
\boldsymbol{J}_{j,t}
\right]
+
\mathbb{E}\left[
\omega_{h,k}^{\mathrm{imp}}\left(
\rho_{\delta}\left(P_{kt}\right)
\right)
\rho_{\delta}'\left(P_{kt}\right)
\boldsymbol{J}_{k,t}
\right],\\
Z_{h,t}^{\mathrm{res}}
&=
\left(
\boldsymbol{B}_{h}^{\mathrm{res}}
\right)^{\prime}
\boldsymbol{L}_t,\\
Z_{h,t}^{\mathrm{dist}}
&=
\sum_{j=1}^{n}
\int_{\mathbb{R}}
\omega_{h,j}^{\mathrm{dist}}\left(
F_{j0}\left(u\right)
\right)
\left[
F_{j0}\left(u\right)
-
\mathbf{1}\left\{U_{jt}\leqslant u\right\}
\right]
\,\mathrm{d}u,\\
Z_{h,t}^{\mathrm{imp}}
&=
\int_{\mathbb{R}}
\omega_{h,k}^{\mathrm{imp}}\left(
\rho_{\delta}\left(
F_{k0}\left(u\right)
\right)
\right)
\rho_{\delta}'\left(
F_{k0}\left(u\right)
\right)
\left[
F_{k0}\left(u\right)
-
\mathbf{1}\left\{U_{kt}\leqslant u\right\}
\right]
\,\mathrm{d}u.
\end{aligned}
\label{eq:density-free-influence}
\end{equation}
\end{proposition}
The result follows by changing variables from ranks to innovation units in the two quantile derivatives. For the generated-residual term, differentiability of \(G\) gives
\[
\boldsymbol{r}_j\left(p\right)
=
\mathbb{E}\left[
\boldsymbol{J}_{j,t}
\mathrel{\big|}
U_{jt}=Q_{j0}\left(p\right)
\right].
\]
Appendix~\ref{app:residual-quantiles} gives the formal argument.
We now construct sample analogues. For a transition estimator defined by the estimating equation
\[
\frac{1}{T}
\sum_{t=1}^{T}
\boldsymbol{S}_t\left(
\boldsymbol{\hat{\beta}}
\right)
=
\boldsymbol{0},
\]
let
\[
\boldsymbol{\hat{H}}
=
\frac{1}{T}
\sum_{t=1}^{T}
\frac{
\partial
\boldsymbol{S}_t\left(
\boldsymbol{\hat{\beta}}
\right)
}{
\partial\boldsymbol{\beta}^{\prime}
},
\qquad
\boldsymbol{\hat{L}}_t
=
-
\boldsymbol{\hat{H}}^{-1}
\boldsymbol{S}_t\left(
\boldsymbol{\hat{\beta}}
\right).
\]
Other regular transition estimators enter through their estimated influence contributions. Holding \(\boldsymbol{\hat{Q}}^{\mathrm{E}}\) fixed, estimate the direct transition derivative by
\[
\boldsymbol{\hat{A}}_h
=
\frac{1}{S}
\sum_{s=1}^{S}
\left.
\frac{\partial}{\partial\boldsymbol{\beta}}
D_h\left(
\boldsymbol{P}_{s,1:h};
\boldsymbol{\beta},
\boldsymbol{\hat{Q}}^{\mathrm{E}},
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right)
\right|_{\boldsymbol{\beta}=\boldsymbol{\hat{\beta}}}.
\]
The derivative and the path sensitivities in equation~\eqref{eq:propagation-weights} are obtained by automatic differentiation through equation~\eqref{eq:paired-paths}, or by the equivalent finite recursive derivatives in Appendix~\ref{app:recursive-derivatives}. To estimate a propagation weight at rank \(p\), fix the relevant simulated rank at \(p\) and average over the remaining ranks. Thus, the conditional expectations in equation~\eqref{eq:propagation-weights} are numerical integrals over inputs controlled by the researcher rather than nonparametric regressions on observed data.
Let \(\hat{R}_{jt}\) be the rank of \(\hat{U}_{jt}\) among the component-\(j\) residuals, set
\[
\hat{P}_{jt}
=
\frac{\hat{R}_{jt}-1/2}{T},
\]
and define
\[
\Delta\hat{U}_{j,r}
=
\hat{U}_{j,\left(r+1\right)}
-
\hat{U}_{j,\left(r\right)}.
\]
Also let
\[
\boldsymbol{\hat{J}}_{j,t}
=
\left.
\frac{\partial}{\partial\boldsymbol{\beta}}
G_j\left(
\boldsymbol{Y}_t,
\boldsymbol{Y}_{t-1};
\boldsymbol{\beta}
\right)
\right|_{\boldsymbol{\beta}=\boldsymbol{\hat{\beta}}}.
\]
Estimate the generated-residual coefficient by
\[
\boldsymbol{\hat{B}}_{h}^{\mathrm{res}}
=
\frac{1}{T}
\sum_{t=1}^{T}
\left[
\sum_{j=1}^{n}
\hat{\omega}_{h,j}^{\mathrm{dist}}\left(
\hat{P}_{jt}
\right)
\boldsymbol{\hat{J}}_{j,t}
+
\hat{\omega}_{h,k}^{\mathrm{imp}}\left(
\rho_{\delta}\left(
\hat{P}_{kt}
\right)
\right)
\rho_{\delta}'\left(
\hat{P}_{kt}
\right)
\boldsymbol{\hat{J}}_{k,t}
\right].
\]
The four estimated contributions are
\begin{equation}
\begin{aligned}
\hat{Z}_{h,t}^{\mathrm{tr}}
&=
\boldsymbol{\hat{A}}_h^{\prime}
\boldsymbol{\hat{L}}_t,
\qquad
\hat{Z}_{h,t}^{\mathrm{res}}
=
\left(
\boldsymbol{\hat{B}}_{h}^{\mathrm{res}}
\right)^{\prime}
\boldsymbol{\hat{L}}_t,\\
\hat{Z}_{h,t}^{\mathrm{dist}}
&=
\sum_{j=1}^{n}
\sum_{r=1}^{T-1}
\Delta\hat{U}_{j,r}
\hat{\omega}_{h,j}^{\mathrm{dist}}\left(
\frac{r}{T}
\right)
\left[
\frac{r}{T}
-
\mathbf{1}\left\{
\hat{R}_{jt}\leqslant r
\right\}
\right],\\
\hat{Z}_{h,t}^{\mathrm{imp}}
&=
\sum_{r=1}^{T-1}
\Delta\hat{U}_{k,r}
\hat{\omega}_{h,k}^{\mathrm{imp}}\left(
\rho_{\delta}\left(
\frac{r}{T}
\right)
\right)
\rho_{\delta}'\left(
\frac{r}{T}
\right)
\left[
\frac{r}{T}
-
\mathbf{1}\left\{
\hat{R}_{kt}\leqslant r
\right\}
\right],\\
\hat{Z}_{h,t}
&=
\hat{Z}_{h,t}^{\mathrm{tr}}
+
\hat{Z}_{h,t}^{\mathrm{res}}
+
\hat{Z}_{h,t}^{\mathrm{dist}}
+
\hat{Z}_{h,t}^{\mathrm{imp}}.
\end{aligned}
\label{eq:estimated-influence-contributions}
\end{equation}
The two spacing sums evaluate the integrals in equation~\eqref{eq:density-free-influence} over the intervals between adjacent residual order statistics.
For a fixed collection of \(M\) horizons or response variables, apply equation~\eqref{eq:estimated-influence-contributions} to each response and stack the results in
\[
\boldsymbol{\hat{Z}}_t
=
\left(
\hat{Z}_{1,t},
\ldots,
\hat{Z}_{M,t}
\right)^{\prime}.
\]
Let
\[
\overline{\boldsymbol{\hat{Z}}}
=
\frac{1}{T}
\sum_{t=1}^{T}
\boldsymbol{\hat{Z}}_t,
\qquad
\boldsymbol{\tilde{Z}}_t
=
\boldsymbol{\hat{Z}}_t
-
\overline{\boldsymbol{\hat{Z}}},
\]
and define
\[
\boldsymbol{\hat{\Gamma}}_{\ell}
=
\frac{1}{T}
\sum_{t=\ell+1}^{T}
\boldsymbol{\tilde{Z}}_t
\boldsymbol{\tilde{Z}}_{t-\ell}^{\prime}.
\]
A feasible long-run covariance estimator is
\begin{equation}
\boldsymbol{\hat{\Omega}}
=
\boldsymbol{\hat{\Gamma}}_0
+
\sum_{\ell=1}^{q_T}
\mathcal{K}\left(
\frac{\ell}{q_T+1}
\right)
\left(
\boldsymbol{\hat{\Gamma}}_{\ell}
+
\boldsymbol{\hat{\Gamma}}_{\ell}^{\prime}
\right),
\label{eq:feasible-long-run-covariance}
\end{equation}
where \(\mathcal{K}\) is a standard HAC weighting function and \(q_T\) satisfies the corresponding lag-truncation conditions. When \(\left\{\boldsymbol{Z}_t\right\}\) is a martingale difference sequence, set \(q_T=0\), and equation~\eqref{eq:feasible-long-run-covariance} reduces to the sample covariance.
\begin{corollary}[Feasible Wald inference]
\label{cor:feasible-wald-inference}
Suppose the conditions of Corollary~\ref{cor:main-gaussian-limit} hold, the nuisance estimates above are consistent, and the numerical integration error is \(o_{\mathbb{P}}\left(1\right)\). Then
\[
\boldsymbol{\hat{\Omega}}
\xrightarrow{\mathbb{P}}
\boldsymbol{\Omega}.
\]
For response \(m\), an asymptotic pointwise \(1-\alpha\) confidence interval is
\begin{equation}
\mathcal{I}_{m,1-\alpha}^{\mathrm{pt}}
=
\left[
\hat{\psi}_{m,S}^{\mathrm{E}}
-
z_{1-\alpha/2}
\sqrt{
\frac{\hat{\Omega}_{mm}}{T}
},
\quad
\hat{\psi}_{m,S}^{\mathrm{E}}
+
z_{1-\alpha/2}
\sqrt{
\frac{\hat{\Omega}_{mm}}{T}
}
\right].
\label{eq:pointwise-wald-interval}
\end{equation}
For simultaneous inference over the fixed collection, let
\[
\boldsymbol{\hat{V}}
=
\operatorname{diag}\left(
\hat{\Omega}_{11},
\ldots,
\hat{\Omega}_{MM}
\right),
\qquad
\boldsymbol{\hat{C}}
=
\boldsymbol{\hat{V}}^{-1/2}
\boldsymbol{\hat{\Omega}}
\boldsymbol{\hat{V}}^{-1/2},
\]
and let \(c_{1-\alpha}\left(\boldsymbol{\hat{C}}\right)\) be the \(1-\alpha\) quantile of
\[
\max_{1\leqslant m\leqslant M}
\left|G_m\right|
\]
for
\[
\boldsymbol{G}
\sim
\mathcal{N}\left(
\boldsymbol{0},
\boldsymbol{\hat{C}}
\right).
\]
Simultaneous intervals are
\begin{equation}
\mathcal{I}_{m,1-\alpha}^{\mathrm{sim}}
=
\left[
\hat{\psi}_{m,S}^{\mathrm{E}}
-
c_{1-\alpha}\left(
\boldsymbol{\hat{C}}
\right)
\sqrt{
\frac{\hat{\Omega}_{mm}}{T}
},
\quad
\hat{\psi}_{m,S}^{\mathrm{E}}
+
c_{1-\alpha}\left(
\boldsymbol{\hat{C}}
\right)
\sqrt{
\frac{\hat{\Omega}_{mm}}{T}
}
\right],
\qquad
m=1,\ldots,M.
\label{eq:simultaneous-wald-intervals}
\end{equation}
\end{corollary}
The first-stage influence contributions and residual-map Jacobians are computed analytically from the transition estimator and structural model. Automatic differentiation is applied only to the smooth recursive paths; the order statistics are handled by equation~\eqref{eq:estimated-influence-contributions}. Expectations over rank paths are numerical integrals. Thus, the influence-based procedure requires no innovation-density estimator or quantile-smoothing bandwidth. Under the general weak-dependence formulation, the HAC lag choice concerns serial covariance in the observation-level contributions.
The joint result applies to a fixed number of horizons and response variables. Growing-horizon inference is outside the present claims. The intervals above also impose \(S/T\rightarrow\infty\); Section~\ref{sec:simulation-error} adds simulation uncertainty when \(S\) is proportional to \(T\).
\section{What smoothing changes under the same target}
\label{sec:smoothing}
The comparison in this section keeps the transition estimate and the simulated rank paths fixed. Thus, the paired difference between the two estimators is generated only by replacing \(\boldsymbol{\hat{Q}}^{\mathrm{E}}\) with \(\boldsymbol{\hat{Q}}_{b_T}^{\mathrm{S}}=\mathcal{S}_{b_T}\boldsymbol{\hat{Q}}^{\mathrm{E}}\), applied componentwise, in equation~\eqref{eq:unified-estimator}. For each innovation component,
\begin{equation}
\hat{Q}_{j,b_T}^{\mathrm{S}}-\hat{Q}_{j}^{\mathrm{E}}
=
\left(\mathcal{S}_{b_T}-\operatorname{Id}\right)Q_{j0}
+
\left(\mathcal{S}_{b_T}-\operatorname{Id}\right)
\left(\hat{Q}_{j}^{\mathrm{E}}-Q_{j0}\right).
\label{eq:quantile-smoothing-decomposition}
\end{equation}
The first term is the deterministic approximation error introduced by smoothing. The second records how smoothing changes the sampling error of the empirical quantile function.
\begin{assumption}[Rank smoothing]
\label{ass:smoothing}
The operator \(\mathcal{S}_{b}\) is a componentwise boundary-corrected Gaussian rank smoother with \(b_T\rightarrow0\). On interior ranks, its population version satisfies
\[
\left(\mathcal{S}_{b}q\right)\left(p\right)
=
\int_{\mathbb{R}}
\phi\left(v\right)
q\left(p+bv\right)
\,\mathrm{d}v,
\]
where \(\phi\) is the standard normal density. For each innovation component, \(Q_{j0}\) is twice continuously differentiable over the ranks used by the response, and
\[
\left\|
\mathcal{S}_{b}Q_{j0}
-
Q_{j0}
-
\frac{b^2}{2}Q_{j0}^{\prime\prime}
\right\|_j
=
o\left(b^2\right)
\]
under the weighted norm in Assumption~\ref{ass:first-order}. Moreover,
\[
\sqrt{T}
\left\|
\left(\mathcal{S}_{b_T}-\operatorname{Id}\right)
\left(\hat{Q}_{j}^{\mathrm{E}}-Q_{j0}\right)
\right\|_j
=
o_{\mathbb{P}}\left(1\right),
\qquad
j=1,\ldots,n.
\]
The boundary correction and tail conditions make the same expansion valid for the ranks reached by \(\tau_{\delta}\).
\end{assumption}
Assumption~\ref{ass:smoothing} separates the local smoothing bias from the first-order empirical-quantile fluctuation. To express its effect on the impulse-response function, use the propagation weights in equation~\eqref{eq:propagation-weights} and define
\begin{equation}
B_h^{\mathrm{S}}
=
\frac{1}{2}
\left[
\sum_{j=1}^{n}
\int_{0}^{1}
\omega_{h,j}^{\mathrm{dist}}\left(p\right)
Q_{j0}^{\prime\prime}\left(p\right)
\,\mathrm{d}p
+
\int_{0}^{1}
\omega_{h,k}^{\mathrm{imp}}\left(p\right)
Q_{k0}^{\prime\prime}\left(\tau_{\delta}\left(p\right)\right)
\,\mathrm{d}p
\right].
\label{eq:smoothing-bias-coefficient}
\end{equation}
The first integral covers the unshifted impact innovations and the common future innovations. The second covers the rank-shifted impact innovation.
\begin{theorem}[Paired smoothing expansion]
\label{thm:smoothing-expansion}
Suppose Assumptions~\ref{ass:first-order} and~\ref{ass:smoothing} hold, and let \(S\rightarrow\infty\). At any fixed horizon,
\begin{equation}
\hat{\psi}_{h,S}^{\mathrm{S}}
-
\hat{\psi}_{h,S}^{\mathrm{E}}
=
b_T^2B_h^{\mathrm{S}}
+
R_{h,T,S}^{\mathrm{S}},
\label{eq:smoothing-expansion}
\end{equation}
where
\begin{equation}
R_{h,T,S}^{\mathrm{S}}
=
o_{\mathbb{P}}\left(T^{-1/2}+b_T^2\right)
+
O_{\mathbb{P}}\left(
S^{-1/2}
\left(T^{-1/2}+b_T^2\right)
\right).
\label{eq:smoothing-remainder}
\end{equation}
The result holds jointly for any fixed collection of horizons, initial states, response variables, shocked components, and shock sizes after stacking the corresponding coefficients \(B_h^{\mathrm{S}}\).
\end{theorem}
The first term in equation~\eqref{eq:smoothing-remainder} contains the stochastic part of equation~\eqref{eq:quantile-smoothing-decomposition} and the higher-order terms from the nonlinear recursion. The second is the conditional Monte Carlo error in the paired estimator difference. Common simulation draws make this error proportional to the distance between the two quantile estimates. Consequently, \(S\rightarrow\infty\) is sufficient for the paired expansion, although inference for either estimator without a simulation correction continues to require \(S/T\rightarrow\infty\).
Because both procedures use the same \(\boldsymbol{\hat{\beta}}\) and the same generated residuals, there is no separate direct first-stage term in equation~\eqref{eq:smoothing-expansion}. The common transition and residual-quantile uncertainty remains in the sampling variance of each estimator; smoothing changes their leading difference through \(b_T^2B_h^{\mathrm{S}}\).
Equation~\eqref{eq:smoothing-bias-coefficient} also distinguishes pointwise quantile approximation from its effect after recursive propagation. The local error is weighted by the sensitivity of the date-\(h\) response to an innovation at that rank. When the relevant norms are finite,
\[
\left|B_h^{\mathrm{S}}\right|
\leqslant
\frac{1}{2}
\left[
\sum_{j=1}^{n}
\left\|\omega_{h,j}^{\mathrm{dist}}\right\|_1
\left\|Q_{j0}^{\prime\prime}\right\|_{\infty}
+
\left\|\omega_{h,k}^{\mathrm{imp}}\right\|_1
\left\|Q_{k0}^{\prime\prime}\right\|_{\infty}
\right].
\]
Thus, a small local error can matter when the recursion is persistent or strongly state dependent, while errors at different ranks may offset one another. If \(B_h^{\mathrm{S}}=0\), the second-order bias cancels after propagation and the next nonzero term determines the relevant bandwidth condition.
\begin{corollary}[Bandwidth regimes]
\label{cor:smoothing-regimes}
Let \(\lambda_T=\sqrt{T}b_T^2\), let \(\boldsymbol{B}^{\mathrm{S}}\) stack the coefficients in equation~\eqref{eq:smoothing-bias-coefficient} for a fixed collection of responses, and suppose \(S/T\rightarrow\infty\).
\begin{enumerate}[label=(\roman*),leftmargin=2.2em]
\item If \(\lambda_T\rightarrow0\), then
\[
\sqrt{T}
\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{S}}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)
\xrightarrow{\mathbb{P}}
\boldsymbol{0}.
\]
The empirical and smoothed estimators have the same first-order distribution and the same feasible covariance matrix.
\item If \(\lambda_T\rightarrow\lambda\in\left(0,\infty\right)\), then
\begin{equation}
\sqrt{T}
\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{S}}
-
\boldsymbol{\psi}
\right)
\xrightarrow{\mathrm{d}}
\mathcal{N}
\left(
\lambda\boldsymbol{B}^{\mathrm{S}},
\boldsymbol{\Omega}
\right).
\label{eq:smoothed-local-bias-limit}
\end{equation}
The smoothed estimator has the same first-order covariance as the empirical estimator, but its limiting distribution is shifted by the propagated smoothing bias.
\item If \(\lambda_T\rightarrow\infty\) and \(B_h^{\mathrm{S}}\neq0\), then
\begin{equation}
b_T^{-2}
\left(
\hat{\psi}_{h,S}^{\mathrm{S}}
-
\psi_h
\right)
\xrightarrow{\mathbb{P}}
B_h^{\mathrm{S}}.
\label{eq:smoothing-bias-dominates}
\end{equation}
The smoothed estimator remains consistent because \(b_T\rightarrow0\), but smoothing bias dominates its \(T^{-1/2}\) sampling error.
\end{enumerate}
\end{corollary}
For the second-order smoother used here, first-order equivalence requires \(b_T=o\left(T^{-1/4}\right)\). A bandwidth satisfying \(b_T\sim cT^{-1/4}\) produces the mean shift in equation~\eqref{eq:smoothed-local-bias-limit} with \(\lambda=c^2\), while \(T^{-1/4}=o\left(b_T\right)\) produces the bias-dominated behavior in equation~\eqref{eq:smoothing-bias-dominates}. For an order-\(r\) smoother, the same argument replaces \(b_T^2\) with \(b_T^r\), so the boundary between negligible and first-order bias is \(T^{-1/(2r)}\).
For one response, let \(\sigma_h^2\) be the corresponding diagonal element of \(\boldsymbol{\Omega}\). Under the local-bias regime, a nominal \(1-\alpha\) Wald interval centered at \(\hat{\psi}_{h,S}^{\mathrm{S}}\) and using the standard error from Section~\ref{sec:main-theory-feasible} has limiting coverage
\begin{equation}
\Phi\left(
z_{1-\alpha/2}
-
\frac{\lambda B_h^{\mathrm{S}}}{\sigma_h}
\right)
-
\Phi\left(
-z_{1-\alpha/2}
-
\frac{\lambda B_h^{\mathrm{S}}}{\sigma_h}
\right).
\label{eq:smoothed-wald-coverage}
\end{equation}
This equals \(1-\alpha\) when the propagated bias vanishes and is smaller otherwise. In the bias-dominated regime, the same uncorrected interval has limiting coverage zero. Bias correction is possible if \(b_T^2B_h^{\mathrm{S}}\) can be estimated with error \(o_{\mathbb{P}}\left(T^{-1/2}\right)\), but it requires additional smoothness estimation that the empirical estimator does not use.
\begin{corollary}[Additive shocks and affine transitions]
\label{cor:smoothing-special-cases}
The following cases verify the two channels in equation~\eqref{eq:smoothing-bias-coefficient}.
\begin{enumerate}[label=(\roman*),leftmargin=2.2em]
\item Suppose the impact intervention adds a fixed \(\xi\boldsymbol{e}_k\) to the structural innovation. Write \(\Lambda_{h,s,j}^{\xi}\) for the path derivative along the shocked path. The leading smoothing coefficient becomes
\[
\frac{1}{2}
\sum_{j=1}^{n}
\int_{0}^{1}
\left[
\mathbb{E}\left[
\Lambda_{h,1,j}^{\xi}-\Lambda_{h,1,j}^{0}
\mathrel{\big|}
P_{j1}=p
\right]
+
\sum_{s=2}^{h}
\mathbb{E}\left[
\Lambda_{h,s,j}^{\xi}-\Lambda_{h,s,j}^{0}
\mathrel{\big|}
P_{js}=p
\right]
\right]
Q_{j0}^{\prime\prime}\left(p\right)
\,\mathrm{d}p.
\]
There is no separate shifted-impact term because the amount added in structural innovation units does not depend on an estimated shifted quantile.
\item Suppose
\[
g\left(\boldsymbol{y},\boldsymbol{u};\boldsymbol{\beta}_0\right)
=
\boldsymbol{c}_0
+
\boldsymbol{A}_0\boldsymbol{y}
+
\boldsymbol{B}_0\boldsymbol{u}.
\]
Under the rank-normalized shock,
\begin{equation}
B_h^{\mathrm{S}}
=
\frac{1}{2}
\boldsymbol{a}^{\prime}
\boldsymbol{A}_0^{h-1}
\boldsymbol{B}_0
\boldsymbol{e}_k
\int_{0}^{1}
\left[
Q_{k0}^{\prime\prime}\left(\tau_{\delta}\left(p\right)\right)
-
Q_{k0}^{\prime\prime}\left(p\right)
\right]
\,\mathrm{d}p.
\label{eq:affine-smoothing-bias}
\end{equation}
Common future innovations cancel from the path difference, so smoothing matters only through the impact innovation. Under a fixed additive shock, the affine response does not depend on \(\boldsymbol{Q}_0\), and the empirical and smoothed estimators coincide exactly when they use the same transition estimate.
\end{enumerate}
\end{corollary}
Alternative rank maps and counterfactual future innovation distributions define different population responses rather than alternative estimators of equation~\eqref{eq:population-irf}. They are therefore outside the same-target smoothing comparison.
\section{Full recursive bootstrap and numerical integration}
\label{sec:inference}
\subsection{Full-model re-estimation}
\label{sec:bootstrap}
The residual bootstrap must reproduce every estimated object that enters the recursive impulse response. Throughout this section, \(\boldsymbol{\hat{Q}}^{\mathrm{E}}\) includes any location or scale normalization used to identify the structural innovations. When no finite-sample adjustment is required, it is the empirical quantile function defined in Section~\ref{sec:estimators}.
\begin{algorithm}[t]
\small
\caption{Full recursive residual bootstrap}
\label{alg:recursive-residual-bootstrap}
\KwIn{The sample \(\left\{\boldsymbol{Y}_t\right\}_{t=0}^{T}\), the fitted transition \(\boldsymbol{\hat{\beta}}\), the residual quantiles \(\boldsymbol{\hat{Q}}^{\mathrm{E}}\), the response specification \(\left(\boldsymbol{y},\boldsymbol{a},k,\delta,h\right)\), the response-simulation ranks \(\left\{\boldsymbol{P}_{s,1:h}\right\}_{s=1}^{S}\), a burn-in length \(\ell_T\), and \(B\) bootstrap replications.}
\For{\(b=1,\ldots,B\)}{
Draw \(W_{jt,b}^{*}\stackrel{\mathrm{iid}}{\sim}\operatorname{Unif}\left(0,1\right)\) independently over \(j=1,\ldots,n\) and \(t=1-\ell_T,\ldots,T\), and set \(\boldsymbol{U}_{t,b}^{*}=\boldsymbol{\hat{Q}}^{\mathrm{E}}\left(\boldsymbol{W}_{t,b}^{*}\right)\).\;
Set \(\boldsymbol{Y}_{-\ell_T,b}^{*}=\boldsymbol{Y}_0\) and generate
\[
\boldsymbol{Y}_{t,b}^{*}
=
g\left(
\boldsymbol{Y}_{t-1,b}^{*},
\boldsymbol{U}_{t,b}^{*};
\boldsymbol{\hat{\beta}}
\right)
\]
recursively for \(t=1-\ell_T,\ldots,T\); retain \(\left\{\boldsymbol{Y}_{t,b}^{*}\right\}_{t=0}^{T}\).\;
Re-estimate the transition by the original procedure to obtain \(\boldsymbol{\hat{\beta}}_{b}^{*}\).\;
Compute
\[
\boldsymbol{\hat{U}}_{t,b}^{*}
=
G\left(
\boldsymbol{Y}_{t,b}^{*},
\boldsymbol{Y}_{t-1,b}^{*};
\boldsymbol{\hat{\beta}}_{b}^{*}
\right)
\]
and reconstruct \(\boldsymbol{\hat{Q}}_{b}^{\mathrm{E},*}\) from the bootstrap residual order statistics.\;
Using the original ranks \(\left\{\boldsymbol{P}_{s,1:h}\right\}_{s=1}^{S}\), compute \(\hat{\psi}_{h,S,b}^{\mathrm{E},*}\) from equation~\eqref{eq:unified-estimator}; when studentization is used, recompute \(\boldsymbol{\hat{\Omega}}_{b}^{*}\) by the procedure in Section~\ref{sec:main-theory-feasible}.\;
}
\KwOut{The bootstrap estimates \(\left\{\hat{\psi}_{h,S,b}^{\mathrm{E},*},\boldsymbol{\hat{\Omega}}_{b}^{*}\right\}_{b=1}^{B}\).}
\end{algorithm}
The ranks used to generate each bootstrap time series are redrawn in every replication. In contrast, the response-simulation ranks are held fixed across the original and bootstrap estimates. Holding them fixed removes avoidable numerical noise from the comparison when \(S/T\rightarrow\infty\). Because the maintained model assumes independence across structural components, the algorithm draws the component ranks independently. Contemporaneously dependent innovations would instead require resampling from an estimated joint distribution.
Let \(\mathcal{F}_T\) contain the observed sample and the fixed response-simulation ranks, and let \(\mathbb{P}^*\) and \(\mathbb{E}^*\) denote probability and expectation conditional on \(\mathcal{F}_T\).
\begin{assumption}[Bootstrap regularity]
\label{ass:bootstrap}
The following conditions hold.
\begin{enumerate}[label=(\roman*),leftmargin=2.2em]
\item The stability, differentiability, and moment conditions in Assumption~\ref{ass:first-order} hold uniformly over a neighborhood of \(\boldsymbol{\beta}_0\) and over marginal innovation laws in a neighborhood of \(\boldsymbol{Q}_0\). The same location and scale normalizations are imposed in the original and bootstrap samples. The initialization and burn-in convention in Algorithm~\ref{alg:recursive-residual-bootstrap} affect the bootstrap estimator by \(o_{\mathbb{P}^*}\left(T^{-1/2}\right)\) in probability.
\item Conditional on \(\mathcal{F}_T\), the re-estimated transition parameter satisfies
\[
\sqrt{T}
\left(
\boldsymbol{\hat{\beta}}^*
-
\boldsymbol{\hat{\beta}}
\right)
=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\boldsymbol{L}_t^*
+
o_{\mathbb{P}^*}\left(1\right),
\]
and its conditional law converges in probability to the same Gaussian limit as the right-hand side of equation~\eqref{eq:beta-linearization}.
\item Jointly with part~\textup{(ii)}, the empirical quantile process formed from the re-estimated bootstrap residuals satisfies the conditional counterpart of equation~\eqref{eq:generated-quantile-expansion}:
\[
\left\|
\sqrt{T}
\left(
\hat{Q}_{j}^{\mathrm{E},*}
-
\hat{Q}_{j}^{\mathrm{E}}
\right)
-
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\left[
\Xi_{j,t}^*
+
\boldsymbol{r}_j^*\left(\cdot\right)^{\prime}
\boldsymbol{L}_t^*
\right]
\right\|_j
=
o_{\mathbb{P}^*}\left(1\right)
\]
for \(j=1,\ldots,n\). The stacked bootstrap influence array satisfies the corresponding conditional central limit theorem, and its long-run covariance converges to \(\boldsymbol{\Omega}\).
\item The number of response specifications is fixed and \(S/T\rightarrow\infty\).
\end{enumerate}
\end{assumption}
Primitive sufficient conditions are stated in Appendix~\ref{app:assumptions}, and the proof is given in Appendix~\ref{app:proofs-inference}. The key requirement is joint reproduction of the first-stage expansion and the residual empirical process.
\begin{theorem}[Validity of the full recursive residual bootstrap]
\label{thm:bootstrap-validity}
Suppose Assumptions~\ref{ass:first-order} and~\ref{ass:bootstrap} hold. For any fixed collection of response specifications,
\begin{equation}
\begin{aligned}
\sqrt{T}
\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E},*}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)
&=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\boldsymbol{Z}_t^*
+
o_{\mathbb{P}^*}\left(1\right),\\
\boldsymbol{Z}_t^*
&=
\boldsymbol{Z}_t^{\mathrm{tr},*}
+
\boldsymbol{Z}_t^{\mathrm{res},*}
+
\boldsymbol{Z}_t^{\mathrm{dist},*}
+
\boldsymbol{Z}_t^{\mathrm{imp},*},
\end{aligned}
\label{eq:bootstrap-linearization}
\end{equation}
where the four terms are the bootstrap analogues of the contributions in equation~\eqref{eq:influence-decomposition}. Conditional on \(\mathcal{F}_T\),
\begin{equation}
\sqrt{T}
\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E},*}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)
\xrightarrow{\mathrm{d}^*}
\mathcal{N}\left(
\boldsymbol{0},
\boldsymbol{\Omega}
\right)
\quad
\text{in probability}.
\label{eq:bootstrap-conditional-limit}
\end{equation}
For one response with \(\Omega_{hh}>0\), let
\[
\hat{\sigma}_h^2
=
\hat{\Omega}_{hh},
\qquad
\left(
\hat{\sigma}_h^*
\right)^2
=
\hat{\Omega}_{hh}^*,
\]
and define
\begin{equation}
R_{h,T}
=
\frac{
\sqrt{T}
\left(
\hat{\psi}_{h,S}^{\mathrm{E}}
-
\psi_h
\right)
}{
\hat{\sigma}_h
},
\qquad
R_{h,T}^*
=
\frac{
\sqrt{T}
\left(
\hat{\psi}_{h,S}^{\mathrm{E},*}
-
\hat{\psi}_{h,S}^{\mathrm{E}}
\right)
}{
\hat{\sigma}_h^*
}.
\label{eq:bootstrap-studentized-root}
\end{equation}
Then
\[
\sup_{x\in\mathbb{R}}
\left|
\mathbb{P}^*\left(
R_{h,T}^*
\leqslant x
\right)
-
\mathbb{P}\left(
R_{h,T}
\leqslant x
\right)
\right|
\xrightarrow{\mathbb{P}}
0.
\]
If every marginal variance is positive, the same conclusion holds for the maximum absolute studentized statistic over any fixed collection of responses.
\end{theorem}
Re-estimation of the transition generates the direct parameter term and the generated-residual term in equation~\eqref{eq:bootstrap-linearization}. Reconstructing the residual quantiles generates the future-distribution and shifted-impact terms. Consequently, the proof establishes a conditional version of the complete expansion in Theorem~\ref{thm:main-linearization}, followed by the bootstrap functional delta method for the finite-horizon recursive map \citep{VanDerVaartWellner2023,BeutnerZaehle2016}. Bootstrap quantiles of \(R_{h,T}^*\) therefore yield valid percentile-\(t\) intervals, and the fixed-dimensional maximum statistic yields simultaneous intervals.
We next formalize what is learned when only the response-simulation paths are redrawn. Define the fitted response with exact integration by
\[
\hat{\psi}_h^{\infty}
=
\Psi_h\left(
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right),
\]
and let
\[
\hat{\psi}_{h,S}^{\dagger}
=
\frac{1}{S}
\sum_{s=1}^{S}
D_h\left(
\boldsymbol{P}_{s,1:h}^{\dagger};
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}},
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right),
\]
where the ranks \(\boldsymbol{P}_{s,1:h}^{\dagger}\) are newly drawn while the fitted transition and residual quantiles remain fixed. Let \(\mathbb{P}^{\dagger}\) denote probability over these ranks conditional on \(\mathcal{F}_T\).
\begin{proposition}[Path-only resampling]
\label{prop:path-only-resampling}
Suppose Assumption~\ref{ass:first-order} holds and the fitted path response has a finite conditional second moment. Conditional on \(\mathcal{F}_T\),
\begin{equation}
\sqrt{S}
\left(
\hat{\psi}_{h,S}^{\dagger}
-
\hat{\psi}_h^{\infty}
\right)
\xrightarrow{\mathrm{d}^{\dagger}}
\mathcal{N}\left(
0,
\mathcal{V}_h
\right)
\quad
\text{in probability},
\label{eq:path-only-limit}
\end{equation}
where
\[
\mathcal{V}_h
=
\operatorname{Var}
\left[
D_h\left(
\boldsymbol{P}_{1:h};
\boldsymbol{\beta}_0,
\boldsymbol{Q}_0,
\boldsymbol{y},
\boldsymbol{a},
k,
\delta
\right)
\right].
\]
Consequently, if \(S/T\rightarrow\infty\), then
\[
\sqrt{T}
\left(
\hat{\psi}_{h,S}^{\dagger}
-
\hat{\psi}_h^{\infty}
\right)
\xrightarrow{\mathbb{P}^{\dagger}}
0.
\]
If \(S/T\rightarrow\kappa\in\left(0,\infty\right)\), then the same quantity converges conditionally to
\[
\mathcal{N}\left(
0,
\frac{\mathcal{V}_h}{\kappa}
\right).
\]
Neither limit contains the sampling variance \(\Omega_{hh}\).
\end{proposition}
Accordingly, path-only intervals quantify numerical integration error conditional on the fitted transition and residual distribution. They omit transition estimation, generated residuals, and estimation of both uses of the innovation quantiles. When \(S/T\rightarrow\infty\), such intervals collapse even though the estimator retains sampling uncertainty of order \(T^{-1/2}\). When \(S/T\) has a finite limit, they recover only the additional simulation component derived in Section~\ref{sec:simulation-error}. Centering path-only draws at the original finite-\(S\) estimate does not restore sampling variation; it combines the numerical errors from two simulated averages. Path-only resampling is therefore a Monte Carlo diagnostic, while sampling intervals require the influence-function procedure or Algorithm~\ref{alg:recursive-residual-bootstrap}.
\subsection{A finite number of simulated paths}
\label{sec:simulation-error}
The estimator in equation~\eqref{eq:unified-estimator} combines sampling uncertainty with numerical error from approximating the population expectation by \(S\) simulated paths. The preceding results impose \(S/T\rightarrow\infty\), which makes the second component negligible. We now allow \(S\) to grow at the same rate as \(T\).
Let \(m=1,\ldots,M\) index a fixed collection of response specifications, let \(H=\max_{1\leqslant m\leqslant M}h_m\), and use the same master rank path \(\boldsymbol{P}_{s,1:H}\) to evaluate all responses in simulation draw \(s\). Collect the corresponding path responses in
\[
\boldsymbol{D}\left(
\boldsymbol{P}_{1:H};
\boldsymbol{\beta},
\boldsymbol{Q}
\right)
=
\left(
D_1\left(
\boldsymbol{P}_{1:h_1};
\boldsymbol{\beta},
\boldsymbol{Q}
\right),
\ldots,
D_M\left(
\boldsymbol{P}_{1:h_M};
\boldsymbol{\beta},
\boldsymbol{Q}
\right)
\right)^{\prime},
\]
where the fixed initial states, response variables, shocked components, and shock sizes are suppressed. Define
\begin{equation}
\boldsymbol{\Omega}_{\mathrm{MC}}
=
\operatorname{Var}\left[
\boldsymbol{D}\left(
\boldsymbol{P}_{1:H};
\boldsymbol{\beta}_0,
\boldsymbol{Q}_0
\right)
\right].
\label{eq:mc-covariance}
\end{equation}
The off-diagonal elements of \(\boldsymbol{\Omega}_{\mathrm{MC}}\) record the covariance across horizons and response variables induced by the common rank paths.
\begin{theorem}[Sampling and simulation uncertainty]
\label{thm:finite-simulation-limit}
Suppose Assumption~\ref{ass:first-order}\textup{(i)--(iv)} holds, the collection of responses is fixed, and
\[
\mathbb{E}\left[
\left\|
\boldsymbol{D}\left(
\boldsymbol{P}_{1:H};
\boldsymbol{\beta}_0,
\boldsymbol{Q}_0
\right)
\right\|^{2+\eta}
\right]
<
\infty
\]
for some \(\eta>0\). If the simulation draws are independent of the observed sample and
\[
\frac{S}{T}
\rightarrow
\kappa
\in
\left(0,\infty\right),
\]
then
\begin{equation}
\begin{aligned}
\sqrt{T}\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
-
\boldsymbol{\psi}
\right)
&=
\frac{1}{\sqrt{T}}
\sum_{t=1}^{T}
\boldsymbol{Z}_t\\
&\quad+
\sqrt{\frac{T}{S}}
\frac{1}{\sqrt{S}}
\sum_{s=1}^{S}
\left[
\boldsymbol{D}\left(
\boldsymbol{P}_{s,1:H};
\boldsymbol{\beta}_0,
\boldsymbol{Q}_0
\right)
-
\boldsymbol{\psi}
\right]
+
o_{\mathbb{P}}\left(1\right).
\end{aligned}
\label{eq:finite-simulation-expansion}
\end{equation}
Consequently,
\begin{equation}
\sqrt{T}\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
-
\boldsymbol{\psi}
\right)
\xrightarrow{\mathrm{d}}
\mathcal{N}\left(
\boldsymbol{0},
\boldsymbol{\Omega}
+
\kappa^{-1}
\boldsymbol{\Omega}_{\mathrm{MC}}
\right).
\label{eq:finite-simulation-limit}
\end{equation}
The sampling and simulation components in equation~\eqref{eq:finite-simulation-expansion} are asymptotically independent.
\end{theorem}
When \(S/T\rightarrow\infty\), the simulation term vanishes and equation~\eqref{eq:finite-simulation-limit} reduces to Corollary~\ref{cor:main-gaussian-limit}. When \(S/T\rightarrow\kappa\in\left(0,\infty\right)\), simulation error contributes at the same order as estimation error. If \(S/T\rightarrow0\) and \(\boldsymbol{\Omega}_{\mathrm{MC}}\) is nonzero, simulation error dominates on the \(\sqrt{T}\) scale, and the leading rate becomes \(\sqrt{S}\).
For each simulated path, define the fitted response vector
\[
\boldsymbol{\hat{D}}_{s}^{\mathrm{E}}
=
\left(
D_m\left(
\boldsymbol{P}_{s,1:h_m};
\boldsymbol{\hat{\beta}},
\boldsymbol{\hat{Q}}^{\mathrm{E}}
\right)
\right)_{m=1}^{M},
\qquad
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
=
\frac{1}{S}
\sum_{s=1}^{S}
\boldsymbol{\hat{D}}_{s}^{\mathrm{E}}.
\]
The simulation covariance is estimated by
\[
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{E}}
=
\frac{1}{S-1}
\sum_{s=1}^{S}
\left(
\boldsymbol{\hat{D}}_{s}^{\mathrm{E}}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)
\left(
\boldsymbol{\hat{D}}_{s}^{\mathrm{E}}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)^{\prime}.
\]
Combining this matrix with the sampling covariance estimator from equation~\eqref{eq:feasible-long-run-covariance} gives
\begin{equation}
\boldsymbol{\hat{\Omega}}_{T,S}^{\mathrm{tot}}
=
\boldsymbol{\hat{\Omega}}
+
\frac{T}{S}
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{E}},
\qquad
\operatorname{Var}\left(
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}
\right)
\approx
\frac{\boldsymbol{\hat{\Omega}}}{T}
+
\frac{\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{E}}}{S}.
\label{eq:finite-simulation-correction}
\end{equation}
Thus, the pointwise standard error for response \(m\) is
\[
\operatorname{se}_{m,T,S}
=
\left(
\frac{\hat{\Omega}_{mm}}{T}
+
\frac{\hat{\Omega}_{\mathrm{MC},mm}^{\mathrm{E}}}{S}
\right)^{1/2}.
\]
Pointwise and simultaneous intervals are obtained by replacing \(\boldsymbol{\hat{\Omega}}\) with \(\boldsymbol{\hat{\Omega}}_{T,S}^{\mathrm{tot}}\) in equations~\eqref{eq:pointwise-wald-interval} and~\eqref{eq:simultaneous-wald-intervals}.
Common random numbers are particularly useful when comparing the empirical and smoothed estimators. For \(r,q\in\left\{\mathrm{E},\mathrm{S}\right\}\), define the cross-covariance
\[
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{rq}
=
\frac{1}{S-1}
\sum_{s=1}^{S}
\left(
\boldsymbol{\hat{D}}_{s}^{r}
-
\boldsymbol{\hat{\psi}}_{S}^{r}
\right)
\left(
\boldsymbol{\hat{D}}_{s}^{q}
-
\boldsymbol{\hat{\psi}}_{S}^{q}
\right)^{\prime}.
\]
Let
\[
\boldsymbol{\hat{\Delta}}_{s}
=
\boldsymbol{\hat{D}}_{s}^{\mathrm{S}}
-
\boldsymbol{\hat{D}}_{s}^{\mathrm{E}},
\qquad
\boldsymbol{\hat{\Delta}}_{S}
=
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{S}}
-
\boldsymbol{\hat{\psi}}_{S}^{\mathrm{E}}.
\]
The simulation covariance of the paired estimator difference is
\begin{equation}
\begin{aligned}
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\Delta}
&=
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{SS}}
+
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{EE}}
-
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{SE}}
-
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{ES}}\\
&=
\frac{1}{S-1}
\sum_{s=1}^{S}
\left(
\boldsymbol{\hat{\Delta}}_{s}
-
\boldsymbol{\hat{\Delta}}_{S}
\right)
\left(
\boldsymbol{\hat{\Delta}}_{s}
-
\boldsymbol{\hat{\Delta}}_{S}
\right)^{\prime}.
\end{aligned}
\label{eq:paired-mc-covariance}
\end{equation}
Its contribution on the \(\sqrt{T}\) scale is
\[
\frac{T}{S}
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\Delta}.
\]
Under Assumption~\ref{ass:smoothing}, both path-response functions converge to the same population function as \(b_T\rightarrow0\). With common rank paths,
\[
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{EE}},
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{SS}},
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{ES}},
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\mathrm{SE}}
\xrightarrow{\mathbb{P}}
\boldsymbol{\Omega}_{\mathrm{MC}},
\qquad
\boldsymbol{\hat{\Omega}}_{\mathrm{MC}}^{\Delta}
\xrightarrow{\mathbb{P}}
\boldsymbol{0}.
\]
Consequently, common random numbers remove first-order simulation noise from the smoothing comparison. With independent rank paths, the cross-covariances are zero and the simulation covariance of the difference converges to \(2\boldsymbol{\Omega}_{\mathrm{MC}}\).
For response \(m\), a useful numerical diagnostic is
\[
\hat{\pi}_{m,\mathrm{MC}}
=
\frac{
\hat{\Omega}_{\mathrm{MC},mm}^{\mathrm{E}}/S
}{
\hat{\Omega}_{mm}/T
+
\hat{\Omega}_{\mathrm{MC},mm}^{\mathrm{E}}/S
}.
\]
This quantity reports the fraction of the estimated variance due to numerical integration. In practice, \(S\) can be increased until this fraction is small for the reported responses and the estimates and standard errors are stable. The sufficient asymptotic condition is \(S/T\rightarrow\infty\), but no fixed path count or fixed ratio applies uniformly because simulation variance depends on the horizon, persistence, nonlinearity, initial state, and shock size.
Theorem~\ref{thm:bootstrap-validity} conditions on fixed response-simulation ranks and imposes \(S/T\rightarrow\infty\). When \(S/T\rightarrow\kappa\in\left(0,\infty\right)\), that bootstrap continues to reproduce the sampling component \(\boldsymbol{\Omega}\), while equation~\eqref{eq:finite-simulation-correction} supplies the additional simulation component.
Appendix~\ref{app:implementation} summarizes implementation and reporting conventions supported by the preceding results.
\section{Conclusion}
\label{sec:conclusion}
This paper studies inference for a recursively defined nonlinear impulse response when the structural innovation distributions are estimated from generated residuals. The population response is fixed by a normal-rank impact shock, a common future innovation law, and paired shocked and unshocked paths. The empirical-quantile estimator then changes only the estimate of the innovation quantile functions. This common target is necessary for interpreting differences between empirical and smoothed procedures as estimation effects.
First, the empirical-residual estimator is \(\sqrt{T}\)-asymptotically linear at fixed horizons. Its observation-level influence contribution separates direct transition estimation, the effect of transition estimation on residual order statistics, estimation of ordinary innovation quantiles, and estimation of the shifted impact quantile. After the quantile terms are projected through the recursive response, they can be evaluated from residual ranks and spacings. Feasible inference therefore does not require estimating the innovation density or selecting a smoothing bandwidth.
Next, the same-target smoothing comparison shows that the relevant bias is the quantile-smoothing error after it has been weighted by recursive propagation. Smoothing is first-order irrelevant only when that propagated bias is smaller than \(T^{-1/2}\). At the boundary it shifts the Gaussian limit, and above the boundary it dominates sampling uncertainty. The result distinguishes a local approximation property of the quantile estimator from its effect on the reported nonlinear response.
Finally, valid residual-bootstrap inference requires regeneration and re-estimation of the complete recursive model. Resampling response paths while holding the fitted transition and residual distribution fixed measures numerical integration error, not sampling uncertainty. When the number of simulated paths is proportional to the sample size, numerical integration contributes a separate covariance component; when it grows faster, this component vanishes at the \(\sqrt{T}\) scale. The analysis is deliberately fixed-horizon and maintains independent structural innovation components. Growing horizons, dependent components, and high-dimensional transition estimates require additional arguments beyond the present results.