Asymptotic Theory for Clustered Samples

Asymptotic Theory for Clustered Samples

$ \medskip} \title{\textbf{Asymptotic Theory for Clustered Samples}} \date{February 2019$\footnote{We thank the Co-Editor Han Hong and two referees for helpful comments on a previous version, and Morten Nielsen and James MacKinnon for valuable conversations and suggestions.}$}

Introduction

Clustered samples are widely used in current applied econometric practice. Despite this dominance, there is little formal large-sample theory for estimation and inference. This paper provides such a foundation. We develop a complete, rigorous, and easily-interpretable asymptotic distribution theory for the \textquotedblleft large number of clusters\textquotedblright \ framework. Our theory allows heterogeneous and growing cluster sizes, but requires that the number of clusters $G$ grows with sample size $n$. Our core theory provides a weak law of large numbers (WLLN), central limit theorem (CLT), and consistent clustered variance estimation for clustered sample means. We also provide uniform laws of large numbers and uniform consistent clustered variance estimation appropriate for the distribution theory of nonlinear econometric estimators.

We apply this core theory to develop large sample distribution theory for standard econometric estimators: linear least-squares, 2SLS, MLE, and GMM. For each, we provide conditions for consistent estimation, asymptotic normality, consistent covariance matrix estimation, and asymptotic distributions for t-ratios and Wald statistics. The theory provided in this paper is the first formal theory for such econometric estimators allowing for clustered dependence.

Our assumptions are minimal, requiring only uniform integrability for the WLLN and squared uniform integrability for the CLT and clustered covariance matrix estimators, plus the requirement that individual clusters are asymptotically negligible. Our results show that there are inherent trade-offs in the conditions between the allowed degree of heterogeneity in cluster sizes and the number of finite moments.\ These trade-offs are least restrictive for the WLLN, are more restrictive for the CLT and consistent cluster covariance matrix estimation, and are strongest for CLTs applied to clustered second moments. These trade-offs do not arise in the independent sampling context.

We show that under clustering the convergence rate depends on the degree of clustered dependence. Convergence rates may equal the square root of the sample size, the square root of the number of clusters, be a rate in between these two, or even slower than both. Our assumptions and theory allow for these possibilities. This is in contrast to the existing literature, which imposes specific rate assumptions. One useful finding is that the rate does not need to be known by the user; the asymptotic distribution of t-ratios and Wald statistics does not depend on the underlying rate of convergence. This generalizes similar results in C. Hansen (2007) and related results in Tabord-Meehan (2018).

This paper makes the following technical contributions. We show that the key to extending the classical WLLN and CLT to cluster-level data is developing uniform integrability bounds for cluster sums. To allow for arbitrary within-cluster dependence, this means that such bounds will be scaled by cluster sizes. This leads to bounds on the degree of cluster size heterogeneity which can be allowed under cluster dependence. Some of the most difficult technical work presented here is the extension of classical results to clustered covariance matrix estimators. These are not sample averages, but rather average across clusters of squared cluster sums. Handling such estimators requires a new technical treatment.

Clustered dependence in econometrics dates to the work of Moulton (1986, 1990), Liang and Zeger (1986), and in particular Arellano (1987), who proposed the popular cluster-robust covariance matrix estimator. The method was popularized by the implementation in Stata by Rogers (1994) and the widely-cited paper of Bertrand, Duflo and Mullainathan (2004). Surveys can be found in Wooldridge (2003), Cameron and Miller (2011, 2015), MacKinnon (2012, 2016), and textbook treatments in Angrist and Pischke (2009) and Wooldridge (2010).

The \textquotedblleft large $G$\textquotedblright\ asymptotic theory develops normal approximations under the assumption that $G\rightarrow\infty$. The earliest treatment appears in White (1984). Wooldridge (2010) asserts a distribution theory under the assumption that the cluster sizes are fixed. C. Hansen (2007) provides two sets of asymptotic results, including both $\sqrt{G}$ and $\sqrt{n}$ convergence rates under two distinct assumptions on the rate of convergence of the estimation variance. His results are derived under the assumption that all clusters are identical in size. Carter, Schnepel and Steigerwald (2017) provided asymptotic results allowing for heterogeneous clusters, but their results are limited by atypical regularity conditions. Independently of this paper, Djogbenou, MacKinnon, and Nielsen (2018) have provided a rigorous asymptotic theory for heterogeneous clusters, with similar but stronger regularity conditions than ours. Their primary focus is theory for regression wild bootstrap, while our focus is regularity conditions for general econometric estimators.

An alternative to the \textquotedblleft large $G$\textquotedblright \ asymptotic is the \textquotedblleft fixed $G$\textquotedblright\ framework, which leads to a non-normal inference theory. Contributions to this literature include C. Hansen (2007), Bester, Conley and C. Hansen (2011), and Ibragimov and M\"{u}eller (2010, 2016). A related paper is Conley and Taber (2011) which provide an asymptotic theory under the assumption of a small number of groups with policy changes. Canay, Romano, and Shaikh (2017) provide approximate randomization tests.

Small sample approaches to cluster robust inference include Donald and Lang (2007), Imbens and Koles\'{a}r (2016), and Young (2016). Bootstrap approaches are provided by Cameron, Gelbach and Miller (2008), and MacKinnon and Webb (2017, 2018).

A recent contribution which develops cluster-robust inference for GMM is Hwang (2017).

