Analytic inference with two-way clustering

Analytic inference with two-way clustering

\@startsection{section}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Introduction}

Applied researchers are often reluctant to assume independence between units, because these units may be affected by common shocks. Moreover, these shocks may be of different nature. For instance, the wages of two individuals could be correlated either because these individuals belong to the same industry or because they live in the same area. This case is referred to as two-way clustering because clustering occurs along two dimensions, industry and geographical area in this example. To account for such possible dependence, researchers routinely apply the variance estimator of miglioretti2007marginal (ChenRao2007, MH hereafter), cameron2011 (2011, CGM hereafter) and thompson2011, denoted by $\widehat{V}_u$ below.

However, these “usual” variance estimators have a major drawback, namely, they may not be positive, neither in finite samples nor even asymptotically. Similarly, variance matrices may not be semidefinite positive. We find evidence of this in applied work (see Appendix (ref) for details). CGM propose a fix to this issue, but the corresponding inference in linear regressions is not invariant to affine transforms of the variables anymore. In simulations, we find that the level of a test can vary from zero to one as we change the location and scale of a regressor.

Theoretically speaking, menzel2021bootstrap showed that in a particular data generating process (DGP), standard variance estimators are negative with probability almost 40% asymptotically, which leads to severely distorted inference, even if we use CGM's fix. Beyond this particular example, failure of usual inference occurs whenever estimators are not asymptotically Gaussian, a situation that arises naturally with multiway clustering, as illustrated in menzel2021bootstrap and discussed multiple times below.

The aim of this paper is to suggest an elementary fix for inference, which eliminates the issue of negative variances, remains invariant to affine transforms of the covariates in linear regressions and is asymptotically valid.\footnote{We develop the Stata package twc_inf, available on SSC, which implements this method for linear, probit, logit and poisson regressions.} Consider a univariate equality test. Let $\widehat{V}_1$ and $\widehat{V}_2$ be the variance estimators obtained assuming that only one of the two dimensions of clustering matters, and let $\text{se}_1$ and $\text{se}_2$ be the associated standard errors. Then, we suggest to use as a standard error the maximum between $\text{se}_1$, $\text{se}_2$ and $\text{se}_u$, where the latter is the standard error associated to $\widehat{V}_u$, with the understanding that $\text{se}_u=0$ if $\widehat{V}_u$ is negative. This modification has also been proposed by mackinnon2024jackknife, though they do not establish its validity in cases where the usual method fails. We suggest a similar construction for joint hypothesis testing, which is new to our best knowledge.

We first focus on equality tests for univariate expectations and establish the asymptotic validity of our procedure both in a pointwise and uniform sense. To do so, we model the data as a dissociated, separately exchangeable array, following in particular ddg2021 and menzel2021bootstrap. Then, we rely on results for such arrays, in particular the so-called Aldous-Hoover-Kallenberg representation Aldous1981,Hoover1979,kallenberg1989. Our main insight is that even if the sample mean may not be asymptotically Gaussian and $\widehat{V}_1$ or $\widehat{V}_2$ may remain random asymptotically (once properly normalized), the distribution of our $t$-statistic is asymptotically more concentrated than a standard Gaussian distribution. As a result, the pointwise validity of our univariate tests holds under no further restriction on the DGP. Moreover, we show that our test is equivalent to the usual test whenever the usual $t$-statistic is asymptotically standard Gaussian. Hence, our method does not lead to any power loss asymptotically in cases where usual inference is justified.

The results on univariate means extend to functionals of GMM estimators, as long as a condition that has been overlooked so far holds. Specifically, if linear combinations of the empirical moments under consideration, evaluated at the true parameter, do not all converge at the same rate, the GMM estimator may not be close to the average of the “standard” influence functions, namely the influence functions we use for i.i.d. data. We illustrate this with a simple linear regression example. This issue has consequences for any analytic inference method, ours and the usual one included.\footnote{In our simple linear regression example, usual inference is highly distorted while our method is not, though our theoretical results do not cover this case.} On the other hand, if all linear combinations of the empirical moments converge at the same rate, this peculiar phenomemon disappears and we show the validity of our testing procedure.

We also consider equality tests for multivariate expectations. Unlike in the i.i.d. setup, this extension is not straightforward, however, since the properly normalized matrices $\widehat{V}_1$ and $\widehat{V}_2$ may converge to random and singular matrices, an issue that also affects standard inference and has not been identified yet, to the best of our knowledge. To handle these challenges, we impose sufficient and (partly) necessary conditions on the DGP and resort to a new result on Gaussian matrices that is of independent interest (see Lemma (ref) in Appendix (ref)).

Finally, we compare in simulations our method with the usual one and the bootstrap method of menzel2021bootstrap. We show in particular that usual inference can be very distorted, while ours seems to perform well even in cases not covered by our theory.

\paragraph{Related literature.}

