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.
111,792 characters · 14 sections · 64 citation commands
Fast and Reliable Jackknife and Bootstrap Methods for Cluster-Robust Inference
\onehalfspacing
In applications of linear regression models to many fields of economics and other disciplines, it is common to divide the sample into disjoint clusters and employ a cluster-robust variance matrix estimator (or CRVE) for inference. These estimators are based on the assumption that the disturbances of the regression model are uncorrelated across clusters, but they allow for arbitrary patterns of dependence and heteroskedasticity within each cluster. The literature on cluster-robust inference has grown rapidly in recent years. CM_2015 is a classic survey article. \citet*{CGH_2018} surveys a broader class of methods for dependent data. \citet*{MNW-guide} provides a guide that explores the implications of key theoretical results for empirical practice, with an emphasis on bootstrap methods.
There are several CRVEs for ordinary least squares (OLS) estimates of linear regression models; see (ref). However, because the one known as CV$_{\kern -0.08333em1}$ is easy to compute and is the default in Stata, almost all empirical work to date has used it. Cluster-robust tests and confidence intervals based on CV$_{\kern -0.08333em1}$ may or may not yield reliable inferences. Whether they do so depends primarily on the number of clusters $G$ and how homogeneous these are. When all clusters are roughly equal in size and approximately balanced, asymptotic inference based on CV$_{\kern -0.08333em1}$ seems to be fairly reliable whenever $G$ is at least moderately large (say 50 or more). However, even when $G$ is very large, cluster-robust $t$-tests and Wald tests are at risk of severe over-rejection, and cluster-robust confidence intervals are at risk of severe under-coverage in at least two situations. The first is when one or a few clusters are much larger than the rest, and the second is when the only “treated” observations belong to just a few clusters; \citet*{DMN_2019} discusses the first case, and \citet*{MW-JAE,MW-EJ} discuss the second.
Alternatives to CV$_{\kern -0.08333em1}$ have been known since BM_2002. The first contribution of the present paper, which is discussed in (ref), is to provide a fast method for computing jackknife\kern 0.04167em-based CRVEs, of which the simplest is generally known as CV$_{\kern -0.08333em3}$. By explicitly using the cluster jackknife for computation in an efficient way, our method makes it feasible to employ CV$_{\kern -0.08333em3}$ for inference even in very large samples with a large number of clusters.
Because CV$_{\kern -0.08333em3}$ standard errors used to be hard to compute, there has been very little work comparing the finite\kern 0.04167em-sample performance of $t$-tests based on CV$_{\kern -0.08333em3}$ with those of similar procedures based on CV$_{\kern -0.08333em1}$; a partial exception is NW_2022. The second contribution of this paper is to compare the finite\kern 0.04167em-sample properties of these tests, and also ones based on CV$_{\kern -0.08333em2}$, by simulation; see (ref). In concurrent work that cites our simulations, Hansen-jack provides important theoretical results which suggest that asymptotic inference based on CV$_{\kern -0.08333em3}$ is generally more reliable, and more conservative, than asymptotic inference based on CV$_{\kern -0.08333em1}$.
Existing bootstrap methods for cluster-robust inference are all based on CV$_{\kern -0.08333em1}$. The best known of these (and until now the best performing one) seems to be the wild cluster restricted (or WCR) bootstrap proposed in \citet*{CGM_2008}. There is also a closely related procedure called the wild cluster unrestricted (or WCU) bootstrap, which generally does not work quite as well. The asymptotic validity of these procedures is proved in DMN_2019, which also analyzes their higher-order asymptotic properties. Until a few years ago, the WCR and WCU bootstraps were computationally expensive for large samples, but that is no longer the case. \citet*{RMNW} describes a remarkably efficient implementation in the Stata package boottest, and JGM-fast discusses other methods for fast computation. The boottest routines are now available as a Julia package which can be also be called from R, Python, and \texttt{Stata}. The package \texttt{fwildclusterboot} implements the \texttt{boottest} method natively in \texttt{R} Fischer_2022.
The third contribution of this paper is to propose several new variants of the wild cluster bootstrap. One modification simply replaces CV$_{\kern -0.08333em1}$ by CV$_{\kern -0.08333em3}$. The other, which requires some new results, involves modifying the bootstrap data-generating process, or DGP. Modern treatments of the wild cluster bootstrap, such as MNW-guide, express the bootstrap DGP as a function of the empirical scores. We show how to make the bootstrap DGP more closely resemble the (unknown) true DGP by transforming the residuals before forming the scores. The transformation we propose is based on the jackknife. Accordingly, it does not actually require any calculations that explicitly involve residuals. This makes it very fast when the number of clusters is small relative to the sample size, even when the latter is extremely large.
The next section establishes notation and briefly reviews the literature on asymptotic cluster-robust inference for the linear regression model. (ref) then discusses a computational method for CV$_{\kern -0.08333em3}$ which is conceptually simple and extremely fast in many cases, as we demonstrate in (ref). Next, (ref) discusses several ways of modifying the wild cluster bootstrap. Simulation results in (ref) suggest that our new versions of the WCR and WCU bootstraps perform better, sometimes very much better, than the original ones. This is particularly true when cluster sizes vary greatly. One modified version of the WCR bootstrap that uses transformed scores seems to work especially well in most settings. (ref) presents three empirical examples in which our methods are likely to be more reliable than existing ones. (ref) concludes with a brief discussion of the methods that we recommend in practice.
Consider the linear regression model $y_i = {\bm{x}}_i^\top\!{\bm\beta} + u_i$. If we divide the data into $G$ disjoint clusters, where the allocation of observations to clusters is assumed to be known, this can be written as
The \th{g} cluster has $N_g$ observations, and the total sample size is $N = \sum_{g=1}^G N_g$. In \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, ${\bm{X}}_g$ is an $N_g\times k$ matrix of regressors, ${\bm\beta}$ is a $k$-vector of coefficients, ${\bm{y}}_g$ is an $N_g$-vector of observations on the regressand, and ${\bm{u}}_g$ is an $N_g$-vector of disturbances (or error terms). Stacking the ${\bm{y}}_g$ yields the $N$-vector ${\bm{y}}$, stacking the ${\bm{X}}_g$ yields the $N\times k$ matrix ${\bm{X}}$\kern -0.08333em, and stacking the ${\bm{u}}_g$ yields the $N$-vector ${\bm{u}}$, so that \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} can be rewritten as ${\bm{y}} = {\bm{X}}\kern -0.08333em{\bm\beta} + {\bm{u}}$.
The OLS estimator of ${\bm\beta}$ is
where the second equality depends on the assumption that the data are actually generated by \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}} with true value ${\bm\beta}_0$. Thus, if ${\bm{s}}_g = {\bm{X}}_g^\top {\bm{u}}_g$ is the score vector for the \th{g} cluster,
Obtaining valid inferences evidently requires assumptions about the score vectors. For a correctly specified model, ${\rm E}({\bm{s}}_g)={\bm{0}}$ for all $g$. We further assume that
where ${\bm{\Sigma}}_g$ is the symmetric, positive semidefinite variance matrix of the scores for the \th{g} cluster. The second assumption in \hyperref[{eq:Sigma_g}]{\tagform@{\ref*{eq:Sigma_g}}} is crucial. It states that the scores for every cluster are uncorrelated with the scores for every other cluster.
From the rightmost expression in \hyperref[{eq:betahat}]{\tagform@{\ref*{eq:betahat}}}, we see that the distribution of $\hat{\bm\beta}$ depends on the disturbance subvectors ${\bm{u}}_g$ only through the distribution of the score vectors ${\bm{s}}_g$. It follows immediately that an estimator of $\operatorname{Var}(\hat{\bm\beta})$ should be based on the usual sandwich formula,
Every CRVE replaces the ${\bm{\Sigma}}_g$ in \hyperref[{eq:trueV}]{\tagform@{\ref*{eq:trueV}}} by functions of the ${\bm{X}}_g$ and the residual subvectors $\hat{\bm{u}}_g$. There is more than one way to do this. Since ${\bm{\Sigma}}_g$ is the expectation of ${\bm{s}}_g{\bm{s}}_g^\top$, the simplest approach is just to replace it by $\hat{\bm{s}}_g\hat{\bm{s}}_g^\top$, where $\hat{\bm{s}}_g = {\bm{X}}_g^\top \hat{\bm{u}}_g$ is the empirical score vector for the \th{g} cluster. If in addition we multiply by a correction for degrees of freedom, we obtain
This is by far the most widely-used CRVE in practice, and it is the default in Stata. The leading scalar is chosen so that, when $G=N$\kern -0.08333em, $\hat{\bm{V}}_1(\hat{\bm\beta})$ reduces to the familiar HC$_1$ estimator MW_1985 that is robust only to heteroskedasticity of unknown form.
Inference typically relies on cluster-robust $t$-statistics and Wald statistics based on \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}. If $\beta_j$ denotes the \th{j} element of ${\bm\beta}$, $\hat\beta_j$ the OLS estimate, and $\beta_{0j}$ its value under the null hypothesis, then the appropriate $t$-statistic is
where $\textrm{se}_1(\hat\beta_j)$ is the square root of the \th{j} diagonal element of $\hat{\bm{V}}_1(\hat{\bm\beta})$. Under extremely strong assumptions, \citet*{BCH_2011} shows that $t_j$ asymptotically follows the $t(G-1)$ distribution. Conventional “asymptotic” inference is based on this distribution.
We should expect inferences based on CV$_{\kern -0.08333em1}$ to be reliable if the sum of the ${\bm{s}}_g$, suitably normalized, is well approximated by a multivariate normal distribution with mean zero, and if the ${\bm{s}}_g$ are well approximated by the $\hat{\bm{s}}_g$. But asymptotic inference can be misleading when either or both of these approximations is poor; see DMN_2019 and MNW-guide. Whether or not the first approximation is a good one depends on the model and the data, and there is not much the investigator can do about it. But the second approximation can, in principle, be improved by using modified empirical score vectors instead of the $\hat{\bm{s}}_g$.
Two CRVEs based on this idea, usually known as CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$, were proposed (under different names) in BM_2002. These are the cluster analogs of the hetero\-skedasticity-consistent variance matrix estimators HC$_2$ and HC$_3$ proposed in MW_1985. All of these estimators are designed to compensate, in different ways, for the shrinkage and intra-cluster correlation of the residuals induced by least squares.
The CV$_{\kern -0.08333em2}$ variance matrix is
where the modified score vectors $\grave{\bm{s}}_g$ are defined as
Here ${\bm{M}}_{gg} = {\bf I}_{N_g} - {\bm{X}}_g({\bm{X}}^\top{\bm{X}})^{-1}\kern -0.08333em{\bm{X}}_g^\top$ is the \th{g} diagonal block of the projection matrix ${\bm{M}}_{\bm{X}}$, which satisfies $\hat{\bm{u}} = {\bm{M}}_{\bm{X}}{\bm{u}}$, and ${\bm{M}}_{gg}^{-1/2}$ is the symmetric square root of its inverse. The CV$_{\kern -0.08333em2}$ estimator has been recommended in Imbens_2016 and PT_2018. Following BM_2002, these papers provide methods for computing critical values based on $t$ and $F$ distributions with computed degrees of freedom.
The CV$_{\kern -0.08333em3}$ variance matrix is very similar to CV$_{\kern -0.08333em2}$, but, as we explain in (ref), it is based on the jackknife. The usual definition is
where now the modified score vectors $\acute{\bm{s}}_g$ are defined as
The rescaling factor $(G-1)/G$ in \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} is the analog of the factor $(N-1)/N$ that occurs in jackknife variance matrix estimators at the individual level. This factor implicitly assumes that all clusters are the same size and perfectly balanced, with disturbances that are independent and homoskedastic; an alternative is proposed in NW_2022.
Computing either CV$_{\kern -0.08333em2}$ or CV$_{\kern -0.08333em3}$ using \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}} or \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} is extremely expensive, or even computationally infeasible, when any of the $N_g$ are large. The problem is that, before computing \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}, we apparently need to rescale the residual vector $\hat{\bm{u}}_g$ for each cluster. This involves storing and inverting the $N_g\times N_g$ matrix ${\bm{M}}_{gg}$. Before computing \hyperref[{eq:bis2}]{\tagform@{\ref*{eq:bis2}}}, we also need to compute the symmetric square roots of the ${\bm{M}}_{gg}$, and this requires calculating their eigenvalues and eigenvectors. Of course, when all clusters are very small, this is not difficult. When $G=N$\kern -0.08333em, CV$_{\kern -0.08333em2}$ reduces to HC$_2$, and CV$_{\kern -0.08333em3}$ reduces to HC$_3$, both of which can be computed very quickly.
NAAMW_2020 has recently proposed a method that is much faster for large clusters. Versions of this method apply to both CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$. Instead of rescaling the residual vectors, it calculates the score vectors $\grave{\bm{s}}_g$ or $\acute{\bm{s}}_g$ directly using equations that do not involve any $N_g\times N_g$ matrices. A revised version of this method, which appears to be new, works as follows. First, form the $k\times k$ matrices
Then, for \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}}, calculate the rescaled score vectors
and, for \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}}, calculate the rescaled score vectors
These rescaled score vectors are used in \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}} and \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} as before. Unless all the clusters are very small, computing CV$_{\kern -0.08333em2}$ and CV$_{\kern -0.08333em3}$ using \hyperref[{eq:bis2f}]{\tagform@{\ref*{eq:bis2f}}} and \hyperref[{eq:bis3f}]{\tagform@{\ref*{eq:bis3f}}} is much faster than computing them using \hyperref[{eq:bis2}]{\tagform@{\ref*{eq:bis2}}} and \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}. In the case of CV$_{\kern -0.08333em3}$, an even faster and more intuitive method is available. This jackknife\kern 0.04167em-based method, which we discuss in the next section, can be extremely fast when $N$ is large and $G$ is much smaller than $N$\kern -0.08333em, so that at least some clusters are large; see (ref).
The jackknife is a simple method for reducing bias and estimating standard errors by omitting observations sequentially. Tukey_1958 suggested using the jackknife to estimate standard errors, and Miller_1974 is a classic reference. The key idea of the cluster jackknife is to compute $G$ sets of parameter estimates, each of which omits one cluster at a time. In this section, we use it to compute two closely related CRVEs.
The OLS estimates of ${\bm\beta}$ when each cluster is omitted in turn are
It is easy to obtain the $\hat{\bm\beta}^{(g)}$ in a computationally efficient manner. We start by calculating the cluster-level matrices and vectors
Unless $G$ is very large, this involves very little cost beyond that of computing $\hat{\bm\beta}$, because we can use the quantities in \hyperref[{eq:subthings}]{\tagform@{\ref*{eq:subthings}}} to construct ${\bm{X}}^\top\!{\bm{X}}$ and ${\bm{X}}^\top{\bm{y}}$ and then use \hyperref[{eq:OLSbeta}]{\tagform@{\ref*{eq:OLSbeta}}} to obtain $\hat{\bm\beta}$. For typical values of $k$, it should then be reasonably inexpensive to calculate $\hat{\bm\beta}^{(g)}$ for every cluster using \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}. The main cost, beyond that of computing $\hat{\bm\beta}$, is that we need to calculate the inverse of a $k\times k$ matrix for each of the $\hat{\bm\beta}^{(g)}$\kern -0.08333em.
The cluster jackknife estimator of $\operatorname{Var}(\hat{\bm\beta})$ is the cluster analog of the usual jackknife variance matrix estimator given in Efron_81, among others. It is defined as
where $\bar{\bm\beta} = G^{-1}\sum_{g=1}^G \hat{\bm\beta}^{(g)}$ is the sample average of the $\hat{\bm\beta}^{(g)}$. Notice that \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} calculates the variance matrix around $\bar{\bm\beta}$. Centering around $\bar{\bm\beta}$ is common in jackknife variance estimation, but it is also common to center around $\hat{\bm\beta}$, as in BM_2002.
There is a very close relationship between $\hat{\bm{V}}_{3{\rm J}} (\hat{\bm\beta})$ and $\hat{\bm{V}}_3(\hat{\bm\beta})$. In fact,
which is just \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} with $\bar{\bm\beta}$ replaced by $\hat{\bm\beta}$. This follows from \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} and \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}} because
Note that the summation in \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}} is unchanged if $\hat{\bm\beta}^{(g)} - \hat{\bm\beta}$ is replaced by $\hat{\bm\beta} - \hat{\bm\beta}^{(g)}$.
Although the second equality in \hyperref[{eq:keyres}]{\tagform@{\ref*{eq:keyres}}} is not new, it will turn out to be very useful in (ref), and so we now prove it. The middle expression in \hyperref[{eq:keyres}]{\tagform@{\ref*{eq:keyres}}} can be written as
Using the Woodbury matrix identity,
$\hat{\bm\beta}^{(g)}$ can be written as the sum of four terms, the first of which is just $\hat{\bm\beta}$. Thus the right-hand side of \hyperref[{eq:keyres}]{\tagform@{\ref*{eq:keyres}}} can be written as
The last term in \hyperref[{eq:int2}]{\tagform@{\ref*{eq:int2}}} is identical to the last term in \hyperref[{eq:int1}]{\tagform@{\ref*{eq:int1}}}. The first two terms in \hyperref[{eq:int2}]{\tagform@{\ref*{eq:int2}}} can be rewritten as
where ${\bm{P}}_{gg} = {\bm{X}}_g ({\bm{X}}^\top{\bm{X}})^{-1}{\bm{X}}_g^\top$ is the \th{g} diagonal block of the matrix ${\bm{P}}_{\bm{X}} = {\bf I} - {\bm{M}}_{\bm{X}}$, so that ${\bm{P}}_{gg} = {\bf I}-{\bm{M}}_{gg}$. Inserting this straightforwardly yields the result that
The right-hand side of \hyperref[{eq:lastline}]{\tagform@{\ref*{eq:lastline}}} is the first term in \hyperref[{eq:int1}]{\tagform@{\ref*{eq:int1}}}, which proves the second equality in \hyperref[{eq:keyres}]{\tagform@{\ref*{eq:keyres}}}. When $N_g=1$ for all $g$, $\hat{\bm{V}}_{3{\rm J}}(\hat{\bm\beta})$ is numerically equal to the original HC$_3$ estimator proposed in MW_1985. The modern version of HC$_3$, which uses $\hat{\bm\beta}$ instead of $\bar{\bm\beta}$ and omits the factor of $N/(N-1)$, is due to DM_1993.
Both cluster jackknife estimators may be used to compute cluster-robust $t$-statistics. Since there are $G$ terms in the summation, it is natural to compare these with quantiles of the $t(G-1)$ distribution, as usual. These procedures should almost always be more conservative than $t$-tests based on CV$_{\kern -0.08333em1}$ Hansen-jack. We expect CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm{J}}}$ to be very similar in most cases. This issue will be investigated in (ref), where we conclude that it is reasonable to focus on CV$_{\kern -0.08333em3}$.
Both CV$_3$ and CV$_{\kern -0.08333em3{\rm J}}$ have been available in Stata for some years by using the options “vce(jackknife,mse)” and “vce(jackknife)”, respectively. However, the implementations discussed here are much more efficient when $G$ is not very small. They are available in Stata and R packages, both named summclust; see MNW_summclust and Fischer_summclust. Both packages also calculate a number of summary statistics that may be used to assess the reliability of cluster-robust inference as described in MNW-influence.
The matrix ${\bm{X}}^\top\!{\bm{X}} - {\bm{X}}_g^\top\!{\bm{X}}_g$ in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} can be singular for one or more values of $g$, so that at least some elements of $\hat{\bm\beta}^{(g)}$ cannot be identified. This can happen in otherwise well-specified models when there are cluster-level fixed effects. In that case, the solution is simply to partial them out before running the regression. In other cases where a singularity occurs, there are two possible courses of action. The first is to modify \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} and \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}} so that the summation is taken only over values of $g$ for which $\hat{\bm\beta}^{(g)}$ can be estimated, and $G$ is replaced by the number of clusters for which that is the case (this is the approach followed in the native Stata implementations; see also (ref)). When there are only a few problematic clusters, this approach may be attractive. But since $\hat{\bm\beta}$ and $\bar{\bm\beta}$ would then be based on different samples, it seems likely that CV$_{\kern -0.08333em3{\rm J}}$ and CV$_{\kern -0.08333em3}$ may differ more than they would usually do, which suggests that it may be safer to use the former.
The second course of action is to replace the inverse in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} by a generalized inverse. In practice, this means replacing coefficients that cannot be identified by zeros. When the elements of $\hat{\bm\beta}^{(g)}$ that are of primary interest can always be identified, this approach may be attractive, especially when there are many problematic clusters, as in the example of (ref).
The CV$_{\kern -0.08333em3}$ estimator can be challenging to compute. Following BM_2002, it is natural to employ what we call the “residual method” based on \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} and \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}. To compute the modified score vector $\acute{\bm{s}}_g$ for the \th{g} cluster, this method uses the $N_g$-vector of residuals $\hat{\bm{u}}_g$ and the $N_g\times N_g$ matrix ${\bm{M}}_{gg}^{-1}$. Unless every $N_g$ is small, storing and inverting the ${\bm{M}}_{gg}$ matrices is computationally expensive. Indeed, for even moderately large values of the $N_g$, this can be effectively impossible.
A much faster method, recently proposed in NAAMW_2020 and revised modestly in (ref), uses \hyperref[{eq:bis3f}]{\tagform@{\ref*{eq:bis3f}}} to obtain the modified score vectors $\acute{\bm{s}}_g$. Since it operates directly on the score vectors $\hat{\bm{s}}_g$, we call it the “score method.” An even faster approach, discussed in (ref), computes the $\hat{\bm\beta}^{(g)}$ using \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} and then calculates their variance matrix as \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}}. For obvious reasons, we refer to this as the “jackknife method.”
To compare timings for the residual, score, and jackknife methods, we generate two datasets with $N = 1,048,576 = 2^{20}$ observations. In one case, there are 20 regressors, and in the other case there are 40. The observations are divided into $G$ equal-sized clusters, where $G$ varies from 16 to $512K$ and $K$ denotes $1024=2^{10}$\kern -0.08333em. Thus the cluster size $M=N/G$ varies from 2 to $64K$.
(ref) shows the time in seconds, on a log$_2$ scale, for each of the three methods and the two datasets as a function of cluster sizes $M=N/G$, which vary from 2 to $64K$\kern -0.08333em. These times include the time required to compute the OLS estimates. For both the jackknife and score methods, there is considerable overlap between the computations needed for the OLS estimates and the ones needed for CV$_{\kern -0.08333em3}$. Thus, for large clusters, the cost of computing the OLS estimates and CV$_{\kern -0.08333em3}$ together using one or both of these methods was sometimes less than the cost of computing the OLS estimates alone. This is probably because of cache congestion, which seems to be alleviated by forming ${\bm{X}}^\top{\bm{X}}$ on a cluster-by-cluster basis. For large clusters, the speed of all methods could almost certainly be increased by using a fast BLAS implementation. However, in the interest of programming ease, we have not done this. The jackknife and score methods are already very fast.
In (ref), the residual method works well for very small values of $M$\kern -0.08333em. It is always the fastest method for $M\le4$. We did not perform any timings for $M=1$, where CV$_{\kern -0.08333em3}$ reduces to HC$_3$, because we would have needed a different program that eliminated the loops within each cluster to obtain optimal results. But the residual method is certainly the fastest one for this case. However, its cost rises very rapidly as $M$ increases. Results for this method are only shown for $M \le 4096$, because using it for larger values would have been prohibitively costly. For the largest values of $M$\kern -0.08333em, the cost of the residual method is almost the same for $k=20$ and $k=40$, because it is dominated by the computations needed to form and invert the ${\bm{M}}_{gg}$ matrices.
In contrast, both the score and jackknife methods become faster as $M$ increases and $G$ consequently decreases, except that, when $k=40$, they are both a bit slower for $M=64K$ than for $M=32K$. This probably occurs because of cache congestion. The jackknife method is always quicker than the score method. For small values of $M$\kern -0.08333em, it seems to be faster by a factor of about 12 when $k=20$ and by a factor of about 26 when $k=40$. However, the advantage of the jackknife method gradually diminishes as $M$ increases. When $M=64K$, so that there are only 16 clusters, the jackknife method is only slightly faster.
It is easy to see that the jackknife method will have a big advantage over the residual method whenever cluster sizes vary much, even if most of them are very small. Imagine a sample with, say, 1000 equal-sized clusters and $M=5$. For such a sample, the residual and jackknife methods will perform about the same. Suppose we then merge 100 of the tiny clusters into one large cluster with 500 observations. Doing this will reduce the cost of the jackknife method slightly, but it will greatly increase the cost of the residual method. Indeed, when there is even a single very large cluster, the latter inevitably becomes extremely slow.
Based on these results, the jackknife method for computing CV$_{\kern -0.08333em3}$ is clearly the procedure of choice unless all clusters are tiny (say, $N_g \leq 5$ for all $g$). For datasets with large clusters, an efficient implementation of this method (such as the one provided by the summclust package mentioned in (ref)), can compute both the OLS estimates and the CV$_{\kern -0.08333em3}$ variance matrix in roughly the same amount of time as a reasonably fast program for the OLS estimates alone.
The existing WCR bootstrap is based on CV$_{\kern -0.08333em1}$ standard errors and the restricted empirical score vectors defined in \hyperref[{rescore}]{\tagform@{\ref*{rescore}}} below. Henceforth, we will refer to this as the classic WCR bootstrap, or WCR-C. It often works well, but not always. We therefore propose three new versions of the WCR bootstrap, along with three corresponding versions of the WCU bootstrap. These are based on two distinct modifications. One involves replacing CV$_{\kern -0.08333em1}$ by CV$_{\kern -0.08333em3}$. The other involves modifying the scores used in the bootstrap DGP, in the hope that the modified bootstrap DGP will provide a better approximation to the unknown process that actually generated the data.
We first discuss the bootstrap DGPs for all versions of the wild cluster bootstrap, expressing them in terms of scores instead of observations. This approach is intuitive and computationally attractive RMNW,JGM-fast. In terms of the $G$ score vectors, a generic wild cluster bootstrap DGP is
where $b$ indexes bootstrap samples, $v_g^{*b}$ is a random variate with mean 0 and variance 1, and the $\ddot{\bm{s}}_g$ are empirical score vectors to be discussed below. In most cases, it seems to be best to generate the $v_g^{*b}$ using the Rademacher distribution, which takes the values 1 and $-1$ with equal probabilities DF_2008,DMN_2019. However, since the number of possible Rademacher bootstrap samples that are distinct from the original sample is only $2^G - 1$, it is better to use a distribution with more mass points, such as the six-point distribution proposed in Webb_2022, when $G$ is less than about 12.
The vector $\ddot{\bm{s}}_g$ in \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} is an empirical score vector for the \th{g} cluster. For the WCU-C bootstrap, it is simply the unrestricted empirical score vector $\hat{\bm{s}}_g={\bm{X}}_g^\top\hat{\bm{u}}_g$. For the WCR-C bootstrap, it is the restricted empirical score vector $\tilde{\bm{s}}_g$ defined as
where $\tilde{\bm\beta}$ is the vector of OLS estimates under the null hypothesis. Like $\hat{\bm\beta}$, $\tilde{\bm{s}}_g$ is a $k$-vector, even though some elements of $\tilde{\bm\beta}$ may equal zero or satisfy other linear restrictions. The bootstrap DGP \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} looks very much like the one for the wild score cluster bootstrap for nonlinear models proposed in KS_2012a. In the context of \hyperref[{eq:lrmodel}]{\tagform@{\ref*{eq:lrmodel}}}, however, it is just a different way of writing the bootstrap DGP for the wild cluster bootstrap.
In order to calculate a bootstrap $P$ value or a bootstrap confidence interval, we need to compute $B$ bootstrap test statistics indexed by $b$. These depend only on the bootstrap scores in \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} and the matrix $({\bm{X}}^\top{\bm{X}})^{-1}$\kern -0.08333em. For each bootstrap sample, we use ${\bm{s}}_g^{*b}$ to obtain a bootstrap estimate, not of ${\bm\beta}$ itself, but of the vector ${\bm{\delta}} = {\bm\beta} - \ddot{\bm\beta}$, where $\ddot{\bm\beta} = \tilde{\bm\beta}$ for the WCR bootstrap and $\ddot{\bm\beta} = \hat{\bm\beta}$ for the WCU bootstrap. This estimate is simply
where ${\bm{s}}^{*b} = \sum_{g=1}^G{\bm{s}}_g^{*b}$. When $v_g^{*b}=1$ for all $g$, the bootstrap sample is the same as the original sample. In this very special case, $\hat{\bm{\delta}}^{*b} = {\bm{0}}$ for the WCU bootstrap, and $\hat{\bm{\delta}}^{*b} = \hat{\bm\beta} - \tilde{\bm\beta}$ for the WCR bootstrap.
If we are testing the hypothesis that $\beta_j=0$, where $\beta_j$ is an element of ${\bm\beta}$, then we just need to multiply the \th{j} row of $({\bm{X}}^\top{\bm{X}})^{-1}$ by ${\bm{s}}^{*b}$ in order to obtain $\hat\delta_j^{*b}$, the \th{j} element of ${\bm{\delta}}^{*b}$\kern -0.08333em. The bootstrap $t$-statistic is then equal to
where ${\rm se}(\cdot)$ denotes the standard error formula used to obtain $t_j$, the original $t$-statistic. We automatically get the correct numerator, which is $\hat\beta_j^{*b}$ for the WCR bootstrap, since $\ddot{\bm\beta}=\tilde{\bm\beta}$, and $\hat\beta_j^{*b} - \hat\beta_j$ for the WCU bootstrap, since $\ddot{\bm\beta}=\hat{\bm\beta}$. As usual, a symmetric bootstrap $P$ value is then given by
where ${\mathbb I}(\cdot)$ denotes the indicator function. The bootstrap $P$ value in \hyperref[{eq:bootp}]{\tagform@{\ref*{eq:bootp}}} is simply the fraction of the bootstrap samples for which $|t_j^{*b}|$ is more extreme than $|t_j|$. The value of $B$ should be chosen so that $\alpha(B+1)$ is an integer, where $\alpha$ is the level of the test RM_2007. It is common to use $B=999$, but $B=9,\kern -0.08333em999$ and (when feasible) $B=99,\kern -0.08333em999$ are better choices.
In the classic versions of the wild cluster bootstrap, the standard error formula in \hyperref[{boottj}]{\tagform@{\ref*{boottj}}} is ${\rm se}_1(\cdot)$, the square root of the \th{j} diagonal element of CV$_{\kern -0.08333em1}$. But the results in (ref) make it equally feasible to use standard errors based on CV$_{\kern -0.08333em3}$, even in large samples. This gives us new versions of both the WCR and WCU bootstraps, which we will refer to as WCR-V and WCU-V, because only the variance matrices have changed. The bootstrap standard errors can be calculated without computing an entire variance matrix for each bootstrap sample. For example, the CV$_{\kern -0.08333em3}$ standard error of $\hat\delta_j^{*b}$ is just
where $\hat\delta_{j(g)}^{*b}$ is the \th{j} element of the vector
Only $\hat\delta_j^{*b}$ and the $\hat\delta_{j(g)}^{*b}$ need to be computed for each bootstrap sample. In \hyperref[{bootbeta}]{\tagform@{\ref*{bootbeta}}} and \hyperref[{eq:deloneboot}]{\tagform@{\ref*{eq:deloneboot}}}, the first terms are invariant across bootstrap samples and only need to be computed once.
We now have two versions of the WCR bootstrap, WCR-C and WCR-V, and two versions of the WCU bootstrap, WCU-C and WCR-V. The two WCR bootstraps use the bootstrap DGP \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} with $\ddot{\bm{s}}_g = \tilde{\bm{s}}_g$, and the two WCU bootstraps use the bootstrap DGP \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} with $\ddot{\bm{s}}_g = \hat{\bm{s}}_g$. The “C” and “V” versions calculate both the actual and bootstrap test statistics using ${\rm se}_1(\cdot)$ and ${\rm se}_3(\cdot)$, respectively. These bootstrap methods use the restricted or unrestricted empirical scores in their raw form. But empirical scores differ from true scores, because residuals differ from disturbances. It therefore seems attractive to replace the empirical score vectors by modified score vectors that implicitly rescale the residuals on a cluster-by-cluster basis. This is analogous to methods discussed in DF_2008 and JGM_2013 for the ordinary wild bootstrap. However, quite a lot more algebra is needed.
For the WCU bootstrap, we can simply replace the vectors $\ddot{\bm{s}}_g$ in \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} with the modified empirical score vectors $\acute{\bm{s}}_g$ defined in \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}. Using \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}} is expensive for large clusters, but the result \hyperref[{eq:keyres}]{\tagform@{\ref*{eq:keyres}}} lets us compute $\acute{\bm{s}}_g$ very rapidly as
For large clusters, using \hyperref[{eq:bis3f}]{\tagform@{\ref*{eq:bis3f}}} to compute the $\acute{\bm{s}}_g$ is much faster than using \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}, but using \hyperref[{eq:jackequal}]{\tagform@{\ref*{eq:jackequal}}} is faster still; see (ref). This yields two new bootstrap methods, which we will refer to as WCU-S and WCU-B, respectively. The WCU-S bootstrap (S for score) employs the modified score vectors $\acute{\bm{s}}_g$ instead of $\hat{\bm{s}}_g$, but it uses the familiar ${\rm se}_1(\cdot)$ standard error. The \mbox{WCU-B} bootstrap (B for both) employs both the modified score vectors and the ${\rm se}_3(\cdot)$ standard error.
Finding the analogous versions of the WCR bootstrap takes a bit more work. We need to specify a restricted wild bootstrap DGP based on modified score vectors. Suppose the restrictions have the usual linear form, ${\bm{R}} {\bm\beta} = {\bm{r}}$, for a given matrix ${\bm{R}}$ and a given vector ${\bm{r}}$. We can write this equivalently in terms of free parameters, $\bm\phi$, as ${\bm\beta} = {\bm{H}}\bm\phi +{\bm{h}}$ for a given matrix ${\bm{H}}$ and a given vector ${\bm{h}}$. Then the modified score vectors are
which are the analogs of the $\acute{\bm{s}}_g$ from \hyperref[{eq:bis3}]{\tagform@{\ref*{eq:bis3}}}. Here $\tilde{\bm{M}}_{gg}$ is the \th{g} diagonal block of the projection matrix $\tilde{\bm{M}} = {\bf I} -\tilde{\bm{X}} (\tilde{\bm{X}}^\top\tilde{\bm{X}} )^{-1} \tilde{\bm{X}}^\top$, where $\tilde{\bm{X}} = {\bm{X}}{\bm{H}}$. However, evaluating \hyperref[{eq:bisr3}]{\tagform@{\ref*{eq:bisr3}}} is computationally infeasible when the clusters are not all small. We need to replace \hyperref[{eq:bisr3}]{\tagform@{\ref*{eq:bisr3}}} by something that is feasible for any sample size.
The first step is to compute $\tilde{\bm\beta} = {\bm{H}}\tilde{\bm\phi} +{\bm{h}}$, where $\tilde{\bm{y}} = {\bm{y}} - {\bm{X}}{\bm{h}}$ and $\tilde{\bm\phi} = (\tilde{\bm{X}}^\top\tilde{\bm{X}})^{-1}\tilde{\bm{X}}^\top\tilde{\bm{y}}$. The corresponding estimates when cluster $g$ is omitted are $\tilde{\bm\beta}^{(g)} = {\bm{H}}\tilde{\bm\phi}^{(g)} +{\bm{h}}$, where
Then it can be shown that
To see that \hyperref[{eq:bisr3}]{\tagform@{\ref*{eq:bisr3}}} and \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} are equal, note that the right-hand side of \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} is
where the equality uses the updating formula \hyperref[{eq:update}]{\tagform@{\ref*{eq:update}}} applied to $\tilde{\bm{X}}$\kern -0.08333em, $\tilde{\bm{X}}_g$, and $\tilde{\bm{M}}_{gg}^{-1}$. Then we use the fact that $\tilde{\bm\phi} = (\tilde{\bm{X}}^\top\tilde{\bm{X}})^{-1}\tilde{\bm{X}}^\top \tilde{\bm{y}}$ together with the relation $\tilde{\bm{X}}_g (\tilde{\bm{X}}^\top \tilde{\bm{X}})^{-1}\tilde{\bm{X}}_g^\top = \tilde{\bm{P}}_{gg} = {\bf I}-\tilde{\bm{M}}_{gg}$ to rewrite the last expression as
Replacing $\tilde{\bm{y}}_g$ by ${\bm{y}}_g -{\bm{X}}_g {\bm{h}}$ and $\tilde{\bm{X}}_g$ by ${\bm{X}}_g {\bm{H}}$, and using the fact that ${\bm{H}}\tilde{\bm\phi} = \tilde{\bm\beta} -{\bm{h}}$, the right-hand side of \hyperref[{eq:intexp}]{\tagform@{\ref*{eq:intexp}}} equals \hyperref[{eq:bisr3}]{\tagform@{\ref*{eq:bisr3}}}.
An important special case is the restriction that $\beta_k = 0$. This is obtained by setting ${\bm{R}} = ( 0,\ldots ,0,1 )$ and ${\bm{r}} =0$, or, equivalently, ${\bm{H}} = ({\bf I}_{k-1}, {\bm{0}})^\top$ and ${\bm{h}} = {\bm{0}}$. In this case, we find that $\tilde{\bm{X}} = {\bm{X}}_1$, which contains the first $k-1$ columns of ${\bm{X}}$\kern -0.08333em, and $\tilde{\bm\phi}=\tilde{\bm\beta}_1 = ({\bm{X}}_1^\top{\bm{X}}_1)^{-1}{\bm{X}}_1^\top{\bm{y}}$. The corresponding estimates when each cluster is omitted in turn are
where ${\bm{X}}_{1g}$ contains the first $k-1$ columns of ${\bm{X}}_g$. Then \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} reduces to
Exactly the same arguments that led to \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} can be applied to the modified unrestricted empirical scores, giving us
Either \hyperref[{eq:jackequal}]{\tagform@{\ref*{eq:jackequal}}} or \hyperref[{eq:jackequ}]{\tagform@{\ref*{eq:jackequ}}} can be used to compute the $\acute{\bm{s}}_g$, and both are computationally attractive. However, in situations where both $\dot{\bm{s}}_g$ and $\acute{\bm{s}}_g$ need to be computed, \hyperref[{eq:jackequ}]{\tagform@{\ref*{eq:jackequ}}} may offer some programming advantages relative to \hyperref[{eq:jackequal}]{\tagform@{\ref*{eq:jackequal}}} due to its similarity to \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}}.
The scalar factors in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}} and \hyperref[{eq:CV3}]{\tagform@{\ref*{eq:CV3}}} do not appear in the bootstrap DGPs that correspond to them, because rescaling all the bootstrap scores by the same factor has no impact on the bootstrap $t$-statistics. From \hyperref[{bootbeta}]{\tagform@{\ref*{bootbeta}}} and \hyperref[{eq:deloneboot}]{\tagform@{\ref*{eq:deloneboot}}}, multiplying all the ${\bm{s}}_g^{*b}$ by a scalar $C$ simply makes $\hat{\bm{\delta}}^{*b}$ and all the $\hat{\bm{\delta}}_{(g)}^{*b}$ larger by a factor of $C$. This also makes the empirical scores for every bootstrap sample larger by the same factor. Therefore, from \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}}, \hyperref[{eq:CV2}]{\tagform@{\ref*{eq:CV2}}}, and \hyperref[{eq:CV3}]{\textup{\tagform@{\ref*{eq:CV3}}}}, the variance matrices become larger by a factor of $C^{\kern 0.04167em2}$ and the standard errors by a factor of $C$. The factors of $C$ in the numerator and denominator of $t_j^{*b}$ cancel out, leaving the bootstrap $t$-statistics unchanged.
However, no cancellation occurs for bootstrap tests of $\beta_j=0$ based directly on $\hat\beta_j$ and its bootstrap analog of $\hat\delta_j^{*b}$. In this case, multiplying the right-hand side of \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} by the square root of $G(N-1)/((G-1)(N-k))$ for methods that use CV$_{\kern -0.08333em1}$ and by the square root of $(G-1)/G$ for methods that use CV$_{\kern -0.08333em3}$ should improve the correspondence between the bootstrap DGP and the unknown true DGP. The usual theory of higher-order refinements for the bootstrap suggests that it is generally better to studentize Hall_1992. However, there may be cases in which unstudentized test statistics are of interest \citep*{CSS_2021}. But since we have eight studentized bootstrap methods to study, we do not consider unstudentized ones further.
To generate the transformed scores needed for the WCR/WCU-S and WCR/WCU-B bootstraps, \hyperref[{eq:jackequal}]{\tagform@{\ref*{eq:jackequal}}} and \hyperref[{eq:jackequ}]{\tagform@{\ref*{eq:jackequ}}} must be used for all $G$ clusters. In the event that $\hat{\bm\beta}^{(h)}$ and $\tilde{\bm\beta}^{(h)}$ cannot be calculated for cluster $h$, we have two choices. The simplest is to replace the inverses in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} and \hyperref[{eq:delrest}]{\tagform@{\ref*{eq:delrest}}} by generalized inverses. Alternatively, we could use $\hat{\bm{s}}_h$ instead of $\acute{\bm{s}}_h$ and $\tilde{\bm{s}}_h$ instead of $\dot{\bm{s}}_h$, along with the transformed scores for the remaining clusters. The latter would be appropriate if we have chosen to omit the problematic clusters when computing the cluster-jackknife variance matrix; see the discussion at the end of (ref).
(ref) provides a convenient summary of the eight wild cluster bootstrap methods that we have discussed. Conceptually, they differ along two dimensions. The horizontal dimension represents the way in which the standard errors for both the actual and bootstrap test statistics are calculated. The vertical dimension represents the score vectors used in the four versions of the bootstrap DGP \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}}. Note that the boottest and fwildclusterboot packages now provide fast implementations of the WCR/WCU-S bootstraps as well as the classic ones. This is possible because, in contrast to the WCR/WCU-V and WCR/WCU-B bootstraps, the former do not involve any jackknife calculations for the bootstrap samples. Once the transformed scores have been computed, the fast bootstrap algorithm proposed in RMNW applies directly to the \mbox{WCR/WCU-S} bootstraps.
It seems highly likely that all the methods discussed in this section are asymptotically valid, in the sense that, under suitable regularity conditions, the rejection frequencies for any test converge to the nominal level of the test as $G\to\infty$. Formal proofs could be obtained by modifying the arguments in DMN_2019. For the WCU bootstrap methods, the key fact is that the modified empirical score vectors $\acute{\bm{s}}_g$ computed using \hyperref[{eq:jackequal}]{\tagform@{\ref*{eq:jackequal}}} or \hyperref[{eq:jackequ}]{\tagform@{\ref*{eq:jackequ}}} are asymptotically equal to the ordinary empirical score vectors $\hat{\bm{s}}_g$. For the WCR bootstrap methods, the key fact is that the modified restricted empirical score vectors $\dot{\bm{s}}_g$ defined in \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} are asymptotically equal to the ordinary restricted empirical score vectors $\tilde{\bm{s}}_g$ in \hyperref[{rescore}]{\tagform@{\ref*{rescore}}}.
Simulation results in MW-JAE,MW-EJ, Brewer_2018, DMN_2019, JGM-fast, and several other papers have shown that the reliability of bootstrap and asymptotic methods for cluster-robust inference depends heavily on the number of clusters, the extent to which cluster sizes vary, and (in the case of treatment effects) both the number of treated clusters and their sizes. Many of our experiments therefore focus on these features.
The model we consider is
where the $u_{gi}$ are generated by a normal random-effects model with intra-cluster correlation $\rho$. The way in which the $k-1$ non-constant regressors are generated varies across the experiments. The hypothesis to be tested is that $\beta_k=0$.
In most of our experiments, there are $N=400G$ observations, which are divided among the $G$ clusters using the formula
where $[x]$ means the integer part of $x$. The value of $N_G$ is then set to $N - \sum_{g=1}^{G-1} N_g$. The key parameter here is $\gamma$, which determines how uneven the cluster sizes are. When $\gamma=0$ and $N/G$ is an integer, \hyperref[{eq:gameq}]{\tagform@{\ref*{eq:gameq}}} implies that $N_g = N/G$ for all $g$. For $\gamma>0$, cluster sizes vary more and more as $\gamma$ increases. The largest value of $\gamma$ that we use is 4. In that case, when $G=24$ and $N=9600$, the largest cluster (1513 observations) is about 47 times as large as the smallest cluster (32 observations). In contrast, when $\gamma=2$, the largest cluster (899 observations) is just under seven times as large as the smallest (130 observations).
The sample sizes that we employ are unusually large for experiments of this type. Since cluster-robust inference is often used with samples that have hundreds of thousands or even millions of observations, we want our results to apply to such cases. In preliminary experiments, we found that the results tended to change slightly, but systematically, as small values of $N/G$ were increased. Results for $N/G>400$ are very similar to ones for $N/G=400$, so we use 400 in all the experiments based on \hyperref[{eq:gameq}]{\tagform@{\ref*{eq:gameq}}}. Because the bootstrap samples are generated using scores, the cost of the experiments increases much less than proportionally with $N/G$.
All experiments use $400,\kern -0.08333em000$ replications. This number is so large that experimental randomness is negligible. The most important determinant of computational cost is $k$, the number of regressors. As can be seen from \hyperref[{eq:wcscore}]{\tagform@{\ref*{eq:wcscore}}} and \hyperref[{eq:jackeqr}]{\tagform@{\ref*{eq:jackeqr}}} or \hyperref[{eq:jackequ}]{\tagform@{\ref*{eq:jackequ}}}, generating each bootstrap sample involves $O(k^2G)$ operations. So does calculating the test statistics using either CV$_{\kern -0.08333em1}$ or CV$_{\kern -0.08333em3}$. Thus the experiments can be somewhat costly when $k$ is large.\footnote{The method of RMNW, which can only be used for the WCR/WCU-C and WCR/WCU-S bootstraps, is usually less expensive when $k$ is not small, but our programs did not use it.} Nevertheless, many of our experiments involve $k\ge10$. We do this because results in JGM-fast suggest that the performance of many methods of inference deteriorates as $k$ increases. Previous Monte Carlo experiments, which often use $k\le3$, may therefore have tended to give too optimistic a picture.
It might seem that substantial savings could be achieved by partialing out all regressors except the one(s) of interest prior to performing the bootstrap. However, this trick only works in certain special cases. For methods based on the jackknife, it is easy to see the problem. If we were to partial out some of the regressors prior to computing the delete\kern 0.04167em-one\kern 0.04167em-cluster estimates in \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}}, then the computed $\hat{\bm\beta}^{(g)}$ would depend on the values of the partialed-out regressors for the full sample, including those in the \th{g} cluster which was supposed to be deleted. Consequently, the values of the delete\kern 0.04167em-one\kern 0.04167em-cluster estimates would be incorrect if we partialed out any regressor that affects more than one cluster (such as industry-level fixed effects with firm-level clustering).
An important exception is when the regressors that are partialed out are cluster fixed effects or fixed effects at a finer level (such as firm-level fixed effects with industry-level clustering), because each of them affects only some or all of the observations within a single cluster. In fact, it is essential to partial out fixed effects of this type if using a generalized inverse is to be avoided.
The experiments in this subsection deal with rejection frequencies under the null hypothesis. We consider both asymptotic tests based on the $t(G-1)$ distribution and the wild cluster bootstrap tests listed in (ref).
The experiments in (ref) focus on variation in cluster sizes. There are always 9600 observations, 24 clusters, and 10 regressors. Cluster sizes vary according to \hyperref[{eq:gameq}]{\tagform@{\ref*{eq:gameq}}}. Regressors 2 through $k-1$ in \hyperref[{eq:simmod}]{\tagform@{\ref*{eq:simmod}}} follow a normal random-effects model that yields intra-cluster correlations of 0.50. The test regressor either follows the same normal distribution as the others (in the three panels on the left), or a $\chi^2(1)$ distribution (in the three panels on the right). In the latter case, it is the square of a normally distributed random variable generated by the same random-effects model as the other regressors. The disturbances are also generated by such a model, but with $\rho=0.10$. We focus on rejection frequencies for a test that $\beta_k=0$ at the 5% level.
The results for asymptotic tests, based on the $t(23)$ distribution and shown in Panels (a) and (b), are striking. Note that a square\kern 0.04167em-root transformation has been applied to the vertical axis to prevent these panels from being too tall. Tests based on CV$_{\kern -0.08333em1}$ over-reject substantially. The extent of the over-rejection increases with $\gamma$, and, except for $\gamma=4$, it is more severe in Panel (b) than in Panel (a). A regressor that follows the $\chi^2(1)$ distribution necessarily has some extreme values. These become points of high leverage, which makes inference more difficult in Panel (b).
Although tests based on CV$_{\kern -0.08333em2}$ always reject considerably less often than ones based on CV$_{\kern -0.08333em1}$, they also over-reject significantly and to an extent that increases with $\gamma$. In contrast, tests based on CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ either under-reject slightly all the time, in Panel (a), or under-reject very slightly for larger values of $\gamma$, in Panel (b). The results for CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ are extremely similar. The latter always rejects more often than the former, because the difference between \hyperref[{eq:jack}]{\tagform@{\ref*{eq:jack}}} and \hyperref[{eq:jack3}]{\tagform@{\ref*{eq:jack3}}} is the positive semi-definite matrix $((G-1)/G) (\hat{\bm\beta}-\bar{\bm\beta})(\hat{\bm\beta}-\bar{\bm\beta})^\top$. Since CV$_{\kern -0.08333em3}$ tends to under-reject slightly in (ref), it might seem that CV$_{\kern -0.08333em3{\rm J}}$ is to be preferred. However, as we shall see, there are also many cases in which CV$_{\kern -0.08333em3}$ over-rejects, and CV$_{\kern -0.08333em3{\rm J}}$ therefore over-rejects slightly more. In practice, it would be perfectly reasonable to report either CV$_{\kern -0.08333em3}$ or CV$_{\kern -0.08333em3{\rm J}}$. We never encountered a case in which it made any real difference.
The results for the WCR bootstrap tests, shown in Panels (c) and (d), are surprising. In the past, WCR-C has been the only variant of the WCR bootstrap, and numerous Monte Carlo experiments have suggested that it is the procedure of choice. But WCR-B performs notably better than WCR-C for every value of $\gamma$, and both WCR-V and WCR-S perform better still. Note that, although these two procedures perform almost the same here, this is not true in general. Oddly, all the WCR procedures perform better in Panel (d), where the test regressor is highly skewed, than they do in Panel (c), where it is Gaussian. The rather mediocre performance of WCR-C must be due, at least in part, to the fact that $k=10$, which is a larger number than has been used in most previous experiments; see (ref) below.
Some of the results for the WCU bootstrap tests, shown in Panels (e) and (f), are also surprising. It is not a surprise that WCU-C rejects more often than WCR-C or that its performance is much worse in Panel (f) than in Panel (e). However, the fact that the other three WCU procedures over-reject much less often than WCU-C may well be surprising. In both panels, WCU-B is clearly the procedure of choice. WCU-V and WCU-S perform much better than \mbox{WCU-C}, but worse than \mbox{WCU-B}. In Panels (c) and (d), the differences between \mbox{WCU-V} and \mbox{WCU-S} are small, but larger than the differences between \mbox{WCR-V} and \mbox{WCR-S}.
(ref) is similar to (ref), but the number of regressors $k$ is now on the horizontal axis, and $\gamma=2$. In Panels (a) and (b), CV$_{\kern -0.08333em1}$ over-rejects to an increasing extent as $k$ increases. So does CV$_{\kern -0.08333em2}$, although it always over-rejects considerably less than CV$_{\kern -0.08333em1}$. In contrast, CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ over-reject modestly for small values of $k$ and under-reject modestly for large ones.
Panels (c) and (d) look a lot like the same panels in (ref), even though what is on the horizontal axis is different. WCR-C performs quite well for very small values of $k$, but it over-rejects more and more severely as $k$ increases. WCR-B performs much better than WCR-C, but WCR-V and WCR-S perform even better. In Panel (d), where the test regressor is highly skewed, they both perform extremely well for all values of $k$.
Panels (e) and (f) also look a lot like the same panels in (ref). WCU-C performs quite poorly, over-rejecting more and more severely as $k$ increases. In contrast, WCU-B performs quite well in Panel (e) and fairly well in Panel (f), and there is no tendency for its performance to deteriorate as $k$ increases. As before, the two other bootstrap methods generally perform much better than WCU-C but slightly worse than WCU-B.
In the next set of experiments, we focus on what happens as $G$ increases. (ref) shows rejection frequencies as functions of $G$, which varies from 12 to 84 by 6, and implicitly also $N$\kern -0.08333em, since $N=400G$. In these experiments, $\gamma=2$ and $k=10$. We report results for only five methods, instead of twelve. We omit CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em2}$, because they never perform very well, and CV$_{\kern -0.08333em3{\rm J}}$ because it is almost identical to CV$_{\kern -0.08333em3}$. Among the restricted bootstrap methods, we report WCR-C, because it was until now the procedure of choice. We also report WCR-S and WCR-B, but we do not report WCR-V, because it yields results nearly identical to those of WCR-S and is harder to compute. Among the unrestricted bootstrap methods, we report only WCU-B, because it always seems to outperform the other WCU methods.
In Panel (a), using CV$_{\kern -0.08333em3}$ with the $t(G-1)$ distribution under-rejects quite noticeably for very small values of $G$, but it performs extremely well for $G \geq 30$. The bootstrap methods always over-reject, with WCR-C always the worst of them. For $G \geq 42$, however, all the bootstrap methods perform very well, with WCR-S the winner by a tiny margin.
Panel (b) is more interesting than Panel (a). The extreme skewness of the $\chi^2(1)$ regressor apparently affects the results quite a bit, even when $G=84$. Using CV$_{\kern -0.08333em3}$ with the $t(G-1)$ distribution under-rejects for small values of $G$ but over-rejects for larger values, where it is the worst method. Note that $G=24$, the value in (ref), is near where the curve for CV$_{\kern -0.08333em3}$ crosses the .05 line in (ref). The best method is WCR-S in every case. It performs remarkably well for $G \geq 30$. However, all three WCR methods perform well for the larger values of $G$. Indeed, by most standards, every method shown in Panel (b) of (ref) works very well, unless $G$ is less than about 30. For $G=84$, CV$_{\kern -0.08333em3}$ is the worst method, but even it rejects only 5.49% of the time. For comparison, CV$_{\kern -0.08333em1}$ rejects 9.04% of the time, and CV$_{\kern -0.08333em2}$ rejects 7.15%. The best method, WCR-S, rejects 4.97% of the time.
Many applications of cluster-robust inference involve treatment at the cluster level, and existing methods generally perform very poorly when either the number of treated clusters or the number of control clusters is small. Using CV$_{\kern -0.08333em1}$ with the $t(G-1)$ distribution or WCU-C leads to severe over-rejection, and using WCR-C leads to severe under-rejection MW-JAE, MW-EJ. Our next set of experiments therefore focuses on the model
where $x_g$ is a treatment dummy, ${\bm{Z}}_{gi}$ is a row vector of other regressors, and $u_{gi}$ is generated by a random-effects model with intra-cluster correlation $\rho$. The treatment dummy equals 1 for $G_1$ of the $G$ clusters and 0 for the remaining $G_0=G-G_1$. The clusters that are treated are chosen at random. The ${\bm{Z}}_{gi}$ consist of eight more dummy variables. For each of these variables and each cluster, a probability $\pi_g$ between 0.25 and 0.75 is chosen at random for each replication. Then each observation for that variable in that cluster equals 1 with probability $\pi_g$ and 0 otherwise. Thus all the regressors are dummies, which vary at the individual level in a way that varies across clusters.
(ref) shows rejection frequencies based on the $t(G-1)$ distribution. In the left-hand column, there are 12 clusters and 4800 observations. In the right-hand one, there are 24 clusters and 9600 observations. The value of $\gamma$ is 0 in the top row, 2 in the middle row, and 4 in the bottom row. The number of treated observations $G_1$ varies between 2 and $G-2$ on the horizontal axes. It would have been impossible to set $G_1=1$ or $G_1=G-1$, because CV$_{\kern -0.08333em2}$, CV$_{\kern -0.08333em3}$, and CV$_{\kern -0.08333em3{\rm J}}$ cannot be computed in those cases. This is obvious from \hyperref[{eq:delone}]{\tagform@{\ref*{eq:delone}}} for the jackknife\kern 0.04167em-based estimators. When the single treated cluster is omitted, the coefficient of interest in $\hat{\bm\beta}^{(g)}$ is not identified.
As previous work has shown, tests that use CV$_{\kern -0.08333em1}$ tend to over-reject severely when either $G_0$ or $G_1$ is small. This is evident in (ref). The over-rejection is worst in Panel (f), where both $\gamma$ and $G$ are largest. CV$_{\kern -0.08333em2}$ over-rejects less than CV$_{\kern -0.08333em1}$, but it still does not work very well, except perhaps for values of $G_1$ near $G/2$ when $\gamma=0$; see Panels (a) and (b). In contrast, CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$, which perform almost identically, are much less prone to over-reject than the other two CRVEs. They actually under-reject for values of $G_1$ fairly near $G/2$ when $\gamma=0$, and they perform very well for values of $G_1$ near $G/2$ when $\gamma=2$. Oddly, CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em3{\rm J}}$ over-reject less seriously for extreme values of $G_1$ when $\gamma$ is large than when $\gamma$ is small.
(ref) shows results for four bootstrap tests for the same set of experiments as in (ref). When $\gamma=0$, all three variants of the WCR bootstrap perform almost identically. However, as $\gamma$ increases, their performance starts to differ. WCR-S seems to reject least frequently, which is a good thing for intermediate values of $G_1$ and a bad thing for extreme values. In contrast, WCR-B under-rejects least severely for extreme values of $G_1$. However, for intermediate values, it over-rejects less than WCR-C but more than WCR-S.
The most surprising results in (ref) are the ones for the unrestricted wild bootstraps. We do not report results for WCU-C or WCU-S, because they would have required a much longer vertical axis. WCU-C rejects almost 28% of the time in its worst case ($G=24$, $G_1=2$, $\gamma=4$), and WCU-S rejects over 12% of the time in its worst case ($G=24$, $G_1=2$, $\gamma=0$). In contrast, WCU-B is arguably the best method overall when $G=12$, and it performs very well for intermediate values of $G_1$ when $G=24$. In addition, it never over-rejects as severely as CV$_{\kern -0.08333em3}$ for extreme values of $G_1$.
Simulations in DMN_2019 suggest that many methods work poorly when one cluster is much bigger than the others. Even when $\gamma=4$, the largest cluster in our experiments is never dramatically larger than all the rest, although this happens quite often in empirical work. For instance, more than half of all incorporations in the United States occur in Delaware HS_2020. This implies that studies of the effects of corporate governance based on changes in state laws, where standard errors are clustered by state of incorporation, are likely to encounter severe errors of inference. To investigate this phenomenon, we create artificial samples with 50 clusters based on data for incorporations by year and state from SW-data. There are $205,\kern -0.08333em566$ observations, of which $108,\kern -0.08333em538$, or 52.80%, are for Delaware. The second-largest cluster is Nevada, with $17,\kern -0.08333em010$ or 8.27%, and the smallest is Montana, with 101 or 0.05%.
We perform a set of experiments similar to the ones in (ref) using these artificial samples. There are 10 regressors, generated in the same way as before, with one exception. Because investigators are surely aware of whether or not the largest cluster (Delaware) is treated, it is always treated in our experiments. The other clusters to be treated (between 1 and 47 of them) are chosen at random. Because the largest cluster is always treated, the rejection frequencies are no longer the same for $G_1$ and $G-G_1$ treated clusters. However, since this is a pure treatment model, the results for $G_1$ treated clusters that include Delaware must be the same as the results for $G-G_1$ treated clusters that exclude Delaware.
The results in (ref) are striking. In Panel (a), using either CV$_1$ or CV$_2$ leads to over-rejection that varies between severe and extreme. Using CV$_3$ and CV$_{3{\rm J}}$ also leads to over-rejection, but it is much less severe. For between 20 and 41 treated clusters, rejection frequencies are less than 0.07. In Panel (b), WCU-C over-rejects severely, and WCR-C can either over-reject or under-reject, often severely. In contrast, our new bootstrap methods work remarkably well. The best of them is WCU-B, which always rejects less than 9% of the time and sometimes rejects just about 5% of the time. WCR-S and WCR-B also perform much better than WCR-C, except when $G_1$ is very large, in which case they under-reject severely.
Even though it is based on real data, the distribution of cluster sizes in the experiments of (ref) is very extreme. The performance of CV$_3$ and three of our new bootstrap methods is far from perfect, but it is generally very much better than that of existing methods. Thus jackknife\kern 0.04167em-based methods seem to be remarkably robust to heterogeneity in cluster sizes.
It is natural to worry that a new test may be less powerful than existing tests, especially when it performs much better under the null hypothesis. In this section, we therefore investigate test power. Studying power is tricky, because it is unreasonable to compare tests that have noticeably different rejection frequencies under the null. If, for example, an asymptotic test rejects 15% of the time under the null and a bootstrap test rejects 6% of the time, then we would expect the asymptotic test to have substantially more power than the bootstrap test. But the additional power may be entirely spurious, simply reflecting the finite\kern 0.04167em-sample over-rejection by the former.
One way to compare tests with different rejection frequencies under the null is to “size\kern 0.04167em-adjust” them. But this approach has two serious conceptual difficulties. First, size\kern 0.04167em-adjusted tests are infeasible. What do we learn by comparing tests that cannot actually be performed? Second, there are often many ways to size\kern 0.04167em-adjust a given test, and they may yield quite different results. The idea of size\kern 0.04167em-adjustment is to base rejection frequencies for tests under the alternative on critical values calculated by simulation under the null. But, in general, there exists an infinite number of DGPs that satisfy the null hypothesis. If they all yield the same critical values, then there is no problem. But if they yield different critical values, as will often be the case, then we have to choose which null DGP to use. It seems natural to make the null DGP used for critical values as close as possible to the alternative DGP. DM_2006 suggests a particular way of doing this, based on the Kullback-Leibler information criterion, but this approach means using a different critical value for each set of values of the parameters under test.
To avoid the difficulties just discussed, we focus on four cases where the tests of interest all perform quite well under the null. They are treatment experiments similar to the ones in (ref), with $G=24$, $N=9600$, and $k=5$. In Panels (a) and (b), $G_1=12$, so that precisely half the clusters are treated. In Panels (c) and (d), $G_1=6$, so that the effects of having few treated clusters are apparent but not severe. In order to avoid excessive power loss, we use $B=999$ for the bootstrap tests. We use $k=5$ instead of $k=10$ partly to reduce computational cost and partly to improve test performance under the null.
(ref) shows rejection frequencies as a function of $\beta_k$, the actual coefficient on the treatment dummy in \hyperref[{eq:treat}]{\tagform@{\ref*{eq:treat}}}, when the null hypothesis is that $\beta_k=0$. In Panels (a) and (c), $\gamma=0$, so that every cluster has exactly 400 observations. In Panel (a), the perfectly balanced case, all five power functions are visually indistinguishable. In Panel (c), where only six clusters are treated, CV$_3$ has noticeably more power than any of the bootstrap methods, which are all but identical.
In Panels (b) and (d), cluster sizes vary from 32 to 1513. All tests are now substantially less powerful than in Panels (a) and (c), because, whenever there is intra-cluster correlation, the information content of a sample declines as the cluster sizes become more variable. The most striking result in both panels is that WCU-B has noticeably less power than any of the other methods. This is especially true in Panel (d), where WCU-B over-rejects modestly under the null but becomes by far the least powerful method for larger values of $\beta_k$.
The pattern for CV$_{\kern -0.08333em3}$ is similar but much less pronounced. Under the null hypothesis, it over-rejects slightly under the null in Panel (b) and noticeably in Panel (d), with rejection frequencies of 0.0612 and 0.0775, respectively. But for large enough values of $\beta_k$, it has less power than WCR-C and WCR-S, especially in Panel (d). The latter two methods also have slightly more power than WCR-B in Panel (b) and noticeably more in Panel (d) for large values of $\beta_k$. Interestingly, WCR-V, which for clarity is not shown in the figure, has somewhat less power than either WCR-C or WCR-S in Panels (b) and (d), where cluster sizes vary a lot. In contrast, it is almost indistinguishable from both these methods in Panels (a) and (c), where cluster sizes are constant.
Based on these results, the procedure of choice appears to be WCR-S. For larger values of $\beta_k$, it is always one of the two most powerful tests. WCR-C has similar power, and it also works well under the null in these experiments, but it is much more prone to over-reject than WCR-S in (ref). Happily, WCR-S is already available in computationally efficient packages for Stata, R, and \texttt{Python}; see (ref).
Cluster-robust standard errors and bootstrap methods are often used to form confidence intervals. Although we do not perform any Monte Carlo experiments explicitly to study the properties of confidence intervals, these can be inferred from (ref) and the results in (ref). Most confidence intervals are implicitly or explicitly obtained by inverting a hypothesis test. When such a test has approximately the correct rejection frequency, the resulting confidence interval must have approximately correct coverage. Similarly, when such a test has high power, the resulting confidence interval must be relatively short.
In many of the experiments in (ref), tests based on CV$_{\kern -0.08333em3}$ and the $t(G-1)$ distribution are much less prone to over-reject than tests based on CV$_{\kern -0.08333em1}$. This suggests that the coverage of confidence intervals based on CV$_{\kern -0.08333em3}$ standard errors will often be much better than the coverage of ones based on CV$_{\kern -0.08333em1}$ standard errors. Even more reliable intervals may often (but not always) be obtained by using the WCR-S or WCR-B bootstraps, which perform much better than the classic WCR-C bootstrap in many cases. The WCU-B bootstrap also performs well in many cases under the null, but the results in Panels (b) and (d) of (ref) suggest that, when cluster sizes vary a lot, intervals based on it may be longer than ones based on WCR-B, which in turn may be slightly longer than ones based on WCR-S.
The WCR-S bootstrap has excellent performance in many of the experiments of (ref), seems to have slightly better power than WCR-B in Panels (b) and (d) of (ref), and is easy to compute. Therefore, we tentatively recommend that confidence intervals should be obtained by inverting WCR-S bootstrap tests. However, using CV$_{\kern -0.08333em3}$ standard errors and the $t(G-1)$ distribution would often lead to very similar intervals.
Of course, it is easier to obtain a confidence interval by using a standard error and the $t(G-1)$ distribution than by inverting a bootstrap test, and it is easier to invert any form of WCU bootstrap test than any form of WCR bootstrap test. However, the computational cost of inverting WCR bootstrap tests can be remarkably small, even for very large samples; see RMNW and JGM-fast.
In this section, we consider three empirical examples. These suggest that the new bootstrap procedures proposed in (ref) may sometimes yield results very similar to those from the existing WCR-C and WCU-C procedures, but they may also yield results which differ noticeably from those and from each other.
Our first example is based on MNW-guide. It exploits differences in the minimum wage across states and years to estimate the impact of minimum wages on hours worked for teenagers.
The data on hours at the individual level from the American Community Survey (ACS) are obtained from IPUMS IPUMS_2020 and cover the years 2005\kern 0.04167em--2019. The minimum wage data come from Neumark_2019 and are collapsed to state\kern 0.04167em-year averages to match the ACS frequency. We restrict attention to teenagers aged 16\kern 0.04167em--19, keeping only individuals who are children of the respondent to the survey and who have never been married. We drop individuals who had completed one year of college by age 16 and those reporting in excess of 60 hours usually worked per week. We also restrict attention to individuals who identify as either black or white. There are $492,\kern -0.08333em827$ observations in 51 clusters, which correspond to all 50 states plus the District of Columbia.
The model we estimate is
where $y_{ist}$ is usual hours worked per week for individual $i$. The parameter of interest is $\beta$, which is the coefficient on $\textrm{mw}_{\kern -0.08333em st}$, the minimum wage in state $s$ at time $t$. The row vector ${\bm{Z}}_{ist}$ collects a large set of individual-level controls, including race, gender, age, and education. There are also state and year fixed effects, denoted by $\delta_s$ and $\eta_t$, respectively.
As MNW-guide discusses, clustering could in principle be done at several different levels. However, the one that is most appealing and seems to be supported by the data is clustering at the state level. The 51 clusters vary considerably in size. The smallest has 258 observations, and the largest has $35,\kern -0.08333em995$. The ratio of these numbers is more than twice as large as for $\gamma=4$ in the experiments of (ref). The mean number of observations per cluster is $9,\kern -0.08333em663$, and the median is $7,\kern -0.08333em082$. This suggests that inference based on CV$_{\kern -0.08333em1}$ and the $t(50)$ distribution may not be reliable. Other measures of cluster heterogeneity, which are discussed in the original paper, lead to the same conclusion.
(ref) presents our key results. As expected, the CV$_{\kern -0.08333em3}$ $t$-statistic is somewhat smaller than the CV$_{\kern -0.08333em1}$ $t$-statistic, and the $P$ value based on the $t(50)$ distribution is therefore somewhat larger. The four WCR $P$ values are larger than either of them, but still below 0.05, and the four WCU $P$ values are notably smaller than the WCR ones. Because $B$ is so large (larger than really needed), the simulation standard errors for the WCR bootstrap $P$ values are about 0.0002.
Based on how similar the four WCR $P$ values are, and on how well many of the WCR methods perform in the experiments of (ref), we tentatively conclude that the “true” $P$ value for the test of $\beta=0$ is probably between 0.034 and 0.039. Thus the null hypothesis can safely be rejected at the 0.05 level but not at the 0.01 level.
The second example comes from \citet*{AMT_2022}. This paper examines the impact of political turnover on the quality of public services. Specifically, it examines several outcomes following close mayoral elections in Brazil. One of these outcomes is the test scores of fourth-grade students. The paper uses a regression discontinuity design to identify the treated and control municipalities, but it conducts the analysis using OLS. We replicate one such regression, found in Table 3, Column 5 of the original paper:
The dependent variable is the test score one year after an election. $\textrm{IVM}_{mt}$ is the incumbent vote margin in the close election which occurs in year $t$. Accordingly, the treatment variable is ${\mathbb I}(\textrm{IVM}_{mt}<0)$, which equals 1 when the incumbent party loses the election and a turnover occurs, and the coefficient of interest is $\beta$. This regression is estimated using a sample which is determined by a selected bandwidth. While the paper considers several bandwidths, we focus on the bandwidth 0.110, as this results in the largest sample.
The paper clusters the standard errors at the municipality level. Since there are 2101 municipalities, many of them located close to each other, it seems possible that this level of clustering is too fine. We therefore consider state-level clustering. However, there are only 26 states in Brazil, and they vary in size from 420 to $64,\kern -0.08333em953$ with partial leverages from 0.000234 to 0.179318 MNW-influence. With this much heterogeneity across clusters, relying on CV$_{\kern -0.08333em1}$ may be risky.
(ref) presents our key results. As expected, the CV$_{\kern -0.08333em1}$ standard error for clustering by state is smaller than the CV$_{\kern -0.08333em3}$ standard error. Contrary to our expectations, however, both are a bit smaller than the CV$_{\kern -0.08333em1}$ standard error for clustering by municipality. The four WCR $P$ values are similar to each other and to the $P$ value based on the CV$_{\kern -0.08333em1}$ $t$-statistic and the $t(25)$ distribution. Surprisingly, the four WCU $P$ values are noticeably larger than the WCR ones. Nevertheless, since every test rejects at the 0.05 level, there is evidence against the null hypothesis.
The third example is taken from Xu_2018, which explores the effect of patronage in the colonial era of Britain on the appointment of governors to colonies. Part of the analysis examines whether the extent to which the current secretary of state and a governor are “connected” led to more desirable colony postings. We replicate the results of one such regression, found in Table 3, Column 3 of the original paper:
Here $\textrm{log(revenue)}_{ist} $ is the initial revenue for colony $s$ when governor $i$ was appointed in year $t$. The main variable of interest is $\textrm{connected}_{it}$, which is a binary variable set equal to 1 when the governor and the secretary share connections such as having attended the same elite boarding school, or Oxford or Cambridge, or both being in the aristocracy, or having shared ancestry. The variable $\textrm{served}_{it}$ is the number colonies in which the governor has served up to the year of appointment. The regression also has fixed effects for governors ($\gamma_i$), years ($\tau_t$), and the duration of the governorship ($\delta_{it}$).
The paper clusters the standard errors at the bilateral pair (or dyad) level between the secretary of state and the governor. However, the dependent variable is observed at multiple times for each colony, so it seems likely that there would be dependence across observations for the same colony. The regression does not include colony fixed effects, which would have reduced this dependence, because, with so many other fixed effects, it was impossible to include them. Thus, it seems plausible that the standard errors should be clustered at the colony level instead of the dyad level, and we investigate this approach. Switching from dyadic clustering to clustering by colony actually reduces the CV$_{\kern -0.08333em1}$ standard error. However, even though there are 70 colonies, they are quite unbalanced; the number of observations per colony ranges from 4 to 104. Partial leverages also vary greatly, and they seem to be roughly proportional to cluster sizes. Perhaps in consequence, the CV$_{\kern -0.08333em3}$ standard error is 47% larger than the CV$_{\kern -0.08333em1}$ one.
In view of the dramatic difference between the CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ $t$-statistics, the various wild bootstrap methods provide valuable information. The bootstrap $P$ values are all somewhat larger than the one for the CV$_{\kern -0.08333em1}$ $t$-statistic based on the $t(69)$ distribution, but they are all much smaller than the corresponding one for the CV$_{\kern -0.08333em3}$ $t$-statistic. The smallest bootstrap $P$ value is the one for the classic WCR-C method. At $0.0535$, it is not much larger than the one based on the $t(69)$ distribution. Surprisingly, every WCU $P$ value is larger than the corresponding WCR $P$ value.
This example is deliberately extreme, because the number of regressors (573) is unusually large relative to the number of observations (3510). Perhaps in consequence, the full coefficient vectors $\hat{\bm\beta}^{(g)}$ are not identified for 61 out of the 70 clusters. However, since the $\hat\beta_1^{(g)}$ coefficients are always identified, we used a generalized inverse to compute both CV$_{\kern -0.08333em3}$ and the bootstrap DGPs for the WCR/WCU-S and WCR/WCU-B bootstraps. The alternative approach of trying to estimate a variance matrix based on only 9 out of 70 clusters seems very dubious, and it yields an implausibly small standard error of just 0.0379. However, the large number of singularities may explain why the CV$_{\kern -0.08333em3}$ and CV$_{\kern -0.08333em1}$ standard errors differ as much as they do.
Because $k$ is so large in this example, we suspect that the $t$-test based on CV$_{\kern -0.08333em3}$ may be prone to under-reject, and that both the $t$-test based on CV$_{\kern -0.08333em1}$ and the WCR-C bootstrap test may be prone to over-reject; see (ref). Nevertheless, the $P$ values for the new WCR bootstrap methods are only modestly larger than the WCR-C $P$ value. The fact that all the bootstrap $P$ values lie between 0.0535 and 0.0738 suggests that the “true” $P$ value probably also lies within, or at least not too far outside, this interval. We conclude that there seems to be only weak evidence against the null hypothesis.
The classic CV$_{\kern -0.08333em1}$ estimator given in \hyperref[{eq:CV1}]{\tagform@{\ref*{eq:CV1}}} is by far the most popular CRVE for linear regression models, but standard errors based on it are often much too small. The cluster jackknife estimator, often called CV$_{\kern -0.08333em3}$, has been known for many years but is much less widely used. In (ref), we discuss how to compute CV$_{\kern -0.08333em3}$ in a computationally efficient fashion. Except when all clusters are tiny, this is the fastest available method for computing it; see (ref). Inference based on CV$_{\kern -0.08333em3}$ and the Student's $t(G-1)$ distribution seems to be much more reliable than inference based on CV$_{\kern -0.08333em1}$ and that distribution; see (ref). This accords with theoretical results in Hansen-jack, which provides no simulations and cites the ones in this paper.
Although combining CV$_{\kern -0.08333em3}$ standard errors and the $t$ distribution often works well, it does not always do so. Bootstrap methods may well perform better, and they also provide a valuable robustness check. In (ref), we prove some simple, but by no means obvious, algebraic results about the relationship between cluster jackknife estimates and score vectors at the cluster level. These results allow us to obtain new and easy-to\kern 0.04167em-compute variants of the wild cluster bootstrap. These typically perform better than the classic variants, now called WCR-C and WCU-C. The eight new and existing variants are summarized in (ref). Of these, the ones that use CV$_{\kern -0.08333em1}$ together with modified bootstrap score vectors, called WCR-S and WCU-S, are particularly easy to compute. They are available in packages for Stata and R.
Prior to this paper, there were already quite a few methods for inference in linear regression models with clustered disturbances MNW-guide, and (ref) has added six new variants of the wild cluster bootstrap. Empiricists may reasonably ask what methods they should use in practice. As discussed in detail in MNW-guide,MNW-influence, the first thing to do is to investigate the clustering structure of the model and dataset. For instance, it is good practice to calculate the effective number of clusters CSS_2017 as well as various measures of leverage and influence at the cluster level MNW-influence. When these measures indicate that clusters are well-balanced, and the (effective) number of clusters is large (say, more than 100), then CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ should yield very similar standard errors. In such cases, it is probably safe to rely on CV$_{\kern -0.08333em3}$ standard errors together with the $t(G-1)$ distribution.
However, the number of clusters will often be much less than 100. Moreover, measures of cluster-level leverage and influence may indicate that clusters are not well-balanced. This can happen, for example, when cluster sizes vary a lot, when there are few treated clusters, or when the distributions of key regressors vary greatly across clusters. In such cases, CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ can yield quite different standard errors. Recall (ref), where the CV$_{\kern -0.08333em3}$ standard error is 47% larger than the CV$_{\kern -0.08333em1}$ standard error, even though there are 70 clusters. Whenever CV$_{\kern -0.08333em1}$ and CV$_{\kern -0.08333em3}$ differ substantially, bootstrap $P$ values or confidence intervals are likely to be more reliable than conventional ones based on either of those CRVEs, and it is probably a good idea to compute both WCR-C and WCR-S $P$ values.
In most cases, it is advisable to compute wild cluster bootstrap $P$ values and/or confidence intervals using at least $9,\kern -0.08333em999$ bootstrap samples. This is usually not computationally difficult. However, there might be exceptions when either the number of clusters or the number of regressors is unusually large. Of course, when the number of clusters is very large, the bootstrap will not be needed unless the clusters are severely unbalanced, but that can happen.
If we had to recommend just one method, it would be the WCR-S bootstrap proposed in (ref). This method uses ordinary CV$_{\kern -0.08333em1}$ standard errors, which makes it easy to compute, but the bootstrap DGP employs restricted scores that have been transformed using the cluster jackknife. In some of our experiments, the WCR-S bootstrap works substantially better than the classic (and popular) WCR-C bootstrap; see, in particular, (ref) and (ref). We generally do not recommend using the unrestricted wild cluster bootstrap, except perhaps as a robustness check or when it is desired to generate a large number of confidence intervals using just one set of bootstrap samples.
{0pt} \addcontentsline{toc}{section}{\refname}