Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
85,369 characters · 14 sections · 58 citation commands
Bootstrap Inference under General Two-way Clustering with Serially and Spatially Dependent Common Effects
\thispagestyle{empty}
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 (bertrand2004much) and Petersen (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 (davezies2025analytic), we show that heterogeneous score components impose a fundamental limit on uniformly consistent inference under two-way clustering. Moreover, based on Menzel (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 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 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 (menzel2021bootstrap). Unlike Menzel (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 (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 menzel2021bootstrap and Juodis juodis2021shock, it does 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) presents the two-way clustering model and the five asymptotic regimes. Section (ref) introduces the PWB procedures and develops the theory for bootstrap validity. Section (ref) reports various simulation results under five regimes. Section (ref) concludes. Technical proofs and additional results are deferred to the appendix.
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
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
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
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). 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 (davezies2021empirical, davezies2022marcinkiewicz), MacKinnon, Nielsen, and Webb (mackinnon2021wild), and Menzel (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 (chiang2023standard), Chen and Vogelsang (chen2023fixed), Hounyo and Lin (hounyo2025jackknife, 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.
Following Conley 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 conley1999gmm; see also Conley and Molinari conley2007spatial; Kelejian and Prucha kelejian2007hac).
We further describe the dependence structure using an Aldous-Hoover-Kallenberg (AHK, Aldous aldous1981representations; Hoover hoover1979relations; Kallenberg kallenberg1989representation) representation, a standard device in the multiway clustering literature.
Assumption (ref) 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
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 (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:
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.
Assumption (ref) 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$.
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:
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 (ref) 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):
It corresponds to strong clustering along at least one dimension and generalizes settings such as Condition (16) of MacKinnon et al. 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.
\subparagraph{(V&N) Vanishing clustering but non-Gaussian interaction:}
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.
\subparagraph{(V&G) Vanishing clustering and Gaussian limit:}
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. 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}$.
\subparagraph{(I&N) Intermediate clustering with non-Gaussian interaction:}
All terms in (ref) 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)$.
\subparagraph{(I&G) Intermediate clustering with Gaussian limit:}
Here the interaction contamination vanishes, and a properly designed procedure can make the estimation error asymptotically negligible.
\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 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.
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.
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 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.\ (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} =
\] 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 (shao2010dependent) and Hounyo (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 (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:
where
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}$:
where the scaling terms $\bm{\vartheta}_{a}$ and $\bm{\vartheta}_{d}$ are define as follows:
and
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.
Assumptions (ref) (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.
Assumption (ref) (i) imposes a standard mixing condition from the time-series literature.\footnote{ It weakens the dependence restriction in Chiang et al.\ chiang2023standard by requiring an \(\alpha\)-mixing condition rather than a \(\beta\)-mixing condition. Chiang et al. 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) (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\}. \]
Assumption (ref) collects standard high-level regularity conditions used in the spatial dependence and spatial HAC literature, see, e.g., Conley 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.
Assumption (ref) requires the aggregate variances of different score components to be of the same order. Importantly, Assumption (ref) 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 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) 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.
Proposition (ref) 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)). 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), 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}\}$.
Theorem (ref) 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}$.
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}$:
and
where $EVC\left(\cdot\right)$ replaces any negative eigenvalue by zero to ensure positive semi-definiteness.
Lemma (ref) 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
where
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 (ref) and (ref), 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 (ref) holds. In practice, we recommend setting
\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 (ref) holds. In practice, we recommend
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). 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}\}. \]
Theorem (ref)(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)(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 (ref) and (ref). 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), 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 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) 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 leeb2008can. The simulation results in Table (ref) 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 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), PWB-H is uniformly asymptotically exact over all regimes except I&N, i.e., when (ref) holds, where uniform consistency is ruled out by Proposition (ref)(a). This supports PWB-H as a general-purpose inference procedure.
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)(c) indicates a fundamental tension: any method that is uniformly valid in all feasible regimes 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.
Figure (ref) (visualizing Proposition (ref)) 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.
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\} . \]
Finally, note that the two non-Gaussian regimes cannot be distinguished by DRC. One might hope to construct a sharper classifier, but Proposition (ref)(b) shows that, absent additional information on the DGP, no procedure can uniformly distinguish these two regimes.
Table (ref) 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 (juodis2021shock), CHS refers to the variance estimator proposed by Chiang et al. (chiang2023standard), and MWCB stands for the multiway cluster bootstrap developed by Hounyo and Lin (hounyo2024wild).
Several inference methods have been proposed under the assumption of no temporal dependence. For example, the variance estimator of Cameron, Gelbach, and Miller (cameron2011robust), the wild bootstrap procedures of MacKinnon et al. (mackinnon2021wild), and several bootstraps of Menzel (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. (cameron2011robust) estimator and MacKinnon et al. (mackinnon2021wild) wild bootstrap methods are similar to the CHS variance estimator and MWCB, respectively, and Menzel (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 menzel2021bootstrap and Juodis juodis2021shock. We summarize several key differences between these approaches and our PWB framework (under no temporal dependence):
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
However, this approach is essentially perturbing residuals, as the first equality in ((ref)) 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}$.
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 fernandez2021low and Chen, Fern{\'a}ndez-Val, Weidner 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:
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:
Such an AR(1) process with $\rho<1$ satisfies Assumption (ref).
The latent components $\left(\alpha_{i,k}^{x},\alpha_{i}^{u}\right)$ are generated by adopting the spatial design of Conley and Molinari 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
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 (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:
Figure (ref) reports the corresponding rejection frequencies. Panel (a) presents the results under DGP (ref), 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 (ref), 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 (ref), 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 (ref), 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 (ref), 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.
We also report simulation results for the dependence-regime classifier (DRC) in Table (ref). 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.
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.