First and foremost, our paper contributes to the literature on analytic inference under multiway clustering. As mentioned above, the variance estimator $\widehat{V}_u$ was proposed by miglioretti2007marginal, cameron2011 and thompson2011. These papers do not show the validity of the corresponding inference. The latter is established by menzel2021bootstrap for sample means of univariate variables if such means are asymptotically Gaussian. menzel2021bootstrap also shows that if sample means are not asymptotically Gaussian, inference based on $\widehat{V}_u$ may not be valid. chiang2023using extend Menzel's results by showing asymptotic Gaussianity for specific drifting sequences of DGPs. Another extension of Menzel's result, to large $T$ panel data where temporal shocks can be dependent both over time and across individuals, is considered by Chiang2024. Yap2025_1 shows the validity of usual inference under the same independence structure as here, but without exchangeability. Compared to these papers, we show that a simple modification of inference based on $\widehat{V}_u$ solely also works in non-Gaussian cases, while being equivalent to it in Gaussian cases. To our knowledge, this is the first analytic inference method for which validity is established in non-Gaussian cases.

Several papers also consider resampling-based inference, and here we just mention a few of them. ddg2021 show the validity of the so-called pigeonhole bootstrap, and a multiplier bootstrap, for “non-degenerate” DGPs, for which the estimator under consideration converges at a slow rate. mackinnon2021wild show the validity of a certain wild bootstrap method in some Gaussian regimes. menzel2021bootstrap develop other wild bootstrap schemes and show that one of them is pointwise valid both in Gaussian and non-Gaussian regimes, while another one controls size over a large set of DGPs but is possibly conservative juodis2021. Our paper complements Menzel's by showing that to some extent, adaptivity is also possible with analytic inference in this set-up. Our approach also has the advantage of being computationally very cheap and not requiring any tuning parameter.

\paragraph{Organization of the paper.}

Section (ref) introduces the setup and the tests on univariate expectations we propose, and presents our pointwise and uniform results for these tests. Section (ref) extends our approach to the GMM case. Section (ref) analyses differences between our method and others, in particular $\widehat{V}_u$, in simulations. Section (ref) concludes. The appendix gathers some extensions, such as equality tests for multivariate expectations, and most of the proofs. The remaining proofs and supporting lemmas can be found in the Online Appendix.

\@startsection{section}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Inference on scalar expectations}

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Set-up and inference method}

We have access to the observed random variables $(Y_{ij})_{1\le i\le C_1,1\le j\le C_2}$. The two indices $i$ and $j$ correspond to the two dimensions of “clustering”, with a dependence structure that will be clarified below. For instance, these dimensions may correspond to industries and geographical areas. We would like to make inference on a scalar parameter $\theta_0 := E[\overline{Y}]$, where, for any array of random variables $(D_{ij})_{1\leq i \leq C_1, 1 \leq j \leq C_2}$, we let $\overline{D}:=(C_1C_2)^{-1}\sum_{i=1}^{C_1}\sum_{j=1}^{C_2} D_{ij}$. Note that for simplicity, we assume in this section to observe a single random variable in each “cell” $(i,j)$. If, instead, we observe $(Y_{ij\ell})_{\ell=1,...,N_{ij}}$ in cell $(i,j)$ and if $\theta_0=E[(C_1C_2)^{-1}\sum_{i,j}\sum_{\ell=1}^{N_{ij}}Y_{ij\ell}]$, we can always define $Z_{ij}$ as $\sum_{\ell=1}^{N_{ij}} Y_{ij\ell}$ and $\theta_0$ as $E[\overline{Z}]$. For other parameters of interest (e.g., $\theta_0=E[\overline{Z}]/E[\overline{N}]$, the estimand corresponding to $\sum_{i,j}\sum_{\ell}^{N_{ij}}Y_{ij\ell}/\sum_{i,j} N_{ij}$), assuming a single observation per cell is not without loss of generality, and we do allow for this case in our GMM setup in Section (ref).

Hereafter, we mostly consider tests of nominal level $\alpha\in (0,1)$ of the null hypothesis that $\theta_0 = \theta$, against $\theta_0\ne \theta$; we also briefly discuss unilateral tests, as well as confidence intervals. As those proposed by MH and CGM, our tests rely on the following three variance estimators:

For simplicity, and since they do not matter asymptotically, we do not consider the degrees-of-freedom corrections suggested by CGM.

We first present the test proposed by MH and CGM. Let $\widehat{V}_u:=\widehat{V}_1+\widehat{V}_2-\widehat{V}_{12}$, where the index “u” refers to “usual”. Then, MH and CGM consider the test $\phi_{u,\alpha} := \mathds{1}\left\{|t_u|>z_{1-\alpha/2}\right\}$, where $z_{1-\alpha/2}$ is the quantile of order $1-\alpha/2$ of a standard normal distribution and $$t_u := \frac{\overline{Y}-\theta}{\widehat{V}_u^{1/2}}.$$ This approach has one major drawback, namely, $\widehat{V}_u$ can be negative, in which case the test above is not defined. This may happen in finite samples even in DGPs for which this test is asymptotically valid and thus $P(\widehat{V}_u>0)$ tends to one. There are also cases for which the latter condition does not hold and the test is asymptotically invalid. Suppose for instance that $Y_{ij}=\theta_0+U_{i0}U_{0j}$ with $(U_{i0})_{i\ge 1}$ and $(U_{0j})_{j\ge 1}$ two independent sequences of centered, i.i.d. variables. menzel2021bootstrap shows that in this example, $P(\widehat{V}_u<0)\to 39.3\%$. This multiplicative structure is a special case of a factor model (in fact, similar results would hold if we replaced $U_{i0}U_{0j}$ by $\sum_{k=1}^K U_{i0,k}U_{0j,k}$) and thus seems relevant in economic contexts.

