Asymptotic results under multiway clustering

Asymptotic results under multiway clustering

Introduction

Taking into account dependence between observations is crucial for making correct inference. Common shocks tend to correlate observations positively, leading to overly optimistic inference when ignored bertrand2004. As a result, estimation of standard errors robust to clustering has become pervasive in applied economics. In particular, following the very influential work of cameron2011,\footnote{According to the Web of Science and Google Scholar, cameron2011 is the most cited paper in econometrics since 2009.} empirical studies now routinely report standard errors accounting for multiway clustering. Perhaps surprisingly however, econometric theory has lagged behind this practice. cameron2011 conduct simulations suggesting the validity of their method but neither prove that their variance estimators are consistent, nor that estimators of parameters of interest are themselves asymptotically normal. And more generally, there are still very few theoretical results under multiway clustering.

The goal of this paper is to fill this gap, by developing general tools for inference on linear but also nonlinear estimators with multiway clustering. We consider for that purpose a fairly general set-up including as particular cases one- and two-way clustering. To understand the key underlying restrictions, let us consider the example of two-way clustering where the first dimension is the sector of activity and the second is the area of residence, e.g. counties or states. Such an example would be appropriate when studying for instance individual wages. We index the two dimensions respectively by $j_1\in \{1,...,C_1\}$ and $j_2\in\{1,...,C_2\}$. We call a cell any pair $(j_1,j_2)$, corresponding therefore to a specific sector of activity and area of residence. Then two units in different cells $(j_1,j_2)$ and $(j'_1,j'_2)$ are assumed independent whenever $j_1\neq j'_1$ and $j_2\neq j'_2$. Otherwise they may be dependent in an unrestricted way. The idea behind is that units sharing at least one cluster may be affected by common shocks, e.g. sectorial shocks or local shocks in the previous example.

Following most of the literature, we consider an asymptotic framework where $\underline{C}=\min(C_1,C_2)$ tends to infinity.\footnote{A growing strand of the literature on one-way clustering has also considered fixed-$\underline{C}$ asymptotics, with cluster sizes tending to infinity. We refer in particular to donald2007inference, ibragimov2010, bester2011inference, ibragimov2016 and canay2018. To our knowledge, no paper has considered such a set-up with multiway clustering yet.} In particular, we allow for random, possibly unbounded, cell sizes. Cell sizes may also be correlated with the data themselves. These features are important to account for cluster heterogeneity, following the terminology of carter2017. Given this set-up, our first contribution is a general weak convergence result on the empirical process. To our knowledge, this weak convergence result is new, even under one-way clustering. When considering processes indexed by a finite class of functions, it is equivalent to a simple multivariate central limit theorem (CLT) on sample averages. But when considering infinite classes of functions, this result is also key for proving asymptotic normality of nonlinear estimators like GMM estimators or smooth functionals of the empirical cumulative distribution function (cdf). Also, up to moment restrictions that have to be slightly adapted, our conditions on the class of functions indexing the empirical process are the same as with i.i.d. data. This means that results on, e.g., GMM estimators already established for i.i.d. data can be extended directly to multiway clustering.

Then, we prove the consistency of three asymptotic variance estimators, including that suggested by cameron2011.\footnote{mackinnon2017 also prove the consistency of the estimator of cameron2011 under two-way clustering, but under restrictions that may not hold in practice, as we argue below.} If this latter estimator is asymptotically valid, it has the drawback of being possibly negative in practice. We develop another simple estimator that is also consistent and avoids this drawback. Our Monte Carlo simulations suggest that this estimator may perform significantly better than that suggested by cameron2011 when $\underline{C}$ is small.

Next, we prove the asymptotic validity of a general bootstrap scheme adapted to multiway clustering, called the pigeonhole bootstrap. This resampling scheme differs from the usual multinomial bootstrap by explicitly taking into account the particular dependence structure implied by multiway clustering. The idea is to sample independently each dimensions of clustering and to select cells (with possible repetitions) that are at the intersection of selected clusters in various dimensions. This bootstrap was suggested by mccullagh2000resampling and studied by owen2007pigeonhole but to our knowledge, no weak convergence result has been obtained on it yet, even for sample averages. Again, we prove a general weak convergence result on the pigeonhole bootstrap process. This result implies the validity of the pigeonhole bootstrap for sample averages but also for GMM or smooth functionals of the cdf.\footnote{Similarly to the usual multinomial bootstrap but contrary to, e.g., the wild bootstrap, this bootstrap has the advantage of being universal. Namely, as a resampling scheme, it can be applied in the same way irrespective of the estimation procedure.} Monte Carlo simulations suggest that the pigeonhole bootstrap may work very well even with $\underline{C}$ as small as 5 (resp. 3) under two-way (resp. three-way) clustering.

As in the i.i.d. setting, weak convergence of the empirical process relies on two main ingredients: a multivariate CLT and the asymptotic equicontinuity of the process. To prove the multivariate CLT, we use the Aldous-Hoover representation for exchangeable arrays Aldous1981,Aldous83,Hoover1979,kallenberg05 and techniques related to U-statistics, in particular H\'ajek projections. For the asymptotic equicontinuity, a key step, as in the i.i.d. setting, is to symmetrize the initial process. We do so by generalizing the standard symmetrization lemma vanderVaartWellner1996, using again the Aldous-Hoover representation and an adaptation to our framework of arguments used in ArconesGine1993. The same kind of strategy is used conditional on the data and combined with ergodicity arguments to establish the consistency of the pigeonhole bootstrap process.

The literature on clustering is vast but has mostly focused on linear models under one-way clustering, following the seminal papers of pfeffermann1981, moulton1986 (moulton1986, moulton1990) liang1986 and arellano1987. Without being exhaustive, we also refer to hansen2007, cameron2008bootstrap, carter2017, mackinnon2017wild and hansen2017 for more recent contributions.

The only papers we are aware of considering multiway clustering are the recent works of menzel2017 and mackinnon2017. menzel2017 focuses on sample averages. Contrary to us, he studies inference both with and without asymptotically normality. He also shows that refinements in asymptotic approximations are possible using the wild bootstrap. mackinnon2017 focus on linear regressions with two-way clustering. For such models, they show asymptotic normality and the consistency of the variance estimator of cameron2011. They also show the validity of a certain wild boostrap in this context.

Compared to these papers, our contributions are the following. First, our empirical process result allows us to consider nonlinear estimators. To our knowledge, we are thus the first to show the asymptotic normality of general GMM estimators with multiway clustering. Second, we propose for linear and nonlinear models a new variance estimator that is always positive, very simple to compute and that seems to perform better in practice than that of cameron2011. Third, we show the general validity of the pigeonhole bootstrap with multiway clustering. Finally, even in linear models, we obtain our results under different conditions from those in menzel2017 and mackinnon2017. Contrary to menzel2017, we do not impose cell sizes equal to one, or i.i.d. units within cells. mackinnon2017 assume, through their Assumption 3, that $\overline{N}$, the average of the cell sizes, satisfies $\overline{N} \underline{C}^{2/(2+\lambda)} \rightarrow 0$ for some $\lambda>0$. In other words, the vast majority of cells has to become empty as $\underline{C}$ tends to infinity. This condition may not hold in applications. In contrast, while our framework allows for empty cells, it implies that $\overline{N}$ converges in probability to a positive constant.

The paper is organized as follows. Section (ref) describes the assumptions we impose on the data generating process and the parameters of interest we consider afterwards. Section (ref) provides our main results on the convergence of the empirical process and the pigeonhole bootstrap empirical process. Section (ref) discusses applications of these results to linear and nonlinear estimators. In particular, we show therein the consistency of various asymptotic variance estimators, and asymptotic normality of GMM and smooth functionals of the cdf. We also show the consistency of the pigeonhole bootstrap for inference on such estimators. Section (ref) explores through simulations the finite-sample properties of inference based on asymptotic normality or the pigeonhole bootstrap. Section (ref) concludes. The appendix gathers extensions, additional details on simulations and all the proofs of our results.