The organization of the paper is as follows. After Section 2, which introduces cluster sampling, Sections 3-8 cover the core asymptotic theory, providing rigorous conditions for the WLLN (Section 3), rates of convergence (Section 4), the CLT (Section 5), cluster-robust covariance matrix estimation (Section 6), the ULLN (Section 7), and the CLT for clustered second moments (Section 8). Following this, we provide the distribution theory for the core econometric estimators, specifically linear regression and 2SLS (Section 9), Maximum Likelihood (Section 10), and GMM\ (Section 11). Each of these latter sections are written self-sufficiently, so they can be used directly by readers. Proofs of the core theorems are provided in the Appendix, and proofs for the applications are provided in the Supplemental Appendix.

Cluster Sampling

The observations are $X_{i}\in \mathbb{R} ^{p}$ for $i=1,...,n$. They are grouped into $G$ mutually independent known clusters, indexed $g=1,...,G$, where the $g^{th}$ cluster has $n_{g}$ observations. The clustering can be due to the sampling scheme, or done by the researcher due to known correlation structures. The number of observations $n_{g}$ per cluster (the \textquotedblleft cluster sizes\textquotedblright) may vary across clusters. The total number of observations are $n=\sum _{g=1}^{G}n_{g}$. It will also be convenient to double-index the observations as $X_{gj}$ for $g=1,...,G$ and $j=1,...,n_{g}$.

As is conventional in the clustering literature, the only dependence assumption we make is that the observations are independent across clusters, while the dependence within each cluster is unrestricted. Furthermore, we do not require that the observations or clusters come from identical distributions. Thus our framework includes i.n.i.d (independent, not necessarily identically distributed) as the special case $n_{g}=1$.

The notation and assumptions allow for linear panel data models with cluster-specific fixed effects. In this case the observations $X_{gj}$ should be viewed as clustered-demeaned observations. Another common application is linear panel data models with both cluster-specific and time-specific fixed effects. Our assumptions do not cover this case as removing the time effects will induce cross-cluster correlations. This is essentially \textquotedblleft multiway\textquotedblright\ clustering and requires different methods. See MacKinnon, Nielsen and Webb (2017).

Our distributional framework is asymptotic as $n$ and $G$ simultaneously diverge to infinity. This is typically referred to as the \textquotedblleft large $G$\textquotedblright\ framework. Our assumptions, however, will allow $G$ to diverge at a rate slower than $n$, by allowing the cluster sizes $n_{g}$ to diverge. This is in contrast to the early asymptotic theory for clustering, which implicitly assumed that the cluster sizes were bounded.

Our theory assumes that the clusters are known, and observations are independent across clusters. This is a substantive restriction. Alternatively, it may be possible to develop a distribution theory which allows weak dependence across clusters, but we do not do so here.

A word on notation. For a vector $a$ let $\left\Vert a\right\Vert =\left( a^{\prime}a\right) ^{1/2}$ denote the Euclidean norm. For a positive semi-definite matrix $A$ let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote its smallest and largest eigenvalue, respectively. For a general matrix $A$ let $\left\Vert A\right\Vert =\sqrt{\lambda_{\max}\left( A^{\prime }A\right) }$ denote the spectral norm. For a positive semi-definite matrix $A$ let $A^{1/2}$ denote the symmetric square root matrix such that $A^{1/2}A^{1/2}=A$. We let $C$ denote a generic positive constant, that may be different in different uses.

Weak Law of Large Numbers

For our core theory (WLLN & CLT), we focus on the sample mean $\overline {X}_{n}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ as an estimator of $E\overline{X}_{n} $. It will be convenient to define the cluster sums \[ \widetilde{X}_{g}=\sum_{j=1}^{n_{g}}X_{gj} \] which are mutually independent under clustered sampling. The sample mean can then be written as \[ \overline{X}_{n}=\frac{1}{n}\sum_{g=1}^{G}\widetilde{X}_{g}. \]

We use the following regularity condition.

The condition ((ref)) states that $X_{i}$ is uniformly integrable\footnote{A referee points out that the sup in ((ref)) could be weakened to an average. However our later results will use uniform integrability conditions similar to ((ref)) so we state all results in this format.}. This condition is identical to the standard condition for the WLLN for independent heterogeneous observations, and thus Theorem (ref) is a direct generalization of the WLLN for i.n.i.d. samples. ((ref)) simplifies to $E\left\Vert X_{i}\right\Vert <\infty$ when the observations have identical marginal distributions. A sufficient condition allowing for distributional heterogeneity is $\sup_{i}E\left\Vert X_{i}\right\Vert ^{r}<\infty$ for some $r>1$.

Assumption (ref) states that each cluster size $n_{g}$ is asymptotically negligible. This implies $G\rightarrow\infty$, so we do not explicitly need to list the latter as an assumption. Assumption (ref) allows for considerable heterogeneity in cluster sizes. It allows the cluster sizes to grow with sample size, so long as the growth is not proportional. For example, it allows clusters to grow at the rate $n_{g}=n^{\alpha}$ for $0\leq\alpha<1$.

Assumption (ref) is necessary for parameter estimation consistency while allowing arbitrary within-cluster dependence. Otherwise a single cluster could dominate the sample average. To see this, suppose that there is a cluster $\ell$ such that all observations within the cluster are identical, so that $X_{\ell j}=Z_{\ell}$ for some non-degenerate random variable $Z_{\ell}$, and that this cluster violates Assumption (ref), so that $n_{\ell}/n\rightarrow c>0$. Suppose for all other clusters that $EX_{gj}=0$ and $n_{g} /n\rightarrow0$. Then $\overline{X}_{n}\xrightarrow{p}Z_{\ell}$ and is inconsistent. Thus Assumption (ref) is necessary for the WLLN ((ref)) if we allow for unstructured cluster heterogeneity.