To solve these issues, we propose the following modification, which was suggested before by mackinnon2024jackknife. Let $\text{se}_k:=\widehat{V}_k^{1/2}$ for $k\in\{1,2\}$, $\text{se}_u:=\max(0,\widehat{V}_u)^{1/2}$ and let $\text{se}:=\max(\text{se}_1,\text{se}_2,\text{se}_u)$. Then, consider the test

with the convention that $\phi_{\alpha}=1$ when $\text{se}=0$. Remark that we simply replace the usual standard error $\text{se}_u$ with the maximum of $\text{se}_u$, $ \text{se}_1$ and $\text{se}_2$. An intuition behind this test is that the second (resp., the first) dimension of clustering may not matter. In such a case, it would be more natural to consider $\text{se}_1$ (resp. $\text{se}_2$) rather than $\text{se}_u$. We then take a conservative approach by picking the maximum of these three standard errors. It turns out, however, that when $\overline{Y}$ is asymptotically Gaussian, our test is not asymptotically conservative. Our test is (potentially) conservative only in non-Gaussian cases, for which $\phi_{u,\alpha}$ may be asymptotically invalid. In the rest of paper, we do not discuss unilateral tests nor confidence intervals on $\theta_0$ (given by $[\overline{Y} \pm \, z_{1-\alpha/2} \, \text{se}]$), but our results on bilateral tests below directly extend to them.

\paragraph{Multivariate extension.} Suppose that $Y_{ij}\in\mathbb R^d$, and we want to test $E[\overline{Y}]=\theta$. In Appendix (ref), we consider a test in a same spirit as above. Specifically, we basically take the minimum of three $F$-tests, with an adjustment in cases some of the involved variance matrices are singular. To our knowledge, this test is new.

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Assumptions}

We obtain our results below under two conditions, which put restrictions on the data generating process and the asymptotic framework. The first assumption clarifies the dependence structure underlying two-way clustering and imposes minimal moment conditions:

The first condition implies that two subsets of $\bm{Y}$ sharing no common cluster are independent. On the other hand, this condition does not impose any restriction on the dependence between $Y_{ij}$ and $Y_{ij'}$ or between $Y_{ij}$ and $Y_{i'j}$. The second condition states that the labels $i$ and $j$ do not carry any information: replacing them by any other labelling (through permutations) leads to the same distribution of the array. This implies in particular that the variables $(Y_{ij})_{i,j \ge 1}$ are identically distributed. Imposing finite second moments and non-zero variance allows us to rule out pathological situations without any variation in the observed data and is required to derive the asymptotic properties of our tests. Finally, allowing the distribution of $\bm{Y}$ to depend on $(C_1,C_2)$ is essential when studying the (asymptotic) uniform validity of our inference method, as in this case we must study the asymptotic behavior of our test under sequences of DGPs rather than a fixed DGP.

Our second assumption pertains to the asymptotic framework. We suppose hereafter that both $C_1$ and $C_2$ tend to infinity, but without restricting their respective rates of convergence.

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{A useful decomposition}

Before showing our results, we present a useful decomposition. We assume that Assumption (ref) holds and introduce the variable $W_{ij} := Y_{ij}-\theta_0$ and vector $\bm{W} := (Y_{ij}-\theta_0)_{(i,j)\in \mathbb N^{*2}}$. First remark that as a dissociated and separately exchangeable array, $\bm{W}$ satisfies a Aldous-Hoover-Kallenberg (AHK for short) representation, see Aldous1981, Hoover1979 and kallenberg1989. Namely, there exist i.i.d. continuously distributed random variables $(U_{i0},U_{0j},U_{ij})_{i,j \ge 1}$ and a function $\tau$ such that almost surely,

We can assume without loss of generality (wlog) that $U_{i0}$, $U_{0j}$ and $U_{ij}$ are centered and admit second-order moments. The variables $U_{i0}$ and $U_{0j}$ may be seen as row and column shocks, respectively, while $U_{ij}$ can be interpreted as a “cell”-specific shock. Then, we consider a similar decomposition as that in menzel2021bootstrap. Specifically, let us define

Observe that by construction,