The set up

In this section, we define and discuss the restrictions we impose on the data generating process, and the parameters of interest. We suppose to have $k$ non-nested partitions of the population, which correspond to the different dimensions of clustering. We denote the index of the first dimension of clustering (e.g. sector of activity) by $j_1$, the second (e.g. area of residence) by $j_2$ etc. Hereafter, the intersection of $k$ given clusters in the different dimensions (e.g., the second sector of activity and the third area of residence if $(j_1,j_2)=(2,3)$) is called a cell. Cells are indexed by the $k$-tuple $\bm{j}=(j_1,...,j_k)$ for $j_i=1,...,C_i$, where $C_i$ denotes the number of clusters in the sample for dimension $i$. With $k=2$, cells may be seen as matrix entries where the dimensions of clustering would be rows and columns. With $k>2$, cells correspond to the entries of a multidimensional array. We let $\bm{j}\geq \bm{j}'$ to mean that $j_i\geq j'_i$ for all $i=1,...,k$. In the following, we let $\boldsymbol{1}=(1,...,1)$ and $\bm{C}=(C_1,...,C_k)$. The number of observations within each cell is denoted by $N_{\bm{j}}$. The random vector corresponding to unit $\ell=1,...,N_{\bm{j}}$ in cell $\bm{j}$ (with $\boldsymbol{1}\leq \bm{j}\leq \bm{C}$) is then denoted $Y_{\ell,\bm{j}}$, with $Y_{\ell,\bm{j}}\in \mathcal{Y} \subset \mathbb R^l$.

The key assumptions of this $k$-way clustering are the following. First, the sequences $(N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell \geq 1})$ are identically distributed, but not necessarily independent across $\bm{j}$, since two cells with at least one common cluster may face common shocks. Second, $(N_{\bm{j}}, (Y_{\ell,\bm{j}})_{\ell\geq 1})$ and $(N_{\bm{j}'}, (Y_{\ell,\bm{j}'})_{\ell\geq 1})$ are independent if $j_i\neq j'_i$ for all $i=1,...,k$. Third, we consider a sample $(Y_{1,\bm{j}},...,Y_{N_{\bm{j}},\bm{j}})_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}}$ where $\underline{C}=\min_{i\in\{1,...,k\}}C_i$ tends to infinity. Assumption (ref) formalizes all these conditions.

To better understand Assumptions (ref).1 and (ref).2, consider first for simplicity two-way clustering with $N_{\bm{j}}=1$ almost surely. The data can then be depicted as follows.

Assumption (ref).1 imposes for instance that $(Y_{1,(1,1)}, Y_{1,(1,2)})$ has the same distribution as $(Y_{1,(2,1)}$, $Y_{1,(2,2)})$. More generally, data of all rows, or data of all columns, are assumed to have the same distribution. Another way to state this is that the DGP is invariant by a relabelling of each dimension of clustering. This assumption is natural in many settings, with the notable exception of time series. Importantly, Assumption (ref).1 does not impose that $(Y_{1,(1,1)},Y_{1,(1,2)})$ has the same distribution as $(Y_{1,(1,1)}, Y_{1,(2,2)})$. This would indeed amount to neglecting possible dependence within a specific row.

Assumption (ref).2 imposes that any two blocks on the diagonal that do not overlap are independent. In particular $Y_{1,(1,1)}$ and $Y_{1,(2,2)}$ are assumed independent, contrary to, e.g., $Y_{1,(1,1)}$ and $Y_{1,(1,2)}$.\footnote{To be precise, Assumption (ref).2 remains silent on the joint distribution of cells sharing at least one cluster: they may or may not be independent. Thus, i.i.d. sampling of cells is compatible with Assumption (ref).} When combined with Assumption (ref).1, it also implies that cells sharing no rows and columns are mutually independent, since they have the same distribution as cells on the diagonal, which themselves are mutually independent (by applying repeatedly Assumption (ref).2).

Let us come back to the general case with possibly $N_{\bm{j}}\neq 1$. Assumption (ref).2 does not impose any restriction on the distribution of $(N_{\bm{j}}, (Y_{\ell,\bm{j}})_{\ell\geq 1})$. Hence, the dependence between $N_{\bm{j}}$ and the $(Y_{\ell,\bm{j}})_{\ell\geq 1}$, and the dependence between the $(Y_{\ell,\bm{j}})_{\ell\geq 1}$ within cell $\bm{j}$, are left unrestricted. This implies for instance that conditional on $N_{\bm{j}}$, the correlation between $Y_{\ell,\bm{j}}$ and $Y_{\ell',j}$ may vary with $N_{\bm{j}}$. In this sense, we allow for cluster heterogeneity, as defined by carter2017. Also, $Y_{\ell,\bm{j}}$ may have a different distribution from $Y_{\ell',\bm{j}}$, for $\ell\neq \ell'$.

Assumption (ref).3 only excludes arrays that are almost surely empty. Assumption (ref).4 states that only the $N_{\bm{j}}$ first units in each cell $\bm{j}$ are observed. It also specifies our asymptotic framework, in which all dimensions of the array grow large. The condition that $\underline{C}/C_i$ tends to $\lambda_i\geq 0$ is very mild since it allows for different rates of convergence along the different dimensions of clustering.

Whereas the data generating process is defined at the cell level, parameters of interest are virtually always defined at the unit level. To see this, consider again the example of wages. When considering average wages, one usually focuses on units (e.g., individuals) rather than cells, which means that the parameter of interest satisfies

This definition differs from $\mathbb E\left(Y_{\ell,\boldsymbol{1}}\right)$ or its symmetrized version $$\theta_{0,w}=\mathbb E\left(\frac{1}{N_{\boldsymbol{1}}} \sum_{\ell=1}^{N_{\boldsymbol{1}}}Y_{\ell,\boldsymbol{1}}\right),$$ which may seem more natural.\footnote{Note though that the three parameters coincide when $N_{\boldsymbol{1}}=1$ or when $N_{\boldsymbol{1}}$ is independent of $(Y_{\ell,\boldsymbol{1}})_{\ell \geq \boldsymbol{1}}$} But $\theta_0$ is actually the right parameter of interest if one wants to weight equally each individual, rather than weighting equally each cell (e.g., each sector $\times$ area of residence). In the latter case, we would put more weight on individuals lying in small cells. The plug-in estimator of $\theta_0$ corresponds to the empirical mean at the individual level:

where $\Pi_C=\prod_{i=1}^k C_i$ denotes the total number of cells. We study inference on $\theta_0$ based on $\widehat{\theta}$ in Section (ref) below.

More generally, we consider parameters of interest that depend on the unit-level distribution of $Y$, defined by

This implies that the median of wages at the individual level is defined by $\theta_0=F_Y^{-1}(1/2)$. We consider smooth functionals of $F_Y$ in Section (ref) below.

Finally, we consider in Section (ref) moment restrictions at the unit level, rather than at the cell level. Namely, we consider a parameter of interest $\theta_0\in \Theta$ satisfying

for a vector-valued function $m(y,\theta)$. The average parameter defined by (ref) is a particular case of (ref), with $m(Y_{\ell,\boldsymbol{1}},\theta)=Y_{\ell,\boldsymbol{1}}-\theta$. This GMM framework also encompasses linear models and pseudo maximum likelihood estimators of nonlinear models such as logit or probit models. These latter estimators are covered by taking $m$ as the score of the model. Such estimators are not usual maximum likelihood estimators, since they ignore potential correlations between observations within cells, and between cells sharing at least one cluster.

We establish below that in the three cases above, namely expectations, smooth functionals of $F_Y$ and GMM, the corresponding estimators are asymptotically normal. We also develop valid inference on the corresponding estimands. To establish such results, we first study in Section (ref) the asymptotic behavior of empirical processes and their bootstrap counterparts.