Assumption (ref) is equivalent to the condition

To see this, first observe that since $\sum_{g=1}^{G}n_{g}=n$, the left-hand-side of ((ref)) is smaller than $\max\limits_{g\leq G} n_{g}/n\rightarrow0$ under Assumption (ref). Thus Assumption (ref) implies ((ref)). Second, \[ \max\limits_{g\leq G}\dfrac{n_{g}}{n}=\left( \max\limits_{g\leq G} \dfrac{n_{g}^{2}}{n^{2}}\right) ^{1/2}\leq\left( \sum_{g=1}^{G}\dfrac {n_{g}^{2}}{n^{2}}\right) ^{1/2}\rightarrow0 \] under ((ref)). Thus ((ref)) implies Assumption (ref), so the two are equivalent.

Rate of Convergence

Under i.i.d. sampling the rate of convergence of the sample mean is $n^{-1/2} $. Clustering can alter the rate of convergence. In this section we explore possible rates of convergence. From the work of C. Hansen (2007) it has been understood that if the dependence within each cluster is weak then the rate of convergence would be the i.i.d. rate $n^{-1/2}$ but if the dependence within each cluster is strong then the rate of convergence would be determined by the number of clusters: $G^{-1/2}$. What we now show is that the rate of convergence can be in between or even slower than these rates.

The convergence rate can be calculated as the standard deviation of the sample mean. For simplicity we focus on the scalar case $p=1$. The standard deviation of $\overline{X}_{n}$ is \[ \text{sd}\left( \overline{X}_{n}\right) =\frac{1}{n}\left( \sum_{g=1} ^{G}\text{var}(\widetilde{X}_{g})\right) ^{1/2}. \]

We now consider several examples. For our first four we take the case where the clusters are all the same size: $n_{g}=n^{\alpha}$ for $0<\alpha<1$. In this case the number of clusters is $G=n^{1-\alpha}$.

We first consider a case where the convergence is the i.i.d. rate $n^{-1/2}$.

\noindentExample 1. The observations are independent within each cluster and $\text{var}(X_{i})=1$. Then \[ \text{var}(\widetilde{X}_{g})=n_{g}=n^{\alpha} \] and \[ \text{sd}\left( \overline{X}_{n}\right) =n^{-1/2}. \]

The $n^{-1/2}$ rate extends to any case where the within-cluster dependence is weak, including autoregressive and moving average dependence.

Our second example is a case where the convergence is determined by the number of clusters.

\noindentExample 2. The observations are identical within each cluster (e.g. perfectly correlated) and $\text{var}(X_{i})=1$. Then \[ \text{var}(\widetilde{X}_{g})=n_{g}^{2}=n^{2\alpha} \] and \[ \text{sd}\left( \overline{X}_{n}\right) =n^{-(1-\alpha)/2}=G^{-1/2}. \]

The assumption that the observations are perfectly correlated is not essential to obtain the $G^{-1/2}$ rate. What is important is that there is a common component to the observations within a cluster.

Our third example is a case where the convergence rate is in between the above two cases. Not surprisingly, it can obtained by constructing strong but decaying within-cluster dependence.

\noindentExample 3. The observations are correlated within each cluster with $\text{var}(X_{i})=1$ and $\text{cov}(X_{gj},X_{gl})=1/|j-l|$. Then \[ \text{var}(\widetilde{X}_{g})\sim n_{g}\log n_{g}\sim n^{\alpha}\log n \] and \[ \text{sd}\left( \overline{X}_{n}\right) \sim\sqrt{\log n/n}. \] Furthermore, $G\text{var}\left( \overline{X}_{n}\right) \rightarrow0.$ Thus $\text{sd}\left( \overline{X}_{n}\right) $ converges at a rate in between $n^{-1/2}$ and $G^{-1/2}$.

Our next two examples are somewhat surprising. They are cases where the convergence rate is slower than both $n^{-1/2}$ and $G^{-1/2}$.

\noindentExample 4. The observations follow random walks within each cluster: $X_{gj}=X_{gj-1}+\varepsilon_{gj}$ with $\varepsilon_{gj}$ i.i.d. $(0,1)$ and $X_{g0}=0.$ Then \[ \text{var}(\widetilde{X}_{g})\sim n_{g}^{3} \] and \[ \text{sd}\left( \overline{X}_{n}\right) \sim n^{\alpha-1/2}. \] Thus $\text{sd}\left( \overline{X}_{n}\right) $ converges at a rate slower than both $n^{-1/2}$ and $G^{-1/2}$.

\noindentExample 5. The clusters are of two sizes, $n_{g}=1$ and $n_{g}=n^{\alpha}$. There are $G_{1}=n/2$ of the first type and $G_{2} =n^{1-\alpha}/2$ of the second type. (So $G=G_{1}+G_{2}=O\left( n\right) $.) Within each cluster the observations are identical and have unit variances. var$(\widetilde{X}_{g})$ for the two types of clusters are $1$ and $n^{2\alpha}$, respectively. Then \[ \text{sd}\left( \overline{X}_{n}\right) =\left( \frac{G_{1}+G_{2} n^{2\alpha}}{n^{2}}\right) ^{1/2}=\left( \frac{1+n^{\alpha}}{2n}\right) ^{1/2}=O\left( n^{-(1-\alpha)/2}\right) . \] Thus $\text{sd}\left( \overline{X}_{n}\right) $ converges at at a rate slower than both $n^{-1/2}$ and $G^{-1/2}$.