Finally, we define $\Omega_1:=V(\alpha_1)$, $\Omega_2:=V(\beta_1)$, $\Omega_3:=V(\gamma_{11})$ and $\Omega_4:=V(\varepsilon_{11})$. Because the AHK decomposition is not unique, it may seem that $(\alpha_i,\beta_j,\gamma_{ij}, \varepsilon_{ij})_{i,j \ge 1}$ and the $(\Omega_k)_{k=1,...,4}$ depend on the choice of the variables $(U_{i0},U_{0j},U_{ij})_{i,j \ge 1}$. The following lemma shows that this is not the case. Let $\mathcal{S}_{1,n}:=\sigma(W_{ij}:j>n, i\geq 1)$, $\mathcal{S}_{2,n}:=\sigma(W_{ij}:i>n, j\geq 1)$ and $\mathcal{S}_{12,n}:=\sigma(W_{ij}:\max(i,j)>n)$ and $\mathcal{S}_1:=\bigcap_{n\geq 1}\mathcal{S}_{1,n}$, $\mathcal{S}_2:=\bigcap_{n\geq 1}\mathcal{S}_{2,n}$ and $\mathcal{S}_{12}:=\bigcap_{n\geq 1}\mathcal{S}_{12,n}$.

To our knowledge, there is no simple expression for $\Omega_3$, though we can still express it as a function of $\bm{W}$ only through the following equality $$\Omega_3 = V\left\{E[W_{11}|\mathcal{S}_{12}] - E[W_{11}|\mathcal{S}_1] - E[W_{11}|\mathcal{S}_2]\right\}.$$

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Pointwise results}

We now state and discuss validity results for our test when the probability distribution of $\bm{Y}$ does not vary with $n$.

Even if Theorem (ref) follows from our uniform result below (Theorem (ref)), let us give some intuition on its proof. Assume first that $\Omega_1+\Omega_2>0$. In that case, we show that $\overline{Y}-\theta_0=O_p\left((\Omega_1/C_1+\Omega_2/C_2)^{1/2}\right)$ and $$\frac{\overline{Y}-\theta_0}{\text{se}_u}=\left[\frac{\overline{\alpha}+\overline{\beta}}{(\Omega_1/C_1 + \Omega_2/C_2)^{1/2}}+\frac{\overline{\gamma}+\overline{\varepsilon}}{(\Omega_1/C_1 + \Omega_2/C_2)^{1/2}}\right] + o_P(1).$$ By the central limit theorem, the first fraction on the right-hand side converges to a standard normal distribution. Also, observing that $\text{Cov}(\gamma_{ij},\gamma_{i'j'})=V(\gamma_{11})$ $\times \mathds{1}\left\{i=i',j=j'\right\}$ and $\text{Cov}(\varepsilon_{ij},\varepsilon_{i'j'})=V(\varepsilon_{11})\mathds{1}\left\{i=i',j=j'\right\}$, we prove that $\overline{\gamma}+\overline{\varepsilon}=O_p((C_1C_2)^{-1/2})$. As a result, $$\frac{\overline{Y}-\theta_0}{\text{se}_u} \stackrel{d}{\longrightarrow} \mathcal{N}(0,1).$$ This proves the asymptotic validity of usual inference, as well as asymptotic normality of $\overline{Y}$, if $\Omega_1+\Omega_2>0$. Moreover, we show that $\text{se}/\text{se}_u\stackrel{p}{\longrightarrow} 1$ (see Eq. (ref) in the appendix), which implies that our test is also asymptotically valid in this case, and in fact equivalent to the usual test.

Next, assume that $\Omega_1+\Omega_2=0$. Then, $\overline{\alpha}=\overline{\beta}=0$. Since we still have $\overline{\gamma}+\overline{\varepsilon}=O_p((C_1C_2)^{-1/2})$, we obtain

Equation (ref) implies the estimator converges at a faster rate when $\Omega_1+\Omega_1=0$. If $\Omega_3=0$, then $\overline{\gamma}=0$ and $(C_1C_2 /\Omega_4)^{1/2} \overline{\varepsilon}\stackrel{d}{\longrightarrow} \mathcal{N}\left(0,1\right)$ ensuring that $\overline{Y}$ is again asymptotically normal. Moreover, we establish that $(C_1C_2/\Omega_4)^{-1} \text{se} \stackrel{p}{\longrightarrow} 1$ and $\text{se}_u/\text{se}\stackrel{p}{\longrightarrow} 1$. Thus, in this case again, the usual test and ours are equivalent. Moreover, they are both asymptotically valid and non-conservative.

Finally, if $\Omega_3>0$, two complications occur. First, $\overline{\gamma}$ is not asymptotically normal and second, the standard errors remain random asymptotically. The key point we establish is that conditional on $(U_{0j})_{j\ge 1}$ (say), we have

Since the limit (Gaussian) distribution in (ref) does not depend on the $(U_{0j})_{j\ge 1}$, we obtain unconditional convergence as well. Combined with (ref), this yields $$\frac{\overline{Y}-\theta_0}{\text{se}_1} \stackrel{d}{\longrightarrow} \mathcal{N}(0,1).$$ We finally obtain (ref) using the fact that $\text{se}\ge \text{se}_1$ and that $\theta_0=\theta$ under the null hypothesis.

