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.
85,369 characters
Bootstrap Inference under General Two-way Clustering with Serially and Spatially Dependent Common Effects
\maketitle
\begin{abstract}
\noindent This paper develops bootstrap procedures for inference in linear regression
models with two-way clustered data. We characterize the estimator's
asymptotic behavior in five mutually exclusive and exhaustive regimes: three
Gaussian and two non-Gaussian. We establish four impossibility results: heterogeneous score components preclude uniform
consistency; uniform consistency also fails in one
non-Gaussian (infeasible) regime; the infeasible regime is not uniformly distinguishable from a
feasible one; and uniform validity over all feasible regimes rules out uniform
conservativeness over the infeasible regime.
To address the feasible regimes, we propose a data-driven regime classifier
and a projection-based wild bootstrap procedure. The procedure delivers
uniformly valid inference across the four feasible regimes while allowing
serial dependence along the second clustering dimension and spatial dependence
along the first. This combination of regime adaptivity and flexible dependence
is new to the two-way clustering literature. Monte Carlo
simulations confirm the accuracy and flexibility of the proposed methods in
settings with complex clustering structures.
\bigskip{}
\textbf{JEL Classification}: C15, C23, C31, C80
\noindent \medskip{}
\textbf{Keywords}: Bootstrap, clustered
data, two-way clustering, robust inference, wild bootstrap.
\end{abstract}
\vspace*{-0.5cm}
\vfill{}
\thispagestyle{empty} \pagebreak{}
\section{Introduction}
Understanding dependence structures is essential for valid statistical
inference. In many empirical settings, the classical assumption of
independence is violated, especially in time series and panel data
contexts. Ignoring such dependence, whether temporal or cross-sectional,
can lead to biased standard errors, inflated Type I error rates, and
invalid inference, as emphasized by Bertrand, Duflo, and Mullainathan~(\citeyear{bertrand2004much})
and Petersen~(\citeyear{petersen2008estimating}).
This paper addresses a key challenge in modern econometrics: conducting
inference under two-way clustering when the common effects may exhibit
arbitrary forms of serial and spatial dependence. This setting is
practically relevant but, to the best of our knowledge, has not been
studied in the existing literature. We first characterize the asymptotic distribution of the estimator
$\widehat{\bm{\beta}}$ under general two-way clustering structures, yielding five mutually exclusive
and exhaustive regimes, classified by the relative growth of variance
along the two clustering dimensions: three Gaussian and two non-Gaussian.
Building on the insight of Davezies,
Haultf{\ae}uille, and Guyonvarch~(\citeyear{davezies2025analytic}), we show that
heterogeneous score components impose a fundamental limit on uniformly
consistent inference under two-way clustering.
Moreover, based on Menzel~(\citeyear{menzel2021bootstrap}), who shows that
no inference procedure can achieve uniform consistency across the
full function space, we extend and sharpen this result by isolating
a particularly challenging regime. In this infeasible regime, consistent
inference is still impossible when the DGP is fully unspecified. We
further demonstrate that lack of knowledge of the true variance along
both dimensions is the \emph{only} obstacle to attaining uniform consistency
across the full function space. We also establish that uniform validity
across the four feasible regimes necessarily rules out uniform conservativeness
in the infeasible one.
Next, we develop a data-driven classification procedure that distinguishes
among all the feasible regimes using two informative discriminants.
This enables researchers to tailor inference procedures to the prevailing
asymptotic environment. However, we also show that the infeasible
scenario is not uniformly distinguishable from one feasible regime,
imposing fundamental limits on what can be learned from the data without
stronger assumptions.
Third, we introduce a family of projection-based wild bootstrap (PWB)
procedures and establish their \emph{uniform} consistency across different
regimes.\footnote{The code used to implement the methods in this paper is available
at:\\ \url{https://jiahaoecon.github.io/webpage/research/}.} Among these, the PWB-H variant delivers uniformly valid inference
across all four feasible regimes. The theory uses joint large-$N$, large-$T$ asymptotics,
but imposes no restriction on the relative growth of $N$ and $T$.
When temporal and spatial dependencies are absent, some PWB variants
closely resemble the approach of Menzel~(\citeyear{menzel2021bootstrap}).
Unlike Menzel (\citeyear{menzel2021bootstrap})'s bootstrap, which
mixes wild and i.i.d. components, our method applies the wild bootstrap
to projected scores, which preserves the correlations for uniform
consistency. Moreover, PWB shares conceptual similarities with the
method of Juodis~(\citeyear{juodis2021shock}). Our approach differs
mainly in two key respects. First, the procedure
detects whether the limiting distribution is asymptotically Gaussian
and, based on this diagnostic, to adaptively select the appropriate
regime, thereby achieving simultaneous validity across multiple scenarios.
Consequently, even though PWB employs the tuning thresholds as in
Menzel~\citeyearpar{menzel2021bootstrap} and Juodis~\citeyearpar{juodis2021shock},
it does \emph{not} inherit the boundary non-uniformity induced by
the indifference region that is intrinsic to purely threshold-based
tuning. Second, we introduce a scaling adjustment that corrects the
estimation error of the variance components along both clustering
dimensions, which is crucial for validity in regimes where the contribution
of the noise term is asymptotically negligible.
Overall, the paper makes six contributions: (i) extends the framework to allow spatial dependence
along the first dimension; (ii) shows heterogeneous score components preclude uniform consistency; (iii) characterizes five
mutually exclusive and exhaustive asymptotic regimes for two-way clustering;
(iv) pinpoints one infeasible regime where uniform consistency fails,
showing that unknown true variance is the sole obstacle to uniform
consistency over the full function space, and further establishes
that uniform validity across the four feasible regimes necessarily
precludes uniform conservativeness in the infeasible regime; (v)
develops a data-driven regime classifier for the four feasible regimes
and shows that the infeasible regime cannot be uniformly distinguished; and
(vi) proposes a projection-based wild bootstrap family and a hybrid
implementation that is uniformly asymptotically valid across all feasible
regimes. To the best of our knowledge, this level of
generality is new to the two-way clustering literature.
The rest of the paper is organized as follows. Section~\ref{sec:model}
presents the two-way clustering model and the five asymptotic regimes.
Section~\ref{sec:method} introduces the PWB procedures and develops
the theory for bootstrap validity. Section~\ref{sec:simulations}
reports various simulation results under five regimes. Section~\ref{sec:conclusion}
concludes. Technical proofs and additional results are deferred to
the appendix.
\begin{comment}
Juodis (2021) introduces the adaptive wild (AdaWild) bootstrap. Our
method is different from the AdaWild bootstrap in three aspects: first,
one of our method applies , while the AdaWild bootstrap is not valid.
first, it applies a data-dependent tuning parameter that minimizing
the MSE, while for the AdaWild bootstrap, it is unclear if there exists
an optimal data-driven way for an ``optimal'' value. The random
variable is Rademacher distributed, which mimic the fourth moment
and kurtosis correction.
\end{comment}
\section{Model Setting and Five Asymptotic Regimes}
\label{sec:model}
\subsection{Two-way Clustering with Serially and Spatially Dependent Common Effects}
We consider a linear regression model with two clustering dimensions.
Let $i=1,\dots,N$ index clusters in the first dimension and $t=1,\dots,T$
index clusters in the second dimension. For each intersection $(i,t)$,
suppose that
\begin{equation}
\bm{y}_{it}=\bm{X}_{it}\bm{\beta}+\bm{u}_{it},
\end{equation}
where $\bm{y}_{it}$ is an $M_{it}\times1$ vector of outcomes, $\bm{X}_{it}$
is an $M_{it}\times K$ matrix of regressors, $\bm{u}_{it}$ is an
$M_{it}\times1$ vector of disturbances, and $\bm{\beta}$ is a $K\times1$
parameter vector. Here, $M_{it}$ denotes the number of observations
in intersection $(i,t)$. This setup allows for unbalanced panels
as well as multiple observations within a given intersection.
Let
\[
M_{i}=\sum_{t=1}^{T}M_{it},\qquad M_{t}=\sum_{i=1}^{N}M_{it},\qquad M=\sum_{i=1}^{N}\sum_{t=1}^{T}M_{it},
\]
where $M_{i}$ and $M_{t}$ denote the numbers of observations in
cluster $i$ and cluster $t$, respectively, and $M$ is the total
sample size.
Stacking all observations yields
\begin{equation}
\bm{y}=\bm{X}\bm{\beta}+\bm{u},\label{eq: linear regression model}
\end{equation}
where $\bm{y}$ is an $M\times1$ vector, $\bm{X}$ is an $M\times K$
matrix, and $\bm{u}$ is an $M\times1$ vector. We use $\bm{y}_{i},\bm{X}_{i},\bm{u}_{i}$
to denote the subvectors and submatrices associated with cluster $i$
in the first dimension, and $\bm{y}_{t},\bm{X}_{t},\bm{u}_{t}$ to
denote those associated with cluster $t$ in the second dimension.
The ordinary least squares estimator of $\bm{\beta}$ is
\begin{equation}
\widehat{\bm{\beta}}=\left(\bm{X}^{\top}\bm{X}\right)^{-1}\bm{X}^{\top}\bm{y}.
\end{equation}
For simplicity, in the main text we focus on the case in which each
intersection contains the same finite number of observations. The
Internet Appendix IC extends the framework to allow for heterogeneous numbers
of observations and missing intersections.
The main feature of the model is its dependence structure, which
is illustrated in Figure~\ref{fig:RV-generate-multiway-clustering-with-TS}.
Panels~(a)-(c) show the classical independence and one-way clustering dependence structures. Panel~(d)
depicts the conventional two-way clustering setup with no additional
dependence; see, e.g., Davezies,
D'Haultf{œ}uille, and Guyonvarch~(\citeyear{davezies2021empirical},
\citeyear{davezies2022marcinkiewicz}), MacKinnon, Nielsen, and Webb~(\citeyear{mackinnon2021wild}),
and Menzel~(\citeyear{menzel2021bootstrap}). Panel~(e) allows for
serial dependence in the time common effect $\{\bm{\xi}_{t}\}_{t=1}^{T}$; see, e.g., Chiang, Hansen, and Sasaki~(\citeyear{chiang2023standard}), Chen and Vogelsang (\citeyear{chen2023fixed}), Hounyo and Lin (\citeyear{hounyo2025jackknife}, \citeyear{hounyo2024wild}).
Panel~(f) further extends the existing literature by allowing for spatial
dependence in the cross-sectional common effect
$\{\bm{\alpha}_{i}\}_{i=1}^{N}$. This yields a more general two-way
clustering environment with both serial and spatial dependence.
\begin{figure}[t!]
\centering \includegraphics[width=1\textwidth]{cluster} \caption{An illustration of two-way clustering with possibly dependent common
effects. A lighter color indicates weaker dependence. In this figure,
spatial dependence decays as $|i-j|$ increases.}
\label{fig:RV-generate-multiway-clustering-with-TS}
\end{figure}
Following Conley~\citeyearpar{conley1999gmm}, we assume that the
spatial process $\{\bm{\alpha}_{i}\}_{i=1}^{N}$ is indexed by locations
$\{\mathfrak{s}_{i}\}_{i=1}^{N}\subset\mathcal{H}$, where $\mathcal{H}$
is a regular lattice in $\mathbb{R}^{2}$.\footnote{The restriction to
\(\mathbb R^{2}\) is imposed for expositional simplicity. The arguments can
be extended to \(\mathbb R^{c}\), for any fixed integer \(c>0\), with
minor modifications.} Let $\mathfrak{d}_{ij}\equiv\|\mathfrak{s}_{i}-\mathfrak{s}_{j}\|$
denote the true Euclidean distance, satisfying $\mathfrak{d}_{ii}=0$,
$\mathfrak{d}_{ij}=\mathfrak{d}_{ji}$, $\mathfrak{d}_{ij}>0$ for
$i\neq j$, and $\mathfrak{d}_{ij}\leq\mathfrak{d}_{ij'}+\mathfrak{d}_{j'i}$.
The spatial
literature also allows $\mathfrak{d}_{ij}$ to be measured with error.
In particular, one may observe $\widetilde{\mathfrak{d}}_{ij}=\mathfrak{d}_{ij}+\varsigma_{ij}$
where $\varsigma_{ij}$ is bounded measurement error; under mild regularity
conditions, replacing $\mathfrak{d}_{ij}$ by $\widetilde{\mathfrak{d}}_{ij}$
preserves consistency as further discussed in the theory (Conley~\citeyear{conley1999gmm};
see also Conley and Molinari~\citeyear{conley2007spatial}; Kelejian
and Prucha~\citeyear{kelejian2007hac}).
We further
describe the dependence structure using an Aldous-Hoover-Kallenberg (AHK, Aldous~\citeyear{aldous1981representations};
Hoover~\citeyear{hoover1979relations}; Kallenberg~\citeyear{kallenberg1989representation})
representation, a standard device in the multiway clustering literature.
\begin{assumption}\label{as: AHS representation} There exists a
Borel measurable function $f$ such that
\begin{equation}
\left(\bm{y}_{it},\bm{X}_{it},\bm{u}_{it}\right)=f\left(\bm{\alpha}_{i},\bm{\xi}_{t},\bm{\varepsilon}_{it}\right),\label{eq: real DGP}
\end{equation}
where $\{\bm{\alpha}_{i}\}$, $\{\bm{\xi}_{t}\}$, and $\{\bm{\varepsilon}_{it}\}$
are mutually independent random elements with uniform marginals on
$[0,1]$. The sequence $\{\bm{\alpha}_{i}\}$ is a strictly stationary spatially dependent process, $\{\bm{\xi}_{t}\}$
is a strictly stationary serially dependent process, and
$\{\bm{\varepsilon}_{it}\}$ is i.i.d. over $(i,t)$. The function
$f$ is allowed to vary with the sample sizes $N$ and $T$. \end{assumption}
Assumption~\ref{as: AHS representation} allows for dependence beyond the
standard two-way clustering structure. Along the first dimension,
observations may be dependent within a cluster and across nearby clusters
through the spatially correlated common effect \(\bm{\alpha}_{i}\). Along
the second dimension, dependence may arise within a cluster and across
nearby clusters through the serially correlated common effect
\(\bm{\xi}_{t}\). Consequently, observations at intersections \((i,t)\)
and \((i',t')\) may be dependent not only when they share the same cluster
in either dimension, but also when \(i\) and \(i'\) are spatially linked or
when \(t\) and \(t'\) are serially linked.\footnote{All serial and spatial dependence considered in this paper operates,
as is standard, through the common effects \(\bm{\alpha}_i\) and/or
\(\bm{\xi}_t\). Individual-specific
dependence in \(\bm{\varepsilon}_{it}\) would introduce an additional
dependence channel and is outside the scope of the present paper.}
The representation is useful as it allows for a Hoeffding type decomposition
as follows:
\[
\bm{X}_{it}^{\top}\bm{u}_{it}=\bm{a}_{i}+\bm{d}_{t}+\bm{w}_{it}+E\left(\bm{X}_{it}^{\top}\bm{u}_{it}\right),
\]
where
\begin{align*}
\bm{a}_{i} & =E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\alpha}_{i})-E(\bm{X}_{it}^{\top}\bm{u}_{it}),\\
\bm{d}_{t} & =E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\xi}_{t})-E(\bm{X}_{it}^{\top}\bm{u}_{it}),\\
\bm{w}_{it} & =\bm{X}_{it}^{\top}\bm{u}_{it}-E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\alpha}_{i})-E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\xi}_{t})+E(\bm{X}_{it}^{\top}\bm{u}_{it})=\bm{v}_{it}+\bm{e}_{it},\\
\bm{v}_{it} & =E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\alpha}_{i},\bm{\xi}_{t})-E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\alpha}_{i})-E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\xi}_{t})+E(\bm{X}_{it}^{\top}\bm{u}_{it}),\\
\bm{e}_{it} & =\bm{X}_{it}^{\top}\bm{u}_{it}-E(\bm{X}_{it}^{\top}\bm{u}_{it}\mid\bm{\alpha}_{i},\bm{\xi}_{t}).
\end{align*}
Observe that $E\left(\bm{e}_{it}\mid\bm{\alpha}_{i},\bm{\xi}_{t}\right)=0$.
The application of the tower property of conditional expectation yields
that $E\left(\bm{e}_{it}\mid\bm{a}_{i},\bm{d}_{t},\bm{v}_{it}\right)=0$.
Similarly, it holds that $E\left(\bm{v}_{it}\mid\bm{a}_{i},\bm{d}_{t}\right)=0$.
The decomposition is based on a score expansion tailored to the OLS
setting; for a generic $M$-estimator, one may replace ${\bm{X}}_{it}^{\top}{\bm{u}}_{it}$
with the corresponding score (influence function).
To separate the effect of $\bm{\alpha}_{i}$ and $\bm{\xi}_{t}$ in
$\bm{v}_{it}$, we follow Menzel (\citeyear{menzel2021bootstrap})
and assume a low-rank approximation. Given that any square-integrable
function of $(\bm{\alpha}_{i},\bm{\xi}_{t})$ admits an expansion
in a tensor-product orthonormal basis, each component of $\bm{v}_{it}$,
$v_{it,k}$, admits a spectral decomposition:
\begin{equation}
v_{it,k}=\sum_{l,l'=1}^{\infty}c_{ll'f,k}\phi_{l,k}\left(\bm{\alpha}_{i}\right)\psi_{l',k}\left(\bm{\xi}_{t}\right)\label{eq: spectral representation}
\end{equation}
under the $L_{2}\left(F_{\alpha\xi}\right)$ norm on the space of
smooth functions of $\left(\alpha,\xi\right)\in\left[0,1\right]^{2}$.
$F_{\alpha\xi}$ is the joint distribution of $\bm{\alpha}_{i},\bm{\xi}_{t}$.
For each $k$, $\left\{ \phi_{l,k}\left(\bm{\alpha}_{i}\right)\right\} _{l}$
and $\left\{ \psi_{l',k}\left(\bm{\xi}_{t}\right)\right\} _{l'}$ are
orthonormal.
\begin{assumption}\label{as: spectral representation} Assume that
there exists a sequence $\left\{ \bar{c}_{ll'}\right\} $ such that
$c_{ll'f,k}\leq\bar{c}_{ll'}$ for each $l$, $l'$ and $k$, and
$\sum_{l,l'=1}^{\infty}\bar{c}_{ll'}^{2}<\infty$. The first three
moments of $\phi_{l,k}\left(\bm{\alpha}_{i}\right)$ and $\psi_{l',k}\left(\bm{\xi}_{t}\right)$
are uniformly bounded by a constant $B>0$ for all $l$, $l'$, and $k$.
\end{assumption} Assumption~\ref{as: spectral representation} restricts
$v_{it,k}$ to be well-approximated
by the leading components of the SVD-type expansion. This is natural in empirical panels where
unit--time dependence is captured by a small number of interactive
effects, as in common-shock/factor-type specifications widely used
in macro--finance and firm-level panels (e.g., interactive fixed
effects, multi-factor return models, and factor-augmented panel regressions).
We now introduce notations that will be used throughout this paper.
Let $\odot$ denote the Hadamard (elementwise) product. For any matrices
$\bm{A}$, write $\mathbf{A}>0$ to indicate that the matrix $\mathbf{A}$
is positive definite. $\lambda_{\max}(\bm{A})$ and $\lambda_{\min}(\bm{A})$
denote the maximum and minimum eigenvalue of $\bm{A}$, respectively.
As is standard in the bootstrap literature, we write $\bm{A}_{n}^{*}\xrightarrow{P^{*}}\bm{A}$
and $\bm{A}_{n}^{*}\xrightarrow{d^{*}}\bm{A}$ to indicate that the
sequence of bootstrap random matrices $\bm{A}_{n}^{*}$ converges
in bootstrap probability and in bootstrap distribution, respectively,
to $\bm{A}$, as $N,T\to\infty$.
\subsection{Five Asymptotic Regimes}
\label{subsec:Limiting-Distribution-five-scenarios}
Define variances in different dimensions as $\bm{\sigma}_{a,f}^{2}=\frac{1}{N}\sum_{i,j=1}^{N}E(\bm{a}_{i}\bm{a}_{j}^{\top})$,
$\bm{\sigma}_{d,f}^{2}=\frac{1}{T}\sum_{t,t'=1}^{T}E\big(\bm{d}_{t}\bm{d}_{t'}^{\top}\big)$,
$\bm{\sigma}_{e,f}^{2}=\frac{1}{NT}\sum_{i,i'=1}^{N}\sum_{t,t'=1}^{T}E\big(\bm{e}_{it}\bm{e}_{i't'}^{\top}\big)$,
and $\bm{\sigma}_{v,f}^{2}=\frac{1}{NT}\sum_{i,i'=1}^{N}\sum_{t,t'=1}^{T}E\big(\bm{v}_{it}\bm{v}_{i't'}^{\top}\big)$.
Moreover, let $\phi_{l}(\bm{\alpha}_{i})\equiv\big(\phi_{l,1}(\bm{\alpha}_{i}),\ldots,\phi_{l,K}(\bm{\alpha}_{i})\big)^{\top}$,
$\psi_{l'}(\bm{\xi}_{t})\equiv\big(\psi_{l',1}(\bm{\xi}_{t}),\ldots,\psi_{l',K}(\bm{\xi}_{t})\big)^{\top}$,
and $\bm{c}_{ll',f}\equiv(c_{ll'f,1},\ldots,c_{ll'f,K})^{\top}$.
Note that the subscript $f$ indicates dependence on the function
$f$ (which may vary with $N$ and $T$); for notational simplicity,
we suppress the explicit $(N,T)$ dependence in the subscript. Assuming
$E\left(\bm{X}_{it}^{\top}\bm{u}_{it}\right)=0$, we can decompose
the mean of scores $\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\bm{s}_{it}\equiv\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\bm{X}_{it}^{\top}\bm{u}_{it}$
as follows:
\begin{align}
\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\bm{s}_{it}= & \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\left(\bm{a}_{i}+\bm{d}_{t}+\bm{e}_{it}+\sum_{l,l'=1}^{\infty}\bm{c}_{ll',f}\odot\phi_{l}\left(\bm{\alpha}_{i}\right)\odot\psi_{l'}\left(\bm{\xi}_{t}\right)\right)\nonumber \\
= & \frac{\bm{\sigma}_{a,f}}{\sqrt{N}}\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\bm{\sigma}_{a,f}^{-1}\bm{a}_{i}+\frac{\bm{\sigma}_{d,f}}{\sqrt{T}}\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\bm{\sigma}_{d,f}^{-1}\bm{d}_{t}+\frac{1}{\sqrt{NT}}\frac{\bm{\sigma}_{e,f}}{\sqrt{NT}}\sum_{i=1}^{N}\sum_{t=1}^{T}\bm{\sigma}_{e,f}^{-1}\bm{e}_{it}\nonumber \\
& +\frac{1}{\sqrt{NT}}\sum_{l,l'=1}^{\infty}\bm{c}_{ll',f}\odot\left(\frac{1}{\sqrt{N}}\sum_{i=1}^{N}\phi_{l}\left(\bm{\alpha}_{i}\right)\right)\odot\left(\frac{1}{\sqrt{T}}\sum_{t=1}^{T}\psi_{l'}\left(\bm{\xi}_{t}\right)\right)\nonumber \\
\equiv & \frac{\bm{\sigma}_{a,f}}{\sqrt{N}}\bm{Z}_{N}^{a}+\frac{\bm{\sigma}_{d,f}}{\sqrt{T}}\bm{Z}_{T}^{d}+\frac{\bm{\sigma}_{e,f}}{\sqrt{NT}}\bm{Z}_{NT}^{e}+\frac{1}{\sqrt{NT}}\sum_{l,l'=1}^{\infty}\bm{c}_{ll',f}\odot\bm{Z}_{N,l}^{\phi}\odot\bm{Z}_{T,l'}^{\psi}.\label{eq:decomposition score}
\end{align}
Here, the correlation $\left\{ Cov\left(\bm{Z}_{N}^{a},\bm{Z}_{N,l}^{\phi}\right)\right\} _{l}$
are not necessarily zero, given that $\bm{Z}_{N}^{a}$ and $\left\{ \bm{Z}_{N,l}^{\phi}\right\} _{l}$
depend on $\left\{ \bm{\alpha}_{i}\right\} _{i}$. The same for $\left\{ Cov\left(\bm{Z}_{T}^{d},\bm{Z}_{T,l'}^{\psi}\right)\right\} _{l'}$.
All remaining components are pairwise uncorrelated. Hence, $\bm{\sigma}_{NT,f}^{2}\equiv Var\left(\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\bm{s}_{it}\right)=\frac{1}{NT}(T\bm{\sigma}_{a,f}^{2}+N\bm{\sigma}_{d,f}^{2}+\bm{\sigma}_{e,f}^{2}+\bm{\sigma}_{v,f}^{2}).$
Heuristically, the decomposition~\eqref{eq:decomposition score}
implies that for each $k$, the limit law of $\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}{s}_{itk}$
is determined by a collection of (possibly dependent) Gaussian primitives,
with relative magnitudes governed by the variance components $T\sigma_{ak,f}^{2}$,
$N\sigma_{dk,f}^{2}$, and $\sigma_{vk,f}^{2}$, where $\sigma_{ak,f}^{2}$,
$\sigma_{dk,f}^{2}$, and $\sigma_{vk,f}^{2}$ denote the $k$th diagonal
elements of $\bm{\sigma}_{a,f}^{2}$, $\bm{\sigma}_{d,f}^{2}$, and
$\bm{\sigma}_{v,f}^{2}$, respectively. Moreover, denote the (matrix)
square roots of $\sigma_{\bullet k,f}^{2}$ and $\bm{\sigma}_{\bullet,f}^{2}$
by $\sigma_{\bullet k,f}$ and $\bm{\sigma}_{\bullet,f}$, respectively.
\paragraph{(D) Strong clustering in at least one dimension.}
This regime occurs when at least one clustering variance diverges
(D):
\begin{equation}
(\mathrm{D}):\quad T\sigma_{ak,f}^{2}\to\infty\quad\text{or}\quad N\sigma_{dk,f}^{2}\to\infty.\label{eq: both diverge}
\end{equation}
It corresponds to strong clustering along at least one dimension and
generalizes settings such as Condition~(16) of MacKinnon et al. \citeyearpar{mackinnon2021wild},
which keeps $({\sigma}_{ak,f}^{2},{\sigma}_{dk,f}^{2})$ fixed. In
this case, the leading behavior is dominated by $\frac{{\sigma}_{ak,f}}{\sqrt{N}}{Z}_{Nk}^{a}+\frac{{\sigma}_{dk,f}}{\sqrt{T}}{Z}_{Tk}^{d}$.
\paragraph{(V) Vanishing clustering effects.}
Suppose both clustering variances vanish:
\[
T\sigma_{ak,f}^{2}=o(1),\qquad N\sigma_{dk,f}^{2}=o(1).
\]
The limit then hinges on the interaction component ${\sigma}_{vk,f}^{2}$,
yielding two subcases.
\smallskip{}
\subparagraph{(V\&N) Vanishing clustering but non-Gaussian interaction:}
\begin{equation}
(\mathrm{V\&N}):\quad T\sigma_{ak,f}^{2}=o(1),\quad N\sigma_{dk,f}^{2}=o(1),\quad\sigma_{vk,f}^{2}>0.\label{eq: vanish nongaussian}
\end{equation}
This regime features no (or weak) clustering along either dimension,
while the non-Gaussian interaction component $\frac{1}{\sqrt{NT}}\sum_{l,l'=1}^{\infty}c_{ll'k}\,Z_{Nk,l}^{\phi}\,Z_{Tk,l}^{\psi}$
remains non-negligible. Such behavior may arise in interactive factor
designs, where multiplicative dependence is intrinsic.
\smallskip{}
\subparagraph{(V\&G) Vanishing clustering and Gaussian limit:}
\begin{equation}
(\mathrm{V\&G}):\quad T\sigma_{ak,f}^{2}=o(1),\quad N\sigma_{dk,f}^{2}=o(1),\quad\sigma_{vk,f}^{2}=o(1).\label{eq: vanish gaussian}
\end{equation}
In this case, the leading term reduces to $\frac{{\sigma}_{ek,f}}{\sqrt{NT}}{Z}_{NTk}^{e}$,
which is Gaussian. This corresponds to ``no clustering'' beyond
the intersection level (e.g., it generalizes Condition~(17) of MacKinnon
et al. \citeyearpar{mackinnon2021wild}).
\paragraph{(I) Intermediate clustering strength.}
The intermediate regime arises when no clustering variance diverges,
yet at least one remains non-negligible:
\[
\max\{T\sigma_{ak,f}^{2},\,N\sigma_{dk,f}^{2}\}\to\varphi_{k}\in(0,\infty).
\]
Again, two subcases are determined by $\sigma_{vk,f}^{2}$.
\smallskip{}
\subparagraph{(I\&N) Intermediate clustering with non-Gaussian interaction:}
\begin{equation}
(\mathrm{I\&N}):\quad\max\{T\sigma_{ak,f}^{2},\,N\sigma_{dk,f}^{2}\}\to\varphi_{k}\in(0,\infty),\qquad\sigma_{vk,f}^{2}>0.\label{eq: converge non-gaussian}
\end{equation}
All terms in~\eqref{eq:decomposition score} may contribute. In fully
unspecified models, the limit may fail to be consistently estimable
due to contamination from $\{{v}_{it,k}\}$: the estimation error
$\widehat{e}_{NT,k}$ can be of order $O_{P}\!\bigl((NT)^{-1/2}{\sigma}_{vk,f}\bigr)$,
which is not negligible relative to the target $O_{P}\!\bigl((NT)^{-1/2}\bigr)$.
\smallskip{}
\subparagraph{(I\&G) Intermediate clustering with Gaussian limit:}
\begin{equation}
(\mathrm{I\&G}):\quad\max\{T\sigma_{ak,f}^{2},\,N\sigma_{dk,f}^{2}\}\to\varphi_{k}\in(0,\infty),\qquad\sigma_{vk,f}^{2}=o(1).\label{eq: converge gaussian}
\end{equation}
Here the interaction contamination vanishes, and a properly designed
procedure can make the estimation error asymptotically negligible.
\begin{remark}[Subsequence reduction and exhaustive regime classification]\label{rem:subsequence-regimes}
Let $\mathcal{B}$ be the class of Borel measurable set of functions
satisfying Assumptions~\ref{as: AHS representation}--\ref{as: same rate},
and allow $f=f_{NT}\in\mathcal{B}$ to vary with $(N,T)$. The induced
parameter vector is compact such that along any sequence $\{f_{NT}\}$
there exists a \emph{convergent} subsequence $\{f_{N_{k}T_{k}}\}$.
Working without loss of generality along an arbitrary convergent subsequence,
the limit necessarily falls into exactly one of the five regimes \eqref{eq: both diverge}-\eqref{eq: converge gaussian}.
Hence, the five regimes are \emph{mutually exclusive and exhaustive}.
For notational simplicity, we keep writing $\{f_{NT}\}$ in place
of the selected subsequence, with the understanding that the proofs
proceed from an arbitrary convergent sequence and use the subsequence
reduction to formally cover all sequences. \end{remark}
\paragraph{Practical relevance of different regimes.}
Among the five regimes, the two Gaussian cases (D) and (V\&G) are
canonical and empirically common, and are largely covered by existing
methods: (D) corresponds to strong clustering in at least one dimension,
whereas (V\&G) corresponds to essentially no clustering beyond the
intersection level. The practical difficulty is that empirical DGPs
need not fall neatly into either extreme. Regime (I\&G) captures the
\emph{transition} between (V\&G) and (D): clustering is present but
not strong enough to behave as in (D), while it is also too strong
to be safely treated as negligible as in (V\&G). In applications,
the boundary separating these regimes is typically unclear as further
demonstrated below, and may be hard to diagnose in finite samples.
Hence, procedures calibrated only for the two endpoints can be sensitive
to local departures from their target regimes, which makes (I\&G)
practically important and motivates inference that is uniformly valid
across the Gaussian continuum (V\&G)--(I\&G)--(D).
Regime (V\&N) is also practically relevant. It arises when main effects
are weak, yet the interaction component remains non-negligible because
common shocks interact with heterogeneous loadings. This structure
can arise in panels with interactive factor features (e.g., $y_{it}=\alpha_{i}\xi_{t}$),
such as asset-pricing panels (common risk shocks with heterogeneous
exposures), firm--time panels (aggregate shocks with heterogeneous
sensitivities), or shift--share type designs (common shocks with
heterogeneous shares). In such environments, treating the limit as
Gaussian may be misleading. Finally, (I\&N) highlights the intrinsic
difficulty of inference when intermediate clustering coexists with
non-vanishing interaction contamination.
\section{Projection-Based Wild Bootstrap (PWB)}
\label{sec:method}
\subsection{Oracle PWB under Known Variance}
For the bootstrap procedure, it is essential to replicate the dependence
structure of the true DGP. We now provide the algorithm procedure
for the projection-based wild bootstrap method.
\begin{description}
\item [{Algorithm 1.}] \textbf{Projection-Based Wild Bootstrap Algorithm}
\end{description}
\begin{enumerate}[label=\textbf{Step \arabic*:}, leftmargin=*, itemsep=1ex]
\item Regress $\bm{y}$ on $\bm{X}$ to obtain the regression estimate $\widehat{\bm{\beta}}$,
the residual $\widehat{\bm{u}}_{it}$, the empirical score $\widehat{\bm{s}}_{it}=\bm{X}_{it}^{\top}\widehat{\bm{u}}_{it}$.
\item Generate the bootstrap score
\begin{equation}
\bm{s}_{it}^{*b}=\widehat{\bm{a}}_{i}\cdot\eta_{i}^{*b}+\widehat{\bm{d}}_{t}\cdot\eta_{t}^{*b}+\widehat{\bm{w}}_{it}\cdot\eta_{i}^{*b}\eta_{t}^{*b}+\bar{\bm{s}}_{NT}.\label{eq: bootstrap score}
\end{equation}
Here, $\widehat{\bm{a}}_{i}$, $\widehat{\bm{d}}_{t}$, $\widehat{\bm{w}}_{it}$,
$\bar{\bm{s}}_{NT}$ are based on projections of $\widehat{\bm{s}}_{it}$
and serve as sample analogs of $\bm{a}_{i}$, $\bm{d}_{t}$, $\bm{w}_{it}$,
$E\left(\bm{X}_{it}^{\top}\bm{u}_{it}\right)$, respectively. These
components are discussed in further detail below. $b$ is an index
for the bootstrap number. $\eta_{i}^{*b}$ and $\eta_{t}^{*b}$ are bootstrap multipliers discussed below.
\item Compute the bootstrap statistics $\widehat{\bm{\beta}}^{*b}-\widehat{\bm{\beta}}=\left(\bm{X}^{\top}\bm{X}\right)^{-1}\sum_{i}\sum_{t}\bm{s}_{it}^{*b}$.
\item Repeat Step 2 to Step 3 for $B$ times. The confidence interval for
$\bm\varrho^\top(\widehat{\bm{\beta}}-\bm{\beta}_{0})$ is then constructed by the
empirical distribution of $\left\{ \bm\varrho^\top(\widehat{\bm{\beta}}^{*b}-\widehat{\bm{\beta}})\right\} _{b}$
under the null hypothesis $\mathcal{H}_{0}:\bm\varrho^\top \bm{\beta}=\bm\varrho^\top\bm{\beta}_{0}$, where $\bm\varrho$ is a known unit vector.
\end{enumerate}
The spatially correlated bootstrap multipliers are generated as
\[
{\bm{\eta}}_{N}^{*b}
=
({\eta}_{i}^{*b})_{i=1}^{N}
=
\mathbb{K}_{N}^{1/2}\widetilde{\bm{\eta}}_{N}^{*b},\quad\widetilde{\bm{\eta}}_{N}^{*b}=(\widetilde{\eta}_{i}^{*b})_{i=1}^{N},
\]
where \(\{\widetilde{\eta}_{i}^{*b}\}_{i=1}^{N}\) are i.i.d. Rademacher random variables.
Given a bandwidth \(\mathfrak d_N>0\), let
\(
\mathbb{K}_{N}
=
\bigl[
\mathcal K(\mathfrak d_{ij}/\mathfrak d_N)
\bigr]_{i,j=1}^{N},
\)
where \(\mathcal K(\cdot)\) is a bounded kernel on \(\mathbb R\). We use the
Wendland \(C^2\) kernel,
\(
\mathcal K(u)=(1-|u|)_{+}^{4}(4|u|+1),
\)
which is compactly supported and ensures positive semi-definite. The compact
support is convenient for the large-sample theory, while positive
semi-definiteness ensures that \(\mathbb K_N^{1/2}\) is well defined.
Hence, by construction, the covariance structure of the spatially
correlated bootstrap multipliers is induced by \(\mathbb K_N\); in
particular,
\(
{Cov}^{*}
\bigl(
{\eta}_{i}^{*b},
{\eta}_{j}^{*b}
\bigr)
=
\mathcal K(\mathfrak d_{ij}/\mathfrak d_N).
\)
Thus, the proposed bootstrap procedure is designed to preserve the
spatial covariance structure of the original data. In the simulations, we set \(\mathfrak d_N\asymp\lfloor N^{1/8}\rfloor\), which
performs well across the designs considered.\footnote{Kim and Sun~\citeyearpar{kim2011spatial} propose an optimal bandwidth
choice for spatially dependent data by modeling spatial dependence as a
linear process. Their bandwidth rule, however, is not directly applicable
under the general mixing framework maintained in this paper. Deriving an
optimal bandwidth choice for the present setting remains an open question
and is left for future research.} This construction is
closely related to the spatially dependent wild bootstrap of
Conley et al.\ (\citeyear{conley2023bootstrap}).
The serially correlated bootstrap multipliers are generated by a simple
two-state Markov construction. Starting from
\(\eta_{0}^{*b}\), drawn from the Rademacher distribution, define, for
\(t\geq 0\),
\[
\eta_{t+1}^{*b}
=
\begin{cases}
\eta_{t}^{*b}, & \text{with probability } (1+q)/2,\\
-\eta_{t}^{*b}, & \text{with probability } (1-q)/2.
\end{cases}
\]
This construction yields serially dependent Rademacher multipliers, with
\({E}^{*}(\eta_t^{*b})=0\) and
\({Cov}^{*}(\eta_t^{*b},\eta_{t-h}^{*b})=q^{h}\) for
\(h\geq 0\).\footnote{The dependent multipliers $\eta_t^{*b}$ can be applied directly to the wild bootstrap in conventional time series settings; see also the dependent wild bootstrap (DWB) of Shao (\citeyear{shao2010dependent}) and Hounyo (\citeyear{hounyo2023wild}).}
The parameter $q$ quantifies serial dependence. The correlation is
$corr^{*}\left(\eta_{t}^{*b},\eta_{t+\iota}^{*b}\right)=q^{\iota}$,
for each $\iota\geq0$ and each $t\geq0$, which is a form of the
Laplacian kernel $k\left(\iota\right)=\exp\left(-\iota/S_{T}\right)$,
with $S_{T}=\left(\ln q\right)^{-1}$.\footnote{The multipliers $\{\eta_{t}^{*b}\}_{t}$ impose geometric decay, much
stronger than the general $\alpha$-mixing class assumed for the time
process. However, the goal here is only to approximate the long-run
autocovariance structure, and this specification is well-suited for that purpose.} The rule of thumb of selecting $q$ follows the plug-in bandwidth
guideline proposed by Andrews (\citeyear{andrews1991heteroskedasticity})
with kernel $k\left(\iota\right)$:
\[
\widehat{q}=\exp\left(-\omega^{-1/3}T^{-1/3}\right),
\]
where $\omega$ represents a measure of autocorrelation, defined as
follows: for $k=1,\ldots,K$, let $\widehat{s}_{t,k}$ denote the
$k$-th element of $\widehat{\bm{s}}_{t}$, and let $\widehat{\rho}_{k}$
be the coefficient obtained by regressing $\widehat{s}_{t,k}$ on
$\widehat{s}_{t-1,k}$. Then, $\omega$ is given by: $\omega=\sum_{k=1}^{K}\frac{\widehat{\rho}_{k}^{2}}{(1-\widehat{\rho}_{k})^{4}}\Biggl/\sum_{k=1}^{K}\frac{(1-\widehat{\rho}_{k}^{2})^{2}}{(1-\widehat{\rho}_{k})^{4}}.$
Given the estimated score $\left\{ \widehat{\bm{s}}_{it}\right\} _{i,t}$,
we now construct an empirical analog of the decomposition. Specifically,
we project the empirical score onto different dimensions:
\begin{equation}
\widehat{\bm{s}}_{it}=\ddot{\bm{a}}_{i}+\ddot{\bm{d}}_{t}+\ddot{\bm{w}}_{it}+\bar{\bm{s}}_{NT},\label{eq: decompose sit}
\end{equation}
where
\begin{align*}
\ddot{\bm{a}}_{i}= & \frac{1}{T}\sum_{t=1}^{T}\widehat{\bm{s}}_{it}-\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\widehat{\bm{s}}_{it}\equiv\bar{\bm{s}}_{iT}-\bar{\bm{s}}_{NT},\\
\ddot{\bm{d}}_{t}= & \frac{1}{N}\sum_{i=1}^{N}\widehat{\bm{s}}_{it}-\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\widehat{\bm{s}}_{it}\equiv\bar{\bm{s}}_{Nt}-\bar{\bm{s}}_{NT},\\
\ddot{\bm{w}}_{it}= & \widehat{\bm{s}}_{it}-\ddot{\bm{a}}_{i}-\ddot{\bm{d}}_{t}-\bar{\bm{s}}_{NT}=\widehat{\bm{s}}_{it}-\bar{\bm{s}}_{iT}-\bar{\bm{s}}_{Nt}+\bar{\bm{s}}_{NT}.
\end{align*}
For the oracle PWB, $\widehat{\bm{a}}_{i}$, $\widehat{\bm{d}}_{t}$,
and $\widehat{\bm{w}}_{it}$ in $(\ref{eq: bootstrap score})$ depend
on the true variances in two dimensions $\bm{\sigma}_{a,f}^{2}$ and
$\bm{\sigma}_{d,f}^{2}$:
\begin{equation}
\widehat{\bm{a}}_{i}=\bm{\vartheta}_{a}\ddot{\bm{a}}_{i},\quad\widehat{\bm{d}}_{t}=\bm{\vartheta}_{d}\ddot{\bm{d}}_{t},\quad\text{and}\quad\widehat{\bm{w}}_{it}=\ddot{\bm{w}}_{it},\label{eq: Oracle PWB}
\end{equation}
where the scaling terms $\bm{\vartheta}_{a}$ and $\bm{\vartheta}_{d}$
are define as follows:
\begin{equation}
\bm{\vartheta}_{a}=\bm{\sigma}_{a,f}\left(\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathcal{K}\!\left(\frac{\mathfrak{d}_{ij}}{\mathfrak{d}_{N}}\right)\ddot{\bm{a}}_{i}\ddot{\bm{a}}_{j}^{\top}\right)^{-1/2}\label{eq: PWB oracle scaling}
\end{equation}
and
\begin{equation}
\bm{\vartheta}_{d}=\bm{\sigma}_{d,f}\left(\dfrac{1}{T}\sum\limits _{t=1}^{T}\sum_{\tau=1}^{T}q^{\left|t-\tau\right|}\ddot{\bm{d}}_{t}\ddot{\bm{d}}_{\tau}^{\top}\right)^{-1/2}.\label{eq:PWB oracle scaling 2}
\end{equation}
Here, each scaling term applies a whitening transformation (i.e., the negative square root matrix) with the
true variance, ensuring that the normalized statistic reproduces the
corresponding component in the limiting distribution. This whitening
is needed here because $\ddot{\bm{a}}_{i}$ is distorted by the contribution
of $\bm{v}_{it}$, which may not be separable.
\begin{assumption} \label{as: moment and variance} For some $\delta>0$
and $\zeta>1$, (i) $E(\bm{X}_{it}^{\top}\bm{u}_{it})=0$,
$\bm{Q}\equiv E(\bm{X}_{it}^{\top}\bm{X}_{it})>0$, $E(\left\Vert \bm{X}_{it}\right\Vert ^{8(\zeta+\delta)})\le C_{1}<\infty$,
$E(\left\Vert \bm{u}_{it}\right\Vert ^{8(\zeta+\delta)})\le C_{2}<\infty$,
and uniformly over $(N,T)$ such that the corresponding (possibly
$(N,T)$-dependent) variance satisfies $\sigma_{\bullet,f}>0$, the
$8(\zeta+\delta)$th moments of $\sigma_{a,f}^{-1}a_{ik}$, $\sigma_{d,f}^{-1}d_{tk}$,
$\sigma_{e,f}^{-1}e_{it,k}$, and $\sigma_{v,f}^{-1}v_{it,k}$ are
bounded, for each $k$. (ii) $\lambda_{\min}(\lim_{N,T\to\infty}(T\bm{\sigma}_{a,f}^{2}+N\bm{\sigma}_{d,f}^{2}+\bm{\sigma}_{v,f}^{2}+\bm{\sigma}_{e,f}^{2}))>0$.
\end{assumption}
Assumptions~\ref{as: moment and variance}~(i) and (ii) impose standard
moment conditions and require the limiting smallest eigenvalue of the
variance of the sum of scores to be bounded away from zero.
\begin{assumption} \label{as: time mixing} For the same
$\zeta$ and $\delta$ as in Assumption \ref{as: moment and variance},
(i) $\bm{\xi}_{t}$ is a $\alpha$-mixing sequence
with a mixing coefficient $\alpha(s)$ such that $\alpha(s)=O(s^{-\lambda})$
for a $\lambda>2\zeta/(\zeta-1)$.
(ii) $q\to1$ as $T\to\infty$, and $(-\ln q)^{-1}=o(T^{1/2})$. \end{assumption}
Assumption~\ref{as: time mixing} (i) imposes a standard mixing
condition from the time-series literature.\footnote{
It weakens the
dependence restriction in Chiang et al.\ \citeyearpar{chiang2023standard}
by requiring an \(\alpha\)-mixing condition rather than a
\(\beta\)-mixing condition. Chiang et al. \citeyearpar{chiang2023standard} project the sum of products of scores
\(\sum_{t,t'}\bm{s}_{it}\bm{s}_{it'}^{\top}\) onto the second dimension, rather than decomposing \(\bm{s}_{it}\) and studying the resulting
product terms. This approach requires a more delicate treatment of
fourth-order summations.} Assumption~\ref{as: time mixing}~(ii) imposes
a condition on the kernel function $q^{\iota}$. Specifically, under
this assumption, the kernel tends to 1 as $T\to\infty$, and tends
to 0 when $\iota\geq T^{1/2}$ as $T\to\infty$, holding all other
factors constant.
Let $\mathcal{F}_{A}\equiv\sigma(\{\bm{\alpha}_{i}:i\in A\})$ for
$A\subset\{1,\dots,N\}$ and define $dist(A,B)\equiv\min\{\mathfrak{d}_{ij}:i\in A,\,j\in B\}$
and the strong mixing coefficient
\[
\alpha_{d_{1},d_{2}}(r)\equiv\sup\Big\{\big|P(G\cap H)-P(G)P(H)\big|:\;G\in\mathcal{F}_{A},\ H\in\mathcal{F}_{B},|A|\leq d_{1},|B|\leq d_{2},\ dist(A,B)\ge r\Big\}.
\]
\begin{assumption}\label{as:spatial} For the same
$\zeta$ and $\delta$ as in Assumption \ref{as: moment and variance}, (i) The sampling region expands
in two non-opposing directions as $N\to\infty$. (ii)
$\alpha_{\infty,\infty}(r)^{1-\frac{1}{2(\zeta+\delta)-1}}=o(r^{-4})$.
(iii) $\mathfrak{d}_{N}\to\infty$, as $N\to\infty$, and $\mathfrak{d}_{N}=o(N^{1/6})$.
(iv) If only $\widetilde{\mathfrak{d}}_{ij}={\mathfrak{d}}_{ij}+\varsigma_{ij}$ is observed, assume
$\sup_{i,j}|\varsigma_{ij}|\le C_{\varsigma}<\infty$ a.s., $\{\varsigma_{ij}\}$
are independent of $\{\bm{\alpha}_{i}\}$, $\{\bm{\xi}_{t}\}$, and
$\{\bm{\varepsilon}_{it}\}$. $\big[\mathcal{K}(\widetilde{\mathfrak{d}}_{ij}/\mathfrak{d}_{N})\big]_{i,j=1}^{N}$
is symmetric and positive semi-definite. \end{assumption}
Assumption~\ref{as:spatial} collects standard high-level regularity
conditions used in the spatial dependence and spatial HAC literature,
see, e.g., Conley \citeyearpar{conley1999gmm}. Part (i) is an increasing-domain
condition. Part (ii) imposes a sufficiently fast decay of strong mixing
to deliver a CLT for the sample mean and to control the variance of
the HAC estimator. Part~(iii) restricts the growth rate of the bandwidth
$\mathfrak{d}_{N}$ on the regular lattice $\mathcal{H}$.
In particular, it is chosen so that the maximal neighborhood size
is $O(\mathfrak{d}_{N}^{2})=o(N^{1/3})$.
Part (iv) allows the use of noisy distances, with uniformly bounded
measurement error independent of the latent spatial component.
\begin{assumption} \label{as: same rate} $\lambda_{\max}(\bm \sigma_{NT,f}^2)/\lambda_{\min}(\bm \sigma_{NT,f}^2)=O(1)$. \end{assumption}
Assumption~\ref{as: same rate} requires the aggregate variances of different
score components to be of the same order. Importantly, Assumption~\ref{as: same rate} does not impose coordinatewise
homogeneity within each componentwise variance term. For instance, the
diagonal entries of \(\bm\sigma_{e,f}^2\) may have heterogeneous orders
across coordinates, and the same is allowed for
\(\bm\sigma_{a,f}^2\), \(\bm\sigma_{d,f}^2\), and
\(\bm\sigma_{v,f}^2\). This allows different components of the score vector to fall into
different asymptotic regimes. Related restrictions have also been imposed in the recent two-way clustering literature; see, for example, Assumption 5 in Davezies,
D'Haultf{œ}uille, and Guyonvarch~\citeyearpar{davezies2025analytic}. In particular, their innovative Example 2 illustrates that, in the absence of such a condition, a standard least-squares approximation may fail under two-way clustered dependence. Proposition~\ref{prop: impossibility heter score} below extends this insight by showing that, without a
restriction of this type, the difficulty is not specific to least-squares
approximations but applies to all data-dependent procedures.
We denote the distribution of the original statistic and bootstrap
statistic using $P_{NT,f}$ and $P_{NT,f}^{*}$, respectively. $\left\Vert \cdot\right\Vert _{\infty}$
is the Kolmogorov metric.
\begin{proposition}[Impossibility due to heterogeneous scores]
\label{prop: impossibility heter score}
Suppose that $\bm{\beta}=\bm{\beta}_0$. Let \(\mathcal D\) denote the collection of all measurable maps of
the observed data \(\{(\bm y_{it},\bm X_{it})\}_{i,t}\) into distribution
functions. Let \(\mathcal B_0\) be a class of DGPs \(f\) satisfying
Assumptions~\ref{as: AHS representation}--\ref{as:spatial}
and for all coordinate \(k\)'s, condition~\eqref{eq: converge non-gaussian} does not hold.
Assume that, for each \(f\in\mathcal B_0\) and $k$, there exists a deterministic
normalizing sequence \(a_{NTk,f}>0\) such that
\(
a_{NTk,f}\bigl(\widehat\beta_{k}-\beta_{k,0}\bigr)
\overset{d}{\to} L_{kf},
\)
where \(L_{kf}\) has a nondegenerate distribution function \(G_{kf}\). Then there
exists \(\varepsilon>0\) and $k$ such that
\[
\liminf_{N,T\to\infty}
\inf_{\widehat D\in\mathcal D}
\sup_{f\in\mathcal B_0}
P_{NT,f}\!\left(
\bigl\|
\widehat D\bigl(\{(\bm y_{it}^{(f)},\bm X_{it}^{(f)})\}_{i,t}\bigr)-G_{kf}
\bigr\|_{\infty}
>\varepsilon
\right)
>0.
\]
\end{proposition}
Proposition~\ref{prop: impossibility heter score} shows that, even after
excluding the infeasible (I\&N) regime, no feasible
data-dependent procedure can uniformly estimate the limiting law of
\(\widehat\beta_k\) when the aggregate score components are allowed to be
heterogeneous (without Assumption \ref{as: same rate}). In particular, uniform estimation may fail when different
coordinates exhibit different stochastic orders, for example when some
coordinates fall in regime (D) while others fall in the remaining regimes. The difficulty arises because the heterogeneous
coordinates of the score vector are mixed through
\((\bm X^\top \bm X)^{-1}\). This reflects
a fundamental limitation of inference with heterogeneous scores under
two-way clustering.
With Assumption \ref{as: same rate}, each coordinate of $\widehat{\bm{\beta}}-\bm{\beta}_{0}$ converges at the same rate. We define the (infeasible) rate of convergence
for $\widehat{\bm{\beta}}-\bm{\beta}_{0}$, $r_{NT,f}=\min\{\sqrt{N}\sigma_{a1,f}^{-1},\sqrt{T}\sigma_{d1,f}^{-1},\sqrt{NT}\}$.
\begin{theorem}\label{thm: main}
For the oracle PWB, under the null hypothesis
$\mathcal{H}_{0}:\bm\varrho^\top\bm{\beta}=\bm\varrho^\top\bm{\beta}_{0}$,
\begin{equation}
\Bigl\|
P_{NT,f}^{*}\!\left(
r_{NT,f}\bm\varrho^\top\bigl(\widehat{\bm{\beta}}^{*}-\widehat{\bm{\beta}}\bigr)
\right)
-
P_{NT,f}\!\left(
r_{NT,f}\bm\varrho^\top\bigl(\widehat{\bm{\beta}}-\bm{\beta}_{0}\bigr)
\right)
\Bigr\|_{\infty}
\xrightarrow{P} 0
\label{eq:theorem 1}
\end{equation}
holds uniformly over the entire function space $\mathcal{B}$ satisfying
Assumptions~\ref{as: AHS representation}--\ref{as: same rate}.
\end{theorem} Theorem~\ref{thm: main} shows the uniform consistency
result for the oracle PWB over the entire function space $\mathcal{B}$.\footnote{It is possible that \(\widehat{\bm{\beta}}-\bm{\beta}_{0}\) converges at rate
\(O_P(N^{-1/2})\), while a particular linear combination
\(\bm{\varrho}^{\top}(\widehat{\bm{\beta}}-\bm{\beta}_{0})\) converges at the faster
rate \(O_P((NT)^{-1/2})\). This can occur in extreme cases where two limiting
components are very close but not identical. In such cases, the normalization
\(r_{NT,f}\) remains of order \(\sqrt{N}\), so both original and bootstrap distributions converge to the same degenerate
limit. Consequently, the present theory does not directly cover
inference on hypotheses such as \(\beta_1>\beta_2\) in such extreme cases.}
However, this procedure is generally infeasible in practice, as it
requires knowledge of the DGP, which is typically unspecified. Interestingly,
it implies that the sole obstacle to achieving uniform consistency
is the lack of knowledge of the true variances in two dimensions $\bm{\sigma}_{a,f}^{2}$
and $\bm{\sigma}_{d,f}^{2}$.
\begin{comment}
\[
\bm{\vartheta}_{a}=\begin{cases}
1, & \text{if }T\bm{\sigma}_{a,f}^{2}\text{ diverges};\\[8pt]
\frac{{\mu}_{a}^{1/2}}{\sqrt{T}}\left(\dfrac{1}{N}\sum\limits _{i=1}^{N}\ddot{\bm{a}}_{i}\ddot{\bm{a}}_{i}^{\top}\right)^{-1/2}, & \text{if }T\bm{\sigma}_{a,f}^{2}\to{\mu}_{a};\\[12pt]
0, & \text{if }T\bm{\sigma}_{a,f}^{2}=o(1),
\end{cases}
\]
\text{and}
\[
\bm{\vartheta}_{d}=\begin{cases}
1, & \text{if }N\bm{\sigma}_{d,f}^{2}\text{ diverges};\\
\frac{{\mu}_{d}^{1/2}}{\sqrt{N}}\left(\dfrac{1}{T}\sum\limits _{t=1}^{T}\ddot{\bm{d}}_{t}\ddot{\bm{d}}_{t}^{\top}+\frac{1}{T}\sum_{\iota=1}^{T-1}q^{\iota}\sum_{t=1}^{T-\iota}\left(\ddot{\bm{d}}_{t}^{\top}\ddot{\bm{d}}_{t+\iota}+\ddot{\bm{d}}_{t+\iota}^{\top}\ddot{\bm{d}}_{t}\right)\right)^{-1/2}, & \text{if }N\bm{\sigma}_{d,f}^{2}\to{\mu}_{d};\\[12pt]
0, & \text{if }N\bm{\sigma}_{d,f}^{2}=o(1).
\end{cases}
\]
\end{comment}
\subsection{Feasible PWBs with Estimated Variance}
In practice, the true values of $\bm{\sigma}_{a,f}^{2}$ and $\bm{\sigma}_{d,f}^{2}$
are generally unknown to the researchers, and hence we propose the
feasible PWB method. Define the variance estimator of $\bm{\sigma}_{a,f}^{2}$
and $\bm{\sigma}_{d,f}^{2}$:
\begin{align}
\widehat{\bm{\sigma}}_{a}^{2}=EVC\left(\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathcal{K}\!\left(\frac{\mathfrak{d}_{ij}}{\mathfrak{d}_{N}}\right)\ddot{\bm{a}}_{i}\ddot{\bm{a}}_{j}^{\top}-\frac{1}{NT^{2}}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{t=1}^{T}\mathcal{K}\!\left(\frac{\mathfrak{d}_{ij}}{\mathfrak{d}_{N}}\right)\ddot{\bm{w}}_{it}\ddot{\bm{w}}_{jt}^{\top}\right)\label{eq: sigma a estimator}
\end{align}
and
\begin{align}
\widehat{\bm{\sigma}}_{d}^{2}= & EVC\Biggl(\frac{1}{T}\sum_{t=1}^{T}\sum_{\tau=1}^{T}q^{\left|t-\tau\right|}\ddot{\bm{d}}_{t}\ddot{\bm{d}}_{\tau}^{\top}-\frac{1}{N^{2}T}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{\tau=1}^{T}q^{\left|t-\tau\right|}\ddot{\bm{w}}_{it}^{\top}\ddot{\bm{w}}_{it}\Biggl),\label{eq: sigma d estimator}
\end{align}
where $EVC\left(\cdot\right)$ replaces any negative eigenvalue by
zero to ensure positive semi-definiteness.
Lemma~\ref{lemma:limit form-1} demonstrates that the magnitude
of $\widehat{\bm{\sigma}}_{a}^{2}$ can be informative about the true
order of $\bm{\sigma}_{a,f}^{2}$, despite some ambiguity arising
from certain alternative cases. Based on this insight, we define
\begin{align}
\widehat{\bm{a}}_{i}= & \widehat{\bm{\vartheta}}_{a}\ddot{\bm{a}}_{i},\quad\widehat{\bm{d}}_{t}=\widehat{\bm{\vartheta}}_{d}\ddot{\bm{d}}_{t},\quad\text{and}\quad\widehat{\bm{w}}_{it}=\ddot{\bm{w}}_{it},\label{eq: PWB}
\end{align}
where
\begin{align}
\widehat{\bm{\vartheta}}_{a} & =\bm{D}_{a}\left(\widehat{\bm{\mu}}_{a}\right)\cdot\widehat{\bm{\sigma}}_{a}\left(\frac{1}{N}\sum_{i=1}^{N}\sum_{j=1}^{N}\mathcal{K}\!\left(\frac{\mathfrak{d}_{ij}}{\mathfrak{d}_{N}}\right)\ddot{\bm{a}}_{i}\ddot{\bm{a}}_{j}^{\top}\right)^{-1/2},\label{eq:PWB DV scaling}\\
\widehat{\bm{\vartheta}}_{d} & =\bm{D}_{d}\left(\widehat{\bm{\mu}}_{d}\right)\cdot\widehat{\bm{\sigma}}_{d}\left(\frac{1}{T}\sum_{t=1}^{T}\sum_{\tau=1}^{T}q^{\left|t-\tau\right|}\ddot{\bm{d}}_{t}^{\top}\ddot{\bm{d}}_{\tau}\right)^{-1/2}.\label{eq:PWB DV scaling 2}
\end{align}
Here, $\bm{D}_{a}(\widehat{\bm{\mu}}_{a})$ denotes the $K\times K$
diagonal matrix with
\[
\bigl[\bm{D}_{a}(\widehat{\bm{\mu}}_{a})\bigr]_{kk}\;=\;\mathbb{I}\!\left\{ \widehat{\sigma}_{ak}^{2}\ge\widehat{\mu}_{ak}\right\} ,\qquad k=1,\ldots,K.
\]
Define $\bm{D}_{d}(\widehat{\bm{\mu}}_{d})$ analogously by $\bigl[\bm{D}_{d}(\widehat{\bm{\mu}}_{d})\bigr]_{kk}=\mathbb{I}\!\left\{ \widehat{\sigma}_{dk}^{2}\ge\widehat{\mu}_{dk}\right\} $.
Relative to the oracle scalings \(\bm{\nu}_{a}\) and \(\bm{\nu}_{d}\) in
\eqref{eq: PWB oracle scaling} and \eqref{eq:PWB oracle scaling 2}, their
feasible counterparts \(\widehat{\bm{\nu}}_{a}\) and
\(\widehat{\bm{\nu}}_{d}\) also implement a whitening transformation, but replace
the unknown population variances with their sample estimates and incorporate
indicator matrices. These indicators serve as safeguards in regimes where the
variance estimators may be unreliable, preventing spurious rescaling of
components whose variances are poorly estimated.
We propose several variants of the PWB method, each based on different
choices for the tuning parameters $\widehat{\bm{\mu}}_{a}$ and $\widehat{\bm{\mu}}_{d}$.
These parameter choices are designed to help detect whether the variance
components $T\bm{\sigma}_{a,f}^{2}$ and $N\bm{\sigma}_{d,f}^{2}$
diverge, vanish, or remain bounded, which in turn affects the validity
of the bootstrap procedure.
\paragraph{PWB-D (Divergence-sensitive).}
The tuning parameters $\widehat{\bm{\mu}}_{a,D}$ and $\widehat{\bm{\mu}}_{d,D}$
are chosen to diagnose divergence of the scaled variance components
$T\bm{\sigma}_{a,f}^{2}$ and $N\bm{\sigma}_{d,f}^{2}$, respectively.
In particular, their role is to detect whether the data fall into
scenario~(D), i.e., whether condition~\eqref{eq: both diverge}
holds. In practice, we recommend setting
\begin{equation}
\widehat{\bm{\mu}}_{a,D}=\mathbf{1}_K\cdot{\log T}/{T}\quad\text{and}\quad\widehat{\bm{\mu}}_{d,D}=\mathbf{1}_K\cdot{\log N}/{N}.\label{eq: tuning pwb-d}
\end{equation}
\paragraph{PWB-V (Vanishing-sensitive).}
In contrast, the tuning parameters $\widehat{\bm{\mu}}_{a,V}$ and
$\widehat{\bm{\mu}}_{d,V}$ are selected to discriminate whether the
scaled variances $T\bm{\sigma}_{a,f}^{2}$ and $N\bm{\sigma}_{d,f}^{2}$,
together with $\bm{\sigma}_{v,f}^{2}$, vanish; that is, whether the
(V\&G) condition~\eqref{eq: vanish gaussian} holds. In practice,
we recommend
\begin{equation}
\widehat{\bm{\mu}}_{a,V}=\mathbf{1}_K\cdot1/(T\log T)\text{ and }\widehat{\bm{\mu}}_{d,V}=\mathbf{1}_K\cdot1/(N\log N).\label{eq: tuning pwb-v}
\end{equation}
The logarithmic terms in the definitions of \(\widehat{\bm{\mu}}_{a}\)
and \(\widehat{\bm{\mu}}_{d}\) reflect the absence of a sharp finite-sample
boundary between the relevant regimes, a feature intrinsic to the problem
and illustrated further in Section~\ref{sec: DRC}. Such thresholding
rules necessarily generate a shrinking but nonempty indifference region,
within which classification may be unstable. For this reason, procedures
that rely directly on a hard classification may be affected near the
boundary. The specific choices of \(\widehat{\bm{\mu}}_{a}\) and
\(\widehat{\bm{\mu}}_{d}\), however, are not DGP-specific. They are used
only to guide regime classification, and the hybrid procedure recommended below is designed precisely to
avoid relying on a sharp boundary classification; hence these threshold
choices do not affect its uniform validity. The simulation evidence also suggests that its finite-sample performance
is not sensitive to the particular threshold choice, provided that the
thresholding rule can reliably identify the relevant regime.
Define the feasible rate of convergence
\[
\widehat{r}_{NT,f}=\min\{\sqrt{N}\widehat{\sigma}_{a1}^{-1},\sqrt{T}\widehat{\sigma}_{d1}^{-1},\sqrt{NT}\}.
\]
\begin{theorem}\label{thm: main-2} Suppose Assumptions~\ref{as: AHS representation}--\ref{as: same rate}
hold. Under the null hypothesis $\mathcal{H}_{0}:\bm\varrho^\top\bm{\beta}=\bm\varrho^\top\bm{\beta}_{0}$,
and with ${{r}}_{NT,f}$ replaced by its feasible counterpart $\widehat{{r}}_{NT,f}$,
the convergence result in~\eqref{eq:theorem 1} holds uniformly in
the following cases:
\begin{enumerate}[label=(\alph*)]
\item \textbf{(PWB-D).} For each $k$, either one of \eqref{eq: vanish nongaussian},
\eqref{eq: vanish gaussian} holds or
\[
T\sigma_{ak,f}^{2}>2\log T\quad\text{or}\quad N\sigma_{dk,f}^{2}>2\log N.
\]
\item \textbf{(PWB-V).} For each $k$, one of \eqref{eq: both diverge},
\eqref{eq: vanish gaussian}, or \eqref{eq: converge gaussian} holds.
\end{enumerate}
\end{theorem}
Theorem~\ref{thm: main-2}(a) implies that PWB-D is uniformly valid
when clustering vanishes (regimes V\&G and V\&N), and also under sufficiently
strong clustering in at least one dimension (a condition strictly
stronger than D). The precise boundary of such ``strictly'' strong
clustering region is determined by the chosen tuning thresholds $\widehat{\bm{\mu}}_{a,D}$ and
$\widehat{\bm{\mu}}_{d,D}$.
Theorem~\ref{thm: main-2}(b) shows that PWB-V is uniformly consistent
on (almost) all class of DGPs that yield a Gaussian limit. Importantly, the tuning thresholds do not impose
a similarly strict restriction as PWB-D: in Gaussian scenario,
when $T\bm{\sigma}_{a,f}^{2}$ (or $N\bm{\sigma}_{d,f}^{2}$) vanishes,
its estimator $T\widehat{\bm{\sigma}}_{a}^{2}$ (or $N\widehat{\bm{\sigma}}_{d}^{2}$)
vanishes as well, so any misclassification induced by the indicator
$\widehat{\bm{D}}_{a}(\widehat{\bm{\mu}}_{a,V})$ (or $\widehat{\bm{D}}_{d}(\widehat{\bm{\mu}}_{d,V})$)
is innocuous in \eqref{eq:PWB DV scaling} and \eqref{eq:PWB DV scaling 2}. In fact, this suggests that the indicator matrices,
and the associated tuning parameters $\widehat{\bm{\mu}}_{a,V}$ and
$\widehat{\bm{\mu}}_{d,V}$, are superfluous for the validity of PWB-V.
They are retained for subsequent use and notational clarity.
\paragraph{PWB-H (Hybrid).}
By Theorem~\ref{thm: main-2}, PWB-V is consistent whenever the limiting
distribution is Gaussian, while PWB-D is valid in some settings with
non-Gaussian limiting distributions. However, neither method works
in both of the following scenarios: (I\&G) and
(V\&N). This is because the variance estimator cannot provide information
to distinguish these two scenarios. These observations motivate a
hybrid procedure that attempts to adapt to both underlying structures
simultaneously based on a second factor.
We define a data-driven combination of PWB-D and PWB-V via the tuning
parameters
\[
\widehat{\bm{\mu}}_{a,H}=\bm{D}^{*}\widehat{\bm{\mu}}_{a,D}+(\bm{I}_{K}-\bm{D}^{*})\widehat{\bm{\mu}}_{a,V}\quad\text{and}\quad\widehat{\bm{\mu}}_{d,H}=\bm{D}^{*}\widehat{\bm{\mu}}_{d,D}+(\bm{I}_{K}-\bm{D}^{*})\widehat{\bm{\mu}}_{d,V},
\]
where $\bm{D}^{*}$ is a diagonal matrix whose $k$-th diagonal element
is constructed based on a Kolmogorov Smirnov (KS) normality test applied
to bootstrap statistics:
\[
D_{k}^{*}=\mathbb{I}\left\{ P_{{\rm KS}}\left(\widehat{t}_{k}^{*b}\right)_{b=1}^{B}<\kappa\right\} .
\]
Here, $P_{{\rm KS}}\left(\widehat{t}_{k}^{*b}\right)_{b=1}^{B}$ returns
the KS $p$-values computed on each group $\left(\widehat{t}_{k}^{*1},\ldots,\widehat{t}_{k}^{*B}\right)$.
The standardized version bootstrap
statistic $\widehat{t}_{k}^{*b}=\sum_{i,t}s_{it,k}^{*b}\Biggl/\sqrt{\frac{1}{B-1}\sum_{b=1}^{B}\left(\sum_{i,t}s_{it,k}^{*b}\right)^{2}}$
is computed under the PWB-V procedure. As we show that PWB-V reproduces
the \emph{form} of the asymptotic law: the PWB-V bootstrap statistic
is asymptotically Gaussian if and only if the original statistic is
asymptotically Gaussian. Observe that the components of $\widehat{\mu}_{ak,H}$
and $\widehat{\mu}_{dk,H}$ may differ across $k$, allowing each component
to have its own limiting behavior. Under the non-Gaussian alternative, this KS $p$-value converges to
zero at an exponential rate, which motivates the choice $\kappa=1/B$.
Since the number of bootstrap replications $B$ can be chosen large
and is not tied to the sample size, this threshold does not induce
the shrinking indifference region associated with the tuning parameters
$\widehat{{\mu}}_{\bullet}$. Proposition \ref{prop: two indicator factors} further shows that the KS diagnostic distinguishes the Gaussian and
non-Gaussian cases with probability approaching one, including transition regimes that allow for drifting parameters. Moreover, such pre-test classifier
does not induce the type of post-selection distortion discussed by Leeb and P{\"o}tscher
\citeyearpar{leeb2008can}. The
simulation results in Table~\ref{tab:average correct rate for DRC varying NT, ols}
are consistent with this theoretical finding.
Note that PWB-H combines the strengths of PWB-V and PWB-D. When $D_{k}^{*}$
indicates a Gaussian limit, PWB-H applies PWB-V, which covers (almost)
all Gaussian regimes, namely D, I\&G, and V\&G. When the rule indicates
a non-Gaussian limit, PWB-H switches to PWB-D, which covers the non-Gaussian
regime V\&N.
More importantly, although PWB-H uses logarithmic tuning thresholds,
it does \emph{not} inherit the boundary non-uniformity that can arise
for PWB-D in certain regimes. This is because PWB-H avoids the thresholds
that may be misleading: in the non-Gaussian region it relies on D-type
thresholds (for which $\widehat{\bm{\mu}}_{a,D}$ remains well behaved
even when $\widehat{\bm{\mu}}_{a,V}$ can be misleading), whereas
in the Gaussian region it adopts the V-type procedure, which is uniformly
valid even when $\widehat{\bm{\mu}}_{a,D}$ may be misleading.
Consequently, as formalized in Theorem~\ref{thm: main-3}, PWB-H
is uniformly asymptotically exact over all regimes except I\&N, i.e.,
when~\eqref{eq: converge non-gaussian} holds, where uniform consistency
is ruled out by Proposition~\ref{prop: impossibility}(a). This supports
PWB-H as a general-purpose inference procedure.
\begin{theorem}\label{thm: main-3} \textbf{(PWB-H).} Suppose Assumptions~\ref{as: AHS representation}--\ref{as: same rate}
hold. Under the null hypothesis $\mathcal{H}_{0}:\bm\varrho^\top\bm{\beta}=\bm\varrho^\top\bm{\beta}_{0}$,
and with ${{r}}_{NT,f}$ replaced by its feasible counterpart $\widehat{{r}}_{NT,f}$,
the convergence result in~\eqref{eq:theorem 1} holds uniformly for
PWB-H, except when condition~\eqref{eq: converge non-gaussian} holds
for some $k$. \end{theorem}
A natural next question is whether one can go one step further: when
uniform validity fails in the I\&N regime, can the procedure at least
be made conservative there. Proposition~\ref{prop: impossibility}(c)
indicates a fundamental tension: any method that is uniformly valid
in all feasible regimes \emph{cannot} be uniformly conservative in the
I\&N regime. Within the two-way clustering framework studied here, this trade-off helps
clarify the robustness of PWB-H: it delivers uniform validity over the feasible
regimes, while the remaining difficulty in I\&N reflects an intrinsic reflects a limitation that cannot be uniformly resolved without sacrificing validity elsewhere.
\subsection{Discussion}
\subsubsection{Dependence Regime Classification}
\label{sec: DRC}
Figure~\ref{fig: variance estimators} (visualizing Proposition~\ref{prop: two indicator factors})
summarizes the relationship between $T\widehat{\sigma}_{ak}^{2}$
and $T\sigma_{ak,f}^{2}$ across Gaussian and non-Gaussian scenarios.
When $D_{k}^{*}=0$ (Gaussian scenario), $T\widehat{\sigma}_{ak}^{2}$
matches the stochastic order of $T\sigma_{ak,f}^{2}$. When $D_{k}^{*}=1$
(non-Gaussian scenario), this order matching fails: under both non-Gaussian
regimes, $T\widehat{\sigma}_{ak}^{2}=O_{P}(1)$. Consequently, these
two non-Gaussian cases may overlap each other.
There exists no clear practical boundary between neighboring regimes, for
instance, between V\&G and I\&G, or between I\&G and D. Accordingly,
we use a $\log T$ or $\log N$ threshold to separate regimes, which
may create an indifference zone near the cutoff. We acknowledge this
limitation, but it is intrinsic: the data cannot, in general, cleanly
distinguish regimes that differ only in such local asymptotic behavior.
\begin{figure}[t!]
\centering \includegraphics[width=0.8\textwidth]{indicators} \caption{\textbf{Stochastic order of the variance estimator across true variance
regimes.}}
\label{fig: variance estimators}
\end{figure}
Notice that relying solely on ``cluster-strength'' diagnostics (building
on the variance estimators) can be misleading, since Gaussian and
non-Gaussian regimes may yield identical diagnostics. Hence, we next
introduce a Dependence Regime Classifier (DRC) to distinguish different
feasible regimes. For each $k$, define the indicators
\[
\widehat{D}_{D,k}=\max\!\left\{ \bigl[\bm{D}_{a}(\widehat{\bm{\mu}}_{a,D})\bigr]_{kk},\;\bigl[\bm{D}_{d}(\widehat{\bm{\mu}}_{d,D})\bigr]_{kk}\right\} ,\qquad\widehat{D}_{V,k}=\max\!\left\{ \bigl[\bm{D}_{a}(\widehat{\bm{\mu}}_{a,V})\bigr]_{kk},\;\bigl[\bm{D}_{d}(\widehat{\bm{\mu}}_{d,V})\bigr]_{kk}\right\} .
\]
\begin{description}
\item [{\textbf{Algorithm 2.}}] \textbf{Dependence Regime Classifier (DRC)}
\end{description}
\begin{enumerate}[label=\textbf{Step \arabic*:}, leftmargin=*, itemsep=1ex]
\item \textbf{Gaussian vs.\ non-Gaussian.} If $D_{k}^{*}=0$, treat the limiting distribution of $\frac{1}{NT}\sum_{i,t}{s}_{itk}$ as Gaussian
and proceed to Step~2. Otherwise, treat it as non-Gaussian and proceed
to Step~3.
\item \textbf{Gaussian branch.} If $\widehat{D}_{D,k}=1$, classify the
component as (D). Else if $\widehat{D}_{V,k}=0$, classify it as (V\&G).
Otherwise, classify it as lying in the \emph{Gaussian transition region}
between (V\&G) and (D) (i.e., one of V\&G, I\&G, or D). For simulation
reporting, we assign this case to (I\&G).
\item \textbf{Non-Gaussian branch.} Classify it as belonging to the \emph{non-Gaussian
region} (i.e., one of V\&N or I\&N).
\end{enumerate}
Finally, note that the two non-Gaussian regimes cannot be distinguished
by DRC. One might hope to construct a sharper classifier, but Proposition~\ref{prop: impossibility}(b)
shows that, absent additional information on the DGP, no procedure
can uniformly distinguish these two regimes.
\subsubsection{Comparison of Existing Methods and Key Differences}
Table~\ref{tab: asymptotic property of methods} provides a summary
of the asymptotic properties of the limiting distributions and the
validity of various inference methods across different regimes. AdaWild
denotes the autoregressive double adaptive wild bootstrap introduced
by Juodis~(\citeyear{juodis2021shock}), CHS refers to the
variance estimator proposed by Chiang et al.~(\citeyear{chiang2023standard}),
and MWCB stands for the multiway cluster bootstrap developed by Hounyo
and Lin~(\citeyear{hounyo2024wild}).
Several inference methods have been proposed under the assumption
of no temporal dependence. For example, the variance estimator of
Cameron, Gelbach, and Miller~(\citeyear{cameron2011robust}), the
wild bootstrap procedures of MacKinnon et al. (\citeyear{mackinnon2021wild}),
and several bootstraps of Menzel~(\citeyear{menzel2021bootstrap})
are all designed without explicitly accounting for autocorrelation.
Nonetheless, their properties can be understood within the general
framework developed here. Particularly when $\{\bm{\xi}_{t}\}_{t}$
do not capture serial dependence, the behavior of Cameron et al.~(\citeyear{cameron2011robust})
estimator and MacKinnon et al.~(\citeyear{mackinnon2021wild}) wild
bootstrap methods are similar to the CHS variance estimator and MWCB,
respectively, and Menzel~(\citeyear{menzel2021bootstrap}) bootstrap
procedure with (without) model selection shares the same asymptotic
validity as our PWB-D (PWB-V) method.
The idea of decomposing the score into three components prior to bootstrapping
also appears in Menzel~\citeyearpar{menzel2021bootstrap} and Juodis~\citeyearpar{juodis2021shock}.
We summarize several key differences between these approaches and
our PWB framework (under no temporal dependence):
\begin{enumerate}
\item Menzel~\citeyearpar{menzel2021bootstrap} employs i.i.d.\ resampling
with wild weights, whereas Juodis~\citeyearpar{juodis2021shock}
and PWB use wild weights to reproduce dependence. This further complicates
the bootstrap joint CLT because the resulting bootstrap components
are no longer independent.
\item Our data-dependent rescaling (introduced to maintain validity in I\&N)
and the thresholding indicators are conceptually related to the bootstrap
with model selection in Menzel~\citeyearpar{menzel2021bootstrap}.
In Juodis~\citeyearpar{juodis2021shock}, thresholding indicators
appear, but without the accompanying rescaling weights.
\item Unlike Menzel~\citeyearpar{menzel2021bootstrap} and Juodis~\citeyearpar{juodis2021shock},
we propose a hybrid bootstrap (PWB-H) that detects whether the
self-normalized statistic is asymptotically Gaussian. The guiding
principle is to deploy the divergence-sensitive and vanishing-sensitive
schemes precisely in the regimes where they are most reliable.
\item Consequently, relative to existing methods, the hybrid bootstrap offers
two main advantages: (i) it is valid on the union of the regimes covered
by the divergence- and vanishing-sensitive schemes, and (ii) it avoids
the boundary non-uniformity induced by the indifference region inherent
in threshold-based tuning.
\end{enumerate}
\begin{table}
\resizebox{\columnwidth}{!}{
{\Huge
\begin{tabular}{llcccccc}
\hline \hline
\multicolumn{2}{l}{\multirow{2}{*}{\makecell{Limiting behavior of \\$\bm{\sigma}_{a,f}^{2}$, $\bm{\sigma}_{d,f}^{2}$, and $\bm{\sigma}_{v,f}^{2}$}}} &
\multirow{2}{*}{\makecell{Divergent, $\left(\ref{eq: both diverge}\right)$}} &
\multicolumn{2}{c}{Vanish}& &
\multicolumn{2}{c}{\makecell[c]{Intermediate}}\\
\cline{4-5}\cline{7-8}
&&& V\&N, \eqref{eq: vanish nongaussian} &
V\&G, \eqref{eq: vanish gaussian} &&
I\&N, \eqref{eq: converge non-gaussian} &
I\&G, \eqref{eq: converge gaussian} \\
\hline
\multirow{3}{*}{\makecell[l]{Asymptotic \\
Properties}} & {Limiting form of $\widehat{\bm{\beta}}$} & Gaussian & non-Gaussian & Gaussian && non-Gaussian & Gaussian \\
&{Order of $\widehat{\bm{\beta}}-\bm{\beta}$} & $O_P\left(\max\left\{\frac{\bm{\sigma}_{a,f}}{\sqrt{N}},\frac{\bm{\sigma}_{d,f}}{\sqrt{T}}\right\}\right)$ & $O_P\left(\frac{1}{\sqrt{NT}}\right)$ & $O_P\left(\frac{1}{\sqrt{NT}}\right)$ & &$O_P\left(\frac{1}{\sqrt{NT}}\right)$ & $O_P\left(\frac{1}{\sqrt{NT}}\right)$ \\
&{Order of $\max\{T\widehat{\bm{\sigma}}_{a}^{2},N\widehat{\bm{\sigma}}_{d}^{2}\}$} & diverge & $O_P\left(1\right)$&$o
_P\left(1\right)$ & &$O_P\left(1\right)$&$O_P\left(1\right)$ \\
\hline
\multirow{7}{*}{\makecell[l]{Methods}} &{CHS CRVE} &$\checkmark$&&$\checkmark$&&&$\checkmark$\\
&
MWCB & $\checkmark$ & & $\checkmark$ & & & $\checkmark$ \\
&AdaWild & $\checkmark$ & $\checkmark$ & $\checkmark$ && & \\
&Oracle PWB & $\checkmark$ & $\checkmark$ & $\checkmark$ & & $\checkmark$ & $\checkmark$ \\
&PWB-D & $\checkmark$ & $\checkmark$ & $\checkmark$ & & & \\
&PWB-V & $\checkmark$ & & $\checkmark$ & & & $\checkmark$ \\
&PWB-H & $\checkmark$ & $\checkmark$ & $\checkmark$ & & & $\checkmark$ \\
\hline \hline
\end{tabular}}}\caption{\textbf{Asymptotic properties of statistics and various methods.} A checkmark is put when the method is consistent in the sub-scenario (ignore the indifference region).
}\label{tab: asymptotic property of methods}
\end{table}
\subsubsection{Other Possible Alternatives}
In this paper, the PWB methods are primarily built on the empirical
score $\widehat{\bm{s}}_{it}$ and employ a non-studentized statistic
designed to accommodate a broad range of settings. When perturbing
the score, we cannot directly compute the bootstrap empirical score
$\widehat{\bm{s}}_{it}^{*}$, which is required for forming a studentized
statistic. One possible workaround is
\begin{equation}
\widehat{\bm{s}}_{it}^{*}=\bm{X}_{it}^{\top}\widehat{\bm{u}}_{it}^{*}=\bm{X}_{it}^{\top}\widehat{\bm{u}}_{it}-\bm{X}_{it}^{\top}\bm{X}_{it}(\widehat{\bm{\beta}}^{*}-\widehat{\bm{\beta}})=\widehat{\bm{s}}_{it}-\bm{X}_{it}^{\top}\bm{X}_{it}(\bm{X}^{\top}\bm{X})^{-1}\sum_{i}\sum_{t}\bm{s}_{it}^{*}.\label{eq: bootstrap empirical score}
\end{equation}
However, this approach is essentially perturbing residuals, as the
first equality in (\ref{eq: bootstrap empirical score}) fixes $\bm{X}_{it}$.
Another alternative is to perturb the score directly, such as $\widehat{\bm{s}}_{i}^{*}=\widehat{\bm{s}}_{i}\eta_{i}^{*}$.
But under Rademacher weights, the resulting bootstrap variance estimator
becomes identical to that of the original $t$-statistic, yielding
no gain relative to the non-studentized approach. We also experimented
with other weight distributions, but found no notable improvement.
Note that we multiply both $\widehat{\bm{a}}_{i}$ and $\widehat{\bm{w}}_{it}$
by the same $\eta_{i}^{*b}$, and both $\widehat{\bm{d}}_{t}$ and
$\widehat{\bm{w}}_{it}$ by the same $\eta_{t}^{*b}$. This design
enables the bootstrap to capture the correlation structures reflected
in $\left\{ Cov\left(\bm{Z}_{N}^{a},\bm{Z}_{N,l}^{\phi}\right)\right\} _{l}$
and $\left\{ Cov\left(\bm{Z}_{T}^{d},\bm{Z}_{T,l}^{\psi}\right)\right\} _{l}$,
which is necessary for the validity of oracle PWB over the entire
parameter function space. By contrast, a pure Efron (or pigeonhole)
bootstrap is ill-suited for reproducing these two underlying correlation
structures of the data, because such resampling tends to overinflate
the variance of $\bm{e}_{it}$.
\begin{comment}
\begin{table}
\resizebox{\columnwidth}{!}{
\begin{tabular}{llcccccc}
\hline
\multicolumn{2}{l}{\multirow{2}{*}{\makecell{Limiting behavior of }}} & & & & & & \tabularnewline
$\bm{\sigma}_{a,f}^{2}$, $\bm{\sigma}_{d,f}^{2}$, and $\bm{\sigma}_{v,f}^{2}$ & $T\bm{\sigma}_{a,f}^{2}$ or $N\bm{\sigma}_{d,f}^{2}$ diverge, $\left(\ref{eq: both diverge}\right)$ & \multicolumn{2}{c}{$T\bm{\sigma}_{a,f}^{2}=o(1)$ and $N\bm{\sigma}_{d,f}^{2}=o(1)$,
$\left(\ref{eq: both diminish}\right)$} & & \multicolumn{2}{c}{\makecell[c]{ $T\bm{\sigma}_{a,f}^{2}\to{\mu}_{a}$ or $N\bm{\sigma}_{d,f}^{2}\to{\mu}_{d}$,
with neither}} & \tabularnewline
term diverging asymptotically, $\left(\ref{eq: both converge}\right)$ & & & & & & & \tabularnewline
\hline
& & & $(i)\bm{\sigma}_{v,f}^{2}>0$ & $(ii)\bm{\sigma}_{v,f}^{2}=o(1)$ & & $(i)\bm{\sigma}_{v,f}^{2}>0$ & $(ii)\bm{\sigma}_{v,f}^{2}=o(1)$ \tabularnewline
\hline
\multicolumn{2}{l}{Limit } & Gaussian & non-Gaussian & Gaussian & & non-Gaussian & Gaussian \tabularnewline
\multicolumn{2}{l}{Convergence rate} & $\sqrt{N}\bm{\sigma}_{a,f}$ or $\sqrt{T}\bm{\sigma}_{d,f}$ & $\sqrt{NT}$ & $\sqrt{NT}$ & & $\sqrt{NT}$ & $\sqrt{NT}$ \tabularnewline
\hline
\multirow{4}{*}{\makecell[l]{Bootstrap}} & & & & & & & \tabularnewline
& MWCB & valid & invalid & valid & & invalid & valid \tabularnewline
& AdaWild & valid & valid & valid & & invalid & invalid \tabularnewline
& Oracle PWB & valid & valid & valid & & valid & valid \tabularnewline
& PWB & valid & valid & valid & & invalid & valid \tabularnewline
\hline
\end{tabular}} \caption{\textbf{Asymptotic properties of various bootstrap methods.}}
Valid means it is uniformly consistent in this case, otherwise invalid.
\end{table}
\end{comment}
\section{Simulation Results}
\label{sec:simulations} In this section, we examine the performance
of various methods and report the most relevant results. The additional
results, including results for less favorable alternatives, heteroskedasticity framework, varying levels of dependence, nonseparable panel models, are provided in the Internet Appendix ID.\footnote{The nonseparable panel DGP is discussed in Fern{\'a}ndez-Val, Freeman, and Weidner \citeyearpar{fernandez2021low} and Chen, Fern{\'a}ndez-Val, Weidner \citeyearpar{chen2021nonlinear}.)} The overall pattern mirrors the findings reported in the main text.
The main goal of our simulation is to support the theoretical results.
We generate data based on the linear model:
\begin{align}
& y_{it}=\beta_{1}+\sum_{k=2}^{K}\beta_{k}X_{it,k}+u_{it},\label{eq: simulation DGP}\\
& X_{it,k}=f_{1k}(\alpha_{i,k}^{x},\xi_{t,k}^{x},\varepsilon_{it,k}^{x}),\text{ and}\\
& u_{it}=f_{2}(\alpha_{i}^{u},\xi_{t}^{u},\varepsilon_{it}^{u}).
\end{align}
We assume that $(\alpha_{i,k}^{x},\alpha_{i}^{u},\xi_{t,k}^{x},\xi_{t}^{u},\varepsilon_{it,k}^{x},\varepsilon_{it}^{u})$
are mutually independent random variables. $(\varepsilon_{it,k}^{x},\varepsilon_{it}^{u})$ are standard Gaussian distributed, independent across $i$ and $t$, and $k$. The latent
components $\left(\xi_{t,k}^{x},\xi_{t}^{u}\right)$ are serially
dependent over $t$, following an AR(1) procedure:
\begin{equation}
\xi_{t}=\rho\xi_{t-1}+\widetilde{\xi}_{t},\text{ where \ensuremath{\widetilde{\xi}_{t}} are independent draws from \ensuremath{\mathcal{N}(0,1-\rho^{2}).}}\label{eq: simulation DGP xi time}
\end{equation}
Such an AR(1) process with $\rho<1$ satisfies Assumption~\ref{as: time mixing}.
The latent
components $\left(\alpha_{i,k}^{x},\alpha_{i}^{u}\right)$ are generated by adopting the spatial design of Conley and Molinari~\citeyearpar{conley2007spatial},
which specifies a stationary, finite-range moving-average field with
geometrically decaying weights. Let $\{\mathfrak{s}_{i}\}_{i=1}^{N}\subset\mathbb{R}^{2}$
denote the (fixed) spatial locations and define the Euclidean distance
$\mathfrak{d}_{ij}\equiv\|\mathfrak{s}_{i}-\mathfrak{s}_{j}\|$. For
each component $k\in\{1,\dots,K\}$, draw i.i.d.\ Gaussian innovations
$\{z_{j,k}\}_{j=1}^{N}$ with $z_{j,k}\stackrel{\text{i.i.d.}}{\sim}\mathcal{N}(0,\sigma_{z}^{2})$,
independent across $j$ and $k$. Fix a range parameter $m>0$ and
a dependence parameter $\rho_{\mathfrak{d}}$, and define
the $m$-neighborhood $\mathcal{N}_{m}(i)\equiv\{j\le N:\mathfrak{d}_{ij}\le m\}$.
We then generate the spatial effect by the truncated, geometrically
weighted average
\begin{align}
\alpha_{i,k}\equiv\sum_{j\in\mathcal{N}_{m}(i)}w_{ij}\,z_{j,k},\qquad w_{ij}\equiv\rho_{\mathfrak{d}}^{\,\mathfrak{d}_{ij}}\mathbf{1}\{\mathfrak{d}_{ij}\le m\}.\label{eq:conley_molinari_dgp}
\end{align}
Because the weights are compactly supported, $\alpha_{i,k}$ and $\alpha_{j,k}$
are independent whenever $\mathfrak{d}_{ij}>2m$, so dependence decays
with distance and vanishes beyond a finite cutoff, matching the short-range
dependence implied by our mixing assumptions.
In our baseline implementation, we set $m=5$,
$\rho_{\mathfrak{d}}=0.10$, and $\sigma_{u}^{2}=1$. We further discuss the effect of varying levels of parameters in the Internet Appendix ID.
We set $\beta_{k}=1$ for all $k=1,\ldots,K$. For the number of regressors
($K$), MacKinnon (\citeyear{mackinnon2021fast}) suggests that performances
of many methods deteriorate with increasing $K$. Choosing a small
value of $K$, such as $K=2$ (a constant term and one regressor),
may yield an overly optimistic assessment. We hence choose $K=5$
and examine the true value of $\beta_{5}$ with various methods.
For functions $f_{1k}$ and $f_{2}$, we consider simulation designs
that cover all five scenarios:
\begin{flalign}
& \bullet\ \text{D, (\ref{eq: both diverge}) holds:} & {x_{it,k}=\alpha_{i,k}^{x}+\xi_{t,k}^{x}+\varepsilon_{it,k}^{x}\text{ and }u_{it}=\alpha_{i}^{u}+\xi_{t}^{u}+\varepsilon_{it}^{u},}\label{eq: dgp 1}\\
& \bullet\ \text{V\&N, (\ref{eq: vanish nongaussian}) holds:} & x_{it,k}=\alpha_{i,k}^{x}\xi_{t,k}^{x}\text{ and }u_{it}=\alpha_{i}^{u}\xi_{t}^{u},\label{eq: dgp 2}\\
& \bullet\ \text{V\&G, (\ref{eq: vanish gaussian}) holds:} & x_{it,k}=\varepsilon_{it,k}^{x}\text{ and }u_{it}=\varepsilon_{it}^{u},\label{eq: dgp 3}\\
& \bullet\ \text{I\&N, (\ref{eq: converge non-gaussian}) holds:} & {x_{it,k}=\left(\alpha_{i,k}^{x}+N^{-1/4}\right)\xi_{t,k}^{x}\text{ and }u_{it}=\left(\alpha_{i}^{u}+N^{-1/4}\right)\xi_{t}^{u},}\label{eq: dgp 4}\\
& \bullet\ \text{I\&G, (\ref{eq: converge gaussian}) holds:} & {x_{it,k}=\left(\varepsilon_{it,k}^{x}+N^{-1/4}\right)\xi_{t,k}^{x}\text{ and }u_{it}=\left(\varepsilon_{it}^{u}+N^{-1/4}\right)\xi_{t}^{u}.}
\label{eq: dgp 5}
\end{flalign}
\begin{figure}[t!]
\centering \begin{subfigure}[t]{0.49\hsize} \subcaption{DGP
(\ref{eq: dgp 1})}\includegraphics[width=1\textwidth]{spatial_dgp_D.png}
\end{subfigure}
\begin{subfigure}[t]{0.49\hsize} \subcaption{DGP
(\ref{eq: dgp 2})}\includegraphics[width=1\textwidth]{spatial_dgp_V_N.png}
\end{subfigure} \begin{subfigure}[t]{0.49\hsize} \subcaption{DGP
(\ref{eq: dgp 3})} \includegraphics[width=1\textwidth]{spatial_dgp_V_G.png}
\end{subfigure}
\begin{subfigure}[t]{0.49\hsize} \subcaption{DGP (\ref{eq: dgp 4})}\includegraphics[width=1\textwidth]{spatial_dgp_I_N.png}
\end{subfigure} \begin{subfigure}[t]{0.49\hsize} \subcaption{DGP
(\ref{eq: dgp 5})} \includegraphics[width=1\textwidth]{spatial_dgp_I_G.png}
\end{subfigure}
\caption{\textbf{Rejection Frequency for DGPs (\ref{eq: dgp 2})-(\ref{eq: dgp 5}).}
For each bootstrap
method, $B=999$. Results are based on 5,000 Monte Carlo replicates.
The predetermined significance level is 5\%.}
\label{fig: rej frequency DGPs 2-5}
\end{figure}
Figure~\ref{tab:average correct rate for DRC varying NT, ols} reports the
corresponding rejection frequencies. Panel~(a) presents the results under
DGP~\eqref{eq: dgp 1}, where there is strong cluster dependence along both
dimensions. As expected, the performance of all methods improves as the
number of clusters increases, which is consistent with the theoretical
prediction that all methods are valid in this regime.
Panel~(b) reports the results under DGP~\eqref{eq: dgp 2}, where the
limiting distribution is non-Gaussian. In line with the theory, PWB-V
does not provide valid inference, whereas the other methods perform well.
Notably, PWB-H performs relatively well even in small samples, although
this may partly reflect finite-sample randomness. Since PWB-H is a hybrid
of PWB-D and PWB-V, its finite-sample performance generally lies between
those of the two benchmark procedures.
Panel~(c) presents the results under DGP~\eqref{eq: dgp 3}, a setting in
which the dependence is primarily driven by intersection-level clustering.
All methods perform adequately in this scenario, again consistent with the
theoretical predictions.
Panel~(d) considers DGP~\eqref{eq: dgp 4}, the most challenging case, in
which the limiting distribution is non-Gaussian and no feasible procedure
can achieve asymptotic validity. The simulation results confirm this theoretical impossibility result, as
none of the methods exhibits further improvement when the number of
clusters becomes sufficiently large. Panel~(e) examines DGP~\eqref{eq: dgp 5}, where the asymptotic
distribution is Gaussian. In this setting, PWB-D is invalid and exhibits
substantial overrejection even when the sample size is large. By contrast,
PWB-V and PWB-H perform well as $N$ and $T$ increase.
Overall, the simulation results across the five scenarios corroborate the
theoretical findings. Among the feasible procedures, PWB-H delivers robust
performance across most scenarios, except in the infeasible regime where no
feasible method can be asymptotically valid. Therefore, in view of both the
theoretical analysis and the simulation evidence, we recommend PWB-H for
practical applications.
\begin{table}[t!]
\centering
\begin{tabular}{lccccccc}
\hline \hline
$N,T$ & 20 & 30 & 50 & 70 & 100 & 150 & 200 \tabularnewline
\hline
D, DGP \eqref{eq: dgp 1} & 0.991 & 0.998 & 0.999 & 0.999 & 0.999 & 0.999 & 0.999\tabularnewline
V\&N, DGP \eqref{eq: dgp 2} & 0.987 & 0.996 & 0.998 & 0.999 & 0.999 & 0.999 & 0.999 \tabularnewline
V\&G, DGP \eqref{eq: dgp 3} & 0.686 & 0.715 & 0.786 & 0.812 & 0.860 & 0.880 & 0.922\tabularnewline
I\&N, DGP \eqref{eq: dgp 4} & 0.972 & 0.990 & 0.990 & 0.991 & 0.997 & 0.999 & 0.999 \tabularnewline
I\&G, DGP \eqref{eq: dgp 5} & 0.538 & 0.631 & 0.735 & 0.810 & 0.855 & 0.933 & 0.964 \tabularnewline
\hline \hline
\end{tabular}\caption{\textbf{Classification accuracy of the dependence-regime classifier
(DRC) across varying \(N\) and \(T\).} Entries report the fraction of
replications in which the DRC selects the population regime; the optimal
classification probability approaches one over distinguishable regimes. For
the two uniformly indistinguishable non-Gaussian regimes, classification is
counted as correct when the DRC selects the non-Gaussian branch. For each bootstrap
method, $B=999$. Results are based on
5,000 Monte Carlo replicates.}
\label{tab:average correct rate for DRC varying NT, ols}
\end{table}
We also report simulation results for the dependence-regime classifier
(DRC) in Table~\ref{tab:average correct rate for DRC varying NT, ols}.
The entries in the table measure the accuracy of regime classification, not
the empirical size of the subsequent confidence interval; in the population
limit, the relevant classification probability is expected to converge to one
over the distinguishable regimes.
For the two indistinguishable non-Gaussian regimes, we record a classification
as correct whenever the procedure flags the component as non-Gaussian
(i.e., assigns it to either of the two non-Gaussian regimes). The
classification accuracy increases with the sample sizes $N$ and $T$,
and all scenarios exceed 90\% accuracy when $N=T=200$. When $N=T=20$, V\&N and I\&N demonstrate very high accuracy, indicating the KS diagnostic remains accurate in small samples as long as $B$ is sufficiently large. By contrast, I\&G and V\&G exhibit lower classification
accuracy because they are more likely to be confused when the variance
components are imprecisely estimated in small samples.
In the unreported results, we also consider
designs with a small numbers of factors ($K=2$), and
obtain similar results.
\section{Conclusion}
\label{sec:conclusion} This paper contributes to the econometrics literature on inference under
two-way clustering with serially and spatially dependent common effects. We
characterize the limiting distribution of the OLS estimator across five
mutually exclusive and exhaustive regimes, determined by the relative
contributions of the two clustering dimensions and the interaction component.
These regimes include both Gaussian and non-Gaussian limits and imply
different requirements for valid inference.
We show that one non-Gaussian regime is intrinsically infeasible: without
additional restrictions on the DGP, no procedure can achieve uniformly
consistent inference in that regime. We further establish two additional
impossibility results. First, the infeasible regime cannot be uniformly
distinguished from one feasible regime. Second, heterogeneous score components
under two-way clustering preclude uniformly consistent inference. Together,
these results identify fundamental limits on what can be learned from the
data in two-way clustered settings.
To address the feasible cases, we propose a family of projection-based wild
bootstrap procedures. The hybrid procedure, PWB-H, combines a data-driven
Gaussianity diagnostic, variance-scaling adjustments, and dependence-adaptive
bootstrap multipliers. It delivers uniformly valid inference across all four
feasible regimes, while the remaining non-Gaussian regime is shown to be
fundamentally beyond the reach of uniformly valid data-driven inference.
Monte Carlo simulations confirm that PWB-H performs well across a range of
dependence structures and sample sizes.
\section*{Appendix}