The final example illustrates the importance of considering heterogeneous cluster sizes. The reason why the convergence rate is slower than both $n^{-1/2}$ and $G^{-1/2}$ is because the number of clusters is determined by the large number of small clusters, but the convergence rate is determined by the (relatively) small number of large clusters.

What we have seen is that the convergence rate $\text{sd}\left( \overline {X}_{n}\right) $ can equal the square root of sample size $n^{-1/2}$, can equal the square root of the number of groups $G^{-1/2}$, can be in between $G^{-1/2}$ and $n^{-1/2}$, or can be slower than both $n^{-1/2}$ and $G^{-1/2}$.

When $\overline{X}_{n}$ is a vector, it is likely that its elements converge at different rates since they can have different within-cluster correlation structures. For example, some variables could be independent within clusters while others could be identical within clusters.

These examples show that under cluster dependence the convergence rate is context-dependent and variable-dependent, and it is therefore important to allow for general rates of convergence and to not impose arbitrary rates in asymptotic analysis.

Central Limit Theory

Under i.i.d. sampling the standard deviation of the sample mean is of order $O(n^{-1/2})$, so $\sqrt{n}$ is the appropriate scaling to obtain the central limit theorem (CLT). As discussed in the previous section, clustering can alter the rate of convergence, so it is essential to standardize the sample mean by the actual variance rather than an assumed rate. The variance matrix of $\sqrt{n}\overline{X}_{n}$ is

We use the following regularity condition.

Theorem (ref) provides a CLT for cluster samples which generalizes the classic CLT for independent heterogeneous samples. The latter holds with $r=2 $, $n_{g}=1$ and $G=n$.

Assumption (ref) and (ref) are stronger than Assumption (ref) and (ref), and thus the conditions for the CLT imply those for the WLLN.

The condition ((ref)) states that $\left\Vert X_{i}\right\Vert ^{r}$ is uniformly integrable. When $r=2$ this is similar to the Lindeberg condition for the CLT under independent heterogeneous sampling. ((ref)) simplifies to $E\left\Vert X_{i}\right\Vert ^{r}<\infty$ when the observations have identical marginal distributions. A sufficient condition allowing for distributional heterogeneity is $\sup_{i}E\left\Vert X_{i}\right\Vert ^{s}<\infty$ for some $s>r\geq2$.

Assumption (ref) ((ref)) is a restriction on the cluster sizes. It involves a trade-off with the number of moments $r$. It is least restrictive for large $r$, and more restrictive for small $r$. As $r\rightarrow\infty$ it approaches $\max_{g\leq G}n_{g}^{2}/n = O(1)$, which is implied by Assumption (ref) ((ref)).

Assumption (ref) allows for growing and heterogeneous cluster sizes. For example, it allows clusters to grow uniformly at the rate $n_{g}=n^{\alpha}$ for $0\leq\alpha\leq(r-2)/2(r-1)$. (Note that this requires the cluster sizes to be bounded if $r=2$.) It also allows for only a small number of clusters to grow. For example, suppose that $n_{g}=\overline{n}$ (bounded) for $G-K$ clusters and $n_{g}=G^{\alpha/2 }$ for $K$ clusters, with $K$ fixed. Then Assumption (ref) holds for any $\alpha<1$ and $r\geq2$.

Assumption (ref) ((ref)) is implied by

and they are equivalent when the cluster sizes are homogeneous. In general, however, ((ref)) is less restrictive than ((ref)). For example, when $r=2$, ((ref)) requires the cluster sizes to be bounded, while ((ref)) does not. (Consider the heterogeneous example given in the previous paragraph. This satisfies ((ref)) but not ((ref)) when $r=2$.)

The condition ((ref)) specifies that $\text{var}\left( \sqrt{n} \alpha^{\prime}\overline{X}_{n}\right) $ does not vanish for any conformable vector $\alpha\neq0$. This excludes degenerate cases and perfect negative within-cluster correlation. In general, if $X_{i}$ is non-degenerate then ((ref)) is not restrictive as there is no reasonable setting where it will be violated. If $\overline{X}_{n}$ converges at rate $n^{-1/2}$ then $\lambda_{n}=O(1)$ but when $\overline{X}_{n}$ converges at rate slower than $n^{-1/2}$ then $\lambda_{n}$ will actually diverge with $n$. It should also be mentioned that condition ((ref)) allows the components of $\Omega_{n}$ to converge at different rates.

Our proof of Theorem (ref) actually uses the conditions

and

instead of ((ref))-((ref)). ((ref))-((ref)) is weaker than ((ref))-((ref)) when $\lambda_{n}$ diverges to infinity (which occurs when $\overline{X}_{n}$ converges at a rate slower than $n^{-1/2}$). Since the sequence $\lambda_{n}$ is unknown in an application it is difficult to interpret the assumptions ((ref))-((ref)). Hence we prefer the assumptions ((ref))-((ref)).

The conditions ((ref))-((ref)) may be stronger than necessary when within-cluster dependence is weak, but are necessary under strong within-cluster dependence. To see this, suppose that all observations within a cluster are identical, so that $X_{gj}=Z_{g}$ and $Z_{g}$ has a finite variance but no higher moments. Then the Lindeberg condition for the CLT can be simplified to \[ \sum_{g=1}^{G}\frac{n_{g}^{2}}{n\lambda_{n}}E\left( \left\Vert Z_{g} \right\Vert ^{2}1\left( \left\Vert Z_{g}\right\Vert ^{2}\geq\frac {n\lambda_{n}\varepsilon}{n_{g}^{2}}\right) \right) \rightarrow0 \] for all $\varepsilon>0$. Each term in the sum must limit to zero, which requires ((ref))-((ref)) with $r=2$.

