Refined Cluster Robust Inference

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.

29,408 characters · 7 sections · 38 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstract\linespread{1.2} It has become standard for empirical studies to conduct inference robust to cluster dependence and heterogeneity. With a small number of clusters, the normal approximation for the $t$-statistics of regression coefficients may be poor. This paper tackles this problem using a critical value based on the conditional Cram\'er-Edgeworth expansion for the $t$-statistics. Our approach guarantees third-order refinement, regardless of whether a regressor is discrete or not, and, unlike the cluster pairs bootstrap, avoids resampling data. Simulations show that our proposal can make a difference in size control with as few as 10 clusters. \begin{description} • Cluster robust inference, Cram\'er-Edgeworth expansion, Asymptotic refinement • C12, C21 \end{description}

Introduction

Cluster robust inference has become a standard practice in applied microeconometrics. Robustness to arbitrary dependence within a cluster often comes at the cost of a reduction in the effective sample size. As a result, despite a large overall sample size, the researchers have to account for the non-Gaussian distribution of the $t$-statistic for significance tests as was the case in classical statistics student1908probable. In particular, it is well known that the standard normality-based inference method may lead to over-rejection when the number of clusters is small cameron2015practitioner,comeron/miller:2026,mackinnon2023cluster. Instead of assuming Gaussian regression errors, one can use only mild moment restrictions and derive corrected critical values using higher order asymptotic expansions of Edgeworth1905 and Cramer1928.

In this paper, we propose a new analytical correction to critical values based on inverting the Cram\'er-Edgeworth expansion of the $t$-statistic null distribution. The resulting inference method is third-order asymptotically accurate and robust against heterogeneous cluster dependence. We first derive the Cram\'er-Edgeworth term up to the second order and then adjust the critical value based on the estimated Cram\'er-Edgeworth term.