Weak convergence results

Empirical processes

Let $\mathcal{F}$ denote a class of real-valued functions. In this section, we study the empirical process $\mathbb{G}_{C}$ defined on $\mathcal{F}$ by $$\mathbb{G}_{C}f=\sqrt{\underline{C}}\left\{\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}\sum_{\ell=1}^{N_{\bm{j}}} f(Y_{\ell,\bm{j}}) - \mathbb E\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}f(Y_{1,\boldsymbol{1}})\right]\right\}.$$ Specifically, we prove that under restrictions on $\mathcal{F}$, $\mathbb{G}_{C}$ converges weakly to a Gaussian process as $\underline{C}$ tends to infinity. While we refer to, e.g., vanderVaartWellner1996 for a formal definition of weak convergence of empirical processes, we recall that this result is stronger than pointwise asymptotic normality of $\mathbb{G}_{C}f$. Our result below will therefore entail central limit theorems for means of the form $$\frac{1}{\Pi_C} \sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}\sum_{\ell=1}^{N_{\bm{j}}} f(Y_{\ell,\bm{j}}),$$ and therefore, by the delta method (considering $f(y)=y$ and $f(y)=1$), for sample averages defined by (ref). But such a result is not sufficient for the asymptotic normality of, e.g., smooth functionals of the empirical cdf or GMM estimators. Convergence of the whole process, on the other hand, allows one to establish such results. To establish this convergence, we cannot apply standard results on the empirical process for two reasons. First, the different cells are potentially dependent rather than i.i.d. Second, even if they were i.i.d., we do not consider the usual empirical process at the cell-level, because the class of functions is defined at the unit level and we sum over a random number of units within each cell.

Before giving our main asymptotic result on $\mathbb{G}_{C}$, we introduce additional notation related to a generic class $\mathcal{G}$. An envelope of $\mathcal{G}$ is a measurable function $G$ satisfying $G(u)\geq\sup_{f\in \mathcal{G}}|f(u)|$. For any $\varepsilon>0$ and any norm $||.||$ on a space containing $\mathcal{G}$, $N(\varepsilon,\mathcal{G},||.||)$ denotes the minimal number of $||.||$-closed balls of radius $\varepsilon$ with centers in $\mathcal{G}$ needed to cover $\mathcal{G}$.\footnote{With a slight abuse of language, we use here the term norms in lieu of seminorms. For instance, Assumption (ref) involves seminorms rather than norms. Also, we use $\|.\|$ for (semi)norms on functions and $|.|$ for norms on finite-dimensional objects. Specifically, for any vector $b$, $|b|$ denotes the Euclidean norm of $b$; and for any matrix $A$, $|A|$ denotes the Frobenius norm of $A$.} The norms we consider hereafter are $\|f\|_{\mu,r}=(\int |f|^rd\mu)^{1/r}$ for any $r\geq 1$ and probability measure $\mu$. Finally, a class of measurable functions $\mathcal{G}$ is pointwise measurable if there exists a countable subclass $\mathcal{H}\subset \mathcal{G}$ such that elements of $\mathcal{G}$ are pointwise limit of elements of $\mathcal{H}$.

We consider the following standard assumptions on the class $\mathcal{F}$ indexing $\mathbb{G}_C$.

Assumption (ref) is not necessary but usually imposed cherno2014,Kato2016 to avoid measurability issues and the use of outer expectations. For further discussion about these classes, we refer to Kosorok2006 (Kosorok2006, pp.137-140). Assumption (ref) imposes a condition on what is usually referred to as the uniform entropy integral, see, e.g., vanderVaartWellner1996. Finiteness of the uniform entropy integral is satisfied by any VC-type class of functions cherno2014, or by the convex hull of such classes under some restrictions. These conditions are nearly the same as those used with i.i.d. data. The only difference lies in the moment conditions. When $\mathcal{F}$ is finite, we require a second moment condition that is the exact analog of the moment condition for usual central limit theorems. When $\mathcal{F}$ is infinite, on the other hand, we require the slightly stronger condition $\mathbb E\left[N_{\boldsymbol{1}}^2\right]<+\infty$ and $\mathbb E\left[\left(N_1 \sum_{\ell=1}^{N_{\boldsymbol{1}}} F^2\left(Y_{\ell,\boldsymbol{1}}\right)\right)\right]<+\infty$. Note however that the two conditions are equivalent whenever $N_{\boldsymbol{1}}$ is bounded.

Theorem (ref) shows the weak convergence of $\mathbb{G}_C$ towards a centered gaussian process $\mathbb{G}$, and gives the form of the covariance kernel of $\mathbb{G}$. The result holds under Assumption (ref), but this is not the only possible restriction on the class of functions. In Appendix (ref), we show the same result under smoothness restrictions on $\mathcal{F}$ instead of Assumption (ref).

Let us summarize the proof of Theorem (ref). Weak convergence of $\mathbb G_C$ holds under two main conditions. First, $(\mathbb G_Cf_1,...,\mathbb G_C f_m)$ should be asymptotically normal for any $(f_1,...,f_m)$ in $\mathcal{F}$ and any $m\geq 1$. Second, one should establish asymptotic equicontinuity. Regarding finite-dimensional convergence, we proceed in several steps. To simplify the discussion, we consider here the case where $m=1$ and two-way clustering. We first exploit the Aldous-Hoover representation Aldous1981,Hoover1979,kallenberg05, which extends de Finetti's theorem to separately exchangeable random sequences. This result ensures the existence of mutually independent variables $(U_{j_1, 0}, U_{0,j_2}, U_{\bm{j}})_{j_1\geq 1, j_2\geq 1, \bm{j}\geq \boldsymbol{1}}$ such that for all $\bm{j}$,

The variable $U_{j_1,0}$ (resp. $U_{0, j_i}$) may be seen as a shock specific to cluster 1 (resp. 2), while $U_{\bm{j}}$ can be interpreted as a shock specific to Cell $\bm{j}$.\footnote{With more than two dimensions of clustering, the representation is similar but we have to include shocks specific to each subset of $(j_1,...,j_k)$. For instance, with $k=3$, we have also to consider shocks such as $U_{j_1,j_2,0}$.}

In the second step, we consider the H\'ajek projection of $\mathbb{G}_C f_1$ on the set $\mathscr{S}$ of random variables depending only on the marginal cluster specific factors, namely $$\mathscr{S}=\left\{\sum_{j_1=1}^{C_1} g_{j_1,0}(U_{j_1,0}) + \sum_{j_2=1}^{C_2} g_{0,j_2}(U_{0,j_2}), \; g_{j_1,0}\in L^2(U_{j_1,0}), g_{0,j_2}\in L^2(U_{0,j_2})\right\}.$$ We prove that $\mathbb{G}_Cf_1$ gets close, in a $L^2$ sense, to its H\'ajek projection as $\underline{C}\rightarrow\infty$. Asymptotic normality then follows by a simple CLT on the H\'ajek projection.

To complete the proof of the theorem, we have to establish asymptotic equicontinuity. Roughly speaking, this means that whenever $f_1$ and $f_2$ are close to each other, $\mathbb G_Cf_1-\mathbb G_Cf_2$ is close to zero vanderVaartWellner1996. For that purpose, we prove a symmetrization lemma similar to Lemma 2.3.1 in vanderVaartWellner1996. To do so, we adapt arguments used in the proofs of Theorem 3.1 in ArconesGine1993 where independent copies of random variables are introduced to control U-statistics. Following this idea, we introduce independent copies of the $(U_{\bm{j}})_{\bm{j}>0}$ that come from the Aldous-Hoover representation. By the symmetrization lemma, we can then bound fluctuations of $\mathbb G_C$ by a function of the entropy of the class $$\widetilde{\mathcal{F}}=\left\{g(n,y_1,...,y_n)=\sum_{i=1}^{n}f(y_i): n\in \mathbb{N}, (y_1,...,y_n) \in \mathcal{Y}^{n}; f\in \mathcal{F} \right\}.$$ Note that this class is related to, but different from $\mathcal{F}$. We have defined the class of function $\mathcal{F}$ at the unit (e.g., individual) level because parameters of interest are defined at this level. But the stochastic model in Assumption (ref) is stated at the cell level, which explains why, intuively, we need to control the complexity of $\widetilde{\mathcal{F}}$. We show that this is possible under Assumption (ref). By what precedes, this implies the asymptotic equicontinuity of $\mathbb G_C$.