\paragraph{Asymptotically exact tests.} Our test $\phi_{\alpha}$ is conservative in non-Gaussian regimes. It is actually possible to consider an asymptotically exact test. To understand how, remark that $t_u$ is asymptotically exact when $\Omega_1+\Omega_2>0$, in which case $\overline{Y}$ has a slow rate of convergence, whereas in view of (ref) and (ref), the test based on $t_1:=(\overline{Y}-\theta)/\text{se}_1$ is asymptotically exact when $\Omega_1+\Omega_2=0$. Moreover, we show in the proof of Theorem (ref) that $(\widehat{V}_1+\widehat{V}_2)/\widehat{V}_{12}$ converges to infinity when $\Omega_1+\Omega_2>0$, whereas $(\widehat{V}_1+\widehat{V}_2)/\widehat{V}_{12} =O_P(1)$ when $\Omega_1+\Omega_2=0$. Now, consider

where $\underline{C}:=\min(C_1,C_2)$ and $s_{\underline{C}}$ is such that $s_{\underline{C}}\to \infty$ and $s_{\underline{C}}/\underline{C} \to 0$. Such conditions ensure that one selects the statistic that is asymptotically exact with probability approaching one.\footnote{In this sense, this construction is related to Menzel's bootstrap with selection.} As a result, the corresponding test is also asymptotically exact. Remark also that to treat the dimensions of clustering symmetrically, one could replace $t_1$ by $t_{j_{\textrm{max}}}$, with $j_{\textrm{max}}:=\arg\max_{k=1,2}C_k$. However, this test suffers from at least two drawbacks. First, the precise choice of the tuning parameter $s_{\underline{C}}$ remains unclear. Second, the test associated with $t_a$ does not have uniform guarantees, contrary to $t$.

\paragraph{Power loss.} Related to the previous point, we explore in Appendix (ref) to what extent our test $\phi_\alpha$ is conservative, by computing the average increase in confidence intervals we obtain when using $\text{se}$ instead of $\text{se}_1$ on asymptotically non-Gaussian DGPs. Across multiple draws of possible DGPs, we obtain an average increase of the length of around 9%, with a maximum of around 25%.

\paragraph{Multivariate extension.} We derive in Appendix (ref) the asymptotic validity of the joint test mentioned above. We obtain therein a similar result as Theorem (ref), under an additional condition ensuring that the limit of estimated variance matrices have full rank. Note that this latter issue is specific to the multivariate set-up.

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Uniform results}

We now consider a uniform version of Theorem (ref). In this context, we have to make some of our previous conditions uniform. Let $\mathcal{P}$ denote the set of probability distributions such that Assumption (ref) holds. By Theorem (ref), we automatically obtain “uniform” asymptotic validity on $\mathcal{P}$ as long as it is finite; but additional restrictions have to be imposed otherwise.

To introduce these restrictions, we index relevant objects such as expectation signs or $\Omega_k$ by $P$. For any $\tau_P$ that satisfies Equation (ref), let us define\footnote{Here and in the proofs, the variables $(U_{ij})_{i,j\ge 0}$ are supposed wlog to be uniformly distributed on $[0,1]$.}

We define $\tau_{2P}$ similarly, just replacing $\Omega_{1P}$ and $U_{1,0}=u_1$ by $\Omega_{2P}$ and $U_{0,1}=u_2$. For any $m>0$ and $H$ compact subset of $L_2([0,1]^3,\mathbb R)$, let us introduce

The compactness restriction states that we can approximate elements of $H$ uniformly well by elements of a finite-dimensional space. We comment on the other restrictions in $\mathcal{P}_{m,H}$ and $\mathcal{P}_{m,H}^G$ below.

Let us sketch the proof of Theorem (ref). First, we show that it suffices to establish the result for any sequence of DGPs $(P_n)_{n\ge 1}$ in $\mathcal{P}_{m,H}$ (or in $\mathcal{P}_{m,H}^G$). The difficulty, then, is that for such a sequence, the four terms in the decomposition (ref) may matter asymptotically. To illustrate this, consider the following sequence of DGPs:

where $C_1=C_2=n$, $(b_1,...,b_4)\in\mathbb R^4$ and the $(U_{ij})_{i,j\ge 0}$ are i.i.d., mean-zero variables. By dividing $b_1 U_{i0} + b_2 U_{0j}$ by $n^{1/2}$, we make the four terms of the decomposition ($\overline{\alpha}$, $\overline{\beta}$, $\overline{\gamma}$ and $\overline{\varepsilon}$) converge at the same rate, namely $n=(C_1C_2)^{1/2}$. The term $\overline{\alpha}+\overline{\beta}+\overline{\varepsilon}$ is asymptotically normal but $\overline{\gamma}$ is not, and it is not asymptotically independent of the first term. The general asymptotic distribution of $\overline{Y}$ for such sequences of DGPs is complicated and given by Lemma (ref) in Appendix (ref). Still, we can explain the logic of our results in the simple example given by (ref). Specifically, Lemma (ref) implies that $$\bigg[(\overline{Y}-\theta_{0P})/\sqrt{V(\overline{Y})}, (\widehat{V}_1, \widehat{V}_2, \widehat{V}_{12})/V(\overline{Y})\bigg] \stackrel{d}{\longrightarrow} (L, V_1, V_2, V_{12}),$$ where, letting $(Z_1, Z_2, Z_4)$ be three i.i.d. standard normal variables,