We now compare our conditions with those of Djogbenou, MacKinnon, and Nielsen (2018). Their Assumption 3 states (in our notation) for $r\geq4$

Equation ((ref)) implies and is stronger than ((ref)). Calculations similar to those in our appendix show that $\lambda_{n}\leq O\left( \max _{g}n_{g}\right) =O(n)$. So ((ref)) also implies \[ \left( \max\limits_{g\leq G}\dfrac{n_{g}^{2}}{n\lambda_{n}}\right) ^{1/2}=\max\limits_{g\leq G}\dfrac{n_{g}}{n^{(r-2)/2(r-1)}\lambda _{n}^{r/2(r-1)}}\left( \frac{\lambda_{n}}{n}\right) ^{1/2(r-1)}=o\left( 1\right) \] which is ((ref)). Thus our conditions ((ref))-((ref)) are less restrictive than their condition ((ref)), and do not require $r\geq4$.

Cluster-Robust Variance Matrix Estimation

We now discuss cluster-robust covariance matrix estimation.

We first consider the case where $X_{i}$ is mean zero (or equivalently that the mean is known). In this case the covariance matrix equals \[ \Omega_{n}=\frac{1}{n}\sum_{g=1}^{G}E\left( \widetilde{X}_{g}\widetilde{X} _{g}^{\prime}\right) . \] In this case a natural estimator is \[ \widetilde{\Omega}_{n}=\frac{1}{n}\sum_{g=1}^{G}\widetilde{X}_{g} \widetilde{X}_{g}^{\prime}. \]

Theorem (ref) shows that the cluster-robust covariance matrix estimator is consistent, and replacing the covariance matrix in the CLT\ with the estimated covariance matrix does not affect the asymptotic distribution. Implications of ((ref)) are that cluster-robust t-ratios are asymptotically standard normal, and that cluster-robust Wald statistics are asymptotically chi-square distributed with $p$ degrees of freedom.

Construction of practical covariance matrix estimators is context-specific, depending on the mean structure. For example, suppose that $\mu=EX_{i}$ does not vary across observations. In this case we can write \[ \Omega_{n}=\frac{1}{n}\sum_{g=1}^{G}E\left( \widetilde{X}_{g}\widetilde{X} _{g}^{\prime}\right) -\frac{1}{n}\sum_{g=1}^{G}n_{g}^{2}\mu\mu^{\prime}. \] The natural estimator for $\mu$ is $\overline{X}_{n}$ and that for $\Omega _{n}$ is \[ \widehat{\Omega}_{n}=\frac{1}{n}\sum_{g=1}^{G}\widetilde{X}_{g}\widetilde{X} _{g}^{\prime}-\frac{1}{n}\sum_{g=1}^{G}n_{g}^{2}\overline{X}_{n}\overline {X}_{n}^{\prime}. \]

Uniform Laws of Large Numbers

Now consider a uniform WLLN. Consider functions $f(x,\theta)\in\mathbb{R}^{k}$ indexed on $\theta\in\Theta$ where $\Theta$ is compact. Define the sample mean \[ \overline{f}_{n}(\theta)=\frac{1}{n}\sum_{i=1}^{n}f(X_{i},\theta). \]

The following result is an application of Theorem 3 of Andrews (1992).

We also consider a uniform law for the clustered variance. Set $\mu (\theta)=Ef(X_{i},\theta)$ so that it does not vary across observations. The variance of $\sqrt{n}\overline{f}_{n}(\theta)$ is

where $\widetilde{f}_{g}(\theta)=\sum_{j=1}^{n_{g}}f(X_{gj},\theta)$ are the cluster sums. An appropriate estimator for $\Omega_{n}(\theta)$ is \[ \widehat{\Omega}_{n}(\theta)=\frac{1}{n}\sum_{g=1}^{G}\widetilde{f}_{g} (\theta)\widetilde{f}_{g}(\theta)-\frac{1}{n}\sum_{g=1}^{G}n_{g}^{2} \overline{f}_{n}(\theta)\overline{f}_{n}(\theta)^{\prime}. \] In practice, a simpler estimator \[ \widetilde{\Omega}_{n}(\theta)=\frac{1}{n}\sum_{g=1}^{G}\widetilde{f} _{g}(\theta)\widetilde{f}_{g}(\theta)^{\prime} \] is often used if $\mu(\theta_{0})=0$ for $\theta_{0}\in interior\left( \Theta\right) $ and $\widehat{\theta}\xrightarrow{p}\theta_{0}$ for some estimator $\widehat{\theta}$.

The following result is an extension of Theorem (ref) to the case of clustered variance estimators. It also relies on Theorem 3 of Andrews (1992).

Central Limit Theorem for Clustered Second Moments

