Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
40,338 characters · 5 sections · 60 citation commands
On Using The Two-Way Cluster-Robust Standard Errors
Multi-way clustering is ubiquitous in empirical studies. For example, market structures by construction induce two-way clustering, where common supply shocks cause cluster dependence within a firm across markets and common demand shocks cause cluster dependence within a market across firms. To account for such forms of cluster dependence, researchers often use the two-way cluster-robust (TWCR) standard errors proposed by CGM2011 and thompson2011simple.
A key to the inference based on the TWCR standard errors of CGM2011 and thompson2011simple is the asymptotic gaussianity. However, the recent econometrics literature menzel2021bootstrap has pointed out the potential non-gaussianity of the limit distribution. In this light, this literature proposes alternative inference procedures that do not rely on asymptotic gaussianity.
With this said, a large number of empirical papers have already reported their TWCR standard errors. Specifically, CGM2011 and thompson2011simple have attracted 3,500 citations and 1,600 citations, respectively,\footnote{We obtained these numbers from Google Scholar in December 2022.} where most of these papers are empirical research papers that actually use their TWCR standard errors. If the non-gaussianity were indeed a common feature, then the whole body of this empirical economics literature would require re-investigation.
A natural question is, therefore, whether the asymptotic gaussianity of the statistics commonly occurs under two-way clustering. Some simulation studies will easily convince us that it is fairly common, and non-gaussianity is rather exceptional. This observation is encouraging and supports the common practice of two-way cluster-robust inference based on the TWCR standard errors.
In this paper, we provide a theoretical justification for this encouraging observation. By taking the asymptotics through the lens of triangular arrays, we show that the gaussian limit distribution is the norm rather than an outlier. Specifically, we show that two-way clustered triangular arrays are guaranteed to have asymptotically gaussian limiting distributions under very mild conditions. What these conditions concern shares a natural resemblance to the number of factors in factor models, and these conditions are therefore easily interpretable from the viewpoint of economic models. Concretely, our theory suggests that a researcher should worry about the potential non-gaussianity only in those peculiar situations in which the data-generating model takes the form of a sum of a very small number of mean-zero two-way interactive factors. In other words, a researcher can in fact enjoy the TWCR standard errors in most situations. We, therefore, hope that this paper provides a theoretical justification for the legitimacy of most, if not all, of the thousands of those empirical papers that use the TWCR standard errors of CGM2011 and thompson2011simple and shed some lights on the empirical practice of TWCR for future empirical papers to come as well.
{\bf Relation to the Literature:} Our way of viewing the data generating processes and asymptotics in terms of triangular arrays and exploiting the asymptotically gaussian degenerate one-sample $U$-statistic structure is closely related to the literature of specification testing hong1995consistent,fan1996consistent,kankanala2022kernel, many weak instruments andrews2007testing,newey2009generalized, small-bandwidth asymptotics cattaneo2014small, regression discontinuity designs porter2015regression, network formation models graham2017econometric, many regressors cattaneo2018alternative, and algorithmic subsampling lee2022least, to list but a few. Notably, non-gaussianity can be safely ruled out in the corresponding degenerate $U$-statistics in all these applications. Unlike these existing papers, however, statistics based on two-way clustered triangular arrays do not take the form of a one-sample degenerate $U$-statistic structure -- it rather takes the form of a two-sample $U$-statistics with unknown kernel and unobserved underlying random variables. This key difference renders the proof strategies in the existing literature inapplicable, and thus motivates us to develop our own new theoretical results. In terms of the setup, our paper is built upon the literature that utilizes exchangeable models for network or two-way dependence considered in bickel2011method, menzel2021bootstrap and DDG2019, to list a few. Other alternative models for asymptotics under this type of dependence structures exist and are studied in, for example, tabord2019inference and verdier2020estimation.
Our main result, a central limit theorem (CLT) for means of two-way clustered triangular arrays, is related to those CLT results for various degenerate $U$-statistics, such as hall1984central, de1987central, and eubank1999central that are based on martingale structures. However, due to the two-way clustering, their martingale construction does not work in our setting. It is also related to the CLT in khashimov1989limit, which is a special case of our CLT, albeit it is shown to rely on a non-martingale-based proof strategy. As the Hoeffding-type decomposition of the two-way clustered mean contains extra components in comparison with the degenerate two-sample $U$-statistics, it remains unclear whether the proof strategy of Khashimov, which relies on approximating characteristic functions directly, can be readily adapted to cover our case. Therefore, we instead take a martingale-based approach to derive our own CLT.
We first provide an informal overview of the practical implications of our main result in Section (ref).
Two-way clustered data $\{D_{it}: 1 \le i \le N, 1 \le t \le T\}$ are generated by $i$-specific factors, $j$-specific factors, and idiosyncratic components. To fix ideas, consider the simple yet generic data generating process (DGP)
where $\{\alpha_{ij}\}_{j=0}^J$ are $i$-specific latent factors, $\{\gamma_{tj}\}_{j=0}^J$ are $t$-specific latent factors, and $\varepsilon_{it}$ is an idiosyncratic component. The reason that this DGP is highly representative will be made clear in Section (ref). Suppose that the factor loadings $\lambda_j$ are non-zero, and $\alpha_{ij}$ and $\gamma_{tj}$ are zero-mean non-degenerate factors for $j \in 1,\cdots,J$. In this simple setup (ref), Table (ref) summarizes the cases in which a researcher can and cannot use the TWCR standard error for $\widehat\theta = (NT)^{-1} \sum_{i=1}^N\sum_{t=1}^T D_{it}$.
This table suggests that a researcher may use the TWCR standard error in most of cases. The only pathetic situation is when the data generating model is too simple in the sense that the number $J$ of interacting economic factors is small and all of additive latent factors, $\alpha_{i0}$, $\gamma_{t0}$, and $\varepsilon_{it}$, are degenerate as in the first row in Table (ref). Section (ref) presents a formal theoretical justification for this practical guidance. Simulation studies in Section (ref) illustrate how large $J$ should be in practice.
For each $N,T\in \mathbb N$, let
where $(\alpha_i)_i$, $(\gamma_t)_t$, and $(\varepsilon_{it})_{it}$ are mutually independent i.i.d. latent Borel-random variables, and $f_{NT}$ is a real-valued Borel-measurable function. The existence of this nonlinear factor-type structure is implied by a symmetry condition known as "separate exchangeability," see the discussions in menzel2021bootstrap and DDG2019. For ease of writing and without loss of generality, we normalize the location to $E[D_{it}]=0$. Define
Our goal is to show the asymptotic gaussianity of $\widehat \theta_{NT}$.
Note that the Hoeffding-type decomposition yields
One can easily verify the following properties.
The $W_{NT}$ component in the decomposition (ref) is the potentially non-gaussian part. Specifically, it is a completely degenerate two-sample $U$-statistic, and has a non-gaussian limit distribution if $f_{NT}$ is fixed over $N,T$ -- see menzel2021bootstrap. We are going to argue that even this potentially non-gaussian component can be, and often is, asymptotically gaussian if we treat the data generating process as a triangular array. Hence, the whole $\widehat\theta_{NT}$ in (ref) is gaussian as well in this framework under some mild extra conditions.
We will write $a_{NT}(\alpha_i)=a_i$, $b_{NT}(\gamma_t)=b_t$, and $w_{NT}(\alpha_i,\gamma_t)=w_{it}$, when we want to emphasize the fact that they are transformations of the underlying latent random variables. The following lemma follows directly from eagleson1979orthogonal.
By the orthonormality of the system,
hold for all $j=1,2,\cdots$.
Lemma (ref) and Equation ((ref)) together imply that the DGP introduced in Equation ((ref)) in Section (ref) for an informal overview is indeed highly representative.
We impose the following conditions. Let us follow the mathematical convention $0/0=0$.
Assumption (ref) (i) is standard in the literature that uses exchangeable arrays to model two-way clustering -- see menzel2021bootstrap and DDG2019, for example, for more discussion on exchangeable models. Assumption (ref) (ii) states that the $U$-statistic of interest is degenerate and has at least two moments. In the setting of two-way clustering, this condition does not impose any restriction, but comes naturally from the property of the Hoeffding decomposition ((ref)). Assumption (ref) (iii) imposes the same growth rate between $N$ and $T$. This condition can be relaxed at the cost of more complicated moment restrictions in the conditional moments of $w_{NT}$. Part (iii) further imposes a standard Lyapunov-type condition ((ref)), as well as ((ref)), a condition that is analogous to the second half of Condition (2.1) in hall1984central. This Hall-type condition restricts the size of the second moment of the conditional cross-products, $E[w_{NT}(\alpha_1,\gamma_1)w_{NT}(\alpha_2,\gamma_1)|\alpha_1,\alpha_2]$ and $E[w_{NT}(\alpha_1,\gamma_1)w_{NT}(\alpha_1,\gamma_2)|\gamma_1,\gamma_2]$, relative to the size of the second moment of $w_{NT}$. See Remark (ref) below for a detailed discussion on its implications.
The first condition in Assumption (ref) imposes a standard Lyapunov condition on the leading terms in the Hoeffding decomposition ((ref)). The second condition in Assumption (ref) is a mild restriction on how the linear term and quadratic terms in the Hoeffding decomposition ((ref)) can correlate in second moments, relatively to the variances of each component. Note that, by construction, the linear and quadratic components are uncorrelated. This condition is analogous to condition (1.6) in eubank1999central.
This theorem implies that even the potentially non-gaussian term $W_{NT}$ alone in the decomposition (ref) can be gaussian. The following example illustrates a case in point.
We now turn to the asymptotic gaussianity of the whole $\widehat \theta_{NT}$. From ((ref)), observe that only the $R_{NT}$ component remains random conditionally on all $\alpha_i$'s and $\gamma_t$'s, and is conditionally independent over $(i,t)$. An application of the Lyapunov CLT therefore implies that $R_{NT}$ is conditionally asymptotically gaussian given $\alpha_i$'s and $\gamma_t$'s. Thus, by the asymptotic gaussianity of $L_{NT}+W_{NT} $ from Theorem (ref), applying Theorem 1 in chen2007asymptotic yields the asymptotic gaussianity of $\widehat\theta_{NT}$. This conclusion is summarized as a corollary below.
We focus on the data generating processes that correspond to the potentially non-gaussian components in menzel2021bootstrap, i.e., the completely degenerate two-sample U-statistics, and their variants. Such data generating processes are derived from the orthonormal representation in Section (ref). Specifically, we generate a sample $\{D_{it}: 1 \le i \le N, 1 \le t \le T \}$ by
where $\alpha_{ij} \sim N(0,1)$, $\gamma_{tj} \sim N(0,1)$, and $\varepsilon_{it} \sim N(0,1)$ are mutually independent and independent across $i \in \{1,\cdots,N\}$, $t \in \{1,\cdots,T\}$ and $j \in \{1,\cdots,J\}$. Note that this DGP is a special case of the generic DGP of Equation ((ref)) with $\lambda_{j}=J^{-1}$. We vary the data generating parameters $\delta$, $J$, and $\phi$ across sets of Monte Carlo simulations. Each set of simulations consists of 10,000 random draws of $\{D_{it}: 1 \le i \le N, 1 \le t \le T \}$. The sample size is set to $N=T=50$ throughout.
For the population mean $ \theta = E[D_{it}] $ as the parameter of interest, consider the two-way sample mean estimator $ \widehat\theta_{NT} = (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T D_{it}. $ Note that our theory from Section (ref) predicts that the gaussian approximation is asymptotically reasonable unless both $|\delta|$ and $J$ are too small. For each draw, we construct the 95% confidence intervals by two existing methods: one is based on the widely used two-way cluster-robust standard error CGM2011,thompson2011simple which presumes the asymptotic gaussianity; and the other is the recent bootstrap method by menzel2021bootstrap which does not require the asymptotic gaussianity. We will hereafter refer to the first method by `CGM' and the second method by `M' for brevity.
Figure (ref) illustrates QQ plots of $\widehat\theta_{NT}$ for $\delta \in \{0.0,0.5,1.0\}$, $J \in \{1,50,100\}$ and $\phi = 0.5$. The top left plot is associated with the most pathetic case with small $|\delta|$ and small $J$, for which our theory cannot guarantee that the gaussian approximation is asymptotically reasonable. This QQ plot shows that the sample quantiles of $\widehat\theta_{NT}$ indeed fail to conform with the theoretical quantiles. On the other hand, when $|\delta|$ or $J$ takes a larger value, the gaussian approximation becomes more reasonable. Specifically, the remaining eight QQ plots in Figure (ref) show that the sample quantiles are sufficiently close to the theoretical quantiles. These results support our theoretical prediction from Section (ref).
Figure (ref) illustrates coverage frequencies as a function of $\phi \in [0.0,1.0]$ for $J \in \{1,50,100\}$ and $\delta \in \{0.0,0.5,1.0\}$, with the nominal coverage probability of 95%. The positions of the nine coverage plots in Figure (ref) correspond to those of the nine QQ plots in Figure (ref). Thus, the top left plot is associated with the most pathetic case in which $\widehat\theta_{NT}$ is far away from gaussian. In this plot, `CGM' suffers from severe under-coverage and `M' in contrast suffers from over-coverage. The under-coverage by `CGM' can be explained by the non-gaussianity of $\widehat\theta_{NT}$.
In the remaining eight plots in Figure (ref), on the other hand, `CGM' achieves fairly precise coverage. Again, as predicted by our theory from Section (ref) and also evidenced by the QQ plots in Figure (ref), these remaining eight cases admit reasonably close gaussian approximations for the distribution of $\widehat\theta_{NT}$. Therefore, these precise coverage results by `CGM' in these eight cases meet our expectations.
Figure (ref) discretely varies $J$ and $\phi$ while it continuously varies $\phi$. We next change our perspectives by continuously varying $\delta$. Figure (ref) illustrates coverage frequencies as a function of $\delta \in [0.0,1.0]$ for $J \in \{1,50,100\}$ and $\phi \in \{0.0,0.5,1.0\}$. As before, the nominal probability is 95%.
On the three plots in the left column of Figure (ref), where we set $J=1$, `CGM' suffers from under-coverage and `M' suffers from over-coverage in the region where $\delta$ takes small values. As $\delta$ increases, however, their coverage accuracy improves. On the remaining six plots in the middle and right columns of Figure (ref), where we set $J=50$ and $100$, respectively, `CGM' achieves fairly accurate coverage regardless of the values taken by $\delta$. These results again conform with our theoretical prediction of asymptotic gaussianity, which holds unless both $|\delta|$ and $J$ are too small.
Thus, far, we have analyzed the simulation results by continuously varying $\phi$ (Figure (ref)) and $\delta$ (Figure (ref)). We next present simulation results with finer variations of $J$. Figure (ref) illustrates coverage frequencies as a function of $J \in \{1,\cdots,100\}$ for $\delta \in \{0.0,0.5,1.0\}$ and $\phi \in \{0.0,0.5,1.0\}$.
On the three plots in the left column of Figure (ref), where we set $\delta=0.00$, `CGM' suffers from under-coverage and `M' suffers from over-coverage in the region where $J$ takes small values. As $J$ increases, however, their coverage accuracy improves. On the remaining six plots in the middle and right columns of Figure (ref), where we set $\delta=0.50$ and $1.00$, respectively, `CGM' achieves fairly accurate coverage regardless of the values of $J$. These results again conform with our theoretical prediction of asymptotic gaussianity, which holds unless both $|\delta|$ and $J$ are too small. We emphasize that gaussian limiting distributions still provide good approximations in most cases even when $\phi=0$ and $\delta=0$, i.e. even when only the degenerate two-sample $U$-statistics are present in the DGP.
In summary, the asymptotic gaussian approximation appears reasonable except for the pathetic case in which both $|\delta|$ and $J$ are too small, as predicted by our theory from Section (ref). Consequently, the conventional two-way cluster-robust standard errors CGM2011,thompson2011simple, which presume the asymptotic gaussianity, perform sufficiently well unless both $|\delta|$ and $J$ are too small. Recall that we have already focused on the data generating processes that correspond to the potentially non-gaussian components in menzel2021bootstrap. Even within this class with the potential non-gaussianity, we confirm that gaussianity is in fact rather common than exceptional.
Thanks to the earlier work in the literature, we now have a more complete understanding of the settings and conditions under which the asymptotic gaussianity is guaranteed. To our knowledge, there are three other scenarios in addition to our conditions in which the gaussianity holds under two-way clustering. Let us maintain $N\sim T$. Table (ref) summarizes them.
The most well-known case is the “non-degeneracy” condition imposed in e.g. DDG2019, which assumes that at least one of $E[D_{it}|\alpha_i]$ and $E[D_{it}|\gamma_t]$ is random with its variance bounded away from zero. This condition imposes that at least one of the latent cluster-specific shocks has an impact on the level of the conditional mean. If this condition is satisfied, then $\widehat\theta_{NT}\approx L_{NT}$ and the asymptotic gaussianity holds regardless of how the $W_{NT}$ term behaves.
The next case is illustrated by the “sparse network” asymptotics, which was first considered in bickel2011method for dyadic network graphs, and studied under two-way clustering by graham2022sparse. In this case, one assumes $D_{it}\in\{0,1\}$ and $p:=P(D_{it}=1)$ is converging to zero or one. Under this setting, one can exploit the binary nature of the outcome variable and show $W_{NT}=O_p(p/N)$ to be asymptotically dominated by $R_{NT}=O_p(p^{1/2}/N)$, and thus $\widehat \theta_{NT}\approx L_{NT}+R_{NT}$. Therefore, the asymptotic gaussianity is guaranteed by an application of a martingale CLT. The third case is related to the second one and considers kernel density estimation for two-way clustered random variables, studied in graham2022kernel.
In light of these conditions for the asymptotic gaussianity along with our CLT result, we propose the decision tree shown in Figure (ref) to researchers considering adoption of the TWCR standard errors for inference.
First, decide whether the researcher is willing to assume the non-degeneracy or a sparse network. If the answer is yes, then one can use the TWCR standard errors. Otherwise, ask if the model under consideration involves only very few latent factors. If the answer is no, then one can still use the TWCR standard errors. In many economic applications, researchers would rather want to suppose that there are many latent factors in economic agents, e.g., many dimensions of ability measures, many dimensions of health measures, and many dimensions of cultural attributes. Our theoretical result endorses the validity of the TWCR standard errors in these applications that are common in the economic research.
In summary, we conclude that the two-way cluster-robust (TWCR) standard errors, which were proposed by CGM2011 and thompson2011simple and have been used by thousands of empirical research papers, are valid under most, if not all, circumstances of the economic research.