We develop the approach of hall1983inverting in the setup of linear regression with heterogeneous clusters.\footnote{In this paper, we consider a two-sided alternative hypothesis and the second-order Cram\'er-Edgeworth expansion. As a result, we do not need to apply the Cram\'er-Edgeworth expansion of the estimated coefficient recursively, as in hall1983inverting.} As special sub-cases, our approach also nests inference on the sample mean of non-identically distributed data and regression coefficients in cross-sectional regression. The resulting inference is third-order asymptotically accurate in the sense that the actual test size differs from the prespecified significance level $\alpha$ only by $o(G^{-1})$ with $G$ clusters.

Our simulation studies support this theoretical size control.

A cluster pairs bootstrap might be a popular choice for asymptotic refinement, but our proposed method has a few advantages over the cluster pairs bootstrap.\footnote{The residual bootstrap cannot be applied for cluster robust inference when the sample sizes vary across clusters.} First and most importantly, existing simulation studies show that the cluster pairs bootstrap does not perform well with a small number of clusters. For example, cameron2008bootstrap use a simulation design based on bertrand2004much and explain that the poor performance of the cluster pairs bootstrap is due to the fact that the resampled values of the Gram matrix $X'X$ are nearly singular. Our approach of analytically inverting the Cram\'er-Edgeworth expansion avoids this problem by not resampling $X'X$ while achieving asymptotic refinement. Second, the cluster pairs bootstrap uses independence across $X_1,\ldots,X_G$, while our expansion does not require it. We allow $X_1,\ldots,X_G$ to be correlated, e.g., through adaptive randomization. Third, the standard proof for the cluster pairs bootstrap's asymptotic refinement liu1988bootstrap,hall2013bootstrap excludes discrete regressors. To apply the results from hall2013bootstrap to the regression framework, we need to impose the Cram\'er's condition on the regressors, but the Cram\'er's condition fails for discrete random variables bhattacharya2010normal.\footnote{We could not find a sufficient condition for the cluster pairs bootstrap's asymptotic refinement that allows for discrete regressors and non-identical distributions. A weaker version of the Cram\'er's condition has been proposed, e.g., in bai1991edgeworth, but the cluster pairs bootstrap's asymptotic refinement without the classical Cram\'er's condition is beyond the scope of this paper.} Last, our proposed critical value has a closed-form expression, and for this reason, its computation is much faster than resampling.

A number of other methods have been proposed for the construction of standard errors and confidence intervals. Our proposal is unique in the following sense: it demonstrates good finite-sample performance in the simulation designs based on bertrand2004much, while simultaneously achieving third-order asymptotic refinement. For example, cameron2008bootstrap propose the wild cluster bootstrap and demonstrate its good finite-sample properties in simulations. djogbenou2019asymptotic provide the asymptotic size control of the wild cluster bootstrap, and canay2021wild show the size control of the wild cluster bootstrap even with a fixed number of the clusters as long as there are a large number of observations per cluster. However, Theorem 5.2 of djogbenou2019asymptotic shows that the wild cluster bootstrap does not achieve asymptotic refinement when the score has non-zero skewness. Our simulation results in Section (ref) confirm this result of djogbenou2019asymptotic, and our proposed method demonstrates better size control than the wild cluster bootstrap in such cases.

This paper is also related to the use of the $t$-distribution for the critical value and an adjustment of degrees of freedom in the $t$-distribution mccaffrey2003bias,bester2011inference,imbens2016robust,young2016improved,hansen2021exact,hansen2022jackknife. These papers assume normal and homoscedastic errors. On the other hand, our approach does not rely on normal errors or a specific covariance structure for the error terms within a cluster.

The paper proceeds as follows. Section (ref) formally introduces the regression model with clustered errors and the new critical value. Section (ref) presents Monte Carlo simulations for the proposed critical value and existing ones. Section (ref) concludes. The appendix collects all the proofs and additional results.

Inverting a Cram\'er-Edgeworth Expansion with Clustered Errors

We have the dataset of $\{(Y_{ig},X_{ig}):i=1,\ldots,N_g,\ g=1,\ldots,G\}$ and consider the regression model $Y_{ig}=X_{ig}'\beta+u_{ig}$ with $E[u_{ig}\mid X_{1g},\ldots,X_{N_gg}]=0$ and $\dim(X_{ig})=k$. We assume the observations $\{(Y_{ig},X_{ig}):i=1,\ldots,N_g\}$ are independent across $g=1,\ldots,G$. To emphasize the fact that we use the independence across $g$, we use the matrix notation with $Y_g=(Y_{1g},\ldots,Y_{N_gg})'$ and $X_g=(X_{1g}',\ldots,X_{N_gg}')'$, and write the regression model succinctly as $$ Y_{g}=X_{g}\beta+u_{g}\mbox{ with }E[u_{g}\mid X_{g}]=0. $$ For a fixed vector $\lambda$ and a hypothesized value $c_0$ for $\lambda'\beta$, we consider the hypothesis testing problem of $$ H_0:\lambda'\beta=c_0\mbox{ vs }H_1: \lambda'\beta\ne c_0 $$ with the significance level $\alpha\in (0,1)$. The OLS estimator for $\beta$ is $$ \hat{\beta}=\left(\frac{1}{G}\sum_{g=1}^GX_g'X_g\right)^{-1}\left(\frac{1}{G}\sum_{g=1}^GX_g' Y_g\right). $$ Consider an asymptotic variance estimator for $\lambda'\hat{\beta}$ defined by $$ \hat\sigma^2=\frac{1}{G}\sum_{g=1}^G(\lambda'\Pi X_g'\hat{u}_g)^2\mbox{ with }\Pi=\left(\frac{1}{G}\sum_{g=1}^GX_{g}'X_{g}\right)^{-1}\mbox{ and }\hat{u}_g=Y_g-X_g\hat\beta. $$ Define the $t$-statistic by $$ t=\sqrt{G}\frac{\lambda'\hat{\beta}-c_0}{\hat\sigma}. $$ From now on, we estimate the null distribution for the above $t$-statistic $t$ and construct a critical value for it. We treat the covariates $\mathbf{X}=\{X_{g}\}^\infty_{g=1}$ as fixed, so we investigate $Pr(|t|\leq z\mid\mathbf{X}=\mathbf{x})$ for a given sequence of constants $\mathbf{x}=\{x_{g}\}^\infty_{g=1}$.

We consider the numerator and denominator of $$ t=\sqrt{G}\frac{(\lambda'\hat{\beta}-c_0)/\sigma}{\hat\sigma/\sigma} $$ under the null hypothesis $H_0$, where $\sigma^2$ is the asymptotic variance for $\lambda'\hat{\beta}$ defined by $$ \sigma^2=\frac{1}{G}\sum_{g=1}^G\sigma_g^2\mbox{ with }\sigma_{g}^2=E\left[(\lambda'\Pi X_g'u_g)^2\mid\mathbf{X}=\mathbf{x}\right]. $$ The numerator has the linear representation of $$ (\lambda'\hat{\beta}-c_0)/\sigma=\frac{1}{G}\sum_{g=1}^G\omega_{1g}\mbox{ with } \omega_{1g}=\sigma^{-1}\lambda'\Pi X_g'u_g. $$ The square of the denominator $\hat\sigma^2/\sigma^2$ has the following quadratic representation. The proof is given in Section (ref).

lem$$ \hat\sigma^2/\sigma^2=1-\left(\frac{1}{G}\sum_{g=1}^G\omega_{2g}\right)'\Gamma\left(\frac{1}{G}\sum_{g=1}^G\omega_{2g}\right)+\frac{1}{G}\sum_{g=1}^G\omega_{3g} $$ where $$ \Gamma=\left( \begin{array}{cc} -\frac{1}{G}\sum_{g=1}^GX_g'X_g\Pi'\lambda\lambda'\Pi X_g'X_g&I_k\\ I_k&O \end{array} \right) $$ and \begin{align*} \omega_{2g}&=\sigma^{-1}\left(\begin{array}{c}I_k\\ X_g'X_g\Pi'\lambda\lambda'\end{array}\right) \Pi X_g'u_g, \\ \omega_{3g}&=\sigma^{-2}((\lambda'\Pi X_g'u_g)^2-\sigma_{g}^2). \end{align*}

To approximate the null distribution for the $t$-statistic, we approximate the distribution of $\frac{1}{\sqrt{G}}\sum_{g=1}^G(\omega_{1g},\omega_{2g}',\omega_{3g})'$ up to $o(G^{-1})$. For this purpose, we use the following moments of $\omega_{1g}$ and $\omega_{2g}$:

align*[align* omitted — 406 chars of source]

We assume these moments are bounded.

assumption$\mu_{1,2}$, $\mu_{2,2}$, $\mu_{1,1,1}$, and $\mu_{1,1,1,1}$ are bounded uniformly in $G$.\footnote{All the results in this paper only require the assumptions hold for sufficiently large $G$. To simplify the exposition, we impose the assumptions for every value of $G$.}

The (second-order) Cram\'er-Edgeworth expansion for the $t$-statistic's null distribution is expressed under the following assumption, for which we provide a sufficient condition in Section (ref).

assumptionUnder $H_0$, there is a neighborhood $\mathcal{N}$ of $\Phi^{-1}(1-\alpha/2)$ such that $$ \sup_{z\in\mathcal{N}}|Pr(|t|\leq z\mid\mathbf{X}=\mathbf{x})-(2 \Phi(z) - 1 +2G^{-1}q_2(z)\phi(z))| = o(G^{-1}), $$ where $He_r$ are the $r$-th order Hermite polynomials and $$ q_2(z)=-\left(\frac{1}{2}(k_2+k_1^2)He_1(z)+\frac{1}{24}(k_4+4k_1k_3)He_3(z)+\frac{1}{72}k_3^2He_5(z)\right). $$ The parameters $k_1,\ldots,k_4$ are defined as \begin{align*} k_{1}&=\nu_1 \\ k_{2}&=\nu_2-\nu_1^2 \\ k_{3}&=\nu_3-3\nu_1 \\ k_{4}&=\nu_4-4\nu_1\nu_3-6\nu_2+12\nu_1^2 \end{align*} where \begin{align*} \nu_1&=-\frac{\mu_{1,1,1}}{2} \\ \nu_2&=2\mu_{1,1,1}^2 + (\mu_{2,2} +2\mu_{1,2}'\Gamma\mu_{1,2}) \\ \nu_3&=-\frac{7}{2}\mu_{1,1,1} \\ \nu_4&=-2\mu_{1,1,1,1} +28\mu_{1,1,1}^2 +6\mu_{2,2} +24\mu_{1,2}'\mu_{1,2}. \end{align*}

The function $q_2(z)$ in the (second-order) Cram\'er-Edgeworth expansion is unknown since we do not know the population objects of $\mu_{1,2}'\mu_{1,2},\mu_{2,2},\mu_{1,1,1},\mu_{1,1,1,1}$. We can estimate them using their sample analogs:

align*[align* omitted — 393 chars of source]

where

align*[align* omitted — 199 chars of source]

We can construct the estimator $\hat{q}_2(z)$ for $q_2(z)$ using these sample analogs. Our proposed critical value is $$ \hat{cv}=\Phi^{-1}(1-\alpha/2)- G^{-1}\hat q_2(\Phi^{-1}(1-\alpha/2)) $$ and the resulting confidence interval for $\lambda'\beta$ is $\lambda'\hat\beta\pm \hat{cv}\sqrt{\hat\sigma/G}$.

We impose the following conditions on moments.

assumption(i) $\sigma^2$ and the minimum eigenvalue of $\frac{1}{G}\sum_{g=1}^Gx_{g}'x_{g}$ are bounded away from zero uniformly in $G$. (ii) $\frac{1}{G}\sum_{g=1}^G\left\| x_g'x_g\right\|^4$ and $\frac{1}{G}\sum_{g=1}^GE[\left\| X_g'u_g\right\|^j\mid\mathbf{X}=\mathbf{x}]$ are bounded uniformly in $G$ for every positive integer $j$.

The second condition makes the characteristic function of $\frac{1}{G}\sum_{g=1}^G\omega_g$ infinitely differentiable and simplifies the proofs. We may weaken it to the bounded moment condition up to a certain order by truncating $\omega_g$ (cf., hall2013bootstrap; bhattacharya1978validity). The bounded higher-order moments are crucial for our inference because we estimate the population objects of $\mu_{1,2},\mu_{2,2},\mu_{1,1,1},\mu_{1,1,1,1}$. For example, $\mu_{1,1,1,1}$ is the average fourth moment of $\omega_{1g}$. In the proof, we use the fourth moment of its estimator and thus require bounded 16th moments.

In Theorem (ref) below, we show the size control for the proposed critical value. Note that even if Assumption (ref) fails, we can still achieve size control with an asymptotic approximation error of the standard rate $O(G^{-1})$, as long as asymptotic normality holds. This point resembles the efficiency gain of the feasible generalized least squares estimation comeron/miller:2026.

thm$Pr(|t|\leq \hat{cv}\mid\mathbf{X}=\mathbf{x})=1-\alpha+o(G^{-1})$ under $H_0$ and Assumption (ref)-(ref).

The proof is provided in Section (ref).

remOur proposed confidence interval does not rely on a particular choice of the asymptotic variance estimator in the following sense. Suppose we use another asymptotic variance estimator for $\lambda'\hat{\beta}$ denoted by $\tilde\sigma^2$. Then the $t$-statistic becomes $$ \sqrt{G}\frac{\lambda'\hat{\beta}-c_0}{\tilde\sigma}=\frac{\hat\sigma}{\tilde\sigma}t. $$ We can accordingly change the critical value as $(\hat\sigma/\tilde\sigma)\hat{cv}$, but the resulting confidence interval is still $\lambda'\hat\beta\pm \hat{cv}\sqrt{\hat\sigma/G}$.

remIn this paper, we treated $\mathbf{X}$ as fixed, so we do not apply the Cram\'er-Edgeworth expansion to $\frac{1}{G}\sum_{g=1}^GX_{g}'X_{g}$. It is different from djogbenou2019asymptotic, who consider the regressor $\mathbf{X}$ as random when they derive the Cram\'er-Edgeworth expansion for the $t$-statistic $t$.\footnote{They do not use the Cram\'er-Edgeworth expansion to construct an asymptotically refined inference, but instead use it to show the wild bootstrap cannot achieve the asymptotic refinement of $o(G^{-1})$. As a result, this remark is not relevant for their analysis on the wild bootstrap because their conclusion relies on the fact that the wild bootstrap cannot replicate the first four cumulants of the Cram\'er-Edgeworth expansion.} First, our expansion does not require the independence of $X_{g}$ across $g$. As a consequence, we allow $X_1,\ldots,X_G$ to be generated from an adaptive randomization. Second, the assumption in djogbenou2019asymptotic for the Cram\'er-Edgeworth expansion excludes empirically relevant cases such as binary regressors. With the random regressor, they apply the Cram\'er-Edgeworth expansion to the term $\frac{1}{G}\sum_{g=1}^GX_{g}'X_{g}$ and assume the Cram\'er's condition on $X_{g}'X_{g}$. This condition fails if $X_g$ includes discrete variables bhattacharya2010normal.\footnote{Consider the case where $X_{ig}$ includes a binary variable $D_g$. The matrix $X_g'X_g$ includes $N_gD_g^2$, which takes the two values of $N_g$ and $0$. By choosing the vector $t$ such that $t'\mathrm{vech}(X_{g}'X_{g})=\|t\|N_gD_g^2$, we have $$ E[\exp(it'\mathrm{vech}(X_{g}'X_{g}))] = Pr(D_g=1)\exp(i\|t\|N_g)+(1-Pr(D_g=1)) = 1\mbox{ if $\|t\|N_g$ is a multiple of }2\pi. $$ It implies that this Cram\'er's condition fails.} In this paper, we do not require the Cram\'er's condition on $X_{g}'X_{g}$. Third, the resulting Cram\'er-Edgeworth expansion has fewer unknown parameters to estimate. Namely, we do not need to estimate the moments of $X_{g}'X_{g}$.

remThe above Cram\'er-Edgeworth expansion for the $t$-statistic in Assumption (ref) is different from the one for the (standardized) sample mean $\sigma^{-1}\frac{1}{\sqrt{G}}\sum_{g=1}^G\lambda'\Pi X_g'u_g$ since we estimate the average variance $\sigma^2$. One of the qualitative differences is that we can estimate the above $q_2(z)$ using the sample analogs for the coefficients, while we cannot estimate the Cram\'er-Edgeworth expansion for $\sigma^{-1}\frac{1}{\sqrt{G}}\sum_{g=1}^G\lambda'\Pi X_g'u_g$ without additional assumptions. In this sense, even if we know the value of $\sigma^2$, we need to use the estimated value of $\sigma^2$ for the inference strategy proposed in this paper.\footnote{The Cram\'er-Edgeworth expansion for $\sigma^{-1}\frac{1}{\sqrt{G}}\sum_{g=1}^G\lambda'\Pi X_g'u_g$ involves the average squared variance $$ \theta = \frac{1}{G}\sum_{g=1}^G(\sigma_{g}^2)^2 = \frac{1}{G}\sum_{g=1}^G \left(\lambda'\Pi X_g'E[u_gu_g'\mid\mathbf{X}=\mathbf{x}]X_g\Pi\lambda \right)^2 $$ and the error variance matrix $E[u_gu_g'\mid\mathbf{X}=\mathbf{x}]$ is not specified in this paper. To illustrate how the average squared variance appears in the Cram\'er-Edgeworth expansion, we assume $\sigma=1$. In order to approximate the distribution of $\frac{1}{\sqrt{G}}\sum_{g=1}^G\lambda'\Pi X_g'u_g$ up to $o(G^{-1})$, we characterize $\log(E[\exp(iz\lambda'\Pi X_g'u_g/\sqrt{G})])$ up to $o(G^{-2})$. We have the Taylor expansion of $$ E[\exp(iz\lambda'\Pi X_g'u_g/\sqrt{G})\mid\mathbf{X}=\mathbf{x}] = 1 +\sigma_{g}^2\frac{(iz/\sqrt{G})^2}{2!} +\cdots. $$ Since $\log(1+z)=z-\frac{1}{2}z^2+O(\|z\|^3)$ as $z\rightarrow 0$, we have $$ \log(E[\exp(iz\lambda'\Pi X_g'u_g/\sqrt{G})\mid\mathbf{X}=\mathbf{x}]) = \sigma_{g}^2\frac{(iz/\sqrt{G})^2}{2!} -\frac{1}{2}\left(\sigma_{g}^2\frac{(iz/\sqrt{G})^2}{2!}\right)^2+\cdots. $$ By the independence across $g$, the log characteristic function of $\frac{1}{\sqrt{G}}\sum_{g=1}^G\lambda'\Pi X_g'u_g$ is the sum of the log characteristic function of $\lambda'\Pi X_g'u_g/\sqrt{G}$ and therefore it has the average squared variance $\frac{1}{G}\sum_{g=1}^G(\sigma_{g}^2)^2$. In the above equation, the squared variance $\sigma_{g}^2$ appears because the logarithm function is nonlinear and we use the quadratic approximation to achieve the remainder term of $o(G^{-1})$.}

Cram\'er-Edgeworth Expansion in Assumption (ref)

In this subsection, we provide a sufficient condition for the Cram\'er-Edgeworth expansion in Assumption (ref).

assumptionThere are random vectors $\eta_1,\ldots,\eta_G$, and matrices $\mathbf{M}_1$ and $\mathbf{M}_2$ such that $\omega_{g}=\mathbf{M}_1\eta_g$ and $\eta_{g}=\mathbf{M}_2\omega_g$ and that the minimum eigenvalue of $\mathbf{V}_G=Var(\frac{1}{\sqrt{G}}\sum_{g=1}^G\eta_g\mid\mathbf{X}=\mathbf{x})$ is bounded away from zero uniformly in $G$.

This assumption removes the redundant or duplicated elements from $\omega_{g}$. This removal is necessary to normalize the random variable $\frac{1}{\sqrt{G}}\sum_{g=1}^G\eta_g$ by using the matrix square root of the variance matrix $\mathbf{V}_G$.

assumptionThere are positive numbers $C,R,b$ with $b<2$ such that $\frac{1}{G}\sum_{g=1}^G\left|E\left[\exp\left(it'\eta_g\right)\mid\mathbf{X}=\mathbf{x}\right]\right|\leq 1-C\|t\|^{-b}$ for every $t$ with $\|t\|>R$ and for sufficiently large $G$.

The assumption is the mean weak Cram\'er's condition proposed in angst2017weak. We provide a sufficient condition for Assumption (ref).

thmSuppose $u_g$ has an absolutely continuous distribution given $\mathbf{X}=\mathbf{x}$ (with respect to the $N_g$-dimensional Lebesgue measure), and there are positive numbers $R,c_1,c_2$ such that $|\mathbf{B}(t)|\geq c_2G$ for every $t$ with $\|t\|=1$ and sufficiently large $G$, where $\mathbf{B}(t)=\{g: \pi^2/(16R^2)\leq Var(t'\eta_g\mid\mathbf{X}=\mathbf{x})\leq c_1\mathrm{med}(f_{t'\eta_g\mid\mathbf{X}=\mathbf{x}}(t'\eta_g)\mid\mathbf{X}=\mathbf{x})^{-2}\}$.\footnote{$f_{t'\eta_g\mid\mathbf{X}=\mathbf{x}}$ is the probability density function of $t'\eta_g$ given $\mathbf{X}=\mathbf{x}$ for every $g$ with $Var(t'\eta_g\mid\mathbf{X}=\mathbf{x})>0$. For every $t$, the distribution of $t'\eta_g$ given $\mathbf{X}=\mathbf{x}$ is either a point mass of zero or absolutely continuous (with respect to the $1$-dimensional Lebesgue measure). Therefore, if $Var(t'\eta_g\mid\mathbf{X}=\mathbf{x})>0$ then the distribution of $t'\eta_g$ given $\mathbf{X}=\mathbf{x}$ has a probability density function.} Then Assumption (ref) holds.

The proof is given in Section (ref).

The above three assumptions constitute a sufficient condition for Assumption (ref).

thmAssumptions (ref)-(ref) imply Assumption (ref).

The proof is given in Section (ref).

Monte Carlo simulations

In this section, we investigate the finite-sample performance of the critical value proposed in Section (ref) using simulated data. We compare our method (denoted by “Analytical” in the figures) with a few existing methods, such as (i) the $t_{G-1}$ critical value (“Student”), (ii) the restricted wild cluster bootstrap with Rademacher weights by cameron2008bootstrap (“CWB”), and (iii) the pairs percentile-$t$ cluster bootstrap (“Pairs”). We use 10,000 simulations and 1000 draws for bootstrap procedures. We consider two designs that are challenging for existing methods but can be accommodated well using our analytic corrections. The first design features binary regressors, which is challenging for the pairs cluster bootstrap, while the second design features skewed errors, which is challenging for the cluster wild bootstrap.

figure[figure omitted — 316 chars of source]

figure[figure omitted — 319 chars of source]

Design with Binary Regressors

In this section, we follow the simulation design from bertrand2004much and cameron2008bootstrap. It uses a state-year panel of excess earnings from 1979 to 1999 based on the Current Population Survey.\footnote{We use the data from the replication package of cameron2015practitioner: \url{https://cameron.econ.ucdavis.edu/research/papers.html} In particular, this simulation exercise uses the variable lnwage from CPS_panel.dta from 1979 to 1999.} For each simulation draw, we randomly select $G$ out of 50 states with replacement. We randomly select the policy change time uniformly from $\{1984,\ldots,1993\}$ and assume that half of the $G$ states experience the policy change after the selected time period. We construct the policy dummy variable accordingly. By definition, this policy dummy variable has a zero coefficient in the population. We regress the excess earnings on the policy dummy, the year dummies, and the state dummies, and conduct the significance test for the coefficient of the policy dummy variable.

In the data generating process of this section, the skewness of the score is close to zero (with $\hat\mu_{1,1,1}=0.02$ for $G=10^4$). Although the Cram\'er condition does not hold for this design because all the variables are discrete bhattacharya2010normal, we can consider whether a critical value accounts for the skewness and kurtosis of the $t$-statistic, which are the key components of the second-order Cram\'er-Edgeworth expansion. Our proposed method matches these moments by estimating them explicitly. At the same time, the wild cluster bootstrap approximates the skewness and kurtosis well because it uses zero skewness and estimates the kurtosis consistently djogbenou2019asymptotic.

Figure (ref) shows the rejection probabilities for different methods. As documented in cameron2008bootstrap, the pairs cluster bootstrap under-rejects for small values of $G$ (e.g., $G=6, 10$), while the inference based on the $t_{G-1}$ critical value over-rejects. All methods control the size approximately when $G$ is sufficiently large (e.g., $G=50$). Our proposed method exhibits comparable performance to the wild cluster bootstrap, even for a small value of $G=10$.

Design with a Skewed Error Distribution

To compare the methods when the error has large skewness, we consider the case with $N_g=1$, $X_{ig}=1$, and $Y_{ig}$ follows the exponential distribution with unit mean. This design is used in Section 3 of hall1983inverting for one-sided tests. Since the error has a skewness of $2$, Theorem 5.2 of djogbenou2019asymptotic implies that the wild cluster bootstrap does not have asymptotic refinement in this case.

Figure (ref) shows the rejection probabilities for different methods with the skewed error distribution. Again, the rejection probabilities of all the methods approach the prespecified significance level ($5\%$) as $G$ increases, which confirms their asymptotic validity. However, the finite-sample performance differs. The rejection probabilities of our proposed method approach $5\%$ faster than the wild cluster bootstrap. This is consistent with the theoretical fact that the wild cluster bootstrap does not have asymptotic refinement in this data distribution. In contrast, our analytical approach and the pairs cluster bootstrap both achieve the asymptotic refinement. It explains why these methods are similar to each other and much closer to nominal size than the other two methods in Figure (ref).

Conclusion

In this paper, we propose an inference method for linear regression with clustered errors, and it achieves third-order asymptotic refinements. Unlike the cluster pairs bootstrap, it does not resample the Gram matrix of $\frac{1}{G}\sum_{g=1}^GX_g'X_g$, thus avoiding the small-sample issues of the cluster pairs bootstrap cameron2008bootstrap. Our simulation results show favorable finite-sample performance of the proposed method. Notably, it works comparably to the wild bootstrap in the simulation design based on bertrand2004much and for some designs (with skewed distributions) it has better size control than the wild bootstrap.