Although our primary focus is the sample mean, the core theory can be extended to statistics which are not sample means. In this section, we focus on the vectorized variance estimators \[ \overline{f}_{G}=\frac{1}{n}\sum_{g=1}^{G}\widetilde{f}_{g} \] where \[ \widetilde{f}_{g}=\widetilde{X}_{g}\otimes\widetilde{X}_{g} \] or \[ \widetilde{f}_{g}=\left( \widetilde{X}_{g}-n_{g}\overline{X}_{n}\right) \otimes\left( \widetilde{X}_{g}-n_{g}\overline{X}_{n}\right) . \] The WLLN for $\overline{f}_{G}$ holds by Theorem (ref) ((ref)) and Theorem (ref) ((ref)), and the ULLN for $\overline{f}_{G}$ holds by Theorem (ref). However, the CLT given in Theorem (ref) cannot be applied to $\overline{f}_{G}$ because $\overline{f}_{G}$ cannot be written as the sample mean over $i$. We provide the CLT for $\overline{f}_{G}$ below. This is useful to establish asymptotic distributions of estimators in a non-standard setting. For example, the asymptotic distribution of the generalized method of moments (GMM) estimators depends on the limiting distribution of the weight matrix when the moment condition is misspecified (Hall and Inoue, 2003; Lee, 2014; Hansen and Lee, 2018).

Similar to the sample mean, the convergence rate of $\overline{f}_{G}$ can vary under cluster dependence. Consider $\widetilde{f}_{g} = \widetilde{X} _{g}\otimes\widetilde{X}_{g}$ and assume $p=1$ for simplicity. The standard deviation of $\overline{f}_{G}$ is \[ \text{sd}\left( \overline{f}_{G}\right) =\frac{1}{n}\left( \sum_{g=1} ^{G}\text{var}\left( \widetilde{X}_{g}\widetilde{X}_{g}\right) \right) ^{1/2} = \frac{1}{n}\left( \sum_{g=1}^{G}\sum_{j=1}^{n_{g}}\sum_{l=1}^{n_{g} }\text{var}\left( X_{gj}X_{gl}\right) \right) ^{1/2}. \] Under i.i.d. sampling $\text{sd}\left( \overline{f}_{G}\right) =O\left( n^{-1/2}\right) $. Under the Examples 1 and 2 in Section (ref), the convergence rate is $G^{-1/2}$.

Define the variance matrix of $\sqrt{n}\overline{f}_{G}$ as

We use the following regularity condition.

Note that Assumption (ref) is a strengthening of Assumption (ref).

Finally we provide a CLT combining the previous results. For $Y_{i} \in\mathbb{R}^{s}$, $i=1,...,n$, obtained by cluster sampling, let $\widetilde{\psi}_{g}$ be the stacked vector \[ \widetilde{\psi}_{g}=\left(

\right) \] or \[ \widetilde{\psi}_{g}=\left(

\right) \] and $\overline{\psi}_{G}=n^{-1}\sum_{g=1}^{G}\widetilde{\psi}_{g}$. Let the variance matrix of $\sqrt{n}\overline{\psi}_{G}$ be \[ \Omega_{n}=E\left( n\left( \overline{\psi}_{G}-E\overline{\psi}_{G}\right) \left( \overline{\psi}_{G}-E\overline{\psi}_{G}\right) ^{\prime}\right) . \] The following Corollary provides the CLT for the joint process. Since it immediately follows from Theorems (ref) and (ref), the proof is omitted.

Linear Regression and Two-Stage Least Squares

It is useful to use cluster-level notation. Let $\boldsymbol{y}_{g} =(y_{g1},...,y_{gn_{g}})^{\prime}$, $\boldsymbol{X}_{g}=(\boldsymbol{x} _{g1},...,\boldsymbol{x}_{gn_{g}})^{\prime}$ and $\boldsymbol{Z} _{g}=(\boldsymbol{z}_{g1},...,\boldsymbol{z}_{gn_{g}})^{\prime}$ denote an $n_{g}\times1$ vector of dependent variables, $n_{g}\times k$ matrix of regressors, and $n_{g}\times l$ matrix of instruments for the $g^{th}$ cluster. A linear model can be written using cluster notation as

where $\boldsymbol{e}_{g}$ is a $n_{g}\times1$ error vector. The case of linear regression holds as the special case where $\boldsymbol{Z} _{g}=\boldsymbol{X}_{g}$ and $l=k$ (so that (ref) becomes identity). Assume $l\geq k$. ((ref)) is the structural equation and ((ref)) is the first-stage equation.

The two-stage least squares (2SLS) estimator for $\boldsymbol{\beta}$ can be written as \[ \widehat{\boldsymbol{\beta}}=\left( \sum_{g=1}^{G}\boldsymbol{X}_{g}^{\prime }\boldsymbol{Z}_{g}\left( \sum_{g=1}^{G}\boldsymbol{Z}_{g}^{\prime }\boldsymbol{Z}_{g}\right) ^{-1}\sum_{g=1}^{G}\boldsymbol{Z}_{g}^{\prime }\boldsymbol{X}_{g}\right) ^{-1}\left( \sum_{g=1}^{G}\boldsymbol{X} _{g}^{\prime}\boldsymbol{Z}_{g}\left( \sum_{g=1}^{G}\boldsymbol{Z} _{g}^{\prime}\boldsymbol{Z}_{g}\right) ^{-1}\sum_{g=1}^{G}\boldsymbol{Z} _{g}^{\prime}\boldsymbol{y}_{g}\right) . \] We first show consistency of $\widehat{\boldsymbol{\beta}}$. Define

Next we provide the asymptotic distribution. Define

The residuals for the $g^{th}$ cluster are \[ \widehat{\boldsymbol{e}}_{g}=\boldsymbol{y}_{g}-\boldsymbol{X}_{g}\widehat{ \boldsymbol{\beta}}. \]

Define

