The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
89,079 characters
Asymptotic results under multiway clustering
\title{Asymptotic results under multiway clustering\thanks{We would like to thank Stéphane Bonhomme, Clément de Chaisemartin, Isabelle Méjean, seminar participants at Bristol, CREST and Yale and attendees of the 2018 IAAE conference for helpful comments.}}
\author{Laurent Davezies\thanks{CREST, [email removed]}
\and Xavier D'Haultf\oe uille
\thanks{CREST. [email removed]}
\and Yannick Guyonvarch
\thanks{CREST. [email removed]}}
\maketitle
\bigskip
\begin{abstract}
If multiway cluster-robust standard errors are used routinely in applied economics, surprisingly few theoretical results justify this practice. This paper aims to fill this gap. We first prove, under nearly the same conditions as with i.i.d. data, the weak convergence of empirical processes under multiway clustering. This result implies central limit theorems for sample averages but is also key for showing the asymptotic normality of nonlinear estimators such as GMM estimators. We then establish consistency of various asymptotic variance estimators, including that of \cite{cameron2011} but also a new estimator that is positive by construction. Next, we show the general consistency, for linear and nonlinear estimators, of the pigeonhole bootstrap, a resampling scheme adapted to multiway clustering. Monte Carlo simulations suggest that inference based on our two preferred methods may be accurate even with very few clusters, and significantly improve upon inference based on \cite{cameron2011}.
\medskip
\textbf{Keywords:} Multiway clustering, Empirical processes, Cluster-robust standard errors, Pigeonhole bootstrap, GMM.
\medskip
\textbf{JEL codes:} C13, C15, C21, C23.
\end{abstract}
\newpage
\section{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 \citep{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 \cite{cameron2011},\footnote{According to the Web of Science and Google Scholar, \cite{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. \cite{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.
\medskip
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.
\medskip
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 \cite{donald2007inference}, \cite{ibragimov2010}, \cite{bester2011inference}, \cite{ibragimov2016} and \cite{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 \cite{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.
\medskip
Then, we prove the consistency of three asymptotic variance estimators, including that suggested by \cite{cameron2011}.\footnote{\cite{mackinnon2017} also prove the consistency of the estimator of \cite{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 \cite{cameron2011} when $\underline{C}$ is small.
\medskip
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 \cite{mccullagh2000resampling} and studied by \cite{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.
\medskip
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 \citep{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 \citep[see for instance Lemma 2.3.1 in][]{vanderVaartWellner1996}, using again the Aldous-Hoover representation and an adaptation to our framework of arguments used in \cite{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.
\medskip
The literature on clustering is vast but has mostly focused on linear models under one-way clustering, following the seminal papers of \cite{pfeffermann1981}, \citeauthor{moulton1986} (\citeyear{moulton1986}, \citeyear{moulton1990}) \cite{liang1986} and \cite{arellano1987}. Without being exhaustive, we also refer to \cite{hansen2007}, \cite{cameron2008bootstrap}, \cite{carter2017}, \cite{mackinnon2017wild} and \cite{hansen2017} for more recent contributions.
\medskip
The only papers we are aware of considering multiway clustering are the recent works of \cite{menzel2017} and \cite{mackinnon2017}. \cite{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. \cite{mackinnon2017} focus on linear regressions with two-way clustering. For such models, they show asymptotic normality and the consistency of the variance estimator of \cite{cameron2011}. They also show the validity of a certain wild boostrap in this context.
\medskip
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 \cite{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 \cite{menzel2017} and \cite{mackinnon2017}. Contrary to \cite{menzel2017}, we do not impose cell sizes equal to one, or i.i.d. units within cells. \cite{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.
\medskip
The paper is organized as follows. Section \ref{sec:setup} describes the assumptions we impose on the data generating process and the parameters of interest we consider afterwards. Section \ref{sec:gen_results} provides our main results on the convergence of the empirical process and the pigeonhole bootstrap empirical process. Section \ref{sec:appli} 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{sec:MC} explores through simulations the finite-sample properties of inference based on asymptotic normality or the pigeonhole bootstrap. Section \ref{sec:conclu} concludes. The appendix gathers extensions, additional details on simulations and all the proofs of our results.
\section{The set up}
\label{sec:setup}
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$.
\medskip
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{as:dgp} formalizes all these conditions.
\begin{hyp}\label{as:dgp}
~
\begin{enumerate}
\item The array $(N_{\bm{j}}, (Y_{\ell,\bm{j}})_{\ell\geq 1})_{\bm{j}\geq \boldsymbol{1}}$ is separately exchangeable. Namely, for any $(\pi_1,...,\pi_k)$ $k$-tuple of permutations of $\mathbb{N}$,
$$(N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\geq 1})_{\bm{j}\geq \boldsymbol{1}}\overset{d}{=}(N_{\pi_1(j_1),...,\pi_k(j_k)},(Y_{\ell,\pi_1(j_1),...,\pi_k(j_k)})_{\ell\geq 1})_{\bm{j}\geq \boldsymbol{1}}.$$
\item For any $\bm{c}\geq \boldsymbol{1}$, $(N_{\bm{j}}, (Y_{\ell,\bm{j}})_{\ell\geq 1})_{\boldsymbol{1} \leq \bm{j} \leq \bm{c}}$ is independent of $(N_{\bm{j}'}, (Y_{\ell,\bm{j}'})_{\ell\geq 1})_{\bm{j}'\geq \bm{c}+\boldsymbol{1}}$.
\item $E(N_{\boldsymbol{1}})>0$.
\item The econometrician observes $(N_{\bm{j}},(Y_{\ell,\bm{j}})_{1\leq \ell \leq N_{\bm{j}}})_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}$, with $\underline{C}\rightarrow \infty$ and for all $i=1...k$, $\underline{C}/C_i\rightarrow \lambda_i\geq 0$.
\end{enumerate}
\end{hyp}
To better understand Assumptions \ref{as:dgp}.1 and \ref{as:dgp}.2, consider first for simplicity two-way clustering with $N_{\bm{j}}=1$ almost surely. The data can then be depicted as follows.
\begin{figure}[H]
$$\begin{array}{|c|c|c|c|c|} \hline
&1&2& \cdots &C_2\\
\hline
1& Y_{1,(1,1)} & Y_{1,(1,2)} & \cdots & Y_{1, (1,C_2)} \\
\hline
2& Y_{1,(2,1)}& Y_{1,(2,2)}& \cdots & Y_{1,(2,C_2)} \\
\hline
\vdots & \vdots & \vdots & \vdots & \vdots \\
\hline
C_1& Y_{1,(C_1,1)} & Y_{1,(C_1,2)}& \cdots & Y_{1,(C_1,C_2)} \\
\hline
\end{array}$$
\end{figure}
Assumption \ref{as:dgp}.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{as:dgp}.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.
\medskip
Assumption \ref{as:dgp}.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{as:dgp}.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{as:dgp}.} When combined with Assumption \ref{as:dgp}.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{as:dgp}.2).
\medskip
Let us come back to the general case with possibly $N_{\bm{j}}\neq 1$. Assumption \ref{as:dgp}.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 \cite{carter2017}. Also, $Y_{\ell,\bm{j}}$ may have a different distribution from $Y_{\ell',\bm{j}}$, for $\ell\neq \ell'$.
\medskip
Assumption \ref{as:dgp}.3 only excludes arrays that are almost surely empty. Assumption \ref{as:dgp}.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.
\medskip
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
\begin{equation}
\theta_0 =\frac{\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}Y_{\ell,\boldsymbol{1}}\right)}{\mathbb E(N_{\boldsymbol{1}})}.
\label{eq:expectation}
\end{equation}
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:
\begin{equation}
\widehat{\theta}= \frac{\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}\sum_{\ell=1}^{N_{\bm{j}}}Y_{\ell,\bm{j}}}{\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}N_{\bm{j}}}= \frac{\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} \sum_{\ell=1}^{N_{\bm{j}}}Y_{\ell,\bm{j}}}{\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} N_{\bm{j}}},
\label{eq:sample_avg}
\end{equation}
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{sub:sample_averages} below.
\medskip
More generally, we consider parameters of interest that depend on the unit-level distribution of $Y$, defined by
\begin{equation}
F_Y(y) = \frac{\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}\mathds{1}\{Y_{\ell,\boldsymbol{1}}\leq y\}\right)}{\mathbb E(N_{\boldsymbol{1}})}.
\label{eq:def_F}
\end{equation}
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{sub:nonlinear_functionals_of_the_distribution} below.
\medskip
Finally, we consider in Section \ref{sub:gmm_estimators} moment restrictions at the unit level, rather than at the cell level. Namely, we consider a parameter of interest $\theta_0\in \Theta$ satisfying
\begin{equation}
\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta_0)\right)=0,
\label{eq:GMM}
\end{equation}
for a vector-valued function $m(y,\theta)$. The average parameter defined by \eqref{eq:expectation} is a particular case of \eqref{eq:GMM}, 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.
\medskip
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{sec:gen_results} the asymptotic behavior of empirical processes and their bootstrap counterparts.
\section{Weak convergence results}
\label{sec:gen_results}
\subsection{Empirical processes}
\label{sec:emp_process}
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., \cite{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 \eqref{eq:sample_avg}. 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.
\medskip
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{as:vc} 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}$.
\medskip
We consider the following standard assumptions on the class $\mathcal{F}$ indexing $\mathbb{G}_C$.
\begin{hyp}\label{as:measurability}
$\mathcal{F}$ is a pointwise measurable class of functions.
\end{hyp}
\begin{hyp}\label{as:vc}
The class $\mathcal{F}$ admits an envelope $F$ with either:\\
- $\mathbb E\left[\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}} F\left(Y_{\ell,\boldsymbol{1}}\right)\right)^2\right]<+\infty$ and $\mathcal{F}$ is finite; \\
- or $\mathbb E\left[N_{\boldsymbol{1}}^2\right]<+\infty$, $\mathbb E\left[N_{\boldsymbol{1}}\sum_{\ell=1}^{N_{\boldsymbol{1}}} F\left(Y_{\ell,\boldsymbol{1}}\right)^2\right]<+\infty$ and
\begin{equation*}
\int_0^{+\infty}\sup_Q\sqrt{\log N\left(\varepsilon||F||_{Q,2},\mathcal{F},||.||_{Q,2}\right)}d\varepsilon<+\infty,
\end{equation*}
where the supremum is taken over the set of probability measures with finite support on $\mathcal{Y}$.
\end{hyp}
Assumption \ref{as:measurability} is not necessary but usually imposed \cite[see, e.g.][]{cherno2014,Kato2016} to avoid measurability issues and the use of outer expectations. For further discussion about these classes, we refer to \citeauthor{Kosorok2006} (\citeyear{Kosorok2006}, pp.137-140). Assumption \ref{as:vc} imposes a condition on what is usually referred to as the uniform entropy integral, see, e.g., \cite{vanderVaartWellner1996}. Finiteness of the uniform entropy integral is satisfied by any VC-type class of functions \citep[see][for a definition]{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.
\begin{thm} \label{thm:unifTCL}
Suppose that Assumptions~\ref{as:dgp}-\ref{as:vc} hold. Then the process $\mathbb{G}_C$ converges weakly to a centered Gaussian process $\mathbb{G}$ on $\mathcal{F}$ as $\underline{C}$ tends to infinity. Moreover, the covariance kernel $K$ of $\mathbb{G}$ satisfies:
$$K(f_1,f_2) = \sum_{i=1}^{k}\lambda_i \mathbb{C}ov\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}f_1(Y_{\ell,\boldsymbol{1}}), \sum_{\ell=1}^{N_{\boldsymbol{2}_i}}f_2(Y_{\ell,\boldsymbol{2}_i})\right),$$
where $\boldsymbol{2}_i$ is the $k$-tuple with $2$ in each entry but 1 in entry $i$.
\end{thm}
Theorem \ref{thm:unifTCL} 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{as:vc}, but this is not the only possible restriction on the class of functions. In Appendix \ref{app:smoothness_class}, we show the same result under smoothness restrictions on $\mathcal{F}$ instead of Assumption \ref{as:vc}.
\medskip
Let us summarize the proof of Theorem \ref{thm:unifTCL}. 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 \citep[][for the almost-sure version]{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}$,
\begin{equation}
(N_{\bm{j}}, (Y_{\ell, j})_{\ell\geq 1})= \tau(U_{j_1, 0},U_{0,j_2}, U_{\bm{j}}).
\label{eq:representation}
\end{equation}
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}$.}
\medskip
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.
\medskip
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 \citep[see, e.g.,][Section 2.1.2, for a formal definition]{vanderVaartWellner1996}. For that purpose, we prove a symmetrization lemma similar to Lemma 2.3.1 in \cite{vanderVaartWellner1996}. To do so, we adapt arguments used in the proofs of Theorem 3.1 in \cite{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{as:dgp} 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{as:vc}. By what precedes, this implies the asymptotic equicontinuity of $\mathbb G_C$.
\medskip
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{thm:unifTCL} implies that the asymptotic variance of $\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} S_{\bm{j}}/\Pi_C$ is
\begin{equation}
\sum_{i=1}^k \lambda_i \mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right).
\label{eq:formula_variance}
\end{equation}
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{as:dgp}, each covariance term is positive. Specifically, and considering for simplicity $k=2$, we establish in the proof of Theorem \ref{thm:unifTCL} 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 \eqref{eq:representation}.
\medskip
Now, let us give some intuitions on \eqref{eq:formula_variance}. 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.
\subsection{Pigeonhole bootstrap processes}
\label{sub:boot}
\label{sec:boot}
We now consider the bootstrap counterpart of the weak convergence result in Theorem \ref{thm:unifTCL}. 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 \cite{mccullagh2000resampling} and studied, in the case of the sample mean and for particular models, by \cite{owen2007pigeonhole}. We are, however, not aware of any result concerning the asymptotic validity of the pigeonhole bootstrap for inference. Theorem \ref{thm:boot_unif} below aims to fill this gap.
\medskip
We first recall the principle of the pigeonhole bootstrap:
\begin{enumerate}
\item For each $i\in\{1,...,k\}$, $C_i$ elements are sampled with replacement and equal probability in the set $\{1,...,C_i\}$. For each $j_i$ in this set, let $W^i_{j_i}$ denote the number of times $j_i$ is selected this way.
\item Cell $\bm{j}=(j_1,...,j_k)$ is then selected $W_{\bm{j}}=\prod_{i=1}^k W^i_{j_i}$ times in the bootstrap sample.
\end{enumerate}
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{as:dgp}, 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.
\medskip
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{thm:unifTCL}.
As discussed in, e.g., \citeauthor{vanderVaartWellner1996} (1996, Chapter 3.6), conditional weak convergence in probability amounts to proving
\begin{equation}
\sup_{h\in \text{BL}_1} \left|\mathbb E\left(h(\mathbb{G}_C^{\ast})|\{N_{\bm{j}},(Y_{\ell,\bm{j}})_{\ell\geq 1}\}_{\bm{j}\geq \boldsymbol{1}}\right)-\mathbb E\left(h(\mathbb{G})\right)\right|\stackrel{\mathbb{P}}{\longrightarrow} 0,
\label{eq:conv_boot_proc}
\end{equation}
where $\text{BL}_1$ is the set of bounded and Lipschitz functions from $\ell^{\infty}(\mathcal{F})$ to $\mathbb{R}$.
\begin{thm}
Suppose that Assumptions \ref{as:dgp}-\ref{as:vc} hold. Then $\mathbb{G}^{\ast}_{C}$ converges weakly to $\mathbb{G}$ in probability, namely \eqref{eq:conv_boot_proc} holds.
\label{thm:boot_unif}
\end{thm}
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{thm:boot_unif} follows the same lines as that of Theorem \ref{thm:unifTCL} 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 \citep{kallenberg05} and we prove Lindeberg-Feller conditions for some statistics defined on the exchangeable array.
\section{Applications}
\label{sec:appli}
\subsection{Simple averages and linear models}
\label{sub:sample_averages}
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
\begin{equation}
\widehat{\theta}=\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}}S_{\bm{j}}.
\label{eq:def_theta_hat}
\end{equation}
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 \eqref{eq:sample_avg} and linear models. Provided that $E(S_{\bm{j}}^2)<+\infty$, we have, by Theorem \ref{thm:unifTCL},
\begin{equation}
\sqrt{\underline{C}}\left(\widehat{\theta}-\theta_0\right) \stackrel{d}{\longrightarrow} \mathcal{N}\left\{0,\sum_{i=1}^{k}\lambda_i \mathbb{C}ov\left(S_{\boldsymbol{1}}, S_{\boldsymbol{2}_i}\right)\right\}.
\label{eq:TCL_simple}
\end{equation}
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 \cite{menzel2017} for the construction and analysis of a wild bootstrap procedure for sample averages on such clustered data.}
\medskip
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$:
\begin{equation}
\label{eq:Vch2}
\widehat{V}_2 = \sum_{i=1}^k \frac{\underline{C}}{C_i} \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)'.
\end{equation}
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.
\medskip
To avoid these issues, we suggest to simply add to $\widehat{V}_2$ pairs sharing more than one cluster. Specifically, we consider
\begin{align}
\widehat{V}_1 & = \sum_{i=1}^k\frac{\underline{C}}{C_i} \frac{1}{C_i\prod_{s\neq i}C_s^2}\sum_{(\bm{j},\bm{j}'): \bm{j}_i=\bm{j}'_i} \left(S_{\bm{j}}-\widehat{\theta}\right) \left(S_{\bm{j}'}-\widehat{\theta}\right)' \nonumber \\
& = \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} S_{\bm{j}}-\widehat{\theta}\right)
\left(\frac{1}{\prod_{s\neq i}C_s}\sum_{\bm{j}:j_i=j'_i} S_{\bm{j}}-\widehat{\theta}\right)'. \label{eq:V1_positive}
\end{align}
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{prop:as_var} 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 \eqref{eq:V1_positive} 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{sec:MC} 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 \eqref{eq:CGM} 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.
\medskip
We now compare our two estimators with that proposed by \cite{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
\begin{align*}
\mathbb V(\widehat{\theta}) & = \frac{1}{\Pi_C^2} \sum_{(\bm{j}, \bm{j}') \in \cup_{i=1}^k \mathcal{B}_i} \mathbb{C}ov\left(S_{\bm{j}}, S_{\bm{j}'}\right) \\
& = \sum_{m=1}^k (-1)^m \sum_{1\leq i_1 < ... < i_m \leq k} \frac{1}{\Pi_C^2} \sum_{(\bm{j}, \bm{j}')\in \mathcal{B}_{i_1,...,i_m}} \mathbb{C}ov\left(S_{\bm{j}}, S_{\bm{j}'}\right).
\end{align*}
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}$
\begin{equation}
\widehat{V}_{\text{cgm}} = \underline{C} \sum_{m=1}^k (-1)^m \sum_{1\leq i_1 < ... < i_m \leq k} \frac{1}{\Pi^2_C} \sum_{(\bm{j}, \bm{j}')\in \mathcal{B}_{i_1,...,i_m}} \left(S_{\bm{j}}-\widehat{\theta}\right) \left(S_{\bm{j}'}-\widehat{\theta}\right)'.
\label{eq:CGM}
\end{equation}
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 \cite{cameron2011} for more details. As with $\widehat{V}_1$, the appeal of Formula \eqref{eq:CGM} 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$.
\medskip
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
\begin{align}
\widehat{V}_1 & = \underline{C} \sum_{i_1=1}^2 \frac{1}{\Pi^2_C} \sum_{(\bm{j}, \bm{j}')\in \mathcal{B}_{i_1}} \left(S_{\bm{j}}-\widehat{\theta}\right) \left(S_{\bm{j}'}-\widehat{\theta}\right)', \nonumber \\
\widehat{V}_{\text{cgm}} & = \widehat{V}_1 - \frac{\underline{C}}{\Pi^2_C} \sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} \left(S_{\bm{j}}-\widehat{\theta}\right) \left(S_{\bm{j}}-\widehat{\theta}\right)'.\label{eq:compar_V1_Vcgm}
\end{align}
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 \eqref{eq:compar_V1_Vcgm}. 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 \eqref{eq:compar_V1_Vcgm} 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.
\medskip
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{prop:as_var} 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.
\begin{prop}
\label{prop:as_var}
Suppose that Assumption \ref{as:dgp} holds and $\mathbb{E}\left[S_{\boldsymbol{1}}S_{\boldsymbol{1}}'\right]<+\infty$. Then $\widehat{V}_1$, $\widehat{V}_2$ and $\widehat{V}_{\text{cgm}}$ are consistent for $V$. Moreover, if $V$ is positive definite, we have, for any $k\in \{1,2,\text{cgm}\}$ and $\alpha\in (0,1)$,
$$\lim_{\underline{C}\rightarrow+\infty} \mathbb{P}\left(R_k^{1-\alpha}\ni\theta_0\right)=1-\alpha$$
\end{prop}
\medskip
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{as:dgp}, $\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 \cite{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
\begin{equation}
S_{\bm{j}}=(U_{j_1,0} - \mathbb E(U_{j_1,0})) (U_{0,j_2} - \mathbb E(U_{0,j_2}))+U_{\bm{j}}.
\label{eq:non_normal_example}
\end{equation}
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 \eqref{eq:non_normal_example} holds. We refer to \cite{menzel2017}, Example 1.6, for more details.
\medskip
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{sub:boot} 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{thm:boot_unif}.
\begin{prop}
\label{prop:boot_mean}
Suppose that Assumption \ref{as:dgp} holds, $\mathbb{E}\left[S_{\boldsymbol{1}}S_{\boldsymbol{1}}'\right]<+\infty$ and $V$ is positive definite. Then $$\lim_{n\rightarrow\infty} \mathbb{P}\left(R_{\text{boot}}^{1-\alpha} \ni \theta_0\right) = 1-\alpha.$$
\end{prop}
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.
\medskip
We now discuss how Propositions \ref{prop:as_var} and \ref{prop:boot_mean} extend to other parameters of interest. First, let us consider $\theta_0=E(S_{\boldsymbol{1}})/E(N_{\boldsymbol{1}})$, as in Section \ref{sec:setup}. Assume that $\mathbb E(S_{\boldsymbol{1}}^2)<\infty$ and $\mathbb E(N_{\boldsymbol{1}}^2)<\infty$. By Theorem \ref{thm:unifTCL} applied to $\mathcal{F}=\{\text{Id},1\}$ and the delta method, we have
\begin{equation}
\sqrt{\underline{C}}\left(\frac{\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}}S_{\bm{j}}}{\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} N_{\bm{j}}}-\theta_0\right) = \frac{\sqrt{\underline{C}}}{\Pi_C}\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} T_{\bm{j}} +o_p(1).
\label{eq:lin_ratio}
\end{equation}
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 \eqref{eq:Vch2}, \eqref{eq:V1_positive} and \eqref{eq:CGM} 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{prop:as_var}, 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 \cite{vanderVaart2000}. More generally, Propositions \ref{prop:as_var} and \ref{prop:boot_mean} 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})$.
\medskip
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
\begin{align}
\theta_0 & = \mathbb E\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}} X_{\ell,\boldsymbol{1}} X_{\ell,\boldsymbol{1}}' \right]^{-1} \mathbb E\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}} X_{\ell,\boldsymbol{1}} \tilde{Y}_{\ell,\boldsymbol{1}} \right] \label{eq:def_theta_lin} \\
\widehat{\theta} & = \left(\frac{1}{\Pi_C}\sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} X_{\ell,\bm{j}} X_{\ell,\bm{j}}' \right)^{-1} \left(\frac{1}{\Pi_C} \sum_{\ell=1}^{N_{\boldsymbol{1}}} \sum_{\boldsymbol{1} \leq \bm{j} \leq \bm{C}} X_{\ell,\bm{j}} \tilde{Y}_{\ell,\bm{j}}\right). \label{eq:def_thetahat_lin}
\end{align}
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$.
\begin{prop}\label{prop:linear}
Suppose that Assumption \ref{as:dgp} holds with $Y_{\ell,\bm{j}}=(\tilde{Y}_{\ell,\bm{j}},X'_{\ell,\bm{j}})'$. Suppose also
$\mathbb E\left((\sum_{\ell=1}^{N_{\boldsymbol{1}}}|Y_{\ell,\boldsymbol{1}}|^2)^2\right)<+\infty$ and $\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}X_{\ell,\boldsymbol{1}} X'_{\ell,\boldsymbol{1}}\right)$ non-singular. Let $\theta_0$ and $\widehat{\theta}$ be defined by \eqref{eq:def_theta_lin} and \eqref{eq:def_thetahat_lin}. Then
$$\sqrt{\underline{C}}\left(\widehat{\theta}-\theta_0\right)\stackrel{d}{\longrightarrow}\mathcal{N}(0,V),$$
with $V=J^{-1}HJ^{-1}$, $J=\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}X_{\ell,\boldsymbol{1}} X'_{\ell,\boldsymbol{1}}\right)$ and
$$H=\sum_{i=1}^k\lambda_i \mathbb E\left[\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}X_{\ell,\boldsymbol{1}} u_{\ell,\boldsymbol{1}}\right)\left(\sum_{\ell=1}^{N_{\boldsymbol{2}_i}}u_{\ell,\boldsymbol{2}_i} X'_{\ell,\boldsymbol{2}_i} \right)\right].$$
\end{prop}
Next, we show similar results as in Propositions \ref{prop:as_var} and \ref{prop:boot_mean}. 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}$.
\begin{prop}\label{prop:linear2}
Under the assumptions of Proposition \ref{prop:linear}, $\widehat{V} \stackrel{\mathbb{P}}{\longrightarrow} V$. If $V$ is positive definite, inference based on either asymptotic normality and $\widehat{V}$, or the pigeonhole bootstrap, is valid.
\end{prop}
Propositions \ref{prop:linear} and \ref{prop:linear2} complement the results of \cite{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{prop:linear2} also shows the consistency of a new, positive, variance estimator and the asymptotic validity of the pigeonhole bootstrap in this context of linear models.
\subsection{Nonlinear functionals of the distribution}
\label{sub: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{sec:gen_results} can be applied to such smooth, nonlinear functionals of the empirical distribution. Let $F_Y$ be defined as in \eqref{eq:def_F} 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 \citep[see, e.g., Example 20.15 in][]{vanderVaart2000} or the nonlinear difference-in-difference estimand of \cite{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..
\medskip
We consider the plug-in estimator $\widehat{\theta}=g(\widehat{F_Y})$ of $\theta_0$, with
\begin{equation}
\widehat{F_Y}(y) = \frac{\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}}\sum_{\ell=1}^{N_{\bm{j}}}\mathds{1}\{Y_{\ell,\bm{j}}\leq y\}}{\sum_{\boldsymbol{1}\leq \bm{j}\leq \bm{C}} N_{\bm{j}}}.
\label{eq:def_Fhat}
\end{equation}
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{prop:functional_dm} 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.
\begin{prop}
\label{prop:functional_dm}
Suppose that $\theta_0=g(F_Y)$ and $\widehat{\theta}=g(\widehat{F_Y})$, where $F_Y$ and $\widehat{F_Y}$ are defined respectively by \eqref{eq:def_F} and \eqref{eq:def_Fhat} and $g$ is Hadamard differentiable at $F_Y$ tangentially to $\mathbb{D}_0$. Suppose also that Assumption \ref{as:dgp} holds and $\mathbb E(N_{1,\boldsymbol{1}}^2)<+\infty$. Then:
\begin{enumerate}
\item $\sqrt{\underline{C}}(\widehat{F_Y}-F_Y)$ converges weakly, as a process indexed by $y$, to a Gaussian process $\mathbb{G}$ with kernel $K$ satisfying
\begin{equation}
K(y_1,y_2)=\frac{1}{\mathbb E(N_{1,\boldsymbol{1}})^2} \sum_{i=1}^{k}\lambda_i \mathbb{C}ov\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}\mathds{1}\{Y_{\ell,\boldsymbol{1}}\leq y_1\}, \sum_{\ell=1}^{N_{\boldsymbol{2}_i}}\mathds{1}\{Y_{\ell,\boldsymbol{2}_i}\leq y_2\}\right).
\label{eq:kernel_had_diff}
\end{equation}
\item If $\mathbb{G} \in \mathbb{D}_0$ with probability one,
$$\sqrt{\underline{C}} \left(\widehat{\theta}-\theta_0\right) \rightarrow \mathcal{N}(0,\mathbb V(g'_{F_Y}(\mathbb{G}))).$$
\item If $\mathbb{G} \in \mathbb{D}_0$ with probability one and $\mathbb V(g'_{F_Y}(\mathbb{G}))$ is positive definite, $$\lim_{n\rightarrow\infty} \mathbb{P}\left(R_{1-\alpha}^{\text{boot}} \ni \theta_0\right) = 1-\alpha.$$
\end{enumerate}
\end{prop}
The first part follows from Theorem \ref{thm:unifTCL}, and a linearization of the ratio akin to \eqref{eq:lin_ratio}. The second follows from the first part and the functional delta method \citep[see, e.g.,][Theorem 20.8]{vanderVaart2000}. The third part is a direct consequence of Theorem \ref{thm:boot_unif} and the functional delta method for the bootstrap \citep[see, e.g.,][Theorem 3.9.11]{vanderVaartWellner1996}.
\medskip
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$ \citep[see, e.g.][Lemma 21.3]{vanderVaart2000}. Moreover, we prove in Appendix \ref{app:example_quantile} 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{prop:functional_dm} 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.
\medskip
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{sub:sample_averages} 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.
\subsection{GMM estimators}
\label{sub:gmm_estimators}
Finally, we consider parameters defined by the moment restrictions \eqref{eq:GMM}, 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.
\begin{hyp}
~
\begin{enumerate}
\item $\theta_0$ belongs to the interior of $\Theta$, a compact subset of $\mathbb{R}^p$.
\item $\mathbb E\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right]=0$ if and only if $\theta=\theta_0$.
\item For any $\theta\in \Theta$ we have $\lim_{\theta'\rightarrow \theta}\mathbb E\left[\left|\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta')-\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right|^2\right]=0$.
\item $\theta\mapsto \mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right)$ is differentiable at $\theta_0$ with a jacobian matrix $J$ of rank $p$.
\item For all $s=1,...,L$ the class $\mathcal{F}_s=\{y\mapsto m_s(y,\theta):\theta \in \Theta\}$ fulfills Assumptions \ref{as:measurability}-\ref{as:vc}.
\item $\widehat{\Xi}$ is a sequence of random symmetric matrices of size $L$ tending in probability to $\Xi$, which is positive definite.
\end{enumerate}
\label{as:GMM}
\end{hyp}
Assumptions \ref{as:GMM}.1, \ref{as:GMM}.2 and \ref{as:GMM}.6 are standard. Assumption \ref{as:GMM}.3, combined with \ref{as:GMM}.1 and \ref{as:GMM}.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{as:GMM}.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{as:GMM}.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 \cite{newey1994large}. Finally, by Theorem \ref{thm:unifTCL}, Assumption \ref{as:GMM}.5 ensures the stochastic equicontinuity condition \citep[e.g., Condition (v) in Theorem 7.2 of][]{newey1994large}, which together with \ref{as:GMM}.4, is key to obtain $\sqrt{\underline{C}}$-asymptotic normality of $\widehat{\theta}$.
\medskip
To illustrate that Assumption \ref{as:GMM} 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{sub:quantile_IV} that Assumptions \ref{as:GMM}.3-\ref{as:GMM}.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{as:quantile_IV} in Appendix \ref{sub:quantile_IV}.}
\begin{thm}
Suppose that Assumptions \ref{as:dgp} and \ref{as:GMM} hold. Then $\widehat{\theta}$ is well-defined with probability approaching one and
$$\sqrt{\underline{C}}\left(\widehat{\theta}-\theta_0\right)\stackrel{d}{\longrightarrow}\mathcal{N}\left(0,V_0\right),$$
where $V_0 = (J'\Xi J)^{-1} J' \Xi H \Xi J (J'\Xi J)^{-1}$ and
$$H = \sum_{i=1}^k \lambda_i \mathbb E\left[\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta_0)\right)\left(\sum_{\ell=1}^{N_{\boldsymbol{2}_i}}m(Y_{\ell,\boldsymbol{2}_i},\theta_0)\right)'\right].$$
\label{thm:AN_GMM}
\end{thm}
To our knowledge, Theorem \ref{thm:AN_GMM} is the first result on the asymptotic normality of GMM estimators under multiway clustering. Theorem \ref{thm:AN_GMM} 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.
\medskip
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.
\begin{hyp}
~
\begin{enumerate}
\item The jacobian matrix $J$ of $\theta\mapsto\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}m(Y_{\ell,\boldsymbol{1}},\theta)\right)$ at $\theta_0$ admits the following representation
$$J=\mathbb E\left(\sum_{\ell=1}^{N_{\boldsymbol{1}}}d(Y_{\ell,\boldsymbol{1}},\theta)\right)$$
for some matrix-valued function $d(.,.)=\left(d_{r,s}(.,.)\right)_{1\leq r\leq p,1\leq s\leq L}$.
\item For all $(r,s)\in\{1,...,p\}\times\{1,...,L\}$, the class $\mathcal{F}_{r,s}=\{y\mapsto d_{r,s}(y,\theta):\theta \in \Theta\}$ fulfills Assumption \ref{as:measurability} and admits an envelope $F_{r,s}$ such that $E[\sum_{\ell=1}^{N_{\boldsymbol{1}}}F_{r,s}(Y_{\ell, \boldsymbol{1}})]<+\infty$, and for any $\varepsilon>0$, $\sup_{Q} N(\varepsilon ||.||_{Q,1},\mathcal{F}_{r,s},||.||_{Q,1})<\infty$ where the supremum is taken over the set of probability measures with finite support on $\mathcal{Y}$.
\item $\lim_{\theta'\to\theta_0}\mathbb{E}\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}d(Y_{\ell,\boldsymbol{1}},\theta')\right]=\mathbb{E}\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}d(Y_{\ell,\boldsymbol{1}},\theta_0)\right]$.
\item For every $i=1,...,k$, $$\lim_{\theta'\to\theta_0}\mathbb{E}\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}\sum_{\ell=1}^{N_{\boldsymbol{2}_i}}m(Y_{\ell,\boldsymbol{1}},\theta')m(Y_{\ell,\boldsymbol{2}_i},\theta')\right]=\mathbb{E}\left[\sum_{\ell=1}^{N_{\boldsymbol{1}}}\sum_{\ell=1}^{N_{\boldsymbol{2}_i}}m(Y_{\ell,\boldsymbol{1}},\theta_0)m(Y_{\ell,\boldsymbol{2}_i},\theta_0)\right].$$
\end{enumerate}
\label{as:varGMM}
\end{hyp}
Assumption \ref{as:varGMM}.1 refines Assumption \ref{as:GMM}.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{as:varGMM}.2, the condition on the classes are of Glivenko-Cantelli type, and weaker than Assumption \ref{as:vc}. The continuity conditions in Points 3 and 4 are similar to Assumption \ref{as:GMM}.3, but are imposed on different functions related to $J$ and $H$ rather than on the moment conditions themselves.
\medskip
We now define our estimator of $V_0$, which is based on estimators of $J$ and $H$. Given Assumption \ref{as:varGMM}.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
\begin{align*}
&\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}}}m\left(Y_{\ell,\bm{j}},\widehat{\theta}\right)\right)
\left(\frac{1}{\prod_{s\neq i}C_s}\sum_{\bm{j}:j_i=j'_i} \sum_{\ell=1}^{N_{\bm{j}}}m\left(Y_{\ell,\bm{j}},\widehat{\theta}\right) \right)'.
\end{align*}
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}$.
\begin{thm}
\label{thm:var_gmm}
Assume that Assumptions \ref{as:dgp} and \ref{as:GMM} hold and $H$ is positive definite. Then:
\begin{enumerate}
\item If Assumption \ref{as:varGMM} holds as well, $\widehat{V}\stackrel{\mathbb{P}}{\longrightarrow} V_0$ and confidence regions and tests on $\theta_0$ based on asymptotic normality and $\widehat{V}$ are asymptotically valid.
\item Confidence regions and tests on $\theta_0$ based on the pigeonhole bootstrap are asymptotically valid.
\end{enumerate}
\end{thm}
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{thm:var_gmm} implies for instance that under the conditions displayed above, the pigeonhole bootstrap is valid for quantile IV regressions under multiway clustering.
\section{Monte Carlo Simulations}
\label{sec:MC}
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:
\begin{enumerate}
\item ``two-way, Gaussian'': our baseline scenario is a two-way balanced design ($C_1=C_2=\underline{C}$) with one observation per cell ($N_{\bm{j}}=1$). Each $Y_{1,\bm{j}}$ is drawn in a standard Gaussian distribution, but the variance due to cell shocks only represent 60\% of the total variance, whereas row and column shocks represent 20\% of the variance each:
\begin{equation}\label{eq:DGPsim}
Y_{1,(j_1,j_2)}=\frac{1}{\sqrt{5}}\left(U_{j_1,0}+U_{0,j_2}+\sqrt{3}U_{j_1,j_2}\right), \quad \left(U_{j_1,0},U_{0,j_2},U_{j_1,j_2}\right)\sim \mathcal{N}(0,I_3).
\end{equation}
The parameter of interest is $\theta_0=\mathbb E(Y_{1,\boldsymbol{1}})$ and we consider $C_1=C_2$ taking values in $\{5,10,30,50,100\}$.
\item ``two-way, w/o adjust'': this scenario is as the baseline, except that we compute variance estimators without including the finite-sample corrections described below. The purpose is to investigate the effects of this correction in practice.
\item ``two-way, binary'': this scenario is as the baseline, except that $\theta_0=\mathbb E(\mathds{1}\left\{Y_{1,\boldsymbol{1}}>0\right\})$. The goal is to investigate whether accuracy of inference in our baseline scenario is driven by the fact that $\widehat{\theta}$ itself is normally distributed.
\item ``two-way, probit'': in this scenario, we consider a simple probit model, with random cell sizes. Namely, we suppose that the $(N_{\bm{j}})_{\bm{j}\geq \boldsymbol{1}}$ are independent, with $N_{\bm{j}}\sim 1+\mathcal{P}(5)$. Our outcome variable is then defined by
$$\tilde{Y}_{\ell, (j_1, j_2)}=\mathds{1}\left\{\beta_0+\theta_0 X_{\ell,(j_1, j_2)}+\frac{1}{\sqrt{6}}\left(U_{(j_1,0)}+U_{(0,j_2)}+U_{(j_1,j_2)}\right) + \frac{U_{\ell,\bm{j}}}{\sqrt{2}}>0\right\},$$
where the $(X_{\ell,\bm{j}})_{\ell\geq 1,\bm{j}\geq\boldsymbol{1}}$, $(U_{\bm{j}})_{\bm{j}\in \mathbb{N}^2}$ and $(U_{\ell,\bm{j}})_{\ell\geq 1,\bm{j}\geq\boldsymbol{1}}$ are mutually independent and standard normal variables (as above, the $(U_{\bm{j}})_{\bm{j}\in \mathbb{N}^2}$ and $(U_{\ell,\bm{j}})_{\ell\geq 1,\bm{j}\geq\boldsymbol{1}}$ are assumed unobserved). $\theta_0$ is again our parameter of interest and $(\beta_0,\theta_0)=(0,1)$. In these simulations $(\widehat{\beta},\, \widehat{\theta})$ is the pseudo-maximum likelihood estimator of $(\beta_0,\, \theta_0)$, i.e. the usual probit estimator obtained on the pooled sample. The first aim of this scenario is to study the sensitivity of inference to the non-linearity of the estimator. The second aim is to investigate the sensitivity of inference to randomness in cell sizes $N_{\bm{j}}$.
\item ``three-way, Gaussian'': this scenario differs from the baseline in that we consider three-way clustering. As in the baseline, $N_{\bm{j}}=1$, $Y_{\bm{j}}$ is Gaussian and 60\% of the variance is due to cell shocks. 6,67\% of the variance is due to shocks specific to dimension 1, 2 or 3 of the clustering. The remaining variance is due to shocks common to dimensions 1 and 2, 2 and 3, and 1 and 3. Specifically,
$$Y_{1,\bm{j}}=\frac{1}{\sqrt{15}}\left(U_{(j_1,0,0)}+U_{(0,j_2,0)}+U_{(0,0,j_3)} + U_{(j_1,j_2,0)}+U_{(j_1,0,j_3)} + U_{(0,j_2,j_3)} + 3U_{\bm{j}}\right),$$
where the $(U_{\bm{j}})_{\bm{j} \in \mathbb{N}^3}$ are independent standard normal variables. We consider $C_1=C_2=C_3$ taking values in $\{3,5,10,30,50\}$. These values were chosen so as to correspond roughly to the same number of cells as in the corresponding cases under our baseline scenario.
\end{enumerate}
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 \eqref{eq:CGM} 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$.
\medskip
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
\begin{align*}
\widehat{\Sigma}_1 & =\frac{C_1}{C_1-1} \frac{1}{\Pi_C^2}\sum_{j_1=1}^{C_1} \left(\sum_{j_2=1}^{C_2}\sum_{\ell=1}^{N_{(j_1,j_2)}}(Y_{\ell,(j_1,j_2)}-\widehat{\theta})\right)^2 \\
\widetilde{\Sigma}_1& =\frac{C_1}{C_1-1} \frac{1}{C_1^2C_2(C_2-1)}\sum_{j_1=1}^{C_1} \left(\sum_{\substack{1\leq j_2,j_2'\leq C_2\\
j_2'\neq j_2}}(Y_{1,(j_1,j_2)}-\widehat{\theta})(Y_{1,(j_1,j_2')}-\widehat{\theta})\right),\\
\widehat{\Sigma}_{12}&=\frac{C_1C_2}{C_1C_2-1} \frac{1}{\Pi_C^2}\sum_{\boldsymbol{1} \leq \bm{j}\leq \bm{C}} \left(\sum_{\ell=1}^{N_{\bm{j}}}(Y_{\ell,\bm{j}}-\widehat{\theta})\right)^2.
\end{align*}
$\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{sec:additional_details_on_the_simulations}. Note that when $\widehat{V}_2$ or $\widehat{V}_{\text{cgm}}$ are negative, we simply set the confidence intervals to the point estimates.
\medskip
We also compute Efron's percentile bootstrap confidence interval based on the pigeonhole bootstrap presented in Section \ref{sec:boot}:
$$\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.
\medskip
The results are displayed in Table \ref{tab:MC_results}. 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 \citep{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.
\medskip
For the same design and still $\underline{C}=5$, inference based on the estimator of \cite{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.
\medskip
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.
\medskip
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.
\medskip
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.
\medskip
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$.
\medskip
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.
\begin{table}[H]
\caption{Coverage rates on the five scenarios (nominal coverage rate: 0.95)}
\label{tab:MC_results}
\begin{center}
\begin{tabular}{c c c c ccccc}
\hline \hline
&& & &\multicolumn{5}{c}{$\underline{C}$}\\
Scenario& & Interval &&5&10&30&50&100 \\
\hline
& & $\text{IC}_{\text{boot}}$ & &0.929& 0.94& 0.948& 0.952& 0.951 \\
two-way, Gaussian && $\text{IC}_1$ & &0.935 &0.939 &0.949& 0.957& 0.955 \\
&& $\text{IC}_2$ & &$\underset{\text{\tiny{[16.2\%]}}}{0.653^*}$& 0.87& 0.926& 0.943& 0.949\\
&& $\text{IC}_{\text{cgm}}$ & &$\underset{\text{\tiny{[0.5\%]}}}{0.875^*}$ &0.916 &0.936& 0.952& 0.952\\[3mm]
\cline{3-9}
&& $\text{IC}_1$ & &0.904& 0.933& 0.945& 0.955& 0.952 \\
two-way, w/o adjust.&& $\text{IC}_2$ & &$\underset{\text{\tiny{[16.2\%]}}}{0.626^*}$ &0.853& 0.922& 0.942& 0.949\\
&& $\text{IC}_{\text{cgm}}$ & &$\underset{\text{\tiny{[1.4\%]}}}{0.816^*}$ &0.897& 0.93& 0.945& 0.949 \\[3mm]
\cline{3-9}
&& $\text{IC}_{\text{boot}}$ & &0.952& 0.97& 0.955& 0.951& 0.953\\
two-way, binary && $\text{IC}_1$ & &0.937 &0.959& 0.957 &0.955& 0.952\\
&&$\text{IC}_2$&&$\underset{\text{\tiny{[31.8\%]}}}{0.561^*}$& $\underset{\text{\tiny{[0.8\%]}}}{0.84^*}$ & 0.925 &0.926 &0.946\\
&& $\text{IC}_{\text{cgm}}$ & &$\underset{\text{\tiny{[0.8\%]}}}{0.837^*}$ &0.921& 0.94& 0.944& 0.948\\[3mm]
\cline{3-9}
&& $\text{IC}_{\text{boot}}$ & &0.938& 0.977& 0.982& 0.976& 0.964\\
two-way, probit && $\text{IC}_1$ & &0.97& 0.977& 0.978& 0.977& 0.959\\
&&$\text{IC}_2$&&$\underset{\text{\tiny{[77.2\%]}}}{0.165^*}$ &$\underset{\text{\tiny{[56.9\%]}}}{0.295^*}$& $\underset{\text{\tiny{[6.8\%]}}}{0.697^*}$& $\underset{\text{\tiny{[0.1\%]}}}{0.829^*}$& 0.898\\
&& $\text{IC}_{\text{cgm}}$ & &$\underset{\text{\tiny{[9.0\%]}}}{0.755^*}$&$\underset{\text{\tiny{[1.5\%]}}}{0.872^*}$ &0.925 &0.935& 0.936\\[3mm]
\cline{3-9}
& &&&\multicolumn{5}{c}{$\underline{C}$}\\
& & & &3&5&10&15&20 \\
\cline{5-9}
&& $\text{IC}_{\text{boot}}$ & &0.958 &0.966 &0.960 &0.956 &0.957\\
three-way, Gaussian && $\text{IC}_1$ & &0.942 &0.956 &0.957 &0.952 &0.958\\
&& $\text{IC}_2$ & &$\underset{\text{\tiny{[85.4\%]}}}{0.096^*}$ &$\underset{\text{\tiny{[25.5\%]}}}{0.484^*}$ &0.86 & 0.914 & 0.925\\
&& $\text{IC}_{\text{cgm}}$ & &$\underset{\text{\tiny{[5.6\%]}}}{0.769^*}$ &$\underset{\text{\tiny{[0.5\%]}}}{0.859^*}$ &0.919 &0.934 &0.937\\
\hline \hline
\end{tabular}
\end{center}
\begin{tablenotes}
Coverage rate estimated on 1000 simulations. The bootstrap confidence intervals are based on 1000 bootstrap samples. ${}^*$ indicates that some estimated variance were negative, in which case the share of negative variance is reported in brackets below. When an estimated variance is negative, we set the corresponding confidence interval to the point estimate.
\end{tablenotes}
\end{table}
\section{Conclusion}
\label{sec:conclu}
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.
\newpage
\bibliography{biblio}
\newpage