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.
98,529 characters
Two-Point Dependent Wild Bootstrap for Weakly Dependent Estimating Equations
\maketitle
\begin{abstract}
\noindent
This paper develops a general two-point dependent wild bootstrap (DWB) for weakly dependent estimating equations.
Its key feature is that the two-point marginal distribution and the latent serial dependence specification can be chosen separately.
The construction combines a normalized two-point distribution with a stationary latent Gaussian process via a Gaussian copula transformation, includes dependent Rademacher and Mammen multipliers, and nests the classical iid two-point wild bootstrap as the serially independent case.
The induced multiplier autocovariances determine the lag weights in a corresponding heteroskedasticity- and autocorrelation-consistent (HAC) covariance estimator, which coincides exactly with the conditional covariance of the bootstrap estimating-equation sum.
We establish first-order bootstrap validity for asymptotically linear estimators by showing that the matched-HAC estimator consistently estimates the long-run covariance and that the bootstrap estimating-equation sum converges conditionally to the same Gaussian limit as its original-sample counterpart, yielding valid HAC-studentized $z$-tests and the corresponding Wald and Lagrange multiplier tests.
Monte Carlo experiments in nonlinear generalized method of moments (GMM) and linear regression show that Rademacher DWB generally provides more accurate finite-sample size control for $z$-tests than the Mammen and Gaussian DWB.
A GMM application to a nonlinear short-rate mean-reversion model illustrates the practical relevance of the proposed two-point DWB.
\end{abstract}
\noindent\textbf{Keywords:} dependent multipliers; two-point wild bootstrap; Rademacher multipliers; Mammen multipliers; dependent wild bootstrap; HAC covariance estimation; estimating equations. \\
\textbf{JEL codes:} C12, C14, C22.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Introduction}
Wild bootstrap methods provide a convenient approach to robust inference under heteroskedasticity when the disturbance distribution is otherwise left unspecified \citep{Wu1986,Liu1988,Mammen1993}.
Bounded two-point auxiliary distributions have played an important role in this literature.
Mammen's two-point distribution reproduces the first three moments relevant for classical higher-order arguments in heteroskedastic regression, while \citet{DavidsonMonticiniPeel2007} study a general class of mean-zero, unit-variance two-point distributions and the trade-off between their higher moments.
The symmetric Rademacher distribution provides the particularly simple $\{-1,1\}$ case and has attractive finite-sample properties in conventional wild bootstrap testing \citep{davidsonflachaire2008}.
Thus, even after imposing mean zero and unit variance, the choice of multiplier distribution leaves scope to vary higher-order marginal features such as skewness and kurtosis.
For dependent data, independent wild bootstrap weights do not reproduce the serial covariance relevant for inference.
The dependent wild bootstrap (DWB) addresses this problem by introducing serial dependence into the auxiliary variables.
A particularly convenient way of generating such dependence is through Gaussian multipliers, but this also fixes the marginal distribution of the multiplier: standardized Gaussian multipliers are unbounded and have fixed higher-order moments.
This removes the freedom to retain bounded two-point marginal distributions such as the Rademacher and Mammen laws.
This paper addresses this problem by showing that the two-point marginal law and the latent serial dependence specification can be chosen separately.
We develop a two-point DWB using a latent stationary Gaussian process and a Gaussian copula threshold transformation.
For a normalized two-point distribution with distribution function $F$, we define
\begin{align}
\xi_t
=
F^{-1}\{\Phi(Z_t)\},
\label{eq:intro_copula_transform}
\end{align}
where $(Z_t)$ is a stationary Gaussian process.
The probability mechanism underlying Gaussian thresholding is classical \citep{EmrichPiedmonte1991}; we use it here to construct a dependent multiplier process while preserving exactly the chosen two-point marginal law.
The marginal distribution of $\xi_t$ is exactly $F$, while its serial dependence is induced by the correlation structure of the latent Gaussian process through the threshold transformation.
Thus the chosen two-point marginal law can be preserved while the latent dependence specification is varied.
The same device accommodates Rademacher, Mammen, and other normalized bounded two-point multipliers, within the same dependence construction.
When the latent Gaussian process is serially independent, the construction reduces exactly to the classical iid two-point wild bootstrap, so the conventional iid Rademacher and Mammen wild bootstraps are nested as special cases.
The construction is related to, but distinct from, existing DWB and dependent multiplier procedures.
The DWB of \citet{Shao2010} introduces serial dependence into the auxiliary variables, and related dependent multiplier procedures have subsequently been developed for other statistics and empirical processes \citep{LeuchtNeumann2013,doukhannemann2015,buecherkojadinovic2016}.
Recent work by \citet{HounyoLin2026} develops multiway wild-cluster bootstrap procedures for linear regression with serially correlated common time effects, including a dependent Rademacher construction tailored to the time dimension.
Their objective is to reproduce the dependence induced by multiway clustering and serially correlated time effects, whereas our construction provides a general device for weakly dependent estimating equations that accommodates arbitrary normalized two-point laws and allows the marginal law and serial dependence profile to be specified separately.
To our knowledge, this is the first general two-point DWB that accommodates an arbitrary normalized two-point marginal law while allowing its serial dependence profile to be specified separately.
The construction also has a natural connection to heteroskedasticity-and-autocorrelation-consistent (HAC) covariance estimation.
The latent Gaussian correlation determines the covariance of the transformed two-point multipliers, which in turn determines the corresponding HAC lag weights.
Using these multiplier autocovariances as lag weights yields a matched-HAC estimator whose finite-sample value coincides exactly with the conditional covariance of the bootstrap score sum. For consistency, we separate the standard population-score HAC approximation from the additional effects of estimating and recentering the score contributions, and show that the latter are asymptotically negligible.
The general connection between DWB and HAC covariance estimation is already present in \citet{Shao2010} and is made explicit by \citet{davidsonmonticini2014}; in the present framework, however, the multiplier process, its autocovariance sequence, and the associated HAC lag weights arise jointly from the specified two-point marginal law and latent dependence structure.
We then establish first-order validity of the resulting two-point DWB in a general weakly dependent estimating-equation framework.
In particular, we show that the bootstrap estimating-equation sum converges conditionally to the same Gaussian limit as its original-sample counterpart, yielding valid bootstrap distributions for studentized statistics and the corresponding Wald and Lagrange multiplier statistics based on asymptotically linear estimators.
The proof requires an additional step beyond existing DWB results because the Gaussian copula dependence profile is not finite-range.
We construct a finite-dependent approximation that preserves the exact two-point marginal distribution, control the resulting approximation error, and establish the required conditional Gaussian limit using a blocking argument.
For linear regression, the validity results also extend to the null-imposed residual bootstrap under primitive conditions.
The Monte Carlo analysis considers nonlinear generalized method of moments (GMM) and linear regression under serial dependence, deterministic heteroskedasticity, and non-Gaussian innovations.
Across the designs considered, Rademacher DWB generally yields more accurate finite-sample size control for $z$-tests than Mammen DWB and Gaussian DWB.
Size distortions are also generally smaller for the restricted $z$-statistic $\tilde z_{LM}$ than for the unrestricted $z$-statistic $\hat z_W$, with particularly accurate size control obtained by combining $\tilde z_{LM}$ with Rademacher DWB.
The simulations further show that these differences across bootstrap procedures cannot be explained primarily by differences among their HAC kernels, since the corresponding asymptotic HAC procedures behave very similarly.
The corresponding two-step GMM results give qualitatively similar conclusions.
We illustrate the method using a nonlinear GMM specification based on the conditional-mean implication of the Cox--Ingersoll--Ross short-rate model.
The application shows that asymptotic and bootstrap approximations can lead to different conclusions at conventional significance levels.
The rest of the paper is organized as follows.
Section~\ref{sec:setup} introduces the estimating-equation framework and motivating examples.
Section~\ref{sec:constructing_multipliers} develops the dependent two-point multiplier construction, its matched-HAC representation, and its Rademacher and Mammen special cases.
Section~\ref{sec:estimating_equations} develops the bootstrap inference procedures for asymptotically linear estimators, including unrestricted and restricted studentized procedures and the regression-specific residual bootstrap.
Section~\ref{sec:asymptotics} establishes the asymptotic validity of the proposed procedures.
Section~\ref{sec:mc} presents the Monte Carlo experiments.
Section~\ref{sec:empirical} gives the empirical illustration.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Setup and First-Order Framework}\label{sec:setup}
To fix ideas, suppose that an estimator $\hat{\boldsymbol{\theta}}$ for a parameter of interest $\boldsymbol{\theta}_0\in\Theta\subset\mathbb{R}^k$ is asymptotically linear:
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
=
\mathbf{B}_0\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)+o_p(1),
\label{eq:asym_linear_theta}
\end{align}
where $\mathbf{g}_t(\boldsymbol{\theta})$ is a $d$-dimensional score, moment, or estimating-equation contribution and $\mathbf{B}_0$ is a $k\times d$ population linearization matrix.
The true value satisfies $\E\mathbf{g}_t(\boldsymbol{\theta}_0)=\mathbf{0}$.
We assume that $\{\mathbf{g}_t(\boldsymbol{\theta}_0)\}$ is covariance-stationary, weakly serially dependent, and satisfies
\begin{align}
\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\boldsymbol{\Omega}_0),
\label{eq:score_clt}
\end{align}
where $\boldsymbol{\Omega}_0
=
\sum_{h=-\infty}^{\infty}\boldsymbol{\Gamma}_0(h)>0$, $
\boldsymbol{\Gamma}_0(h)=\E\{\mathbf{g}_t(\boldsymbol{\theta}_0)\mathbf{g}_{t-h}(\boldsymbol{\theta}_0)'\}$.
Then
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{V}_{\theta}),
\label{eq:theta_limit_generic}
\end{align}
where $\mathbf{V}_{\theta}
:=
\mathbf{B}_0\boldsymbol{\Omega}_0\mathbf{B}_0'$.
Under serial dependence, feasible inference therefore requires estimation of the long-run covariance matrix $\boldsymbol{\Omega}_0$.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{HAC inference under serial dependence}\label{subsec:hac_inference}
Let $\hat{\mathbf{g}}_t=\mathbf{g}_t(\hat{\boldsymbol{\theta}})$ and $\overline{\hat{\mathbf{g}}}=T^{-1}\sum_{t=1}^T\hat{\mathbf{g}}_t$.
For $h\geq0$, define the sample autocovariance matrix
\begin{align}
\hat{\boldsymbol{\Gamma}}(h)
&=
\frac{1}{T}\sum_{t=h+1}^T
(\hat{\mathbf{g}}_t-\overline{\hat{\mathbf{g}}})
(\hat{\mathbf{g}}_{t-h}-\overline{\hat{\mathbf{g}}})'.
\label{eq:generic_sample_autocovariance}
\end{align}
Let $w_T(h)$ denote symmetric lag weights satisfying $w_T(0)=1$.
A generic HAC estimator of $\boldsymbol{\Omega}_0$ is
\begin{align}
\hat{\boldsymbol{\Omega}}_{\mathrm{HAC}}
&=
\hat{\boldsymbol{\Gamma}}(0)
+
\sum_{h=1}^{T-1}
w_T(h)
\left\{
\hat{\boldsymbol{\Gamma}}(h)+\hat{\boldsymbol{\Gamma}}(h)'
\right\}.
\label{eq:generic_hac}
\end{align}
For example, the Bartlett-type \citet{NeweyWest1987} estimator uses $w_T(h)=(1-|h|/b_T)\mathbf{1}\{|h|\leq b_T\}$, where $b_T$ is the bandwidth or truncation parameter.
Let $\hat{\mathbf{B}}$ be a consistent sample analogue of $\mathbf{B}_0$ and define
\begin{align}
\hat{\mathbf{V}}_{\theta,\mathrm{HAC}}
&=
\hat{\mathbf{B}}\hat{\boldsymbol{\Omega}}_{\mathrm{HAC}}\hat{\mathbf{B}}'.
\label{eq:generic_hac_vtheta}
\end{align}
When $\hat{\boldsymbol{\Omega}}_{\mathrm{HAC}}\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0$ and $\hat{\mathbf{B}}\stackrel{p}{\longrightarrow}\mathbf{B}_0$, replacing $\mathbf{V}_{\theta}$ by $\hat{\mathbf{V}}_{\theta,\mathrm{HAC}}$ gives conventional first-order HAC inference for $\boldsymbol{\theta}_0$.
Such asymptotic inference can nevertheless be sensitive in finite samples to long-run variance estimation and to the normal or chi-square approximation.
This motivates a dependent wild bootstrap that reproduces the serial covariance relevant for HAC inference while retaining bounded two-point multipliers.
The construction below chooses the HAC lag weights from the autocovariances of the bounded dependent multipliers and reproduces exactly the same weights in the bootstrap.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Examples}\label{subsec:examples}
\paragraph{Linear regression.}
Consider
\begin{align}
y_t=\mathbf{x}_t'\boldsymbol{\theta}_0+u_t,
\qquad
\E(\mathbf{x}_tu_t)=\mathbf{0},
\end{align}
where $\mathbf{x}_t$ is a $k$-dimensional regressor and $\boldsymbol{\theta}_0$ is the parameter of interest.
Define
\begin{align}
\mathbf{g}_t(\boldsymbol{\theta})=\mathbf{x}_t(y_t-\mathbf{x}_t'\boldsymbol{\theta}).
\end{align}
If $\mathbf{Q}=\E(\mathbf{x}_t\mathbf{x}_t')$ is nonsingular, the ordinary least squares (OLS) estimator satisfies
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
=
\mathbf{Q}^{-1}\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{x}_tu_t+o_p(1).
\end{align}
Thus $\mathbf{B}_0=\mathbf{Q}^{-1}$ and $\mathbf{g}_t(\boldsymbol{\theta}_0)=\mathbf{x}_tu_t$.
The estimated contribution is $\hat{\mathbf{g}}_t=\mathbf{g}_t(\hat{\boldsymbol{\theta}})=\mathbf{x}_t\hat u_t$, $\hat u_t=y_t-\mathbf{x}_t'\hat{\boldsymbol{\theta}}$.
\paragraph{Maximum likelihood and smooth $M$-estimation.}
Let $\ell_t(\boldsymbol{\theta})$ be a log-likelihood, quasi-log-likelihood, or smooth criterion contribution and define the score
\begin{align}
\mathbf{g}_t(\boldsymbol{\theta})=\partial \ell_t(\boldsymbol{\theta})/\partial\boldsymbol{\theta}.
\end{align}
Suppose that $\hat{\boldsymbol{\theta}}$ solves $T^{-1}\sum_{t=1}^T\mathbf{g}_t(\hat{\boldsymbol{\theta}})=\mathbf{0}$.
Let $\mathbf{G}_0=\E\left\{\partial\mathbf{g}_t(\boldsymbol{\theta}_0) /\partial\boldsymbol{\theta}'\right\}$.
If $\mathbf{G}_0$ is nonsingular and the standard Taylor expansion is valid, then
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
=
-\mathbf{G}_0^{-1}\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)+o_p(1),
\end{align}
so $\mathbf{B}_0=-\mathbf{G}_0^{-1}$.
Under correct likelihood specification, $\mathbf{G}_0$ is the negative information matrix; under quasi-likelihood or general smooth $M$-estimation, it is the corresponding sensitivity matrix.
\paragraph{GMM.}
Let $\boldsymbol{\theta}_0$ be defined by $\E\mathbf{g}_t(\boldsymbol{\theta}_0)=\mathbf{0}$, where $\mathbf{g}_t(\boldsymbol{\theta})$ is a $d$-dimensional moment vector.
Consider a GMM estimator that minimizes
\begin{align}
\overline{\mathbf{g}}_T(\boldsymbol{\theta})'\mathbf{W}_T\overline{\mathbf{g}}_T(\boldsymbol{\theta}),
\qquad
\overline{\mathbf{g}}_T(\boldsymbol{\theta})=T^{-1}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}),
\end{align}
where $\mathbf{W}_T$ is a symmetric positive-definite weighting matrix satisfying $\mathbf{W}_T\stackrel{p}{\longrightarrow}\mathbf{W}_0$ for some symmetric positive-definite matrix $\mathbf{W}_0$.
Define
$\mathbf{D}_0
=
\E\left\{
\partial\mathbf{g}_t(\boldsymbol{\theta}_0)/\partial\boldsymbol{\theta}'
\right\}$.
If $\mathbf{D}_0'\mathbf{W}_0\mathbf{D}_0$ is nonsingular, then
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
=
-(\mathbf{D}_0'\mathbf{W}_0\mathbf{D}_0)^{-1}\mathbf{D}_0'\mathbf{W}_0
\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)
+o_p(1),
\end{align}
so
$\mathbf{B}_0
=
-(\mathbf{D}_0'\mathbf{W}_0\mathbf{D}_0)^{-1}\mathbf{D}_0'\mathbf{W}_0$.
\paragraph{Panel data regression.}
The framework also covers panel settings in which cross-sectional score or moment contributions are aggregated into a time-indexed estimating-equation sequence. Provided that the resulting sequence satisfies the first-order conditions above, the dependent two-point multiplier can be applied along the time dimension in the same way as for the preceding examples.
Panel-specific theory and methods for large panels are developed by Dai, Matsushita, and Yamagata~(\citeyear{DaiMatsushitaYamagata2026}).
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Constructing Dependent Wild Bootstrap Multipliers}\label{sec:constructing_multipliers}
This section constructs the multiplier process used for score and moment bootstrap inference and links its covariance to HAC estimation.
The construction first specifies a bounded two-point marginal distribution and then introduces serial dependence through a latent Gaussian copula while preserving the exact marginal law.
The resulting covariance kernel defines the original-sample HAC estimator and, by construction, the conditional covariance of the bootstrap score sum.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Bounded two-point wild bootstrap distributions}\label{sec:twopoint}
Fix $p\in(0,1)$. Let $\xi$ take two values $a_p<0<b_p$ with probabilities
\begin{align}
\Pro(\xi=a_p)=p,
\qquad
\Pro(\xi=b_p)=1-p.
\end{align}
Imposing $\E\xi=0$ and $\E\xi^2=1$ gives
\begin{align}
a_p=-\sqrt{\frac{1-p}{p}},
\qquad
b_p=\sqrt{\frac{p}{1-p}}.
\label{eq:two_point_values}
\end{align}
The corresponding distribution function is denoted by $F_p$.
The third and fourth moments are
\begin{align}
\E\xi^3
=
\frac{2p-1}{\sqrt{p(1-p)}},
\qquad
\E\xi^4
=
\frac{(1-p)^2}{p}+\frac{p^2}{1-p}.
\label{eq:two_point_moments}
\end{align}
Two choices are central.
\paragraph{Davidson--Flachaire/Rademacher two-point multiplier.}
For $p=1/2$, $a_p=-1$ and $b_p=1$. This is the Rademacher two-point
wild bootstrap multiplier. It is the symmetric bounded choice emphasized by
\citet{davidsonflachaire2008}. It satisfies $\E\xi=0$, $\E\xi^2=1$,
$\E\xi^3=0$, and $\E\xi^4=1$.
\paragraph{Mammen skewness-replicating two-point multiplier.}
Mammen's two-point distribution is obtained by choosing $p$ so that
$\E\xi^3=1$. This gives
\begin{align}
p_M=\frac{\sqrt{5}+1}{2\sqrt{5}},
\qquad
a_M=\frac{1-\sqrt{5}}{2},
\qquad
b_M=\frac{1+\sqrt{5}}{2}.
\label{eq:mammen_distribution}
\end{align}
It satisfies $\E\xi=0$, $\E\xi^2=\E\xi^3=1$, and $\E\xi^4=2$. Thus Mammen's distribution is the skewness-replicating bounded two-point choice.
The distinction between these two cases is useful. The Rademacher choice is symmetric and maximally simple. The Mammen choice preserves the classical third-moment matching condition. Both are bounded and both are special cases of the same two-point family.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Copula-based bounded dependent two-point multipliers}
\label{sec:construction}
Let $\{Z_t\}_{t\in\mathbb{Z}}$ be a stationary standard Gaussian sequence satisfying
\begin{align}
\E Z_t=0,
\qquad
\E Z_t^2=1,
\qquad
\Corr(Z_t,Z_{t-h})=\rho_h,
\qquad
h\in\mathbb{Z}.
\label{eq:latent_gaussian_general}
\end{align}
By stationarity, $\rho_{-h}=\rho_h$.
For a sample of size $T$, define the $T\times1$ latent Gaussian vector
\begin{align}
\mathbf{z}_T
=
(Z_1,\ldots,Z_T)',
\qquad
\mathbf{z}_T
\sim
N(\mathbf{0},\mathbf{R}_T),
\qquad
(\mathbf{R}_T)_{t,s}
=
\rho_{|t-s|}.
\end{align}
Apply the Gaussian copula transformation coordinatewise:
\begin{align}
U_t
=
\Phi(Z_t),
\qquad
\xi_t
=
F_p^{-1}(U_t),
\qquad
t=1,\ldots,T.
\end{align}
Equivalently, with $q_p=\Phi^{-1}(p)$,
\begin{align}
\xi_t
=
a_p\,\mathbf{1}\{Z_t\leq q_p\}
+
b_p\,\mathbf{1}\{Z_t>q_p\}.
\label{eq:threshold_twopoint}
\end{align}
Thus each $\xi_t$ has exactly the prescribed two-point marginal distribution $F_p$.
Writing
\begin{align}
\boldsymbol{\xi}_T
=
(\xi_1,\ldots,\xi_T)',
\end{align}
the transformation maps the Gaussian vector $\mathbf{z}_T$ into a $T$-dimensional random vector whose coordinates each take the two values $a_p$ and $b_p$, with dependence inherited from the joint Gaussian distribution of $\mathbf{z}_T$.
For any two coordinates separated by lag $h$, $(Z_t,Z_{t-h})$ is bivariate standard normal with correlation $\rho_h$.
Since $\xi_t
=
F_p^{-1}\{\Phi(Z_t)\}$ and
$\xi_{t-h}
=
F_p^{-1}\{\Phi(Z_{t-h})\}$,
their correlation is determined by $\rho_h$, given $p$.
We denote the resulting correlation transformation by $K_p$, so that $K_p(\rho_h)$ is the correlation between $\xi_t$ and $\xi_{t-h}$.
The following proposition gives the explicit form of this correlation transformation and establishes that the resulting multiplier autocovariance sequence is valid.
\begin{prop}[Copula-induced correlation map]
\label{prop:copula_kernel}
For the normalized two-point distribution $F_p$,
\begin{align}
K_p(\rho)
=
\frac{\Phi_2(q_p,q_p;\rho)-p^2}{p(1-p)},
\qquad
-1\leq\rho\leq1,
\label{eq:general_copula_kernel}
\end{align}
where $\Phi_2(\cdot,\cdot;\rho)$ is the bivariate standard normal distribution function with correlation $\rho$.
Consequently, since $\{\xi_t\}$ is stationary, its autocovariance sequence is given by
\begin{align}
\E(\xi_t\xi_{t-h})
=
K_p(\rho_h),
\qquad
h\in\mathbb{Z}.
\end{align}
Moreover,
$\mathbf{K}_{p,T}
:=
\E(\boldsymbol{\xi}_T\boldsymbol{\xi}_T')$
is positive semi-definite for every finite $T$.
\end{prop}
\begin{proof}
Let $I_t=\mathbf{1}\{Z_t\leq q_p\}$.
Using $\xi_t=a_pI_t+b_p(1-I_t)$, $pa_p+(1-p)b_p=0$, and the definitions of $a_p$ and $b_p$,
\begin{align}
\xi_t
=
\frac{p-I_t}{\sqrt{p(1-p)}}.
\end{align}
Therefore
\begin{align}
\E(\xi_t\xi_{t-h})
&=
\frac{\E\{(p-I_t)(p-I_{t-h})\}}{p(1-p)}
\nonumber\\
&=
\frac{\E(I_tI_{t-h})-p^2}{p(1-p)}.
\end{align}
Since
$\E(I_tI_{t-h})
=
\Phi_2(q_p,q_p;\rho_h)$,
the expression for $K_p(\rho_h)$ follows.
Finally, for any $\mathbf{a}=(a_1,\ldots,a_T)'\in\mathbb{R}^T$,
$\mathbf{a}'\mathbf{K}_{p,T}\mathbf{a}
=
\E[
(
\sum_{t=1}^T a_t\xi_t
)^2
]
\geq0$; thus, $\mathbf{K}_{p,T}$ is positive semi-definite for every finite $T$.
\end{proof}
For implementation, we specify the symmetric autocorrelation sequence of the latent Gaussian process as
\begin{align}
\rho_{T,h}
&=
\exp\left\{
-\left(\frac{h}{b_T}\right)^2
\right\},
\qquad
h\in\mathbb{Z},
\label{eq:latent_gaussian_corr}
\end{align}
where $b_T>0$ is the dependence-scale parameter.
For $h\geq0$, Proposition~\ref{prop:copula_kernel} gives
\begin{align}
\E(\xi_t\xi_{t-h})
=
K_p(\rho_{T,h})
=
K_p\left[
\exp\left\{
-\left(\frac{h}{b_T}\right)^2
\right\}
\right].
\label{eq:induced_kernel_general}
\end{align}
For notational simplicity, the dependence of $Z_t$ and $\xi_t$ on $T$ is suppressed when \eqref{eq:latent_gaussian_corr} is used.
The induced multiplier autocovariance $K_p(\rho_{T,h})$ will also serve as the HAC lag weight at lag $h$.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{From multiplier dependence to matched-HAC inference}\label{sec:matched}
\paragraph{Original-sample HAC.}
In the generic HAC estimator \eqref{eq:generic_hac}, choose
\begin{align}
w_T(h)
&=
K_p(\rho_{T,h}).
\label{eq:matched_hac_weight}
\end{align}
The resulting original-sample HAC estimator is
\begin{align}
\hat{\boldsymbol{\Omega}}_p
&:=
\hat{\boldsymbol{\Gamma}}(0)
+
\sum_{h=1}^{T-1}
K_p(\rho_{T,h})
\left\{
\hat{\boldsymbol{\Gamma}}(h)+\hat{\boldsymbol{\Gamma}}(h)'
\right\}.
\label{eq:matched_hac}
\end{align}
Equivalently,
\begin{align}
\hat{\boldsymbol{\Omega}}_p
&=
\frac{1}{T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
(\hat{\mathbf{g}}_t-\overline{\hat{\mathbf{g}}})
(\hat{\mathbf{g}}_s-\overline{\hat{\mathbf{g}}})'.
\label{eq:double_sum_hac}
\end{align}
The corresponding original-sample covariance estimator for the linearized estimator is
\begin{align}
\hat{\mathbf{V}}_{\theta}
&=
\hat{\mathbf{B}}\hat{\boldsymbol{\Omega}}_p\hat{\mathbf{B}}'.
\label{eq:vtheta_hat}
\end{align}
Thus $\hat{\boldsymbol{\Omega}}_p$ is the matched-HAC estimator targeting the
long-run covariance $\boldsymbol{\Omega}_0$, with lag weight at $h$ given by the
multiplier autocovariance $K_p(\rho_{T,h})$.
Section~\ref{sec:asymptotics} establishes conditions under which
$\hat{\boldsymbol{\Omega}}_p\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0$.
\paragraph{Bootstrap match.}
For each bootstrap replication, generate $\{\xi_t^*\}_{t=1}^T$ independently of the data from the same Gaussian copula threshold construction, so that
\begin{align}
\E^*(\xi_t^*)
&=
0,
\qquad
\E^*(\xi_t^*\xi_{t-h}^*)
=
K_p(\rho_{T,h}),
\label{eq:bootstrap_multiplier_covariance}
\end{align}
and define
\begin{align}
\hat{\mathbf{g}}_t^*
&=
\xi_t^*(\hat{\mathbf{g}}_t-\overline{\hat{\mathbf{g}}}).
\label{eq:score_bootstrap_theta}
\end{align}
Conditional on the data,
\begin{align}
\E^*(\hat{\mathbf{g}}_t^*\hat{\mathbf{g}}_{t-h}^{*'})
&=
K_p(\rho_{T,h})
(\hat{\mathbf{g}}_t-\overline{\hat{\mathbf{g}}})
(\hat{\mathbf{g}}_{t-h}-\overline{\hat{\mathbf{g}}})'.
\label{eq:bootstrap_score_pair_covariance}
\end{align}
Consequently, the conditional covariance of the bootstrap score sum exactly matches the original-sample HAC estimator:
\begin{align}
\Var^*\left(
\frac{1}{\sqrt{T}}\sum_{t=1}^T\hat{\mathbf{g}}_t^*
\right)
&=
\frac{1}{T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
(\hat{\mathbf{g}}_t-\overline{\hat{\mathbf{g}}})
(\hat{\mathbf{g}}_s-\overline{\hat{\mathbf{g}}})'\\
&=
\hat{\boldsymbol{\Omega}}_p.
\label{eq:exact_hac_match}
\end{align}
Hence the original-sample HAC estimator and the conditional covariance of the bootstrap score sum coincide exactly.
This exact equality is the matching property: the HAC weight $K_p(\rho_{T,h})$ at lag $h$ is also the covariance of two bootstrap multipliers $h$ periods apart.
Consequently, once matched-HAC consistency is established, the conditional covariance of the bootstrap score sum also converges to $\boldsymbol{\Omega}_0$.
The asymptotic results below further show that the bootstrap score sum converges conditionally to the same $N(\mathbf{0},\boldsymbol{\Omega}_0)$ limit as its original-sample counterpart.
For draw-specific studentization, define the bootstrap analogue of $\hat{\boldsymbol{\Omega}}_p$ from the bootstrap score array:
\begin{align}
\hat{\boldsymbol{\Omega}}_p^*
&=
\hat{\boldsymbol{\Gamma}}^*(0)
+
\sum_{h=1}^{T-1}
K_p(\rho_{T,h})
\left\{
\hat{\boldsymbol{\Gamma}}^*(h)+\hat{\boldsymbol{\Gamma}}^*(h)'
\right\},
\label{eq:bootstrap_matched_hac}
\end{align}
where, for $h\geq0$,
$\hat{\boldsymbol{\Gamma}}^*(h)
=
\frac{1}{T}\sum_{t=h+1}^T
(\hat{\mathbf{g}}_t^*-\overline{\hat{\mathbf{g}}^*})
(\hat{\mathbf{g}}_{t-h}^*-\overline{\hat{\mathbf{g}}^*})'$, $
\overline{\hat{\mathbf{g}}^*}
=
\frac{1}{T}\sum_{t=1}^T\hat{\mathbf{g}}_t^*$.
The associated draw-specific covariance estimator of $\hat{\boldsymbol{\theta}}^*$ is
\begin{align}
\hat{\mathbf{V}}_{\theta}^*
&=
\hat{\mathbf{B}}\hat{\boldsymbol{\Omega}}_p^*\hat{\mathbf{B}}'.
\label{eq:bootstrap_theta_covariance}
\end{align}
Note that $\hat{\boldsymbol{\Omega}}_p^*$ and $\hat{\mathbf{V}}_{\theta}^*$ vary across bootstrap replications.
\begin{rem}[Positive semi-definiteness]
Equation \eqref{eq:exact_hac_match} shows that $\hat{\boldsymbol{\Omega}}_p$ is a conditional covariance matrix and is therefore positive semi-definite for every finite sample.
At the same time, \eqref{eq:matched_hac} identifies it as a HAC estimator with lag weights $K_p(\rho_{T,h})$.
\end{rem}
\subsubsection{The Rademacher special case}
When $p=1/2$, $q_p=0$, $a_p=-1$, and $b_p=1$.
Then $\xi_t=\sgn(Z_t)$ up to the value assigned at zero.
The Gaussian sign-correlation identity gives
\begin{align}
K_{1/2}(\rho)=\frac{2}{\pi}\arcsin(\rho).
\label{eq:arcsine_special_case}
\end{align}
Writing $x=h/b_T$, the corresponding lag-weight function is
\begin{align}
K_{1/2}\{\exp(- x^2)\}
=
\frac{2}{\pi}\arcsin\{\exp(- x^2)\}.
\label{eq:arcsine_kernel}
\end{align}
Near the origin,
\begin{align}
\frac{2}{\pi}\arcsin\{\exp(- x^2)\}
=
1-\frac{2\sqrt{2}}{\pi}|x|+o(|x|),
\qquad x\to0.
\label{eq:kernel_cusp}
\end{align}
Thus the Rademacher case has an explicit arcsine lag-weight function with a Bartlett-like cusp.
\subsubsection{The Mammen special case}
When $p=p_M$ in \eqref{eq:mammen_distribution}, the multiplier process has exact Mammen marginals at every date.
The covariance transformation is
\begin{align}
K_{p_M}(\rho)=\frac{\Phi_2(q_{p_M},q_{p_M};\rho)-p_M^2}{p_M(1-p_M)}.
\label{eq:mammen_kernel}
\end{align}
There is no arcsine simplification in this asymmetric case, but the transformation is explicit and can be evaluated numerically. For a specified multiplier correlation at a single lag, the corresponding latent Gaussian correlation can be obtained numerically by inverting $K_{p_M}$. In our construction, however, the latent Gaussian correlation sequence is specified directly through \eqref{eq:latent_gaussian_corr}, so no pointwise inversion of a prespecified sequence of multiplier correlations across lags is required.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Dependent Wild Bootstrap for Estimating Equations}\label{sec:estimating_equations}
The setup in Section~\ref{sec:setup} gives
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
&=
\mathbf{B}_0\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)+o_p(1),
\end{align}
and $\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)\stackrel{d}{\longrightarrow} N(\mathbf{0},\mathbf{V}_{\theta})$ with $\mathbf{V}_{\theta}=\mathbf{B}_0\boldsymbol{\Omega}_0\mathbf{B}_0'$.
The dependent wild bootstrap applies the bootstrap contributions $\{\hat{\mathbf{g}}_t^*\}$ defined in \eqref{eq:score_bootstrap_theta} to this first-order representation.
The bootstrap analogue of the asymptotic linear representation is
\begin{align}
\sqrt{T}(\hat{\boldsymbol{\theta}}^*-\hat{\boldsymbol{\theta}})
&=
\hat{\mathbf{B}}
\frac{1}{\sqrt{T}}\sum_{t=1}^T
\hat{\mathbf{g}}_t^*,
\label{eq:bootstrap_theta_update}
\end{align}
where $\hat{\mathbf{B}}$ is the sample analogue of $\mathbf{B}_0$.
The notation $\hat{\boldsymbol{\theta}}^*$ is symbolic in this general formulation: the procedure resamples the first-order representation directly and does not require re-estimation from bootstrap data.
In the linear regression implementation considered in Subsection~\ref{subsec:reg_residual}, where the observed regressors are held fixed across bootstrap draws, this direct bootstrap update coincides exactly with the unrestricted OLS estimator computed from the generated bootstrap sample.
For likelihood scores or just-identified smooth estimating equations, $\hat{\mathbf{B}}=-\hat{\mathbf{G}}^{-1}$ with $\hat{\mathbf{G}}=T^{-1}\sum_{t=1}^T\partial\mathbf{g}_t(\hat{\boldsymbol{\theta}})/\partial\boldsymbol{\theta}'$.
For efficient two-step GMM,
$\hat{\mathbf{B}}=-(\hat{\mathbf{D}}'\mathbf{W}_T\hat{\mathbf{D}})^{-1}
\hat{\mathbf{D}}'\mathbf{W}_T$, where $\hat{\mathbf{D}}=T^{-1}\sum_{t=1}^T\partial\mathbf{g}_t(\hat{\boldsymbol{\theta}})/\partial\boldsymbol{\theta}'$.
The weighting matrix $\mathbf{W}_T$ is constructed from a first-step GMM estimator and satisfies $\mathbf{W}_T\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0^{-1}$.
When a matched-HAC estimator is used in this construction,
$\mathbf{W}_T-\hat{\boldsymbol{\Omega}}_p^{-1}=o_p(1)$.
The matched original-sample covariance estimator $\hat{\mathbf{V}}_{\theta}$ is given in \eqref{eq:vtheta_hat}, and \eqref{eq:exact_hac_match} shows that its moment covariance component is reproduced exactly by the bootstrap.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Unrestricted and restricted studentized implementations}\label{subsec:wald_lm_implementations}
We consider a fixed number $\ell$ of linear restrictions, $H_0:\mathbf{A}\boldsymbol{\theta}=\mathbf{a}_0$, where $\mathbf{A}\in\mathbb{R}^{\ell\times k}$ has rank $\ell$ and $\mathbf{a}_0\in\mathbb{R}^\ell$.
It is useful to formulate inference first in terms of feasible studentized $\ell$-dimensional vectors and then obtain the conventional Wald and LM statistics by taking quadratic forms.
The unrestricted Wald implementation uses the unrestricted estimator, whereas the restricted LM implementation imposes the null before forming the bootstrap moments.
We write $\hat{\mathbf{z}}_{W}$ for the feasible unrestricted studentized vector and $\tilde{\mathbf{z}}_{LM}$ for the feasible restricted studentized vector.
For $\ell=1$, we refer to the scalar quantities $\hat z_W$ and $\tilde z_{LM}$ as the unrestricted and restricted $z$-statistics, respectively.
For scalar restrictions, retaining the signed studentized statistic is
particularly useful because it preserves the direction of departures from
the null, which can be economically meaningful in regression and structural
parameter applications. It also permits the lower and upper tails of the
sampling approximation to be assessed separately, whereas the corresponding
quadratic Wald and LM statistics fold the two directions together.
All inverse square roots below denote the symmetric positive-definite inverse square root of the corresponding covariance matrix.
For the unrestricted implementation, define the feasible studentized restriction vector
\begin{align}
\hat{\mathbf{z}}_{W}
&=
(\mathbf{A}\hat{\mathbf{V}}_{\theta}\mathbf{A}')^{-1/2}
\sqrt{T}(\mathbf{A}\hat{\boldsymbol{\theta}}-\mathbf{a}_0),
\qquad
W_T
=
\hat{\mathbf{z}}_{W}'\hat{\mathbf{z}}_{W}.
\label{eq:wald_stat_generic}
\end{align}
Thus the usual Wald statistic is the squared Euclidean norm of the studentized restriction vector.
To approximate the distribution of $\hat{\mathbf{z}}_{W}$, use the bootstrap analogue in \eqref{eq:bootstrap_theta_update} and, for recomputed studentization, define
\begin{align}
\hat{\mathbf{z}}_{W}^*
&=
(\mathbf{A}\hat{\mathbf{V}}_{\theta}^*\mathbf{A}')^{-1/2}
\sqrt{T}\mathbf{A}(\hat{\boldsymbol{\theta}}^*-\hat{\boldsymbol{\theta}}),
\qquad
W_T^*
=
\hat{\mathbf{z}}_{W}^{*'}\hat{\mathbf{z}}_{W}^*.
\label{eq:wald_bootstrap_generic}
\end{align}
For fixed studentization, replace $\hat{\mathbf{V}}_{\theta}^*$ in \eqref{eq:wald_bootstrap_generic} by the original-sample matrix $\hat{\mathbf{V}}_{\theta}$.
The restricted implementation imposes the null before the moment contribution is formed.
Since $\operatorname{rank}(\mathbf{A})=\ell$, use an invertible linear reparameterization $\boldsymbol{\eta}=(\boldsymbol{\eta}_1',\boldsymbol{\eta}_2')'$ such that $\boldsymbol{\eta}_2=\mathbf{A}\boldsymbol{\theta}$, where $\boldsymbol{\eta}_2\in\mathbb{R}^\ell$.
Under the null, $\boldsymbol{\eta}_2=\mathbf{a}_0$, while $\boldsymbol{\eta}_1$ contains the $k-\ell$ nuisance parameters.
Write $\boldsymbol{\eta}_{1,0}$ for the true nuisance parameter, $\tilde{\boldsymbol{\eta}}_1$ for its restricted estimate, and $\tilde{\boldsymbol{\theta}}$ for the corresponding restricted estimate in the original parameterization, and define $\tilde{\mathbf{g}}_t=\mathbf{g}_t(\tilde{\boldsymbol{\theta}})$ and $\overline{\tilde{\mathbf{g}}}=T^{-1}\sum_{t=1}^T\tilde{\mathbf{g}}_t$.
The matched-HAC estimator based on the restricted moments is
\begin{align}
\tilde{\boldsymbol{\Omega}}_p
&=
\frac{1}{T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
(\tilde{\mathbf{g}}_t-\overline{\tilde{\mathbf{g}}})
(\tilde{\mathbf{g}}_s-\overline{\tilde{\mathbf{g}}})'.
\label{eq:restricted_matched_hac}
\end{align}
Let $\tilde{\mathbf{G}}=(\tilde{\mathbf{G}}_1,\tilde{\mathbf{G}}_2)$ be the sample Jacobian with respect to $(\boldsymbol{\eta}_1',\boldsymbol{\eta}_2')'$ at $\tilde{\boldsymbol{\theta}}$, where $\tilde{\mathbf{G}}_1\in\mathbb{R}^{d\times(k-\ell)}$ and $\tilde{\mathbf{G}}_2\in\mathbb{R}^{d\times\ell}$.
Following \citet{NeweyMcFadden1994}, define
\begin{align}
\tilde{\mathbf{G}}_{2\cdot1}
&=
\tilde{\mathbf{G}}_2
-
\tilde{\mathbf{G}}_1
(\tilde{\mathbf{G}}_1'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_1)^{-1}
\tilde{\mathbf{G}}_1'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_2.
\label{eq:g2_partial_generic}
\end{align}
By construction, $\tilde{\mathbf{G}}_1'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_{2\cdot1}=\mathbf{0}$.
The corresponding $\ell$-dimensional restricted studentized vector is
\begin{align}
\tilde{\mathbf{z}}_{LM}
&=
-\left(
\tilde{\mathbf{G}}_{2\cdot1}'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\tilde{\mathbf{G}}_{2\cdot1}
\right)^{-1/2}
\tilde{\mathbf{G}}_{2\cdot1}'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\sqrt{T}\,\overline{\tilde{\mathbf{g}}},
\qquad
LM_T
=
\tilde{\mathbf{z}}_{LM}'\tilde{\mathbf{z}}_{LM}.
\label{eq:lm_stat_generic}
\end{align}
For a scalar restriction, the leading minus sign is chosen so that, to first order, the sign of the restricted $z$-statistic $\tilde z_{LM}$ agrees with that of the unrestricted $z$-statistic $\hat z_W$, while leaving the LM statistic unchanged.
Thus $LM_T$ is the usual LM statistic, given by the squared Euclidean norm of the restricted studentized vector.
For the direct bootstrap, define the restricted bootstrap moments and their sample average by
\begin{align}
\tilde{\mathbf{g}}_t^*
&=
\xi_t^*(\tilde{\mathbf{g}}_t-\overline{\tilde{\mathbf{g}}}),
\qquad
\overline{\tilde{\mathbf{g}}^*}
=
T^{-1}\sum_{t=1}^T\tilde{\mathbf{g}}_t^*.
\label{eq:restricted_moment_bootstrap}
\end{align}
Note that, conditional on the data,
\begin{align*}
\Var^*\left(
\frac{1}{\sqrt{T}}\sum_{t=1}^T\tilde{\mathbf{g}}_t^*
\right)
&=
\tilde{\boldsymbol{\Omega}}_p
\end{align*}
holds exactly.
For draw-specific studentization, define
\begin{align}
\tilde{\boldsymbol{\Omega}}_p^*
&=
\frac{1}{T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
(\tilde{\mathbf{g}}_t^*-\overline{\tilde{\mathbf{g}}^*})
(\tilde{\mathbf{g}}_s^*-\overline{\tilde{\mathbf{g}}^*})'.
\label{eq:restricted_bootstrap_matched_hac}
\end{align}
The original-sample projection $\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}$ is held fixed across bootstrap draws, and the recomputed feasible bootstrap restricted studentized vector is
\begin{align}
\tilde{\mathbf{z}}_{LM}^*
&=
-\left\{
\tilde{\mathbf{G}}_{2\cdot1}'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\tilde{\boldsymbol{\Omega}}_p^*
\tilde{\boldsymbol{\Omega}}_p^{-1}
\tilde{\mathbf{G}}_{2\cdot1}
\right\}^{-1/2}
\tilde{\mathbf{G}}_{2\cdot1}'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\sqrt{T}\,\overline{\tilde{\mathbf{g}}^*},
\qquad
LM_T^*
=
\tilde{\mathbf{z}}_{LM}^{*'}\tilde{\mathbf{z}}_{LM}^*.
\label{eq:lm_bootstrap_generic}
\end{align}
For fixed studentization, replace $\tilde{\boldsymbol{\Omega}}_p^*$ in the covariance matrix of \eqref{eq:lm_bootstrap_generic} by $\tilde{\boldsymbol{\Omega}}_p$, leaving the same original-sample projection fixed.
For restricted efficient two-step GMM, the role of the projection $\tilde{\mathbf{G}}_{2\cdot1}$ can also be seen from the first-order conditions.
The weighting matrix constructed from the first-step GMM estimator satisfies $\mathbf{W}_T-\tilde{\boldsymbol{\Omega}}_p^{-1}=o_p(1)$, so the first-order conditions imply $\tilde{\mathbf{G}}_1'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\overline{\tilde{\mathbf{g}}}
=
o_p(T^{-1/2})$. Hence, to first order,
$\tilde{\mathbf{G}}_{2\cdot1}'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\sqrt{T}\,\overline{\tilde{\mathbf{g}}}
=
\tilde{\mathbf{G}}_2'
\tilde{\boldsymbol{\Omega}}_p^{-1}
\sqrt{T}\,\overline{\tilde{\mathbf{g}}}$;
see \citet[][p.~2230]{NeweyMcFadden1994}.
However, the direct bootstrap does not re-estimate the restricted nuisance parameter in each draw, so the corresponding bootstrap first-order condition is not available.
The projection $\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}$ is therefore retained in \eqref{eq:lm_bootstrap_generic} to account for the first-order effect of restricted nuisance-parameter estimation.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Regression-specific residual resampling}\label{subsec:reg_residual}
Regression models permit an additional implementation in which the same fixed $\ell$-dimensional linear restriction is imposed and the multiplier is applied to residuals rather than directly to the score.
\subsubsection{Original-sample statistics}
Consider
\begin{align}
y_t
=
\mathbf{x}_t'\boldsymbol{\theta}+u_t,
\qquad
H_0:\mathbf{A}\boldsymbol{\theta}=\mathbf{a}_0,
\qquad
t=1,\ldots,T,
\label{eq:reg_model}
\end{align}
where $\mathbf{x}_t=(1,\mathbf{z}_t')'\in\mathbb{R}^k$.
We assume throughout the regression framework that the intercept is included under both the null and alternative and is left unrestricted by $\mathbf{A}$.
Let $\hat{\boldsymbol{\theta}}$ denote the unrestricted OLS estimator and define
\begin{align}
\hat{\mathbf{Q}}_x
&=
T^{-1}\sum_{t=1}^T\mathbf{x}_t\mathbf{x}_t',
\\
\hat u_t
&=
y_t-\mathbf{x}_t'\hat{\boldsymbol{\theta}},
\qquad
\hat{\mathbf{g}}_t
=
\mathbf{x}_t\hat u_t.
\end{align}
The unrestricted normal equations imply $\overline{\hat{\mathbf{g}}}=\mathbf{0}$.
Let $\hat{\boldsymbol{\Omega}}_p$ denote the matched-HAC estimator based on $\{\hat{\mathbf{g}}_t\}$ and set
\begin{align}
\hat{\mathbf{V}}_{\theta}
=
\hat{\mathbf{Q}}_x^{-1}\hat{\boldsymbol{\Omega}}_p\hat{\mathbf{Q}}_x^{-1}.
\label{def:vcov_wald_linear}
\end{align}
The regression Wald statistic is therefore \eqref{eq:wald_stat_generic} with $\hat{\mathbf{B}}=\hat{\mathbf{Q}}_x^{-1}$.
For the restricted implementation, choose a fixed nonsingular matrix $\mathbf{C}=(\mathbf{C}_1,\mathbf{C}_2)$ such that $\mathbf{A}\mathbf{C}_1=\mathbf{0}$ and $\mathbf{A}\mathbf{C}_2=\mathbf{I}_\ell$, and write $\boldsymbol{\theta}=\mathbf{C}_1\boldsymbol{\eta}_1+\mathbf{C}_2\boldsymbol{\eta}_2$ so that $\boldsymbol{\eta}_2=\mathbf{A}\boldsymbol{\theta}$.
Let $\tilde{\boldsymbol{\eta}}_1$ denote the restricted OLS estimator under $\boldsymbol{\eta}_2=\mathbf{a}_0$ and define
\begin{align}
\tilde{\boldsymbol{\theta}}
&=
\mathbf{C}_1\tilde{\boldsymbol{\eta}}_1+\mathbf{C}_2\mathbf{a}_0,
\\
\tilde u_t
&=
y_t-\mathbf{x}_t'\tilde{\boldsymbol{\theta}},
\qquad
\tilde{\mathbf{g}}_t
=
\mathbf{x}_t\tilde u_t,
\qquad
\overline{\tilde{\mathbf{g}}}
=
T^{-1}\sum_{t=1}^T\tilde{\mathbf{g}}_t.
\end{align}
Because the intercept is an unrestricted nuisance parameter, the restricted normal equations imply $\sum_{t=1}^T\tilde u_t=0$.
Thus the restricted residuals themselves require no additional recentering before multiplication, although the restricted moment array is centered in the direct bootstrap and matched-HAC studentizer.
For the matched-HAC implementation, let $\tilde{\boldsymbol{\Omega}}_p$ be defined from $\{\tilde{\mathbf{g}}_t-\overline{\tilde{\mathbf{g}}}\}$ as in \eqref{eq:restricted_matched_hac}.
For $\mathbf{g}_t(\boldsymbol{\theta})=\mathbf{x}_t(y_t-\mathbf{x}_t'\boldsymbol{\theta})$,
\begin{align}
\tilde{\mathbf{G}}
=
(\tilde{\mathbf{G}}_1,\tilde{\mathbf{G}}_2)
=
-\hat{\mathbf{Q}}_x(\mathbf{C}_1,\mathbf{C}_2).
\end{align}
Using these regression-specific derivative matrices, define $\tilde{\mathbf{G}}_{2\cdot1}$ by \eqref{eq:g2_partial_generic}.
The regression LM statistic is therefore \eqref{eq:lm_stat_generic} with these regression-specific objects.
\subsubsection{Residual resampling}
For a Wald-type residual bootstrap, set $u_t^\dagger=\xi_t^*\hat u_t$ and generate
\begin{align}
y_t^*
=
\mathbf{x}_t'\hat{\boldsymbol{\theta}}+u_t^\dagger.
\end{align}
Here $u_t^\dagger$ is the generated bootstrap disturbance.
Let $\hat{\boldsymbol{\theta}}^*$ denote the unrestricted OLS estimator from the bootstrap sample.
As shown below, this estimator coincides exactly with the direct-score bootstrap update denoted by $\hat{\boldsymbol{\theta}}^*$ in \eqref{eq:bootstrap_theta_update}, so we use the same notation.
Define the post-estimation bootstrap residual and residual-bootstrap moment contribution by
\begin{align}
\hat u_t^*
&=
y_t^*-\mathbf{x}_t'\hat{\boldsymbol{\theta}}^*,
\qquad
\hat{\mathbf{g}}_{t,\mathrm{res}}^*
=
\mathbf{x}_t\hat u_t^*.
\end{align}
The generated disturbance $u_t^\dagger$ is generally different from the post-estimation residual $\hat u_t^*$.
The OLS normal equations give exactly
\begin{align}
\sqrt T\mathbf{A}(\hat{\boldsymbol{\theta}}^*-\hat{\boldsymbol{\theta}})
=
\mathbf{A}\hat{\mathbf{Q}}_x^{-1}
T^{-1/2}\sum_{t=1}^T\xi_t^*\hat{\mathbf{g}}_t.
\label{eq:reg_resid_wald_numerator}
\end{align}
Thus the Wald numerator from refitting the bootstrap sample is exactly the same as that obtained by direct score resampling in \eqref{eq:bootstrap_theta_update}, since
$\overline{\hat{\mathbf{g}}}=\mathbf{0}$.
For fixed studentization, replace $\hat{\mathbf{V}}_{\theta}^*$ in \eqref{eq:wald_bootstrap_generic} by the original-sample matrix
$\hat{\mathbf{V}}_{\theta}$ defined in \eqref{def:vcov_wald_linear}; for recomputed studentization, construct $\hat{\mathbf{V}}_{\theta}^*$
from the draw-specific moment array $\{\hat{\mathbf{g}}_{t,\mathrm{res}}^*\}$.
For an LM-type residual bootstrap, impose the restriction before resampling by setting $\tilde u_t^\dagger=\xi_t^*\tilde u_t$ and generate
\begin{align}
y_t^*
=
\mathbf{x}_t'\tilde{\boldsymbol{\theta}}+\tilde u_t^\dagger.
\end{align}
For each bootstrap sample, re-estimate $\boldsymbol{\eta}_1$ under $\boldsymbol{\eta}_2=\mathbf{a}_0$, let $\tilde{\boldsymbol{\eta}}_1^*$ denote the resulting nuisance estimate, and define
\begin{align}
\tilde{\boldsymbol{\theta}}^*
&=
\mathbf{C}_1\tilde{\boldsymbol{\eta}}_1^*+\mathbf{C}_2\mathbf{a}_0,
\\
\tilde u_t^*
&=
y_t^*-\mathbf{x}_t'\tilde{\boldsymbol{\theta}}^*,
\qquad
\tilde{\mathbf{g}}_{t,\mathrm{res}}^*
=
\mathbf{x}_t\tilde u_t^*.
\end{align}
The definition of $\tilde{u}_t^*$ and $\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_1=\mathbf{0}$ give the exact decomposition\footnote{This decomposition does not require the regression to contain an intercept; the intercept condition above is used only to ensure that the restricted residuals have zero sample mean.}
\begin{align}
\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}
T^{-1/2}\sum_{t=1}^T\tilde{\mathbf{g}}_{t,\mathrm{res}}^*
=
\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}
T^{-1/2}\sum_{t=1}^T
\xi_t^*(\tilde{\mathbf{g}}_t-\overline{\tilde{\mathbf{g}}})
+
\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}
\overline{\tilde{\mathbf{g}}}\,
T^{-1/2}\sum_{t=1}^T\xi_t^* .
\label{eq:reg_resid_lm_numerator}
\end{align}
The first term on the right-hand side is the centered direct-bootstrap projected moment, while the second is a centering correction that is asymptotically negligible under $H_0$.
For recomputed studentization, construct $\tilde{\boldsymbol{\Omega}}_p^*$ from $\{\tilde{\mathbf{g}}_{t,\mathrm{res}}^*\}$ after subtracting its bootstrap
sample average, and use it as the middle covariance matrix of \eqref{eq:lm_bootstrap_generic}, while $\tilde{\mathbf{G}}_{2\cdot1}$ and $\tilde{\boldsymbol{\Omega}}_p^{-1}$ remain fixed at their original restricted-sample values.
Primitive conditions for residual resampling are given in Subsection~\ref{subsec:reg_validity}.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Asymptotic Results}\label{sec:asymptotics}
In this section, we derive the asymptotic properties of the dependent wild bootstrap.
The primary inferential results are Gaussian approximations for the studentized estimating-equation vector and the feasible unrestricted and restricted studentized vectors; the familiar chi-square limits for the Wald and LM statistics then follow by continuous mapping.
\begin{ass}
\begin{itemize}
\item[(i)] The asymptotic expansion \eqref{eq:asym_linear_theta} and the score central limit theorem (CLT) \eqref{eq:score_clt} hold, with $\boldsymbol{\Omega}_0$ positive definite.
\item[(ii)] $\sup_t\E\|\mathbf{g}_t(\boldsymbol{\theta}_0)\|^2<\infty$.
\item[(iii)] $\sup_t\E\sup_{\boldsymbol{\theta}\in\Theta}\|(\partial/\partial\boldsymbol{\theta}')\mathbf{g}_t(\boldsymbol{\theta})\|^2<\infty$.
\item[(iv)] $\boldsymbol{\theta}_0$ is an interior point of $\Theta$, and $\mathbf{g}_t(\boldsymbol{\theta})$ is continuously differentiable in a neighborhood of $\boldsymbol{\theta}_0$.
\end{itemize}
\label{ass:moment_weakdepend}
\end{ass}
Assumption~\ref{ass:moment_weakdepend}(i) restates the first-order framework in Section~\ref{sec:setup}.
Assumption~\ref{ass:moment_weakdepend}(iii) is a standard condition used to establish HAC consistency and asymptotic normality of smooth estimators \citep{Andrews1991,NeweyMcFadden1994}.
\begin{ass}[Bandwidth]
The bandwidth $b_T$ used in \eqref{eq:double_sum_hac} satisfies $b_T\to\infty$ and $b_T/T\to0$.
\label{ass:bandwidth}
\end{ass}
Using the same lag weights $K_p(\rho_{T,h})$ as in the feasible matched-HAC estimator,
define the corresponding infeasible finite-sample HAC matrix based on $\mathbf{g}_t(\boldsymbol{\theta}_0)$ by
\begin{align} \boldsymbol{\Omega}_{p,T} = \boldsymbol{\Gamma}_T(0)+ \sum_{h=1}^{T-1}K_p(\rho_{T,h})\left\{\boldsymbol{\Gamma}_T(h)+\boldsymbol{\Gamma}_T(h)'\right\}, \end{align}
where $\boldsymbol{\Gamma}_T(h) = T^{-1}\sum_{t=h+1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)\mathbf{g}_{t-h}(\boldsymbol{\theta}_0)'$.
Thus $\boldsymbol{\Omega}_{p,T}$ separates the standard HAC approximation problem from the additional effects of estimating and recentering the score contributions.
Consistency of the feasible matched-HAC estimator, established in Theorem~\ref{thm:matched_hac_consistency} below, therefore has two ingredients: consistency of $\boldsymbol{\Omega}_{p,T}$ for $\boldsymbol{\Omega}_0$, and asymptotic negligibility of replacing the true-parameter scores by their estimated and recentered counterparts. The first ingredient is stated in the next assumption.
Throughout this section, the latent correlation sequence $\{\rho_{T,h}\}$ is given by \eqref{eq:latent_gaussian_corr} and is used both to construct the bootstrap multipliers and to define the matched-HAC weights.
The population-score HAC consistency condition is as follows.
\begin{ass}
$\boldsymbol{\Omega}_{p,T}\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0$.
\label{ass:consistency_infHAC}
\end{ass}
Assumption~\ref{ass:consistency_infHAC} holds under standard weak-dependence conditions.
For example, Assumption A of \citet{Andrews1991} implies this condition when the lag-weight function $x\mapsto K_p\{\exp(-x^2)\}$ belongs to the $\mathcal{K}_1$ class.
Lemma~\ref{lem:kernel_radem} verifies this property.
\begin{lem}
The function $x\mapsto K_p\{\exp(-x^2)\}$ is continuous, equals one at $x=0$, is symmetric in $x$, and satisfies $\int_{-\infty}^{\infty}|K_p\{\exp(-x^2)\}|^qdx<\infty$ for each $p\in(0,1)$ and finite $q>0$.
Therefore $x\mapsto K_p\{\exp(-x^2)\}$ belongs to the $\mathcal{K}_1$ class of \citet{Andrews1991}.
\label{lem:kernel_radem}
\end{lem}
Lemma~\ref{lem:kernel_radem} shows that the HAC lag-weight function induced by the two-point multiplier construction satisfies the standard kernel regularity conditions. Hence, under conventional weak-dependence conditions, Assumption~\ref{ass:consistency_infHAC} follows from existing HAC theory. The next result shows that estimation and recentering do not affect this consistency asymptotically.
\begin{thm}[matched-HAC consistency]
Suppose that Assumptions~\ref{ass:moment_weakdepend}--\ref{ass:consistency_infHAC} hold.
If $b_T/T^{1/2}\to0$, then $\hat{\boldsymbol{\Omega}}_p\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0$.
\label{thm:matched_hac_consistency}
\end{thm}
Matched-HAC consistency provides the covariance approximation needed for studentization, but bootstrap validity additionally requires a conditional Gaussian approximation to the bootstrap score sum.
For this purpose, we strengthen the weak-dependence and moment conditions as follows.
\begin{ass}
\begin{itemize}
\item[(i)] $\mathbf{g}_t(\boldsymbol{\theta}_0)$ is strongly mixing with mixing coefficients $\alpha(h)$.
\item[(ii)] There exists $\delta\geq2$ such that $\sup_t\E\|\mathbf{g}_t(\boldsymbol{\theta}_0)\|^{2+\delta}<\infty$ and $\sum_{h=1}^{\infty}\alpha(h)^{\delta/(2+\delta)}<\infty$.
\end{itemize}
\label{ass:moment_stronger}
\end{ass}
Assumption~\ref{ass:moment_stronger} is the same as Assumption 3.1 of \citet{Shao2010}.
Under this stronger condition, the next theorem establishes the central bootstrap approximation at the level of the estimating-equation sum.
\begin{thm}[Conditional Gaussian approximation]
Suppose that Assumptions~\ref{ass:moment_weakdepend}--\ref{ass:moment_stronger} hold.
If $b_T/T^{\delta/(2+2\delta)}\to0$, then
${T}^{-\frac{1}{2}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\boldsymbol{\Omega}_0)$,
${T}^{-\frac{1}{2}}\sum_{t=1}^T\hat{\mathbf{g}}_t^*
\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\boldsymbol{\Omega}_0)$ in probability,
and
\begin{align}
\sup_{\mathbf{x}\in\mathbb{R}^d}
\left|
\Pro\left(\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)\leq\mathbf{x}\right)
-
\Pro^*\left(
\frac{1}{\sqrt{T}}\sum_{t=1}^T\hat{\mathbf{g}}_t^*
\leq\mathbf{x}
\right)
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:score_bootstrap_cdf_validity}
\end{align}
\label{thm:score_boot_consistent}
\end{thm}
Theorem~\ref{thm:score_boot_consistent} is stated at the score level and is the key bootstrap result.
For asymptotically linear estimators, it transfers directly through the corresponding linear representation, yielding the following result.
\begin{cor}[Gaussian approximation for estimators]
Suppose that Assumptions~\ref{ass:moment_weakdepend}--\ref{ass:moment_stronger} hold, $\operatorname{rank}(\mathbf{B}_0)=k$,
$b_T/T^{\delta/(2+2\delta)}\to0$, and $\hat{\mathbf{B}}\stackrel{p}{\longrightarrow}\mathbf{B}_0$.
Then
$\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
\stackrel{d}{\longrightarrow} N(\mathbf{0},\mathbf{V}_{\theta})$,
$\sqrt{T}(\hat{\boldsymbol{\theta}}^*-\hat{\boldsymbol{\theta}})
\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\mathbf{V}_{\theta})$
in probability,
and
\begin{align}
\sup_{\mathbf{x}\in\mathbb{R}^k}
\left|
\Pro\left(\sqrt{T}(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)\leq\mathbf{x}\right)
-
\Pro^*\left(\sqrt{T}(\hat{\boldsymbol{\theta}}^*-\hat{\boldsymbol{\theta}})\leq\mathbf{x}\right)
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:theta_bootstrap_cdf_validity}
\end{align}
\label{cor:theta_hat_boot_consistent}
\end{cor}
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Gaussian validity of the unrestricted and restricted studentized vectors}\label{subsec:wald_lm_validity}
The Gaussian approximation is the primary first-order result for inference, while the conventional Wald and LM chi-square laws follow by taking squared norms.
\begin{cor}[Unrestricted bootstrap validity]\label{cor:wald}
Assume that $H_0:\mathbf{A}\boldsymbol{\theta}=\mathbf{a}_0$ holds, where $\mathbf{A}\in\mathbb{R}^{\ell\times k}$ has rank $\ell$, $\mathbf{a}_0\in\mathbb{R}^\ell$, and $\ell$ is fixed.
Suppose that $\mathbf{A}\mathbf{V}_{\theta}\mathbf{A}'$ is positive definite and $\hat{\mathbf{V}}_{\theta}\stackrel{p}{\longrightarrow}\mathbf{V}_{\theta}$.
Under the conditions of Corollary~\ref{cor:theta_hat_boot_consistent}, let the bootstrap vector use either the fixed studentizer $\hat{\mathbf{V}}_{\theta}$ or a recomputed studentizer satisfying $\hat{\mathbf{V}}_{\theta}^*\stackrel{p^*}{\longrightarrow}\mathbf{V}_{\theta}$ in probability.
Then
$\hat{\mathbf{z}}_{W}
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell)$,
$\hat{\mathbf{z}}_{W}^*
\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell)$ in probability,
and
\begin{align}
\sup_{\mathbf{x}\in\mathbb{R}^\ell}
\left|
\Pro^*(\hat{\mathbf{z}}_{W}^*\leq\mathbf{x})-\Pro(\hat{\mathbf{z}}_{W}\leq\mathbf{x})
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:wald_bootstrap_gaussian_validity}
\end{align}
Consequently, $W_T=\hat{\mathbf{z}}_{W}'\hat{\mathbf{z}}_{W}\stackrel{d}{\longrightarrow}\chi_\ell^2$, $W_T^*=\hat{\mathbf{z}}_{W}^{*'}\hat{\mathbf{z}}_{W}^*\stackrel{d^*}{\longrightarrow}\chi_\ell^2$ in probability, and
\begin{align}
\sup_{x \geq 0}
\left|
\Pro^*(W_T^*\leq x)-\Pro(W_T\leq x)
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:wald_bootstrap_validity}
\end{align}
\end{cor}
For the restricted-score LM result, let
$\mathbf{G}=(\mathbf{G}_1,\mathbf{G}_2)$ denote the probability limit of
$\tilde{\mathbf{G}}=(\tilde{\mathbf{G}}_1,\tilde{\mathbf{G}}_2)$ under $H_0$, and let
$\mathbf{G}_{2\cdot1}
=
\mathbf{G}_2
-
\mathbf{G}_1
(\mathbf{G}_1'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_1)^{-1}
\mathbf{G}_1'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_2$.
\begin{cor}[Restricted bootstrap validity]\label{cor:lm}
Assume $H_0:\mathbf{A}\boldsymbol{\theta}=\mathbf{a}_0$ holds as in Corollary~\ref{cor:wald}.
Suppose that the conditions of Theorem~\ref{thm:score_boot_consistent} hold, $\tilde{\mathbf{G}}_j\stackrel{p}{\longrightarrow}\mathbf{G}_j$ for $j\in\{1,2\}$, $\tilde{\boldsymbol{\Omega}}_p\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0$, $\mathbf{G}_1'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_1$ is nonsingular, and $\mathbf{G}_{2\cdot1}'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_{2\cdot1}$ is positive definite.
Suppose that the restricted score sum admits the following expansion with respect to $\boldsymbol{\eta}$ under the null:
\begin{align}
\sqrt{T}\,\overline{\tilde{\mathbf{g}}}
=
\frac{1}{\sqrt{T}}\sum_{t=1}^T\mathbf{g}_t(\boldsymbol{\theta}_0)
+
\mathbf{G}_1\sqrt{T}(\tilde{\boldsymbol{\eta}}_1-\boldsymbol{\eta}_{1,0})
+
o_p(1),
\qquad
\sqrt{T}(\tilde{\boldsymbol{\eta}}_1-\boldsymbol{\eta}_{1,0})=O_p(1).
\label{eq:lm_restricted_moment_linearization}
\end{align}
For recomputed studentization, suppose additionally that $\tilde{\boldsymbol{\Omega}}_p^*\stackrel{p^*}{\longrightarrow}\boldsymbol{\Omega}_0$ in probability.
Then the fixed-studentizer and recomputed-studentizer versions both satisfy
$\tilde{\mathbf{z}}_{LM}
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell)$,
$\tilde{\mathbf{z}}_{LM}^*
\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell)$
in probability,
and
\begin{align}
\sup_{\mathbf{x}\in\mathbb{R}^\ell}
\left|
\Pro^*(\tilde{\mathbf{z}}_{LM}^*\leq\mathbf{x})-\Pro(\tilde{\mathbf{z}}_{LM}\leq\mathbf{x})
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:lm_bootstrap_gaussian_validity}
\end{align}
Consequently, $LM_T=\tilde{\mathbf{z}}_{LM}'\tilde{\mathbf{z}}_{LM}\stackrel{d}{\longrightarrow}\chi_\ell^2$, $LM_T^*=\tilde{\mathbf{z}}_{LM}^{*'}\tilde{\mathbf{z}}_{LM}^*\stackrel{d^*}{\longrightarrow}\chi_\ell^2$ in probability, and
\begin{align}
\sup_{x \geq 0}
\left|
\Pro^*(LM_T^*\leq x)-\Pro(LM_T\leq x)
\right|
\stackrel{p}{\longrightarrow}0.
\label{eq:lm_bootstrap_validity}
\end{align}
\end{cor}
For recomputed studentization, Corollaries~\ref{cor:wald} and \ref{cor:lm} impose high-level conditions requiring consistency of the recomputed studentizers under the bootstrap law.
Section~\ref{subsec:reg_validity} provides primitive sufficient conditions for these requirements in the regression residual-bootstrap setting.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Primitive validity of regression residual resampling}\label{subsec:reg_validity}
The preceding results concern direct score resampling, and we now give primitive conditions for the regression-specific residual bootstrap in Subsection~\ref{subsec:reg_residual}.
Under $H_0$, write
\begin{align}
y_t=\mathbf{x}_t'\boldsymbol{\theta}_0+u_t,
\qquad
\mathbf{A}\boldsymbol{\theta}_0=\mathbf{a}_0,
\end{align}
and retain the fixed matrix $\mathbf{C}=(\mathbf{C}_1,\mathbf{C}_2)$ from Subsection~\ref{subsec:reg_residual}.
Throughout this subsection, $\mathbf{x}_t=(1,\mathbf{z}_t')'$ and the intercept is left unrestricted by $\mathbf{A}$.
Write $\mathbf{Q}_x=\E(\mathbf{x}_t\mathbf{x}_t')$ and $\mathbf{g}_t=\mathbf{g}_t(\boldsymbol{\theta}_0)=\mathbf{x}_tu_t$.
For least squares, $\mathbf{B}_0=\mathbf{Q}_x^{-1}$, so $\mathbf{V}_{\theta}=\mathbf{Q}_x^{-1}\boldsymbol{\Omega}_0\mathbf{Q}_x^{-1}$.
Let
\begin{align}
\mathbf{G}_1=-\mathbf{Q}_x\mathbf{C}_1,
\qquad
\mathbf{G}_2=-\mathbf{Q}_x\mathbf{C}_2,
\end{align}
and define
\begin{align}
\mathbf{G}_{2\cdot1}
=
\mathbf{G}_2
-
\mathbf{G}_1(\mathbf{G}_1'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_1)^{-1}\mathbf{G}_1'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_2.
\end{align}
For this regression model,
\begin{align}
\left(\mathbf{G}_{2\cdot1}'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_{2\cdot1}\right)^{-1}
=
\mathbf{A}\mathbf{V}_{\theta}\mathbf{A}'.
\label{eq:reg_lm_wald_variance_relation}
\end{align}
Let $g_{t,j}$ denote the $j$th component of $\mathbf{g}_t$, and let $\operatorname{cum}(g_{s,a},g_{t,b},g_{u,c},g_{v,d})$ denote the fourth-order cumulant of $(g_{s,a},g_{t,b},g_{u,c},g_{v,d})$.
The following conditions provide primitive sufficient conditions for the original-sample regression results.
\begin{ass}[Primitive regression conditions]
\label{ass:reg_additional}
\begin{itemize}
\item[(i)] $\E(\mathbf{x}_tu_t)=\mathbf{0}$ and $\mathbf{Q}_x$ is positive definite.
\item[(ii)] The process $\{(\mathbf{x}_t,u_t)\}_{t\in\mathbb{Z}}$ is strictly stationary and strongly mixing with coefficients $\alpha(h)$, and for some $\delta\geq2$, $\E\|\mathbf{x}_tu_t\|^{2+\delta}<\infty$ and $\sum_{h=1}^{\infty}\alpha(h)^{\delta/(2+\delta)}<\infty$.
\item[(iii)] $\E\|\mathbf{x}_t\|^4<\infty$ and $\boldsymbol{\Omega}_0$ is positive definite.
\item[(iv)] For every $j_1,j_2,j_3,j_4\in\{1,\ldots,k\}$,
$\sum_{h_1,h_2,h_3\in\mathbb{Z}}
\left|
\operatorname{cum}(g_{0,j_1},g_{h_1,j_2},g_{h_2,j_3},g_{h_3,j_4})
\right|
<\infty$.
\item[(v)] $b_T\to\infty$ and $b_T/T^{1/2}\to0$.
\end{itemize}
\end{ass}
For fixed $p\in(0,1)$, define the observed matched two-point HAC estimators
\begin{align}
\hat{\boldsymbol{\Omega}}_p
&=
{1\over T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
\hat{\mathbf{g}}_t\hat{\mathbf{g}}_s',
\\
\tilde{\boldsymbol{\Omega}}_p
&=
{1\over T}\sum_{t=1}^T\sum_{s=1}^T
K_p(\rho_{T,t-s})
(\tilde{\mathbf{g}}_t-\overline{\tilde{\mathbf{g}}})
(\tilde{\mathbf{g}}_s-\overline{\tilde{\mathbf{g}}})'.
\label{eq:reg_observed_twopoint_hac}
\end{align}
We use the regression Wald and LM statistics defined in Subsection~\ref{subsec:reg_residual}.
The residual bootstrap requires the additional conditions below.
Condition (i) is needed for the conditional multiplier approximation, including fixed studentization, while conditions (ii)--(iii) are imposed only when the bootstrap HAC studentizer is recomputed.
\begin{ass}[Residual-bootstrap conditions]
\label{ass:reg_bootstrap}
\begin{itemize}
\item[(i)] With $\delta$ as in Assumption~\ref{ass:reg_additional}(ii), $b_T/T^{\delta/(2+2\delta)}\to0$.
\item[(ii)] For a recomputed bootstrap HAC studentizer, $\E\|\mathbf{x}_t\|^{4+2\delta}+\E|u_t|^{4+2\delta}<\infty$.
\item[(iii)] If $b_T^*$ denotes the bootstrap HAC bandwidth, then $b_T^*\to\infty$ and $b_T^*=O(b_T)$.
\end{itemize}
\end{ass}
Let $\hat{\boldsymbol{\Omega}}_p^*$ denote the matched-HAC estimator recomputed from the unrestricted residual-bootstrap moment array $\{\hat{\mathbf{g}}_{t,\mathrm{res}}^*\}$ using bandwidth $b_T^*$, and set $\hat{\mathbf{V}}_{\theta}^*=\hat{\mathbf{Q}}_x^{-1}\hat{\boldsymbol{\Omega}}_p^*\hat{\mathbf{Q}}_x^{-1}$.
Let $\tilde{\boldsymbol{\Omega}}_p^*$ denote the matched-HAC estimator recomputed from the restricted residual-bootstrap moment array $\{\tilde{\mathbf{g}}_{t,\mathrm{res}}^*\}$ after subtracting its bootstrap sample average, using bandwidth $b_T^*$.
For the recomputed residual-bootstrap LM statistic, retain the original-sample projection $\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}$ and use the draw-specific projected covariance $\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\boldsymbol{\Omega}}_p^*\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_{2\cdot1}$ in \eqref{eq:lm_bootstrap_generic}.
\begin{prop}[Regression and residual-bootstrap Gaussian validity]
\label{prop:reg_validity}
Suppose Assumption~\ref{ass:reg_additional} holds.
The studentized vectors below are the feasible statistics defined in Subsection~\ref{subsec:wald_lm_implementations}.
\begin{itemize}
\item[(i)] The unrestricted OLS estimator satisfies
\begin{align}
\sqrt T(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
=
\mathbf{Q}_x^{-1}{1\over\sqrt T}\sum_{t=1}^T\mathbf{g}_t+o_p(1),
\qquad
\sqrt T(\hat{\boldsymbol{\theta}}-\boldsymbol{\theta}_0)
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{V}_{\theta}).
\label{eq:reg_theta_primitive_expansion}
\end{align}
Moreover,
\begin{align}
\hat{\boldsymbol{\Omega}}_p&\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0,
\qquad
\tilde{\boldsymbol{\Omega}}_p\stackrel{p}{\longrightarrow}\boldsymbol{\Omega}_0,
\qquad
\hat{\mathbf{V}}_{\theta}\stackrel{p}{\longrightarrow}\mathbf{V}_{\theta},
\\
\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_{2\cdot1}
&\stackrel{p}{\longrightarrow}
\mathbf{G}_{2\cdot1}'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_{2\cdot1}.
\label{eq:reg_observed_hac_consistency}
\end{align}
Hence under $H_0$,
\begin{align}
\hat{\mathbf{z}}_{W}
&\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell),
\qquad
\tilde{\mathbf{z}}_{LM}
\stackrel{d}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell),
\label{eq:reg_observed_gaussian_limits}
\end{align}
and consequently $W_T\stackrel{d}{\longrightarrow}\chi_\ell^2$ and $LM_T\stackrel{d}{\longrightarrow}\chi_\ell^2$.
\item[(ii)] If Assumption~\ref{ass:reg_bootstrap}(i) also holds, then the residual-bootstrap unrestricted and restricted studentized vectors with the observed matched-HAC studentizers held fixed satisfy
\begin{align}
\hat{\mathbf{z}}_{W}^*
&\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell),
\qquad
\tilde{\mathbf{z}}_{LM}^*
\stackrel{d^*}{\longrightarrow}
N(\mathbf{0},\mathbf{I}_\ell)
\label{eq:reg_boot_fixed_clt}
\end{align}
in probability.
\item[(iii)] If Assumption~\ref{ass:reg_bootstrap}(ii)--(iii) also holds, then
\begin{align}
\hat{\boldsymbol{\Omega}}_p^*&\stackrel{p^*}{\longrightarrow}\boldsymbol{\Omega}_0,
\qquad
\tilde{\boldsymbol{\Omega}}_p^*\stackrel{p^*}{\longrightarrow}\boldsymbol{\Omega}_0,
\qquad
\hat{\mathbf{V}}_{\theta}^*\stackrel{p^*}{\longrightarrow}\mathbf{V}_{\theta},
\\
\tilde{\mathbf{G}}_{2\cdot1}'\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\boldsymbol{\Omega}}_p^*\tilde{\boldsymbol{\Omega}}_p^{-1}\tilde{\mathbf{G}}_{2\cdot1}
&\stackrel{p^*}{\longrightarrow}
\mathbf{G}_{2\cdot1}'\boldsymbol{\Omega}_0^{-1}\mathbf{G}_{2\cdot1}.
\label{eq:reg_boot_recomputed_hac_consistency}
\end{align}
Hence the recomputed residual-bootstrap studentized vectors also satisfy the conditional Gaussian limits in \eqref{eq:reg_boot_fixed_clt}.
Consequently,
\begin{align}
\sup_{\mathbf{x}\in\mathbb{R}^\ell}
\left|
\Pro^*(\hat{\mathbf{z}}_{W}^*\leq\mathbf{x})-\Pro(\hat{\mathbf{z}}_{W}\leq\mathbf{x})
\right|
&\stackrel{p}{\longrightarrow}0,
\\
\sup_{\mathbf{x}\in\mathbb{R}^\ell}
\left|
\Pro^*(\tilde{\mathbf{z}}_{LM}^*\leq\mathbf{x})-\Pro(\tilde{\mathbf{z}}_{LM}\leq\mathbf{x})
\right|
&\stackrel{p}{\longrightarrow}0.
\label{eq:reg_boot_cdf_validity}
\end{align}
In particular, $W_T^*\stackrel{d^*}{\longrightarrow}\chi_\ell^2$ and $LM_T^*\stackrel{d^*}{\longrightarrow}\chi_\ell^2$ in probability for both fixed and recomputed studentization.
\end{itemize}
\end{prop}
Proposition~\ref{prop:reg_validity} gives primitive Gaussian validity for the same fixed $\ell$ linear restrictions considered in the general theory, with the Wald and LM chi-square limits following as consequences.
The proof is given in Appendix~\ref{app:reg_primitive}.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Bandwidth choice}\label{sec:bandwidth}
The validity results above impose rate conditions on the bandwidth $b_T$.
For the Rademacher case, the arcsine kernel in \eqref{eq:arcsine_kernel} has a first-order cusp at the origin, as shown in \eqref{eq:kernel_cusp}, suggesting the standard first-order HAC bandwidth rate $b_T=\lceil c_bT^{1/3}\rceil$.
This choice is permitted under Assumption~\ref{ass:moment_stronger} when $\delta>2$.
In the simulations, we set $c_b=1$.
Sensitivity to the dependence window can be examined by varying the bandwidth scale $c_b$.
For a general two-point law, the same latent Gaussian correlation profile induces the lag-weight function $x\mapsto K_p\{\exp(-x^2)\}$, with the transformation $K_p$ depending on the chosen marginal distribution.
In particular, the simulations use this common latent correlation profile for both the Rademacher and Mammen DWB procedures.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Monte Carlo Experiments}\label{sec:mc}
This section examines the finite-sample performance of the proposed dependent two-point multiplier procedures in nonlinear GMM and linear regression models.
To keep the main comparison focused, we report results for a demanding design combining serial dependence, deterministic heteroskedasticity, and skewed innovations.
Results for alternative innovation distributions, together with two-step GMM results and a comparison with conventional LM and Wald bootstrap tests, are reported in Appendix~\ref{app:additional_mc}.
In both experiments, sample sizes are $T\in\{100,200,400\}$.
For each design, $R=5{,}000$ Monte Carlo replications are conducted with $B=499$ bootstrap draws per replication.
All reported tests concern a single restriction, $\ell=1$, and use either the restricted $z$-statistic $\tilde z_{LM}$ or the unrestricted $z$-statistic $\hat z_W$.
We use the signed $z$-statistics as the primary finite-sample diagnostic because they preserve the direction of departures from the null and allow the lower and upper tails of the sampling approximation to be assessed separately.
The main results therefore report equal-tail two-sided $z$-tests: rejection occurs below the $\alpha/2$ critical value or above the $1-\alpha/2$ critical value, using the standard Normal reference distribution for asymptotic inference and the corresponding conditional bootstrap distribution for bootstrap inference.
For $\ell=1$, the conventional LM and Wald statistics are obtained by squaring $\tilde z_{LM}$ and $\hat z_W$, respectively, and hence fold the two directions together.
For comparison, the corresponding conventional LM and Wald bootstrap tests are reported in Appendix~\ref{app:symmetric_tests}.
The common bandwidth is $b_T=\lceil T^{1/3}\rceil$, which is also used as the moving-block-bootstrap block length.
For the dependent two-point multipliers, the latent Gaussian correlation at lag $h\geq0$ is $\rho_{T,h}=\exp\{-(h/b_T)^2\}$.
The asymptotic comparisons use the Bartlett HAC kernel and the kernels induced by the dependent Rademacher and Mammen constructions.
The bootstrap comparisons use the dependent Rademacher and Mammen procedures, a Gaussian DWB with Bartlett-correlated Gaussian multipliers, and a moving-block-bootstrap benchmark.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Nonlinear GMM}\label{subsec:mc_nonlinear_gmm}
The first experiment considers the nonlinear model
\begin{align}
y_t
&=
\theta_{1,0}+\exp(\theta_{2,0}x_t)
+
\sigma_t u_t,
\qquad t=1,\ldots,T,
\label{eq:mc_nonlinear_model}
\end{align}
where $\sigma_t^2=1/2$ for $t\leq T/2$ and $\sigma_t^2=3/2$ for $t>T/2$.
We set $\theta_{1,0}=0$ throughout and consider the test $H_0:\theta_{2,0}=0.5$, with $\theta_{2,0}=0.5$ for size and $\theta_{2,0}=0.7$ for power.
The regressor follows
\begin{align}
x_t
&=
\rho_x x_{t-1}
+
\sqrt{1-\rho_x^2}\,v_t,
\qquad
v_t\sim N(0,1),
\label{eq:mc_x_process}
\end{align}
with $\rho_x=0.8$.
The serially dependent component $u_t$ is generated as
\begin{align}
u_t
&=
\rho_u u_{t-1}
+
\sqrt{1-\rho_u^2}\,\varepsilon_t,
\label{eq:mc_u_process}
\end{align}
with $\rho_u=0.5$.
The autoregressive processes $x_t$ and $u_t$ are initialized at zero, and the first 300 observations are discarded before retaining the sample of size $T$.
We consider standard Normal, standardized Student-$t_5$, and standardized centered $\chi_1^2$ innovations $\varepsilon_t$.
The main text reports the centered $\chi_1^2$ results, while the Normal and Student-$t_5$ results are given in Appendix~\ref{app:additional_mc}.
Let $\mathbf{z}_t=(1,x_t,x_{t-1})'$ and define the moment contribution $\mathbf{g}_t(\boldsymbol{\theta})=\mathbf{z}_t\{y_t-\theta_1-\exp(\theta_2x_t)\}$.
There are three moment conditions and two parameters.
We consider equal-tail two-sided $z$-tests of $H_0:\theta_2=0.5$ versus $H_1:\theta_2\neq0.5$.
The one-step GMM estimator uses the identity weighting matrix.
For the restricted $z$-statistic $\tilde{z}_{LM}$, $\theta_2$ is fixed at its null value and the nuisance parameter $\theta_1$ is estimated subject to this restriction.
The unrestricted $z$-statistic $\hat{z}_{W}$ is based on the unrestricted one-step GMM estimator.
For the bootstrap version of the restricted $z$-statistic $\tilde{z}_{LM}$, centered restricted moment contributions are resampled and the HAC studentizer is recomputed in each bootstrap draw.
For the bootstrap version of the unrestricted $z$-statistic $\hat{z}_{W}$, the corresponding unrestricted centered moment contributions are used.
The dependent Rademacher and Mammen procedures use their matched-HAC kernels, while the Gaussian DWB uses Bartlett-correlated Gaussian multipliers and Bartlett studentization.
As the moving-block-bootstrap (MBB) benchmark, we use the Hall--Horowitz recentered procedure \citep{HallHorowitz1996}.
Throughout the nonlinear GMM experiment and the empirical application, MBB refers to this recentered version.
The recentering subtracts the conditional bootstrap mean induced by unequal representation of observations near the sample boundaries, thereby imposing the bootstrap analogue of the moment condition.
Further details are reported in Appendix~\ref{app:mbb_recentering}.
\begin{table}[!ht]
\centering
\caption{Empirical size and power of equal-tail two-sided $z$-tests based on the GMM estimator under serially dependent heteroskedastic $\chi_1^2$ innovations}
\label{tab:nlgmm-main-chi2}
\begin{threeparttable}
\small
\setlength{\tabcolsep}{5.0pt}
\begin{tabular}{
l
@{\hspace{1.0em}}
S[table-format=2.2]
S[table-format=2.2]
S[table-format=2.2]
@{\hspace{1.5em}}
S[table-format=2.2]
S[table-format=2.2]
S[table-format=2.2]
}
\toprule
Method
& \multicolumn{3}{c}{Size (\%)}
& \multicolumn{3}{c}{Power (\%)} \\
\cmidrule(lr){2-4}
\cmidrule(lr){5-7}
\multicolumn{1}{r}{$T$}
& {100} & {200} & {400}
& {100} & {200} & {400} \\
\midrule
\multicolumn{7}{l}{\textit{Panel A: Restricted $z$-statistic, $\tilde{z}_{LM}$}} \\
\multicolumn{7}{l}{\textit{Asymptotic}} \\
\quad ASY Bartlett
& 9.04 & 8.58 & 6.46
& 47.42 & 67.72 & 89.22 \\
\quad ASY Rad.\ kernel
& 9.08 & 8.60 & 6.62
& 47.34 & 67.02 & 88.92 \\
\quad ASY Mam.\ kernel
& 9.12 & 8.68 & 6.56
& 47.58 & 67.26 & 89.08 \\
\multicolumn{7}{l}{\textit{Bootstrap}} \\
\quad Dep.\ Rademacher
& 6.12 & 6.68 & 5.50
& 40.38 & 63.08 & 87.14 \\
\quad Dep.\ Mammen
& 11.16 & 11.28 & 8.08
& 52.42 & 72.26 & 90.58 \\
\quad Gaussian DWB
& 6.68 & 7.26 & 5.66
& 42.36 & 65.12 & 88.06 \\
\quad MBB
& 10.18 & 10.04 & 7.20
& 50.60 & 70.38 & 89.70 \\
\midrule
\multicolumn{7}{l}{\textit{Panel B: Unrestricted $z$-statistic, }$\hat{z}_{W}$} \\
\multicolumn{7}{l}{\textit{Asymptotic}} \\
\quad ASY Bartlett
& 11.72 & 9.82 & 7.12
& 64.44 & 80.94 & 94.66 \\
\quad ASY Rad.\ kernel
& 11.82 & 9.96 & 7.18
& 64.64 & 80.80 & 94.62 \\
\quad ASY Mam.\ kernel
& 11.92 & 9.96 & 7.18
& 64.66 & 80.84 & 94.66 \\
\multicolumn{7}{l}{\textit{Bootstrap}} \\
\quad Dep.\ Rademacher
& 9.28 & 8.08 & 6.18
& 60.16 & 78.22 & 93.70 \\
\quad Dep.\ Mammen
& 12.44 & 10.48 & 7.56
& 63.34 & 79.78 & 93.76 \\
\quad Gaussian DWB
& 10.06 & 8.76 & 6.54
& 61.72 & 79.34 & 94.04 \\
\quad MBB
& 9.76 & 8.90 & 6.60
& 58.58 & 76.72 & 92.10 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item \textit{Notes:} All reported tests are equal-tail two-sided $z$-tests for $H_0:\theta_2=0.5$ at the nominal $\alpha=0.05$ level based on the one-step GMM estimator. The size is the rejection frequency when $\theta_2=0.5$, while the power is the rejection frequency when $\theta_2=0.7$. Critical values are given by the $\alpha/2$ and $1-\alpha/2$ quantiles of the respective reference distributions. ASY Bartlett, ASY Rad.\ kernel, and ASY Mam.\ kernel denote asymptotic procedures based on the Bartlett kernel and the kernels induced by the dependent Rademacher and Mammen constructions, respectively. Dep.\ Rademacher, Dep.\ Mammen, and Gaussian DWB denote the corresponding dependent wild bootstrap procedures, and MBB denotes the moving-block bootstrap.
\end{tablenotes}
\end{threeparttable}
\end{table}
The corresponding results are summarized in Table~\ref{tab:nlgmm-main-chi2}, which shows a clear advantage of the bootstrap approximation for the restricted $z$-statistic $\tilde{z}_{LM}$.
At $T=100$, the three asymptotic procedures have rejection frequencies of about $9\%$, whereas Dep.\ Rademacher and Gaussian DWB reduce them to $6.12\%$ and $6.68\%$, respectively.
The improvement persists as the sample size increases, with rejection frequencies of $5.50\%$ and $5.66\%$ at $T=400$.
In contrast, Dep.\ Mammen and MBB remain appreciably oversized, particularly at $T=100$ and $T=200$.
The three asymptotic procedures give almost identical results.
Hence, the substantial differences among the bootstrap procedures cannot be attributed primarily to differences between the Bartlett and matched two-point HAC kernels.
Size distortions are generally larger for the unrestricted $z$-statistic $\hat{z}_{W}$ than for the restricted $z$-statistic $\tilde{z}_{LM}$.
Nevertheless, among the bootstrap procedures, Dep.\ Rademacher generally provides the most successful size control.
Raw rejection frequencies under the alternative differ across procedures, but these differences largely track the corresponding differences in size distortion.
We therefore do not interpret the raw power rankings independently of size.
The results for Normal and Student-$t_5$ innovations reported in Appendix~\ref{app:additional_mc} are qualitatively very similar:
Dep.\ Rademacher continues to provide the most accurate size control for the restricted $z$-statistic, while size distortions are generally larger for the unrestricted $z$-statistic.
The corresponding two-step GMM results are reported in Appendix~\ref{app:twostep_mc} and give qualitatively very similar conclusions: size distortions are generally smaller for the restricted $z$-statistic $\tilde z_{LM}$ than for the unrestricted $z$-statistic $\hat z_W$, and Dep.\ Rademacher again provides the most successful size control for the restricted $z$-statistic.
\@startsection {subsection}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\normalsize\bf}{Linear regression}\label{subsec:mc_regression_lm}
The second experiment considers the linear regression
\begin{align}
y_t
&=
\alpha + \theta x_t + \sigma_t u_t,
\qquad
t=1,\dots,T.
\label{eq:mc_regression}
\end{align}
We set $\alpha=0$. The regressor $x_t$ and error process $u_t$ are generated by \eqref{eq:mc_x_process} and \eqref{eq:mc_u_process} as in the nonlinear GMM experiment with standardized $\chi_1^2$ innovations.
A two-sided test, $H_0:\theta=0$ versus $H_1:\theta\neq0$, is implemented. Empirical size is evaluated at $\theta=0$, while power is evaluated at $\theta=0.20$.
For the restricted $z$-statistic $\tilde{z}_{LM}$, the regression is estimated under $H_0$ and the corresponding score is constructed using the restricted residuals.
The unrestricted $z$-statistic $\hat{z}_{W}$ is based on the unrestricted OLS estimator and its HAC covariance.
The bootstrap procedures use residual resampling, with the HAC studentizer recomputed in each bootstrap draw.
The dependent Rademacher and Mammen procedures use their respective matched-HAC kernels, while Gaussian DWB and MBB use Bartlett studentization.
Results for Normal and Student-$t_5$ innovations are reported in Appendix~\ref{app:additional_mc}.
\begin{table}[!ht]
\centering
\caption{Empirical size and power of equal-tail two-sided $z$-tests based on OLS estimator under serially dependent heteroskedastic $\chi_1^2$ innovations}
\label{tab:reg-main-chi2}
\begin{threeparttable}
\small
\setlength{\tabcolsep}{5.0pt}
\begin{tabular}{
l
@{\hspace{1.0em}}
S[table-format=2.2]
S[table-format=2.2]
S[table-format=2.2]
@{\hspace{1.5em}}
S[table-format=2.2]
S[table-format=2.2]
S[table-format=2.2]
}
\toprule
Method
& \multicolumn{3}{c}{Size (\%)}
& \multicolumn{3}{c}{Power (\%)} \\
\cmidrule(lr){2-4}
\cmidrule(lr){5-7}
\multicolumn{1}{r}{$T$}
& {100} & {200} & {400}
& {100} & {200} & {400} \\
\midrule
\multicolumn{7}{l}{\textit{Panel A: Restricted $z$-statistic, $\tilde{z}_{LM}$}} \\
\multicolumn{7}{l}{\textit{Asymptotic}} \\
\quad ASY Bartlett
& 9.42 & 7.82 & 7.44
& 37.64 & 54.28 & 78.96 \\
\quad ASY Rad.\ kernel
& 9.30 & 7.56 & 7.52
& 37.46 & 54.18 & 78.60 \\
\quad ASY Mam.\ kernel
& 9.38 & 7.60 & 7.56
& 37.72 & 54.32 & 78.72 \\
\multicolumn{7}{l}{\textit{Bootstrap}} \\
\quad Dep.\ Rademacher
& 5.72 & 5.34 & 5.60
& 28.92 & 46.54 & 73.32 \\
\quad Dep.\ Mammen
& 11.76 & 10.50 & 10.32
& 41.12 & 56.68 & 79.86 \\
\quad Gaussian DWB
& 7.26 & 6.14 & 6.34
& 32.56 & 49.42 & 76.04 \\
\quad MBB
& 6.32 & 5.64 & 6.14
& 30.46 & 48.02 & 75.22 \\
\midrule
\multicolumn{7}{l}{\textit{Panel B: Unrestricted $z$-statistic, }$\hat{z}_{W}$} \\
\multicolumn{7}{l}{\textit{Asymptotic}} \\
\quad ASY Bartlett
& 11.88 & 9.32 & 8.40
& 42.42 & 58.08 & 80.62 \\
\quad ASY Rad.\ kernel
& 12.46 & 9.40 & 8.36
& 42.64 & 57.90 & 80.60 \\
\quad ASY Mam.\ kernel
& 12.34 & 9.42 & 8.42
& 42.74 & 58.14 & 80.68 \\
\multicolumn{7}{l}{\textit{Bootstrap}} \\
\quad Dep.\ Rademacher
& 7.88 & 6.16 & 6.42
& 32.94 & 49.40 & 75.00 \\
\quad Dep.\ Mammen
& 11.12 & 9.00 & 8.70
& 38.38 & 53.36 & 76.62 \\
\quad Gaussian DWB
& 8.58 & 7.18 & 6.92
& 35.96 & 52.06 & 77.18 \\
\quad MBB
& 6.44 & 5.70 & 6.18
& 31.64 & 48.88 & 75.62 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item \textit{Notes:} All reported tests are equal-tail two-sided $z$-tests for $H_0:\theta=0$ at the nominal $\alpha=0.05$ level. The size is the rejection frequency when $\theta=0$, while the power is the rejection frequency when $\theta=0.20$. Critical values are given by the $\alpha/2$ and $1-\alpha/2$ quantiles of the respective reference distributions. ASY Bartlett, ASY Rad.\ kernel, and ASY Mam.\ kernel denote asymptotic inference based on the Bartlett kernel and the kernels induced by the dependent Rademacher and Mammen constructions, respectively. Dep.\ Rademacher, Dep.\ Mammen, and Gaussian DWB denote the corresponding dependent wild bootstrap procedures, and MBB denotes the moving-block bootstrap.
\end{tablenotes}
\end{threeparttable}
\end{table}
Table~\ref{tab:reg-main-chi2} summarizes the regression results, which reinforce the main findings from the nonlinear GMM experiment.
For the restricted $z$-statistic $\tilde{z}_{LM}$, Dep.\ Rademacher provides the most accurate size control, followed closely by MBB and then Gaussian DWB, while Dep.\ Mammen is clearly the most oversized.
As in the nonlinear GMM experiment, size distortions are generally larger for the unrestricted $z$-statistic $\hat{z}_{W}$ than for $\tilde{z}_{LM}$, with Dep.\ Mammen being the main exception.
The notable difference is that MBB provides slightly better size control than Dep.\ Rademacher for $\hat{z}_{W}$.
Overall, however, the best finite-sample size control is obtained by combining the restricted $z$-statistic $\tilde{z}_{LM}$ with Dep.\ Rademacher, giving the same main conclusion as in the nonlinear GMM experiment.
The corresponding results for Normal and Student-$t_5$ innovations reported in Appendix~\ref{app:additional_mc} are qualitatively very similar and preserve the main conclusion that the restricted $z$-statistic combined with Dep.\ Rademacher provides particularly accurate finite-sample size control.
The power results are broadly similar across procedures, and the relatively small differences in rejection frequencies appear to reflect mainly the corresponding differences in size distortion.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Empirical Application: Mean Reversion in Short-Term Interest Rates}
\label{sec:empirical}
We illustrate the proposed bootstrap procedures using a nonlinear GMM specification motivated by the Cox--Ingersoll--Ross (CIR) model \citep{CoxIngersollRoss1985}.
In continuous time, measured in years, the short rate $r(\tau)$ follows
$dr(\tau)
=
\kappa\{\mu-r(\tau)\}\,d\tau
+
\sigma\sqrt{r(\tau)}\,dW(\tau)$,
where $\mu$ is the long-run mean level, $\kappa>0$ is the speed of mean reversion, $\sigma>0$ governs the instantaneous volatility, and $W(\tau)$ is a standard Brownian motion.
The model implies the conditional mean
$\E\{r(\tau+\Delta)\mid\mathcal{F}_\tau\}
=
\mu
+
\{r(\tau)-\mu\}\exp(-\kappa\Delta)$,
where $\mathcal{F}_\tau$ denotes the information available at time $\tau$ and $\Delta$ is the horizon measured in years.
Thus, $\mu$ determines the level toward which the short rate reverts, while $\kappa$ determines the speed at which deviations from this level decay, with the associated mean-reversion half-life given by $\log(2)/\kappa$ years.
We use only this conditional-mean implication and do not impose the CIR conditional-variance specification.
The empirical analysis uses the monthly three-month Treasury bill secondary-market rate as an observable proxy for the short rate.
The sample runs from January 1985 to August 2026; see Appendix~\ref{app:cir_data} for further details on the data and their construction.
Let $r_t$ denote the observed rate in month $t$, and let $h$ denote the horizon measured in months.
Evaluating the CIR conditional mean at the monthly observation dates gives
$\E_t(r_{t+h})
=
\mu
+
(r_t-\mu)\exp(-\kappa h/12)$,
where $\E_t(\cdot):=\E(\cdot\mid\mathcal{F}_t)$.
Accordingly, we construct unconditional moment restrictions from this conditional-mean implication using functions of the current short rate as instruments.
Let $x_t=r_t/10$, where this normalization defines the relative weighting of the moments under the identity-weighted criterion. We retain the same normalization in the two-step comparison.
The implied moment restrictions are $\E\mathbf{g}_t(\boldsymbol{\theta}_0)=\mathbf{0}$, where
\begin{align}
\mathbf{g}_t(\boldsymbol{\theta})
&=
\mathbf{z}_t
\left[
r_{t+h}
-
\mu
-
(r_t-\mu)\exp\left(-\kappa h/12\right)
\right],
\qquad
\boldsymbol{\theta}=(\mu,\kappa)',
\label{eq:cir_moment}
\end{align}
with $\mathbf{z}_t=(1,x_t,x_t^2)'$.
This gives three moment conditions for the two parameters $\mu$ and $\kappa$.
We examine the model for $h=6,9,12$ months and test the null hypothesis of a two-year mean-reversion half-life, $H_0:\kappa=\log(2)/2$.
For the empirical comparison, we focus on the restricted $z$-statistic $\tilde z_{LM}$, which delivered the most accurate size control in the Monte Carlo experiments, particularly for the dependent Rademacher procedure.
Table~\ref{tab:cir-h9} reports the main specification with $h=9$ months, while the corresponding results for $h=6$ and $h=12$ are reported in Appendix~\ref{app:cir_additional}.
The one-step GMM estimate implies a mean-reversion half-life substantially longer than the two-year value under the null.
All three asymptotic procedures reject the null at the $5\%$ level.
Among the bootstrap procedures, however, only Gaussian DWB rejects, and its result lies very close to the $5\%$ rejection boundary.
The dependent Rademacher procedure gives the largest bootstrap $p$-value and does not reject the null, consistent with its stronger finite-sample size control in the Monte Carlo experiments.
Thus, in this application, the empirical conclusion at the conventional $5\%$ level depends on the distributional approximation used.
Appendix~\ref{app:cir_additional} reports the one-step results for $h=6$ and $h=12$, together with a two-step GMM robustness check for
the main $h=9$ specification.
The evidence varies with the horizon, while the two-step results preserve the main qualitative comparison between the asymptotic and bootstrap procedures.
\begin{table}[!htb]
\centering
\caption{Equal-tail two-sided $z$-tests of a two-year mean-reversion half-life, $h=9$ months}
\label{tab:cir-h9}
\small
\begin{threeparttable}
\begin{tabular}{lrrrr}
\toprule
Method & $\tilde z_{LM}$ & $p$-value & $q_{0.025}$ & $q_{0.975}$ \\
\midrule
\multicolumn{5}{l}{\textit{Asymptotic}} \\
\quad ASY Bartlett & -2.27 & \multicolumn{1}{l}{$0.023^{**}$} & -1.96 & 1.96 \\
\quad ASY Rad.\ kernel & -2.11 & \multicolumn{1}{l}{$0.035^{**}$} & -1.96 & 1.96 \\
\quad ASY Mam.\ kernel & -2.13 & \multicolumn{1}{l}{$0.033^{**}$} & -1.96 & 1.96 \\
\multicolumn{5}{l}{\textit{Bootstrap}} \\
\quad Dep.\ Rademacher & -2.11 & \multicolumn{1}{l}{0.089} & -2.41 & 2.40 \\
\quad Dep.\ Mammen & -2.13 & \multicolumn{1}{l}{0.074} & -2.32 & 2.28 \\
\quad Gaussian DWB & -2.27 & \multicolumn{1}{l}{$0.050^{**}$} & -2.27 & 2.26 \\
\quad MBB & -2.27 & \multicolumn{1}{l}{0.061} & -2.38 & 2.38 \\
\bottomrule
\end{tabular}
\begin{tablenotes}[flushleft]
\footnotesize
\item \textit{Notes:} The sample is January 1985--August 2026 ($T=491$).
The one-step GMM estimates are $\hat\mu=2.65$ and $\hat\kappa=0.18$, implying a mean-reversion half-life of 3.77 years.
All reported tests are equal-tail two-sided $z$-tests based on the restricted $z$-statistic $\tilde z_{LM}$; the test procedures are defined as in the nonlinear GMM Monte Carlo experiment in Table~\ref{tab:nlgmm-main-chi2}.
Bootstrap $p$-values and quantiles are based on 9,999 draws.
$^{**}$ denotes rejection at the 5\% level.
\end{tablenotes}
\end{threeparttable}
\end{table}
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}{Conclusion}
In this paper, we have developed a general two-point dependent wild bootstrap (DWB) for weakly dependent estimating equations.
Its key feature is that the two-point marginal law and the latent dependence profile can be specified separately.
The construction combines a normalized two-point distribution with a stationary latent Gaussian process through a Gaussian copula transformation, includes dependent Rademacher and Mammen multipliers as special cases, and nests the classical iid two-point wild bootstrap under serial independence.
The induced multiplier autocovariances determine the corresponding HAC lag weights, and the resulting matched-HAC estimator coincides exactly with the conditional covariance of the bootstrap estimating-equation sum.
We establish consistency of this covariance estimator and a conditional Gaussian approximation for the bootstrap score sum, yielding first-order validity for HAC-studentized $z$-statistics and the corresponding Wald and LM statistics.
The Monte Carlo results complement these first-order validity results by showing substantial differences in finite-sample performance across multiplier distributions.
Rademacher DWB generally yields more accurate finite-sample size control for \(z\)-tests than Mammen DWB and Gaussian DWB, with particularly good performance when combined with the restricted \(z\)-statistic \(\tilde z_{LM}\).
At the same time, the corresponding asymptotic HAC procedures behave very similarly, indicating that these differences cannot be attributed primarily to the associated HAC kernels.
The empirical application further illustrates that bootstrap and asymptotic approximations can lead to practically meaningful differences in inference, including different conclusions at conventional significance levels.
The generality of the two-point construction motivates several extensions already being developed in companion manuscripts.
Matsushita, Nishi and Yamagata (\citeyear{MatsushitaNishiYamagata2026}) study higher-order properties of the dependent wild bootstrap, while
Dai, Matsushita and Yamagata (\citeyear{DaiMatsushitaYamagata2026}) extend the approach to large panels, where cross-sectional aggregation yields time-indexed score processes, building on the panel DWB framework of \citet{GaoPengYan2024}.
A distinct extension, motivated by matrix-based spatial multiplier constructions such as \citet{ConleyGoncalvesKimPerron2023},
Nishi and Yamagata (\citeyear{NishiYamagata2026CrossTemporal}) develop bounded two-point multipliers under both cross-sectional and temporal dependence, with applications to spatial and related network models.
Together, these developments illustrate the broader applicability of the proposed construction beyond the weakly dependent time-series setting considered here.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}*{Acknowledgments}
We are grateful to Yukitoshi Matsushita for helpful discussions and useful comments.
This work was supported by JSPS KAKENHI (grant numbers 25KJ0041, 25K00625, 25K05036 and 25H00544).
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}*{Data Availability and Disclosure Statement}
All data used in this study are publicly available. Data sources and replication code are provided in the Online Appendix and supplementary material.
The authors declare no conflicts of interest.
\@startsection {section}{1}{\z@}{-3.5ex plus -1ex minus-.2ex}{2.3ex plus .2ex}{\large\bf}*{Generative AI disclosure}
Generative AI tools, including ChatGPT (OpenAI) and Claude (Anthropic), were used for language editing and coding assistance in the Monte Carlo experiments and empirical analysis.
All research design, methodology, analysis, interpretation, and verification were undertaken by the authors, who take full responsibility for the paper and code.
\bibliographystyle{apalike}
\bibliography{DWB}
\newpage