The variance estimator is \[ \widehat{V}_{n}=d_{n}\left( \widehat{Q}_{n}^{\prime}\widehat{W}_{n} ^{-1}\widehat{Q}_{n}\right) ^{-1}\widehat{Q}_{n}^{\prime}\widehat{W}_{n} ^{-1}\widehat{\Omega}_{n}\widehat{W}_{n}^{-1}\widehat{Q}_{n}\left( \widehat{Q}_{n}^{\prime}\widehat{W}_{n}^{-1}\widehat{Q}_{n}\right) ^{-1}. \] with $d_{n}$ a possible finite-sample degree-of-freedom adjustment. For example, C. Hansen (2007) proposed $d_{n}=G/(G-1)$ for the regression case (under homogeneous cluster sizes), and Stata sets \[ d_{n}=\left( \frac{n-1}{n-k}\right) \left( \frac{G}{G-1}\right) \] for the OLS and 2SLS estimators under cluster option.

The standard errors for $R_{n}^{\prime}\widehat{\boldsymbol{\beta}}$ can be obtained by taking the square roots of the diagonal elements of $n^{-1} R_{n}^{\prime}\widehat{V}_{n}R_{n}$.

(Pseudo) Maximum Likelihood

Suppose that we observe a sequence of random vectors $X_{i}\in\mathbb{R}^{p}$, $i=1,...,n$ with the same marginal distributions from a density $f(x,\boldsymbol{\theta})$ where $\boldsymbol{\theta}\in\mathbb{R}^{k}$. Let $\boldsymbol{X}_{g}=(X_{g1},...,X_{gn_{g}})^{\prime}$ be a $n_{g}\times p$ matrix for each cluster. For the observations in the cluster $g$, let $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}_{0})$ be the joint density. Since the observations within the same cluster need not be independent, $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}_{0})\neq\prod_{i=1}^{n_{g} }f(X_{gi},\boldsymbol{\theta}_{0})$ in general. This also implies that $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}_{0})\neq f_{h}(\boldsymbol{X} _{h},\boldsymbol{\theta}_{0})$ for $g\neq h$. Given specification of $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}_{0})$, the maximum likelihood estimator (MLE) can be obtained as the maximizer of \[ \sum_{g=1}^{G}\log f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}). \] However, the joint density $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta})$ may be difficult to specify in practice. A simpler alternative is to use a pseudo-likelihood $\prod_{i=1}^{n_{g}}f(X_{gi},\boldsymbol{\theta}_{0})$ for the joint density $f_{g}(\boldsymbol{X}_{g},\boldsymbol{\theta}_{0})$, and specify the log likelihood function as \[ L_{n}(\theta)=\sum_{g=1}^{G}\sum_{j=1}^{n_{g}}\log f(X_{gj},\boldsymbol{\theta }). \] Define the pseudo-MLE as \[ \widehat{\boldsymbol{\theta}}=\operatorname*{arg\,max}_{\theta\in\Theta}L_{n} (\boldsymbol{\theta}). \] This estimator is also called the partial (or pooled) MLE (Wooldridge, 2010).

This estimator is the standard implementation of MLE under clustered dependence. To our knowledge there is no existing distribution theory for this standard estimator.

We first show consistency of $\widehat{\boldsymbol{\theta}}$. The following is based on Theorem 2.1 of Newey and McFadden (1994).

Next we show the asymptotic distribution. Define

Define the sample versions

The variance estimator is \[ \widehat{V}_{n} = \widehat{H}_{n}(\widehat{\boldsymbol{\theta}})^{-1} \widehat{\Omega}_{n}(\widehat{\boldsymbol{\theta}})\widehat{H}_{n} (\widehat{\boldsymbol{\theta}})^{-1}. \] Note that the information matrix equality does not hold because $\sum _{j=1}^{n_{g}}\log f(X_{gj},\boldsymbol{\theta}_{0})\neq f_{g}(\boldsymbol{X} _{g},\boldsymbol{\theta}_{0})$ in general.

The standard errors for $R_{n}^{\prime}\widehat{\boldsymbol{\beta}}$ can be obtained by taking the square roots of the diagonal elements of $n^{-1} R_{n}^{\prime}\widehat{V}_{n}R_{n}$.

Generalized Method of Moments

Suppose that we observe a sequence of random vectors $X_{i}\in\mathbb{R}^{p}$, $i=1,...,n$ from cluster sampling. A known moment function is given by $m(X_{i},\boldsymbol{\theta})$ where $m(\cdot,\cdot)$ is $l\times1$ and $\boldsymbol{\theta}$ is $k\times1$. Define the cluster sum as \[ \widetilde{m}_{g}(\boldsymbol{\theta})=\sum_{j=1}^{n_{g}}m(X_{gj} ,\boldsymbol{\theta}). \] An unconditional moment model in cluster notation is given by

We assume that $\boldsymbol{\theta}_{0}$ is identified and $l>k$ so the moment model is over-identified. Write the sample mean of the moment function as \[ \overline{m}_{n}(\boldsymbol{\theta})=\frac{1}{n}\sum_{i=1}^{n}m(X_{i} ,\boldsymbol{\theta}). \] Since (ref) holds for all $g=1,...,G$, the usual unconditional moment condition $E\overline{m}_{n}(\boldsymbol{\theta}_{0})=0$ follows. The generalized method of moments (GMM) estimator is given by

where $\widehat{W}_{n}^{-1}$ is an $l\times l$ positive definite weight matrix, which may or may not depend on an estimated parameter. Typically, the weight matrix is obtained by plugging in a preliminary consistent estimator, $\widetilde{\boldsymbol{\theta}}$, so that $\widehat{W}_{n}^{-1} =\widehat{W}_{n}(\widetilde{\boldsymbol{\theta}})^{-1}$.