with $(c_1,...,c_4)$ a vector that is related to $(\Omega_{1P},...,\Omega_{4P})$ (see Lemmas (ref) and (ref) for details) and also to $(b_1,...,b_4)$. In particular, if $b_2=0$, which implies $\Omega_{2P}=0$, we have $c_2=0$. Then, $L|Z_2 \sim \mathcal{N}(0, c_4^2 + (c_1 + c_3 Z_2)^2)$. As a result, $L/V_1^{1/2} \sim \mathcal{N}(0,1)$ and thus, as in the non-normal, pointwise case,

Similarly, if $b_1=0$, so that $\Omega_{1P}=0$, we can show that $(\overline{Y}-\theta_{0P})/\text{se}_2 \stackrel{d}{\longrightarrow} \mathcal{N}\left(0,1\right)$. The conclusion on (ref) follows as in the pointwise case.

To what extent are the conditions in $\mathcal{P}_{m,H}$ necessary? Without fully answering this question, we can at least ascertain that the last condition in $\mathcal{P}_{m,H}$, namely $\Omega_{1P} \wedge \Omega_{2P} = 0$ or $\Omega_{3P} \le m^{-1} \left(\Omega_{1P}+\Omega_{2P}\right)$, cannnot be omitted. To see this, let us consider the following particular case of (ref):

for some $\zeta\in\mathbb R$ and an i.i.d. sequence $(U_{ij})_{i,j\ge 0}$ with standard normal distribution. Remark that no $\mathcal{P}_{m,H}$ includes the full sequence $(P_n)_{n\ge 1}$, since for all $n\ge 1$ $\Omega_{1P_n} \wedge \Omega_{2P_n} > 0$, $\Omega_{3P_n}=\zeta$ and $\Omega_{1P_n}+\Omega_{2P_n}\to 0$. Now, using Lemma (ref), we are able to simulate the asymptotic distribution of the test statistic $t$ in this case, for any $\zeta\in\mathbb R$. It appears that the test is not asymptotically valid for $\zeta\in(0,1.16]$, with an asymptotic level peaking at around $11\%$ for $\zeta\simeq 0.65$. Mathematically, the expressions for $L$ and $V_1$ above show that $L|Z_2 \sim \mathcal{N}(c_2 Z_2, V_1)$. Moreover, $c_2\ne 0$ and thus we do not obtain (ref) anymore.

\@startsection{section}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Inference based on GMM estimators}

We extend our results on scalar expectations to tests based on smooth functionals of GMM estimators. Because we consider possibly nonlinear estimators, we have to explicitly account for the fact that several units may be observed in each cell $(i,j)$. We thus assume to observe $(A_{ij\ell})_{1\le i\le C_1,1\le j\le C_2, 1 \le \ell \le N_{ij}}$, where $N_{ij}\in \mathbb N$ is a random variable denoting the number of units observed in cell $(i,j)$. We are interested in testing $\theta_0=\theta$ against $\theta_0\neq \theta$, where $\theta_0:=\varphi(\beta_0)\in \mathbb R$ for some function $\varphi$ and $\beta_0\in \Theta \subseteq \mathbb R^p$ satisfies

with $\psi(a,\beta) \in \mathbb R^q$ $(q\ge p)$, and we use the convention $\sum_{\ell=1}^0 t_{\ell}=0$ for any sequence $(t_{\ell})_{\ell\ge 1}$. We let $\widehat{\theta}:=\varphi(\widehat{\beta})$, where $\widehat{\beta}$ is the GMM estimator:

for some symmetric, positive matrix $\Upsilon_n$ and with $\psi_{ij}(\beta) :=\sum_{\ell=1}^{N_{ij}} \psi(A_{ij\ell},\beta)$. To build our test, let $J:=E\left[\partial \psi_{11}(\beta_0)/\partial \beta\right]'$ and $\Psi_{ij} := -\partial \varphi(\beta_0)/\partial \beta'\left(J'\Upsilon J \right)^{-1} J'\Upsilon$ $\psi_{ij}(\beta_0)$ (we assume below the existence of the derivatives). As with i.i.d. data, we expect to have

Then, let $\widehat{\Psi}_{ij} := -\partial \varphi(\widehat{\beta})/\partial \beta'\left(\widehat{J}'\Upsilon_n \widehat{J} \right)^+ \widehat{J}'\Upsilon_n \psi_{ij}(\widehat{\beta})$, with $A^+$ the Moore-Penrose inverse of $A$ and $\widehat{J} := (C_1C_2)^{-1}\sum_{i=1}^{C_1}\sum_{j=1}^{C_2}\partial \psi_{ij}(\widehat{\beta})/\partial \beta'$. We consider a similar test as in (ref), namely

We compute $\text{se}$ as in Section (ref) except that in view of (ref), we replace $Y_{ij}$ with $\widehat{\Psi}_{ij}$ in $\widehat{V}_1, \widehat{V}_2$ and $\widehat{V}_{12}$.

