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.
53,474 characters · 12 sections · 48 citation commands
Genuinely Robust Inference for Clustered Data
\allowdisplaybreaks \onehalfspacing
{7.45mm} \sloppy
Cluster-robust (CR) standard errors are designed to account for within-cluster correlations. Such correlations often arise by construction, for example, within an industry He1998 or within a state BeDuMu2004. Today, even when a model does not inherently induce cluster dependence, the application of CR methods using observable group identifiers has become a common practice.
The foundational theory \citep*[][]{Wh84,LiZe86,Ar87} for CR inference methods assumes small cluster sizes $N_g$ (uniformly bounded above by $\overline{N} < \infty$) with a large number of clusters, $G \rightarrow \infty$, where $N_g$ denotes the number of entities in the $g$-th cluster for $g \in \{1,2,\dots,G\}$. Procedures based on this theory are implemented through the `cluster()' and `vce(cluster)' options in Stata, and they are utilized in nearly all, if not all, empirical studies that report CR standard errors.
It has been recognized that large cluster sizes $N_g$ can result in inflated CR standard errors \citep*[e.g.,][p. 324]{CaMi15}. Recent theoretical advancements carter2017asymptotic,DjMaNi19,HaLe19,Ha22,BuCaShTa2022 accommodate larger cluster sizes $N_g$, eliminating the requirement that $N_g \leq \overline{N}$ and thereby broadening the applicability of the `cluster()' and `vce(cluster)' options, among others. With this said, they still impose the restriction $\max_g N^2_g / N \rightarrow 0$ of vanishing maximum cluster size relative to the square root of the whole sample size $N = \sum_{g=1}^G N_g$ as $G \rightarrow \infty$.
A natural question is whether the relaxed condition $\max_g N_g^2 / N \rightarrow 0$ accommodates a wide range of data sets. To answer this, we analyze empirical papers published in top journals.\footnote{We studied all articles published in American Economic Review and Econometrica between 2020 and 2021. Among them, we extracted a list of papers reporting estimation and inference results based on regressions, IV regressions, and their variants. Furthermore, we focus on articles using publicly available data sets for replication. See Section (ref) for further details of this study.} All of these articles employ the aforementioned Stata options for CR standard errors, thereby implicitly assuming $\max_g N_g^2 / N \rightarrow 0$. Table (ref) summarizes the number of articles with $\max_g N_g^2 / N$ falling into each bin on a logarithmic scale. Notably, 55 percent (respectively, 39, 29, and 16 percent) of the articles use data sets where $\max_g N_g^2 / N \geq 1$ (respectively, $\geq 10$, $\geq 100$, and $\geq 1000$). In other words, the condition $\max_g N_g^2 / N \rightarrow 0$, required for the validity of conventional CR inference, may not hold for a nontrivial portion of these published articles.
The condition $\max_g N_g^2 / N \to 0$ is sufficient but not necessary for asymptotic normality, implying that the adequacy of normality-based confidence intervals and tests cannot be evaluated solely by assessing the plausibility of this condition. To address this, we establish a necessary and sufficient condition for the validity of conventional cluster-robust (CR) inference -- see Theorem (ref). Specifically, the limiting distribution is normal if and only if the score of the largest cluster is ignorable. When clusters are unignorably large, regression estimates exhibit non-Gaussian limiting distributions, as illustrated in Figure (ref).\footnote{Details on these non-Gaussian distributions are provided in Section (ref).} Using this characterization, formal statistical tests based on sasaki2023diagnostic reject the null hypothesis of normality in 24 of the 31 papers (77 percent) reported in Table (ref)--see Table (ref).
Non-Gaussian limiting distributions invalidate conventional critical values, such as “1.96,” as well as bootstrap critical values. For example, using 1.96 results in sizes of 0.053, 0.087, and 0.250 (instead of the desired 0.050) when the nuisance parameter $\alpha$ equals 1.75, 1.50, and 1.25, respectively, as shown in Figure (ref). The empirical bootstrap fails in these cases of infinite variance, and the widely used wild cluster bootstrap and pairs cluster bootstrap are also inconsistent. Later, we formally establish these negative results as Proposition (ref).
To address this issue, we introduce the cluster score (CS) bootstrap, a novel inferential method for clustered data that extends the $m$-out-of-$n$ and score bootstraps. This method provides valid critical values adaptively across all limiting distributions depicted in Figure (ref). In addition, we provide a data-driven choice of the tuning parameter and justify its theoretical validity.
{\bf Relation to the Literature:} The literature on cluster-robust inference has a long history, dating back to Wh84, LiZe86, and Ar87. For a comprehensive review, we refer readers to CaMi15,cameronreview and mackinnon2023cluster. More recently, sampling frameworks in which cluster sizes are treated as random variables have been investigated by bai2022inference, BuCaShTa2022, and cavaliere2022econometrics. We adopt a model-based perspective with an increasing number of clusters and unrestricted intra-cluster dependence, or, asymptotically equivalently, a sampling-based perspective where the growing number of sampled clusters represents a negligible fraction of the superpopulation. This framework is well suited to the empirical contexts encountered in most applications.\footnote{ An alternative framework assumes a fixed number of clusters with growing cluster sizes, where asymptotic normality can be derived under additional assumptions of weak intra-cluster dependence, as in Ha07, IbMu10, canay2021wild. ibragimov2016inference and hansen2022jackknife consider inference under gaussian assumptions. Another line of research advances the integration of design-based and sampling-based asymptotics, particularly under explicit treatment assignment schemes such as randomized experiments, as considered by AbAtImWo23. Extending our method to this design-based framework is an avenue for future research. }
In an insightful recent work, kojevnikov2021some establish an impossibility result for consistent estimation of the asymptotic variance when the sample contains a single large cluster under a triangular array setup. They further provide a necessary and sufficient condition on the cluster structure for the asymptotic variance to be consistently estimable. Our findings complement their result by showing that normal approximation for $t$-statistics fails in the presence of unignorably large clusters. Furthermore, our proposed procedure overcomes this limitation, as it does not rely on consistent variance estimation. We demonstrate that the self-normalized statistic converges in distribution and formally derive its limiting stable distribution in such settings. Importantly, the implementation of the bootstrap inference procedure does not require knowledge of the unknown rate or consistent variance estimation, owing to the self-normalizing nature of the test statistics.
The aspect of our paper that establishes non-Gaussianity under certain conditions connects to a branch of the econometrics literature exploring the possibility of non-Gaussian limiting distributions and, in some cases, impossibility results. For instance, hirano2012impossibility show that regular estimation, and hence Gaussian asymptotics, are impossible for a class of estimators characterized by the maximum. More directly related, menzel2021bootstrap highlight the potential non-Gaussianity of estimators under two-way clustering and establish an impossibility result on the uniform consistency of tests. In contrast, we demonstrate that non-Gaussianity can arise even under one-way clustering. Moreover, such negative results may be even more pervasive in this widely used empirical setting.
Our key distributional approximation results build on logan1973limit, lepage1981convergence, and gine1997student. For theoretical foundations of probability and statistics with heavy-tailed distributions, we refer readers to resnick1987extreme,resnick2008extreme and samorodnitsky1994stable. Also see ibragimov2015heavy,chernozhukov2016extremal,chernozhukov2017extremal for applications in economics and finance. Our inference procedure builds on the resampling theory developed in arcones1989bootstrap,arcones1991additions and bickel2008choice. For the related discussions on the inconsistency of empirical bootstrap for means of random variables with infinite variance, see, e.g., athreya1987bootstrap and knight1989bootstrap.
While the idea extends to a general class of econometric models, we consider the linear model\footnote{The assumption $\mathbb{E}[U_{gi}|X_g]=0$, while standard in the literature, is stronger than required for our asymptotic results. It can be relaxed to $\mathbb{E}[\sum_{i=1}^{N_g}X_{gi}U_{gi}]=0$.}
for ease of exposition as well as its popular use in practice, where $X_g = (X_{g1},\ldots,X_{gN_g})'$, $U_g = (U_{g1},\ldots,U_{gN_g})'$, $g \in \{1,\ldots,G\}$ indexes clusters, and $N_g$ denotes the size of the $g$-th cluster. Define the OLS estimator and its cluster-robust (CR) variance estimator by
respectively, for some finite sample adjustment factor $a_G$ such that $a_G \to 1$ as $G \to \infty$, where $S_g = \sum_{i=1}^{N_g} X_{gi} U_{gi}$, $\widehat S_g = \sum_{i=1}^{N_g} X_{gi} \widehat U_{gi}$, and $\widehat U_{gi}=Y_{gi}-X_{gi}'\widehat\theta$. For simplicity of writing, we set $a_G=1$ throughout as it does not affect our asymptotic arguments.
Consider a linear transformation $\delta=r'\theta$ of the regression coefficient vector $\theta$, such that $r\in \mathbb{R}^{\dim(\theta)}$ and $\|r\|=1$, as the parameter of interest. Let the corresponding estimator and its CR standard error be denoted by
respectively. We are interested in conducting inference for $\delta$ using the t-statistic
based on the CR standard error.
To state our assumption, we introduce a few definitions. A random variable $\eta$ is said to be stable if it has a domain of attraction in that there exists a sequence of i.i.d. random variables $\xi_1,\xi_2,\ldots$ and sequences of positive numbers $A_G$ and real numbers $D_G$ such that
A function $L(\cdot)$ is said to be slowly varying at $\infty$ if $\lim_{t\to \infty}L(yt)/L(t)=1$ for all $y>0$. If $\eta$ is stable, then $A_G$ takes the form of $G^{1/\alpha}L(G)$ for some $\alpha\in (0,2]$ and some slowly varying function $L(\cdot)$ at $\infty$ (cf. Proposition 2.2.13 in embrechts1997modelling). If $\alpha \in (1,2]$, then $D_G$ can be chosen to be $G\cdot \mathbb{E}[\xi_g]$. The number $\alpha$ is called the index of stability, and $\eta$ is said to be $\alpha$-stable. In such a case, $\xi_g$ is said to belong to the domain of attraction of an $\alpha$-stable distribution. Although this concept may look esoteric to some readers, it essentially states that a sum of i.i.d. random variables, after being suitably centered and normalized, converges in distribution to a limiting random variable, and it, in particular, encompasses the standard cases where central limit theorems (CLTs) hold. In other words, econometricians and economists adopting the standard inference (e.g., the conventional critical value of 1.96) implicitly make this (and even stronger) assumption.
This assumption is arguably general. It is significantly weaker than requiring a central limit theorem to hold and even covers scenarios where the asymptotic normality fails. It also encompasses two notable cases frequently considered in economics and econometrics.
First, the case of $\alpha = 2$ encompasses the conventional setting in which $r'(\widehat\theta - \theta)$ attains the standard convergence rate of $\sqrt{G}$ by the central limit theorem. In this case, the limiting $\alpha$-stable distribution is necessarily Gaussian Geluk2000. This scenario also includes certain non-standard cases with a Gaussian limiting distribution but without finite variance, such as a Pareto random variable with Pareto exponent equal to $2$. The vast majority of econometric papers deriving asymptotic normality implicitly rely on this high-level condition in Assumption (ref), or on even stronger ones.
Second, the case of $\alpha < 2$ entails a power law pena2009self, i.e.,
for some slowly varying functions $L_1(\cdot)$ and $L_2(\cdot)$, where $L_2(\cdot)$ may depend on $u_1$ and $u_2$. In this case, the index $\alpha$ of stability coincides with the Pareto exponent\footnote{Specifically, the Pareto distribution has CDF $F(t) = 1 - t^{-\beta}$ for $t \geq 1$.} $\beta$, in the sense that $\alpha = \min\{\beta,2\}$. Thus, when $\alpha < 2$, the score has infinite variance. For more precise details, see Theorem (ref) in Appendix (ref). In this case, unignorably large clusters are indeed unignorable, since the sample sum of the (scaled) scores becomes asymptotically proportional to the (scaled) score of the largest cluster; see Remark (ref) in Appendix (ref) for further discussion. Hence, the asymptotic distribution cannot be Gaussian in this case.
The literature in urban economics and economic geography establishes that (truncated) city size distributions frequently exhibit a power-law behavior in the upper tail, with estimated exponents typically in the range of $1$ to $1.5$; see eeckhout2004gibrat,ioannides2013us and references therein. Consequently, when observations are strongly correlated at the city level, this implies $\alpha < 2$, yielding a non-Gaussian limiting distribution.
The i.i.d. requirement across clusters in Assumption (ref) is standard in this literature (e.g., BuCaShTa2022,cavaliere2022econometrics,bai2022inference). This assumption is mild because: (1) the conditional distributions of $S_g$ and $X'_g X_g$ given $N_g = n_g$ may vary across $n_g$; and (2) the distributions of individuals within each cluster need not be identical. Moreover, $S_g$ and $X_g$ may be arbitrarily correlated with the cluster size $N_g$, provided that the regression exogeneity condition is satisfied.
To simplify the exposition, we focus on the case where $v'S_g$ and $u_1' X_g' X_g u_2$ share a common stability index $\alpha$. This simplification is rationalized if the tail behavior of their distributions is driven by the tail behavior of the cluster-size distribution $N_g$; see Section (ref) for an illustrative example. However, this simplification is adopted only for notational simplicity and can be relaxed at the cost of substantially more cumbersome exposition.
This section shows that conventional methods of cluster-robust (CR) inference are valid if and only if $\alpha = 2$. In other words, they necessarily fail when $\alpha < 2$. We begin with heuristic discussions in Section (ref) and then develop formal results in Section (ref). We further report how frequently cases with $\alpha < 2$ arise in empirical economic applications.
The intuition behind the fragility of conventional CR methods is as follows. When $\alpha < 2$, the cluster size $N_g$ does not have a finite variance. If intra-cluster dependence is non-trivial, this infinite variance of $N_g$ is inherited by the score $S_g$, causing the CLT for OLS (and also other estimators) to fail. Cases with $\alpha < 2$ are quite plausible in empirical data, and we show that this is indeed the case for the majority of recent empirical papers published in Econometrica and the American Economic Review, as discussed in Section (ref) in detail.
To provide a simple and transparent illustration, consider the sample average
which is a special case of the OLS estimator (ref) with $X_{gi} = 1$. The true parameter is the mean $\theta = \mathbb{E}[Y_{gi}]$, which is normalized to $\theta = 0$ without loss of generality. For clarity, suppose the extreme case of perfect intra-cluster dependence, i.e., $Y_{gi} \equiv Y_{g}$ for all $i \in \{1, \ldots, N_g\}$ within each cluster $g$. For simplicity, also assume that $N_g$ is independent of $Y_g$. These assumptions are made purely for expositional clarity and are not essential for our results.
In this case, we obtain
The denominator converges to $\sqrt{\mathbb{E}[N_{g}]}$ as long as $\alpha > 1$. For the numerator, note that $ \mathrm{Var}\!\left(G^{-1/2}\sum_{g=1}^{G} N_{g} Y_{g}\right) = \mathrm{Var}[N_g] \cdot \mathrm{Var}[Y_g], $ which is infinite when $\alpha < 2$. Indeed, Theorem 1 of Geluk2000 implies that, if the distribution of $N_g Y_g$ is $\alpha$-stable, then the limiting distribution
for some sequences $a_G \simeq G^{1/\alpha} \to \infty$ and $b_G \in \mathbb{R}$, has the characteristic function
which differs from the Gaussian characteristic function. Thus, the CLT for $\widehat{\theta}$ fails. Further discussion of this example is provided in Appendix (ref).
In summary, the stability index $\alpha$ determines both the convergence rate and the limiting distribution. When intra-cluster correlation is non-trivial, the tail heaviness of $N_g$ carries over to that of $S_g$. Consequently, the $t$-ratio of the conventional CR method is asymptotically normal if and only if $\alpha = 2$. While we currently present heuristic arguments in a simplified setting, we formalize and generalize this claim in Theorem (ref) and Proposition (ref) below.
The following theorem formalizes and generalizes the discussion from the previous subsection.
A proof is provided in Appendix (ref).
The theorem implies that conventional CR inference based on common variance estimators, such as CR1, CR2, CR3, and the jackknife, together with normal critical values (e.g., $\approx 1.96$ for the 97.5th percentile) fails whenever $\alpha < 2$.
With Theorem (ref), we now characterize the curves shown in Figure (ref) from Section (ref). The left, middle, and right panels of Figure (ref) display the limiting distributions of the $t$-statistic under $p = 0.25$, $0.50$, and $0.75$, respectively, where
for $v$ as given in Assumption (ref), measures the limiting asymmetry of tail probabilities. Each panel in Figure (ref) depicts three non-Gaussian limiting distributions corresponding to $\alpha = 1.25$, $1.50$, and $1.75$ with distinct line styles, together with the normal reference case ($\alpha = 2.00$). The key takeaway is that conventional methods of CR inference, which rely on normal approximation, become increasingly size-distorted as $\alpha$ decreases and as $p$ deviates from $0.5$.
Another class of conventional approaches consists of cluster bootstraps. Two main bootstrap-based CR inference methods are commonly used in the literature: the pairs cluster bootstrap and the wild cluster bootstrap cameron2008bootstrap. It is well established that the empirical bootstrap is inconsistent when the variance of the score is infinite athreya1987bootstrap,knight1989bootstrap. In light of the power-law characterization (ref), the pairs cluster bootstrap, essentially the empirical bootstrap applied to cluster-wise sums treated as independent units, is inconsistent under Assumption (ref) with $\alpha < 2$. Moreover, Theorem (ref) in Appendix (ref) shows that the wild cluster bootstrap is likewise inconsistent under Assumption (ref) with $\alpha < 2$. The following proposition summarizes these results.
Given that the case of $\alpha < 2$ invalidates all conventional methods of CR inference, a natural question is how frequently such cases arise in empirical economics. To address this, we examined all articles published in two leading journals, the American Economic Review and Econometrica, during 2020–2021. From these, we extracted the subset of papers reporting estimation and inference results based on regressions, IV regressions, and related variants. We further restricted attention to articles that employ publicly available datasets due to replicability.
For these articles, we test the null hypothesis $H_0: \alpha = 2$ against the alternative $H_1: \alpha < 2$ for the score. Such a test can be implemented via the likelihood ratio test of sasaki2023diagnostic, which considers the surrogate null hypothesis $H_0: \beta \ge 2$ against the alternative $H_1: \beta < 2$ in light of (ref), where $\beta$ denotes the tail exponent of the score.\footnote{The test of the null hypothesis $H_0: \beta \ge 2$ against the alternative $H_1: \beta < 2$ is implemented using the Stata command testout y x1 x2 ..., cluster(cid) for least-squares estimation and testout y x1 x2 ..., iv(z) cluster(cid) for instrumental variables estimation, both following sasaki2023diagnostic.}
Table (ref) summarizes the set of papers included in our study. The first two columns report the journals and years of publication. The next column, “All \#,” indicates the total number of eligible articles according to the selection criteria described above. The column group labeled “Cluster” contains articles in which CR inference is applied to at least one regression result. Within this group, the column “\#” reports the number of such articles, while the column “Test $\alpha < 2$” reports the fraction of these articles for which the test rejects the null hypothesis in one or more regression specifications. The final row presents the column totals.
During 2020--2021, the American Economic Review published 30 articles that met our selection criteria. Of these, 21 reported CR standard errors. The null hypothesis is rejected in 16 of these 21 papers. In other words, inference based on the conventional CR method may be misleading in approximately $76\%$ of the articles employing it.
During 2020--2021, Econometrica published 14 articles that met our selection criteria. Of these, 10 reported CR standard errors. The null hypothesis is rejected in 8 of these 10 papers. In other words, inference based on the conventional CR method may be misleading in approximately $80\%$ of the articles employing it.
Combining the two journals, we find that inference may be misleading in as many as $77\%$ of the 31 articles that employ the conventional CR method. Thus, problematic practice appears to be prevalent even in these highly influential outlets.\footnote{Spreadsheets of all the test results with specific papers and specific equations are available upon request.} All of the above issues with conventional CR methods motivate our proposed approach: the cluster score (CS) bootstrap, which accommodates non-Gaussian limiting distributions, to be presented in Section (ref).
In light of the limitations of conventional CR inference methods discussed in the previous section, we introduce a novel cluster score (CS) bootstrap procedure to approximate the limiting distribution of $(\widehat{\delta} - \delta)/\widehat{\sigma}$. This procedure remains valid whether the limiting distribution is Gaussian or non-Gaussian. Section (ref) describes the proposed method, and Section (ref) provides its theoretical justification.
Our objective is to conduct statistical inference for $\delta$ using the $t$-statistic defined in (ref). Let the CDF $J_G^\ast$ of the $t$-statistic be
We will show that, under suitable regularity conditions, $J_G^*$ converges to the CDF $J^*$ of the corresponding limiting distribution.
Let $b$ denote the number of resampled clusters, chosen according to Algorithm (ref) (to be presented in Section (ref)). For a large positive integer \(M\), draw \(M\) i.i.d. multinomial random vectors \((w_1^j,\ldots,w_G^j)_{j=1}^M\), each with \(b\) trials and uniform cell probability \(1/G\), independently of the data. This is equivalent to sampling \(b\) clusters with replacement uniformly from the \(G\) clusters in the data, and \(w_g^j\) records the counts of how many times the $g$-th cluster is selected in the $j$-th bootstrap sample. Define the CS bootstrap estimator and its associated variance estimator by
where $ \widehat{S}_{g,j} = X_g' \big( Y_g - X_g \widehat{\theta}_{b,j} \big). $
Note that the inverse factor $\left(\sum_{g=1}^G X_g'X_g \right)^{-1}$ is computed from the full unweighted sample, whereas the linear component and its variance are constructed from the bootstrap sample. Practical motivations for this feature will be discussed in Remark (ref).
Define the bootstrapped empirical distribution function $\widehat{L}_{G,b}$ of $(\widehat{\delta}_{b,j} - \widehat{\delta})/\widehat{\sigma}_{b,j}$ by
For any $a \in (0,1)$, define the corresponding critical value as
As will be formally established in Section (ref), this critical value is guaranteed to satisfy
as $G \to \infty$. Hence, the CS bootstrap method provides asymptotically valid inference.
{\bf Practical Implication:} For the $t$-statistic, one may continue to use the conventional CR “standard error” $\widehat{\sigma}$ even though it may diverge.\footnote{Note that the “standard error” $\widehat{\sigma}$ does not converge in probability when $\alpha < 2$.} However, rather than relying on conventional Gaussian critical values (e.g., $\Phi^{-1}(0.025) \approx -1.96$ and $\Phi^{-1}(0.975) \approx 1.96$), one should instead employ the bootstrap-based critical values $\widehat{c}_{G,b}(0.025)$ and $\widehat{c}_{G,b}(0.975)$ obtained from the CS bootstrap procedure. These critical values, for example, can be used to construct a 95% confidence interval for $\delta$. $\blacktriangle$
We now provide theoretical support for our proposed CS bootstrap method. The following theorem provides a formal justification of its robust asymptotic validity.
The proof is non-trivial and proceeds by considering two distinct cases. The first, corresponding to $\alpha < 2$, is formalized in Lemma (ref) in Appendix (ref) and proved in detail in Appendix (ref). The second case, $\alpha = 2$, is analyzed in Appendix (ref), which also synthesizes the two regimes to establish Theorem (ref). Here, asymptotics are taken with respect to $G \to \infty$ for a given DGP; further studies on uniformity can be found in Appendix (ref).
While the theory requires $b \to \infty$ and $b/G = o(1)$ as $G \to \infty$, in practice the researcher must select a finite value of $b$. This choice should be neither too large nor too small. Intuitively, if $b$ is chosen too close to $G$, the largest clusters are sampled too frequently in the bootstrapped $t$-statistics, which prevents the procedure from adequately reflecting the heavy-tailed nature of the DGP when $\alpha < 2$. Conversely, if $b$ is too small, the bootstrapped $t$-statistics become excessively noisy. Thus, one seeks a value of $b$ that lies in a stable range, such that small perturbations of $b$ (e.g., increasing or decreasing it by one) have only minimal impact on the bootstrap distribution. In the context of the $m$-out-of-$n$ bootstrap, bickel2008choice formalized this idea and proposed a data-driven algorithm with theoretical guarantees for its validity. Here, we introduce a modified version of their algorithm tailored to our setting.
The following theorem provides theoretical guarantees for the data-driven choice of $\widehat{b}$ in the CS bootstrap. Specifically, it shows that Theorem (ref) continues to hold with our data-driven choice of $\widehat{b}$.
A proof is found in Appendix (ref). Hence, following the selection of $\widehat b$, asymptotically valid inference can be carried out using the corresponding critical value $\widehat c_{G,\widehat b}(1-a)$.
In this section, we present simulation studies evaluating the finite-sample performance of our proposed method of genuinely robust CR inference, based on the cluster score (CS) bootstrap, in comparison with conventional CR methods.
The data-generating design is defined as follows. We consider the cluster treatment model with individual covariates $$ Y_{gi} = \theta_0 + \theta_1 T_g + \sum_{j=1}^K \theta_j X_{g,i,j+1} + U_{gi}, $$ following \citet*[][Equation (40)]{MaNiWe2022fast}, among others. The binary treatment variable $T_g$ equals one for $\lceil 0.2G \rceil$ clusters and zero for the remaining $G - \lceil 0.2G \rceil$ clusters, where $\lceil a \rceil$ denotes the smallest integer greater than or equal to $a$. Cluster sizes are drawn independently as $N_g \sim \lceil \text{Pareto}(1,\alpha) \rceil$ for $g \in \{1,\ldots,G\}$. For each $g \in \{1,\ldots,G\}$, we independently draw $N_g$-variate random vectors $(\widetilde X_{g1j},\ldots,\widetilde X_{gN_gj})' \sim \mathcal{N}(0,\Omega)$ for $j \in \{1,\ldots,K\}$ and $(\widetilde U_{g1},\ldots,\widetilde U_{gN_g})' \sim \mathcal{N}(0,\Omega)$ in the baseline design, where $\Omega$ is an $N_g \times N_g$ covariance matrix with $\Omega_{ii}=1$ for all $i \in \{1,\ldots,N_g\}$ and $\Omega_{ii'}=1/2$ whenever $i \neq i'$. The controls are constructed as $X_{gij} = 0.2 F_{\text{Beta}(2,2)}^{-1} \circ \Phi(\widetilde X_{gij})$, where $F_{\text{Beta}(2,2)}$ and $\Phi$ denote the CDFs of the $\text{Beta}(2,2)$ and standard normal distributions, respectively. The errors are constructed heteroskedastically as $U_{gi} = 0.2 \widetilde U_{gi}$ if $T_g=0$ and $U_{gi} = \widetilde U_{gi}$ if $T_g=1$.
We vary the exponent parameter $\alpha \in \{1.1,1.2,\ldots,1.9,2.0\}$ across simulation sets. The regression coefficients are fixed at $(\theta_0,\theta_1,\theta_2,\ldots,\theta_{K+1})' = (1,1,1,\ldots,1)'$ throughout, while the covariate dimension varies as $K \in \{5,10\}$. The sample size (i.e., the number of clusters) is set to $G = 50$ across all simulations, which is roughly comparable to the number of U.S.\ states. Each simulation set consists of 10,000 Monte Carlo iterations.
Figure (ref) reports the Monte Carlo coverage frequencies. The horizontal axis represents the value of $\alpha$, and the vertical axis represents the coverage frequency. In the legend, `CSB' denotes the CS bootstrap, while `WCB' and `CR1' denote the wild cluster bootstrap and the CR1 standard error with normal critical values, respectively. The nominal coverage probability of 95% is indicated by the horizontal gray line at 0.95.
The CS bootstrap performs best, followed by the WCB and the CR1. Overall, the CS bootstrap consistently delivers coverage frequencies closest to the nominal 95% level across the range of $\alpha$. By contrast, both conventional methods suffer from under-coverage, particularly for small values of $\alpha$.
akhtari2022political study the effects of political turnover on various outcomes measuring the quality of public services in Brazil. In their original paper (Table 3, Column 5), they estimate the following linear model by OLS:
The dependent variable, $\text{score}_{gi+1}$, is the test score of fourth-grade students in the year following an election. The main explanatory variable, $\text{IVM}_{g}$, is the incumbent vote margin. Thus, $\mathds{1}\{\text{IVM}_{g}<0\}$ equals one when the incumbent party loses the election. The parameter of interest is $\delta = \theta_1$, which measures the effect of political turnover on test scores.\footnote{This effectively implements a sharp regression discontinuity design, although the original paper estimates the effect by OLS using this linear specification.} While the original paper considers alternative `bandwidths,' we focus on the bandwidth 0.110 to maximize the sample size, following prior work MaNiWe2022fast, which replicates this regression.
The original paper clusters standard errors at the municipality level, and we follow this definition of the cluster unit. There are $G = 2101$ municipalities in the data, with $\max_{1 \leq g \leq G} N_g^2 / N \approx 26$. Thus, the assumption $\max_{1 \leq g \leq G} N_g^2 / N \to 0$, under which conventional CR inference methods are guaranteed to work, is difficult to justify in this application.
Table (ref) reports the $p$-values for $\delta = \theta_1$ based on alternative inference methods. Column HC1 reports the $p$-value using conventional inference without clustering, i.e., the HC1 standard error. Columns CR1 and WCB report the $p$-values using conventional CR inference methods, namely the CR1 standard error and the wild cluster bootstrap,\footnote{While there are four variants of WCB MaNiWe2022fast, they yield identical $p$-values up to the reported digits, and hence we summarize them in a single column.} respectively, with normal approximation. Finally, column CSB reports the $p$-value based on our proposed inference method using the CS bootstrap. We employ the same code as in the simulation studies in Section (ref), including the choice of $b$ based on the minimum volatility method.
The $p$-value is zero up to the third digit when standard errors are not clustered (HC1). Conventional CR inference methods (CR1 and WCB) with normal approximation yield larger $p$-values, but the statistical significance remains unchanged. By contrast, our proposed CS bootstrap method produces a much larger $p$-value, rendering the effect $\delta = \theta_1$ statistically insignificant, unlike any of the conventional methods. These results highlight that failing to account for potential non-Gaussianity in the limiting distributions, particularly in the presence of unignorably large clusters, can lead to erroneous statistical conclusions.
Conventional methods for cluster-robust inference often fail to provide consistent results in the presence of unignorably large clusters. In this paper, we formalize this limitation by deriving a necessary and sufficient condition for consistency. We document that 77% of empirical research articles published in the American Economic Review and Econometrica during 2020--2021 contain model specifications fail to satisfy this condition.
To address this challenge, we propose the CS bootstrap and establish its size control across a wide class of data-generating processes where conventional methods break down. Our simulation studies confirm the reliability and effectiveness of the proposed method, underscoring its practical value in overcoming the limitations of existing cluster-robust inference techniques. We further demonstrate the failure of the wild cluster bootstrap in Section (ref) and discuss the related uniformity issues in Section (ref) in the appendix, reinforcing the need for our proposed approach. Finally, we demonstrate that correctly accounting for potential non-Gaussianity can overturn empirical conclusions. We conclude the paper by modifying a well-known haiku by hiranohaiku.