We now comment on the asymptotic kernel $K$ of $\mathbb G_C$. For simplicity, let $f_1(y)=f_2(y)=y$ and define $S_{\bm{j}}=\sum_{\ell=1}^{N_{\boldsymbol{1}}} Y_{\ell,\bm{j}}$. Theorem (ref) implies that the asymptotic variance of $\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} S_{\bm{j}}/\Pi_C$ is

This formula may seem surprising, because it is not obvious at first glance that it is positive. But it turns out that under Assumption (ref), each covariance term is positive. Specifically, and considering for simplicity $k=2$, we establish in the proof of Theorem (ref) that $$\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\bm{2}_1}\right) = \mathbb V\left(\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}f(Y_{\ell,\boldsymbol{1}})\bigg| U_{\boldsymbol{1}, 0}\right)\right),\; \mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\bm{2}_2}\right) = \mathbb V\left(\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}f(Y_{\ell,\boldsymbol{1}})\bigg| U_{0, \boldsymbol{1}}\right)\right),$$ where $U_{\boldsymbol{1}, 0}$ and $U_{0,\boldsymbol{1}}$ appear in the representation (ref).

Now, let us give some intuitions on (ref). This formula involves cells sharing exactly one common cluster, namely cluster 1 in dimension $i$. To better understand why only such terms appear, consider $$\mathbb V\left(\frac{\sqrt{\underline{C}}}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} S_{\bm{j}}\right).$$ This variance is complicated because of the particular dependence structure due to multiway clustering. To simplify it, we can write it as the sum of covariances between cells sharing no common cluster, cells sharing one common cluster... and finally the covariances of cells with themselves. The number of pairs of cells sharing no common cluster is $\Pi_C \times \prod_{i=1}^k (C_i-1)$, which is of the order $\Pi^2_C$ as $\underline{C}$ tends to infinity. The number of pairs of cells sharing one common cluster is $\Pi_C \sum_{i=1}^k \prod_{j\neq i} (C_j-1)$, which is smaller than $k\Pi^2_C/\underline{C}$. Hence, the number of such pairs of cells is negligible compared to the number of pairs of cells sharing no common cluster. Similarly, we can prove that the number of cells sharing more than one common cluster is negligible compared with the number of cells sharing one common cluster. Hence, intuitively, the variance will be equivalent to the sum of only covariances between cells sharing either no or just one common cluster. But by independence, the covariance between cells sharing no common cluster is actually zero. Hence, at the end of the day, we only get covariances between cells sharing just one common cluster.

Pigeonhole bootstrap processes

We now consider the bootstrap counterpart of the weak convergence result in Theorem (ref). Bootstrap offers several advantages over usual inference based on asymptotic normality. First, it avoids the computation of theoretical formulas of asymptotic variances, which can be difficult with, e.g., multistep estimators. Second, it often exhibits a better behavior than normal approximations in finite samples. Still, a consistent bootstrap scheme in our clustering setting needs to reproduce the dependence between cells. We consider for that purpose the “pigeonhole bootstrap”, suggested by mccullagh2000resampling and studied, in the case of the sample mean and for particular models, by owen2007pigeonhole. We are, however, not aware of any result concerning the asymptotic validity of the pigeonhole bootstrap for inference. Theorem (ref) below aims to fill this gap.

We first recall the principle of the pigeonhole bootstrap:

By construction, any bootstrap sample consists of exactly $\Pi_C$ cells. Also, dependence between cells sharing cluster $i$ is achieved through the term $W^i_{j_i}$. Actually, one can check that conditional on the data $(N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\geq 1})_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}$, the bootstrap weights $(W_{\bm{j}})_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}}$ satisfy the first condition in Assumption (ref), and the second asymptotically.\footnote{As in the i.i.d. setting where bootstrap weights are asymptotically independent, the weights of cells sharing no common cluster become independent as $\underline{C}\rightarrow\infty$.} This suggests that the pigeonhole bootstrap could be asymptotically valid.

We now consider the bootstrap counterpart of the empirical process $\mathbb G_C$. For any $f\in\mathcal{F}$, let us define $$\mathbb{G}^{\ast}_C(f)=\frac{\sqrt{\underline{C}}}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}} \left(W^{\bm{C}}_{\bm{j}}-1\right) \sum_{\ell=1}^{N_{\bm{j}}}f(Y_{\ell,\bm{j}}).$$

The asymptotic validity of the pigeonhole bootstrap amounts to showing that conditional on the data $\{N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\geq 1}\}_{\bm{j}\geq \boldsymbol{1}}$, $\mathbb{G}_C^{\ast}$ converges weakly in probability to the process $\mathbb{G}$ defined in Theorem (ref). As discussed in, e.g., vanderVaartWellner1996 (1996, Chapter 3.6), conditional weak convergence in probability amounts to proving

where $\text{BL}_1$ is the set of bounded and Lipschitz functions from $\ell^{\infty}(\mathcal{F})$ to $\mathbb{R}$.

As we shall see below, this theorem ensures the asymptotic validity of the pigeonhole bootstrap not only for sample means, but also for smooth functionals of the empirical cdf and GMM estimators. The proof of Theorem (ref) follows the same lines as that of Theorem (ref) to get weak convergence of the bootstrap process conditionally on the original data. To ensure the unconditional boostrap consistency we also use some ergodicity arguments kallenberg05 and we prove Lindeberg-Feller conditions for some statistics defined on the exchangeable array.

Applications

Simple averages and linear models

As before, let $S_{\bm{j}}=\sum_{\ell=1}^{N_{\bm{j}}} Y_{\ell,\bm{j}}$. We first investigate here how inference can be conducted on $\theta_0=E(S_{\boldsymbol{1}})$ based on the plug-in estimator

We focus first on $\widehat{\theta}$ for simplicity but show at the end of the section how our reasoning extends to sample averages as defined by (ref) and linear models. Provided that $E(S_{\bm{j}}^2)<+\infty$, we have, by Theorem (ref),

A first strategy to make inference on $\theta_0$ is therefore to use the normal approximation and a consistent estimator of the asymptotic variance. A second strategy is to rely on the pigeonhole bootstrap.\footnote{A third strategy is to rely on other bootstrap schemes. We refer to menzel2017 for the construction and analysis of a wild bootstrap procedure for sample averages on such clustered data.}

First, let us consider inference based on asymptotic normality. The asymptotic variance $V$ depends on $\lambda_i$ and $\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right)=\mathbb E\left((S_{\boldsymbol{1}}-\theta_0)(S_{\boldsymbol{2}_i}-\theta_0)\right)$, for $i=1,...,k$. $\lambda_i$ can simply be approximated by $\underline{C}/C_i$. Regarding the covariance term, observe that $\boldsymbol{1}$ and $\boldsymbol{2}_i$ share exactly one cluster. It is then natural to consider the estimator $$\widehat\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right) = \frac{1}{C_i\prod_{s\neq i}C_s(C_s-1)}\sum_{(\bm{j},\bm{j}')\in\mathcal{A}_i}\left(S_{\bm{j}}-\widehat{\theta}\right)\left(S_{\bm{j}'}-\widehat{\theta}\right)',$$ where $\mathcal{A}_i:=\left\{(\bm{j},\bm{j}'):j_i=j_i',\quad j_s\neq j_s'\quad\forall s\neq i\right\}$. This estimator is the average of cross products between clusters sharing just one common cluster, the denominator $C_i\prod_{s\neq i}C_s(C_s-1)$ corresponding to the number of such pairs. This leads to the following estimator for $V$:

