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.
86,913 characters · 20 sections · 52 citation commands
Estimation and exclusion restrictions in clustered linear models
Clustered datasets---such as panel, network, spatial, or group-structured data in which individuals or firms are observed repeatedly or nested within groups---have become increasingly common in empirical research. These settings present several methodological challenges. One is the need to accommodate heterogeneity, often addressed by including fixed effects and flexible time or group trends, which typically requires many control variables. A more fundamental difficulty is within-cluster dependence: observations in the same cluster may be correlated due to, e.g., spatial and network interference, spillover effects, and time-series dependence. Such dependencies complicate statistical inference and limit the effective aggregation of information across observations. \looseness=-1
The concept of exogeneity becomes more nuanced in the presence of clustered data. As we demonstrate, assuming only a per-observation exclusion restriction---uncorrelatedness of the error with the regressor within an observation---may deliver no identifying variation for consistent estimation of structural parameters. By contrast, strict exogeneity---each error term uncorrelated with all regressors in the cluster---is often implausible in many empirical contexts. \looseness=-1
To bridge strict exogeneity and unrestricted dependence, we assume the error term is uncorrelated with a subset of within-cluster regressors; the subset is application-specific. In panel data, it is common to assume errors are mean-zero conditional on current and past (but not future) regressors, allowing feedback from shocks to future policies. In spatial contexts, localized spillovers motivate analogous restrictions: regressors corresponding to sufficiently distant units within a cluster may plausibly be uncorrelated with a given error term, even if nearby regressors are not. Examples include spatial leakage in deforestation jayachandran2017cash, spillovers in policing blattman2021place, and fiscal intervention experiments in rural Kenya egger2022general. In network settings, treatment effects may propagate along social ties, so regressors associated with nodes that are unconnected—or sufficiently distant in the network—may be uncorrelated with a given error term paluck2016changing. These partial exogeneity assumptions more accurately reflect the dependence structures present in many empirical applications.\looseness=-1
When only a subset of exclusion restrictions holds, it becomes inappropriate to treat regressors as fixed or to condition on them. Their stochastic variation may correlate with the errors, undermining standard arguments that rely on conditioning on regressors. As a result, even estimators like OLS must be interpreted as ratios of stochastic quadratic forms, which introduces several challenges. \looseness=-1
The most well-known is Nickell bias nickell1981biases, which occurs when the expected value of the OLS error numerator is nonzero, leading to asymptotic bias. While prominent in dynamic panels with fixed effects and lagged outcomes, analogous bias arises more broadly under clustered dependence. A second challenge is inference: standard variance estimators often fail to account for dependence among cross-cluster terms in the quadratic form. Finally, the denominator may remain asymptotically random, unlike in classical settings, leading to weak identification and further complicating inference. \looseness=-1
This paper studies the estimation of structural parameters in linear regression models with clustered data, high-dimensional controls, and researcher-specified exclusion restrictions. A key feature of the framework is the grouping of observations into disjoint clusters, permitting within-cluster dependence and independence across clusters. The setting accommodates unbalanced panels with multiple high-dimensional fixed effects, as well as spatial, network, and group-structured data. Our approach extends dynamic panel methods to a broader class of models, addressing the central challenges of bias, inference, and identification. \looseness=-1
Our first contribution is to characterize a class of correctly centered internal instrument estimators that remove the leading asymptotic bias of OLS. The approach accommodates high-dimensional exogenous controls and adapts to the exclusion structure specified in the application. The resulting estimator is an internal instrument IV estimator, defined as the one closest to OLS in a specified norm. The estimator is easy to implement and asymptotically efficient under some conditions. \looseness=-1
The procedure has a simple interpretation: for each observation, controls are partialled out using only those observations whose errors are uncorrelated with the observation's regressor, yielding an observation-specific leave-out projection. A just-identified IV regression is then performed on the transformed equation, using the original regressor as the instrument. \looseness=-1
We further characterize the efficiency loss from varying the strength of exclusion restrictions. In particular, assuming exogeneity only at the observation level eliminates all identifying variation when fixed effects are present. This observation underscores the importance of carefully specifying the exogeneity structure in applied work. \looseness=-1
Our second contribution addresses inference. We show that, outside special cases such as (i) strict exogeneity and (ii) block-diagonal residualization (e.g., models with only cluster fixed effects), the estimator's numerator is a nontrivial quadratic form in the errors. In these more general settings, standard cluster-robust variance estimators may fail, and valid inference requires central limit theorems for clustered quadratic forms. These issues are especially acute in models with many high-dimensional controls, such as two-way fixed effects verdier2018estimation. We propose a new central limit theorem and variance estimator that account for this structure. \looseness=-1
Our third contribution addresses weak identification, which arises when the internal instrument---though valid---captures little identifying variation due to many controls or weak exclusion restrictions. As in the classical dynamic panel literature blundell1998initial, bun2010weak, weak identification can compromise standard inference. We develop inference procedures that remain valid even when the estimator's denominator exhibits substantial sampling variability. \looseness=-1
We propose a comprehensive empirical strategy that combines a correctly centered estimator with identification-robust inference and confidence sets. We apply it to a prominent randomized evaluation of a large-scale fiscal intervention in rural Kenya egger2022general. A key challenge in this setting is spatial interference: treatment in one village affects outcomes in neighboring villages, complicating the choice of plausible exclusion restrictions. We show that weakening the exclusion restrictions yields less precise estimates and wider confidence sets. \looseness=-1
This paper contributes to several strands of the econometrics literature. It relates to the extensive work on linear dynamic panel data models---a mature field with too many important contributions to list exhaustively. For recent overviews on correcting Nickell bias, see the surveys by okui20211 and bun2015dynamic. Our estimation approach is more closely aligned with the internal instrument framework developed in anderson1981estimation, arellano1991some, ahn1995efficient, alvarez2003time, and blundell1998initial, and we extend the framework to accommodate more general models. In addition, we directly address weak identification concerns noted in the dynamic panel literature, particularly in bun2010weak. We also contribute to the emerging literature on quadratic central limit theorems for clustered data structures. These results have recently been employed in many-instrument panel settings, such as ligtenberg2023inference, and are crucial for understanding the distribution of estimators that are quadratic functionals of the errors. \looseness=-1
The remainder of the paper is organized as follows. Section (ref) describes the data structure, a plausible set of exclusion restrictions, and two modeling perspectives---outcome-based and design-based. Section (ref) characterizes correctly centered internal instrument estimators, introduces our proposed estimator as the solution to an efficiency optimization problem, and interprets it as a leave-out internal instrument procedure. Section (ref) addresses uncertainty quantification, and Section (ref) establishes a central limit theorem for quadratic forms with clustered data. Section (ref) develops inference procedures robust to weak identification and a corresponding variance estimator. Section (ref) illustrates the full strategy in an empirical application, showing how both the estimator and, especially, its uncertainty depend on the maintained set of assumptions. \looseness=-1
\paragraph{Notation.} We use \(0< c < 1 \) and \( C > 0 \) to denote generic finite constants that may vary across equations but do not depend on the sample size. For a matrix \( A \), we let \( A^+ \) denote its Moore–Penrose generalized inverse, \( \operatorname{tr}(A) \) its trace, \( \|A\|_F^2 = \operatorname{tr}(A'A) \) its Frobenius norm squared, and \( \|A\| \) its operator (spectral) norm. \looseness=-1
We consider the estimation of a structural parameter $\beta$ in the linear regression model
where \(\ell\) indexes the \(n\) observations in the dataset. The parameter of interest, \(\beta\), is a scalar, though most of our results extend to vectors of fixed dimension. The vector of strictly exogenous controls \(w_\ell\) has dimension \(K\), which may be large but is assumed to satisfy \(K < n\). We denote by \(W = [w_1, \dots, w_n]'\) the matrix of control variables. For simplicity, we treat \(W\) as non-random; equivalently, all results may be interpreted conditional on \(W\). Throughout, we assume that \(W\) has full rank and define the projection matrix $M = I_n - W(W'W)^{-1}W',$ which projects out the variation associated with the controls \(W\).\looseness=-1
In contrast to standard cross-sectional settings with independent observations, we allow for dependence within clusters. Specifically, we assume the data can be partitioned into \(N\) disjoint clusters, with independence across clusters and arbitrary dependence within them. Formally, let \(\{ S_i \}_{i=1}^N\) be a partition of the index set \(\{1, \dots, n\}\), so that \(\cup_{i=1}^N S_i = \{1, \dots, n\}\) and \(S_i \cap S_j = \emptyset\) for \(i \ne j\). Let \(i(\ell)\) denote the index of the cluster containing observation \(\ell\), and define \(T_i = |S_i|\) as the size of cluster \(i\), so that \(n = \sum_{i=1}^N T_i\). This cluster structure reflects the sampling design and is essential for the application of laws of large numbers and central limit theorems. \looseness=-1
Assume that random regressors $x_\ell$ and errors $e_\ell$ have uniformly bounded second moments. There are a variety of exclusion assumptions one may make in regression (ref). The weakest of these is the contemporaneous exclusion restriction, \( \mathbb{E}[x_\ell e_\ell] = 0 = \mathbb{E}[e_\ell] \) for all \( \ell \), which defines the parameter \( \beta \). When all observations are independent, this restriction (under mild regularity conditions) ensures that ordinary least squares (OLS) yields a consistent estimator of \( \beta \). However, as we show later, in the presence of clustered data, when the regression includes cluster fixed effects, the contemporaneous exclusion assumption alone fails to yield a consistent estimator of \( \beta \).
At the other extreme is the strict exogeneity condition: \( \mathbb{E}[e_\ell \!\mid\! x] = 0 \) where \( x = [x_1, \dots, x_n]' \), which implies \( \mathbb{E}[x_{\tilde\ell}e_\ell] = 0 \) for all \( \tilde\ell \) and \( \ell \). Under strict exogeneity, OLS is unbiased conditional on $x$. In particular, by the Frisch-Waugh-Lovell theorem, we have:
where \( y = [y_1, \dots, y_n]' \), \( e = [e_1, \dots, e_n]' \) and \( \tilde x = Mx \) denotes the residual from regressing \( x \) on \( W \). Conditionally on \( x \), the regression errors \( \{e_\ell\} \) are the only source of randomness. Since \( \hat\beta^\mathrm{LS} \) is linear in the errors, standard results---such as unbiasedness, the law of large numbers, and the central limit theorem---hold. The only complication introduced by the clustered data structure is the need to adjust standard errors appropriately: cluster-robust standard errors are required because the relevant level of independence is at the cluster level. \looseness=-1
Strict exogeneity is often too strong to hold in applications with clustered data, where unmodeled dependence between regressors and outcomes within clusters is likely. Researchers should tailor the exclusion restrictions to reflect institutional or structural features of the setting. We assume that the researcher can specify a matrix of exclusion restrictions, denoted by $\mathcal{E}$, which is an $n\times n$ indicator matrix encoding the moment conditions the researcher is willing to impose. Specifically, \( \mathcal{E}_{\tilde\ell\ell} = 1 \) encodes the assumption \( \mathbb{E}[x_{\tilde\ell} e_\ell ] = 0 \), while \( \mathcal{E}_{\tilde\ell\ell} = 0 \) encodes the absence of a restriction and allows \( \mathbb{E}[x_{\tilde\ell} e_\ell ] \ne 0 \).\looseness=-1
Due to the assumed independence across clusters, we automatically have \( \mathbb{E}[x_{\tilde\ell}e_\ell] = 0 \) whenever \( i(\ell) \ne i(\tilde\ell) \). For notational simplicity, we assume \( \mathcal{E}_{\tilde\ell\ell} = 1 \) in such cases. Thus, zero entries in \( \mathcal{E} \) may only occur within the blocks corresponding to clusters. Zeros within these blocks represent the absence of exclusion restrictions. We consider the problem of estimating \( \beta \) in regression (ref) under the exclusion restrictions encoded in the indicator matrix \( \mathcal{E} \). \looseness=-1
Assumption (ref) ensures contemporaneous exogeneity, \( \mathbb{E}[ x_\ell e_\ell] = 0 \), which underlies the desirable properties of OLS in standard cross-sectional settings.
When strict exogeneity fails, the OLS estimator may suffer from asymptotic bias and, in some cases, may even be inconsistent. \looseness=-1
This lemma is a relatively straightforward extension of the well-known Nickell bias nickell1981biases. Nickell’s original work focused on a more restricted setting commonly referred to as a “dynamic panel data” model. \looseness=-1
The source of Nickell’s bias is well understood: differencing out fixed effects mixes observations from different periods, and future values of the regressor may correlate with current errors. Our framework is more general than the standard dynamic panel specification in (ref). It is worth emphasizing that the result established by nickell1981biases---and generalized in Lemma (ref)---is best understood as an asymptotic bias: a term which, when subtracted from the OLS estimator, yields a consistent (or asymptotically unbiased) estimator. Asymptotic bias is distinct from finite sample bias, as it ignores the randomness in the denominator of the OLS expression. While this randomness becomes negligible asymptotically, it is present in finite samples.\looseness=-1
Model (ref) combined with Assumption (ref) can be viewed as an instance of outcome modeling, since it specifies the outcome equation. In cross-sectional settings, however, OLS estimators perform well when either the outcome or the treatment equation is correctly specified. The alternative approach—based on the treatment equation—is referred to as design-based and is commonly used when researchers understand the random assignment of treatment, such as in randomized controlled trials (RCTs). In this paper, we define a design-based model by specifying the following treatment equation, along with Assumption (ref):
The asymptotic bias expression from Lemma (ref) can equivalently be restated in the design-based framework:
Figure (ref) illustrates this asymptotic bias in a simulation designed to mimic the network-interference setting in Example (ref). We consider a DGP based on a stochastic block model where clusters represent school classes and within-cluster edges represent friendships. The researcher observes the friendship network but not the intensity of each connection, denoted as weight $G_{\tilde\ell\ell}=G_{\ell\tilde\ell}\sim\mathrm{Exp}(1)$ for connected pairs and zero otherwise. Treatment is binary and assigned independently across units with probability $\mu_{i(\ell)}$, which varies across clusters. The outcome $y_\ell$ includes a spillover component depending on the treatment of friends weighted by unobserved friendship intensities: $ y_\ell = \beta \, x_\ell + \alpha \sum_{\tilde\ell} G_{\ell \tilde\ell} \, x_{\tilde\ell} + e_\ell. $ The least-squares estimator controlling for cluster fixed effects exhibits substantial bias in this setting. Full details of the simulation design are provided in Appendix (ref). Despite random assignment of treatment, the OLS estimator can be highly misleading when interference is present. The bias arises from the correlation induced by spillovers within clusters, illustrating the Nickell-type bias in clustered regressions with network interference.
The simulation illustrates the usefulness of the design-based perspective in settings where the treatment assignment mechanism is credible but the outcome equation may be difficult to specify. The primary exposition in this paper is framed within the outcome modeling approach. However, all results can equivalently be reformulated under the design-based perspective with only minor modifications. We provide such restatements in the remarks following the main theorems. Additionally, where appropriate, we discuss the possibility of a doubly robust formulation, in which the researcher assumes that either the outcome equation or the treatment equation is correctly specified, without committing to which one.
\looseness=-1
This subsection shows that unbiased estimation is generally impossible once regressors are random, even under standard orthogonality conditions, and motivates a weaker requirement that we call correct centering. Assume we have clustered data \( \{y_\ell, x_\ell\} \) generated according to model (ref) under Assumption (ref). The set of controls \( W \), the cluster assignments \( S_i \), and the exclusion restriction matrix \( \mathcal{E} \) are given to the econometrician and treated as fixed. We consider the problem of estimating the parameter \( \beta \). Let \( \mathcal{F} \) denote the class of all joint distributions \( F \) over \( (x, y) \) that satisfy model (ref) under Assumption (ref).
In fact, the proof of Lemma (ref) establishes something even stronger: in a cross-sectional setting with an independent sample $\{y_i, x_i\}_{i=1}^n$ from a correctly specified regression model $y_i = \beta x_i + e_i$, where the regressor $x_i$ is random and $\mathbb{E}[e_i] = \mathbb{E}[x_i e_i] = 0$, no unbiased estimator of $\beta$ exists. The OLS estimator is not unbiased in this setting because the randomness of the regressors induces randomness in the OLS denominator, making it impossible to pass the expectation operator through the ratio. This observation motivates the need to consider a broader and more appropriate class of estimators, in which the normalization is separated from the centering.
\looseness=-1 Correct centering is a weaker condition than unbiasedness. The key distinction lies in the treatment of randomness in the denominator. When the denominator is non-random—so that $\mathbb{E}[C_2(x)] = C_2(x)$—the two concepts coincide. In a correctly specified cross-sectional regression with random regressors, OLS is correctly centered but not unbiased.
For many of the estimators we consider, we can establish a form of the law of large numbers (under broadly applicable technical conditions) that ensures $C_2(x) - \mathbb{E}[C_2(x)] \xrightarrow{p} 0$ and $C_1(x, y) - \mathbb{E}[C_1(x, y)] \xrightarrow{p} 0$—with any necessary normalization absorbed into the functions themselves. For instance, Lemma (ref) establishes that $C_2(x) = \frac{1}{n}x'Mx = \frac{1}{n}\mathbb{E}[x'Mx] + o_p(1)$ and $C_1(x, y) = \frac{1}{n}\mathbb{E}[x'My] + o_p(1)$. Regarding the denominator, what matters is the size of its uncertainty relative to its signal—that is, the strength of identification. If $\frac{C_2(x)}{\mathbb{E}[C_2(x)]} \xrightarrow{p} 1$ and the above law of large numbers holds, then correct centering implies consistency (or asymptotic unbiasedness) of the estimator.
To better understand the concept, it is helpful to recall that in a standard cross-sectional overidentified IV regression, the two-stage least squares (2SLS) estimator is not correctly centered. When the number of instruments is large, the resulting bias can be of first-order importance and is known as the many-instrument problem. In contrast, a one-step just-identified IV estimator—as well as jackknife or leave-one-out IV estimators in overidentified cross-sectional settings—remains correctly centered whenever the instruments are exogenous.
Here, (POP) stands for the partialling-out property, and (CC) refers to correct centering. Under condition (POP), the estimator satisfies: \[ \hat\beta^A = \frac{x' A y}{x' A x} = \frac{x' A (\beta x + W \delta + e)}{x' A x} = \beta + \frac{x' A e}{x' A x}. \] Condition (CC) guarantees that \( \mathbb{E}[x' A e] = 0 \), ensuring correct centering of the estimator. \looseness=-1
The class of estimators described in Lemma (ref) coincides with the class of just-identified linear IV estimators that use linear internal instruments. Specifically, define the instrument \( z = A'x \), a linear transformation of the regressors. Condition (CC) implies the exogeneity condition \( \mathbb{E}[z_\ell e_\ell] = 0 \) for all \( \ell \). If we estimate model (ref) via a just-identified IV regression using \( z \) as the instrument for \( x \), and including \( W \) as controls, then a straightforward generalization of the Frisch-Waugh-Lovell theorem for IV estimators yields: \[ \hat\beta^{\mathrm{IV}} = \frac{z' M y}{z' M x} = \frac{x' A M y}{x' A M x} = \frac{x' A y}{x' A x} = \hat\beta^A, \] where the final equality follows from the partialling-out property (POP). \looseness=-1
Our approach is closely related to methods based on internal instruments for dynamic panel data, such as anderson1981estimation, arellano1991some, and ahn1995efficient. Both arellano1991some and ahn1995efficient stack all available exclusion restrictions and estimate the model using overidentified 2SLS or GMM. In contrast, we construct a just-identified linear combination of internal instruments and implement a one-step IV procedure. This delivers a correctly centered estimator. Overidentified procedures generally do not enjoy this property, as the combination of many moment conditions together with an estimated weighting matrix introduces additional randomness. When the degree of overidentification is large, this can result in substantial many-instrument bias alvarez2003time. Different choices of the matrix $A$ satisfying restrictions (CC) and (POP) span a broad class of linear internal instruments and, in particular, include the infeasible optimally weighted instrument of arellano1991some as a special case. \looseness=-1
Let $\mathcal{A}$ denote the set of $n \times n$ matrices satisfying conditions (POP) and (CC). This set forms a linear subspace within the space of $n \times n$ matrices. The full space of $n \times n$ matrices has dimension $n^2$. Condition (POP), which is equivalent to $AW = 0_{n \times K}$, imposes $nK$ linear constraints. Condition (CC) imposes an additional $L$ zero restrictions, where $L$ is the number of entries $A_{\tilde\ell \ell}$ required to be zero due to (CC). Therefore, the dimension of $\mathcal{A}$ is at least $n(n - K) - L$. If the dimension of $\mathcal{A}$ is positive, then our estimator class contains infinitely many candidates, raising the natural question: Which matrix $A$ should we use? As a default, we propose the choice
The equality holds due to condition (POP). The Frobenius norm defines a Euclidean geometry on the space of matrices, with the corresponding inner product given by $\langle A, B \rangle_F = \operatorname{tr}(A'B)$. The optimization problem (ref) admits a unique solution: the orthogonal projection of $M$ onto the subspace $\mathcal{A}$ with respect to the Frobenius inner product. In the lemma below, we show that this choice of $A^*$ can be motivated by asymptotic efficiency considerations.
The function $V(A)$ represents the asymptotic variance of the IV estimator under strong identification, i.e., when ${x'Ax}/{\mathbb{E}[x'Ax]} \xrightarrow{p} 1$. This efficiency criterion follows the classical approach of chamberlain1987asymptotic. Lemma (ref) shows that, under ideal conditions, our proposed estimator using $A^*$ from (ref) achieves asymptotic efficiency within the class of correctly centered estimators. The two ideal conditions are: first, conditional homoskedasticity, $\mathbb{E}[ee' \!\mid\!x] = \sigma^2 I_n$, which implies that the optimal instrument $z^* $ minimizes the first stage prediction error over all $z = A'x$ with $A \in \mathcal{A}$:
Second, the assumption $\mathbb{E}[xx'] = \lambda I_n$ imposes stochastic homogeneity of the regressors and ensures the equivalence between the optimization problems (ref) and (ref).\looseness=-1
It is natural that the efficiency of an estimator depends on the distributional properties of the regressor $x$, especially in settings like ours where $x$ must be treated as random. Any prior information---such as the scale of $x$ or its dynamic dependence---can, in principle, be used to improve efficiency, for example, through cluster-specific weighting or time transformations (as in GLS). However, caution is needed: the matrix $A$ must remain independent of the randomness in either $x$ or $e$ to avoid bias, such as the many-instrument bias observed in the two-step GMM estimator of arellano1991some. If additional structure is available, such as $\mathbb{E}[xx'] = \Omega$, then asymptotic efficiency is attained by solving: $$ A_{\Omega}^* = \arg\min_{A\in \mathcal{A}}\operatorname{tr}\left((A-I_n)'\Omega(A-I_n)\right) = \arg\min_{A\in \mathcal{A}} \|\Omega^{1/2}(A-I_n)\|_F. $$ Most results in this paper extend naturally to this weighted setting, except for the leave-one-out projection characterization discussed below in Theorem 1(ref).\looseness=-1
This subsection derives and discusses the solution to the optimization problem (ref).\looseness=-1
We begin by introducing what we refer to as leave-out projections. For each index $\tilde{\ell}$, define $\mathcal{E}_{\tilde{\ell},\cdot}$ to be the $\tilde{\ell}$th row of the matrix $\mathcal{E}$, which contains indicators for those indices $\ell$ such that $\mathbb{E}[x_{\tilde{\ell}} e_\ell] = 0$. In other words, $\mathcal{E}_{\tilde{\ell},\cdot}$ identifies all observations for which $x_{\tilde{\ell}}$ is uncorrelated with their error term. Let $\mathbf{i}_{\tilde{\ell}} = \operatorname{diag}(\mathcal{E}_{\tilde{\ell},\cdot})$ denote the $n \times n$ diagonal selection matrix with ones on the diagonal corresponding to entries where $\mathcal{E}_{\tilde{\ell} \ell} = 1$, and zeros elsewhere. Then define $W_{\tilde{\ell}} = \mathbf{i}_{\tilde{\ell}} W,$ an $n \times K$ matrix that retains the rows $w_\ell$ for which $\mathcal{E}_{\tilde{\ell} \ell} = 1$, and sets the remaining rows to zero. The corresponding leave-out projection matrix is given by $M^{(\tilde\ell)}=I-W_{\tilde \ell}(W_{\tilde \ell}'W_{\tilde \ell})^{+}W_{\tilde \ell}',$ which partials out controls using only observations whose errors are uncorrelated with $x_{\tilde\ell}$. Unlike the full projection $M$, it excludes data points that violate exogeneity.\looseness=-1
Characterization (ref) in Theorem (ref) highlights the analytical tractability of $A^*$ by noting its connection with linear projection.\looseness=-1
Characterization (ref) shows that $A^*$ closely resembles the OLS projection matrix $M$, with the difference given by $A^* - M = -BM$, where $B$ is a sparse matrix. Specifically, $B_{\tilde\ell\ell}$ is nonzero only if $\mathbb{E}[x_{\tilde\ell} e_\ell] \ne 0$. As a result, $B$ is block-diagonal (with blocks corresponding to clusters) and has zeros on the main diagonal. Thus, $A^*$ can be computed separately within each cluster, which is computationally appealing. \looseness=-1
Finally, part (ref) of Theorem (ref) offers an intuitive interpretation of the solution. For each observation $\tilde\ell$, we perform an observation-specific transformation of the original model (ref) by partialling out the control $w_{\tilde\ell}$ using only observations uncorrelated with $x_{\tilde\ell}$. Define the residualized outcome as $ y_{\tilde\ell}^*=y_{\tilde\ell}-w_{\tilde\ell}^\prime\hat\delta^y_{{\tilde\ell}}$, where $\hat\delta^y_{{\tilde\ell}}$ is obtained from regressing $y$ on $w$ using only observations with errors uncorrelated with $x_{\tilde\ell}$. Similarly define $x_{\tilde\ell}^*$ and $e_{\tilde\ell}^*$. These definitions yield the transformed equation: \looseness=-1
Since $e^*_{\tilde\ell}$ is constructed using only observations uncorrelated with $x_{\tilde\ell}$, we have $\mathbb{E}[x_{\tilde\ell}e_{\tilde\ell}^*]=0$. Hence, $x_{\tilde\ell}$ is a valid instrument in (ref), and the IV estimator is $\hat\beta=\frac{\sum_{\tilde\ell}x_{\tilde\ell}y^*_{\tilde\ell}}{\sum_{\tilde\ell}x_{\tilde\ell}x^*_{\tilde\ell}}=\frac{x'A^*y}{x'A^*x}=\hat\beta^{A^*}$.\looseness=-1
\paragraph{Price of robustness.} Our estimator can be seen as a robust alternative to OLS when strict exogeneity, $\mathbb{E}[e_\ell \!\mid\!x] = 0$, is not credible. The exclusion matrix $\mathcal{E}$ determines the types of violations it is robust to: Example (ref) allows $\mathbb{E}[x_{i,t+1}e_{it}] \ne 0$, while Example (ref) allows broader future dependence. A natural question is whether this robustness reduces efficiency when strict exogeneity holds. Under the assumptions of Lemma (ref), we have $ V(A^*)=\frac{\sigma^2}{\lambda\operatorname{tr}(A^*)}=\frac{\sigma^2}{\lambda\|A^*\|_F^2}, $ where $\operatorname{tr}(A^*)$ reflects the effective sample size. In the standard panel model with strict exogeneity and homoskedasticity, OLS achieves the efficiency bound with $\operatorname{tr}(M)=\sum_{i=1}^N(T_i-1)=n-N$, as one degree of freedom is lost per cluster due to demeaning. Each block $M_{ii}$ contributes $\operatorname{tr}(M_{ii}) = T_i - 1$.\looseness=-1
In Example (ref), we have $\operatorname{tr}(A_{ii}) = T_i - 1 - \frac{1}{T_i}$; the term $-\frac{1}{T_i}$ reflects the loss of information due to robustness against one-period feedback. In Example (ref), $\operatorname{tr}(A_{ii})=T_i-1-\sum_{t=2}^{T_i}\frac{1}{t}$. With $\sum_{t=2}^{T_i}\frac{1}{t} \ge \log( T_i+1) - \log(2)$, this shows a substantially greater efficiency loss when guarding against dependence on any future regressor. Finally, in a panel model with fixed effects, assuming contemporaneous exogeneity only $\mathbb{E}[x_\ell e_\ell] = 0$ leads to $A^* = 0$. In this case, the model contains no usable identifying variation. \looseness=-1
Let $x_i$ denote the $T_i \times 1$ vector of regressor values within cluster $S_i$---that is, all $x_\ell$ for which $i(\ell) = i$---and define $e_i$ analogously. Due to the cluster structure, the pairs $(x_i, e_i)$ are independent across $i$. Let $A_{ij}$ denote the $T_i \times T_j$ blocks of the matrix $A\in\mathcal{A}$, corresponding to observations in clusters $S_i$ and $S_j$, respectively. The numerators of the estimation errors take the form of quadratic expressions: $$ \quad \hat\beta^A-\beta=\frac{x'Ae}{x'Ax} \quad \text{and} \quad x'Ae=\sum_{i,j}x_i'A_{ij}e_j. $$ As these are quadratic forms in random $(x_i,e_i)$, the variance of the estimators has a non-standard structure. In particular, establishing asymptotic normality may require a CLT for quadratic forms. \looseness=-1
\paragraph{Special case: block-diagonal structure.} It is worth highlighting a special case ---where the numerator reduces to linear forms. If the off-diagonal blocks $A_{ij}$ are equal to zero for all $i\neq j$, then $x'Ae=\sum_i x_i'A_{ii}e_i$. Since $(x_i, e_i)$ are independent across clusters, standard inference follows directly from the law of large numbers and central limit theorem, provided the number of clusters $N \to \infty$ and that the clusters satisfy a negligibility condition. This scenario arises in Examples (ref) and (ref) when only cluster fixed effects $\alpha_i$ are included in $W$. In that case, usual clustered standard errors can be used. \looseness=-1
However, this structure breaks down when more complex controls are included. \looseness=-1
Example (ref) highlights several key points. First, the difference between $e_i^*$ and $e_i$ reflects estimation error from recovering the teacher effects $\gamma_g$. Usage of the full sample to estimate these effects mixes observations and induces complex dependence, complicating inference. If $\gamma_g$ estimates are very precise, then $\sum_i x_{i1} e_i^*$ will be close to $\sum_i x_{i1} e_i$, and ignoring the dependence may still yield valid inference.\looseness=-1
Second, the OLS transformation for a student in group $(1,3)$ involves not only students with shared teachers ($g=1$ or $g=3$), but also those indirectly connected, e.g., students in $(2,4)$, who share a teacher with others in $(1,4)$ or $(2,3)$. More generally, we can view students as nodes in a network, where edges represent shared teachers. OLS demeaning propagates information across this network, creating strong dependencies within connected components and complicating standard inference methods.\footnote{While OLS demeaning is efficient under the benchmark conditions of Lemma (ref), one may prefer a less efficient estimator that induces less dependence. In Example (ref), instead of using the full sample to estimate $\gamma_3 - \gamma_1$, one could rely only on students with the same teacher history, e.g., $\overline{\Delta y}_{(1,3)}.$ This results in independent “super-clusters” defined by the teacher pairs: $(1,3)$, $(1,4)$, $(2,3)$, and $(2,4)$. One can go further by forming pairs of students with identical teacher histories and differencing them, creating many small, independent clusters, then using clustered standard errors on that level. However, this comes at a cost: the resulting estimator may be highly inefficient due to poor estimation of fixed effects. This strategy is explored in verdier2018estimation. \looseness=-1}
Denote $\lambda_i = \mathbb{E}[x_i]$ and $v_i = x_i - \lambda_i$. Then, $$ x'Ae=\sum_{j=1}^N\left(\sum_{i=1}^N\lambda_i^\prime A_{ij}+v_j'A_{jj}\right)e_j+\sum_{j=1}^N\left(\sum_{i\neq j}v_i'A_{ij}\right)e_j=\sum_{j=1}^N\omega_j e_j+\sum_{j=1}^NQ_je_j. $$ where $Q_j = \sum_{i \ne j} v_i' A_{ij}$ is mean-zero and independent of $(e_j, v_j)$. The above is the so-called Hoeffding decomposition hoeffding1948class; the two sums are the first two U-statistics, which are uncorrelated with each other by construction. Thus, the variance decomposes as: \looseness=-1 $$ \mathbb{V}(x'Ae)=\mathbb{V}\left(\sum_{j=1}^n\omega_je_j\right)+\mathbb{V}\left(\sum_{j=1}^nQ_je_j\right). $$ The first term is the first order U-statistic with $\mathbb{V}\left(\sum_{j=1}^n\omega_je_j\right)=\sum_{j=1}^n\mathbb{V}\left(\omega_je_j\right)$. The second term has correlated summands:
Standard cluster-robust variance formulas omits the final term $\sum_{j=1}^N\sum_{i\neq j}\mathbb{E}[v_i'A_{ij}e_jv_j'A_{ji}e_i]$, which captures cross-cluster dependence. This term vanishes in two important cases: (i) under strict exogeneity, $\mathbb{E}[e_j v_j'] = 0$ for all $j$; or (ii) when $A$ is block-diagonal, i.e., $A_{ij} = 0$ for $i \ne j$. These are the settings where standard clustered standard errors yield valid inference. \looseness=-1
We establish a result on the asymptotic Gaussianity of the numerator of the estimation error, $x'Ae$. This result enables the construction of a hypothesis test for $H_0: \beta = \beta_0$, as well as the development of a confidence set for $\beta$ in the presence of clustered data. Implementation details are provided in Section (ref). Subsection (ref) presents a non-technical exposition and discusses the assumptions required to establish Gaussianity. Subsection (ref) is intended for theoretical econometricians and contains the technical details. Readers whose focus is more applied may skip it.
Two main challenges arise when establishing Gaussianity. First, the statistic includes not only a linear term but also a quadratic form in random shocks. Second, intra-cluster dependence introduces complexities, raising the question of how to restrict such dependence and what trade-offs arise between dependence and cluster size. Our asymptotic framework requires the number of clusters to grow, i.e., $N \to \infty$.
The first challenge is addressed by developing a new Central Limit Theorem (Lemma (ref)) for quadratic forms of independent, normalized random vectors. These vectors represent observations within each cluster. The theorem establishes that, under certain negligibility conditions, a quadratic form involving matrix-weighted sums of such vectors converges in distribution to a normal distribution. The core requirement is that the contribution of any single cluster—as measured, for instance, by the Frobenius norm of the corresponding submatrix—becomes negligible relative to the total variance. A key strength of this result is that it avoids direct control over cluster sizes; the cluster size enters only implicitly through restrictions on the weight matrices.
Our main result---Gaussianity of $x'Ae$---is formalized in Theorem (ref) for a general matrix $A$ and later specialized to our proposed estimator matrix $A^*$. Theorem (ref) combines a Lyapunov-type CLT for linear forms with the new quadratic CLT. Two sets of assumptions are required.
The first set (Assumption (ref)) imposes conditions on the distribution of errors, including moment bounds and tail behavior within clusters. Specifically, any normalized linear combination of errors within a cluster (i.e., one with unit variance) must have a uniformly bounded fourth moment. Likewise, normalized combinations of $e_i v_i'$ must have uniformly bounded $(2 + 2\delta)$ moments. These assumptions ensure that intra-cluster dependence is well-approximated by the covariance matrix and rule out extreme tail dependence, such as aligned outliers not reflected in the covariance. Further, we require that certain correlations are bounded away from one and that variances are bounded away from zero, thereby excluding cases of near-perfect predictability.
The second set (Assumption (ref)) concerns the interaction between the weight matrix $A$ and the operator norms of the cluster-specific covariance matrices $\Sigma_i$. These assumptions ensure that each cluster's contribution to the total variance remains asymptotically negligible. By working with $\|\Sigma_i\|$, we allow for a trade-off between cluster size and intra-cluster dependence. Two stylized examples clarify the implications for the quadratic form (Assumption (ref)(ref)) that typically requires stronger restrictions :
Let $\Sigma_i$ denote the covariance matrix of $(e_i', v_i')'$, with submatrices $\Sigma_{e,i}$ and $\Sigma_{v,i}$. Define the normalized shock vector for cluster $i$ as $\xi_i = \Sigma_i^{-1/2} (e_i', v_i')'$, which is a $2T_i \times 1$ vector. Also define the cross-covariance matrix $\Psi_i = \mathbb{E}[e_i (v_i \otimes e_i)']$ of dimension $T_i \times T_i^2$, and the variance matrix $\Phi_i = \mathbb{V}(v_i \otimes e_i)$ of dimension $T_i^2 \times T_i^2$. \looseness=-1
The assumptions above ensure two linear CLTs: one for sums of $e_i$ and one for sums of $v_i \otimes e_i$. Condition (ref) guarantees that $(\Sigma_{v,i}^{-1/2} v_i) \otimes (\Sigma_{e,i}^{-1/2} e_i)$ has eigenvalues bounded away from zero, preventing near-degeneracy. Condition (ref) ensures that correlations between $e_i$ and $v_i \otimes e_i$ remain bounded away from one. Condition (ref) rules out heavy-tailed behavior and strong higher-order dependence that could dominate the covariance-based approximation.\looseness=-1
Part (ref) guarantees that the normalization in the CLT is no slower than $1/\sqrt{n}$. Part (ref) ensures cluster negligibility for the linear CLT, Part (ref) (b) does so for the quadratic CLT, while Part (ref) (a) assumes that the quadratic part is negligible in comparison to the linear one. \looseness=-1
We now discuss how Assumption (ref) applies to the matrix $A^*$ defined in Theorem (ref). According to Lemma (ref), the Frobenius norm $\|A^*\|_F^2$ captures the effective sample size. Assumption (ref)(ref) holds provided $A^*$ retains sufficient variation relative to $n$. From Theorem (ref)(c), each row of $A^*$ corresponds to a projection, and we have: $$\sum_{\tilde\ell=1}^n(A_{\ell\tilde\ell}^*)^2=A_{\ell\ell}^*\mbox{ and }\|A^*\|_F^2=\operatorname{tr}(A^*)=\sum_{\ell}A_{\ell\ell}^*.$$ Since $A^*$ is observable, $\|A^*\|_F^2$ can be computed and directly compared to $n$. A simple sufficient condition for Assumption (ref)(ref) is illustrated below:
Theorem (ref) shows that $A^* = (I - B)M$, where $B$ is a sparse matrix with relatively few nonzero entries. Assumption (ref) ensures that $A^*$ remains close to $M$ under appropriate matrix norms.\looseness=-1
For linear CLTs, the condition $\max_i T_i^2 / n \to 0$ suffices hansen2019asymptotic, sasaki2022non, while the quadratic CLT (Assumption (ref)(ref)(b)) requires the stricter condition $\max_i T_i^4 / n \to 0$.\looseness=-1
In high-dimensional settings with many regressors, an analog of the second part of Assumption (ref)(ref)—e.g., $\frac{\|M \lambda\|_\infty^2}{n} \to 0$—is often imposed as a high-level assumption because direct verification is difficult; see, for example, Assumption 3 in cattaneo2018inference. Lemma (ref) offers low-level primitive conditions that imply this assumption. The bound $\|\lambda\|_\infty<\infty$ is mild since $\lambda_\ell = \mathbb{E}[x_\ell]$, and for a fixed effects projection matrix, we have $\|M\|_\infty = 2$.\looseness=-1
Suppose we test the null hypothesis $H_0: \beta = \beta_0$, and define $U_i = y_i - x_i' \beta_0$. Even under $H_0$, these differ from the true errors $e_i$ due to the inclusion of a predictable component: $U_i = e_i + w_i' \delta$. By the partialling-out property of $A^*$, we have $x'A^*U = x'A^*e.$ \looseness=-1
The jackknife variance estimator, introduced by efron1981jackknife for symmetric statistics based on independent and identically distributed data, has also been advocated as a useful tool in settings involving non-symmetric statistics and non-identically distributed data bousquet2004concentration,mackinnon2023fast. In our setting, we define the statistic of interest as $\mathcal{Z} = x'A^*U = x'A^*e = f(D_1, \dots, D_N)$, where $D_i = (x_i, U_i)$ denotes the observed data for cluster $i$. Let $\mathcal{Z}_{(i)} = f(D_1, \dots, D_{i-1}, 0, D_{i+1}, \dots, D_N)$ denote the statistic computed when cluster $i$ is removed (i.e., replaced by zero). The jackknife estimator of the variance is then given by \[ \hat V_\textrm{JK} = \sum_{i=1}^N (\mathcal{Z} - \mathcal{Z}_{(i)})^2 = \sum_{i=1}^N \left(\sum_{j \ne i} x_j' A^*_{ji} U_i + \sum_{j=1}^N x_i' A^*_{ij} U_j\right)^2. \] \looseness=-1
The following lemma establishes that $\hat V_\textrm{JK}$ is a conservative estimator, potentially overestimating the true variance.
efron1981jackknife showed that the jackknife variance estimator tends to be conservative for symmetric statistics derived from independent and identically distributed data. While the jackknife correctly captures the variance of first-order $U$-statistics, it double-counts second-order terms. In our setting, $\mathcal{Z}$ is non-symmetric and based on independent but non-identically distributed data. The jackknife is conservative for two main reasons: it gives double weight to the quadratic term (as noted in efron1981jackknife), and $\mathcal{Z}_{(i)}$ is improperly centered since $\mathbb{E}[U_i - e_i] \neq 0$. The jackknife variance estimator is unbiased when $A^*_{ij} = 0$ for all $i \ne j$, in which case the removal of a single cluster does not affect the residualization of others. Consequently, $\mathcal{Z}_{(i)}$ is properly centered, and no cross-cluster covariances arise. In general, the degree of conservativeness is small when $\sum_{j=1}^N\sum_{i\neq j}\|A_{ij}^*\|_F^2$ is much smaller than $\sum_{j=1}^N\|A_{jj}^*\|_F^2$.\looseness=-1
Our estimator, $\hat\beta^{A^*} = \frac{x' A^* y}{x' A^* x} = \frac{z' y}{z' x}$, belongs to the class of just-identified internal IV estimators. Similar to the IV estimator for dynamic panels proposed by arellano1991some, our approach may be subject to weak identification concerns (see bun2010weak). To construct inference procedures that are robust to weak identification, we propose using the Anderson–Rubin (AR) test. The AR test delivers valid inference under both strong and weak identification.
To test the null hypothesis $H_0: \beta = \beta_0$, we employ the test statistic \[ AR(\beta_0) = \frac{(x' A^* (y - x \beta_0))^2}{\hat V(\beta_0)}, \] and reject the null if this statistic exceeds the $1 - \alpha$ quantile of the $\chi^2_1$ distribution. The variance estimator $\hat V(\beta_0)$ is given by $\hat V_\textrm{JK}$, which is based on the implied errors $U = U(\beta_0) = y - x \beta_0$. By Theorem (ref), the AR test has asymptotic size no greater than the nominal level. \looseness=-1
A confidence set is obtained by inverting the AR test. Since both the numerator and the denominator are quadratic in $\beta_0$, the inversion reduces to solving a quadratic inequality. The resulting confidence set is always non-empty and always includes the estimator $\hat\beta^{A^*}$. \looseness=-1
When identification is strong, the AR-based test and confidence set are asymptotically equivalent to those based on the $t$-statistic, using $\hat\beta^{A^*}$ as the point estimator and $\frac{\hat V(\hat\beta^{A^*})}{(x' A^* x)^2}$ as the squared standard error. Nevertheless, we advocate using the AR test, as the validity of $t$-statistic-based inference is not generally guaranteed.
Policy interventions with spatial effects present a compelling setting for applying our method. In such contexts, researchers often distinguish between the direct effect of treatment---on the treated units themselves---and the total effect, which includes indirect impacts arising from spillovers. A prominent example is provided by egger2022general (hereafter EHMNW), who study a large-scale fiscal intervention in rural Kenya. In their experiment, one-time cash transfers in total equivalent to approximately 15 percent of local GDP were randomly assigned to 653 villages. We apply our method to this setting to estimate the direct effect of the intervention on treated villages.
The villages (our units of observation) were divided into 68 sub-locations (clusters), the median sub-location contains $7$ villages. Treated villages were randomly selected through a two-stage randomization procedure. Within each treated village, all eligible households received a one-time cash transfer of USD 1,871 (PPP-adjusted). Eligibility was determined through a poverty-based means test. The median village includes $98$ households, approximately one-third of whom met the eligibility criteria. EHMNW examine a range of outcomes, including household consumption, firm behavior, and local prices. Here, we focus exclusively on village-aggregated outcomes for total household consumption. Results for other outcomes are similar and available upon request.
The parameter of interest in our application is the village-level direct treatment effect on average household consumption. Although treatment was randomly assigned, a simple OLS regression is unreliable due to anticipated interference from neighboring villages. EHMNW assume that spatial spillovers do not extend beyond a radius of $R$ kilometers, with their preferred specification (selected through BIC criteria) being $R = 2$.
This setting fits naturally within our design-based framework, which we illustrate through estimation under two specifications. In both cases, we define the exogeneity matrix $\mathcal{E}$ using cutoffs based on pairwise distances between villages. We assume the exclusion restriction $E[v_{\tilde\ell}({y}_\ell - \beta x_\ell)] = 0$ holds if and only if $\ell=\tilde\ell$ or the distance between villages $\tilde\ell$ and $\ell$ is at least $R$ kilometers. In practice, researchers should choose a distance cutoff below which exogeneity may be questionable. Here, we illustrate the behavior of the $\hat\beta^{A^*}$ estimator using a range of cutoffs from 0 to 3 kilometers. In the data, the average pairwise distance between villages within a sub-location is $2.2$ kilometers. On average, each village has $1.5$, $4.6$, and $6.7$ neighbors within the sub-location and within 1, 2, and 3 kilometers\footnote{Based on GPS coordinates. We dropped nine observations for which we do not have GPS coordinates.}, respectively. Our estimates under different choices of $R$ are shown in Figure (ref).
For specification (a), let ${y}_\ell$ denote the average consumption of households in village $\ell$, and let $x_\ell$ denote the binary treatment status of the village. We assume the correctly specified treatment equation: $x_\ell = p_{i(\ell)} + v_\ell,$ where $p_i$ is a cluster (sub-location) fixed effect capturing treatment propensity. We include sub-location fixed effects to flexibly absorb the realized treatment share, which may differ substantially from the ex ante assignment probability. In practice, the realized propensity can depend on cluster size or arise from re-randomization or re-balancing procedures implemented during the experiment.
The results are reported in Panel (a), alongside the baseline estimator from the specification used in EHMNW. The baseline estimator is reported in Table 1 (row 1, column 1) of EHMNW. We find that the estimates remain relatively stable for distance cutoffs below 2 km, suggesting that the bias in estimating the direct effect due to inter-village spillovers is limited. As larger values of $R$ correspond to more relaxed exogeneity assumptions, the effective sample size decreases, leading to larger standard errors. Figure (ref), Panel (a), also displays confidence intervals based on the jackknife variance estimator described in Section (ref). In this specification, since the only controls are cluster fixed effects, the $A^*$ matrix is block-diagonal, and the jackknife estimator coincides with the standard cluster-robust variance estimator.
Specification (b), reported in Figure (ref), Panel (b), makes use of the full variation in treatment intensity. Here, we define ${y}_\ell$ as the average consumption of eligible households in village $\ell$, and let $x_\ell$ denote the total transfer allocated to the village. The assumed treatment assignment equation in this case is $x_\ell = W_\ell' \delta + v_\ell,$ where $W_\ell$ includes the number of eligible households in village $\ell$ along with sub-location fixed effects. The parameter of interest, $\beta$, again represents the village-level direct treatment effect. This parameter is somewhat comparable to the total effect on eligible households reported in EHMNW, which was estimated using an instrumental variables strategy and reported in Table 1 (row 1, column 2) of the paper. Estimates under different exogeneity assumptions are presented along with both jackknife and cluster-robust standard errors. In this specification, the controls $W_\ell$ include a covariate that varies within clusters and is not absorbed by the sub-location fixed effects. As a result, matrix $M$ is no longer block-diagonal, and neither is $A^*$, so the two variance estimators do not coincide.
\paragraph{Discussion.} Both panels illustrate our main takeaway: the point estimates—and especially their precision—are highly sensitive to the researcher’s maintained exogeneity assumptions. Relaxing these assumptions, for example by allowing spillovers to extend within 3 km rather than 2 km, leads to substantially wider confidence intervals. This reflects a reduction in the effective sample size captured by the trace of the $A^*$ matrix under less restrictive exogeneity assumptions.
The structure of the matrix $A^*$ in the continuous treatment specification (same as in Figure (ref), Panel (b)) is shown in Figure (ref), Panel (a) for $R = 1$ and Panel (b) for $R = 3$. A blue dot indicates that the absolute value of the corresponding element of $A^*$ exceeds 0.01, while a red dot indicates that it exceeds 0.1. Villages are sorted by cluster (sub-location), so that block structure, if present, is visible. The figure reveals that the matrix $A^*$ is far from block-diagonal: some villages contribute substantially to the residualization of controls across many clusters. As the cutoff radius $R$ increases, the trace of the matrix $A^*$ decreases. This pattern reflects a reduction in the effective sample size. The ratio of Frobenius norms for off-diagonal blocks to diagonal blocks increases with distance, but remains small. When $R=1$, we have $\operatorname{tr}(A^*)=539$, while the relative Frobenius norm of off-diagonal blocks is $0.043$. When $R=3$, $\operatorname{tr}(A^*)=273$ while the Frobenius norm of off-diagonal blocks is $0.048$. Due to the relative smallness of off-diagonal blocks, the confidence sets using clustered and jackknife variance estimators are very close to each other.