The validity of this test is obtained under the following conditions:

Assumption (ref) is simply Assumption (ref) (without moment constraints) imposed on $\bm{A}^\infty$. Assumption (ref) includes classical regularity conditions that are not specific to our setup with multiway clustering. But we also impose Assumption (ref), which is specific to our setup: it automatically holds with i.i.d. data under Assumption (ref)-(iii). We discuss this condition further below.

The proof of Theorem (ref) can be found in the Supplemental Appendix. In line with Section (ref), we could strengthen our pointwise results on GMMs to uniform ones, by basically imposing uniform versions of Assumptions (ref) and (ref).

While Assumption (ref) is a standard regularity assumption, Theorem (ref) also relies on Assumption (ref). The following example illustrates that without this assumption, the usual linear approximation (ref) may not be valid with two-way clustered data. Note that this issue affects $\widehat{\theta}$ and is thus not specific to our inference method.

\paragraph{Other tests.} With GMMs, the naive test based on $\widehat{V}_u=\widehat{V}_1+\widehat{V}_2-\widehat{V}_{12}$ has the same drawbacks as with sample means: $\widehat{V}_u$ can be negative, and inference may be invalid asymptotically because $\widehat{\theta}$ may not be asymptotically Gaussian. The first issue may actually be common in practice: it arises in 9 out of 15 papers published in the {\it American Economic Review} and that use two-way clustering (see Appendix (ref) for details). To solve this issue, CGM consider the following alternative estimator. First, they estimate the asymptotic variance of $\widehat{\beta}$ by the matrix $\widehat{V}^\beta_u$, obtained as $\widehat{V}_u$ but with the $((\widehat{J}'\Upsilon_n \widehat{J})^+\widehat{J}'\Upsilon_n$ $\psi_{ij}(\widehat{\beta}))_{1\le i\le C_1, 1\le j\le C_2}$ instead of the $(\widehat{\Psi}_{ij})_{1\le i\le C_1, 1\le j\le C_2}$ as inputs. Then, they consider the eigendecomposition of $\widehat{V}^\beta_u$, $P'\Delta P$, and replace the negative eigenvalues in $\Delta$ by 0, and let $\widetilde{V}^\beta_u$ denote the corresponding matrix. Finally, their variance estimator for $\widehat{\theta}$ is $[\partial \varphi(\widehat{\beta})/\partial \beta]' \,\widetilde{V}^\beta_u\, [\partial \varphi(\widehat{\beta})/\partial \beta]$.

This solution has two drawbacks, however. First, it does not restore valid inference when $\widehat{\theta}$ is not asymptotically Gaussian. Second, the test is not invariant to affine transforms of variables in linear regressions, contrary to ours. For instance, we show in Subsection (ref) that adding a constant or changing the scale of a regressor can make the rejection rate vary from 0 to 1. Similarly, changing the reference of a binary regressor affects inference.

\@startsection{section}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Monte Carlo simulations}

We illustrate the performance of our test in two cases: univariate means and linear regressions. In these two cases, we let $C_1=C_2=n \in\{10, 20, 40\}$.

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Univariate sample means}

We first consider $\theta_0= E[Y_{11}]$ estimated by the sample mean $\overline{Y}$, where $$Y_{ij}=\delta_{1n} U_{i0} + \delta_{2n} U_{0j} + U_{i0} U_{0j} + \frac{1}{2} U_{ij},$$ and the $(U_{ij})_{i,j\ge 0}$ are all independent, standard normal variables and $(\delta_{1n}, \delta_{2n})$ are possibly varying with $n$. We consider four DGPs, depending on the values of $(\delta_{1n}, \delta_{2n})$:

We compute rejection rates under the null, by testing for $\theta_0=0$, and under the alternative, by testing for $\theta_0=\theta\ne 0$, with $\theta=0.5$ in DGP1 and $\theta=0.15$ in DGP2 to DGP4. This choice of $\theta$ ensures that power is nontrivial with our sample sizes. We compare our test (“DDG” in the table) with usual inference (“Usual” in the table). Recall that $\text{se}_u=\max(0,\widehat{V}_u)^{1/2}$, so that we automatically reject the null hypothesis with usual inference when $\widehat{V}_u \le 0$. We also consider the bootstrap with selection (BS-S) developed by menzel2021bootstrap. This bootstrap requires a tuning parameter $\kappa_0$: we consider both $\kappa_0=0.05$, as in the programs accompanying menzel2021bootstrap, and a much larger value, $\kappa_0=1.25$.\footnote{In fact, there are two parameters appearing in Menzel's bootstrap with selection, namely $\kappa_a$ and $\kappa_g$ menzel2021bootstrap. But in his simulations, he makes both depend on a single parameter $\kappa_0$, by setting $\kappa_a = \kappa_0 \log(C_1)/C_1$ and $\kappa_g = \kappa_0 \log(C_2)/C_2$. We do not report here the results of his conservative bootstrap (BS-C), which is very conservative in our simulations.}

The results are displayed in Table (ref). As predicted by theory, DDG and usual inference are very close in DGP1, for which the estimator is asymptotically Gaussian and usual inference is valid. For this DGP, the results of the four methods are very similar. In DGP2, on the other hand, the usual variance estimator is negative in around 30% of the samples. Accordingly, the test is highly distorted. Our test is conservative, but less than BS-S with $\kappa=0.05$; its power is similar to that of BS-S with $\kappa=1.25$. In DGP3, our test is again conservative but has higher power than BS-S with $\kappa=0.05$. The bootstrap with $\kappa=1.25$ is the most powerful but slightly overrejects. Usual inference is still distorted, though less so than in DGP2. Finally, in the last DGP, for which we do not have any theoretical guarantee, our test turns out to have a level close to the nominal one. Again, it has slightly larger power than BS-S with $\kappa=0.05$. BS-S with $\kappa=1.25$ slightly overrejects, and usual inference is quite distorted.

The bottom line is that our method compares well in terms of level and power with the bootstrap and has the advantage of not requiring the choice of a tuning parameter, which may be difficult to choose appropriately and does affect rejection rates.

\@startsection{subsection}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Linear regressions}