We will show that $\widehat{V}_2$ is consistent for $V$. A major drawback of this estimator, however, is that it is not necessarily positive. Also, $V$ is the variance of a H\'ajek projection, as explained above. As such, it is likely to underestimate $\mathbb V(\sqrt{\underline{C}}\widehat{\theta})$. Because $\widehat\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right)$ itself slightly underestimates $\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right)$ ($\mathbb E[\widehat\mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right)] = \mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right) - \mathbb V(\widehat{\theta})$), we can expect the corresponding confidence regions to undercover in practice. This intuition is confirmed in our simulations below.

To avoid these issues, we suggest to simply add to $\widehat{V}_2$ pairs sharing more than one cluster. Specifically, we consider

From an asymptotic point of view, the additional terms in $\widehat{V}_1$ correspond to pairs sharing more than one cluster. We show in the proof of Proposition (ref) below that such terms are negligible, implying that $\widehat{V}_1$ is consistent, just as $\widehat{V}_2$. In finite samples, on the other hand, $\widehat{V}_1$ has several advantages over $\widehat{V}_2$. First, $\widehat{V}_1$ is positive, as (ref) shows. Also, it is likely to overestimate $V$, but this may somewhat compensate the fact that $V$ itself underestimates $\mathbb V(\sqrt{\underline{C}}\widehat{\theta})$. And indeed, in the simulations considered in Section (ref) below, inference is more accurate when using $\widehat{V}_1$ rather than $\widehat{V}_2$, in particular when $\underline{C}$ is small. A last advantage is computational. Equation (ref) shows that we can compute this estimator using variance estimators of $\widehat{\theta}$ assuming only one-way clustering along dimensions $i\in\{1,...,k\}$, and then summing these different variances. $\widehat{V}_2$, on the other hand, cannot be obtained as easily. For all these reasons, we recommend using $\widehat{V}_1$ rather than $\widehat{V}_2$ in practice.

We now compare our two estimators with that proposed by cameron2011. Their estimator relies on a reformulation of $\mathbb V(\widehat{\theta})$. For any $m\in\{1,...,k\}$ and $1\leq i_1<...<i_m\leq k$, let $\mathcal{B}_{i_1,...,i_m}=\left\{(\bm{j},\bm{j}'):j_{i_1}=j'_{i_1},...,j_{i_m}=j'_{i_m}\right\}$. Then

The first line follows because if $(\bm{j},\bm{j}')$ share no cluster, $\mathbb{C}ov\left(S_{\bm{j}}, S_{\bm{j}'}\right)=0$. The second line follows from the inclusion-exclusion principle. This leads to the following estimator for the asymptotic variance of $\widehat{\theta}$

We can consider various finite sample adjustments where $1/(\Pi_C)^2$ is replaced by $c_{i_1,...,i_m}/(\Pi_C)^2$, with $c_{i_1,...,i_m}$ tending to one as $\underline{C}$ tends to infinity. We refer to cameron2011 for more details. As with $\widehat{V}_1$, the appeal of Formula (ref) is that we can compute this estimator using variance estimators of $\widehat{\theta}$ assuming only one-way clustering along dimensions $(i_1,...,i_m)$, for all $1\leq i_1 < ... < i_m \leq k$. The estimator $\widehat{V}_{\text{cgm}}$ is still slightly more complicated to compute than $\widehat{V}_1$, as the latter only requires the computation of one-way clustering variances along dimensions $i=1,...,k$.

To further understand the links and differences between $\widehat{V}_1$ and $\widehat{V}_{\text{cgm}}$, it is instructive to consider the case $k=2$. Then the formulas simplify to

In other words, $\widehat{V}_1$ estimates $V$ by counting twice the pairs $(\bm{j},\bm{j}')$ sharing two clusters, or equivalently the pairs $(\bm{j},\bm{j})$, $\boldsymbol{1} \leq \bm{j}\leq \bm{C}$. $\widehat{V}_{\text{cgm}}$ counts such pairs only once, whence the correction in (ref). The cost of this correction is that $\widehat{V}_{\text{cgm}}$ is not always positive. Finally, note that there are only $\Pi_C$ pairs $(\bm{j},\bm{j})$. Thus, the second term in (ref) is of order $\underline{C}/\Pi_C$ and tends to 0 as $\underline{C}$ tends to infinity. We can therefore expect $\widehat{V}_1$ and $\widehat{V}_{\text{cgm}}$ to be asymptotically equivalent.

Finally, for any $k\in \{1,2,\text{cgm}\}$ and $\alpha\in (0,1)$, we consider confidence regions $R_k^{1-\alpha}$ for $\theta_0$ defined by $$R_k^{1-\alpha}= \left\{\theta: \, \underline{C}(\theta-\widehat{\theta})' \widehat{V}_k^{-1} (\theta-\widehat{\theta}) \leq \chi^2_L(1-\alpha)\right\},$$ where $\chi^2_L(1-\alpha)$ is the quantile of order $1-\alpha$ of a $\chi^2_L$ distribution. Proposition (ref) shows that $\widehat{V}_1$, $\widehat{V}_2$ and $\widehat{V}_{\text{cgm}}$ are all consistent estimators of $V$, implying that the confidence regions are also asymptotically valid, as long as $V$ is positive definite.

The condition that $V$ is positive definite basically states that at least one of the dimension of clustering matters, in the sense that for at least one $i\in\{1,...,k\}$, $\mathbb{C}ov(S_{\boldsymbol{1}},S_{\boldsymbol{2}_i})$ is positive definite. Note that under Assumption (ref), $\mathbb{C}ov(S_{\boldsymbol{1}},S_{\boldsymbol{2}_i})$ is necessarily positive; but it may not be positive definite. For instance, consider two-way clustering and $S_{\bm{j}}=U_{j_1,0}+U_{0,j_2}+U_{\bm{j}}\in \mathbb R$, where the $(U_{j_1,0})_{j_1}$, $(U_{0,j_2})_{j_2}$ and $(U_{\bm{j}})_{\bm{j}}$ are all mutually independent. Then $\mathbb{C}ov(S_{\boldsymbol{1}},S_{\boldsymbol{2}_i})>0$ if and only if $\mathbb V(U_{j_1,0})+\mathbb V(U_{0,j_2})>0$. As discussed in menzel2017, $\mathbb{C}ov(S_{\boldsymbol{1}},S_{\boldsymbol{2}_i})$ may not be positive definite even when $S_{\boldsymbol{1}}$ and $S_{\boldsymbol{2}_i}$ are dependent. This is the case if we modify the example above by assuming instead

If $V$ is not positive definite, standard tests and confidence regions are not valid in general. When $V=0$, $\widehat{\theta}$ actually converges at a rate faster than $1/\sqrt{\underline{C}}$ and its asymptotic distribution may be non-normal. This is the case for instance if (ref) holds. We refer to menzel2017, Example 1.6, for more details.

We now turn to the pigeonhole bootstrap. Let $$\widehat{\theta}^* = \frac{1}{\Pi_C} \sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} W_{\bm{j}} S_{\bm{j}},$$ where $W_{\bm{j}}$ is defined in Section (ref) above. Then let $q_{1-\alpha}^*$ denote the quantile of order $1-\alpha$ of the distribution of $|\widehat{\theta}^*-\widehat{\theta}|$ conditional on the data. We consider the confidence region $R_{\text{boot}}^{1-\alpha}$ for $\theta_0$ defined by $$R_{\text{boot}}^{1-\alpha}= \left\{\theta: \, |\widehat{\theta}-\theta|\leq q_{1-\alpha}^*\right\}.$$ The asymptotic validity of $R_{\text{boot}}^{1-\alpha}$ is an immediate consequence of Theorem (ref).

When $\theta_0\in \mathbb R$, an alternative, popular confidence region is the percentile bootstrap. This amounts to considering $[q_{\alpha/2}(\widehat{\theta}^*), q_{1-\alpha/2}(\widehat{\theta}^*)]$. This interval is also valid asymptotically, since the asymptotic distribution of $\widehat{\theta}-\theta_0$ is normal, and therefore symmetric.

We now discuss how Propositions (ref) and (ref) extend to other parameters of interest. First, let us consider $\theta_0=E(S_{\boldsymbol{1}})/E(N_{\boldsymbol{1}})$, as in Section (ref). Assume that $\mathbb E(S_{\boldsymbol{1}}^2)<\infty$ and $\mathbb E(N_{\boldsymbol{1}}^2)<\infty$. By Theorem (ref) applied to $\mathcal{F}=\{\text{Id},1\}$ and the delta method, we have

where $T_{\bm{j}} = (S_{\bm{j}} - N_{\bm{j}} \theta_0)/\mathbb E(N_{\boldsymbol{1}})$. We can then estimate the asymptotic variance of the sample average as previously, by simply replacing $S_{\bm{j}}-\widehat{\theta}$ in (ref), (ref) and (ref) by $$\widehat{T}_{\bm{j}} = \frac{S_{\bm{j}} - N_{\bm{j}} \widehat{\theta}}{\frac{1}{\Pi_C} \sum_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}} N_{\bm{j}}}.$$ Consistency follows as in the proof of Proposition (ref), using consistency of $\widehat{\theta}$ and $\sum_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}} N_{\bm{j}}/\Pi_C$. The pigeonhole bootstrap is also valid for $\mathbb E(S_{\boldsymbol{1}})/\mathbb E(N_{\boldsymbol{1}})$ by applying the simple delta method for the bootstrap, see e.g. Theorem 23.5 in vanderVaart2000. More generally, Propositions (ref) and (ref) extend to parameters of the form $g(\theta_{01},...,\theta_{0R})$, where $\theta_{0r}=\mathbb E(\sum_{\ell=1}^{N_{\boldsymbol{1}}} q_r(N_{\boldsymbol{1}},Y_{\ell,\boldsymbol{1}}))$ ($r=1,...,R$), provided that $g$ is continuously differentiable at $(\theta_{01},...,\theta_{0R})$.

