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.
105,117 characters · 17 sections · 65 citation commands
Jackknife Inference with Two0.04167em--0.08333em Way Clustering
\onehalfspacing
The use of two\kern 0.04167em-way cluster-robust variance estimators for linear regression models was independently proposed by MH_2006, \citet*{CGM_2011}, and Thompson_2011. Although two\kern 0.04167em-way clustering has been widely used in empirical work, the asymptotic theory to justify it is quite recent. See, among others, \citet*{Davezies_2021,Davezies_2025}, \citet*{MNW_2021}, Menzel_2021, \citet*{Chiang_2022-JBES}, \citet*{CKS_2023}, \citet*{CHS_2024}, and yap_2025. The finite\kern 0.04167em-sample properties of statistical inference are much less well understood for two\kern 0.04167em-way clustering than for one\kern 0.04167em-way clustering. For an up-to\kern 0.04167em-date discussion of the latter, with recommendations for empirical practice, see \citet*{MNW-guide}.
The jackknife variance estimator has been around for a very long time Tukey_1958,Efron_81,Efron-Stein. The cluster-jackknife CRVE (sometimes called the CV$_{\kern -0.08333em3}$ estimator) for linear regression models with one\kern 0.04167em-way clustering was proposed in BM_2002 and has been available in Stata for many years. Nevertheless, it has not been studied or applied much until very recently. In part, this is because BM_2002 followed MW_1985 by computing the CV$_{\kern -0.08333em3}$ estimator in a way that is efficient when all clusters are very small but extremely inefficient when any clusters are large; see \citet*{MNW-bootknife}. This seems to have given many investigators the erroneous impression that CV$_{\kern -0.08333em3}$ is very expensive to compute, even though the Stata implementation uses a method that is reasonably efficient when the number of clusters is not too large. An even more efficient method is discussed in \citet*{MNW-influence} and implemented in the Stata package summclust \citep*{MNW_summclust}.
In (ref), we discuss the linear regression model with two\kern 0.04167em-way clustering. Two existing CRVEs are discussed, along with their theoretical and practical deficiencies. For the CRVE that is theoretically soundest, the chief deficiency is that it may not be positive definite in finite samples. We discuss two ways to overcome this problem. One is the eigen-decomposition method suggested in CGM_2011. The other is a new and extremely simple procedure which can readily be implemented using existing software.
In (ref), we show how to extend the cluster-jackknife CRVEs discussed in MNW-bootknife and Hansen-jack to two\kern 0.04167em-way clustering. Two alternative approaches, based on different jackknife constructions, have been proposed very recently; see \citet*{CMO_2025}, which uses empirical likelihood, and Hounyo-jack. In (ref), we prove that our two\kern 0.04167em-way cluster-jackknife CRVE is consistent. Based on what is known about the finite\kern 0.04167em-sample performance of cluster-jackknife CRVEs for one\kern 0.04167em-way clustering, it seems very likely that inference based on our CRVE will be more conservative, and usually more reliable, than conventional inference in the two\kern 0.04167em-way case as well. Some theoretical arguments to support this conjecture are provided in (ref).
In (ref), we use simulation experiments to study the finite\kern 0.04167em-sample performance of several procedures for inference. Using the cluster-jackknife methods of (ref) in combination with either of the procedures discussed in (ref) often performs much better than existing methods for cluster-robust inference. In (ref), we apply several methods to two empirical examples. The results that we obtain are entirely in accord with the simulations in (ref). We conclude that, while conventional methods probably do not yield reliable inferences for these examples, our preferred methods based on the cluster jackknife probably do. Finally, (ref) concludes.
Consider the linear regression model
where ${\bm{y}}$ and ${\bm{u}}$ are $N\times 1$ vectors of observations and disturbances, ${\bm{X}}$ is an $N\times k$ matrix of covariates, and ${\bm\beta}$ is a $k\times1$ parameter vector. The model is assumed to have two dimensions of clustering, where the numbers of clusters in the two dimensions are $G$ and $H$\kern -0.08333em, respectively. It is illuminating to rewrite \hyperref[{model}]{\tagform@{\ref*{model}}} in terms of the intersections of the two clustering dimensions:
Here the vectors ${\bm{y}}_{gh}$ and ${\bm{u}}_{gh}$ and the matrix ${\bm{X}}_{gh}$ contain, respectively, the rows of ${\bm{y}}$, ${\bm{u}}$, and ${\bm{X}}$ that correspond to both the \th{g} cluster in the first clustering dimension and the \th{h} cluster in the second clustering dimension. Similarly, we use ${\bm{y}}_g$, ${\bm{X}}_g$, and ${\bm{u}}_g$ to denote vectors that contain the rows of ${\bm{y}}$, ${\bm{X}}$\kern -0.08333em, and ${\bm{u}}$ for the \th{g} cluster in the first dimension, and ${\bm{y}}_h$, ${\bm{X}}_h$, and ${\bm{u}}_h$ to denote the corresponding rows for the \th{h} cluster in the second dimension. For example, the vector ${\bm{y}}_g$ contains the subvectors ${\bm{y}}_{g1}$ through ${\bm{y}}_{gH}$.
We use $N_g$ to denote the number of observations in cluster $g$ for the first dimension, $N_h$ to denote the number of observations in cluster $h$ for the second dimension, and $N_{gh}$ to denote the number of observations in the intersection of cluster $g$ in the first dimension with cluster $h$ in the second dimension. We assume that $N_g\ge1$ and $N_h\ge1$. Thus, the number of observations in the entire sample is
Note that some of the intersections may be empty, so that $N_{gh}$ might well equal 0 for some values of $g$ and $h$. The number of non-empty intersections is $I \le GH$.
Various score vectors play key roles in cluster-robust inference. The score vector for the entire sample is ${\bm{s}} = {\bm{X}}^\top{\bm{u}}$. The score subvector for cluster $g$ in the first dimension is ${\bm{s}}_g = {\bm{X}}_g^\top{\bm{u}}_g$, and the score subvector for cluster $h$ in the second dimension is ${\bm{s}}_h = {\bm{X}}_h^\top{\bm{u}}_h$. Thus there are $G$ score vectors ${\bm{s}}_g$ and $H$ score vectors ${\bm{s}}_h$. The score subvector for intersection $gh$ is ${\bm{s}}_{gh}= {\bm{X}}_{gh}^\top{\bm{u}}_{gh}$.
The variance matrix of the scores can always be written as
Under two\kern 0.04167em-way clustering, it must be the case that
but the covariances may be arbitrary when either $g=g'$ or $h=h'$. The variance matrices for the score subvectors ${\bm{s}}_g$, ${\bm{s}}_h$, and ${\bm{s}}_{gh}$ are respectively denoted
From \hyperref[{def Omega}]{\tagform@{\ref*{def Omega}}} and \hyperref[{var matrices}]{\tagform@{\ref*{var matrices}}}, it is evident that
This follows from the inclusion-exclusion principle. The third term in \hyperref[{truesig}]{\tagform@{\ref*{truesig}}} is essential to avoid double\kern 0.04167em-counting, but, as we shall see, it causes practical difficulties for estimating ${\bm{\Sigma}}$.
As usual, the OLS estimator of ${\bm\beta}$ is $\hat{\bm\beta}=({\bm{X}}^\top {\bm{X}})^{-1} {\bm{X}}^\top {\bm{y}}$, and the OLS residual vector is $\hat{\bm{u}}$. The subvectors of $\hat{\bm{u}}$ for cluster $g$, cluster $h$, and the intersection $gh$ are denoted $\hat{\bm{u}}_g$, $\hat{\bm{u}}_h$, and $\hat{\bm{u}}_{gh}$, respectively. From standard arguments for sandwich variance matrices,
where the component matrices are
The empirical analog of \hyperref[{Vtrue}]{\tagform@{\ref*{Vtrue}}} is the three\kern 0.04167em-term two\kern 0.04167em-way CRVE
where the estimators on the right-hand side of \hyperref[{3mat}]{\tagform@{\ref*{3mat}}} correspond to \hyperref[{VG}]{\tagform@{\ref*{VG}}}, \hyperref[{VH}]{\tagform@{\ref*{VH}}}, and \hyperref[{VI}]{\tagform@{\ref*{VI}}} and will be defined shortly. The subscript “1” in $\hat{\bm{V}}_1^{(3)}$ identifies this as a CV$_{\kern -0.08333em1}$ estimator, by analogy with the HC$_1$ estimator of MW_1985. The three component estimators in \hyperref[{3mat}]{\tagform@{\ref*{3mat}}} are based on the empirical score subvectors $\hat{\bm{s}}_g$, $\hat{\bm{s}}_h$, and $\hat{\bm{s}}_{gh}$, which take the same form as the actual score subvectors, but with $\hat{\bm{u}}$ replacing ${\bm{u}}$. Thus they are all CV$_{\kern -0.08333em1}$ estimators:
The leading scalar factors here are analogous to the scalar factor for the usual one\kern 0.04167em-way CRVE. Since some of the intersections may contain no observations, some of the $\hat{\bm{s}}_{gh}$ may not exist. In practice, it may therefore be advisable to replace the double summation in \hyperref[{VhatI}]{\tagform@{\ref*{VhatI}}} with a single summation over all non-empty intersections.
The superscript “(3)” on $\hat{\bm{V}}_1^{(3)}$ in \hyperref[{3mat}]{\tagform@{\ref*{3mat}}} emphasizes that this estimator has three terms, which correspond to the three terms in \hyperref[{truesig}]{\tagform@{\ref*{truesig}}}. Because $\hat{\bm{V}}_I$ is subtracted from the sum of $\hat{\bm{V}}_G$ and $\hat{\bm{V}}_H$, the matrix $\hat{\bm{V}}_1^{(3)}$ is not necessarily positive definite in finite samples. This problem is not trivial, and there is more than one way to deal with it.
One approach, suggested in CGM_2011 and implemented in Stata, Version 18 and later, is to compute the eigenvalues of $\hat{\bm{V}}_1^{(3)}$, say $\lambda_1,\ldots,\lambda_k$. When any of them is not positive, $\hat{\bm{V}}_1^{(3)}$ is replaced by the eigen-decomposition $\hat{\bm{V}}_1^{(3+)} = {\bm{U}}{\bm{\Lambda}}^{\!+}{\bm{U}}^\top$\kern -0.08333em\kern -0.08333em, where ${\bm{U}}$ is the $k\times k$ matrix of eigenvectors and ${\bm{\Lambda}}^{\!+}$ is a diagonal matrix with typical diagonal element $\lambda_j^+ = \max \{ \lambda_j,0 \}$. In practice, it may be numerically safer to compare the eigenvalues with a very small positive number, say $\eta$, and define $\lambda_j^+$ as $\max \{\lambda_j,\eta\}$. In our programs, we use $\eta=10^{-12}$. Doing this ensures that $\hat{\bm{V}}_1^{(3+)}$ is positive definite, albeit just barely so.
This approach is not entirely satisfactory. Wald statistics and $t$-statistics based on $\hat{\bm{V}}_1^{(3+)}$ are computable, but they may be extremely large. Even when this does not happen, and all quantities of interest can be computed using $\hat{\bm{V}}_1^{(3)}$, replacing $\hat{\bm{V}}_1^{(3)}$ by $\hat{\bm{V}}_1^{(3+)}$ can change all the standard errors. Moreover, the standard error of any element of $\hat{\bm\beta}$, say $\hat\beta_j$, is not invariant to nonsingular transformations of the remaining columns of the matrix ${\bm{X}}$\kern -0.08333em. Thus, for example, precisely how fixed effects or other dummy variables are specified may affect the standard error of $\hat\beta_j$, even though $\hat\beta_j$ itself is invariant to such reparametrizations. For instance, if one wanted to control for American state fixed effects, the choice of using either Texas or California as the reference group can change the estimated standard error for the treatment regressor of interest.
A simpler way to avoid the problem that $\hat{\bm{V}}_1^{(3)}$ may not be positive definite is to replace it by the two\kern 0.04167em-term estimator
This estimator has been studied in Davezies_2018. It omits the third term in \hyperref[{3mat}]{\tagform@{\ref*{3mat}}} and therefore involves double\kern 0.04167em-counting. The justification for omitting $\hat{\bm{V}}_I$ is that, under additional regularity conditions, it becomes asymptotically negligible as both $G$ and $H$ tend to infinity. Because
is positive definite, it follows that a Wald statistic or $t$-statistic based on $\hat{\bm{V}}_1^{(2)}$ will always be smaller than the same statistic based on $\hat{\bm{V}}_1^{(3)}$, so that the former is more conservative.
The conditions for consistency of $\hat{\bm{V}}_1^{(2)}$ are stronger than the ones needed for $\hat{\bm{V}}_1^{(3)}$. For example, MNW_2021 shows that whenever the scores are actually independent, or whenever they are only correlated at the intersection level, $\hat{\bm{V}}_1^{(2)}$ yields test statistics that are asymptotically too small. In this case, $\hat{\bm{V}}_G \approx \hat{\bm{V}}_H \approx \hat{\bm{V}}_I$. Therefore,
whereas
Thus, in this case, $\hat{\bm{V}}_1^{(2)}$ is approximately twice as large as $\hat{\bm{V}}_1^{(3)}$, and twice as large as it should be. The use of “$\approx$” in \hyperref[{V2hatx}]{\tagform@{\ref*{V2hatx}}} and \hyperref[{V3hatx}]{\tagform@{\ref*{V3hatx}}} is deliberately informal, since we did not take limits or introduce any factors of the sample size in \hyperref[{Vtrue}]{\tagform@{\ref*{Vtrue}}}. For a rigorous treatment, see MNW_2021. The result \hyperref[{V2hatx}]{\tagform@{\ref*{V2hatx}}} suggests that $\hat{\bm{V}}_1^{(2)}$ is also likely to perform poorly in finite samples when most of the intra-cluster correlation is at the intersection level.
We now propose a third way to avoid cases in which test statistics based on the three\kern 0.04167em-term estimator $\hat{\bm{V}}_1^{(3)}$ are not positive. Our proposal is simply to compute three test statistics and use the one that takes the smallest positive value. For the hypothesis that ${\bm{R}}{\bm\beta}={\bm{r}}$, the three Wald statistics are
The statistic we propose to use is
where ${\rm pos}(W_3)$ equals $W_3$ whenever $W_3$ is positive and $+\infty$ whenever it is either negative or undefined, as it can be when $\hat{\bm{V}}_1^{(3)}$ is not positive definite. By using $W_{\min}$ defined in \hyperref[{Wmin}]{\tagform@{\ref*{Wmin}}}, we not only avoid Wald statistics that are not positive numbers but also Wald statistics that are misleadingly large. In a particular sample, one or more diagonal elements of $\hat{\bm{V}}_I$ may randomly happen to be just a little smaller than the sum of the corresponding elements of $\hat{\bm{V}}_G$ and $\hat{\bm{V}}_H$. Thus $\hat{\bm{V}}_1^{(3)}$ can yield extremely large test statistics which are completely misleading.
In most cases, it is not necessary to calculate the entire $\hat{\bm{V}}_1^{(3)}$ matrix. Only the rows and columns needed for the Wald statistic have to be calculated. Of course, when there is only one restriction, we can use a $t$-statistic instead of a Wald statistic. In this case, we just need to find the largest of the three standard errors and calculate a $t$-statistic using that standard error. We will refer to our procedure as the “max-se” procedure because the case of just one restriction is by far the most common one. The max-se procedure has recently been studied in Davezies_2025, which cites an earlier version of this paper.
Henceforth, we denote the variance and standard error estimators based on $\hat{\bm{V}}_1^{(2)}$ and $\hat{\bm{V}}_1^{(3)}$ as CV$_{\kern -0.08333em1}^{\kern 0.08333em(2)}$ and CV$_{\kern -0.08333em1}^{\kern 0.08333em(3)}$ estimators, respectively, the ones based on $\hat{\bm{V}}_1^{(3+)}$ as CV$_{\kern -0.08333em1}^{\kern 0.08333em(3+)}$ estimators, and the ones implicit in \hyperref[{Wmin}]{\tagform@{\ref*{Wmin}}} as CV$_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$ estimators. In the scalar case, $\text{CV}_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$ is $\hat V_1^{(\max)} = \max \{\hat V_1^{(3)}, \hat V_G , \hat V_H \}$. This explains the “$(\max )$” superscript and also makes it clear that, asymptotically, the CV$_{\kern -0.08333em1}^{\kern 0.08333em(3)}$, CV$_{\kern -0.08333em1}^{\kern 0.08333em(3+)}$, and CV$_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$ estimators must be identical whenever the scores are positively correlated in either or both of the $G$ and $H$ dimensions.
In most cases where it makes sense to specify ${\bm{\Sigma}}$ as in \hyperref[{Sigma}]{\tagform@{\ref*{Sigma}}}, the CV$_{\kern -0.08333em1}^{\kern 0.08333em(3)}$, CV$_{\kern -0.08333em1}^{\kern 0.08333em(3+)}$, and CV$_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$ estimators will have exactly the same asymptotic properties. They may or may not be identical in practice. In fact, there are cases where they may differ greatly. This seems to be most common when there is very little intra-cluster correlation and/or the number of clusters is small, and/or the number of regressors is large, as we shall see in (ref).
The component CRVEs defined in \hyperref[{VhatG}]{\tagform@{\ref*{VhatG}}}, \hyperref[{VhatH}]{\tagform@{\ref*{VhatH}}}, and \hyperref[{VhatI}]{\tagform@{\ref*{VhatI}}} all have the form of the widely-used CV$_{\kern -0.08333em1}$ estimator. However, recent work by MNW-bootknife and Hansen-jack strongly suggests that, in the one\kern 0.04167em-way case, it is better to use a CRVE based on the cluster jackknife, which is analogous to the HC$_3$ estimator of MW_1985. The key idea of the cluster jackknife is to compute $G$ (or $H$ or $I$) sets of parameter estimates, each of which omits one cluster at a time, and then compute a CRVE using the variation among these estimates.
Let $J \in \{ G, H, I \}$, and let $j$ denote the corresponding lower-case letter. In the intersection dimension ($J=I$), the summation $\sum_{j=1}^J {\bm{Z}}_j$ should be interpreted as $\sum_{g=1}^G\sum_{h=1}^H {\bm{Z}}_{gh}$ for any ${\bm{Z}}$, where $I$ denotes the number of non-empty intersections, which may be smaller than $GH$. The OLS estimates of ${\bm\beta}$ when each cluster in the $J$ dimension is omitted in turn are
Then the component cluster-jackknife variance matrix estimators are
Thus the three\kern 0.04167em-term jackknife CRVE is
which is analogous to \hyperref[{3mat}]{\tagform@{\ref*{3mat}}}. The subscript “3” here follows the usual notation for jackknife variance matrices; see MNW-bootknife. There is also a two\kern 0.04167em-term jackknife CRVE and, more interestingly, one that is analogous to the CV$_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$ estimator. We refer to the three CRVEs based on the cluster jackknife as CV$_{\kern -0.08333em3}^{\kern 0.08333em(2)}$, CV$_{\kern -0.08333em3}^{\kern 0.08333em(3)}$, and CV$_{\kern -0.08333em3}^{\kern 0.08333em(\max)}$.
The CRVEs defined in \hyperref[{jackj}]{\tagform@{\ref*{jackj}}} are not the only cluster-jackknife variance matrix estimators. Instead of computing variances around $\hat{\bm\beta}$, one can instead compute them around (in the two\kern 0.04167em-way case) the three sample averages, $\bar{\bm\beta}^J = J^{-1}\sum_{j=1}^J \hat{\bm\beta}^{(j)}$. This makes the alternative CRVEs a little smaller than the ones given in \hyperref[{jackj}]{\tagform@{\ref*{jackj}}}. Because simulation experiments in BM_2002 and MNW-bootknife suggest that, in the one\kern 0.04167em-way case with $G$ clusters, inferences based on the alternative jackknife CRVE are almost identical to ones based on $\hat{\bm{V}}_G^{\rm JK}$\kern -0.08333em, we do not study the former in this paper.
Computing the component CRVEs in \hyperref[{jackj}]{\tagform@{\ref*{jackj}}} that are needed for CV$_{\kern -0.08333em3}^{\kern 0.08333em(3)}$, CV$_{\kern -0.08333em3}^{\kern 0.08333em(3+)}$\kern -0.08333em, and CV$_{\kern -0.08333em3}^{\kern 0.08333em(\max)}$ is somewhat more work than computing the ones in \hyperref[{VhatG}]{\tagform@{\ref*{VhatG}}}, \hyperref[{VhatH}]{\tagform@{\ref*{VhatH}}}, and \hyperref[{VhatI}]{\tagform@{\ref*{VhatI}}} that are needed for CV$_{\kern -0.08333em1}^{\kern 0.08333em(3)}$, CV$_{\kern -0.08333em1}^{\kern 0.08333em(3+)}$\kern -0.08333em, and CV$_{\kern -0.08333em1}^{\kern 0.08333em(\max)}$, especially when the number of non-empty intersections, $I$\kern -0.08333em, is large. The first thing is to calculate the cluster-level matrices and vectors
These quantities can be computed for the intersections with a single pass over the $N$ observations. The ones for the $G$ and $H$ dimensions are just summations of the ones for the appropriate intersections. The three sets of $\hat{\bm\beta}^{(j)}$ can then be computed using \hyperref[{delone}]{\tagform@{\ref*{delone}}} for the three clustering dimensions. Unfortunately, this may be expensive when both $k$ and $I$ are large, because computing the omit-\kern 0.04167em one\kern 0.04167em-cluster estimates for the intersections involves inverting $I$ different $k\times k$ matrices.
When computational cost is a concern, it can be reduced significantly by replacing $\hat{\bm{V}}_I^{\rm JK}$ in \hyperref[{3jack}]{\tagform@{\ref*{3jack}}} with $\hat{\bm{V}}_I$, yielding the mixed three-term estimator
Because $\hat{\bm{V}}_I$ is almost always smaller than $\hat{\bm{V}}_I^{\rm JK}$, $\hat{\bm{V}}_{3,1}^{(3)}$ will generally be larger than $\hat{\bm{V}}_3^{(3)}$. However, unless $I$ is small (which can only happen if both $G$ and $H$ are small or $I$ is much smaller than $GH$), the matrices $\hat{\bm{V}}_I$ and $\hat{\bm{V}}_I^{\rm JK}$ tend to be very similar. Thus the difference between \hyperref[{3jack}]{\tagform@{\ref*{3jack}}} and \hyperref[{31jack}]{\tagform@{\ref*{31jack}}} is negligible in most cases, as discussed near the end of (ref). However, it can be noticeable when there are many empty intersections; see (ref).
In many cases, the regression model \hyperref[{model}]{\tagform@{\ref*{model}}} will include fixed effects in the $G$ and $H$ dimensions; that is, two\kern 0.04167em-way fixed effects. If so, it may be rewritten as
Here the matrix ${\bm{Z}}$, which has $p$ columns, corresponds to the actual explanatory variables, and ${\bm\beta}_p$ contains the elements of ${\bm\beta}$ for those variables. The matrices ${\bm{D}}^G$ and ${\bm{D}}^H$ contain dummy variables for the fixed effects in dimensions $G$ and $H$, respectively. Collectively, these have $G+H-1$ columns, say $G$ for ${\bm{D}}^G$ and $H-1$ for ${\bm{D}}^H$\kern -0.08333em. Thus ${\bm{X}}=[{\bm{Z}}\;\;{\bm{D}}^G\;\;{\bm{D}}^H]$, and $k=p+G+H-1$.
For the model \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}}, there is an important computational issue. It is impossible to invert the matrices ${\bm{X}}^\top{\bm{X}} - {\bm{X}}_g^\top{\bm{X}}_g$ and ${\bm{X}}^\top{\bm{X}} - {\bm{X}}_h^\top{\bm{X}}_h$ in \hyperref[{subthings}]{\tagform@{\ref*{subthings}}}, because for each of them the row and column corresponding to the fixed effect for cluster $g$ or cluster $h$ contains only zeros. There are three ways to deal with this issue. The first is just to drop the subsamples in which the inversion is not possible. This is the default in many standard software routines, such as the prefix jackknife in Stata. However, it is not viable for \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}}, because every jackknife replication would have to be dropped. The second approach is to replace the inverse in \hyperref[{delone}]{\tagform@{\ref*{delone}}} by a generalized inverse. Then all of the coefficients except the fixed effect for the omitted cluster can be computed, and the latter is set to zero. Thus, whenever there are two\kern 0.04167em-way fixed effects, $\hat{\bm{V}}_3^{(3)}$ in \hyperref[{3jack}]{\textup{\tagform@{\ref*{3jack}}}} is effectively being defined as a $p\times p$ variance matrix for $\hat{\bm\beta}_p$ instead of a $k\times k$ variance matrix for $\hat{\bm\beta}$.
The third approach is to to partial out the cluster fixed effects before computing the one\kern 0.04167em-way CRVEs. However, this must be done with great care. It is valid to partial out cluster fixed effects in the $G$ dimension when calculating $\hat{\bm{V}}_G^{\rm JK}$\kern -0.08333em, but it is invalid to partial them out when calculating either $\hat{\bm{V}}_H^{\rm JK}$ or $\hat{\bm{V}}_I^{\rm JK}$\kern -0.08333em. The problem is that, after the cluster fixed effects in the $G$ dimension have been partialed out, the observations for every cluster in the $H$ and $I$ dimensions generally depend on observations in some or all of the other clusters in those dimensions. Thus $\hat{\bm\beta}^{(h)}$ and $\hat{\bm\beta}^{(i)}$ would not actually be vectors of omit-\kern 0.04167em one\kern 0.04167em-cluster estimates. Similarly, it is invalid to partial out fixed effects in the $H$ dimension when calculating either $\hat{\bm{V}}_G^{\rm JK}$ or $\hat{\bm{V}}_I^{\rm JK}$\kern -0.08333em. The $I$ dimension is always the most expensive one to deal with, because it involves the largest number of clusters, and it is not valid to partial out fixed effects in either the $G$ or $H$ dimensions when calculating $\hat{\bm{V}}_I^{\rm JK}$\kern -0.08333em. This makes it particularly attractive to use \hyperref[{31jack}]{\tagform@{\ref*{31jack}}} instead of \hyperref[{3jack}]{\tagform@{\ref*{3jack}}} when there are two\kern 0.04167em-way fixed effects.
It is conventional to employ the Student's $t$ distribution with $\min \{ G,H \} -1$ degrees of freedom to obtain $P$ values or critical values for $t$-statistics based on CV$_{\kern -0.08333em1}^{\kern 0.08333em(3)}$. As in the one\kern 0.04167em-way case, it seems reasonable to use the same distribution for $t$-statistics based on CV$_{\kern -0.08333em3}^{\kern 0.08333em(3)}$ as well, and this is the approach that we take.
However, at least two other methods could in principle be used. For one\kern 0.08333em-way clustering, BM_2002 proposes a way to obtain approximate critical values for $t$-tests based on CV$_{\kern -0.08333em1}$ by using a $t$ distribution with a calculated degrees\kern 0.04167em-of-freedom parameter; see also Imbens_2016. Another method based on the same idea is proposed in Hansen-jack,Hansen_2025. For the two\kern 0.04167em-way case, one could in principle use the same sort of approximate critical value. However, we are not aware of any method for obtaining such a critical value for $t$-statistics based on two\kern 0.04167em-way clustering. This is an area for future research.
Another possibility is to use bootstrap methods. The pigeonhole bootstrap of Owen_2007 was studied in Menzel_2021 and found to be conservative in general. That paper also proposed some new and rather complicated bootstrap procedures for inference about the sample mean. The wild cluster bootstrap \citep*{CGM_2008,DMN_2019} has been widely used for inference with one\kern 0.04167em-way clustering, and MNW_2021 suggested using it for two\kern 0.04167em-way clustering as well. In that paper, the usual wild cluster bootstrap for one of the $G$, $H$, or $I$ dimensions is used to generate the bootstrap samples. This procedure is not entirely satisfactory, because the bootstrap samples cannot reproduce the intra-cluster covariances among the residuals. Nevertheless, this wild bootstrap routine is conveniently and efficiently coded in the boottest package in Stata using the bootclust option; see \citet*{RMNW} for details. Recently, Hounyo-boot proposes a wild bootstrap DGP that gives positive weight to both dimensions.
In the absence of any satisfactory alternative, we currently recommend using the cluster jackknife together with critical values based on the Student's $t$ distribution with $\min\{ G,H \} -1$ degrees of freedom. As we shall see in (ref), this approach often works remarkably well. Whether combining the jackknife with a bootstrap procedure would perform even better is a topic for future research; see MNW-bootknife for evidence on this with one\kern 0.04167em-way clustering.
Computing the three\kern 0.04167em-term cluster-jackknife estimator for the two\kern 0.04167em-way fixed-effects model \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}} can be costly when $G$ and $H$ are not fairly small. The cost of forming the ${\bm{X}}_j^\top{\bm{X}}_j$ matrices and the ${\bm{X}}_j^\top{\bm{y}}_j$ vectors is roughly $O(Nk^2) = O(N(G+H+p-1)^2)$, because ${\bm{X}}$ has $k = p+G+H-1$ columns. Since \hyperref[{delone}]{\tagform@{\ref*{delone}}} has to be computed $G+H+I \approx G+H+GH$ times, the cost of computing the cluster-jackknife estimates after the ${\bm{X}}_j^\top{\bm{X}}_j$ matrices and ${\bm{X}}_j^\top{\bm{y}}_j$ vectors have been formed is roughly $O(GHk^2)=O(GH(G+H+p-1)^2)=O(G^4)$ if $G\approx H$\kern -0.08333em.
Most of the computational cost of the two\kern 0.04167em-way cluster jackknife arises from the need to deal with the $I \le GH$ intersections. When $I<\!<GH$, the cost can be greatly reduced if the empty intersections are skipped when calculating the omit-\kern 0.04167em one\kern 0.04167em-cluster estimates using \hyperref[{delone}]{\tagform@{\ref*{delone}}}. An additional reduction is possible by using \hyperref[{31jack}]{\tagform@{\ref*{31jack}}} instead of \hyperref[{3jack}]{\tagform@{\ref*{3jack}}}.
Properties of classic-jackknife variance estimators are well known. However, for the cluster jackknife, the only analysis of theoretical properties that we are aware of is in Hansen-jack. In the context of the linear regression model with one\kern 0.04167em-way clustering, it shows that a certain cluster-jackknife variance estimator (which is not quite the same as $\hat{\bm{V}}_3$, but should usually be very similar) is never downward biased. Moreover, the associated $t$-tests and confidence intervals have worst-case size, or coverage, that is controlled by the Cauchy distribution. In contrast, variance estimators based on CV$_{\kern -0.08333em1}$ can be severely downward biased, the associated $t$-tests have worst-case size of 1, and the associated confidence intervals have worst-case coverage of 0.
In this subsection, we prove consistency of the two\kern 0.04167em-way cluster jackknife CRVE. We will need the following two assumptions.
(ref) guarantees existence of the omit-\kern 0.04167em one\kern 0.04167em-cluster estimators in \hyperref[{delone}]{\tagform@{\ref*{delone}}}, and hence the cluster jackknife. Note that this rules out cluster fixed effects and other cases where a regressor is non-zero for only one cluster (in any dimension). Practical ways to deal with this situation were discussed in the previous section. Existing proofs for the validity of other two\kern 0.04167em-way CRVEs, such as Davezies_2025 and yap_2025, also rule out cluster fixed effects.
(ref) is identical to Assumption 3 in yap_2025, but adapted to our notation. Part (a) is a standard moment condition. Parts (b) and (c) restrict the heterogeneity of cluster sizes and rule out degenerate non-Gaussian cases analyzed in Menzel_2021. They can be viewed as generalizations of Assumptions 2 and 3 in DMN_2019 to allow two\kern 0.04167em-way clustering. Part (d) strengthens \hyperref[{def Omega}]{\tagform@{\ref*{def Omega}}} to independence, and part (e) is a version of the usual rank condition for OLS. See yap_2025 for a detailed discussion. An alternative asymptotic framework is that of separately exchangeable arrays Davezies_2018, under which our results could also be proven by the same arguments yap_2025.
A proof of (ref) is given in (ref). Under (ref), yap_2025 shows that $( \operatorname{Var} ( \hat{\bm\beta} ))^{-1/2}(\hat{\bm\beta} - {\bm\beta}_0 ) \overset{d} \longrightarrow {\rm N}({\bm{0}} , {\bf I}_k)$. Combined with our (ref)(i), this implies that $t$-tests or $F$-tests based on the three jackknife CRVEs, CV$_3^{(3)}$, CV$_3^{(3+)}$, and CV$_{3,1}^{(3)}$, have correct asymptotic size, and associated confidence intervals have correct asymptotic coverage.
The results in (ref)(ii),(iii) require additional conditions to obtain asymptotically valid inference. For part (ii), we have stated the result for linear contrasts. Sufficient primitive conditions for this part were studied in Davezies_2025, and under such conditions it also holds by (ref)(ii) that inference based on the jackknife CV$_3^{(\max)}$ is asymptotically valid, at least for linear contrasts. For part (iii), ${\bm{V}}_I$ is asymptotically negligible when there is sufficient intra-cluster correlation in either of the two main clustering dimensions, $G$ or $H$ Davezies_2018,MNW_2021, and then (ref)(iii) shows that inference based on the jackknife CV$_3^{(2)}$ is asymptotically valid. On the other hand, for example, if there is clustering only at the intersection level, then ${\bm{V}}_I$ is not asymptotically negligible, and $\hat{\bm{V}}_1^{(2)}$ is not consistent; see our discussion around \hyperref[{2mat crve}]{\tagform@{\ref*{2mat crve}}}--\hyperref[{V3hatx}]{\tagform@{\ref*{V3hatx}}}. In that case, $\hat{\bm{V}}_3^{(2)}$ is not consistent either.
Simulation results MNW-guide,MNW-bootknife, MNW-influence,Hansen-jack,Hansen_2025 have shown that one\kern 0.04167em-way cluster-jackknife CRVEs perform better in finite samples than conventional CRVEs in terms of coverage of confidence intervals and size of tests. In this subsection, we discuss possible reasons underlying those results in the context of two\kern 0.04167em-way clustering. The key reason seems to be that cluster-jackknife CRVEs handle cluster size variation, and heterogeneity more generally, better than do conventional CRVEs. This is particularly important for three\kern 0.04167em-term estimators, as we explain.
It is known that CV$_{\kern -0.08333em3}$ estimators are less (downward) biased than CV$_{\kern -0.08333em1}$ ones Efron-Stein,Hansen-jack. The reason for this can be seen intuitively as follows. The one\kern 0.04167em-way cluster-jackknife CRVEs in \hyperref[{jackj}]{\tagform@{\ref*{jackj}}} can be rewritten as
where the modified score vectors $\ddot{\bm{s}}_j$ are defined as
and ${\bm{M}}_{jj}$ denotes the \th{(j,j)} block of ${\bm{M}}_{{\bm{X}}} = {\bf I}_N - {\bm{X}} ({\bm{X}}^\top{\bm{X}} )^{-1}{\bm{X}}^\top$. For a proof of equality of \hyperref[{jackj}]{\tagform@{\ref*{jackj}}} and \hyperref[{jackjMgg}]{\tagform@{\ref*{jackjMgg}}}, see MNW-bootknife. Normalizing the modified score vectors in \hyperref[{ddotbis}]{\tagform@{\ref*{ddotbis}}} by the factor ${\bm{M}}_{jj}^{-1}$ undoes some of the shrinkage caused by least squares. Since the ${\bm{M}}_{jj}$ are inversely related to cluster leverage MNW-influence, the cluster-jackknife CRVE puts more weight on clusters with high leverage compared with the CV$_1$ estimator. This accounts for the smaller bias of the former relative to the latter, because high-leverage clusters are relatively more important in determining the actual variance of the estimator.
In many two\kern 0.04167em-way designs, clusters vary greatly in size and/or leverage in one or both dimensions, or there are few clusters in one dimension. Thus, one or both of $\hat{\bm{V}}_G$ and $\hat{\bm{V}}_H$ is likely to be seriously downward biased. However, $\hat{\bm{V}}_I$ is usually based on a much larger number of clusters, often as many as $G \times H$. As a result, its downward bias is likely to be comparatively moderate. In consequence, when $\hat{\bm{V}}_I$ is subtracted from the sum of $\hat{\bm{V}}_G$ and $\hat{\bm{V}}_H$ to form $\hat{\bm{V}}_{\kern -0.08333em1}^{(3)}$\kern -0.08333em, there is a good chance that the latter will be very severely biased. In contrast, the arguments above and results in Hansen-jack (for one\kern 0.04167em-way clustering) suggest that $\hat{\bm{V}}_G^{\rm JK}$ and $\hat{\bm{V}}_H^{\rm JK}$ are never downward biased, although either or both may be upward biased. It is possible that $\hat{\bm{V}}_I^{\rm JK}$ may be upward biased in this case, but since it is normally based on a much larger number of clusters, any such bias is likely to be modest, and subtracting it is not likely to cause much downward bias in $\hat{\bm{V}}_{\kern -0.08333em3}^{(3)}$ itself. These arguments suggest that $\hat{\bm{V}}_{\kern -0.08333em3}^{(3)}$ is more likely to be positive definite than $\hat{\bm{V}}_{\kern -0.08333em1}^{(3)}$ and that tests based on $\hat{\bm{V}}_{\kern -0.08333em3}^{(3)}$ should be more reliable than ones based on $\hat{\bm{V}}_{\kern -0.08333em1}^{(3)}$.
The above arguments imply that, if the sample is heterogeneous in only one dimension, so that only one of $\hat{\bm{V}}_G$ and $\hat{\bm{V}}_H$ is severely downward biased, then the downward bias in $\hat{\bm{V}}_{\kern -0.08333em1}^{(3)}$ is likely to be relatively moderate. This case probably occurs quite often in panel settings, where samples (and cluster sizes) are often heterogeneous across cross-sectional units but homogeneous across time periods. For instance, a commonly used dataset like the Current Population Survey (CPS) will be unbalanced in terms of the number of observations per state, but strongly balanced in terms of the number of observations per year. In such cases, we would still expect CV$_{\kern -0.08333em3}$-based estimators to be more accurate than CV$_{\kern -0.08333em1}$-based ones, but probably by a smaller margin than in cases with double heterogeneity.
In empirical research, it is very commonly found that some intersections of the two clustering dimensions contain no observations. The possibility of empty intersections can be important for two\kern 0.04167em-way clustering, but it cannot arise for one\kern 0.04167em-way clustering. To examine the importance of empty intersections, consider two hypothetical samples, each with $G=H=10$. Thus there are $100$ intersections. For one sample, no intersections are empty, but 70 of them contain just 1 observation. For the other sample, there are 70 empty intersections. Now consider the cluster-jackknife estimator, $\hat{\bm{V}}_I^{\rm JK}$. In the first sample, it is based on 100 terms. Since dropping just one observation should not change $\hat{\bm\beta}^{(i)}$ very much, the terms in the summation in \hyperref[{jackj}]{\tagform@{\ref*{jackj}}} corresponding to the tiny intersections must all be very small. In the second sample, the cluster-jackknife estimate is based on just 30 terms. The terms that were small in the first sample have vanished, which seems to be a small difference. The only other difference between the two samples is that the leading factor in $\hat{\bm{V}}_I^{\rm JK}$ will be 99/100 in the first sample and 29/30 in the second, which seems inconsequential. Thus the cluster-jackknife estimator handles empty intersections in a reasonable fashion.
Most of our experiments deal with the two\kern 0.04167em-way fixed-effects model \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}}. The number of coefficients is $k=p+G+H-1$, but we focus on tests of a single coefficient, say $\beta_1$. Although \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}} is very widely used, many existing simulation experiments for two\kern 0.04167em-way clustering do not include cluster fixed effects. This is probably because, when the intra-cluster correlations are generated by a random-effects model, cluster fixed effects absorb all of them. For example, the experiments in CGM_2011 and MNW_2021 do not include fixed effects. In contrast, the placebo\kern 0.04167em-regression experiments in Section 3.2 of the former paper use actual data instead of a random-effects model, and they do include two\kern 0.04167em-way fixed effects.
In order to generate data for the model \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}}, the disturbances must be generated in a way that allows for two\kern 0.04167em-way intra-cluster correlation that is not removed by cluster fixed effects. We use factor models of the form
Here $\xi^1_g$ and $\xi^2_g$ are random effects, distributed as ${\rm N}(0,1)$, which apply respectively to the odd-numbered and even-numbered observations within the \th{g} cluster in the $G$ dimension. Similarly, $\xi^1_h$ and $\xi^2_h$ are ${\rm N}(0,1)$ random effects which apply to the odd-numbered and even-numbered observations within the \th{h} cluster in the $H$ dimension. The $\zeta_{ghi}$ are independent standard normals.
The values of $\sigma_g$, $\sigma_h$, and $\sigma_\epsilon$ determine the amount of correlation for the odd-numbered and even-numbered observations within each cluster, and hence the correlations within and across the clusters in the $G$, $H$, and $I$ dimensions. There will be no correlation for observations that belong to different clusters in the $G$ and $H$ dimensions. Specifically, the intra-cluster correlations are $\rho_g$ and $\rho_h$, with $\rho_j=\sigma_j^2$ for $j=g,h$. To ensure that the $z_{ghi}$ have variance unity, the value of $\sigma_\epsilon^2 = 1 - \sigma_g^2 - \sigma_h^2$. This constrains $\rho_g$ and $\rho_h$ not to be too large.
The factor model \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} provides a simple way to generate data for a model with two\kern 0.04167em-way fixed effects. It is based on a one\kern 0.04167em-way DGP used in \citet*{MNW-testing} and can be interpreted in a variety of ways, depending on the nature of the data. The idea is that there are two types of observations within each cluster in each dimension, and all the intra-cluster correlation is within each type. For example, with clustering at the geographical level, there might be two sub\kern 0.04167em-regions. With clustering at the industry level, there might be two types of firm. The key assumption is that the researcher knows which cluster an observation belongs to in each dimension, but not which type. Including cluster fixed effects explains some of the intra-cluster correlation by estimating averages of $\xi^1_g$ and $\xi^2_g$ for each $G$ cluster and averages of $\xi^1_h$ and $\xi^2_h$ for each $H$ cluster, but it does not explain all of it. Thus cluster-robust inference is still needed.
In several of the experiments, we focus on cluster size variation. Following MW-JAE and DMN_2019, the cluster sizes in the $G$ dimension are given by
where $[x]$ denotes the integer part of $x$. The value of $N_G$ is then set to $N - \sum_{g=1}^{G-1} N_g$. The formula \hyperref[{gameq}]{\tagform@{\ref*{gameq}}}, perhaps with a different value of $\gamma$, is also used in the $H$ dimension. Assuming that the distributions are independent, $N_{gh} \approx N_g N_h / N$. In a final step, the cluster sizes are adjusted to ensure that they are all integers with $N = \sum_{g=1}^G N_g = \sum_{h=1}^H N_h = \sum_{g=1}^G\sum_{h=1}^H N_{gh}$.
The way in which the regressors are generated inevitably affects the finite\kern 0.04167em-sample properties of every cluster-robust test statistic. The differences between asymptotic and finite\kern 0.04167em-sample distributions arise mainly from the discrepancies between the disturbance vector ${\bm{u}}$ and the residual vector $\hat{\bm{u}} = {\bm{M}}_{\bm{X}}{\bm{u}}$. We use \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}} to generate the regressor matrix ${\bm{Z}}$ in \hyperref[{TWFE}]{\tagform@{\ref*{TWFE}}} as well as the disturbances. In most experiments, we set $\rho_g^x=\rho_h^x=0.2$ for the regressors and $\rho_g=\rho_h=0.1$ for the disturbances. We use these values because, in practice, regressors often display more intra-cluster correlation than residuals. For this base case, we deliberately avoid situations, to be discussed below, in which the amount of intra-cluster correlation is very small.
Throughout, when a standard error was either non-positive or undefined for a particular replication, which can happen for CV$_{\kern -0.08333em1}^{(3)}$ or CV$_{\kern -0.08333em3}^{(3)}$, we counted it as a rejection.
The first set of experiments focuses on cluster size variation, determined by the parameter $\gamma$ in \hyperref[{gameq}]{\tagform@{\ref*{gameq}}}. (ref) shows rejection frequencies for eight different $t$-tests from 18 experiments with $G=15$, $H=12$, and $N=9,\kern -0.08333em000$. In Panel (a), the value of $\gamma$ is varied simultaneously from 0.0 to 4.0 in both dimensions. In Panel (b), $\gamma=0$ for the $H$ dimension, and $\gamma$ varies from 0.0 to 4.0 for the $G$ dimension. In these experiments, the number of regressors is $p=10$. All tests would have performed better if this had been a smaller number. The effects of varying $p$ will be investigated below.
In Panel (a), when all clusters are the same size (the leftmost point on the horizontal axes), $t$-tests based on the classic CV$_{\kern -0.08333em1}^{(3)}$ variance matrix estimator over-reject noticeably, as do those based on CV$_{\kern -0.08333em1}^{(\max)}$. Rejection frequencies are considerably lower for CV$_{\kern -0.08333em1}^{(3+)}$, and lower still for CV$_{\kern -0.08333em1}^{(2)}$. In contrast, $t$-tests based on the CV$_{\kern -0.08333em3}^{(3)}$ and CV$_{\kern -0.08333em3}^{(\max)}$ estimators are very close to nominal size, while those based on CV$_{\kern -0.08333em3}^{(3+)}$ and CV$_{\kern -0.08333em3}^{(2)}$ under-reject substantially. As the value of $\gamma$ increases, all the CV$_{\kern -0.08333em1}$ rejection frequencies rise sharply, while those for the CV$_{\kern -0.08333em3}$ tests hardly change. Using the max-se procedure has almost no effect when cluster sizes vary little, but it modestly reduces rejection frequencies for CV$_{\kern -0.08333em1}$ tests when they vary a lot.
In Panel (b), the overall patterns are similar. However, as predicted in (ref), rejection frequencies increase less rapidly when $\gamma$ just increases in the $G$ dimension than when it increases in both dimensions. In both panels, as must be the case, $t$-tests based on two\kern 0.04167em-term variance estimators always reject less often than $t$-tests based on three\kern 0.04167em-term ones. This is a good thing for the CV$_{\kern -0.08333em1}$ tests, but not for the CV$_{\kern -0.08333em3}$ ones.
One possibly surprising feature of (ref) is how much the 3+ tests based on the eigen-decomposition differ from the ordinary three\kern 0.04167em-term tests. This happens because, with 36 coefficients to estimate (26 of them fixed effects), the three\kern 0.04167em-term variance matrices are not positive definite, leading to negative eigenvalues in many cases. Whether or not the 3+ tests differ substantially from the ordinary three\kern 0.04167em-term tests depends on how the coefficient being tested loads on the problematic eigen-directions. For this design, the loading seems to be quite large. For CV$_{\kern -0.08333em1}$, the 3+ variant performs better than the usual three\kern 0.04167em-term test, but for CV$_{\kern -0.08333em3}$, it performs worse, under-rejecting about as much as the two\kern 0.04167em-term test.
Except for quite small values of $\gamma$, the intersections in these experiments vary greatly in size. For example, when both values of $\gamma$ equal 2, which is the base case for many of our subsequent experiments, the smallest intersection contains 6 observations, and the largest contains 253. The sizes of the intersections vary much more than those of the $G$ clusters, which range from 223 to 1443, or the $H$ clusters, which range from 282 to 1769. Although these numbers depend on the way in which we generate cluster sizes, it is inevitable that, when the cluster sizes vary in both dimensions, the sizes of the intersections vary more dramatically.
(ref) does not report results for the mixed variance estimator \hyperref[{31jack}]{\tagform@{\ref*{31jack}}}, or its max-se or 3+ variants, because the results are very similar to the corresponding cluster-jackknife variants. For smaller values of $\gamma$, the differences are negligible. The largest differences occur in Panel (a) when $\gamma=4$. In that case, the test based on $\hat{\bm{V}}_3^{(3)}$ rejects 6.12% of the time, while the one based on $\hat{\bm{V}}_{3,1}^{(3)}$ rejects 5.55%. The max-se tests differ much less, rejecting 5.35% and 5.08%, respectively, and the 3+ tests barely differ. In most of our simulations, except when there were a great many empty intersections ((ref)), the differences were much smaller than in this case.
As MNW_2021 shows, test statistics based on the two\kern 0.04167em-term variance estimator are asymptotically too small whenever the scores are asymptotically uncorrelated beyond the intersection level. This suggests that they are likely to under-reject severely when the amount of intra-cluster correlation in either the disturbances or regressors is very small. In (ref), we vary both values of $\rho$ for the disturbances from 0.000 to 0.200. For clarity, the horizontal axis uses a square root transformation. The numbers of clusters, observations, and regressors are larger in Panel (b) than in Panel (a); see the notes to the figure.
Several results stand out in (ref). For small values of $\rho$, the rejection frequencies of the three\kern 0.04167em-term tests are much higher than those of the corresponding max-se tests. This is particularly true for the CV$_{\kern -0.08333em1}$ tests. For the smallest values of $\rho$, the two\kern 0.04167em-term tests under-reject to an extreme extent, as the theory in MNW_2021 predicts. Interestingly, so do the eigen-decomposition tests. In fact, in both panels, the two-term and 3+ tests perform very similarly. In both panels and for all values of $\rho$, the CV$_{\kern -0.08333em3}$-based tests reject less than the corresponding CV$_{\kern -0.08333em1}$-based tests. Except for the smallest values of $\rho$, the CV$_{\kern -0.08333em3}^{(3)}$ and CV$_{\kern -0.08333em3}^{(\max)}$ tests are very similar in Panel (a) and identical in Panel (b), and they perform very well.
For the smallest values of $\rho$ in these experiments, there were a number of replications for which the three\kern 0.04167em-term variance of $\hat\beta_1$ was negative. This happened more often for $G=15$ than for $G=30$, and more often for CV$_{\kern -0.08333em1}$ than for CV$_{\kern -0.08333em3}$. Since we could not calculate the $t$-statistic for these replications, we classified them as rejections. In the most extreme case, when $\rho=0.000$ for $G=15$ ($G=30$), this happened 2.60% (0.24%) of the time for CV$_{\kern -0.08333em1}^{(3)}$ and 2.1% (0.16%) for CV$_{\kern -0.08333em3}^{(3)}$. These numbers declined sharply as the value of $\rho$ increased.
It is not only the correlations of the disturbances that matter. In (ref), we vary the correlations of the regressors, either in both dimensions, in Panel (a), or just in the $H$ dimension, in Panel (b). Although they may seem small, the largest values of the $\rho^x$ parameters here are not far short of the largest possible values; see the discussion below \hyperref[{facDGP}]{\tagform@{\ref*{facDGP}}}. The horizontal axis does not use a square-root scale as (ref) did, because the dependence on $\rho^x$ for small values is not as extreme as the dependence on $\rho$ in that figure.
It is clear from (ref) that the way in which the regressors are distributed can have substantial effects on rejection frequencies. Every test except CV$_{\kern -0.08333em1}^{(3)}$ and CV$_{\kern -0.08333em3}^{(3)}$ can either over-reject or under-reject, depending on the values of the two $\rho^x$ parameters. The two three\kern 0.04167em-term tests always over-reject, although only very slightly for CV$_{\kern -0.08333em3}^{(3)}$ in Panel (a) for $\rho^x>0.05$. The most reliable tests are the ones based on CV$_{\kern -0.08333em3}^{(3)}$ and, especially, CV$_{\kern -0.08333em3}^{(\max)}$. This is particularly the case for larger values of the $\rho^x$ parameters, where all the CV$_{\kern -0.08333em1}$-based tests over-reject substantially. Panels (a) and (b) are quite similar when both $\rho^x$ parameters, or just $\rho^x_h$, are large, but the two panels differ substantially when the intra-cluster correlations are small.
The number of regressors inevitably matters. This was analyzed in the context of heteroskedasticity by \citet*{CJN_2018}. In (ref), $p$ varies from 1 to 16. In Panels (a) and (c), $G=15$, $H=12$, and $N=9,\kern -0.08333em000$. In Panels (b) and (d), $G=30$, $H=24$, and $N=36,\kern -0.08333em000$. As $p$ increases, the rejection rates for the CV$_{\kern -0.08333em1}$ tests increase, but those for the CV$_{\kern -0.08333em3}$ tests decrease slightly. In all panels, the max-se tests perform nearly the same as the three\kern 0.04167em-term tests. Throughout (ref), the CV$_{\kern -0.08333em3}^{(3)}$ and CV$_{\kern -0.08333em3}^{(\max)}$ tests perform very well.
In the lower two panels, the 26 or 53 fixed effects are replaced by a constant term. Rejection frequencies for the CV$_{\kern -0.08333em1}$ tests still increase with $p$, but more slowly, while those for the CV$_{\kern -0.08333em3}$ tests still decrease, at about the same slow rate. In these two panels, the 3+ tests are nearly identical to the ordinary three\kern 0.04167em-term tests. These are the only experiments in which we omit the fixed effects. Their presence evidently has a large impact on the performance of the 3+ tests but a fairly modest effect on that of the other tests.
(ref) suggests that rejection frequencies for all the CV$_{\kern -0.08333em1}$ tests increase fairly rapidly with $p$, the number of regressors that are not fixed effects, while those for all the CV$_{\kern -0.08333em3}$ tests decrease quite slowly. We conjecture that this is happening because all the regressors are correlated within clusters in one or both dimensions. Thus, as the number of regressors increases, more and more of the intra-cluster correlation in the disturbances is explained by the regressors, so that less of it remains in the residuals.
In order to investigate this conjecture, we modify the way in which we generate the regressors. The first one (the test regressor) is generated as before, but then we generate an additional $q$ binary regressors which equal 0 or 1 with probability 0.5. In one set of simulations, they are completely independent across observations. In a second set, they are generated at the intersection level, identical within each intersection and independent across intersections.
(ref) shows rejection frequencies as a function of $q$, which varies from 0 to 30. In both panels, the test based on CV$_{\kern -0.08333em3}^{(\max)}$ performs best, over-rejecting slightly for all values of $q$. The tests based on CV$_{\kern -0.08333em3}^{(3)}$, CV$_{\kern -0.08333em1}^{(3+)}$, and CV$_{\kern -0.08333em1}^{(2)}$ perform nearly as well, with the former two over-rejecting slightly and the latter under-rejecting slightly. The value of $q$ has very little effect on most of the tests.
The differences between (ref) and (ref) are striking. In the former, all the regressors are correlated within both the $G$ and $H$ clusters. We saw there that adding more regressors with this property can substantially increase rejection frequencies for CV$_{\kern -0.08333em1}$ tests and slightly decrease them for CV$_{\kern -0.08333em3}$ tests. In contrast, adding more regressors that are uncorrelated across observations or across intersections has almost no effect on rejection frequencies. Empirical applications of two\kern 0.04167em-way clustering often involve many controls. Whether or not these controls exhibit substantial correlation in either dimension can evidently be important.
All the results can be expected to improve as the number of clusters increases in either or both dimensions. We saw this in (ref). Therefore, in (ref), $G$ varies from 10 to 45 by 5, $H$ is always equal to $4G/5$, and $N$ is proportional to $GH$, so that the sizes of the intersections are roughly constant. In Panel (a), $p=5$. Here, CV$_{\kern -0.08333em1}^{(3)}$ and CV$_{\kern -0.08333em1}^{(\max)}$, and to a lesser extent CV$_{\kern -0.08333em3}^{(3+)}$, always over-reject, but they improve steadily as $G$ (and $H$) increase. CV$_{\kern -0.08333em1}^{(\max)}$ over-rejects less severely than CV$_{\kern -0.08333em1}^{(3)}$ for the smallest values of $G$, but the former is almost indistinguishable from the latter for $G\ge15$. In contrast, CV$_{\kern -0.08333em3}^{(3)}$ always works almost perfectly, with CV$_{\kern -0.08333em3}^{(\max)}$ yielding virtually identical results for $G\ge15$. By what seems to be coincidence, CV$_{\kern -0.08333em1}^{(2)}$ also works well.
In Panel (b) of (ref), $p$ is increased to 15. The CV$_{\kern -0.08333em1}$-based tests now over-reject much more severely, but tests based on CV$_{\kern -0.08333em3}^{(3)}$ and CV$_{\kern -0.08333em3}^{(\max)}$ perform extremely well. In contrast, tests based on CV$_{\kern -0.08333em3}^{(2)}$ and CV$_{\kern -0.08333em3}^{(3+)}$ are almost identical and always under-reject. Clearly, omitting the intersection term or using the eigen-decomposition is helpful for CV$_{\kern -0.08333em1}$, because the three\kern 0.04167em-term tests over-reject, but harmful for CV$_{\kern -0.08333em3}$, because the three\kern 0.04167em-term cluster-jackknife tests are approximately sized correctly.
Up to this point, the data for all of our experiments have been generated in such a way that $I=GH$. In other words, there have been no datasets with empty intersections. But empirical examples with two\kern 0.04167em-way clustering often involve empty intersections. In the next set of experiments, we therefore change the DGP so that intersections can be empty. The details are somewhat complicated and are therefore omitted. There is a parameter that determines the fraction of the intersection clusters, starting with the smallest ones, to be made empty by reallocating their observations proportionally to other clusters. We perform two sets of experiments. In both of them, the 15 clusters in the $G$ dimension are generated from \hyperref[{gameq}]{\tagform@{\ref*{gameq}}} with $\gamma=2$. In the first set $H=10$, and in the second set $H=20$. The maximum observed number of empty intersections is 95 (out of 150) in the first set and 270 (out of 300) in the second set.
(ref) shows rejection frequencies as functions of the fraction of empty intersections. This fraction evidently matters, especially for the CV$_{\kern -0.08333em1}$-based tests, although not dramatically so in these experiments. As usual, $t$-tests based on CV$_{\kern -0.08333em3}^{(\max)}$ always perform best, and in fact they perform extremely well. Some of the other tests perform quite poorly. As a rule, tests that over-reject or under-reject when there are no empty intersections do the same thing to a greater extent when there are many empty intersections.
When the number of empty intersections is large, $I$ is not a great deal larger than $G$ or $H$, and the mixed variance estimators based on \hyperref[{31jack}]{\tagform@{\ref*{31jack}}} are no longer almost the same as the CV$_{\kern -0.08333em 3}$ ones. For the most extreme case, which is the rightmost point in Panel (b) of (ref), $I=30$ and $H=20$. In this case, the test based on $\hat{\bm{V}}_3^{(3)}$ rejects 6.20% of the time, while the one based on $\hat{\bm{V}}_{3,1}^{(3)}$ rejects 4.82%. For the max-se tests, the corresponding values are 4.53% and 4.08%.
Since some tests tend to over-reject and others tend to under-reject under the null hypothesis, it is inevitable that the former will appear to have more power than the latter. (ref) shows power functions for all eight tests for a particular case. The functions never cross, so there is nothing surprising here. For every alternative, the ranking of the tests by power is identical to their ranking by rejection frequencies under the null hypothesis. Thus the fact that all of the CV$_{\kern -0.08333em1}$ tests appear to be more powerful than any of the CV$_{\kern -0.08333em3}$ tests simply occurs because the former are over-sized under the null. In this experiment, the power functions for CV$_1^{\kern 0.08333em(\max)}$ and CV$_3^{\kern 0.08333em(\max)}$ are indistinguishable from those for CV$_1^{\kern 0.08333em(3)}$ and CV$_3^{\kern 0.08333em(3)}$, respectively. All the tests evidently reject with probability one when $\beta_1$ is sufficiently large.
We did not explicitly study coverage of confidence intervals, because inverting tests that over-reject must yield intervals that under-cover, and vice versa. Similarly, inverting tests with higher power will produce shorter intervals than inverting tests with lower power. In these experiments, only the tests based on CV$_3^{\kern 0.08333em(3)}$ and CV$_3^{\kern 0.08333em(\max)}$ have approximately correct rejection frequencies, so only intervals based on them will have approximately correct coverage. Although it would be possible to obtain shorter intervals by using other standard errors, any such intervals would be misleadingly short.
In this section, we study two empirical examples. These involve different types of two\kern 0.04167em-way clustering and have clusters that behave in different ways across the two clustering dimensions.
In a fascinating paper, Alsan_2015 studies the effects of the tsetse fly on African development. The key explanatory variable is the “tsetse suitability index,” or TSI, which measures the extent to which climate (temperature and humidity) is suitable for the tsetse fly to thrive. There are seven dependent variables, which measure various aspects of economic and political development. Each of these is regressed on the TSI, whose coefficient is denoted $\beta$, and on eleven other variables in the columns labeled “(4)” in Table 1 and “(8)” in Table 3 of Alsan_2015. The former uses one\kern 0.04167em-way clustering by “cultural province” and the latter uses two\kern 0.04167em-way clustering by cultural province and country. There are 44 countries and either 43 or 44 cultural provinces, depending on the regressand. Since the total number of observations varies between 315 and 485, most clusters are quite small, and there are many empty intersections. The number of non-empty intersections varies between 112 and 142.
In (ref), we report $P$ values based on sixteen different standard error estimates, eight using conventional standard errors (Panel A) and eight using jackknife ones (Panel B). Because the ordinary three\kern 0.04167em-term and eigen-decomposition three\kern 0.04167em-term standard errors are identical in all cases (to the number of digits reported), we only report the former.
Even though the clusters are quite small (the largest is 63, which is for clustering by country when the dependent variable is the log of population density), the way in which we cluster often makes a substantial difference. Not clustering at all sometimes leads to extremely small $P$ values, as does one\kern 0.04167em-way clustering by intersection. Clustering in two dimensions often, but not always, leads to larger $P$ values than clustering in just one dimension. For two\kern 0.04167em-way clustering, the cluster-jackknife $P$ values are never smaller than the conventional ones, and they are mostly considerably larger. We also calculated $P$ values based on the mixed three-term estimator \hyperref[{31jack}]{\tagform@{\ref*{31jack}}}, but we do not report them because, with so many empty intersections, it is hard to justify using the mixed estimator. They are very similar to the ones in Panel B, although slightly larger in all but one case.
Panel C of (ref) presents a number of the summary statistics calculated by twowayjack for this example. Specifically, it presents coefficients of variation for the partial leverages of the TSI variable and for the $\hat\beta^{(g)}$ for clustering by cultural province, country, and intersection. It also displays the effective number of clusters $G^*=G^*(0)$ for the two primary dimensions CSS_2017,MNW-influence. These diagnostics can help to explain why some of the $P$ values in Panels A and B differ by more than others. The notable $P$ value differences between CV$_{\kern -0.08333em3}^{(\max)}$ and CV$_{\kern -0.08333em1}^{(\max)}$ occur for `plow use,' `indigenous slavery,' and `centralization.' For these three variables, we see the largest coefficients of variation for the omit-\kern 0.04167em one\kern 0.04167em-country estimates, and, to a slightly lesser extent, for the omit-\kern 0.04167em one\kern 0.04167em-culture ones.
The results of (ref) suggest that CV$_{\kern -0.08333em 3}^{(\max)}$ yields the most reliable $P$ values. The CV$_{\kern -0.08333em 3}^{(\max)}$ $P$ value for TSI is less than 0.05 for four of the seven dependent variables. In contrast, the CV$_{\kern -0.08333em1}$ $P$ values for one\kern 0.04167em-way clustering by cultural province used in Alsan_2015 are less than 0.05 for all seven variables in Table 1 (4), and the ones for two\kern 0.04167em-way CV$_{\kern -0.08333em1}^{(3)}$ clustering are less than 0.05 for six of them in Table 3 (8). Thus, although there is still a good deal of evidence that the TSI matters for a variety of outcomes, the evidence is not quite as strong as it originally seemed to be.
Our second example examines the relationship between minimum wages in Canada and the log of hourly earnings. We focus attention on men between 18 and 24 years of age who immigrated to Canada less than ten years ago. Our sample contains 28,599 observations for the years 2008 to 2019. Except for a few federally-regulated industries, minimum wages in Canada are set at the provincial level. They tend to change infrequently, and they never go down. In fact, although our sample contains observations for 1440 province-month pairs, the minimum wage variable takes on only 63 unique values.
The equation we estimate is
where logearn$_{ipmt}$ is the log of hourly earnings for individual $i$ in province $p$ in month $m$ of year $t$, logmw$_{\kern -0.08333em pmt}$ is the log of the minimum wage, bigcity$_{\kern -0.08333em ipmt}$ is a dummy for being in one of nine large cities, age$_{ipmt}$ is a dummy for being 22 to 24, and the remaining regressors are year fixed effects, calendar month fixed effects, and province fixed effects. The total number of regressors, including the constant term, is 35.
This example is one for which reliable cluster-robust inference is likely to be difficult. We cluster by year and province, but there are only 12 years and 10 provinces. The year clusters are reasonably homogeneous in size; they vary from 2051 to 2723 observations. But the province clusters are very heterogeneous; they vary from 163 (P.E.I.) to 6554 (Ontario). Although there are no empty intersections, the smallest contains just 3 observations, and the largest contains 710.
(ref) contains three panels. Panel C presents some cluster diagnostics, calculated using twowayjack. The coefficients of variation are quite revealing. For partial leverage, there is considerable variation across provinces and intersections, but very little across years. For the $\hat\beta^{(g)}$, there is modest variation when leaving out a province or a year, but very little when leaving out an intersection cluster.
These features of the sample suggest that many methods, perhaps all methods, will not yield reliable inferences. In order to investigate this conjecture, we employ placebo\kern 0.04167em-regression simulations as advocated by MNW-guide. These are similar in spirit to the “placebo laws” simulations of \citet*{BDM_2004}. For each of 100,\kern 0.08333em000 simulations and each province, we generate a sequence of values of a placebo regressor that resembles the actual minimum wage sequences: The value tends to stay constant for a while and then rise by a random amount from time to time in a fashion that is correlated across provinces. This placebo regressor is then added to regression \hyperref[{minwage}]{\tagform@{\ref*{minwage}}}, and we calculate sixteen $P$ values for its coefficient based on all sixteen standard errors used for the actual regression. If regression \hyperref[{minwage}]{\tagform@{\ref*{minwage}}} is correctly specified and any particular way of obtaining $P$ values is valid for our sample, then the fraction of the time that the placebo\kern 0.04167em-regression $P$ value is less than 0.05 should be very close to 0.05, subject to experimental error.
Panels A and B of (ref) show both the actual $P$ values and rejection frequencies for the placebo regressions for all sixteen methods. The conventional $P$ values in Panel A imply that the minimum wage is significant at the 0.01 level for all the one\kern 0.04167em-way clustering methods and at the 0.02 level for all the two\kern 0.04167em-way methods. However, the placebo\kern 0.04167em-regression rejection frequencies vary from 15% to 89%, suggesting that none of the conventional $P$ values should be believed.
In contrast, the jackknife $P$ values in Panel B are greater than 0.05 for one\kern 0.04167em-way clustering by province and for the two\kern 0.04167em-way clustering methods. The placebo\kern 0.04167em-regression rejection frequencies for the one\kern 0.04167em-way methods vary between 9% and 89%, suggesting that they should not be trusted. For the two\kern 0.04167em-term two\kern 0.04167em-way estimator CV$_{\kern -0.08333em3}^{(2)}$, the placebo rejection frequency is just 2.5%. This is in line with existing theory ((ref)) and some of our simulations (e.g.\ (ref)), which both show that CV$_{\kern -0.08333em3}^{(2)}$ can under-reject. For the three\kern 0.04167em-term cluster-jackknife estimators, the placebo rejection frequencies are between 5.7% and 6.5%, which is remarkably good in view of the small numbers of clusters and the cluster diagnostics. For these methods, the $P$ value is 0.081.
Because the number of clusters in each of the primary dimensions is small, and cluster sizes vary greatly in one of them, it seems plausible that even the three-term CV$_{\kern -0.08333em3}$ standard errors are too small. The modest over-rejection rates for the placebo regressions in the second row of the last three columns of Panel B support this conjecture. Thus it could makes sense to use the mixed three-term estimator \hyperref[{31jack}]{\tagform@{\ref*{31jack}}} in this case. The placebo regression rejection rates for CV$_{3,1}^{(3)}$ and CV$_{3,1}^{(\rm max)}$ are 0.0506 and 0.0479, respectively. Like that of CV$_{3}^{(\rm max)}$, these are extraordinarily close to 0.05. The $P$ values for both CV$_{3,1}^{(3)}$ and CV$_{3,1}^{(\rm max)}$ are 0.0954, which is also close to that of CV$_{3}^{(\rm max)}$ at 0.0808.
We conclude that, even though all the conventional one\kern 0.04167em-way and two\kern 0.04167em-way standard errors yield strongly significant results (in the first row of Panel A), all the two\kern 0.04167em-way jackknife standard errors (in the first row of Panel B and in the previous paragraph) yield results that are not even close to significant at the 0.05 level. Thus, we cannot be confident that the minimum wage positively affects hourly earnings based on the evidence from this sample.
It is common to assume that the disturbances in linear regression models are clustered in two dimensions. Unless the regressor(s) of interest are uncorrelated in one or both dimensions, it is therefore necessary to employ a cluster-robust variance estimator that allows for two\kern 0.04167em-way clustering. Unfortunately the most widely-used cluster-robust variance matrix estimator (CRVE), CV$_{\kern -0.08333em1}^{(3)}$, due to CGM_2011, is not guaranteed to be positive definite. Inferences based on it are known to be seriously unreliable in finite samples MNW_2021.
In (ref), we discuss several ways to avoid, or at least ameliorate, the problem of undefined standard errors when a CRVE is not positive definite. Most importantly, we propose a new and simple solution to this problem. For tests of a single restriction, it just involves using whichever of three standard errors is the largest. Two of these are based on one\kern 0.04167em-way clustering in each of the two dimensions, and the third is a three\kern 0.04167em-term two\kern 0.04167em-way standard error. Asymptotically, the latter should always be the largest of the three when there really is two\kern 0.04167em-way clustering, but it may not be the largest (and may indeed not be defined) in finite samples. In many cases, our so\kern 0.04167em-called max-se procedure yields results identical to those from the corresponding three\kern 0.04167em-term two\kern 0.04167em-way CRVE, but it can yield substantially lower (and more accurate) rejection frequencies in some cases.
The second, and in our view more widely applicable, contribution of the paper is to propose and study two\kern 0.04167em-way cluster-jackknife CRVEs. Recent work on the cluster-jackknife, or CV$_{\kern -0.08333em3}$, CRVE for one\kern 0.04167em-way clustering Hansen-jack,Hansen_2025,MNW-bootknife,MNW-influence suggests that it can perform much better in finite samples than the usual CV$_{\kern -0.08333em1}$ CRVE. It therefore seems attractive to extend it to the two\kern 0.04167em-way case. This is remarkably simple. We just need to perform three sets of cluster-jackknife calculations, one for each of the two dimensions, and then a third one for their intersections. In many cases, this is straightforward, although cluster fixed effects do raise some computational issues ((ref)), and the calculations can be costly when the number of intersections is large, especially when there are cluster fixed effects. We provide a Stata package called twowayjack that implements our methods and also calculates some cluster diagnostics MNW-influence; see (ref).
In (ref), we study rejection frequencies for $t$-tests based on eight different cluster-robust standard errors. Four of them are of the usual CV$_{\kern -0.08333em1}$ type, and the other four are of the CV$_{\kern -0.08333em3}$ type. In most cases, tests based on the CV$_{\kern -0.08333em3}$ max-se standard error yield the most reliable inferences. Even when they do not, they only perform slightly worse than whatever procedure(s) perform better, and they are usually much more reliable than all of the CV$_{\kern -0.08333em1}$-based tests.
Because most of our simulations involve two\kern 0.04167em-way cluster fixed effects, three\kern 0.04167em-term variance matrices based on either CV$_{\kern -0.08333em1}$ or CV$_{\kern -0.08333em3}$ tend not to be positive definite, so the versions that use an eigen-decomposition ((ref)) can differ greatly from the versions that do not. This always reduces rejection frequencies, which is a good thing for CV$_{\kern -0.08333em1}$ tests but usually a bad thing for CV$_{\kern -0.08333em3}$ tests. Tests based on two\kern 0.04167em-term variance matrices usually reject even less frequently than tests based on three\kern 0.04167em-term variance matrices with the eigen-decomposition. Thus we do not recommend using tests based on either CV$_{\kern -0.08333em3}^{(2)}$ or CV$_{\kern -0.08333em3}^{(3+)}$.
Our simulations show that precisely how the data are generated can have a large effect on finite\kern 0.04167em-sample performance. All the tests are most likely to perform poorly when the number of clusters in either dimension is small, cluster sizes vary greatly, there are many empty intersections, the number of regressors that are clustered in one or both dimensions is large, or either the disturbances or the regressor(s) of interest are only weakly correlated in both dimensions. In many of these cases, alternative test statistics tend to perform quite differently. For example, weak intra-cluster correlation is particularly problematic for the CV${\kern -0.08333em1}$-based, two-term, and 3+ variants. The CV$_{\kern -0.08333em3}^{(\max)}$ tests generally stay close to their nominal size except in the most extreme weak-correlation designs or when combined with other adverse features.
In practice, it can often be illuminating to employ placebo regression simulations, as in (ref). These will show how well alternative tests perform for the particular model and dataset under study. It is probably safe to rely on CV$_{\kern -0.08333em3}^{\kern 0.08333em(\max)}$-based tests if they perform well in these simulations, or perhaps on some other tests if they perform better.