Second, we consider inference in linear regressions, a simple instance of the GMM models discussed in Section (ref). Specifically, we consider the following: $$Y_{ij} = X_{ij}'\beta_0 + \varepsilon_{ij}, \; E[X_{ij}\varepsilon_{ij}]=0,$$ where we wish to conduct inference on $\theta_0$, the second coefficient of $\beta_0$ (corresponding to the first non-constant element of $X_{ij}$). We assume $\theta_0=0$ and consider again four DGPs (where, as above, $(U_{ij})_{i,j\ge 0}$ and $(\widetilde{U}_{ij})_{i,j\ge 0}$ are two independent families of i.i.d. standard normal variables):

We compute the rejection rates under the null and under the alternative by testing for $\theta_0=\theta\ne 0$, with $\theta=0.3$ in DGP1, $\theta=0.15$ in DGP2 and 3 and $\theta=0.13$ in DGP4. Apart from our test and the usual one, we consider CGM's fix detailed in Section (ref). We also consider Menzel's bootstrap with selection (BS-S). As with univariate sample means, this bootstrap requires a tuning parameter, which we also call $\kappa_0$: we consider both $\kappa_0=10$, as in the programs accompanying menzel2021bootstrap, and a smaller value, $\kappa_0=1$.\footnote{As above, the two tuning parameters $\kappa_a$ and $\kappa_g$ are defined in Menzel's programs as $\kappa_a = \kappa_0 \mu_{4e} \log(C_1)/C_1$ and $\kappa_g = \kappa_0 \mu_{4e}\log(C_2)/C_2$, with $\mu_{4e}=[2\max(1/100,\overline{\widehat{\varepsilon}^4})]^{1/2}$, where $\widehat{\varepsilon}$ denotes the residual of the regression. The choice $\kappa_0=1$ also appears in the programs but is commented.}

The results are displayed in Table (ref). Interestingly, in DGP1 for which the usual inference is asymptotically valid, our test leads to substantial improvements when $n=10$, also over CGM. With $\kappa_0=10$, BS-S does not seem to work properly in this DGP, but using $\kappa_0=1$ yields results broadly similar to those of DDG.

Usual inference is highly distorted in DGP2 to DGP4, with in particular a rejection rate of 1 in DGP4. CGM is less distorted but still rejects between 14% and 56% in these three DGPs. Also, as indicated above, inference based on CGM's fix is not invariant to linear change in the regressors. For instance, we obtain a very conservative test, with a rejection rate of 0 under the null, when adding 2 to the first regressor. Conversely, multiplying this regressor by a constant approaching 0 makes the rejection rate tend to 1. Though our theoretical results do not apply for DGP2 to 4, our test seems to behave well in these cases, with rejection rates below 5% under the null for all sample sizes. The two bootstraps differ in DGP2, with $\kappa_0=10$ leading to conservative inference, but behave very similarly for DGP3 and DGP4. They also appear to slightly overreject with DGP3.

\@startsection{section}{2}{0mm}{-1\baselineskip}{1\baselineskip}{\normalfont}{Conclusion}

We have shown that suitable, elementary changes in the usual inference with two-way clustering may result in pointwise valid tests even in non-Gaussian regimes. With sample means, this holds under the same moment condition as with i.i.d. data. For GMM estimators, a condition on the rates of convergence of the different moments is required. This condition is specific to the two-way clustering setup and ensures that the estimator is asymptotically close to the average of its influence functions. We also show uniform validity of the tests over suitable classes of DGPs.

We leave a few questions for future research. The first is whether we can still obtain asymptotically valid inference under weaker restrictions than those we have imposed. The second is whether our proposal extends to multiway clustering with three or more dimensions of clustering. The third is whether simple, analytic inference for dyadic data is possible, including in non-Gaussian regimes. This may not be straightforward: we show in Appendix (ref) that the fix we use with two-way clustering does not lead to valid pointwise inference in this setup.

\linespread{1.3}\selectfont