Finally, let us consider linear models. Then $Y_{\ell,\bm{j}}=(\tilde{Y}_{\ell,\bm{j}},X'_{\ell,\bm{j}})'$, with $\tilde{Y}_{\ell,\bm{j}}$ the outcome variable and $X_{\ell,\bm{j}}$ a vector of covariates. Then the parameter of interest $\theta_0$ and its estimator satisfy

We first show that $\widehat{\theta}$ is asymptotically normal and characterize its asymptotic variance. Hereafter, we define $u_{\ell,\bm{j}}=\widetilde{Y}_{\ell,\bm{j}}-X'_{\ell,\bm{j}}\theta_0$.

Next, we show similar results as in Propositions (ref) and (ref). For conciseness, we focus on an estimator of $V$ similar to $\widehat{V}_1$ rather than $\widehat{V}_2$ and $\widehat{V}_{\text{cgm}}$. Let $\widehat{J}=\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}\sum_{\ell=1}^{N_{\bm{j}}}X_{\ell,\bm{j}} X'_{\ell,\bm{j}}$, $\widehat{u}_{\ell,\bm{j}}=\widetilde{Y}_{\ell,\bm{j}}-X'_{\ell,\bm{j}}\widehat{\theta}$ and $$\widehat{H}=\sum_{i=1}^k\frac{\underline{C}}{C_i} \frac{1}{C_i}\sum_{j'_i=1}^{C_i} \left(\frac{1}{\prod_{s\neq i}C_s}\sum_{\bm{j}:j_i=j'_i} \sum_{\ell=1}^{N_{\bm{j}}}X_{\ell,\bm{j}}\widehat{u}_{\ell,\bm{j}}\right) \left(\frac{1}{\prod_{s\neq i}C_s}\sum_{\bm{j}:j_i=j'_i} \sum_{\ell=1}^{N_{\bm{j}}} \widehat{u}_{\ell,\bm{j}} X'_{\ell,\bm{j}}\right)$$ Then let $\widehat{V}=\widehat{J}^{-1}\widehat{H}\widehat{J}^{-1}$.

Propositions (ref) and (ref) complement the results of mackinnon2017 by showing asymptotic normality and the validity of two inference methods without assuming that the average of the cell sizes $\overline{N}$ satisfies $\overline{N} \underline{C}^{2/(2+\lambda)} \rightarrow 0$ for some $\lambda>0$. Proposition (ref) also shows the consistency of a new, positive, variance estimator and the asymptotic validity of the pigeonhole bootstrap in this context of linear models.

Nonlinear functionals of the distribution

Simple central limit theorems and the usual delta method are not sufficient to yield the asymptotic normality of estimators such as the sample median. We now show how the results in Section (ref) can be applied to such smooth, nonlinear functionals of the empirical distribution. Let $F_Y$ be defined as in (ref) and let $\theta_0=g(F_Y)$. To take examples related to income inequalities (so that here the support of $Y$ is $\mathbb R^+$), we may consider for instance quantiles, interquantile ratios and poverty rates, for which we have respectively $g(F_Y)=F_Y^{-1}(\tau)$ for any $\tau\in (0,1)$, $g(F_Y)=F_Y^{-1}(\tau)/F_Y^{-1}(1-\tau)$ for $q\in (1/2,1)$ and $g(F_Y)=F_Y(\alpha F_Y^{-1}(\beta))$ for $(\alpha,\beta)\in (0,1)^2$. Other examples include the Kaplan-Meier functional vanderVaart2000 or the nonlinear difference-in-difference estimand of Athey06, for which $\theta_0=\int y dF_1(y) - \int \left[F_2^{-1}\circ F_3(y)\right] dF_4(y)$, where $(F_1,...,F_4)$ are the cdf's of $Y$ on four distinct subpopulations..

We consider the plug-in estimator $\widehat{\theta}=g(\widehat{F_Y})$ of $\theta_0$, with

To state the smoothness condition on $g$, we need additional notation and definitions. Let $\mathbb{D}$ denote a subset of the set of all cumulative distribution functions on $\mathbb R^\ell$ and suppose that $g:\mathbb{D}\mapsto \mathbb R^r$. We consider for simplicity here vector-valued functions $g$, but could easily extend our result below to functions taking values in normed spaces. We say that $g$ is Hadamard differentiable at $F_Y$ tangentially to $\mathbb{D}_0$ if there exists a continuous, linear map $g'_{F_Y}:\mathbb{D}\mapsto \mathbb R^r$ such that for every $(h_t)_{t\in \mathbb R^+}$ such that $h_t\rightarrow h\in\mathbb{D}_0$ as $t\downarrow 0$, $$\lim_{t\downarrow 0} \left|\frac{g(F_Y+th_t) - g(F_Y)}{t} - g'_{F_Y}(h)\right|=0.$$

Proposition (ref) shows that if $g$ is Hadamard differentiable at $F_Y$, $g(\widehat{F_Y})$ will be asymptotically normal. We also consider confidence regions based on the bootstrap. As before, we let $R_{1-\alpha}^{\text{boot}}= \left\{\theta: \, |\theta-\widehat{\theta}|\leq q_{1-\alpha}^*\right\}$, where $q_{1-\alpha}^*$ denotes the quantile of order $1-\alpha$ of the distribution of $|\widehat{\theta}^*-\widehat{\theta}|$ conditional on the data.

The first part follows from Theorem (ref), and a linearization of the ratio akin to (ref). The second follows from the first part and the functional delta method vanderVaart2000. The third part is a direct consequence of Theorem (ref) and the functional delta method for the bootstrap vanderVaartWellner1996.

As an illustration, let us consider the example of a quantile, $\theta_0=F_Y^{-1}(\tau)$ for some $\tau\in (0,1)$. Suppose that $F_Y$ is differentiable at $\theta_0$. Then the function $g(F_Y)=F_Y^{-1}(\tau)$ is Hadamard differentiable at $F_Y$, tangentially to the set of functions that are continuous at $\theta_0$ vanderVaart2000. Moreover, we prove in Appendix (ref) that if $\mathbb E[N_{\boldsymbol{1}}^{2+\zeta}]<+\infty$ for some $\zeta>0$, $\mathbb{G}$ is almost surely continuous at $\theta_0$. Hence, Proposition (ref) ensures that $\widehat{\theta}$ is asymptotically normal. Moreover, by Point 3 of the proposition, inference based on the bootstrap is valid, as long as its asymptotic variance is strictly positive.

The third part ensures the consistency of the pigeonhole bootstrap. in principle, one could also use the normal approximation and a consistent estimator of $\mathbb V(g'_{F_Y}(\mathbb{G}))$ to make inference on $\theta_0$. $g'_{F_Y}(\mathbb{G})$ is a linear functional of $\mathbb{G}$, so the same ideas as in Section (ref) above can be applied. This linear functional may however depend on complicated functions of $F_Y$ that must be estimated. For instance, when $g(F)=F^{-1}(\tau)$, $g'_{F_Y}(\mathbb{G})$ depends on the derivative of $F_Y$ taken at $F^{-1}(\tau)$. As a result, additional restrictions may be necessary to achieve the consistency of the variance estimator. We do not explore this avenue further here, as it depends very much on the functional $g$, but consider explicitly this approach in the following section on GMM.

GMM estimators

Finally, we consider parameters defined by the moment restrictions (ref), with possibly nonsmooth moments. We suppose that $m(y,\theta)\in \mathbb R^L$, with $m(y,\theta)=(m_1(y,\theta),...,m_L(y,\theta))'$. We show the asymptotic normality of $\widehat{\theta}$ under the following condition.

Assumptions (ref).1, (ref).2 and (ref).6 are standard. Assumption (ref).3, combined with (ref).1 and (ref).2, ensures that the minimum of $$\theta\mapsto \mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right)'\Xi \, \mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right)$$ is well-separated on $\Theta$, thus ruling out possible inconsistency of the GMM estimator. Note that Assumption (ref).3 is weaker than the standard continuity assumption of $\theta \mapsto m(Y_{\ell,\boldsymbol{1}},\theta)$, which may fail if $m$ includes for instance indicator functions. Assumption (ref).4 is standard in GMM with nonsmooth moments where $\theta \mapsto m(Y_{\ell,\boldsymbol{1}},\theta)$ is not differentiable, see e.g. Condition (ii) in Theorem 7.2 of newey1994large. Finally, by Theorem (ref), Assumption (ref).5 ensures the stochastic equicontinuity condition newey1994large, which together with (ref).4, is key to obtain $\sqrt{\underline{C}}$-asymptotic normality of $\widehat{\theta}$.

To illustrate that Assumption (ref) can handle nonsmooth moments, let us consider the example of quantile IV regressions. Let $Y_{\ell, \bm{j}}=(W_{\ell, \bm{j}}, X'_{\ell, \bm{j}}, Z'_{\ell, \bm{j}})'$, where $W_{\ell, \bm{j}}\in \mathbb R$ denotes the outcome variable, $X_{\ell, \bm{j}}\in \mathbb R^p$ denotes the explanatory, potentially endogenous variable and $Z_{\ell, \bm{j}}\in \mathbb R^L$ denotes the set of instruments ($Z_{\ell, \bm{j}}$ may include some components of $X_{\ell, \bm{j}}$). The moment functions are then $$m(Y_{\ell,\boldsymbol{1}},\theta)=Z_{\ell, \bm{j}}\left(\tau - \mathds{1}\{W_{\ell, \bm{j}}- X'_{\ell, \bm{j}}\theta\leq 0\}\right).$$ Let us assume for simplicity that the $(Y_{\ell, \boldsymbol{1}})_{\ell\geq 1}$ are identically distributed. We show in Appendix (ref) that Assumptions (ref).3-(ref).5 hold if, basically, $\mathbb E[N_{\boldsymbol{1}}^2 |Z_{1,\boldsymbol{1}}|^2]<+\infty$, $X$ is in a compact set, the conditional cdf $F_{W_{1,\boldsymbol{1}}|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}}}(\cdot|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}})$ is continuous everywhere and admits a bounded derivative $f_{W_{1,\boldsymbol{1}}|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}}}(\cdot|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}})$ in a neighborhood of $X_{1,\boldsymbol{1}}'\theta_0$ and the rank of $\mathbb E\left[N_{\boldsymbol{1}} X_{1,\boldsymbol{1}} Z'_{1,\boldsymbol{1}}f_{W_{1,\boldsymbol{1}}|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}}}(X'_{\ell, \bm{j}}\theta_0|X_{1,\boldsymbol{1}},Z_{1,\boldsymbol{1}})\right]$ is equal to $p$.\footnote{For the exact conditions, see Assumption (ref) in Appendix (ref).}

To our knowledge, Theorem (ref) is the first result on the asymptotic normality of GMM estimators under multiway clustering. Theorem (ref) gives also the expression of the asymptotic variance $V_0$. This matrix takes the usual form, except that the matrix $H$, which would simply be $E[m(Y_{\ell,\boldsymbol{1}},\theta_0)m(Y_{\ell,\boldsymbol{1}},\theta_0)']$ without clustering, takes a more complicated form here. This form is in line with our result on the covariance kernel of the empirical process considered above.

We now turn to inference on $\theta_0$. As for sample averages, we consider inference based on asymptotic normality and a consistent estimator of $V_0$, or the pigeonhole bootstrap. To ensure the consistency of our estimator of $V_0$, we impose the following additional regularity condition.

Assumption (ref).1 refines Assumption (ref).4 by imposing some structure on the Jacobian matrix of $\theta\mapsto\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right)$. In Assumption (ref).2, the condition on the classes are of Glivenko-Cantelli type, and weaker than Assumption (ref). The continuity conditions in Points 3 and 4 are similar to Assumption (ref).3, but are imposed on different functions related to $J$ and $H$ rather than on the moment conditions themselves.

We now define our estimator of $V_0$, which is based on estimators of $J$ and $H$. Given Assumption (ref).1, $\widehat{J}$ is the simple plug-in estimator $$\widehat{J}=\frac{1}{\Pi_C}\sum_{\boldsymbol{1}\leq \bm{j} \leq \bm{C}} \sum_{\ell=1}^{N_{\bm{j}}}d(Y_{\ell,\bm{j}},\widehat{\theta}).$$ To estimate $H$, we adapt our previous estimator $\widehat{V}_1$ to this context by considering

Our variance estimator is then $\widehat{V}=(\widehat{J}'\widehat{\Xi}\widehat{J})^{-1}\widehat{J}'\widehat{\Xi} \widehat{H}\widehat{\Xi} \widehat{J}(\widehat{J}'\widehat{\Xi}\widehat{J})^{-1}$.

Note that the pigeonhole bootstrap does not require any additional condition, above that ensuring the $\sqrt{\underline{C}}$-asymptotic normality of the GMM estimator and the fact that $H$ is positive definite. Hence, Theorem (ref) implies for instance that under the conditions displayed above, the pigeonhole bootstrap is valid for quantile IV regressions under multiway clustering.

Monte Carlo Simulations

We now investigate the finite sample properties of the different inference strategies we have considered. We study the coverage rate of confidence intervals based on either asymptotic normality or the pigeonhole bootstrap. We consider the following different cases:

For each scenario, we compute four confidence intervals. The first three are based on the asymptotic normality of $\widehat{\theta}$ and the consistent estimators $\widehat{V}_1$, $\widehat{V}_2$ and $\widehat{V}_{\text{cgm}}$ of the asymptotic variance. The fourth is based on the pigeonhole bootstrap. As explained above, the variance estimator $\widehat{V}_1$ is very easy to compute with popular econometric softwares such as Stata or R, as it satisfies: $$\widehat{V}_1=\underline{C}\sum_{i=1}^k\widehat{\Sigma}_i,$$ where $\widehat{\Sigma}_{i}$ is the clustered estimated variance with respect to the $i$-th dimension of clustering.\footnote{The term $\underline{C}$ accounts for the fact that $\widehat{V}_1$ estimates the asymptotic variance of $\widehat{\theta}$ rather than its variance.} In other words, one has only to compute $k$ variances under one-way clustering and add them to get a consistent estimate of variance under multiway clustering. $\widehat{V}_{\text{cgm}}$ takes a similar form except that one has to consider additional terms, since it is based on the inclusion-exclusion principle (see (ref) above). For instance, with $k=2$, $\widehat{V}_{\text{cgm}}= \widehat{V}_1 - \underline{C}\widehat{\Sigma}_{12}$, where $\widehat{\Sigma}_{12}$ corresponds to the variance under one-way clustering with clusters defined by the intersection of dimensions 1 and 2 (namely, cells with $k=2$). Finally, $\widehat{V}_2$ can also be written as $\widehat{V}_2=\underline{C}\sum_{i=1}^k\widetilde{\Sigma}_i$, but $\widetilde{\Sigma}_i$ does not correspond to the usual estimator of variance under one-way clustering along dimension $i$.

The small-sample correction $C_i/(C_i-1)$ is often used by default for the computation of the clustered variance $\widehat{\Sigma}_i$, and we follow this practice hereafter (also for $\widehat{\Sigma}_{12}$ and $\widetilde{\Sigma}_i$), except in Scenario 2. Hence, in the baseline scenario, we have

$\widehat{\Sigma}_2$ and $\widetilde{\Sigma}_2$ satisfy the same formulas, up to inverting the roles of $C_1$ and $C_2$, and $j_1$ and $j_2$ in the summations. The formulas are identical for the third scenario. In the second scenario, the formulas remain also the same, except that we remove the correction terms $C_i/(C_i-1)$ and $C_1C_2/(C_1C_2-1)$. Finally, the corresponding formulas for the fourth and fifth scenarios are detailed in Appendix (ref). Note that when $\widehat{V}_2$ or $\widehat{V}_{\text{cgm}}$ are negative, we simply set the confidence intervals to the point estimates.

We also compute Efron's percentile bootstrap confidence interval based on the pigeonhole bootstrap presented in Section (ref): $$\text{IC}_{\text{boot}}=\left[q^{*}_{0.025};q^{*}_{0.975}\right],\text{ with } q^*_{\alpha} \text{ the quantile or order } \alpha \text{ of } \theta^{*}|(N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\leq N_{\bm{j}}})_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}}.$$ This confidence interval is valid since the asymptotic distribution of $\widehat{\theta}$ is symmetric. To simulate the distribution of $\theta^{*}|(N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\leq N_{\bm{j}}})_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}}$, we use 1,000 bootstrap replications for each initial sample we draw.

The results are displayed in Table (ref). In the first scenario, the actual coverage when using our preferred estimator of variance and the pigeonhole bootstrap is always very close to the nominal one, even for $\underline{C}$ as small as 5. It is often considered that between 30 and 50 clusters are necessary with one-way clustering to get reliable confidence intervals bertrand2004, cameron2015. Here, we find that even when 40% of the variance of the cells is related to cluster shocks, 25 cells resulting from a $5\times5$ design are sufficient to get reliable inference, at least with fixed cell sizes and when the estimator $\widehat{\theta}$ is Gaussian.

For the same design and still $\underline{C}=5$, inference based on the estimator of cameron2011 leads to an actual coverage of around $88\%$, for a nominal coverage of $95\%$. The confidence intervals based on $\widehat{V}_2$ perform poorly with small samples. In a $5\times 5$ design, the actual coverage is only around 62%. This is partly but not entirely due to the fact that for 16% of the simulations, we get a negative estimator of the variance, implying that we do not cover $\theta_0$. On the other hand, and in line with the theory, we do observe that the coverage rate of IC$_2$ converges to 95% as $\underline{C}$ grows.

The results corresponding to the second scenario show that excluding the small-sample correction deteriorates significantly the coverage rates for $\underline{C} = 5$ and also, when considering IC$_2$ and IC$_{\text{cgm}}$, for $\underline{C} = 10$. IC$_1$ is less sensitive to the adjustment than IC$_2$ and IC$_{\text{cgm}}$. The correction does not have a notable influence when $C_1,C_2\geq 30$ for any of the confidence intervals. But overall, our simulations suggest that this small-sample correction is desirable.

Results with binary outcomes are qualitatively similar to our baseline simulations: the coverage rates of $\text{IC}_1$ and $\text{IC}_{\text{boot}}$ are closer to the nominal rate than those of $\text{IC}_2$ and $\text{IC}_{\text{cgm}}$. But quantitatively, the coverage rate of $\text{IC}_2$ and $\text{IC}_{\text{cgm}}$ are even further away from the nominal rate, falling respectively under 57% and 84% in the $5\times 5$ design. On the other hand, the coverage rate of $\text{IC}_1$ and $\text{IC}_{\text{boot}}$ remain close to the nominal rate (93,5% and 95,2% in the $5 \times 5$ design). Contrary to the baseline case, we observe that $\text{IC}_{\text{boot}}$ and $\text{IC}_1$ (for $\underline{C}\geq 10$) are slightly conservative here.

The results of the probit model give rise to similar conclusions. $\text{IC}_2$ performs even worse with $\underline{C}=5$, with more than 75% of the variance estimates being negative, but somewhat better than previously with $\underline{C}\geq 10$. The coverage rates of $\text{IC}_{\text{cgm}}$ are very close to those observed in the third scenario. Finally, $\text{IC}_1$ and $\text{IC}_{\text{boot}}$ are again closer to the nominal coverage rates, even if they tend to be slightly more conservative than previously.

Finally, our results with three-way clustering show again the very good performance of IC$_1$ and IC$_{\text{boot}}$ for $\underline{C}$ as small as 3. They also emphasize that even with a normal sample average, IC$_2$ and IC$_{\text{cgm}}$ can still severely undercover with a small number of clusters. In particular, neglecting asymptotically negligible terms as done in $\widehat{V}_2$ leads to 85% of negative estimates with $\underline{C}=3$.

Overall, our results suggest that contrary to IC$_2$ and, to a lesser extent, IC$_{\text{cgm}}$, IC$_1$ and IC$_{\text{boot}}$ may be generally reliable, even with few clusters. As explained above, another advantage of IC$_1$ is that it is even simpler to compute than IC$_{\text{cgm}}$. IC$_{\text{boot}}$ may also be useful in particular for estimators whose asymptotic variance takes a complicated form.

Conclusion

In this paper, we have shown two weak convergence results under multiway clustering. The first implies not only simple central limit theorems, but also asymptotic normality of various nonlinear estimators. The second implies the general validity of the pigeonhole bootstrap under multiway clustering. We also establish the consistency of three variance estimators. Inference based on either the pigeonhole bootstrap or asymptotic normality and our preferred variance estimator works very well in simulations, with coverage rates close to their nominal values for no more than five clusters in each dimension with two-way clustering, or even three with three-way clustering.