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.
84,384 characters · 15 sections · 67 citation commands
Cluster-Robust Inference for Quadratic Forms
\setcounter{oldtocdepth}{\value{tocdepth}} \addtocontents{toc}{\setcounter{tocdepth}{-10}}
We study inference for quadratic forms of linear regression coefficients from two regressions with clustered data, given by
where $g$ indexes clusters, $i$ indexes observations within a cluster, and $W_{i, g}$ is a $d$-dimensional vector of covariates, which is treated as fixed. The object of interest,
combines the two sets of regression coefficients using a known nonrandom $d\times d$ matrix $A_0$. Two key features of this setting complicate inference. First, the covariates may be high-dimensional, with $d$ allowed to grow as fast as the sample size. Second, the error terms $(U_{i, g}, V_{i, g})$ may be heteroskedastic and exhibit arbitrary within-cluster correlation, with cluster sizes that may diverge with the sample size.
This framework covers several important applications: if we set $Y=X$, variance components can be written as quadratic forms in regression coefficients KSS2020; tests of many linear restrictions can likewise be cast as a restriction on a quadratic form AS23. The framework also covers inference in \ac{IV} regression models with many potentially weak instruments or controls, and heterogeneous treatment effects---in this case (ref) correspond to reduced-form and first-stage equations, $W$ contains both IVs and controls, and the effect of $X$ on $Y$ may be heterogeneous.
In the IV setting, $W$ is commonly high-dimensional in the popular IV leniency design, which leverages quasi-random assignment of judges or other decision-makers who differ in their leniency---the propensity to grant treatment CFL2024,GHK25. Since these leniency measures are unknown, they must be estimated in a first-stage regression of treatment on decision-maker indicators. Furthermore, since their assignment is typically only random conditional on time and location, time-by-location fixed effects need to be included. Commonly, these designs feature a large number of decision-makers and fine-grained fixed effects, resulting in a high-dimensional first-stage regression. More generally, to avoid restrictive first-stage assumptions under nonrandom instrument assignment, ensuring causal interpretation requires the inclusion of interactions between the instruments and controls, as well as a flexible control specification BBMT22. This can generate a high-dimensional first-stage regression even if the number of baseline instruments and covariates is small. Whenever cases are assigned in batches, this induces clustering in the data.
Similarly, variance component estimation commonly features high-dimensional specifications. For instance, consider the popular two-way worker--firm fixed effects model of AKM99, where $Y$ denotes log wages, $X=Y$, and $W$ contains both worker and firm dummies. In this setting, firms' contributions to wages, measured by the variance component of firm effects, can be formulated as a quadratic form in firm fixed effects, and the sorting of high-wage workers to high-wage firms can be measured by a quadratic form involving worker and firm fixed effects. Unless one has access to a long panel, the number of workers is necessarily large relative to the sample size; often the firm dimension is large as well. To allow for match-specific unobservables, it is desirable to cluster at the worker--firm match level kline2024. Consequently, valid inference must account for both the high dimensionality of the regressors and the clustered dependence structure of the data.
Finally, tests of linear restrictions, such as testing whether a subset of the $\gamma$ coefficients is zero, feature high-dimensional specifications in several applications. For example, value-added models can be validated by testing whether a set of regressors is excluded from the structural regression model ahpw24. Likewise, tests for endogenous peer effects can be cast as testing whether peer characteristics have zero coefficients JL26. As discussed in AS23, testing for the presence of heterogeneity amounts to testing whether a set of fixed effects is zero krw22,fgw16. These zero restrictions can be high-dimensional, and it is desirable to cluster the standard errors to allow for unobserved common shocks: depending on the application, one may want to cluster at the job level krw22 or by classroom (in value-added or peer-effects applications).
Although least squares estimators of linear regression coefficients are unbiased under general conditions, their quadratic forms are biased (the bias is positive by Jensen's inequality if $A_{0}$ is positive semidefinite). The bias scales with the coefficient dimension, so that plugging in least squares estimates of $\gamma$ and $\pi$ into (ref) yields an estimator with non-negligible asymptotic bias if the coefficient dimension is large. In the IV setting, this bias corresponds to the well-known many instrument bias of two-stage least squares bjb95,bekker94. The bias can be purged by using a \ac{L1CO} estimator proposed by KSS2020.
The contribution of this paper is threefold. First, we establish the asymptotic normality of the \ac{L1CO} estimator. Second, we develop a \ac{L3CO} cluster-robust variance estimator and derive conditions for its consistency. Third, we introduce a novel \ac{L2CO} variance estimator that yields valid, but conservative inference under a weaker set of assumptions. These results are derived under primitive conditions that allow for growing cluster sizes, both weak and strong within-cluster dependence, and only impose weak conditions on the regressors, allowing their dimension $d$ to grow potentially as fast as the sample size $n$. In particular, under suitable further regularity conditions, the proposed inference based on the \ac{L3CO} variance estimator is valid provided that the largest cluster size \(n_G\) satisfies \(n_G^{\max\left(\frac{2q}{q-1},3\right)} = o(n)\) and \(n_G^{5} = o(n)\) under weak and strong within-cluster dependence, respectively, where \(2q\) denotes the number of finite moments of the regression errors \((U_{i,g}, V_{i,g})\) for some $q\geq 2$. In contrast, the rate requirements associated with the \ac{L2CO} variance estimator can be relaxed to \(n_G^{\frac{2q}{q-1}} = o(n)\) under weak within-cluster dependence. We also show that when $d$ is not very large and cluster sizes are fixed, inference based on the \ac{L2CO} variance estimator is exact.
These results generalize and unify existing results from two strands of literature. The first strand studies inference in high-dimensional linear regressions. \Textcite{KSS2020} propose the \ac{L1CO} estimator we study. However, their formal results focus on the case with independent data, and their variance estimator is based on sample splitting. While least-squares estimators are unbiased for linear functions of regression coefficients, estimating their asymptotic variance boils down to estimating a quadratic form with a particular $A_{0}$ matrix. The Eicker-Huber-White estimator is a plug-in estimator of this quadratic form, and is thus biased in high-dimensional settings. \Textcite{CJN18_ET,CJN18} show consistency of a Hadamard-type estimator of the quadratic form, first proposed by hrk69, for the case with independent data or clustered data with bounded cluster sizes. \Textcite{CGJN22} extend its consistency to settings with diverging cluster sizes. In this context, under independent sampling, J22 establishes the consistency of the \ac{L1CO} unbiased estimator we study, and AnNg26 generalize the consistency result to the case with clustering.
In contrast, this paper is concerned with inference on quadratic forms, not just their consistent estimation. The variance of the \ac{L1CO} estimator we study depends on products of second moments, and can be written as a quartic form. As a result, unbiased variance estimation using leave-out methods requires leaving three clusters out. \Textcite{MSJ25} study estimation in linear models with high-dimensional controls, clustered data, and weak exogeneity. The weak exogeneity motivates formulating the estimation problem as a quadratic form problem, and they use this representation to construct a just-identified internal-IV estimator. They propose a jackknife variance estimator and show that its expectation is conservative.
The \ac{L3CO} variance estimator we study is inspired by AS23, who provide primitive sufficient conditions for its consistency in the special case with independent data, and $Y=X$. Extending these results to clustered data is substantially more involved, as it requires handling multilevel sums of products of block matrices, where matrix multiplication is non-commutative. We address this challenge through two key technical innovations: (1) a new representation of the L3CO projection matrix, and (2) a decomposition of sums of products of two L3CO projection matrices, each leaving out different clusters. These innovations yield primitive consistency conditions for the L3CO variance estimator in terms of the sample size $n$, the largest cluster size $n_{G}$, the number of moments of the residuals $2q$, the regressor dimension $d$, and the properties of $A_{0}$. The \ac{L2CO} variance estimator that we consider is, to our knowledge, new.
The second strand concerns inference in linear IV models with many weak instruments; see MS24 for a recent survey. Among studies assuming independent data, our results are most closely related to Yap24, who shows that ignoring the additional variation generated by $V_{i, g}$ under heterogeneous treatment effects can lead to underestimation of the asymptotic variance, and proves consistency of the L3CO variance estimator we study under high-level assumptions. The \ac{L1CO} unbiased estimator we study can be thought of as a generalization of the UJIVE estimator in K13 to the clustered setting. \Textcite{BN24,EK2018} also study IV inference with heterogeneous treatment effects, but focus on independent data.
We are only aware of a few IV studies that allow for clustering. \Textcite{CNT23} study many-IV regression under clustering but impose strong restrictions, including a correctly specified structural equation (homogeneous treatment effects), few controls (of order $o(\sqrt{n})$, excluding cluster fixed effects), bounded cluster sizes, and independence of errors within clusters. \Textcite{FLM23} allow within-cluster dependence and propose a cluster-jackknife estimator, but they do not accommodate heterogeneous treatment effects, weak identification, or many controls, and they do not provide a consistent cluster-robust variance estimator for high-dimensional settings. \Textcite{L23} considers the Anderson--Rubin test in many-IV regression under clustering, and LW24 extend it to multidimensional clustering, adapting bias corrections from CNT23 but without formal distributional theory when there are many controls.
We sample data from $G$ clusters indexed by $g\in[G]$, where $[G]$ is a shorthand for $\{1, \dotsc, G\}$. Each cluster contains $n_g$ units indexed by $i\in [n_g]$, such that the total sample size is given by $n=\sum_{g \in [G]}n_g$. Without loss of generality, we assume the last cluster is the largest, so that $n_{G}$ gives the largest cluster size. For each unit, we observe a vector of covariates $W_{i, g} \in \Re^{d}$, and two scalar outcome variables, $Y_{i, g} \in \Re$, and $X_{i, g} \in \Re$. We treat the covariates as deterministic (or, equivalently, we condition on them), and assume the regressions of $Y_{i, g}$ and $X_{i, g}$ onto $W_{i, g}$, defined in (ref), are linear, so that the error terms $U_{i, g}$ and $V_{i, g}$ are mean zero. The matrix $W_g \in \Re^{n_g \times d}$ stacks the row vectors $\{W_{i, g}'\}_{i \in [n_g]}$, and the matrix $W$ stacks the matrices $\{W_g\}_{g\in[G]}$. The matrices $Y$, $X$, $U$, $V$, as well as their cluster-specific counterparts $Y_g$, $X_g$, $U_g$, $V_g$ are defined analogously.
A natural estimator of the quadratic form $\theta$, defined in (ref), is the plug-in estimator that replaces $\pi$ and $\gamma$ with the \ac{OLS} estimates $\hat{\pi}=(W' W)^{-1}W' X$ and $\hat{\gamma}=(W' W)^{-1}W' Y$,
This plug-in estimator displays an overfitting bias given by
The bias arises because the estimation error in $\hat{\pi}$ is correlated with $Y_{g}$. If we replace $\hat{\pi}$ in (ref) with an \ac{OLS} estimator that leaves out cluster $g$, $\hat{\pi}_{-g}=(W'W-W_{g}'W_{g})^{-1}(W'X-W_{g}'X_{g})$, we obtain the leave-out estimator
that is unbiased by construction. This estimator was first proposed by KSS2020 in the setting where $X=Y$, though their formal results focus on the setting with independent data. As discussed in KSS2020, the leave-out estimator can alternatively be thought of as a debiased version of the plug-in estimator. In particular, let \(P = W (W'W)^{-1} W'\) and \(M = I_n - P\) denote the projection and annihilator matrices, respectively. Let \(P_{g,h} \in \mathbb{R}^{n_g \times n_h}\) and \(M_{g,h} \in \mathbb{R}^{n_g \times n_h}\) denote the submatrices of \(P\) and \(M\) corresponding to rows and columns associated with clusters \(g\) and \(h\). We may then equivalently write
where, $B = A - M D$ and $D = \operatorname{Bdiag}(M_{g,g}^{-1}A_{g,g}) \in \Re^{n \times n}$ is the block diagonal matrix with blocks $M_{g, g}^{-1}A_{g, g}$ on its diagonal (so that the diagonal blocks $B_{g, g}$ are all zero). Here, $A_{g, g}$ is defined analogously to $M_{g, g}$, and $(MX)_g \in \Re^{n_g}$ denotes the subvector of $MX$ corresponding to cluster $g$. (ref) follows from the definition of $\hat \theta_{\rm LO}$ since, by the Woodbury identity, $\hat{\pi}_{-g}=\hat{\pi}-(W'W)^{-1}W_{g}'M_{g, g}^{-1}(MX)_{g}$. Since
(ref) implies that $\hat{\theta}_{LO}$ obtains by subtracting off a bias estimate from the plug-in estimator.
The bias correction is not unique: any estimator of the form $X'(A-C)Y$, where $C$ is a matrix such that $W'CW=0$ and $\operatorname{Bdiag}(C_{g,g})=\operatorname{Bdiag}(A_{g,g})$ will be unbiased. To motivate the bias correction matrix $C_{\rm LO}=M\operatorname{Bdiag}(M_{g,g}^{-1}A_{g,g})$ we consider, (ref) below generalizes the classic results in rao70 showing that $C_{\rm LO}$ has certain optimality properties. To state the result, let $Q\ast S$ denote the Khatri-Rao product of two matrices $Q,S$, so that the $(g, h)$ block of the matrix is given by the Kronecker product $(Q\ast S)_{g, h}=Q_{gh}\otimes S_{gh}$. Under independent sampling, the Khatri-Rao product reduces to the Hadamard product, $Q\ast S=Q\odot S$. Let $\operatorname{bvec}(Q)=(\operatorname{vec}(Q_{1,1})',\dotsc,\operatorname{vec}(Q_{G,G}))'$ vectorize the diagonal blocks of an $n\times n$ matrix $Q$.
To interpret this result, note that property (i) is necessary for unbiasedness, while properties (ii) and (iii) ensure that the estimator is invariant to the signals $W\gamma$ and $W\pi$, respectively, in the regressions in (ref). If these conditions hold and the errors $U,V$ are homoskedastic, then the variance of the bias correction is proportional to the square the Frobenius norm, $\text{tr}(C'C)$. Thus, the bias-correction $C_{KR}$ is variance-minimizing among all bias-corrections that are invariant to the signal in both regressions. The resulting estimator $\hat{\theta}_{\rm KR}=X'(B-C_{\rm KR})Y$ (or its version under independent sampling) has been previously studied by, among others, hrk69,CJN18,CGJN22,CNT23,L23. If the invariance requirement is weakened, so that we only require the bias-correction to be invariant to the signal in one of the regressions, the minimum norm unbiased estimator is given by $C_{\rm LO}$ or $C_{\rm LO}'$, which yields the leave-out estimator in (ref) that this paper focuses on (using $C_{\rm LO}'$ yields a numerically equivalent expression if $X=Y$, but not in general).
While the invariance properties of the estimator $\hat{\theta}_{KR}$ are clearly desirable, they come with three limitations relative to the estimator $\hat{\theta}_{LO}$. First, since $\text{tr}(C_{\rm LO}'C_{\rm LO})\leq \text{tr}(C_{\rm KR}'C_{\rm KR})$, it can display greater variability. Second, computing the estimator $\hat{\theta}_{KR}$ is not feasible in large datasets or in datasets with large clusters, as it requires solving a linear system defined by the matrix $M\ast M$, which has dimension $\sum_{g} n_{g}^2$, and may not even be storable in memory. In contrast, the leave-out estimator only requires computation of least squares estimators. Third, the estimator may not exist. A sufficient condition for the existence of $\hat{\theta}_{KR}$ is that the matrix $M\ast M$ be invertible. Under independent sampling, the matrix $M\ast M$ corresponds to the Hadamard product $M\odot M$ (so that $(M\ast M)_{ij}=M_{ij}^{2}$), and HoHo75 show that a sufficient condition for its invertibility is that $\min_i M_{ii} > 1/2$, which holds if $d < n/2$ and the design is sufficiently balanced. While this is only a sufficient condition, simulations reported in J22 show that, with independent sampling, the estimator $\hat{\theta}_{KR}$ often fails to exist when the regression is high-dimensional. Under clustered sampling, general sufficient conditions for the invertibility of the Khatri-Rao product are, to our knowledge, unknown. For these reasons, we focus on the leave-out estimator.
To test a particular null hypothesis $\theta=\theta_0$ at the significance level $\alpha$, we reject if the $t$-statistic based on the leave-out estimator exceeds the usual critical value ${z}_{1-\alpha/2}$, the $1-\alpha/2$ quantile of a standard normal distribution. That is, we reject if
where $\hat \omega_n^2$ is either a consistent or a conservative estimator of the variance of the estimator $\omega_n^2 = \mathbb V(X' B Y)$. We consider variance estimation in (ref); asymptotic validity of this test follows from the fact that $\hat{\theta}_{\rm LO}$ is asymptotically normal, as shown in (ref). Before giving these results, the remainder of this section fleshes out three important special cases that fit this general setup.
To see how our setup covers instrumental variables regressions, decompose the covariate vector $W_{i, g}=(\mathcal W_{i, g}', Z_{i, g}')'$ into a vector of controls $\mathcal{W}_{i, g}\in \Re^{d_w}$ and a vector of instruments $Z_{i, g}\in \Re^{d_z}$, with $d= d_w + d_z$. We allow both $d_w$ and $d_z$ to diverge to infinity, potentially as fast as $n$. Let $\mathcal{Y}_{i, g}$ denote an outcome variable of interest, with $X_{i, g}$ denoting the treatment. Letting
denote the reduced-form regression and (ref) the first-stage regression, the \ac{TSLS} estimand is given by
where $P_{\tilde Z}=\tilde Z(\tilde Z' \tilde Z)^{-1}\tilde Z'$ is the projection matrix of $\tilde Z$, $\tilde Z = M_{\mathcal W} Z$, and $M_{\mathcal W} = I_n - \mathcal W (\mathcal W' \mathcal W)^{-1}\mathcal W'$ is the annihilator matrix of $\mathcal W$. For testing a particular value $\beta_0$ of $\beta$, let $Y=\mathcal{Y}-X\beta_0$ denote the treatment-adjusted outcome, so that (ref) corresponds to the treatment-adjusted outcome regression with $\gamma=\pi_Y -\pi\beta_0$ and $U_{i, g}=e_{i, g}-V_{i, g}\beta_0$. Our setup allows for testing the null value $\beta_0$ by testing whether $\pi'A_0 \gamma=0$, with
Hypothesis tests on $\beta$ are of interest because if the effect of $X_{i, g}$ on $\mathcal{Y}_{i, g}$ is constant, and the instrument satisfies an exclusion restriction, then $\beta$ identifies this constant treatment effect, provided that the regression of $\mathcal{Y}_{i, g}-X_{i, g}\beta$ on $\mathcal{W}_{i, g}$ is linear. Our setting allows for heterogeneous treatment effects and, in this case, a causal interpretation of $\beta$ requires an instrument monotonicity assumption, and the assumption that the reduced-form and first-stage regressions are (approximately) correctly specified. We refer to K13,EK2018,BN24,Yap24,BBMT22 for further discussion and statement of the precise conditions. The requirement that the reduced-form and first-stage regressions be approximately linear gives one motivation for why the covariate vector $W_{i, g}$ may be high-dimensional even in settings where a baseline set of instruments $Q_{i, g}$ and covariates $C_{i, g}$ is low dimensional (as discussed in the introduction, in leniency IV designs, these vectors may already be high-dimensional). In particular suppose we set $\mathcal W_{i, g} = \mathcal W(C_{i, g})$, and $Z_{i, g} = Z(Q_{i, g}, C_{i, g})$ to correspond to technical transformations of these baseline variables ensuring that the linear specifications $W_{i, g}'\pi_Y$ and $W_{i, g}'\pi$ provide good approximations to the reduced-form $\mathbb{E}(\mathcal Y_{i, g} \mid Q_{i, g}, C_{i, g})$ and first stage $\mathbb{E}(X_{i, g} \mid Q_{i, g}, C_{i, g})$, respectively. Typically, such technical transformations involve interactions and series expansions, which can lead to high dimensionality. Our analysis remains valid if the model in (ref) contains an asymptotically negligible approximation error, as in, for example, CJN18. For notational simplicity, we abstract from any such approximation errors and assume that the linearity holds exactly.
The variance $\omega_n^2$ for $X' B Y$ and its estimator $\hat{\omega}_n^2$ implicitly depend on $\beta_0$ through the relation $Y = \mathcal{Y} - X \beta_0$, so the resulting $t$-statistic corresponds to the weak-identification-robust Lagrange multiplier statistic. Our regularity conditions require that either the concentration parameter $\pi' W' P_{\tilde{Z}} W \pi$ or the number of instruments $d_z$ diverges to infinity. This framework encompasses both the strong-identification case with a finite number of instruments and the strong- and weak-identification cases with many instruments, including situations with many weak IVs in which $\pi' W' P_{\tilde{Z}} W \pi / \sqrt{d_z}$ remains bounded.
If the instruments are collectively strong, so that $\pi' W' P_{\tilde{Z}} W \pi / \sqrt{d_z}$ diverges, the leave-out estimator of $\beta$,
will be consistent, even in the presence of many instruments (this stands in contrast with the inconsistency of TSLS, which can be thought of as a plug-in estimator of $\beta$). In this case, we do not have to impose the null when computing the variance of $X'BY$, and can instead base inference on the Wald test, rejecting whenever the absolute value of the $t$-statistic
exceeds $z_{1-\alpha/2}$, where $\hat \omega_n^2(\hat \beta)$ is just $\hat \omega_n^2$ with $\beta_0$ in $Y = \mathcal Y - X \beta_0$ replaced by $\hat \beta$.
The construction of our estimator $\hat{\beta}$ (and its corresponding weak-IV robust Wald test) can be thought of as a cluster-robust version of the UJIVE estimator studied by K13 and further advocated by GHK25, which employs a leave-one-out technique to eliminate the overfitting bias of TSLS\@.
Consider a two-way fixed effect model of log wage determination proposed by AKM99, in which the log-wage $Y_{\ell, t}$ of a worker $\ell\in[L]$ in year $t \in [T_{\ell}]$ is given by
Here $\alpha_{\ell}$ is a worker fixed effect, $\mathcal J(\ell, t) \in [J]$ returns the identity of the firm that employs $\ell$ in year $t$, $\psi_{j}$ is a firm fixed effect, $\mathcal{W}_{\ell, t}$ contains time-varying control variables, and $U_{\ell, t}$ is a time-varying noise term. The error $U_{\ell, t}$ is assumed to be mean zero, which imposes strict exogeneity: workers and firms are allowed to match based on firm and worker fixed effects, but not on time-varying factors influencing wages. As discussed in a recent survey by kline2024, it is desirable to allow this noise to exhibit arbitrary correlation within a worker-firm match to allow for match-specific unobservables, so that $U_{\ell, t}$ and $U_{\ell, t'}$ may be correlated if $J(\ell, t)=J(\ell, t')$. To map this to our setup, following kline2024, let $g$ index worker-firm matches, with $j(g)$ returning the firm and $\ell(g)$ the worker that form the match. Then $n_{g}$ corresponds to the duration of the match (so if a worker $\ell(g)$ spends 5 years at a firm $j(g)$, say, then $n_{g}=5$). Then we may rewrite (ref) as
where $i\in[n_{g}]$ indexes years during the match, $D_{\ell}\in\Re^{L}$ is the $\ell$-th basis vector, and $F_{j}\in\Re^{J-1}$ is the $j$-th basis vector (we normalize the last firm's effect to zero). This maps to our general framework in (ref) by setting $Y_{i, g}=X_{i, g}$, $W_{i, g}=(F_{j(g)}, D_{\ell(g)}, \mathcal{W}_{i, g})$, and $\gamma=\pi=(\psi', \alpha', \delta')'$.
The key parameters of interest in this model are the variances of firm and worker effects, as well as their covariance. A person-year weighted variance of the firm effects is given by $\sigma^{2}_{\psi}=\frac{1}{n}\sum_{g\in[G]}n_{g}(\psi_{j(g)}-\bar{\psi})^{2}$, where $\bar{\psi}=\frac{1}{n}\sum_{g\in[G]}n_{g} \psi_{j(g)}$, and measures the direct contribution of firms to wage inequality. This can be written as a quadratic form in (ref), $\sigma_{\psi}^2 = \gamma'A_\psi \gamma$, with the matrix $A_{\psi}$ given by
where $\bar{F} = \frac{1}{n}\sum_{g \in [G]}n_{g}F_{g}$, and $1_{n}$ is a vector of ones. Variance of worker effects, given by $\sigma^{2}_{\alpha}=\frac{1}{n}\sum_{g\in[G]}n_{g}(\alpha_{\ell(g)}-\bar{\alpha})^{2}$, maps to (ref) analogously.
The covariance between worker and firm effects is given by $\sigma_{\alpha,\psi} = \frac{1}{n}\sum_{g \in [G]}n_{g}(\psi_{j(g)} - \bar\psi) \alpha_{\ell(g)}=\gamma' A_{\alpha,\psi} \gamma$, with
and $\bar{D} = \frac{1}{n} \sum_{g \in [G]}n_{g}{D}_{\ell(g)}$, so that $\sigma_{\alpha, \psi}$ again takes the form of (ref). The covariance $\sigma_{\alpha, \psi}$ measures the contribution of systematic sorting of high wage workers to high wage firms to wage inequality.
While we focus on the AKM99 model for concreteness, variance and covariance components are of interest in numerous other settings that exhibit potentially high-dimensional fixed effects as well. Examples include determining the importance of neighborhoods for intergenerational mobility ChHe18i, the importance of geography for healthcare utilization fgw16, or the importance of classroom assignment in determining student outcomes chetty2011. In some cases, covariance components involve fixed effects estimated from different regressions (so that $Y\neq X$): for instance, chetty2011 are interested in the covariance between classroom fixed effects in an earnings regression and that in a test score regression.
Consider the linear regression (ref), so that \( X_{i,t} = Y_{i,t} \) in (ref). As in AS23, we are interested in testing the linear restriction
where \( q \) may be high-dimensional. This restriction implies
which fits into our framework with \( Y = X \), and $A_0 = R' (R(W' W)^{-1} R')^{-1} R$, and the hypothesis of interest is that $\theta = q' (R(W' W)^{-1} R')^{-1} q$. Note that the plug-in estimator $X'A_{0}X$ corresponds to the numerator of the classic $F$-statistic under homoskedasticity.
By letting $q=0$ and letting $R$ correspond to a matrix that selects a subset of the coefficients $\gamma$, this setup nests testing for the presence of fixed effects, which is how tests of heterogeneity can be cast krw22,fgw16.
The setup also covers validation of value-added models. In particular, value added models are typically estimated by a linear regression specification, where $i$ indexes students within clusters $g$ (such as classrooms), and the vector of regressors $\mathcal{W}^{0}_{i, g}=(\mathcal{W}_{i,g}',\mathcal{D}_{i, g}')'$, consists of a set of school dummies $\mathcal{W}_{i, g}$ indicating which school $i$ attends, and a set of controls $\mathcal{D}_{i, g}$ (that may include lagged test scores). The outcome $Y_{i, g}=X_{i, g}$ denotes student test scores. Provided that $\mathcal{W}^{0}_{i, g}$ is exogenous, the coefficients on the school dummies (i.e., estimates of school fixed effects), can be interpreted as school value added estimates. \Textcite{ahpw24} point out that the exogeneity assumption is testable if we have a set of covariates $\mathcal{Z}_{i, g}$ that affect school attendance, but not test scores directly (such as indicators for winning a lottery to attend oversubscribed schools). \Textcite{ahpw24} then develop a test of exogeneity by adapting the classic Sargan test by viewing $\mathcal{Z}_{i, g}$ as a set of instruments. Their framework requires the instrument and covariate dimension to be fixed, and the sampling to be independent. However, the exogeneity assumption is equivalent to the assumption that the coefficients on $\mathcal{Z}_{i,g}$ in the long regression with $\mathcal{W}_{i, g}=(\mathcal{W}_{i,g}',\mathcal{D}_{i, g}',\mathcal{Z}_{i, g})'$ as the set of regressors all equal zero. This fits our framework by letting $R$ denote the selector matrix that selects the coefficients on $\mathcal{Z}_{i, g}$, and setting $q=0$. Doing so allows us to accommodate high-dimensional instruments and many schools, as well as clustering.
This framework also nests testing for endogenous peer effects. In particular, JL26 point out that in a panel data setting, one can test for peer effects without specifying the network structure by using an Anderson--Rubin test in an \ac{IV} regression of own outcomes on outcomes of potential peers, instrumenting with peer characteristics. However, implementing the test is complicated by the presence of individual and time fixed effects. The implementation in JL26 assumes the regression errors are homoskedastic and independent across both time and individuals. However, an Anderson--Rubin test of a zero effect of endogenous variables is equivalent to an $F$ test on the reduced form. Thus, the null hypothesis is equivalent to testing whether coefficients on the instruments in a reduced-form regression of own outcomes on controls (that may include individual and time effects) and instruments equal zero. This again fits the above framework, if we set $R$ to be a matrix that selects the regressor coefficients on the instruments. Doing so allows us to accommodates both heteroskedasticity and cluster dependence, either in the time dimension or in the cross-section (e.g., by classroom).
This section shows that the estimator $\hat{\theta}_{\textrm{LO}}$ is asymptotically normal. To state the assumptions needed for this result, let $\lambda_{\max}(Q)$ and $\lambda_{\min}(Q)$ denote the largest and smallest eigenvalues of a symmetric matrix $Q$, and for any matrix $Q$, let $\norm{Q}_{op}=\lambda_{\max}(Q'Q)^{1/2}$, and $\norm{Q}_{F}=\text{tr}(Q'Q)^{1/2}$. Finally, $\norm{q}_{2}$ denotes the Euclidean norm of a vector $q$.
To make use of (ref) for inference, we need a consistent or a conservative estimator of the asymptotic variance $\omega_{n}$. We now consider two such estimators. The first estimator we consider, a \acf{L3CO} estimator, is shown to be consistent, while the second one, a \acf{L2CO} estimator, is shown to be conservative.
The L3CO variance estimator we consider generalizes the leave-three-observations out estimators considered in AS23,Yap24 to the case with clustered data, and is defined as
where
Here $\tilde Y_{g,-hk} = Y_g - W_g \hat \gamma_{-ghk}$, and $\tilde X_{g,-hk} = X_g - W_g \hat \pi_{-ghk}$ denote leave-three-clusters-out OLS residuals,
denote leave-three-clusters-out OLS estimators, and
The definition of $\tilde M_{{g,k}, -gh}$ implies that when $g \neq h$,
and thus $\tilde M_{g,k,-gh}$ can partial out $W$ in the sense that $\sum_{k \in [G]} \tilde M_{g,k,-gh} W_k = 0_{n_g \times d}$. The key property of this variance estimator is that $\mathbb E \hat \omega^2_{n,\rm L3CO} = \omega_n^2$, which leads to consistency even when $d$ is large. To establish this formally, we impose the following assumption:
The next theorem establishes the unbiasedness and consistency for the L3CO variance estimator.
The \ac{L3CO} variance estimator may be computationally expensive in large datasets, since computing it involves looping over three cluster indices. This motivates an alternative \ac{L2CO} variance estimator, given by
where
denote leave-two-clusters-out OLS residuals (if $g=h$, we set $\tilde Y_{g,-h} = 0$ and $\tilde X_{g,-h} = 0$; the leave-two-clusters out estimators $\gamma_{-gh}$ and $\pi_{-gh}$ are defined analogously to (ref)). Since its computation only involves looping over two cluster indices, it is computationally cheaper than the \ac{L3CO} variance estimator. The second advantage of this estimator is that it is guaranteed to be non-negative. Third, as (ref) shows, valid inference based on $\hat{\omega}_{n, \rm L2CO}$ can be conducted under less restrictive assumptions. The inference is exact when $d$ is not very large and cluster sizes are fixed, but it can be conservative otherwise.
With large cluster sizes and/or high-dimensional covariates, the upward bias $\mathbb E \hat \omega_{n,\rm L2CO}^2 - \omega_n^2$ could be large. In such settings, we recommend using $\hat \omega^2_{n,\rm L2CO}$ only when computing $\hat \omega_{n,\rm L3CO}^2$ is infeasible.
In practice, it is not necessary to run OLS regression when computing $\hat \gamma_{-ghk}$ and $\hat \pi_{-ghk}$ for each combination of $g,h,k \in [G]^3$, which can be computationally demanding for large $d$. Instead, the leave-out algebra from (ref) implies that the leave-out residuals may be computed directly as
where $M_{ghk,ghk}$ is the block diagonal matrix of $M$ corresponding to clusters $(g,h,k)$. Similarly,
Another potential issue is the numerical instability when solving the system $M_{ghk,ghk}\tilde Y_{ghk,-hk} = (MY)_{ghk}$ that defines $\tilde Y_{g,-hk}$ above (and similarly for the other leave-out residuals). We propose the following approach: whenever the minimum eigenvalue of $M_{ghk,ghk}$ is below a threshold $t_n$, such that $t_n \downarrow 0$ as $n$ goes to infinity, replace it with a ridge regularizer, and instead solve $(M_{ghk,ghk}+t_{n}I_{n_{g}+n_{h}+n_{k}})\tilde Y_{ghk,-hk} = (MY)_{ghk}$. We find that $t_n = 1/\log(n^2)$ performs well in our simulation. Our theory can be adapted to deal with such shrinking regularizer, but for ease of exposition, we do not pursue this extension.
In this section, we compare our test statistic with L3CO variance estimator with several existing methods in the literature: the two-stage least squares estimator (TSLS) with cluster-robust variance estimator; the cluster jackknife instrumental variable estimator (CJIVE) proposed by FLM23 with cluster-robust variance estimator; and the three test statistics (CSW-JIV, CSW-LIM and CSW-FUL) proposed in CNT23. For all methods, we impose the null when computing their corresponding asymptotic variances to guard against weak identification. All the results below are based on $1,000$ simulations.
We first consider the following panel IV regression adapted from CNT23, which assumes a homogeneous treatment effect:
where $\alpha_g$ and $\xi_g$ are cluster-level fixed effects. The fixed effects are generated by $\alpha_g = u_{1g} + g/G, \xi_g = u_{2g} + g/G, g=1, \dots, G$ where $u_{1g}$ and $u_{2g}$ are independent standard normal random variables. These fixed effects are partialled out from the model by demeaning at the cluster level. The instruments in $Z$ are normally distributed with mean $0$ and we allow for cluster-level dependence: for each cluster the covariance matrix is given by
and across clusters these instruments are independent from each other; we set $\theta_1 = 0.5$ in our simulation. The controls in $W$ are generated by $\mathcal W_{i,g} = (z_{i,g}, z_{i,g}^2-1, z_{i,g}^3-3z_{i,g}, z_{i,g}^4-6z_{i,g}^2+3, z_{i,g}(D_{i,g}^{(1)}-0.5), \dots, z_{i,g}(D_{i,g}^{(d_w-4)}-0.5))$ where $z_{i,g}$ are independent standard normal random variables and $D_{i,g}^{(k)}, k=1,\dots,d_w-4$ are independent Bernoulli random variables with success probability $0.5$. The error terms are generated by
with $\sigma_{i,g} = \sqrt{(0.2+z_{i,g}^2)/2.4}$ and $\rho=0.5$, where $\varepsilon_{i,g}$, $\eta_{i,g}$ and $v_{i,g}$ are independent standard normal random variables. We further collect $\tilde{U}_{i,g}$ and $\tilde{V}_{i,g}$ for each cluster and left multiply them by
to generate $U_{i,g}$ and $V_{i,g}$; we set $\theta_2 = 0.7$ in our simulation. For the parameters, we set $\beta = 0.3$, $\gamma = \delta = (1/\sqrt{d_w}) \times \iota_{d_w}$ where $\iota_{d_w}$ is a $d_w \times 1$ vector of ones, and $\pi = t_n \times \iota_{d_z}$ where $\iota_{d_z}$ is a $d_z \times 1$ vector of ones and $t_n = \sqrt{30/(\sqrt{d_z} \times n)}$.
Finally, the number of observations is $n = 600$ and the number of clusters is $G = 150$, where all the clusters have equal cluster size. We consider three different designs for the number of instruments and the number of controls: $d_z = d_w = 50$, $d_z = d_w = 100$ and $d_z = d_w = 150$. We report the empirical size as well as the power curve for each design below.
We consider a simulation setup similar to that in Yap24 where we have many fixed effects, the effect of $X$ on $Y$ is heterogeneous, and the linear regressions for $Y$ and $X$ are correctly specified. Let $t = 1, \dots, K$ index the state and suppose that individuals within the same cluster belong to the same state. We have in total $K = 48$ states, where each state contains $4$ clusters and each cluster contains $4$ individuals. There is also another binary exogenous variable $B \in \{0,1\}$, which is equally distributed within each cluster. The structural equations are given by
where the control variables $\mathcal W_{i,g}$ contain indicators for states with $d_w = 48$ and the instrumental variables $Z_{i,g}$ contain indicators for $k = t \times B$ with $d_z = 48$ (the baseline instrument for $k=0$ is dropped to avoid multicollinearity). For the parameters, we set $\beta = 0.5$, $\pi(k) = 0$ if $k = 0$, $\pi(k) = \sqrt{15\sqrt{K}/n}$ for half of the states and $\pi(k) = -\sqrt{15\sqrt{K}/n}$ for the other half. In addition, we set $\gamma(t) = \delta(t) = 1/\sqrt{K}$ for half of the states and $\gamma(t) = \delta(t) = -1/\sqrt{K}$ for the other half.
The error terms are generated as follows: $V_{i,g}$ is generated as
where $\rho = 0.4$, $u_g$ and $v_{i,g}$ are independent standard normal random variables, and $\Phi$ is the CDF of the standard normal distribution so that marginally $V_{i,g}$ is uniformly distributed on $[-1,1]$. Given $V_{i,g}$, we generate
where we set $\mu = 0.4$ and $\sigma_{\varepsilon} = 0.2$. Lastly, given $V_{i,g}$, $\xi_{i,g}$ is generated by
where we use $\mathcal B(a,b,p)$ to denote the binary random variable that equals to $a$ with probability $p$ and $b$ with probability $1-p$, and we set $p = 2/3$. Here $\sigma_{\xi}^{(k)} = 0$ if $k = 0$, and for those states with $\pi(k) = \sqrt{15\sqrt{K}/n}$, we set $\sigma_{\xi}^{(k)} = \sqrt{30\sqrt{K}/n}$ for half of them and $\sigma_{\xi}^{(k)} = -\sqrt{30\sqrt{K}/n}$ for the other half; the same procedure is applied to those states with $\pi(k) = -\sqrt{15\sqrt{K}/n}$.
As shown in EK2018, the TSLS estimand admits a valid (conditional) causal interpretation in this context, since the monotonicity condition is satisfied. In addition, following the same steps as in Yap24, it can be shown that this estimand equals to $\beta$ under the design above. We report the empirical size as well as the power curve below.
We consider a simulation setup similar to the judge design in which the treatment effect of $X$ on $Y$ is heterogeneous and the linear regressions for $X$ and $Y$ are approximately correctly specified due to the binning method. In this setup, individuals within the same cluster are assigned to the same judge. We have in total $G = 200$ clusters, where each cluster contains 4 individuals (so $n = 800$) and the clusters are distributed evenly across $J$ judges.
If cluster $g$ is assigned to judge $j$ for $j = 1, \dots, J$, then the first- and second-stage equations are given by
where $S_{1,i,g} \sim \text{Uniform}[0,1]$ is an individual-level exogenous variable, and $S_{2,g} \in \{0, 1, \ldots, K-1\}$ is a cluster-level discrete exogenous variable with $K = 5$ levels. The error terms are specified as
where $\{\xi_{g}^{(1)}\}_{g \in [G]}$, $\{\xi_{g}^{(2)}\}_{g \in [G]}$, $\{\{\varepsilon_{i,g}^{(1)}\}_{i \in [n_g]}\}_{g \in [G]}$, and $\{\{\varepsilon_{i,g}^{(2)}\}_{i \in [n_g]}\}_{g \in [G]}$ are independent sequences of standard normal random variables. We set $\alpha_0 = -0.8$, $\alpha_1 = 1.0$, $\alpha_2 = 0.6$, $\beta_0 = -0.8$, $\beta_1 = 1$, $\beta_2 = 0.6$, $\beta_3 = 1.0$, $\rho_1 = 0.8$, $\rho_2 = 0.8$ and $\rho_3 = 0.3$ in our simulation.
We generate piecewise constant basis functions by first defining cells based on $(j, S_{2,g})$, then within each cell, we distribute observations evenly into $n_{\text{bins}} = 6$ bins based on $S_{1,i,g}$. Specifically, within each $(j, S_{2,g})$ cell, we sort observations by $S_{1,i,g}$ and assign them to bins such that each bin has approximately the same number of observations. The control variables $\mathcal{W}_{i,g}$ are then constructed as all interactions between these $S_{1,i,g}$ bin indicators and the dummy variables of $S_{2,g}$, i.e., $\{\mathbf{1}\{S_{2,g} = k\}\}_{k=0}^{K-1}$, resulting in up to $6 \times 5 = 30$ control variables.
Next, we consider $J=4$ judges so that all clusters are distributed equally among the four judges. We create the instrumental variables (IVs) by interacting the dummy for the first three judges, $\{\mathbf{1}\{j(i,g)=l\}\}_{l=1,2,3}$, with the control variables, i.e., $Z_{i,g} = \{\mathbf{1}\{j(i,g)=l\} \mathcal{W}_{i,g}\}_{l=1,2,3}$. Therefore, the entire set of regressors $W_{i,g}$ is defined as $W_{i,g} = (Z_{i,g}', \mathcal{W}_{i,g}')'$.
Given that $(\alpha_1 + \alpha_2 S_{1,i,g}) > 0$, the IV monotonicity condition is satisfied. The linear reduced form regressions are approximately correctly specified due to the use of piecewise constant basis functions and full interactions. Consequently, the TSLS estimand admits a valid (conditional) causal interpretation, as established by EK2018. The value of this causal estimand in our simulation is computed via numerical integration. We report the empirical size as well as the power curve below.
Based on the simulation results in Designs \uppercase\expandafter{\romannumeral1}--\uppercase\expandafter{\romannumeral3}, only our proposed inference method controls asymptotic size in settings with many instruments and controls, heterogeneous treatment effects, and clustered data. Existing methods fail for different reasons. TSLS does not correct for the many-instrument bias, and both TSLS and CJIVE fail to address the many-control bias. The three CSW-type methods correct for both many-instrument and many-control biases, but their corrections rely on independence of observations, even within clusters. Moreover, the variance estimators used by TSLS, CJIVE, and the CSW-type methods are not asymptotically valid in our designs for several reasons: (i) they are developed for homogeneous treatment effect models, whereas Designs \uppercase\expandafter{\romannumeral2} and \uppercase\expandafter{\romannumeral3} feature heterogeneous treatment effects; (ii) TSLS and CJIVE do not account for many controls in variance estimation; (iii) CSW-type methods allow a diverging number of controls only at rates slower than \(\sqrt n\); and (iv) CSW-type variance estimators ignore within-cluster dependence. As a result, existing methods exhibit substantial size distortions, while our method performs well. Since our procedure is the only one that controls size, power comparisons are not particularly meaningful. Nevertheless, the power curves indicate that our method retains substantial power. Developing alternative tests that maintain size control while improving on the power of our test represents an interesting direction for future research.
\printbibliography