We consider two forms of GMM estimator. The first one is based on a non-clustered weight matrix, which takes the form of

for some $l\times1$ vector $v(x,\boldsymbol{\theta})$. This includes the conventional one-step and two-step GMM estimators. For 2SLS, $v(X_{i} ,\boldsymbol{\theta})=Z_{i}$ where $Z_{i}$ is an $l\times1$ vector of instruments. The efficient two-step GMM uses $v(X_{i},\boldsymbol{\theta })=m(X_{i},\boldsymbol{\theta})$ or $v(X_{i},\boldsymbol{\theta} )=m(X_{i},\boldsymbol{\theta})-\overline{m}_{n}(\boldsymbol{\theta})$. The conventional efficient weight matrix, however, does not provide efficiency anymore under cluster sampling because a weight matrix of the form of (ref) is not consistent for the variance matrix of $\sqrt{n} (\overline{m}_{n}(\boldsymbol{\theta})-E\overline{m}_{n}(\boldsymbol{\theta }))$ in general.

The second is based on the clustered efficient weight matrix, which leads to the two-step efficient GMM under cluster sampling. The weight matrix takes the form of

Alternatively, the uncentered version of $\widehat{W}_{n}(\boldsymbol{\theta })$ and $\widehat{\Omega}_{n}(\boldsymbol{\theta})$ can be used to obtain the efficient two-step GMM estimator but the centered version is generally recommended. For more discussion, see Hansen (2018).

Since we assume that the weight matrix depends on a consistent preliminary estimator, we exclude the continuously updating (CU) GMM estimator in our analysis. Whenever possible, we omit the dependence of the weight matrices on $\widetilde{\boldsymbol{\theta}}$ and write $\widehat{W}_{n} = \widehat{W} _{n}(\widetilde{\boldsymbol{\theta}})$. Define $W_{n} = E\widehat{W} _{n}(\boldsymbol{\theta}_{0})$.

We first show consistency of the GMM estimator. The following is based on Theorem 2.1 of Newey and McFadden (1994).

Primitive conditions under which Condition 6 of Theorem (ref) holds can be found given the choice of the weight matrix. For simplicity, we assume that if the conventional weight matrix is used then either $v(X_{i},\boldsymbol{ \theta}) =m(X_{i},\boldsymbol{\theta})$ or $v(X_{i},\boldsymbol{\theta}) = m(X_{i},\boldsymbol{\theta})-\overline{m}_{n}(\boldsymbol{\theta})$. If the clustered weight matrix is used then it takes the form of (ref). The conditions of Theorem (ref) are sufficient for Condition 6 of Theorem (ref) to hold.

To show the asymptotic distribution of the GMM estimator, define

where $Q_{n}=Q_{n}(\boldsymbol{\theta}_{0})$ and $\Omega_{n}=\Omega _{n}(\boldsymbol{\theta}_{0})$. If the clustered efficient weight matrix (ref) is used, then the asymptotic variance matrix simplifies to \[ V_{n}=(Q_{n}^{\prime}\Omega_{n}^{-1}Q_{n})^{-1}. \] Define the sample versions as

and let $\widehat{Q}_{n}=\widehat{Q}_{n}(\widehat{\boldsymbol{\theta}})$ and $\widehat{\Omega}_{n}=\widehat{\Omega}_{n}(\widehat{\boldsymbol{\theta}})$. The variance estimator is \[ \widehat{V}_{n}=(\widehat{Q}_{n}^{\prime}\widehat{W}_{n}^{-1}\widehat{Q} _{n})^{-1}\widehat{Q}_{n}^{\prime}\widehat{W}_{n}^{-1}\widehat{\Omega} _{n}\widehat{W}_{n}^{-1}\widehat{Q}_{n}(\widehat{Q}_{n}^{\prime} \widehat{W}_{n}^{-1}\widehat{Q}_{n})^{-1}, \] if $\widehat{W}_{n}$ is given by (ref) and \[ \widehat{V}_{n}=(\widehat{Q}_{n}^{\prime}\widehat{\Omega}_{n}^{-1} \widehat{Q}_{n})^{-1}, \] if $\widehat{W}_{n}$ is given by (ref), i.e., $\widehat{W} _{n}=\widehat{\Omega}_{n}$.

The over-identifying restrictions test (the J test, hereinafter) is a test based on the GMM criterion to test whether the moment model is correctly specified or not, i.e., $E\widetilde{m}_{g}(\boldsymbol{\theta}_{0})=0$. An implication of cluster sampling is that the conventional J test statistic will not have a standard chi-square asymptotic distribution because the conventional efficient weight matrix is not consistent for the inverse of the variance matrix of the moment function. The GMM criterion (ref) based on the clustered efficient weight matrix (ref) evaluated at the estimator is the robust J test statistic. Define \[ J_{n}(\widehat{\boldsymbol{\theta}})=n\cdot\overline{m}_{n} (\widehat{\boldsymbol{\theta}})^{\prime}\widehat{W}_{n}^{-1}\overline{m} _{n}(\widehat{\boldsymbol{\theta}}). \]

The standard errors for $R_{n}^{\prime}\widehat{\boldsymbol{\beta}}$ can be obtained by taking the square roots of the diagonal elements of $n^{-1} R_{n}^{\prime}\widehat{V}_{n}R_{n}$.