The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
84,384 characters
Cluster-Robust Inference for Quadratic Forms
\maketitle
\begin{abstract}
This paper studies inference for quadratic forms of linear regression
coefficients with clustered data and many covariates. Our framework covers
three important special cases: instrumental variables regression with many
instruments and controls, inference on variance components, and testing
multiple restrictions in a linear regression. Na\"{\i}ve plug-in estimators
are known to be biased. We study a leave-one-cluster-out estimator that is
unbiased, and provide sufficient conditions for its asymptotic normality. For
inference, we establish the consistency of a leave-three-cluster-out variance
estimator under primitive conditions. In addition, we develop a novel
leave-two-cluster-out variance estimator that is computationally simpler and
guaranteed to be conservative under weaker conditions. Our analysis allows
cluster sizes to diverge with the sample size, accommodates strong
within-cluster dependence, and permits the dimension of the covariates to
diverge with the sample size, potentially at the same rate.\bigskip
\noindent \textbf{Keywords:} Many instruments, many covariates, clustered data, cross-fit, judge design \bigskip
\noindent \textbf{JEL codes:} C12, C36, C55
\end{abstract}
\clearpage
\setcounter{oldtocdepth}{\value{tocdepth}}
\addtocontents{toc}{\protect\setcounter{tocdepth}{-10}}
\section{Introduction}
We study inference for quadratic forms of linear regression coefficients from two regressions with clustered data, given by
\begin{align}\label{eq:model_outcome}
Y_{i, g} &= W_{i, g}' \gamma \;+\; U_{i, g},
& \mathbb{E}[U_{i, g}] &= 0,\\
\label{eq:model_treatment}
X_{i, g} &= W_{i, g}' \pi \;+\; V_{i, g},
& \mathbb{E}[V_{i, g}] &= 0,
\end{align}
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,
\begin{equation}\label{eq:mom}
\theta \;=\; \pi' A_0 \gamma,
\end{equation}
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
\parencite[e.g.,][]{KSS2020}; tests of many linear restrictions can likewise be
cast as a restriction on a quadratic form \parencite{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
\cref{eq:model_outcome,eq:model_treatment} 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
\parencite[see][for recent surveys]{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
\parencite[e.g.,][]{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 \textcite{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 \parencite{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 \parencite{ahpw24}. Likewise, tests for endogenous peer effects can be
cast as testing whether peer characteristics have zero coefficients
\parencite{JL26}. As discussed in \textcite{AS23}, testing for the presence of
heterogeneity amounts to testing whether a set of fixed effects is zero
\parencite[e.g.,][]{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 \parencite{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 \cref{eq:mom} 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 \parencite[e.g.,][]{bjb95,bekker94}. The bias can be purged by using a
\ac{L1CO} estimator proposed by \textcite{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 \emph{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 \textcite{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, \textcite{J22} establishes
the consistency of the \ac{L1CO} unbiased estimator we study, and
\textcite{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 \textcite{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 \textcite{MS24} for a recent survey. Among studies assuming
independent data, our results are most closely related to \textcite{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 \textcite{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 \textcite{LW24} extend it to multidimensional clustering,
adapting bias corrections from \textcite{CNT23} but without formal
distributional theory when there are many
controls.
\section{Setup}\label{sec:setup}
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
\cref{eq:model_outcome,eq:model_treatment}, 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 \cref{eq:mom}, 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$,
\begin{equation}\label{eq:theta_pi}
\hat{\theta}_{\rm PI} = \hat{\pi}' A_0 \hat{\gamma}
=\sum_{g\in[G]} \hat{\pi}' A_0 (W' W)^{-1} W_g' Y_g .
\end{equation}
This plug-in estimator displays an overfitting bias given by
\begin{equation*}
\mathbb{E} \hat{\theta}_{\rm PI}-\theta = \sum_{g \in [G]} \mathbb E V_{g}' A_{g, g} U_{g}, \qquad\text{where}\quad A=W(W' W)^{-1} A_0 (W' W)^{-1} W' \in \Re^{n \times n}.
\end{equation*}
The bias arises because the estimation error in $\hat{\pi}$ is correlated with $Y_{g}$. If we replace $\hat{\pi}$ in \cref{eq:theta_pi} 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
\begin{equation*}
\hat \theta_{\rm LO} =\sum_{g\in[G]} \hat{\pi}_{-g}' A_0 (W' W)^{-1} W_g' Y_g
\end{equation*}
that is unbiased by construction. This estimator was first proposed by \textcite{KSS2020} in the setting where $X=Y$, though their formal results focus on the setting with independent data. As discussed in \textcite{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
\begin{equation}\label{eq:theta_ub}
\hat \theta_{\rm LO} = \hat{\theta}_{\rm PI}
-\sum_{g\in[G]} (M X)_g' M_{g, g}^{-1} A_{g, g} Y_{g}=X' B Y,
\end{equation}
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$. \Cref{eq:theta_ub} 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
\begin{equation*}
\mathbb{E}\sum_{g\in[G]} (MX)_g' M_{g, g}^{-1} A_{g, g} Y_{g}=\mathbb{E}\sum_{g\in[G]} V_{g}' A_{g, g} U_g,
\end{equation*}
\cref{eq:theta_ub} 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, \Cref{lemma:rao} below generalizes the classic results in
\textcite{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$.
\begin{lem}\label{lemma:rao}
The bias-correction matrix $C$ that minimizes $\text{tr}(C'C)$ subject to (i)
$\operatorname{Bdiag}(C_{g,g})=\operatorname{Bdiag}(A_{g,g})$ (ii) $CW=0$, and (iii) $W'C=0$ is given by
$C_{\rm KR}=M\operatorname{Bdiag}(\Lambda)M$, where $\operatorname{bvec}(\Lambda)$ solves
$\operatorname{bvec}(A)=(M\ast M)\operatorname{bvec}(\Lambda)$. If property (ii) is not imposed, the
solution is given by $C_{\rm LO}=M\operatorname{Bdiag}(M_{g,g}^{-1}A_{g,g})$, while if
property (iii) is not imposed, it is given by $C_{\rm LO}'$.
\end{lem}
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~\cref{eq:model_outcome,eq:model_treatment}. 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, \textcite{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
\cref{eq:theta_ub} 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 \textcite{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 \textcite{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
\begin{equation*}
\left \vert \frac{\hat \theta_{\rm LO} - \theta_0}{\hat \omega_n} \right \vert \geq {z}_{1-\alpha/2},
\end{equation*}
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 \Cref{sec:variance_estimation}; asymptotic validity of this test follows from the fact that $\hat{\theta}_{\rm LO}$ is asymptotically normal, as shown in \Cref{sec:asymtptic_normality}. Before giving these results, the remainder of this section fleshes out three important special cases that fit this general setup.
\subsection{Instrumental Variables Regression with Many Instruments and Controls}\label{sec:IV}
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
\begin{equation*}
\mathcal{Y}_{i, g} = W_{i, g}\pi_{Y} + e_{i, g}, \quad \mathbb E [e_{i, g}]=0
\end{equation*}
denote the reduced-form regression and \cref{eq:model_treatment} the first-stage regression, the \ac{TSLS} estimand is given by
\begin{equation}\label{eq:beta}
\beta = \frac{\pi' W' P_{\tilde Z} W \pi_Y}{\pi' W' P_{\tilde Z} W\pi},
\end{equation}
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 \cref{eq:model_outcome} 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
\begin{equation*}
A_0=W' P_{\tilde Z} W.
\end{equation*}
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
\textcite{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
\cref{eq:model_outcome,eq:model_treatment} contains an asymptotically negligible
approximation error, as in, for example, \textcite{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$,
\begin{equation*}
\hat \beta = \frac{X' B \mathcal Y}{X' B X},
\end{equation*}
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
\begin{equation*}
\frac{X' B X (\hat \beta - \beta)}{\hat \omega_n(\hat \beta)}=\frac{X' B (\mathcal Y - \beta X)}{\hat \omega_n(\hat \beta)}
\end{equation*}
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 \textcite{K13} and further advocated by \textcite{GHK25}, which employs a leave-one-out technique to eliminate the overfitting bias of TSLS\@.
\begin{rem}
Validity of the Wald test requires consistency of the variance estimator
$\hat{\omega}_{n}^2(\hat{\beta})$ in the sense that
$\hat \omega_n^2(\hat \beta)/\omega_n^2(\beta) \stackrel{p}{\longrightarrow} 1$, where
$\omega_n^2(\beta) := Var \left(X' B (\mathcal Y - \beta X) \right)$. Our
results below show that $\hat \omega_n^2(\beta) /\omega_n^2(\beta) \stackrel{p}{\longrightarrow} 1$.
Since $\hat{\beta}$ is consistent, a stochastic equicontinuity argument can
then be used to show consistency of the variance estimator
$\hat \omega_n^2(\hat \beta)$, using the fact that
$\hat \omega_n^2(\beta_0)/\omega_n^2(\beta)$ is Lipschitz continuous in
$\beta_0$ over a neighborhood of $\beta$.
\end{rem}
\subsection{Variance and Covariance Components in Linear Regressions}\label{sec:var_comp}
Consider a two-way fixed effect model of log wage determination proposed by
\textcite{AKM99}, in which the log-wage $Y_{\ell, t}$ of a worker $\ell\in[L]$ in year
$t \in [T_{\ell}]$ is given by
\begin{equation}\label{eq:akm}
Y_{\ell, t} =
\alpha_{\ell} + \psi_{\mathcal J(\ell, t)} + \mathcal{W}_{\ell, t}'\delta +
U_{\ell, t}.
\end{equation}
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 \textcite{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 \textcite[Section
4.2]{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 \cref{eq:akm} as
\begin{equation*}
Y_{i, g} = F_{j(g)}'\psi + D_{\ell(g)}'\alpha + \mathcal{W}_{i,g}'\delta + U_{i, g},
\end{equation*}
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
\cref{eq:model_outcome,eq:model_treatment} 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 \cref{eq:mom}, $\sigma_{\psi}^2 = \gamma'A_\psi \gamma$, with
the matrix $A_{\psi}$ given by
\begin{equation*}
A_\psi = \frac{1}{n}
\begin{pmatrix}
(F - 1_{n}\bar{F}')'(F - 1_{n}\bar{F}') & 0 \\
0 & 0
\end{pmatrix},
\end{equation*}
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 \cref{eq:mom} 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
\begin{equation*}
A_{\alpha,\psi} = \frac{1}{2 n }
\begin{pmatrix}
0 & (F - 1_{n}\bar{F}')'(D - 1_{n}\bar{D}') & 0 \\
(D - 1_{n}\bar{D}')'(F - 1_{n}\bar{F}') & 0 & 0 \\
0 & 0 & 0
\end{pmatrix},
\end{equation*}
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 \cref{eq:mom}. 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 \textcite{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
\parencite{ChHe18i}, the importance of geography for healthcare utilization
\parencite{fgw16}, or the importance of classroom assignment in determining
student outcomes \parencite{chetty2011}. In some cases, covariance components
involve fixed effects estimated from different regressions (so that $Y\neq X$):
for instance, \textcite{chetty2011} are interested in the covariance between
classroom fixed effects in an earnings regression and that in a test score
regression.
\subsection{Testing Many Linear Restrictions}\label{sec:testing-many-linear}
Consider the linear regression~\eqref{eq:model_outcome}, so that
\( X_{i,t} = Y_{i,t} \) in \cref{eq:model_treatment}. As in \textcite{AS23}, we
are interested in testing the linear restriction
\begin{equation*}
R \gamma = q,
\end{equation*}
where \( q \) may be high-dimensional. This restriction implies
\begin{equation*}
\gamma' \left[ R' (R(W' W)^{-1} R')^{-1} R \right] \gamma
= q' (R(W' W)^{-1} R')^{-1} q,
\end{equation*}
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
\parencite[e.g.,][]{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,
\textcite{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 \textcite{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).
\section{Asymptotic Normality}\label{sec:asymtptic_normality}
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$.
\begin{ass}\label{ass:dgp}
\begin{enumerate}
\item\label{item:ass_clust} \Cref{eq:model_outcome,eq:model_treatment} hold with $W$ nonrandom and
$\{(U_{g},V_{g})\}_{g \in [G]}$ independent across $g \in [G]$.
\item\label{item:ass_lo} $ \min_{g \in [G]} \lambda_{\min} \left( M_{g,g} \right) \geq c$ for
some constant $c>0$.
\item\label{item:ass_moments}
$\max_{g \in [G]}\max_{i \in [n_{g}]} \left( \mathbb E U_{i,g}^{2q}+\mathbb E
V_{i,g}^{2q} \right) \leq C$ for some constants $C<\infty$ and $q\geq 2$.
In addition, there exists a constant $c>0$ and a sequence $u_{n}<\infty$ such that
\begin{equation*}
c \leq \min_{g \in [G]}\lambda_{\min}\left( \Omega_{g}\right)
\leq \max_{g \in [G]}\lambda_{\max}\left( \Omega_{g}\right) \leq u_n,
\end{equation*}
where
\begin{equation*}
\Omega_{g} = \begin{pmatrix}
\Omega_{U,g} & \Omega_{U,V,g} \\
\Omega_{U,V,g}' & \Omega_{V,g}
\end{pmatrix} = \begin{pmatrix}
\mathbb E U_{g} U_{g}' & \mathbb E U_{g} V_{g}' \\
\mathbb E V_{g} U_{g}' & \mathbb E V_{g} V_{g}'
\end{pmatrix}.
\end{equation*}
\end{enumerate}
\end{ass}
\begin{rem}
\Cref{ass:dgp}.\ref{item:ass_clust} describes the data-generating process for
clustered observations. \Cref{ass:dgp}.\ref{item:ass_lo} ensures that the
model remains estimable after leaving out any particular cluster; otherwise
the estimator $\hat{\theta}_{\textrm{LO}}$ would not be well-defined. When the
observations are independent, the condition is equivalent to
$\min_{i \in [n]} M_{i, i} \geq c$ for some small constant $c>0$, which is
necessary for a leave-observation out regression to be feasible, and standard
in the literature \parencite[e.g.,][]{KSS2020,J22,AS23}. In the context of IV
leniency designs or estimation of variance components, one can ensure this
condition holds by ``pruning'' the sample and dropping observations associated
with singleton decision-makers or firms \parencite[see, e.g.,][for
details]{KSS2020}.
Finally, \Cref{ass:dgp}.\ref{item:ass_moments} imposes mild restrictions on
the moments and the within-cluster covariance structure of the error terms.
The upper bound $u_{n}$ captures the strength of within-cluster dependence.
While we can always take $u_{n}$ to equal a constant times the largest cluster
size $n_{G}$, which is the best bound when the errors share a common
cluster-specific component, a tighter upper bound $u_{n}$ obtains when the
within-cluster errors do not exhibit such strong dependence. For example, in a
panel data setting, where $i$ indexes the time periods for which we observe
individuals indexed by $g$, the errors may satisfy weak dependence,
$cov(U_{i,g},U_{i',g}) = \rho^{|i'-i|}$, in which case $u_n$ is bounded. A
tighter upper bound $u_{n}$ and the existence of higher-order moments (higher
$q$) allow us to weaken the conditions on cluster sizes (see \Cref{rem:rate1}
below).
\end{rem}
\begin{ass}\label{ass:reg}
Let $\Pi = W \pi $, $\Gamma = W \gamma$, $r_n = \operatorname{rank}(A)$, $h_n = \norm{(W'W)^{-1/2} A_0 (W'W)^{-1/2}}_{op}$,
\begin{align*}
H & = B' \Pi/h_n,& \tilde H& = B \Gamma/h_n, & \kappa_n& = \norm{B}_{F}^{2}/h_n^2, \\
\lambda_n& = \max_{g \in [G]} \norm{P_{g,g}}_{op},& \zeta_{H,n}& = \max_{g \in [G]}\norm{H_{g}}_2^2, & \zeta_{\tilde H,n}& = \max_{g \in [G]}
\norm{\tilde H_{g}}_2^2, \\
\eta_{n}& = \max \left\{1, \log^2 (r_n) + \log^2 (n) \lambda_n^2 \right\}&
\mu_n^2 &= \norm{H}_2^2, &\text{and}\quad \tilde \mu_n^2 &= \norm{\tilde H}_2^2,
\end{align*}
\begin{enumerate}
\item\label{item:signal_bound} $\max_{g \in [G]}||\Pi_{g}||_2^2 + \max_{g \in [G]}||\Gamma_{g}||_2^2 \lesssim n_G$.
\item\label{item:clt_rate} $G \to \infty$, and the following condition holds
\begin{multline*}
u_n^3 \eta_n (\mu_n^2 + \tilde \mu_n^2)+u_n^{\frac{3q-4}{q-1}}
n_G^{\frac{q}{q-1}} \lambda_n \kappa_n +u_n^4 \eta_n \kappa_n+u_n^{\frac{2(q-2)}{q-1}} n_G^{\frac{2q}{q-1}} \lambda_n^2 \kappa_n\\
+u_n^{\frac{q-2}{q-1}} n_G^{\frac{q}{q-1}} (\zeta_{H,n}\mu_n^2 + \zeta_{\tilde H,n} \tilde \mu_n^2 + \kappa_n)
= o( (\mu_n^2 + \tilde \mu_n^2 + \kappa_n)^2).
\end{multline*}
\end{enumerate}
\end{ass}
\begin{rem}
A sufficient condition for \Cref{ass:reg}.\ref{item:signal_bound}
to hold is that the signal in \cref{eq:model_outcome,eq:model_treatment},
$\Pi_{i, g}$ and $\Gamma_{i, g}$, is bounded, which is mild.
\end{rem}
\begin{rem}\label{rem:kappa}
Across the three examples, many-IV regression, variance components, and many restrictions, the rank \(r_n\) of \(A\) equals the number of instruments, firms, and restrictions, respectively. In the first and third cases, \(\kappa_n \asymp r_n\). In the variance-components case, \textcite{KSS2020} show \(\kappa_n \ge \tfrac{1}{4}\, J\, \dot{\lambda}_J^{\,2}\) for a stochastic block model for the firm connectivity network,\footnote{Two firms are considered connected if at least one worker moves from one firm to the other.} where $J$ is the number of firms, \(\dot{\lambda}_J\) is the smallest eigenvalue of \(E^{1/2}\mathcal{L}E^{1/2}\), \(\mathcal{L}\) is the normalized graph Laplacian of the employer mobility network, and \(E\) collects employer-specific churn rates. If the network's Cheeger constant is bounded away from zero, then \(\kappa_n \asymp J\).
\end{rem}
\begin{rem}\label{rem:lambda}
The parameter \(\lambda_n\) represents the maximum leverage in the cluster
setting, and we naturally have \(\lambda_n \le 1\). When the design matrix
\(W\) is well balanced, we should have \(\lambda_n \lesssim d/n\). In the
first two examples in \Cref{sec:setup}, \(r_n\) denotes, respectively, the
number of instruments and the number of firms, and \(d\) equals \(r_n\) plus
the number of controls. If the number of controls does not dominate \(r_n\),
then \(d \lesssim r_n\), which, combined with the discussion in
\Cref{rem:kappa}, implies $\lambda_n \lesssim \kappa_n/n$.
\end{rem}
\begin{rem}\label{rem:eta}
To establish the asymptotic distribution of the quadratic form, we bound the
operator norm of the upper-triangular part of \(A\), denoted \(\nabla(A)\).
When \(A\) is a projection matrix, \textcite{Chao12} show that
\(\lVert \nabla(A)\rVert_{\mathrm{op}} \lesssim r_n^{1/4}\), which may be
loose when \(r_n \asymp n\). Instead, we provide a sharper bound,
\begin{equation*}
\norm{\nabla(A)}_{\mathrm{op}} \lesssim \log(r_n)\,\norm{A}_{\mathrm{op}},
\end{equation*}
for a general matrix \(A\), which motivates the definition of \(\eta_n\).
\end{rem}
\begin{rem}\label{rem:linear+quadratic}
In the proof, we show that the estimation error
\(\hat{\theta}_{\mathrm{LO}} - \theta\) is decomposed into linear and
quadratic terms. The variability of the linear component is governed by
\(\mu_n^2 + \tilde{\mu}_n^2\); in the many-IV setting, \(\mu_n^2\) coincides
with the usual concentration parameter in IV regression (with \(Y\) the
outcome and \(X\) the endogenous regressor), while \(\tilde{\mu}_n^2\) is the
concentration parameter for the reverse IV regression of \(X\) on \(Y\). The
variability of the quadratic component is captured by \(\kappa_n\), which is
proportional to \(r_n\) in many cases (as discussed above).
\Cref{ass:reg}.\ref{item:ass_moments} implies that at least one of
\(\mu_n^2 + \tilde{\mu}_n^2\) and \(\kappa_n\) must diverge.
The first condition \((\mu_n^2 + \tilde{\mu}_n^2 \to \infty)\) holds in our
three leading examples when, respectively: (i) the concentration parameter
diverges (many-IV), (ii) the total variation across firms diverges (variance
components), and (iii) the total variation of \(W\) projected onto
\((W'W)^{-1/2}R'\) diverges (many restrictions). If the number of IVs is
fixed, then in the first example, this corresponds to strong identification.
Indeed, it is well known that the usual Lagrange multiplier test under weak
identification is not asymptotically normal \parencite{SS97}. The second
condition \((\kappa_n \to \infty)\) holds when the number of instruments
(many-IV), firms (variance components), or restrictions (many restrictions)
diverges, i.e., \(r_n \to \infty\). This accommodates many-weak-IV settings in
which no consistent test for \(\theta\) exists when
\(\mu_n^2/\sqrt{\kappa_n}\) is bounded \parencite{MS22}. Moreover, because
$\kappa_{n}$ is normalized by the operator norm $h_{n}$, this second condition
also entails the Lindeberg-type condition in \textcite[Section~5]{KSS2020},
yielding a Gaussian limit for \(\hat{\theta}_{\mathrm{LO}}\).
Analyzing the limiting distribution of \(\hat{\theta}_{\mathrm{LO}}\) when the
Lindeberg-type condition fails encounters the challenges highlighted by
\textcite[Section~6]{KSS2020}, especially when \(d \asymp n\) (see
\textcite{LWZ24b} for a valid bootstrap procedure when $d= o(n)$). This issue
is further complicated by within-cluster dependence in our setting, and we
leave a full treatment to future work.
\end{rem}
\begin{rem}\label{rem:rate1}
\Cref{ass:reg}.\ref{item:clt_rate} accommodates several scenarios provided that the regularization parameter satisfies $\lambda_n \lesssim \kappa_n / n$, and $r_n \asymp \kappa_n$, as discussed in \Cref{rem:kappa,rem:lambda}.
\textbf{Scenario 1.} Suppose that \(\kappa_n \asymp n\) and $\zeta_{H,n} + \zeta_{\tilde H,n} \lesssim n_G$. If the within-cluster dependence is weak (i.e., \(u_n \lesssim 1\)), then \Cref{ass:reg}.\ref{item:clt_rate} holds provided \(n_G^{\frac{2q}{q-1}} = o(n)\), regardless of the order of \(\mu_n^2 + \tilde{\mu}_n^2\). If the errors have moments of all orders, a sufficient condition is \(n_G^2 = o(n)\). Under strong within-cluster dependence (i.e., \(u_n \lesssim n_G\)), it suffices that \(n_G^4 \log^2 n = o(n)\), again irrespective of \(\mu_n^2 + \tilde{\mu}_n^2\). The same rate requirements apply under strong identification, in the sense that \(\mu_n^2 + \tilde{\mu}_n^2 \asymp n\), irrespective of the order of \(\kappa_n\).
\textbf{Scenario 2.} Suppose instead that \(\kappa_n \asymp G\), which holds by
construction if IVs are assigned at the cluster level or the number of firms is
less than the number of workers, and $H_g$ and $\tilde H_g$, are balanced across
clusters, i.e.,
\[
\frac{\zeta_{H,n}}{\mu_n^2} \lesssim \frac{1}{G}
\qquad \text{and} \qquad
\frac{\zeta_{\tilde H,n}}{\tilde \mu_n^2} \lesssim \frac{1}{G}.
\]
Then, under weak within-cluster dependence, \Cref{ass:reg}.\ref{item:clt_rate} holds if \(n_G^{\frac{2q}{q-1}} G = o(n^2)\) and \(n_G^{\frac{q}{q-1}} = o(G)\), regardless of the order of \(\mu_n^2 + \tilde{\mu}_n^2\). Under strong within-cluster dependence, a sufficient condition is \(n_G^{4} \log^2 n = o(G)\), again regardless of the order of \(\mu_n^2 + \tilde{\mu}_n^2\).
\textbf{Scenario 3.} If the clusters have a bounded size such that $n_G$ is bounded and $G \asymp n$, then \Cref{ass:reg}.\ref{item:clt_rate} holds as long as \(\mu_n^2 + \tilde{\mu}_n^2 + \kappa_n \to \infty \) and $\zeta_{H,n} + \zeta_{\tilde H,n} = o(\mu_n^2 + \tilde{\mu}_n^2 + \kappa_n)$. This setting includes independent observations as a special case with $n_G = 1$. Under this configuration, our condition is sufficient to guarantee the Lindeberg conditions imposed in Theorems~1 and~2 of \textcite{KSS2020}, which in turn imply the asymptotic normality of the proposed estimator.
\end{rem}
\begin{thm}\label{thm:clt}
Suppose \Cref{ass:dgp,ass:reg} hold. Then
\begin{equation*}
\frac{\hat \theta_{\rm LO} - \theta}{ \omega_n} \rightsquigarrow \mathcal{N}(0,1),
\end{equation*}
where
\begin{multline*}
\omega_n^2 = \sum_{g,h \in [G]^2} tr\left( \Omega_{V,g} B_{g,h} \Omega_{U,h} B_{g,h}'\right) + \sum_{g,h \in [G]^2} tr\left( \Omega_{U,V,g} B_{g,h} \Omega_{U,V,h} B_{h,g} \right)\\
+ \sum_{g \in [G]} H_g' \Omega_{U,g} H_g + \sum_{g \in [G]} \tilde H_g' \Omega_{V,g} \tilde H_g + 2 \sum_{g \in [G]} H_g' \Omega_{U,V,g} \tilde H_g .
\end{multline*}
\end{thm}
\section{Variance Estimator}\label{sec:variance_estimation}
To make use of \Cref{thm:clt} 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.
\subsection{Leave-three-clusters-out Variance Estimator}
The L3CO variance estimator we consider generalizes the leave-three-observations
out estimators considered in \textcite{AS23,Yap24} to the case with clustered
data, and is defined as
\begin{equation*}
\hat \omega^2_{n,\rm L3CO} = \hat \omega^2_{n,\rm L3CO,1} + 2\hat \omega^2_{n,\rm L3CO,2} + \hat \omega^2_{n,\rm L3CO,3} - (\hat \omega^2_{n,\rm L3CO,4} + \hat \omega^2_{n,\rm L3CO,5}),
\end{equation*}
where
\begin{align*}
\hat \omega^2_{n,\rm L3CO,1} & = \sum_{g, h, k \in [G]^3} \left(X_h ' B_{h,g} Y_g \right) \left(X_k ' B_{k,g} \tilde Y_{g,-hk} \right), \\
\hat \omega^2_{n,\rm L3CO,2} & = \sum_{g, h, k \in [G]^3} \left(X_h ' B_{h,g} Y_g \right) \left(Y_k ' B_{g,k}' \tilde X_{g,-hk} \right), \\
\hat \omega^2_{n,\rm L3CO,3} & = \sum_{g, h, k \in [G]^3} \left(Y_h ' B_{g,h}' X_g \right) \left(Y_k ' B_{g,k}' \tilde X_{g,-hk} \right), \\
\hat \omega^2_{n,\rm L3CO,4} & = \sum_{g, h, k \in [G]^3} \left(\tilde{Y}_{h, -gk}' B_{g,h}' X_g \right) \left(Y_h' B_{g,h}' \tilde M_{g,k,-gh} X_k \right), \\
\hat \omega^2_{n,\rm L3CO,5} & = \sum_{g, h, k \in [G]^3}
\left(\tilde{X}_{h, -gk}' B_{h,g} Y_g \right) \left(Y_h' B_{g,h}' \tilde M_{g,k,-gh} X_k \right).
\end{align*}
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,
\begin{align}\label{eq:gamma-ghk}
\hat{\gamma}_{-ghk}&=(W'W-\sum_{l \in (g,h,k)} W_{l}'W_{l})^{-1}(W'Y-\sum_{l
\in (g,h,k)} W_{l}'Y_{l}) \quad\text{and}\\
\hat{\pi}_{-ghk}&=(W'W-\sum_{l \in (g,h,k)} W_{l}'W_{l})^{-1}(W'X-\sum_{l \in (g,h,k)} W_{l}'X_{l})\label{eq:pi-ghk}
\end{align}
denote leave-three-clusters-out OLS estimators, and
\begin{equation*}
\tilde M_{{g,k}, -gh} = \left(M_{g,g} - M_{g,h} M_{h,h}^{-1} M_{h,g} \right)^{-1} \left(M_{g,k} - M_{g,h} M_{h,h}^{-1} M_{h,k} \right).
\end{equation*}
The definition of $\tilde M_{{g,k}, -gh}$ implies that when $g \neq h$,
\begin{equation*}
\tilde M_{g,k,-gh} = \begin{cases}
I_{n_g}, \quad k = g; \\
0_{n_g \times n_h}, \quad k = h; \\
- W_{g} (W'W-\sum_{l \in (g,h)} W_{l}'W_{l})^{-1} W_{k}', \quad k\neq g, k\neq h
\end{cases}
\end{equation*}
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:
\begin{ass}\label{ass:var_L3O}
For $g \neq h \neq k$, define
\begin{align*}
S_{k,g}& = M_{k,k} - M_{k,g} M_{g,g}^{-1} M_{g,k}, \\
\tilde S_{k,gh}& = S_{k,g} - \left(M_{k,h} - M_{k,g} M_{g,g}^{-1} M_{g,h}\right) S_{h,g} ^{-1} \left( M_{h,k} - M_{h,g} M_{g,g}^{-1} M_{g,k} \right),
\end{align*} and let \(\phi_n = \max_{g \in [G]}\sum_{h \in [G]} ||P_{g,h}||_{op}^2\).
\begin{enumerate}
\item\label{item:l3o_well_defined} There exists a finite constant $C>0$ such that
\begin{equation*}
\max_{k,g \in [G]^2, k \neq g} \norm{S_{k,g}^{-1}}_{op} + \max_{g,h,k \in [G]^3, k \neq g, h \neq g, k \neq h}
\norm{\tilde S_{k,gh}^{-1}}_{op} \leq C.
\end{equation*}
\item\label{item:l3o_rates} We have
\begin{multline*}
u_n n_G^3 (\phi_n \lambda_n^2 + \phi_n^2 \lambda_n + \phi_n^3 \lambda_n^2 + \lambda_n^2 )\kappa_n +
u_n^2 n_G^2 \phi_n \lambda_n \kappa_n+
u_n n_G^2 \phi_n \lambda_n^3 (\mu_n^2 + \tilde \mu_n^2)\\
+u_n^{\frac{2q-3}{q-1}} n_G^{\frac{q}{q-1}} (\zeta_{H,n} + \zeta_{\tilde H,n}) \kappa_{n}
+u_n^2 n_G \phi_n (\zeta_{H,n} + \zeta_{\tilde H,n}) \kappa_n
+ u_n^4 n_G \lambda_n^2 \kappa_n = o\left( (\mu_n^2 + \tilde \mu_n^2 + \kappa_n)^2\right).
\end{multline*}
\end{enumerate}
\end{ass}
\begin{rem}
\Cref{ass:var_L3O}.\ref{item:l3o_well_defined} ensures that L3CO least-squares
estimator is well-defined and numerically stable. It is also sufficient, but
not necessary, for the L2CO estimator, defined in \Cref{sec:cons-vari-estim}
below, to be well-defined.
\end{rem}
\begin{rem}\label{rem:phi}
The quantity \(\phi_n\) is naturally bounded by \(\lambda_n n_G\). If cluster
sizes are uniformly bounded, or if
\(\lambda_n \lesssim \kappa_n/n \lesssim G/n\) as discussed in
\Cref{rem:kappa,rem:lambda}, then \(\phi_n \lesssim 1\). In addition, if the
projection matrix \(P\) is sparse in the sense that there are only finitely
many nonzero blocks in the row block matrix \([P_{g,1},\ldots,P_{g,G}]\) for
each \(g \in [G]\)---which occurs when \(W\) partitions the clusters---then
\(\phi_n \lesssim 1\). Finally, if the eigenvalues of \(P_{g,h}P_{h,g}\) are
well balanced for all \(g,h \in [G]^2\) in the sense that
\begin{equation*}
\lVert P_{g,h}P_{h,g}\rVert_{\mathrm{op}} \lesssim \frac{\operatorname{tr}(P_{g,h}P_{h,g})}{n_g},
\end{equation*}
then \(\phi_n \lesssim \lambda_n \le 1\) provided that cluster sizes are well
balanced (i.e., \(n_g \asymp n/G\) for all \(g \in [G]\)). These rate
conditions are directly verifiable, since $\phi_n$ is a known function of the
data.
\end{rem}
\begin{rem}\label{rem:rate2}
We now derive the implications of the rate requirements in
\Cref{ass:var_L3O}.\ref{item:l3o_rates} for the three scenarios in
\Cref{rem:rate1}, assuming $\phi_n \lesssim 1$.
\textbf{Scenario 1.} Suppose that \(\kappa_n \asymp n\) and $\zeta_{H,n} + \zeta_{\tilde H,n} \lesssim n_G$. If the within-cluster dependence is weak (i.e., \(u_n \lesssim 1\)), then \Cref{ass:var_L3O}.\ref{item:l3o_rates} holds provided \(n_G^{3} = o(n)\), regardless of the order of \(\mu_n^2 + \tilde{\mu}_n^2\). Under strong within-cluster dependence (i.e., \(u_n \lesssim n_G\)), it suffices that \(n_G^{5} = o(n)\), again irrespective of \(\mu_n^2 + \tilde{\mu}_n^2\). The same rate requirements apply under strong identification, in the sense that \(\mu_n^2 + \tilde{\mu}_n^2 \asymp n\), regardless of the order of \(\kappa_n\).
\textbf{Scenario 2.} Following Scenario 2 in \Cref{rem:rate1}, if the
within-cluster dependence is weak, then \Cref{ass:var_L3O}.\ref{item:l3o_rates} holds if
\(n_G^{3} = o(n)\) and \(n_G^{\frac{q}{q-1}} = o(G) \), regardless of the order
of \(\mu_n^2 + \tilde{\mu}_n^2\). Under strong within-cluster dependence, a
sufficient condition is \(n_G^{5} G = o(n^2)\), again regardless of
\(\mu_n^2 + \tilde{\mu}_n^2\).
\textbf{Scenario 3.} If the clusters have a bounded size such that \(n_G\) is bounded and \(G \asymp n\), then \Cref{ass:var_L3O}.\ref{item:l3o_rates} holds provided \(\mu_n^2 + \tilde{\mu}_n^2 + \kappa_n \to \infty\) and $\zeta_{H,n} + \zeta_{\tilde H,n} = o(\mu_n^2 + \tilde{\mu}_n^2 + \kappa_n)$.
\end{rem}
The next theorem establishes the unbiasedness and consistency for the L3CO variance estimator.
\begin{thm}\label{thm:var_l3co}
Suppose \Cref{ass:dgp,,ass:reg,ass:var_L3O}.\ref{item:l3o_well_defined} hold. Then
$\mathbb E \hat \omega^2_{n,\rm L3CO} = \omega_n^2$. If, in addition,
\Cref{ass:var_L3O}.\ref{item:l3o_rates} holds, then
\begin{equation*}
\hat \omega_{n,\rm L3CO}^2/\omega_n^2 \stackrel{p}{\longrightarrow} 1.
\end{equation*}
\end{thm}
\begin{rem}
We establish the consistency of \(\hat{\omega}_{n,\mathrm{L3CO}}^{2}\) by showing \(\mathbb{V}(\hat{\omega}_{n,\mathrm{L3CO}}^{2}) = o(\omega_n^{4})\). For this purpose, we need to evaluate sums of the form
\[
\sum_{l \in[G]} P_{g,l,-ghk}\, P_{l,g',-g'h'k'},
\]
where \(P_{g,l,-ghk}\) is the \((g,l)\) block of the leave-\((g,h,k)\)-clusters-out projection matrix for \(W\). Evaluating the summation is not trivial because the sets of clusters omitted in the two projection matrices (i.e., $(g,h,k)$ and $(g',h',k')$) need not coincide. It is further complicated in the clustered setting because \(P_{g,l,-ghk}\) is a block matrix rather than a scalar, and matrix multiplication is non-commutative. There are two key technical innovations in our proof:
(1) a new representation of $P_{g,l,-ghk}$ that enables the exact calculation of the sum, and
(2) a detailed decomposition of the summand into matrices that depend on $(g,g',h,h',k,k')$ only through
\[
s \in \{g, gh, gk, ghk\} \times \{g', g'h', g'k', g'h'k'\},
\]
which implies that they are invariant with respect to the remaining indices $(g,g',h,h',k,k')\setminus s$.\footnote{For example, if a matrix is indexed by $s = (g,h,g',k')$, then it is invariant to the remaining indices $(g,g',h,h',k,k')/s = (k,h')$.}
This invariance allows the summation over the remaining indices to pass through these matrices.
\end{rem}
\begin{rem}
The main difficulty in the variance estimation stems from the leave-out construction, which ensures that the estimators of the linear coefficients $\gamma$ and $\pi$ are independent of the other observations appearing in the same summand. Although such independence could also be achieved by sample splitting, this approach does not directly apply in our setting.
To illustrate, consider the estimation of $\omega^2_{n,1}$. Suppose that the clusters are split into two subsets, $I_1$ and $I_2$, and $(\gamma,\pi)$ is estimated using clusters in $I_1$. A natural sample-splitting estimator is
\begin{equation*}
\hat{\omega}^2_{n,\mathrm{SS},1}
= \sum_{g,h,k \in I_2^3}
\left( X_h' B_{h,g} Y_g \right)
\left( X_k' B_{k,g} \bigl( Y_g - W_g \hat{\gamma}_{I_1} \bigr) \right).
\end{equation*}
However, the summation over $(g,h,k)\in I_2^3$ excludes cross-split terms (e.g., $g,h\in I_2$ but $k\in I_1$). This issue remains even if the roles of $I_1$ and $I_2$ are reversed, showing that simple sample splitting is inadequate here. \Textcite{KSS2020} propose a more complicated sample-splitting variance estimator for independent data. Extending their construction to clustered settings is nontrivial, both theoretically and computationally.
\end{rem}
\begin{rem}
The L3CO variance estimator $\hat \omega_{n,\rm L3CO}^2$ is not necessarily nonnegative in finite samples. However, it is possible to construct a variant of the L3CO estimator that is guaranteed to be nonnegative and consistent. Specifically, let
\begin{align*}
& \tilde \omega^2_{n,\rm L3CO,1} = \hat \omega^2_{n,\rm L3CO,1} + 2\hat \omega^2_{n,\rm L3CO,2} + \hat \omega^2_{n,\rm L3CO,3} - 2(\hat \omega^2_{n,\rm L3CO,4} + \hat \omega^2_{n,\rm L3CO,5}),\\
& \tilde \omega^2_{n,\rm L3CO,2} = \hat \omega^2_{n,\rm L3CO,4} + \hat \omega^2_{n,\rm L3CO,5}.
\end{align*}
We note that \(\mathbb{E}\tilde{\omega}^2_{n,\mathrm{L3CO},1}\) corresponds to the variance of the linear component of \(X'BY\), while \(\mathbb{E}\tilde{\omega}^2_{n,\mathrm{L3CO},2}\) corresponds to the variance of its quadratic component. Since
\[
\hat{\omega}^2_{n,\mathrm{L3CO}}
=
\tilde{\omega}^2_{n,\mathrm{L3CO},1}
+
\tilde{\omega}^2_{n,\mathrm{L3CO},2},
\]
and both \(\tilde{\omega}^2_{n,\mathrm{L3CO},1}\) and \(\tilde{\omega}^2_{n,\mathrm{L3CO},2}\) are asymptotically nonnegative, this motivates the following variant of the L3CO estimator:
\[
\tilde{\omega}^2_{n,\mathrm{L3CO}}
=
\bigl|\tilde{\omega}^2_{n,\mathrm{L3CO},1}\bigr|
+
\bigl|\tilde{\omega}^2_{n,\mathrm{L3CO},2}\bigr|.
\]
\end{rem}
\begin{cor}\label{cor:L3CO}
Suppose \Cref{ass:dgp,,ass:reg,ass:var_L3O} hold. Then
\begin{equation*}
\tilde \omega^2_{n,\rm L3CO}/\omega_n^2 \stackrel{p}{\longrightarrow} 1.
\end{equation*}
\end{cor}
\subsection{Leave-two-clusters-out Variance Estimator}\label{sec:cons-vari-estim}
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
\begin{equation*}
\hat \omega^2_{n, \rm L2CO} = \sum_{g \in [G]} \left( \sum_{h \in [G]} \left( X_h ' B_{h,g} \tilde Y_{g,-h} + Y_h ' B_{g,h}' \tilde X_{g,-h} \right) \right)^2,
\end{equation*}
where
\begin{align*}
\tilde Y_{g,-h}& = Y_g - W_g \hat \gamma_{-gh},&
\tilde X_{g,-h}& = X_g - W_g \hat \pi_{-gh}
\end{align*}
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
\cref{eq:gamma-ghk,eq:pi-ghk}). 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 \Cref{thm:var_l2co} 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.
\begin{ass}\label{ass:var_l2co}
For $S_{k,g}$ defined in \Cref{ass:var_L3O}.\ref{item:l3o_well_defined}, the
following holds:
\begin{enumerate}
\item\label{item:l2o_well_defined} There exists a finite constant $C > 0$
such that
$\max_{k,g \in [G]^2, k \neq g} \norm{S_{k,g}^{-1}}_{op}
\leq C$.
\item\label{item:l2o_rates} We have
\begin{multline*}
u_n^{\frac{2q-3}{q-1}} n_G^{\frac{q}{q-1}}(\zeta_{H,n} +
\zeta_{\tilde{H}, n}) \kappa_n +u_n^2 n_G \phi_n \lambda_n^3
(\zeta_{H,n} + \zeta_{\tilde{h}, n}) \kappa_n +
u_n^{\frac{2q-3}{q-1}} n_G^{\frac{2q-1}{q-1}} \phi_n^2 \lambda_n^2 \kappa_n\\
+u_n^3 n_G(\phi_n^2 \lambda_n + \phi_n \lambda_n^2) \kappa_n+
u_n^4 n_G \lambda_n^2 \kappa_n
= o\left( (\mu_n^2 + \tilde \mu_n^2 + \kappa_n)^2\right),
\end{multline*}
where $\phi_n$ is defined in \Cref{ass:var_L3O}.
\end{enumerate}
\end{ass}
\begin{thm}\label{thm:var_l2co}
Suppose \Cref{ass:dgp,,ass:reg,ass:var_l2co}.\ref{item:l2o_well_defined} hold. Then
$\mathbb E \hat \omega^2_{n, \rm L2CO} \geq \omega_n^2$. If, in addition,
\Cref{ass:var_l2co}.\ref{item:l2o_rates} holds, then
\begin{equation}\label{eq:l2o_conservative}
\limsup_{n \to \infty } \mathbb P \left( \left \vert
\frac{\hat \theta_{\rm LO} - \theta}{\hat \omega_{n,\rm L2CO}} \right \vert \geq {z}_{1-\alpha/2} \right) \leq \alpha.
\end{equation}
Suppose further that there exists a finite constant $C > 0$ such that
$\max_{k,g \in [G]^2, k \neq g} \norm{S_{k,g}^{-1}}_{op} \leq
C$, and that $n_G = O(1)$, $\lambda_n = o(1)$ and
$\kappa_n = o (\mu_n^2 + \tilde \mu_n^2)$. Then~\cref{eq:l2o_conservative}
holds with equality.
\end{thm}
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.
\begin{rem}\label{rem:rate3}
One can easily check that the rate conditions in \Cref{ass:var_l2co}.\ref{item:l2o_rates} are weaker than those in \Cref{ass:var_L3O}.\ref{item:l3o_rates}. In particular, consider again the three scenarios in \Cref{rem:rate2}.
\textbf{Scenario 1.} Suppose that \(\kappa_n \asymp n\) and
$\zeta_{H,n} + \zeta_{\tilde H,n} \lesssim n_G$. If the within-cluster
dependence is weak, then \Cref{ass:var_l2co}.\ref{item:l2o_rates} holds provided
\(n_G^{\frac{2q-1}{q-1}} = o(n)\), which is weaker than \(n_G^{3} = o(n)\)
if $q > 2$. Under strong within-cluster dependence, for \Cref{ass:var_l2co}.\ref{item:l2o_rates}
to hold, we still require \(n_G^{5} = o(n)\).
\textbf{Scenario 2.} Following the setting of Scenario~2 in
Remark~\ref{rem:rate1}, Assumption~\ref{ass:var_l2co}.\ref{item:l2o_rates} holds under weak
within-cluster dependence if \( n_G^{\frac{q}{q-1}} = o(G)\) and
\(n_G^{\frac{2q-1}{q-1}} G = o(n^2)\), which is weaker than the requirement
\(n_G^3 = o(n)\) imposed by \Cref{ass:var_L3O}.\ref{item:l3o_rates} in this scenario,
provided that $q>2$. Under strong within-cluster dependence, both
\Cref{ass:var_L3O}.\ref{item:l3o_rates} and \Cref{ass:var_l2co}.\ref{item:l2o_rates} require
$n_G^{5} G = o(n^2)$.
\textbf{Scenario 3.} If clusters have bounded size so that \(n_G\) is
bounded and \(G \asymp n\), then both
Assumption~\ref{ass:var_L3O}.\ref{item:l3o_rates} and
Assumption~\ref{ass:var_l2co}.\ref{item:l2o_rates} hold provided
\(\mu_n^2 + \tilde{\mu}_n^2 + \kappa_n \to \infty\) and
$\zeta_{H,n} + \zeta_{\tilde H,n} = o(\mu_n^2 + \tilde{\mu}_n^2 +
\kappa_n)$.
\end{rem}
\begin{rem}
Note that \(\hat{\omega}^2_{n,\mathrm{L2CO}}\) remains computable even when
leaving three clusters out is infeasible, in the sense that some of the
matrices \(\tilde S_{k,gh}\) in \Cref{ass:var_L3O}.\ref{item:l3o_well_defined}
are not invertible. One may also construct conservative variance estimators
based on ``HC3''-type residuals \parencite[see, e.g.,][]{CJN18}, which can be
computed even when leaving two clusters out is infeasible, that is, when some
of the matrices \(S_{k,g}\) in \Cref{ass:var_L3O}.\ref{item:l3o_well_defined}
are not invertible. Indeed, the computation of HC3-type residuals only
requires \Cref{ass:dgp}.\ref{item:ass_lo}. Establishing formal theoretical
guarantees for such estimators is left for future research.
\end{rem}
\subsection{Practical Guidance}
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 \Cref{sec:setup} implies that the
leave-out residuals may be computed directly as
\begin{align*}
\tilde Y_{g,-hk}& = \left[M_{ghk,ghk}^{-1} (MY)_{ghk} \right]_g,& \text{and}&&
\tilde X_{g,-hk}& = \left[M_{ghk,ghk}^{-1} (MX)_{ghk} \right]_g,
\end{align*}
where $M_{ghk,ghk}$ is the block diagonal matrix of $M$ corresponding to
clusters $(g,h,k)$. Similarly,
\begin{align*}
\tilde Y_{g,-h}& = \left[M_{gh,gh}^{-1} (MY)_{gh} \right]_g & \text{and}&& \tilde X_{g,-h}& = \left[M_{gh,gh}^{-1} (MX)_{gh} \right]_g.
\end{align*}
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.
\section{Simulation}
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 \textcite{FLM23} with cluster-robust variance estimator; and the three test statistics (CSW-JIV, CSW-LIM and CSW-FUL) proposed in \textcite{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.
\subsection{Design \texorpdfstring{\uppercase\expandafter{\romannumeral1}}{I}: Homogeneous Treatment Effect}
We first consider the following panel IV regression adapted from \textcite{CNT23}, which assumes a homogeneous treatment effect:
\begin{align*}
\mathcal Y_{i,g} &= X_{i,g}\beta + \mathcal W_{i,g}'\gamma + \alpha_g + U_{i,g}, \\
X_{i,g} &= Z_{i,g}'\pi + \mathcal W_{i,g}'\delta + \xi_g + V_{i,g},
\end{align*}
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
\begin{equation*}
\Omega_{1g}= \begin{bmatrix}
1 & \theta_1 & \cdots & \theta_1 \\
\theta_1 & 1 & \cdots & \theta_1 \\
\vdots & \vdots & \ddots & \vdots \\
\theta_1 & \theta_1 & \cdots & 1
\end{bmatrix}_{n_g \times n_g} \;\; g=1, \dots, G
\end{equation*}
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
\begin{align*}
\tilde{U}_{i,g} &= \rho \varepsilon_{i,g} + \sqrt{1-\rho^2} \sigma_{i,g} v_{i,g}, \\
\tilde{V}_{i,g} &= \rho \eta_{i,g} + \sqrt{1-\rho^2} \sigma_{i,g} v_{i,g},
\end{align*}
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
\begin{equation*}
\Omega_{2g}= \begin{bmatrix}
1 & 0 & \cdots & 0 \\
\theta_2 & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
\theta_2^{n_g-1} & \theta_2^{n_g-2} & \cdots & 1
\end{bmatrix}_{n_g \times n_g} \;\; g=1, \dots, G,
\end{equation*}
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.
\begin{table}[tp]
\centering
\begin{tabular}{ccccccc}
& \multicolumn{3}{c}{5\% significance level} & \multicolumn{3}{c}{10\% significance level} \\
\cmidrule(rl){2-4} \cmidrule(rl){5-7}
& $d_z=d_w=50$ & $100$ & $150$ & $d_z=d_w=50$ & $100$ & $150$ \\
\midrule
\text{TSLS} & 53.7\% & 92.3\% & 99.7\% & 66.5\% & 96.8\% & 99.9\% \\
\text{CJIVE} & 7.9\% & 22.3\% & 53.4\% & 15.5\% & 31.4\% & 63.3\% \\
\text{CSW-JIV} & 7.8\% & 13.3\% & 22.7\% & 14.2\% & 20.8\% & 31.0\% \\
\text{CSW-LIM} & 7.1\% & 10.5\% & 16.6\% & 14.6\% & 16.8\% & 25.8\% \\
\text{CSW-FUL} & 7.6\% & 10.5\% & 16.4\% & 14.7\% & 16.9\% & 25.6\% \\
\text{L3CO} & 6.3\% & 5.0\% & 5.0\% & 10.4\% & 9.1\% & 9.4\% \\
\bottomrule
\end{tabular}
\caption{Empirical size at 5\% and 10\% significance levels for the \textcite{CNT23} simulation.}\label{tab:size_CSW_combined}
\end{table}
\begin{figure}[tp]
\centering
\begin{tabular}{ccc}
\multicolumn{3}{c}{\textbf{5\% significance level}} \\[2mm]
\includegraphics[width=0.31\textwidth]{res_tstat_50_50_5pct.png} &
\includegraphics[width=0.31\textwidth]{res_tstat_100_100_5pct.png} &
\includegraphics[width=0.31\textwidth]{res_tstat_150_150_5pct.png} \\
$d_z=50,\ d_w=50$ & $d_z=100,\ d_w=100$ & $d_z=150,\ d_w=150$ \\[4mm]
\multicolumn{3}{c}{\textbf{10\% significance level}} \\[2mm]
\includegraphics[width=0.31\textwidth]{res_tstat_50_50_10pct.png} &
\includegraphics[width=0.31\textwidth]{res_tstat_100_100_10pct.png} &
\includegraphics[width=0.31\textwidth]{res_tstat_150_150_10pct.png} \\
$d_z=50,\ d_w=50$ & $d_z=100,\ d_w=100$ & $d_z=150,\ d_w=150$
\end{tabular}
\caption{Power curves for the \textcite{CNT23} simulation. Top row: 5\% significance level. Bottom row: 10\% significance level.}\label{fig:power_CSW_combined}
\end{figure}
\subsection{Design \texorpdfstring{\uppercase\expandafter{\romannumeral2}}{II}: Heterogeneous Treatment Effect, Saturated}
We consider a simulation setup similar to that in \textcite{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
\begin{align*}
\mathcal Y_{i,g} &= X_{i,g} (\beta + \xi_{i,g}) + \mathcal W_{i,g}'\gamma + U_{i,g} \\
X_{i,g} &= \mathbf{1} \left\{ Z_{i,g}' \pi + W_{i,g}'\delta \geq V_{i,g} \right\}
\end{align*}
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
\begin{align*}
V_{i,g} = 2 \Phi (\rho u_g + \sqrt{1-\rho^2} v_{i,g}) - 1
\end{align*}
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
\begin{align*}
U_{i,g} | V_{i,g} \sim
\begin{cases}
\mathcal N (\mu, \sigma_{\varepsilon}^2) & V_{i,g} \geq 0 \\
\mathcal N (-\mu, \sigma_{\varepsilon}^2) & V_{i,g} < 0
\end{cases}
\end{align*}
where we set $\mu = 0.4$ and $\sigma_{\varepsilon} = 0.2$. Lastly, given $V_{i,g}$, $\xi_{i,g}$ is generated by
\begin{align*}
\xi_{i,g} | V_{i,g} \sim
\begin{cases}
\mathcal B(\sigma_{\xi}^{(k)}, -\sigma_{\xi}^{(k)}, p) & V_{i,g} \geq 0 \\
\mathcal B(\sigma_{\xi}^{(k)}, -\sigma_{\xi}^{(k)}, 1-p) & V_{i,g} < 0
\end{cases}
\end{align*}
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 \textcite{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 \textcite{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.
\begin{table}[tp]
\centering
\begin{tabular}{ccc}
& 5\% & 10\% \\ \midrule
\text{TSLS} & 59.6\% & 71.7\% \\
\text{CJIVE} & 10.0\% & 15.8\% \\
\text{CSW-JIV} & 13.8\% & 21.5\% \\
\text{CSW-LIM} & 33.0\% & 45.8\% \\
\text{CSW-FUL} & 30.4\% & 43.0\% \\
\text{L3CO} & 4.8\% & 9.4\% \\ \bottomrule
\end{tabular}
\caption{Empirical size at 5\% and 10\% significance levels for the \textcite{Yap24} simulation.}\label{tab:size_LY_combined}
\end{table}
\begin{figure}[tp]
\centering
\begin{tabular}{cc}
\includegraphics[width=0.47\textwidth]{res_LY_5pct.png} &
\includegraphics[width=0.47\textwidth]{res_LY_10pct.png} \\
5\% significance level & 10\% significance level
\end{tabular}
\caption{Power curves for the \textcite{Yap24} simulation.}\label{fig:power_LY_combined}
\end{figure}
\subsection{Design \texorpdfstring{\uppercase\expandafter{\romannumeral3}}{III}: Heterogeneous Treatment Effect, Approximated}
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
\begin{align*}
X_{i,g} &= \mathbf{1} \left\{U_{i,g} \leq \alpha_0 + \left(\alpha_1 + \alpha_2 S_{1,i,g} \right) \times \frac{j}{J} \right\}, \\
\mathcal Y_{i,g} &= \mathbf{1} \left\{V_{i,g} \leq \beta_0 + \beta_1 X_{i,g} + \beta_2 S_{1,i,g} + \beta_3 S_{2,g} \right\},
\end{align*}
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
\begin{align*}
U_{i,g} &= \rho_1 \xi_{g}^{(1)} + \sqrt{1 - \rho_1^2} \, \varepsilon_{i,g}^{(1)}, \\
V_{i,g} &= \rho_2 U_{i,g} + \sqrt{1 - \rho_2^2} \left( \rho_3 \xi_{g}^{(2)} + \sqrt{1 - \rho_3^2} \, \varepsilon_{i,g}^{(2)} \right),
\end{align*}
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 \textcite{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.
\begin{table}[tp]
\centering
\begin{tabular}{ccc}
& 5\% & 10\% \\ \midrule
\text{TSLS} & 47.5\% & 60.7\% \\
\text{CJIVE} & 11.3\% & 19.7\% \\
\text{CSW-JIV} & 21.8\% & 31.2\% \\
\text{CSW-LIM} & 16.5\% & 23.3\% \\
\text{CSW-FUL} & 16.7\% & 23.4\% \\
\text{L3CO} & 4.1\% & 8.0\% \\
\bottomrule
\end{tabular}
\caption{Empirical size at 5\% and 10\% significance levels for the judge design simulation.}\label{tab:size_judge_combined}
\end{table}
\begin{figure}[tp]
\centering
\begin{tabular}{cc}
\includegraphics[width=0.48\textwidth]{res_judge_const_5pct.png} &
\includegraphics[width=0.48\textwidth]{res_judge_const_10pct.png} \\
5\% significance level & 10\% significance level
\end{tabular}
\caption{Power curves for the judge design simulation.}\label{fig:power_judge_combined}
\end{figure}
\subsection{Remarks}
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
\newpage