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.
140,002 characters
Leave-out estimation of variance components
\maketitle
\begin{abstract}
We propose leave-out estimators of quadratic forms designed for the study of linear models with unrestricted heteroscedasticity. Applications include analysis of variance and tests of linear restrictions in models with many regressors.
An approximation algorithm is provided that enables accurate computation of the estimator in very large datasets.
We study the large sample properties of our estimator allowing the number of regressors to grow in proportion to the number of observations. Consistency is established in a variety of settings where plug-in methods and estimators predicated on homoscedasticity exhibit first-order biases.
For quadratic forms of increasing rank, the limiting distribution can be represented by a linear combination of normal and non-central $\chi^2$ random variables, with normality ensuing under strong identification. Standard error estimators are proposed that enable tests of linear restrictions and the construction of uniformly valid confidence intervals for quadratic forms of interest.
We find in Italian social security records that leave-out estimates of a variance decomposition in a two-way fixed effects model of wage determination yield substantially different conclusions regarding the relative contribution of workers, firms, and worker-firm sorting to wage inequality than conventional methods. Monte Carlo exercises corroborate the accuracy of our asymptotic approximations, with clear evidence of non-normality emerging when worker mobility between blocks of firms is limited.
\end{abstract}
\small{Keywords: variance components, heteroscedasticity, fixed effects, leave-out estimation, many regressors, weak identification, random projection}
\newpage{}
\onehalfspacing
As economic datasets have grown large, so has the number of parameters employed in econometric models. Typically, researchers are interested in certain low dimensional summaries of these parameters that communicate the relative influence of the various economic phenomena under study. An important benchmark comes from \citet{fisher1925statistical}'s foundational work on analysis of variance (ANOVA) which he proposed as a means of achieving a ``separation of the variance ascribable to one group of causes, from the variance ascribable to other groups.''\footnote{See \citet{cochran1980fisher} for a discussion of the intellectual development of this early work.}
A large experimental literature \citep{sacerdote2001peer,graham2008identifying,chetty2011does,angrist2014perils} employs variants of Fisher's ANOVA approach to infer the degree of variability attributable to peer or classroom effects. Related methods are often used to study heterogeneity across firms, workers, and schools in their responsiveness to exogenous regressors with continuous variation \citep{raudenbush1986hierarchical,raudenbush2002hierarchical,arellano2011identifying,graham2012identification}.
In labor economics, log-additive models of worker and firm fixed effects are increasingly used to study worker-firm sorting and the dispersion of firm specific pay premia \citep{abowd1999high,card2013workplace,card2016firms,song2015firming,sorkin2018ranking} and analogous methods have been applied to settings in health economics \citep{finkelstein2016sources,silver2016essays} and the economics of education \citep{arcidiacono2012estimating}.
This paper considers estimation of and inference on \emph{variance components}, which we define broadly as quadratic forms in the parameters of a linear model. Notably, this definition yields an important connection to the recent literature on testing linear restrictions in models with many regressors \citep{anatolyev2012inference,chao2014testing,cattaneo2017inference}. Traditional variance component estimators are predicated on the assumption that the errors in a linear model are identically distributed draws from a normal distribution. Standard references on this subject \citep[e.g.,][]{searle2009variance} suggest diagnostics for heteroscedasticity and non-normality, but offer little guidance regarding estimation and inference when these problems are encountered. A closely related literature on panel data econometrics proposes variance component estimators designed for fixed effects models that either restrict the dimensionality of the underlying group means \citep{bonhomme2019distributional} or the nature of the heteroscedasticity governing the errors \citep{andrews2008high,jochmans2016fixed}.
Our first contribution is to propose a new variance component estimator designed for unrestricted linear models with
heteroscedasticity of unknown form. The estimator is finite sample unbiased and can be written as a naive ``plug-in'' variance component estimator plus a bias correction term that involves ``cross-fit'' \citep{newey2018cross} estimators of observation-specific error variances. We also develop a representation of the estimator in terms of a covariance between outcomes and a ``leave-one-out'' generalized prediction \cite[e.g., as in][]{powell1989semiparametric}, which allows us to apply recent results on the behavior of second order U-statistics. Building on work by \cite{achlioptas2001database}, we propose a random projection method that enables computation of our estimator in very large
datasets with little loss of accuracy.
We study the asymptotic behavior of the proposed leave-out estimator in an environment where the number of regressors may be proportional to the sample size: a framework that has alternately been termed ``many covariates'' \citep{cattaneo2017inference} or ``moderate dimensional'' \citep{lei2016asymptotics} asymptotics. Verifiable design requirements are provided under which the estimator is consistent and we show in an Appendix that these conditions are weaker than those required by jackknife bias correction procedures \citep{quenouille1949approximate,hahn2004jackknife,dhaene2015split}.
A series of examples is discussed where the leave-out estimator is consistent, while estimators relying on jackknife or homoscedasticity-based bias corrections are not.
We present three sets of theoretical results that enable inference based upon our estimator in a variety of settings. The first result concerns inference on quadratic forms of fixed rank, a problem which typically arises when testing a few linear restrictions in a model with many covariates \citep{cattaneo2017inference}. Familiar examples of such applications include testing that particular regressors are significant in a fixed effects model and conducting inference on the coefficients from a projection of fixed effects onto a low dimensional vector of covariates. Extending classic proposals by \cite{horn1975estimating} and \cite{mackinnon1985some}, we show that our leave-out approach can be used to construct an Eicker-White style variance estimator that is unbiased in the presence of unrestricted heteroscedasticity and that enables consistent inference on linear contrasts under weaker design restrictions than those considered by \citet{cattaneo2017inference}.
Next, we derive a result establishing asymptotic normality of quadratic forms of growing rank. Such quadratic forms typically arise when conducting analysis of variance but also feature in tests of model specification involving a large number of linear restrictions \citep{anatolyev2012inference,chao2014testing}. The large sample distribution of the estimator is derived using a variant of the arguments in \citet{chatterjee2008new} and \citet{soelvsten2017robust} and a standard error estimator is proposed that utilizes sample splitting formulations of the sort considered by \cite{newey2018cross}. This standard error estimator is shown to enable consistent inference on quadratic forms of growing rank in the presence of unrestricted heteroscedasticity when the regressor design allows for sample splitting and to provide conservative inference otherwise.
Finally, we present conditions under which the large sample distribution of our estimator is non-pivotal and can be represented by a linear combination of normal and non-central $\chi^{2}$ random variables, with the non-centralities of the $\chi^{2}$ terms serving as weakly identified nuisance parameters. This distribution arises in a two-way fixed effects model when there are ``bottlenecks'' in the mobility network.
Such bottlenecks are shown to emerge, for example, when worker mobility is governed by a stochastic block model with limited mobility between blocks. To construct asymptotically valid confidence intervals in the presence of nuisance parameters, we propose inverting a minimum distance test statistic. Critical values are obtained via an application of the procedure of \cite{andrews2016geometric}. The resulting confidence interval is shown to be valid uniformly in the values of the nuisance parameters and to have a closed form representation in many settings, which greatly simplifies its computation.
We illustrate our results with an application of the two-way worker-firm
fixed effects model of \citet{abowd1999high} to Italian social security records.
The proposed leave-out estimator
finds a substantially smaller contribution of firms to wage inequality
and much more assortativity in the matching of workers to firms than
either the uncorrected plug-in estimator originally considered by \citet{abowd1999high}
or the homoscedasticity-based correction procedure of \citet{andrews2008high}.
When studying panels of length greater than two, we allow for
serial correlation in the errors by employing a generalization of our estimator
that leaves out all the observations in a worker-firm match. Failing to account for this dependence
is shown to yield over-estimates of the variance of firm effects.
Projecting firm effect estimates onto measures of worker age and firm size, we find
that older workers tend to be employed at firms offering higher firm wage effects;
however, this phenomenon is largely explained by the tendency of older
workers to sort to bigger firms. Leave-out standard errors for the coefficients
of these linear projections are found to be several times larger than
a naive standard error predicated on the assumption that the estimated fixed effects are independent of each other.
Stratifying our analysis by birth cohort, we formally reject the null hypothesis
that older and younger workers face identical vectors of firm effects. However, the two sets of firm effects are
estimated to have a correlation coefficient of nearly 0.9, while the plug-in estimate of correlation is only 0.54.
To assess the accuracy of our asymptotic
approximations, we conduct a series of Monte Carlo exercises utilizing the realized mobility patterns of workers between firms. Clear evidence of non-normality arises in the sampling distribution
of the estimated variance of firm effects
in settings where the worker-firm mobility network is weakly connected. The proposed confidence regions are shown to provide reliable size control in both strongly and weakly identified settings.
\section{Unbiased Estimation of Variance Components}
\label{sec:unbiased}
Consider the linear model
\begin{align}\label{eq:model}
y_{i} = x_{i}'\beta + \varepsilon_{i} && (i = 1,\ldots,n)
\end{align}
where the regressors $x_{i}\in\mathbb{R}^{k}$ are non-random and the design matrix $S_{xx}=\sum_{i=1}^{n}x_{i}x_{i}'$ has full rank. The unobserved errors $\{ \varepsilon_{i} \}_{i=1}^{n}$ are mutually independent and obey $\mathbb{E}[\varepsilon_{i}]=0$, but may possess observation specific variances $\mathbb{E}[\varepsilon_{i}^2]=\sigma_{i}^{2}$.
Our object of interest is a quadratic form $\theta=\beta'A\beta$ for some known non-random symmetric matrix $A\in\mathbb{R}^{k\times k}$ of rank $r$. Following \cite{searle2009variance}, when $A$ is positive semi-definite $\theta$ is a \emph{variance} component, while when $A$ is non-definite $\theta$ may be referred to as a \emph{covariance} component. Note that linear restrictions on the parameter vector $\beta$ can be formulated in terms of variance components: for a non-random vector $v$, the null hypothesis $v'\beta=0$ is equivalent to the restriction $\theta=0$ when $A=vv'$. Examples from the economics literature where variance components are of direct interest are discussed in Section \ref{sec:examples}.
\subsection{Estimator}
A naive plug-in estimator of $\theta$ is given by the quadratic form $\hat \theta_{\text{PI}}= \hat{\beta}'A\hat{\beta}$, where $\hat{\beta}=S_{xx}^{-1}\sum_{i=1}^{n}x_{i}y_{i}$ denotes the Ordinary Least Squares (OLS) estimator of $\beta$. Estimation error in $\hat \beta$ leads the plug-in estimator to exhibit a bias involving a linear combination of the unknown variances $\{ \sigma_{i}^{2}\} _{i=1}^{n}$. Specifically, standard results on quadratic forms imply that $\mathbb{E} [ \hat \theta ]=\theta + \text{trace}( A \mathbb{V}[\hat \beta])$, where
\begin{align}\label{eq:bias}
\text{trace}\left( A \mathbb{V}[\hat \beta]\right) = \sum_{i=1}^{n} B_{ii} \sigma_{i}^{2} & \quad \mbox{and} \quad B_{ii}=x_{i}' S_{xx}^{-1} A S_{xx}^{-1} x_{i}.
\end{align}
As discussed in Section \ref{sec:examples}, this bias can be particularly severe when the dimension of the regressors $k$ is large relative to the sample size.
A bias correction can be motivated by observing that an unbiased estimator of the $i$-th error variance is
\begin{align}\label{eq:sigmahat}
\hat{\sigma}_{i}^{2}=y_{i}\left(y_{i}-x_{i}'\hat{\beta}_{-i}\right),
\end{align}
where $\hat \beta_{-i}=\left(S_{xx} - x_i x_i' \right)^{-1}\sum_{\ell\neq i}x_{\ell}y_{\ell}$ denotes the leave-$i$-out OLS estimator of $\beta$. This insight suggests the following bias-corrected estimator of $\theta$:
\begin{align}\label{eq:estimator}
\hat{\theta}=\hat{\beta}'A\hat{\beta}-\sum_{i=1}^{n} B_{ii}\hat{\sigma}_{i}^{2}.
\end{align}
While \cite{newey2018cross} observe that ``cross-fit'' covariances relying on sample splitting can be used to remove bias of the sort considered here, we are not aware of existing estimators involving the leave-one-out estimators $\{\hat{\sigma}_{i}^{2}\} _{i=1}^{n}$.
One can also motivate $\hat{\theta}$ via a change of variables argument. Letting $\tilde{x}_{i}=A S_{xx}^{-1} x_{i}$ denote a vector of ``generalized'' regressors, we can write
\begin{align}
\theta = \beta' A \beta = \beta' S_{xx} S_{xx}^{-1} A \beta = \sum_{i=1}^{n}\beta' x_{i} \tilde {x}_{i}'\beta = \sum_{i=1}^{n}\mathbb{E}\left[y_{i}\tilde{x}_{i}'\beta\right].
\end{align}
This observation suggests using the unbiased \emph{leave-out} estimator
\begin{align}\label{eq:cov}
\hat \theta = \sum_{i=1}^{n}y_{i}\tilde{x}_{i}'\hat{\beta}_{-i}.
\end{align}
Note that direct computation of $\hat \beta_{-i}$ can be avoided by exploiting the representation
\begin{align}\label{eq:Pii}
y_{i}-x_{i}'\hat{\beta}_{-i}=\dfrac{y_i - x_i'\hat \beta}{1-P_{ii}},
\end{align}
where $P_{ii}= x_i' S_{xx}^{-1} x_i$ gives the leverage of observation $i$. Applying the Sherman-Morrison-Woodbury formula \citep{woodbury1949stability,sherman1950adjustment}, this representation also reveals that \eqref{eq:estimator} and \eqref{eq:cov} are numerically equivalent:
\begin{align}
y_i \tilde x_i'\hat \beta_{-i} &= \underbrace{y_i \tilde x_i' S_{xx}^{-1} \sum_{\ell\neq i}x_{\ell}y_{\ell}}_{=y_i \tilde x_i'\hat \beta - B_{ii} y_i^2} + \underbrace{\frac{y_i \tilde x_i' S_{xx} ^{-1} x_i x_i' S_{xx}^{-1}}{1 - P_{ii}} \sum_{\ell\neq i}x_{\ell}y_{\ell}}_{=B_{ii} y_i x_i'\hat \beta_{-i}} = y_i \tilde x_i'\hat \beta - B_{ii} \hat \sigma_i^2.
\end{align}
A similar combination of a change of variables argument and a leave-one-out estimator was used by \cite{powell1989semiparametric} in the context of weighted average derivatives. The JIVE estimators proposed by \cite{phillips1977bias} and \cite{angrist1999jackknife} also use a leave-one-out estimator, though without the change of variables.\footnote{ The object of interest in JIVE estimation is a \emph{ratio} of quadratic forms $\beta_1'S_{xx}\beta_2/\beta_2'S_{xx}\beta_2$ in the two-equation model $y_{ij} = x_{i}'\beta_j + \varepsilon_{ij}$ for $j=1,2$. When no covariates are present, using leave-out estimators of both the numerator and denominator of this ratio yields the JIVE1 estimator of \cite{angrist1999jackknife}.}
\begin{rem}
The $\{\hat\sigma_i^2 \}_{i=1}^n$ can also be used to construct an unbiased variance estimator
\begin{align}
\hat{\mathbb{V}}[\hat{\beta}] = S_{xx}^{-1}\left(\sum_{i=1}^{n} x_{i} x_{i}' \hat{\sigma}_{i}^{2} \right)S_{xx}^{-1}.
\end{align}
Section \ref{sec:InfLow} shows that $\hat{\mathbb{V}}[\hat{\beta}]$ can be used to perform asymptotically valid inference on linear contrasts in settings where existing Eicker-White estimators fail. Specifically, $\hat{\mathbb{V}}[\hat{\beta}]$ leads to valid inference under conditions where the MINQUE estimator of \cite{rao1970estimation} and the MINQUE-type estimator of \cite{cattaneo2017inference} do not exist \cite[see, e.g.,][]{horn1975estimating,verdier2016estimation}.
\end{rem}
\begin{rem}
The quantity $\hat{\mathbb{V}}[\hat{\beta}]$ is closely related to the HC2 variance estimator of \cite{mackinnon1985some}. While the HC2 estimator employs observation specific variance estimators $\hat \sigma^2_{i,\text{HC2}}=\frac{(y_i - x_i'\hat \beta)^2}{1-P_{ii}}$, $\hat{\mathbb{V}}[\hat{\beta}]$ relies instead on $\hat \sigma^2_i=\frac{y_i(y_i - x_i'\hat \beta)}{1-P_{ii}}$
\end{rem}
\begin{rem}\ensuremath{^{\textrm{th}}}\label{rem:cluster}
In some cases it may be important to allow dependence in the errors in addition to heteroscedasticity. A common case arises when the data are organized into mutually exclusive and independent ``clusters'' within which the errors may be dependent \citep{moulton1986random}. The same change of variables argument implies that an estimator of the form $\sum_{i=1}^{n}y_{i}\tilde{x}_{i}'\hat{\beta}_{-c\left(i\right)}$
will be unbiased in such settings, where $\hat{\beta}_{-c\left(i\right)}$ is the OLS estimator obtained after leaving out all observations in the cluster to which observation $i$ belongs.
\end{rem}
\subsection{Large Scale Computation} \label{sec:comp}
From \eqref{eq:estimator} and \eqref{eq:Pii}, computation of $\hat \theta$ relies on the values $\{B_{ii},P_{ii}\}_{i=1}^n$. Section \ref{sec:examples} provides some canonical examples where these quantities can be computed in closed form. When closed forms are unavailable, a number of options exist for accelerating computation. For example, in the empirical application of Section \ref{sec:application}, we make use of a preconditioned conjugate gradient algorithm suggested by \cite{koutis2011combinatorial} to compute exact leave-out variance decompositions in a two-way fixed effects model involving roughly one million observations and hundreds of thousands of parameters (see Appendix \ref{sec:Lambda} for details). However, in very large scale applications involving tens or hundreds of millions of parameters, exact computation of $\{B_{ii},P_{ii}\}_{i=1}^n$ is likely to become infeasible. Fortunately, it is possible to quickly approximate $\hat \theta$ in such settings using a variant of the random projection method introduced by \cite{achlioptas2001database}. We refer to this method as the Johnson-Lindenstrauss approximation (JLA) for its connection to the work of \cite{johnson1984extensions}.
JLA can be described by the following algorithm: fix a $p \in \mathbb{N}$ and generate the matrices $R_B, R_P \in \mathbb{R}^{p \times n}$, where $(R_B,R_P)$ are composed of mutually independent Rademacher random variables that are independent of the data, i.e., their entries take the values $1$ and $-1$ with probability $1/2$. Next decompose $A$ into $A=\frac{1}{2}( A_1'A_2 + A_2'A_1)$ for $A_1,A_2 \in \mathbb{R}^{n \times k}$ where $A_1=A_2$ if $A$ is positive semi-definite.\footnote{Interpretable choices of $A_1$ and $A_2$ are typically suggested by the structure of the problem; see, for instance, the discussion in \ensuremath{^{\textrm{th}}}\ref{ex:AKM} of Section \ref{sec:examples}.} Let
\begin{align}
\hat P_{ii} = \frac{1}{p}\norm*{ R_P X S_{xx}^{-1} x_i}^2
\quad \text{and} \quad
\hat B_{ii} =\frac{1}{p} \left(R_B A_1 S_{xx}^{-1} x_i\right)'\left(R_B A_2 S_{xx}^{-1} x_i\right)
\end{align}
where $X =(x_1,\dots,x_n)'$.
The Johnson-Lindenstrauss approximation to $\hat \theta$ is
\begin{align}
\hat \theta_{JLA} = \hat \beta'A\hat \beta - \sum_{i=1}^n \hat B_{ii} \hat \sigma_{i,JLA}^2,
\end{align}
where $\hat \sigma_{i,JLA}^2 = \frac{y_{i}\left(y_{i}-x_{i}'\hat{\beta}\right)}{1-\hat P_{ii}}\left(1- \frac{1}{p}\frac{3\hat P_{ii}^3 + \hat P_{ii}^2}{1-\hat P_{ii}}\right)$. The term $\frac{1}{p}\frac{3\hat P_{ii}^3 + \hat P_{ii}^2}{1-\hat P_{ii}}$ removes a non-linearity bias introduced by approximating $P_{ii}$.
Section \ref{sec:conspaper} establishes asymptotic equivalence between $\hat \theta_{JLA}$ and $\hat \theta$. Appendix \ref{sec:Lambda} discusses implementation details and numerically illustrates the trade-off between computation time and the bias introduced by JLA for different choices of $p$ under a range of sample sizes. Notably, we show
that JLA allows us to accurately compute a variance decomposition in a two-way fixed effects model with roughly 15 million parameters -- a scale comparable to the study of \cite{card2013workplace} -- in under an hour. A MATLAB package \citep{kssCODE} implementing both the exact and JLA versions of our estimator in the two-way fixed effects model is available online.
\subsection{Relation to Existing Approaches}
As discussed in Section \ref{sec:examples}, several literatures make use of bias corrections nominally predicated on homoscedasticity. A common ``homoscedasticity-only'' estimator takes the form
\begin{align}\label{eq:HO}
\hat \theta_{\text{HO}} = \hat{\beta}'A\hat{\beta}-\sum_{i=1}^{n} B_{ii} \hat{\sigma}_{\text{HO}}^{2}
\end{align}
where $\hat{\sigma}_{\text{HO}}^{2} = \frac{1}{n-k} \sum_{i=1}^n (y_i - x_i' \hat \beta)^2$ is the degrees-of-freedom corrected variance estimator. A sufficient condition for unbiasedness of $\hat \theta_{\text{HO}}$ is that there be no empirical covariance between $\sigma^{2}_{i}$ and $(B_{ii},P_{ii})$. This restriction is in turn implied by the special cases of \emph{homoscedasticity} where $\sigma_i^2$ does not vary with $i$ or \emph{balanced design} where $(B_{ii},P_{ii})$ does not vary with $i$. In general, however, this estimator will tend to be biased \citep[see, e.g.,][chapter 10, or Appendix \ref{app:rel}]{scheffe1959analysis}.
A second estimator, closely related to $\hat{\theta}$, relies upon a jackknife bias-correction \citep{quenouille1949approximate} of the plug-in estimator. This estimator can be written
\begin{align}
\hat \theta_{\text{JK}} = n \hat \theta_\text{PI} - \frac{n-1}{n} \sum_{i=1}^n \hat \theta_{\text{PI},-i}
\quad \text{where} \quad
\hat \theta_{\text{PI},-i} = \hat \beta_{-i}' A \hat \beta_{-i}.
\end{align}
In Appendix \ref{app:rel} we illustrate that jackknife bias-correction tends to over-correct and produce a first order bias in the opposite direction of the bias in the plug-in estimator. This is analogous to the upward bias in the jackknife estimator of $\mathbb{V}[\hat \beta]$ which was derived by \cite{efron1981} and shown by \cite{el2018can} to be of first order importance for inference with many Gaussian regressors.
There are several proposed adaptations of the jackknife to long panels that can decrease bias under stationarity restrictions on the regressors. Letting $t(i)\in\{1,...,T\}$ denote the time period in which an observation is observed, we can write the panel jackknife of \cite{hahn2004jackknife} as
\begin{align}
\hat \theta_\text{PJK} = T \hat \theta_{\text{PI}} - \frac{T-1}{T} \sum_{t=1}^T \hat \theta_{\text{PI},-t}
\quad \text{where} \quad
\hat \theta_{\text{PI},-t} = \hat \beta_{-t}' A\hat \beta_{-t}
\end{align}
and $\hat \beta_{-t} =(\sum_{i : t(i) \neq t} x_i x_i')^{-1} \sum_{i : t(i) \neq t} x_i y_i$
is the OLS estimator that excludes all observations from period $t$. \cite{dhaene2015split} propose a closely related split panel jackknife
\begin{align}
\hat \theta_\text{SPJK} = 2\hat \theta_{\text{PI}} - \frac{\hat \theta_{\text{PI},1} + \hat \theta_{\text{PI},2}}{2}
\quad \text{where} \quad
\hat \theta_{\text{PI},j} = \hat \beta_{j}' A\hat \beta_{j}
\end{align}
and $\hat \beta_{1}$ (and $\hat \beta_{2}$) are OLS estimators based on the first half (and the last half) of an even number of time periods. In Appendix \ref{app:rel}, we illustrate how short panels can lead these adaptations of the jackknife to produce first order biases in the opposite direction of the bias in the plug-in estimator.
\subsection{Finite Sample Properties}\label{fsample}
We now study the finite sample properties of the leave-out estimator $\hat \theta$ and its infeasible analogue $\theta^* = \hat \beta'A\hat \beta - \sum_{i=1}^n B_{ii}\sigma_i^2$, which uses knowledge of the individual error variances. First, we note that $\hat \theta$ is unbiased whenever each of the leave-one-out estimators $\hat \beta_{-i}$ exists, which can equivalently be expressed as the requirement that $\max_i P_{ii} <1$. This condition turns out to also be necessary for the existence of unbiased estimators, which highlights the need for additional restrictions on the model or sample whenever some leverages equal one.
\begin{lem}
\hangindent\leftmargini
\ensuremath{^{\textrm{th}}}\label{lem:unbiased}
\text{1.} If $\max_i P_{ii} <1$, then $\mathbb{E}[\hat \theta] = \theta$.
\begin{enumerate}
\setcounter{enumi}{1}
\item Unbiased estimators of $\theta = \beta'A\beta$ exist for all $A$ if and only if $\max_{i} P_{ii} < 1$.
\end{enumerate}
\end{lem}
Next, we show that when the errors are normal, the infeasible estimator $\theta^*$ is a weighted sum of a series of non-central $\chi^2$ random variables. This second result provides a useful point of departure for our asymptotic approximations and highlights the important role played by the matrix
\begin{align}
\tilde A = S_{xx}^{-1/2} A S_{xx}^{-1/2},
\end{align}
which encodes features of both the target parameter (which is defined by $A$) and the design matrix $S_{xx}$.
Let $\lambda_1,\dots,\lambda_r$ denote the non-zero eigenvalues of $\tilde A$, where $\lambda_1^2 \ge\dots \ge \lambda_r^2$ and each eigenvalue appears as many times as its algebraic multiplicity. We use $Q$ to refer to the corresponding matrix of orthonormal eigenvectors so that $\tilde A = Q D Q'$ where $D = \text{diag}(\lambda_1,\dots,\lambda_r)$. With these definitions we have
\begin{align}
\hat \beta'A\hat \beta = \sum_{\ell = 1}^r \lambda_\ell \hat b_\ell^2,
\end{align}
where $\hat b = (\hat b_1,\dots,\hat b_r)' = Q'S_{xx}^{1/2} \hat \beta$ contains $r$ linear combinations of the elements in $\hat \beta$.
The random vector $\hat b$ and the eigenvalues $\lambda_1,\dots,\lambda_r$ are central to both the finite sample distribution provided below in \ensuremath{^{\textrm{th}}}\ref{lem:fin} and the asymptotic properties of $\hat \theta$ as studied in Sections \ref{sec:InfLow}--\ref{sec:weak}. Each eigenvalue of $\tilde A$ can be thought of as measuring how strongly $\theta$ depends on a particular linear combination of the elements in $\beta$ relative to the difficulty of estimating that combination (as summarized by $S_{xx}^{-1}$). As discussed in Section \ref{sec:weak}, when a few of these eigenvalues are large relative to the others, a form of weak identification can arise.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:fin}
If $\varepsilon_{i} \sim \mathcal{N}(0, \sigma_i^2)$, then
\begin{enumerate}
\item $\hat b \sim \mathcal{N}\left( b , \mathbb{V}[\hat b ] \right)$ where $b = Q' S_{xx}^{1/2} \beta$,
\item $\theta^* = \sum_{\ell = 1}^{r} \lambda_\ell\left(\hat b_\ell^2-\mathbb{V}[\hat b_\ell ]\right)$
\end{enumerate}
\end{lem}
The distribution of $\theta^*$ is a sum of $r$ potentially dependent non-central $\chi^2$ random variables with non-centralities $b = (b_1,\dots,b_r)'$. In the special case of homoscedasticity $(\sigma_i^2 = \sigma^2)$ and no signal $(b=0)$ we have that $\hat b \sim \mathcal{N}\left( 0, \sigma^2 I_r \right)$, which implies that the distribution of $\theta^*$ is a weighted sum of $r$ \emph{independent} central $\chi^2$ random variables. The weights are the eigenvalues of $\tilde A$, therefore consistency of $\theta^*$ follows whenever the sum of the squared eigenvalues converges to zero. The next subsection establishes that the leave-out estimator remains consistent when a signal is present ($b\neq0$) and the errors exhibit unrestricted heteroscedasticity.
\subsection{Consistency}\label{sec:conspaper}
We now drop the normality assumption and provide conditions under which $\hat \theta$ remains consistent. To accommodate high dimensionality of the regressors we allow all parts of the model to change with $n$:
\begin{align}
y_{i,n} = x_{i,n}'\beta_n + \varepsilon_{i,n} && (i = 1,\ldots,n)
\end{align}
where $x_{i,n}\in\mathbb{R}^{k_n}$, $S_{xx,n}=\sum_{i=1}^{n}x_{i,n}x_{i,n}'$, $\mathbb{E}[\varepsilon_{i,n}]=0$, $\mathbb{E}[\varepsilon_{i,n}^2]=\sigma_{i,n}^{2}$ and $\theta_n=\beta_n'A_n\beta_n$ for some sequence of known non-random symmetric matrices $A_n\in\mathbb{R}^{k_n\times k_n}$ of rank $r_n$. By treating $x_{i,n}$ and $A_n$ as sequences of constants, all uncertainty derives from the disturbances $\left\{ \varepsilon_{i,n} : 1\le i \le n, n \ge 1 \right\}$. This \emph{conditional} perspective is common in the statistics literatures on ANOVA \citep{scheffe1959analysis,searle2009variance} and allows us to be agnostic about the potential dependency among the $\{x_{i,n}\}_{i=1}^n$ and $A_n$.\footnote{An \emph{unconditional} analysis might additionally impose distributional assumptions on $A_n$ and consider $\bar \theta = \beta'\mathbb{E}_{A_n}[A_n] \beta$ as the object of interest. The uncertainty in $\hat \theta - \bar \theta$ can always be decomposed into components attributable to $\hat \theta - \theta$ and $\theta - \bar \theta$. Because the behavior of $\theta - \bar \theta$ depends entirely on model choices, we leave such an analysis to future work.} Following standard practice we drop the $n$ subscript in what follows. All limits are taken as $n$ goes to infinity unless otherwise noted.
Our analysis makes heavy use of the following assumptions.
\begin{assumption}\ensuremath{^{\textrm{th}}}\label{ass:reg}
(i) $\max_i \left( \mathbb{E}[\varepsilon_i^4] + \sigma^{-2}_i \right) = O(1)$, (ii) there exist a $c <1$ such that $\max_i P_{ii} \le c$ for all $n$, and (iii) $\max_i (x_i'\beta)^2 = O(1)$.
\end{assumption}
Part $(i)$ of this condition limits the thickness of the tails in the error distribution, as is typically required for OLS estimation \citep[see, e.g.,][page 10]{cattaneo2017inference}. The bounds on $(x_i'\beta)^2$ and $P_{ii}$ imply that $\hat \sigma_i^2$ has bounded variance. Part $(iii)$ is a technical condition that can be relaxed to allow $\max_i (x_i'\beta)^2$ to
increase slowly with sample size as discussed further in Section \ref{sec:verify}.
From $(ii)$ it follows that $\frac{k}{n} \le c < 1$ for all $n$.
The following \ensuremath{^{\textrm{th}}}\nameref{lem:cons} establishes consistency of $\hat \theta$.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:cons}
If \ensuremath{^{\textrm{th}}}\ref{ass:reg} and one of the following conditions hold, then $\hat{\theta} - \theta \overset{p}{\rightarrow} 0$.
\begin{enumerate}
\item[(i)] $A$ is positive semi-definite, $\theta = \beta'A\beta= O(1)$, and $\text{trace}(\tilde A^2) = \sum_{\ell =1}^r \lambda_{\ell}^2= o(1)$.
\item[(ii)] $A=\frac{1}{2}( A_1'A_2 + A_2'A_1)$ where $\theta_1 = \beta'A_1'A_1\beta$ and $\theta_2 = \beta'A_2'A_2\beta$ satisfy (i).
\end{enumerate}
\end{lem}
The first condition of \ensuremath{^{\textrm{th}}}\ref{lem:cons} establishes consistency of variance components given boundedness of $\theta$ and a joint condition on the design matrix $S_{xx}$ and the matrix $A$. The second condition shows that consistency of covariance components follows from consistency of variance components that dominate them via the Cauchy-Schwarz inequality, i.e., $\theta^2 = (\beta'A_1'A_2\beta)^2 \le \theta_1 \theta_2$. In several of the examples discussed in the next section, $\text{trace}(\tilde A^2)$ is of order $r/n^2$, which is necessarily small in large samples. A more extensive discussion of primitive conditions that yield $\text{trace}(\tilde A^2) = o(1)$ is provided in Section \ref{sec:verify}.
We conclude this section by establishing asymptotic equivalence between the leave-out estimator $\hat \theta$ and its approximation $\hat \theta_{JLA}$ under the condition that $p^4$ is large relative to sample size.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:JLA}
If \ensuremath{^{\textrm{th}}}\ref{ass:reg} is satisfied, $n/p^4 = o(1)$, $\mathbb{V}[\hat \theta]^{-1} = O(n)$, and one of the following conditions hold, then $\mathbb{V}[\hat \theta]^{-1/2}{(\hat \theta_{JLA} - \hat \theta - \mathrm{B}_p)}{} = o_p(1)$ where $\abs{\mathrm{B}_p} \le \frac{1}{p} \sum_{i=1}^n P_{ii}^2 \abs{B_{ii}} \sigma_i^2$.
\begin{enumerate}
\item[(i)] $A$ is positive semi-definite and $\mathbb{E}[\hat \beta'A\hat \beta] - \theta = \sum_{i=1}^n B_{ii} \sigma_i^2 = O(1)$.
\item[(ii)] $A=\frac{1}{2}( A_1'A_2 + A_2'A_1)$ where $\theta_1 = \beta'A_1'A_1\beta$ and $\theta_2 = \beta'A_2'A_2\beta$ satisfy (i) and $\frac{\mathbb{V}[\hat \theta_1]\mathbb{V}[\hat \theta_2]}{n\mathbb{V}[\hat \theta]^2} = O(1)$.
\end{enumerate}
\end{lem}
\ensuremath{^{\textrm{th}}}\ref{lem:JLA} requires that $\hat \theta$ is not super-consistent and that the bias in the plug-in estimator is asymptotically bounded, assumptions which can be shown to be satisfied in the examples introduced in the next section. For variance components, the \ensuremath{^{\textrm{th}}}\nameref{lem:JLA} characterizes an approximation bias $\mathrm{B}_p$ in $\hat \theta_{JLA}$ of order $1/p$ and provides an interpretable bound on $\mathrm{B}_p$: the approximation bias is at most $1/p$ times the bias in the plug in estimator $\hat \beta 'A \hat \beta$. For covariance components, asymptotic equivalence follows when the variance components defined by $A_1'A_1$ and $A_2'A_2$ do not converge at substantially slower rates than $\hat \theta$. Under this condition, the approximation bias is at most $1/p$ times the average of the biases in the plug in estimators $\hat \beta 'A_1'A_1 \hat \beta$ and $\hat \beta 'A_2'A_2 \hat \beta$.
These bounds on the approximation bias suggests that a $p$ of a few hundred should suffice for point estimation. However, unless $n/p^2 = o(1)$, the resulting approximation bias needs to be accounted for when conducting inference. Specifically, one can lengthen the tails of the confidence sets proposed in Sections \ref{sec:dist} and \ref{sec:variance} by $\frac{1}{p} \sum_{i=1}^n \hat P_{ii}^2 \abs{\hat B_{ii}} \hat \sigma_{i,JLA}^2$ when relying on JLA.
\section{Examples}\label{sec:examples}
We now consider four commonly encountered empirical examples where our proposed estimation strategy provides an advantage over existing methods.
\begin{example}[Coefficient of determination]\ensuremath{^{\textrm{th}}}\label{ex:R2}
\textit{}
Sewall \citet{wright1921correlation} proposed measuring the explanatory power of a linear model using the coefficient of determination. When $x_i$ includes an intercept, the object of interest and its corresponding plug-in estimator can be written
\begin{align}
R^2 &= \frac{\beta' A \beta }{\beta' A \beta + \frac{1}{n}\sum_{i=1}^{n} \sigma_i^2} = \frac{\sigma^2_{X\beta}}{\sigma^2_y}
\quad \text{and} \quad
\hat R_\text{PI}^2 = \frac{\hat \beta' A \hat \beta}{\frac{1}{n}\sum_{i=1}^{n} (y_i - \bar y)^{2}} = \frac{\hat \sigma^2_{X\beta,\text{PI}}}{\hat \sigma^2_y}
\shortintertext{where}
A &= \frac{1}{n}\sum_{i=1}^{n} (x_i-\bar x)(x_i-\bar x)', \quad \bar x = \frac{1}{n} \sum_{i=1}^n x_i, \quad \bar y = \frac{1}{n} \sum_{i=1}^n y_i.
\end{align}
\cite{theil1961economic} noted that the plug-in estimator of $\sigma^2_{X\beta}$ is biased and proposed an adjusted $R^2$ measure that utilizes the homoscedasticity-only estimator in \eqref{eq:HO}. The above choice of $A$ yields $B_{ii} = \frac{1}{n} (P_{ii} - \frac{1}{n}),$ which implies $\sum_{i=1}^n B_{ii} = \frac{k-1}{n}$. Hence, Theil's proposal can be written
\begin{align}
\hat R^2_\text{adj} = \frac{\hat \sigma^2_{X\beta,\text{HO}}}{\hat \sigma^2_y} = \frac{\hat \beta' A \hat \beta - \frac{k-1}{n}\hat{\sigma}_{\text{HO}}^{2}}{\hat \sigma^2_y}.
\end{align}
A rearrangement gives the familiar representation $\frac{1-\hat R^2_\text{adj}}{1-\hat R_{\text{PI}}^2} = \frac{n-1}{n-k}$ which highlights that the adjusted estimator of $R^2$ relates to the unadjusted one through a degrees-of-freedom correction.
The leave-out estimator of $\sigma^2_{X\beta}$ allows for unrestricted heteroscedasticity and can be found by noting that $\tilde x_i = A S_{xx}^{-1} x_i = \frac{1}{n}(x_i - \bar x)$, which yields
\begin{align}
\hat R^2 = \frac{\hat \sigma_{X\beta}^2}{\hat \sigma^2_y}
\quad \text{where} \quad
\hat \sigma^2_{X\beta}=\frac{1}{n}\sum_{i=1}^{n}y_i (x_{i} - \bar x)'\hat\beta_{-i}.
\end{align}
In general, this estimator does not have an interpretation in terms of degrees-of-freedom corrections. Instead, the explanatory power of the linear model is assessed using the empirical covariance between leave-one-out predictions $(x_{i} - \bar x)'\hat \beta_{-i}$ and the left out observation $y_i$.
\end{example}
\begin{example}[Analysis of covariance]\ensuremath{^{\textrm{th}}}\label{ex:ANOVA}
\textit{}
Since the work of \citet{fisher1925statistical}, it has been common to summarize the effects of experimentally assigned treatments on outcomes with estimates of variance components. Consider a dataset comprised of observations on $N$ groups with $T_{g}$ observations in the ${g}$-th group. The ``analysis of covariance'' model posits that outcomes can be written
\begin{align}\label{eq:panel}
y_{{g} t} =\alpha_{{g}} + x_{{g} t}'\delta + \varepsilon_{{g} t} && ({g}=1,\dots,N, \ t = 1,\dots,T_{g} \ge 2),
\end{align}
where $\alpha_g$ is a group effect and $x_{gt}$ is a vector of strictly exogenous covariates.
A prominent example comes from \cite{chetty2011does} who study the adult earnings $y_{gt}$ of $n = \sum_{{g}=1}^N T_{g}$ students assigned experimentally to one of $N$ different classrooms. Each student also has a vector of predetermined background characteristics $x_{gt}$. The variability in student outcomes attributable to classrooms can be written:
\begin{align}
\sigma_{\alpha}^{2}=\frac{1}{n}\sum_{{g}=1}^{N} T_{g} \left(\alpha_{g}-\bar{\alpha}\right)^{2}
\end{align}
where $\bar{\alpha}=\frac{1}{n}\sum_{{g}=1}^{N} T_{g} \alpha_{{g}}$ gives the (enrollment-weighted) mean classroom effect.
This model and object of interest can written in the notation of the preceding section ($y_i = x_i'\beta + \varepsilon_{i}$ and $\sigma_{\alpha}^{2} = \beta'A\beta$) by letting $i = i(g,t)$ where $i(\cdot,\cdot)$ is bijective with inverse denoted $(g(\cdot),t(\cdot))$, $y_i= y_{{g} t}$, $\varepsilon_{i} = \varepsilon_{{g} t}$,
\begin{align}
x_i &= (d_i',x_{{g} t}')',
\quad
\beta = (\alpha',\delta')', \quad \alpha = (\alpha_1,\dots,\alpha_N)',
\quad
d_i = (\mathbf{1}_{\{{g}=1\}},\dots,\mathbf{1}_{\{{g}=N\}})',
\shortintertext{and}
A &= \mbox{\scriptsize $\begin{bmatrix} A_{dd} & 0 \\ 0 & 0 \end{bmatrix} $}
\quad \text{where} \quad
A_{dd}=\frac{1}{n}\sum_{i=1}^{n} (d_{i} - \bar d) (d_{i} - \bar d)', \quad \bar d=\frac{1}{n} \sum_{i=1}^{n} d_{i}.
\end{align}
\cite{chetty2011does} estimate $\sigma_{\alpha}^{2}$ using a random effects ANOVA estimator \cite[see e.g.,][]{searle2009variance} which is of the homoscedasticity-only type given in \eqref{eq:HO}. As discussed in Section \ref{sec:unbiased} and Appendix \ref{app:rel}, this estimator is in general first order biased when the errors are heteroscedastic and group sizes are unbalanced.
\noindent \textbf{\emph{Special Case: No Common Regressors}}
When there are no common regressors ($x_{{g} t}=0$ for all ${g},t$), the leave-out estimator of $\sigma_{\alpha}^{2}$ has a particularly simple representation:
\begin{align}\label{eq:AP}
\hat{\sigma}_{\alpha}^{2}
&= \dfrac{1}{n}\sum_{g=1}^{N} \left( T_g \left(\hat \alpha_{g}-\hat{\bar{\alpha}}\right)^{2} - \left(1-\frac{T_g}{n}\right) \hat \sigma_g^2 \right)
\quad \text{for} \quad \hat \sigma_g^2 = \frac{1}{T_g-1} \sum_{t=1}^{T_g} (y_{gt} - \hat \alpha_g)^2,
\end{align}
where $\hat \alpha_{g}= \frac{1}{T_g}\sum_{t=1}^{T_g} y_{gt}$, and $\hat{\bar{\alpha}}= \frac{1}{n}\sum_{{g}=1}^{N} T_{g} \hat \alpha_{{g}}$. This representation shows that if the model consists only of group specific intercepts, then the leave-out estimator relies on group level degrees-of-freedom corrections. The statistic in \eqref{eq:AP} was analyzed by \cite{akritas2004heteroscedastic} in the context of testing the null hypothesis that $\sigma_\alpha^2=0$ while allowing for heteroscedasticity at the group level.
\noindent \textbf{\emph{Covariance Representation}}
Another instructive representation of the leave-out estimator is in terms of the empirical covariance
\begin{align}
\hat{\sigma}_{\alpha}^{2} &= \sum_{i=1}^{n} y_i \tilde{d}_i'\hat \alpha_{-i}
\quad \text{where} \quad \hat \beta_{-i} = (\hat \alpha_{-i}',\hat \delta_{-i}').
\end{align}
The generalized regressor $\tilde{d}_i$ can be described as follows: if there are no common regressors then $\tilde{d}_i = \frac{1}{n} (d_i - \bar d)$, which is analogous to \ensuremath{^{\textrm{th}}}\ref{ex:R2}. If the model includes common regressors then $\tilde d_i = \frac{1}{n} \left( (d_i - \bar d) - \hat \varGamma'(x_{g(i)t(i)} - \bar x_{g(i)}) \right)$ where $\bar x_g = \frac{1}{T_{g}} \sum_{t=1}^{T_g} x_{gt}$ and $\hat \varGamma$ is the coefficient vector from an instrumental variables (IV) regression of $d_i - \bar d$ on $x_{g(i)t(i)} - \bar x_{g(i)}$ using $x_{g(i)t(i)}$ as an instrument. The IV residual $\tilde{d}_i$ is uncorrelated with $x_{{g}(i) t(i)}$ and the covariance between ${d_i}$ and $\tilde d_i$ is $A_{dd}$, which ensures that the empirical covariance between $y_i = d_i'\alpha + x_{{g}(i) t(i)}'\delta + \varepsilon_i$ and the generalized prediction $\tilde{d_i}'\hat \alpha_{-i}$ is an unbiased estimator of ${\sigma}_{\alpha}^{2}$.
\end{example}
\begin{example}[Random coefficients]\ensuremath{^{\textrm{th}}}\label{ex:RC}
\textit{}
Group memberships are often modeled as influencing slopes in addition to intercepts \citep{kuh1959validity,hildreth1968some,raudenbush2002hierarchical,arellano2011identifying,graham2012identification,graham2016quantile}. Consider the following ``random coefficient'' model:
\begin{align}\label{random_coeff}
y_{{g} t} = \alpha_{g} + z_{{g} t}\gamma_{g} + \varepsilon_{{g} t} && ({g}=1,\dots,N, \ t=1,\dots,T_{g} \ge 3).
\end{align}
An influential example comes from \cite{raudenbush1986hierarchical}, who model student mathematics
scores as a ``hierarchical'' linear function of socioeconomic status (SES) with school-specific intercepts $(\alpha_{g} \in \mathbb{R})$ and slopes $(\gamma_{g} \in \mathbb{R})$. Letting $\bar{\gamma}=\frac{1}{n}\sum_{{g}=1}^{N} T_{g} \gamma_{{g}}$ for $n = \sum_{{g}=1}^N T_{g}$, the student-weighted variance of slopes can be written:
\begin{align}
\sigma_{\gamma}^{2}=\frac{1}{n}\sum_{{g}=1}^N T_{g}\left(\gamma_{{g}}-\bar{\gamma}\right)^{2}.
\end{align}
In the notation of the preceding section we can write $y_i = x_i'\beta + \varepsilon_{i}$ and $\sigma_\gamma^2 = \beta'A\beta$ where
\begin{align}
x_i = (d_i', d_i' z_{{g} t})', \qquad \beta = (\alpha',\gamma')', \qquad \gamma = (\gamma_1,\dots,\gamma_N)',
\qquad A = \mbox{\scriptsize $\begin{bmatrix} A_{dd} & 0 \\ 0 & 0 \end{bmatrix} $}
\end{align}
for $y_i$, $\varepsilon_i$, $d_i$, $A_{dd}$, and $\alpha$ as in the preceding example.
\cite{raudenbush1986hierarchical} use a maximum likelihood estimator of $\sigma_{\gamma}^{2}$ predicated upon normality and homoscedastic errors. \cite{swamy1970efficient} considers an estimator of $\sigma^{2}_{\gamma}$ that relies on group-level degrees-of-freedom corrections and is unbiased when the error variance is allowed to vary at the group level, but not with the level of $z_{{g} t}$. By contrast, the leave-out estimator is unbiased under arbitrary patterns of heteroscedasticity.
\noindent \textbf{\emph{Covariance Representation} }
The leave-out estimator can be represented in terms of the empirical covariance
\begin{align}
\hat \sigma_{\gamma}^{2}=\sum_{i=1}^{n} y_i \tilde z_i \tilde d_i'\hat \gamma_{-i}
\quad \text{where} \quad \tilde d_i = \frac{1}{n}(d_i - \bar d), \quad \tilde z_i = \frac{z_{g(i)t(i)}-\bar z_{g(i)}}{\sum_{t=1}^{T_{g(i)}}( z_{g(i)t}-\bar z_{g(i)} )^2},
\end{align}
and $\bar z_{g} = \frac{1}{T_g}\sum_{t=1}^{T_{g}} z_{gt}$. Demeaning $z_{g(i)t(i)}$ at the group level makes $\tilde d_i \tilde z_i$ uncorrelated with $d_i$ and scaling by the group variability in $z_{g(i)t}$ ensures that the covariance between $\tilde d_i \tilde z_i$ and $d_i z_{g(i)t(i)}$ is $A_{dd}$. This implies that the empirical covariance between $y_i = d_i'\alpha + z_{g(i)t(i)} d_i'\gamma + \varepsilon_i$ and the generalized prediction $\tilde z_i \tilde d_i'\hat \gamma_{-i}$ is an unbiased estimator of $\sigma_{\gamma}^{2}$.
\end{example}
\begin{example}[Two-way fixed effects]\ensuremath{^{\textrm{th}}}\label{ex:AKM}
\textit{}
Economists often study settings where units possess two or more group memberships, some of which can change over time. A prominent example comes from \cite{abowd1999high} (henceforth AKM) who propose a panel model of log wage determination that is additive in worker and firm fixed effects. This so-called ``two-way'' fixed effects model takes the form:
\begin{align}\label{eq:AKM}
y_{{g} t} = \alpha_{{g}} + \psi_{j({g},t)} + x_{{g} t}'\delta + \varepsilon_{{g} t} && ({g}=1,\dots,N, \ t=1,\dots,T_{g} \ge 2)
\end{align}
where the function $j(\cdot,\cdot):\{ 1,\dots,N\}\times\{1,\dots, \max_g T_{g}\} \rightarrow \{ 0,\dots,J\} $ allocates each of $n = \sum_{{g} =1}^N T_{g}$ person-year observations to one of $J+1$ firms. Here $\alpha_{{g}}$ is a ``person effect'', $\psi_{j({g},t)}$ is a ``firm effect'', $x_{{g} t}$ is a time-varying covariate, and $\varepsilon_{{g} t}$ is a time-varying error. In this context, the mean zero assumption on the errors $\varepsilon_{{g} t}$ can be thought of as requiring both the common covariates $x_{gt}$ and the firm assignments $j(\cdot,\cdot)$ to obey a strict exogeneity condition.
Interest in such models often centers on understanding how much of the variability in log wages is attributable to firms \citep[see, e.g.,][]{card2013workplace,song2015firming}. AKM summarize the firm contribution to wage inequality via the following two parameters:
\begin{align}
\sigma_{\psi}^{2}=\frac{1}{n}\sum_{{g}=1}^{N}\sum_{t=1}^{T_{g}}\left(\psi_{j\left({g},t\right)}-\bar{\psi}\right)^{2}
\quad \text{and} \quad
\sigma_{\alpha,\psi}=\frac{1}{n}\sum_{{g}=1}^{N}\sum_{t=1}^{T_{g}}\left(\psi_{j\left({g},t\right)}-\bar{\psi}\right)\alpha_{{g}}
\end{align}
where $\bar \psi = \frac{1}{n}\sum_{{g}=1}^{N} \sum_{t=1}^{T_{g}} \psi_{j({g},t)}$. The variance component $\sigma_{\psi}^{2}$ measures the contribution of firm wage variability to inequality, while the covariance component $\sigma_{\alpha,\psi}$ measures the additional contribution of systematic sorting of high wage workers to high wage firms.
To represent this model and the corresponding objects of interest in the notation of the preceding section ($y_i = x_i'\beta + \varepsilon_{i}$, $\sigma_{\psi}^{2}= \beta'A_\psi \beta$, and $\sigma_{\alpha,\psi}= \beta'A_{\alpha,\psi}\beta$), let
\begin{align}
x_i = (d_i',f_i',x_{{g} t}')', \ \beta = (\alpha',\psi',\delta')', \ \alpha = (\alpha_1,\dots,\alpha_N)' + \mathbf{1}_N'\psi_0, \ \psi = (\psi_1\,\dots,\psi_J)' - \mathbf{1}_J'\psi_0,
\end{align}
for $y_i$, $\varepsilon_i$, and $d_i$ as in the preceding examples, $f_i = (\mathbf{1}_{\{j({g},t)=1\}},\dots,\mathbf{1}_{\{j({g},t)=J\}})',$
\begin{align}
A_\psi
&= \mbox{\scriptsize $\begin{bmatrix} 0 & 0 & 0 \\ 0 & A_{ff} & 0 \\ 0 & 0 & 0 \end{bmatrix}$ }
\quad \text{where} \quad
A_{ff}=\frac{1}{n}\sum_{i=1}^{n} (f_{i}-\bar f)(f_{i} - \bar f)', \quad
\bar f=\frac{1}{n}\sum_{i=1}^{n} f_{i} ,
\shortintertext{and}
A_{\alpha,\psi}
&=\frac{1}{2} \mbox{\scriptsize $\begin{bmatrix} 0 & A_{df} & 0 \\ A_{df}' & 0 & 0 \\ 0 & 0 & 0\end{bmatrix}$}
\quad \text{where} \quad
A_{df}=\frac{1}{n}\sum_{i=1}^{n} (d_{i} - \bar d)(f_{i} - \bar f)'.
\end{align}
Computation of the Johnson-Lindenstrauss approximation can be facilitated using the representations $A_\psi = A_f'A_f$ and $A_{\alpha,\psi} = \frac{1}{2}(A_d'A_f + A_f'A_d)$ where
\begin{align}
A_f' = \mbox{\scriptsize $ \frac{1}{\sqrt{n}}\begin{bmatrix} 0 & 0 & 0 \\ f_1 - \bar f & \dots & f_n-\bar f \\ 0 & 0 & 0 \end{bmatrix}$} \quad \text{and} \quad A_d' = \mbox{\scriptsize$\frac{1}{\sqrt{n}}\begin{bmatrix} d_1 - \bar d & \dots & d_n-\bar d \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}$ }.
\end{align}
Addition and subtraction of $\psi_0$ in $\beta$ amounts to the normalization, $\psi_0 = 0$, which has no effect on the variance components of interest. As \cite{abowd1999high,abowd2002computing} note, least squares estimation of \eqref{eq:AKM} requires one normalization of the $\psi$ vector within each set of firms connected by worker mobility.
For simplicity, we assume all firms are connected so that only a single normalization is required.\footnote{\cite{bonhomme2019distributional} study a closely related model where workers and firms each belong to one of a finite number of types and each pairing of worker and firm type is allowed a different mean wage. These mean wage parameters are shown to be identified
when each worker type moves between each firm type with positive probability, enabling estimation even when many firms are not connected.}
\noindent \textbf{\emph{Covariance Representation}}
\cite{abowd1999high} estimated $\sigma_{\psi}^{2}$ and $\sigma_{\alpha,\psi}$ using the naive plug-in estimators $\hat{\beta}'A_{\psi}\hat{\beta}$ and $ \hat{\beta}'A_{\alpha,\psi}\hat{\beta}$ which are, in general, biased. \cite{andrews2008high} proposed the ``homoscedasticity-only'' estimators of \eqref{eq:HO}. These estimators are unbiased when the errors $\varepsilon_i$ are independent and have common variance. \cite{bonhomme2019distributional} propose a two-step estimation approach that is consistent in the presence of heteroscedasticity when the support of firm wage effects is restricted to a finite number of values and each firm grows large with the total sample size $n$. Our leave-out estimators, which avoid both the homoscedasticity requirement on the errors and any cardinality restrictions on the support of the firm wage effects, can be written compactly as covariances taking the form
\begin{align}\label{eq:varest}
\hat{\sigma}_{\psi}^{2} =\sum_{i=1}^{n} y_{i} x_{i}' S_{xx}^{-1} A_{\psi} \hat{\beta}_{-i},
\qquad
\hat{\sigma}_{\alpha,\psi} =\sum_{i=1}^{n} y_{i} x_{i}' S_{xx}^{-1} A_{\alpha,\psi} \hat{\beta}_{-i}.
\end{align}
Notably, these estimators are unbiased whenever the leave out estimator $\hat \beta_{-i}$ can be computed, regardless of the distribution of firm sizes.
\noindent \textbf{\emph{Special Case: Two time periods}}
A simpler representation of $\hat{\sigma}_{\psi}^{2}$ is available in the case where only two time periods are available and no common regressors are present ($T_{g} =2$ and $x_{{g} t} =0$ for all ${g},t$). Consider this model in first differences
\begin{align}
\label{fd_model}
\Delta y_{{g}} = \Delta f_{{g}}' \psi + \Delta \varepsilon_{{g}}
&&
({g}=1,\dots,N)
\end{align}
where $\Delta y_{{g}}=y_{{g} 2}-y_{{g} 1}$, $\Delta \varepsilon_{{g}} = \varepsilon_{{g} 2} - \varepsilon_{{g} 1}$, and
$\Delta f_{{g}} = f_{i(g,2)} - f_{i(g,1)}$. The leave-out estimator of $\sigma_{\psi}^2$ applied to this differenced representation of the model is:
\begin{align}
\hat \sigma_{\psi}^2 =\sum_{{g} =1}^N \Delta y_{{g} } \Delta \tilde f_g'\hat{\psi}_{-{g}}
\quad \text{where} \quad
\Delta \tilde f_g = A_{ff} S_{\Delta f\Delta f}^{-1} \Delta f_g.
\end{align}
Note that the quantities $S_{\Delta f\Delta f}$ and $\hat{\psi}_{-{g}}$ correspond respectively to $S_{xx}$ and $\hat{\beta}_{-i}$ in the first differenced model.
\begin{rem}\ensuremath{^{\textrm{th}}}\label{rem:nocluster}The leave-out representation above reveals that $\hat{\sigma}_{\psi}^{2}$ is not only unbiased under arbitrary heteroscedasticity and design unbalance, but also under arbitrary correlation between $\varepsilon_{{g} 1}$ and $\varepsilon_{{g} 2}$. The same can be shown to hold for $\hat{\sigma}_{\alpha,\psi}$. Furthermore, this representation highlights that $\hat{\sigma}_{\psi}^{2}$ only depends upon observations with $\Delta f_{{g} } \neq 0$ (i.e., firm ``movers'').
\end{rem}
\end{example}
\section{Inference on Quadratic Forms of Fixed Rank}\label{sec:InfLow}
While the previous section emphasized variance components where the rank $r$ of $A$ was increasing with sample size, we first study the case where $r$ is fixed. Problems of this nature often arise when testing a few linear restrictions or when conducting inference on linear combinations of the regression coefficients, say $v'\beta$. In the case of two-way fixed effects models of wage determination, the quantity $v'\beta$ might correspond to the difference in mean values of firm effects between male and female workers \citep{card2015bargaining} or to the coefficient from a projection of firm effects onto firm size \citep{bloom2018disappearing}.\fxnote{does bloom actually run such a regression?} A third use case, discussed at length by \cite{cattaneo2017inference}, is where $v'\beta$ corresponds to a linear combination of a few common coefficients in a linear model with high dimensional fixed effects that are regarded as nuisance parameters.
To characterize the limit distribution of $\hat \theta$ when $r$ is small, we rely on a representation of $\theta$ as a weighted sum of squared linear combinations of the data: $\hat \theta = \sum_{\ell = 1}^{r} \lambda_\ell\left(\hat b_\ell^2-\hat{\mathbb{V}}[\hat b_\ell ]\right)$ where
\begin{align}
\hat b &= \sum_{i=1}^n w_i y_i
\quad \text{and} \quad
\hat{\mathbb{V}}[\hat b] = \sum_{i=1}^n w_{i} w_{i}' \hat \sigma_{i}^2
\end{align}
for $w_i = (w_{i1},\dots,w_{ir})'= Q' S_{xx}^{-1/2} x_i$. The following theorem characterizes the asymptotic distribution of $\hat \theta$ while providing conditions under which $\hat b$ is asymptotically normal and $\hat{\mathbb{V}}[\hat b]$ is consistent.
\begin{thm}\ensuremath{^{\textrm{th}}}\label{thm2}
If \ensuremath{^{\textrm{th}}}\ref{ass:reg} holds, $r$ is fixed, and $\max_i w_i'w_i = o(1)$, then
\begin{enumerate}
\item $\mathbb{V}[\hat b]^{-1/2}(\hat b - b) \xrightarrow{d} \mathcal{N}\left( 0 , I_r \right)$ where $b = Q' S_{xx}^{1/2} \beta$,
\item ${\mathbb{V}}[\hat b]^{-1} \hat{\mathbb{V}}[\hat b] \xrightarrow{p} I_r$,
\item $\hat \theta = \sum_{\ell = 1}^{r} \lambda_\ell\left(\hat b_\ell^2-\mathbb{V}[\hat b_\ell]\right) + o_p(\mathbb{V}[\hat \theta]^{1/2})$,
\end{enumerate}
\end{thm}
The high-level requirement of this theorem that $\max_i w_i'w_i = o(1)$ is a Lindeberg condition ensuring that no observation is too influential. One can think of $\max_i w_i'w_i$ as measuring the inverse effective sample size available for estimating $b$: when the weights are equal across $i$, the equality $\sum_{i=1}^n w_i w_i' = I_r$ implies that $w_{i\ell}^2 = \frac{1}{n}$. Since $\frac{1}{n} \sum_{i=1}^n w_i'w_i = \frac{r}{n}$, the requirement that $\max_i w_i'w_i = o(1)$ is implied by a variety of primitive conditions that limit how far a maximum is from the average \cite[see, e.g.,][Appendix A.1]{anatolyev2012inference}. Note that \ensuremath{^{\textrm{th}}}\ref{thm2} does not apply to settings where $r$ is proportional to $n$ because $\max_i w_i'w_i \ge \frac{r}{n}$.
In the special case where $A = vv'$ for some non-random vector $v$, \ensuremath{^{\textrm{th}}}\ref{thm2} establishes that the variance estimator $\hat{\mathbb{V}}[\hat \beta] = S_{xx}^{-1} \left(\sum_{i=1}^n x_i x_i' \hat \sigma_i^2\right)S_{xx}^{-1}$ enables consistent inference on the linear combination $v'\beta$ using the approximation
\begin{align}
\frac{v'(\hat \beta - \beta)}{\sqrt{v'\hat{\mathbb{V}}[\hat \beta]v}} \xrightarrow{d} \mathcal{N}(0,1). \label{lincom}
\end{align}
To derive this result we assumed that $\max_i P_{ii} \le c$ for some $c<1$, whereas standard Eicker-White variance estimators generally require that $\max_i P_{ii} \rightarrow 0$ and \cite{cattaneo2017inference} establish an asymptotically valid approach to inference in settings where $\max_i P_{ii} \le 1/2$. Thus $\hat{\mathbb{V}}[\hat \beta]$ leads to valid inference under weaker conditions than existing versions of Eicker-White variance estimators.
\begin{rem}
\ensuremath{^{\textrm{th}}}\ref{thm2} extends classical results on hypothesis testing of a few linear restrictions, say, $H_0 : R\beta=0$, to allow for many regressors and heteroscedasticity. A convenient choice of $A$ for testing purposes is $\frac{1}{r} R' (R S_{xx}^{-1} R')^{-1} R$ where $r$, the rank of $R \in \mathbb{R}^{r \times k}$, is fixed. Under $H_0$, the asymptotic distribution of $\hat \theta$ is an equally weighted sum of $r$ central $\chi^2$ random variables. This distribution is known up to $\mathbb{V}[\hat b ]$ and a critical value can be found through simulation. For a recent contribution to this literature, see \cite{anatolyev2012inference} who allows for many regressors but considers the special case of homoscedastic errors.
\end{rem}
\section{Inference on Quadratic Forms of Growing Rank}\label{sec:dist}
We now turn to the more challenging problem of conducting inference on $\theta$ when $r$ increases with $n$, as in the examples discussed in Section \ref{sec:examples}. These results also enable tests of many linear restrictions. For example, in a model of gender-specific firm effects of the sort considered by \cite{card2015bargaining}, testing the hypothesis that men and women face identical sets of firm fixed effects entails as many equality restrictions as there are firms.
\subsection{Limit Distribution}
In order to describe the result we introduce $\check x_i = \sum_{\ell=1}^n M_{i\ell} \frac{B_{\ell\ell}}{1-P_{\ell \ell}} x_\ell$ where $M_{i\ell} = \mathbf{1}_{\{i=\ell\}} - x_i S_{xx}^{-1} x_\ell$. Note that $\check x_i $ gives the residual from a regression of $\frac{B_{ii}}{1-P_{ii}} x_i$ on $x_i$. Therefore, $\check x_i =0$ when the regressor design is balanced. The contribution of $\check x_i$ to the behavior of $\hat \theta$ is through the estimation of $\sum_{i=1}^n B_{ii} \sigma_i^2$, which can be ignored in the case where the rank of $A$ is bounded. When the rank of $A$ is large, as implied by condition $(ii)$ of \ensuremath{^{\textrm{th}}}\ref{thm3} below, this estimation error can resurface in the asymptotic distribution. One can think of the eigenvalue ratio in $(ii)$ as the inverse effective rank of $\tilde A$: when all the eigenvalues are equal $\frac{ \lambda_1^2}{\sum_{\ell=1}^r \lambda_\ell^2} = \frac{1}{r}$.
\begin{thm}\ensuremath{^{\textrm{th}}}\label{thm3}
Recall that $\tilde x_i = A S_{xx}^{-1} x_i$ where $\hat \theta = \sum_{i=1}^n y_i \tilde x_i'\hat \beta_{-i}$. If \ensuremath{^{\textrm{th}}}\ref{ass:reg} holds and the following conditions are satisfied
\begin{align}
(i) \ \mathbb{V}[\hat \theta]^{-1} \max_i \left( (\tilde x_i'\beta)^2 + (\check x_i'\beta)^2 \right) = o(1), \quad (ii) \ \frac{ \lambda_1^2}{\sum_{\ell=1}^r \lambda_\ell^2} = o(1),
\end{align}
then $\mathbb{V}[\hat \theta]^{-1/2} (\hat \theta - \theta ) \xrightarrow{d} \mathcal{N}(0,1)$.
\end{thm}
The proof of \ensuremath{^{\textrm{th}}}\ref{thm3} relies on a variation of Stein's method developed in \cite{soelvsten2017robust} and a representation of $\hat \theta$ as a second order U-statistic, i.e.,
\begin{align}\label{eq:Ustat}
\hat \theta = \sum_{i=1}^n \sum_{\ell \neq i} C_{i\ell} y_i y_\ell
\end{align}
where $C_{i\ell} = B_{i\ell} - 2^{-1} M_{i\ell}\left(M_{ii}^{-1} B_{ii} + M_{\ell\ell}^{-1} B_{\ell\ell}\right)$ and $B_{i\ell} = x_i' S_{xx}^{-1} A S_{xx}^{-1} x_\ell$. The proof shows that the ``kernel'' $C_{i \ell}$ varies with $n$ in such a way that $\hat \theta$ is asymptotically normal whether or not $\hat \theta$ is a degenerate U-statistic (i.e., whether or not $\beta$ is zero).
One representation of the variance appearing in \ensuremath{^{\textrm{th}}}\ref{thm3} is
\begin{align}
\mathbb{V}[\hat \theta] &= \sum_{i=1}^n \left( 2\tilde x_i'\beta - \check x_i'\beta \right)^2 \sigma_i^2 + 2 \sum_{i=1}^n \sum_{\ell \neq i} C_{i\ell}^2 \sigma_i^2 \sigma_\ell^2.
\end{align}
Note that this variance is bounded from below by $\min_i \sigma_i^2 \sum_{i=1}^n ( 2\tilde x_i'\beta)^2 + (\check x_i'\beta)^2$ since $\sum_{i=1}^n \tilde x_i'\beta \check x_i'\beta = 0$. Therefore $(i)$ will be satisfied whenever $\max_i \left( (\tilde x_i'\beta)^2 + (\check x_i'\beta)^2 \right)$ is not too large compared to $\sum_{i=1}^n (\tilde x_i'\beta)^2 + (\check x_i'\beta)^2$. As in \ensuremath{^{\textrm{th}}}\ref{thm2}, $(i)$ is implied by a variety of primitive conditions that limit how far a maximum is from the average, but since $(i)$ involves a one dimensional function of $x_i$ it can also be satisfied when $r$ is large. A particularly simple case where $(i)$ is satisfied is when $\beta=0$; further cases are discussed in Section \ref{sec:verify}.
\begin{rem}\ensuremath{^{\textrm{th}}}\label{rem:testbig}
\ensuremath{^{\textrm{th}}}\ref{thm3} can be used to test a large system of linear restrictions of the form $H_0: R\beta=0$ where $r \rightarrow \infty$ is the rank of $R \in \mathbb{R}^{r \times k}$. Under this null hypothesis, choosing $A=\frac{1}{r} R' (R S_{xx}^{-1} R')^{-1} R$ implies ${\mathbb{V}}[\hat \theta]^{-1/2} \hat \theta \xrightarrow{d} \mathcal{N}(0,1)$ since all the non-zero eigenvalues of $\tilde A$ are equal to $\frac{1}{r}$. The existing literature allows for either heteroscedastic errors and moderately few regressors \citep[][$k^3/n \rightarrow 0$]{donald2003empirical} or homoscedastic errors and many regressors \citep[][$k/n \le c <1$]{anatolyev2012inference}. When coupled with the estimator of ${\mathbb{V}}[\hat \theta]$ presented in the next subsection, this result enables tests with heteroscedastic errors and many regressors.
\end{rem}
\begin{rem}
\ensuremath{^{\textrm{th}}}\ref{thm3} extends some common results in the literature on many and many weak instruments \cite[see, e.g.,][]{chao2012asymptotic} where the estimators are asymptotically equivalent to quadratic forms. The structure of that setting is such that $\tilde A = I_r/r$ and $r \rightarrow \infty$, in which case condition $(ii)$ of \ensuremath{^{\textrm{th}}}\ref{thm3} is automatically satisfied.
\end{rem}
\subsection{Variance Estimation}
\label{sec:variance_estimation}
In order to conduct inference based on the normal approximation in \ensuremath{^{\textrm{th}}}\ref{thm3} we now propose an estimator of $\mathbb{V}[\hat \theta]$. The U-statistic representation of $\hat \theta$ in \eqref{eq:Ustat} implies that the variance of $\hat \theta$ is
\begin{align}
\mathbb{V}[\hat \theta] = 4\sum_{i=1}^n \left(\sum_{\ell \neq i} C_{i\ell} x_\ell'\beta \right)^2 \sigma_i^2 + 2\sum_{i=1}^n \sum_{\ell \neq i} C_{i\ell}^2 \sigma_i^2 \sigma_\ell^2.
\end{align}
Naively replacing $\{x_i'\beta,\sigma_i^2\}_{i=1}^n$ with $\{y_i,\hat \sigma_i^2\}_{i=1}^n$ in the above formula to form a plug-in estimator of $\mathbb{V}[\hat \theta]$ will, in general, lead to invalid inferences as $\hat \sigma_i^2 \hat \sigma_\ell^2$ is a biased estimator of $\sigma_i^2 \sigma_\ell^2$. For this reason, we consider estimators of the error variances that rely on leaving out more than one observation. Since this approach places additional restrictions on the design, Appendix \ref{sec:cons} describes a simple adjustment which leads to conservative inference in settings where these restrictions do not hold.
\noindent \textbf{\emph{Sample Splitting}} Our specific proposal is an estimator that exploits two independent unbiased estimators of $x_i'\beta$ that are also independent of $y_i$. We denote these estimators $\widehat{x_i'\beta}_{-i,s} = \sum_{\ell \neq i}^n P_{i\ell,s} y_\ell$ for $s = 1,2$, where $P_{i\ell,s}$ does not (functionally) depend on the $\{y_i\}_{i=1}^n$. To ensure independence between $\widehat{x_i'\beta}_{-i,1}$ and $\widehat{x_i'\beta}_{-i,2}$, we require that $P_{i\ell,1} P_{i\ell,2} =0$ for all $\ell$. Employing these split sample estimators, we create a new set of unbiased estimators for $\sigma_i^2$:
\begin{align}
\tilde \sigma_{i}^2 = \left(y_i - \widehat{x_i'\beta}_{-i,1}\right)\left(y_i - \widehat{x_i'\beta}_{-i,2}\right)
\quad \text{and} \quad
\hat \sigma_{i,-\ell}^2 = \begin{cases}
y_i(y_i - \widehat{x_i'\beta}_{-i,1}), & \text{if } P_{i\ell,1}= 0, \\
y_i(y_i - \widehat{x_i'\beta}_{-i,2}), & \text{if } P_{i\ell,1}\neq 0,
\end{cases}
\end{align}
where $\hat \sigma_{i,-\ell}^2 $ is independent of $y_\ell$ and $\tilde \sigma_{i}^2$ is a cross-fit estimator of the form considered in \cite{newey2018cross}. These cross-fit estimators can be used to construct an estimator of $\sigma_i^2 \sigma_\ell^2$ that, under certain design conditions, will be unbiased.
Letting $P_{im,-\ell} = P_{im,1}1_{\{P_{i\ell,1}= 0\}} + P_{im,2}1_{\{P_{i\ell,1} \neq 0\}}$ denote the weight observation $m$ receives in $\hat \sigma_{i,-\ell}^2$ and $\tilde C_{i\ell} = C_{i\ell}^2 + 2\sum_{m =1 }^n C_{mi} C_{m\ell} (P_{mi,1} P_{m\ell,2} + P_{mi,2} P_{m\ell,1})$, we define
\begin{align}
\widehat{\sigma_i^2 \sigma_\ell^2} &=
\begin{cases}
\hat \sigma_{i,-\ell}^2 \cdot \hat \sigma_{\ell,-i}^2, & \text{if } P_{im,-\ell}P_{\ell m,-i}=0 \text{ for all } m, \\
\tilde \sigma_{i}^2 \cdot \hat \sigma_{\ell,-i}^2, & \text{else if } P_{i\ell,1} + P_{i\ell,2}=0, \\
\hat \sigma_{i,-\ell}^2 \cdot \tilde \sigma_{\ell}^2, & \text{else if } P_{\ell i,1} + P_{\ell i,2}=0, \\
\hat \sigma_{i,-\ell}^2 \cdot (y_\ell - \bar y)^2 \cdot 1_{\{\tilde C_{i\ell} < 0\}}, & \text{otherwise.}
\end{cases}
\end{align}
The first three cases in the above definition correspond respectively to pairs where (i) $\hat \sigma_{i,-\ell}^2$ and $\hat \sigma_{\ell,-i}^2$ are independent, (ii) $\widehat{x_i'\beta}_{-i,1}$ and $\widehat{x_i'\beta}_{-i,2}$ are independent of $y_\ell$, and (iii) $\widehat{x_\ell'\beta}_{-\ell,1}$ and $\widehat{x_\ell'\beta}_{-\ell,2}$ are independent of $y_i$. When any of these three cases apply, we obtain an unbiased estimator of $\sigma_i^2 \sigma_\ell^2$. For the remaining set of pairs $\mathcal{B} = \{(i,\ell) : P_{im,-\ell}P_{\ell m,-i}\neq 0 \text{ for some } m, \ P_{i\ell,1} + P_{i\ell,2} \neq 0, \ P_{\ell i,1} + P_{\ell i,2} \neq 0 \}$ that comprise the fourth case we rely on an unconditional variance estimator which leads to a biased estimator of $\sigma_i^2 \sigma_\ell^2$ and conservative inference.
\noindent \textbf{\emph{Design Requirements}}
Constructing the above split sample estimators places additional requirements on the design matrix $S_{xx}$. We briefly discuss these requirements in the context of \ensuremath{^{\textrm{th}}}\ref{ex:ANOVA,ex:RC,ex:AKM}. In the ANOVA setup of \ensuremath{^{\textrm{th}}}\ref{ex:ANOVA}, leave-one-out estimation requires a minimum group size of two, whereas existence of $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$ requires groups sizes of at least three. Conservative inference can be avoided (i.e., the set $\mathcal{B}$ will be empty) when the minimum group size is at least four. In the random coefficients model of \ensuremath{^{\textrm{th}}}\ref{ex:RC}, minimum group sizes of three and five are sufficient to ensure feasibility of leave-one-out estimation and existence of $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$, respectively. Conservativeness can be avoided with a minimum group size of seven.
In the first differenced two-way fixed effects model of \ensuremath{^{\textrm{th}}}\ref{ex:AKM}, the predictions $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$ are associated with particular paths in the worker-firm mobility network and independence requires that these paths be edge-disjoint. Menger's theorem \citep{menger1927allgemeinen} implies that $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$ exists if the design matrix has full rank when any two observations are dropped. Menger's theorem also implies that conservativeness can be avoided if the design matrix has full rank when any three observations are dropped. In our application, we use Dijkstra's algorithm to find the paths that generate $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$ (see Appendix \ref{sec:samplesplitalg} for further details).
\noindent \textbf{\emph{Consistency}}
The following lemma shows that $\widehat{\sigma_i^2 \sigma_\ell^2} $ can be utilized to construct an estimator of $\mathbb{V}[\hat \theta]$ that delivers consistent inference when sufficiently few pairs fall into $\mathcal{B}$ and provides conservative inference otherwise.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:varCS}
For $s =1,2$, suppose that $\widehat{x_i'\beta}_{-i,s}$ satisfies (unbiasedness) $\sum_{\ell \neq i}^n P_{i\ell,s} x_\ell'\beta=x_i'\beta$, (sample splitting) $P_{i\ell,1} P_{i\ell,2} =0$ for all $\ell$, and (projection property) $\lambda_{\max}(P_sP_s') = O(1)$ where $P_s = (P_{i\ell,s})_{i,\ell}$ is the hat-matrix corresponding to $\widehat{x_i'\beta}_{-i,s}$. Let
\begin{align}
\hat{\mathbb{V}}[\hat \theta] &= 4\sum_{i=1}^n \left( \sum_{\ell \neq i} C_{i\ell} y_\ell \right)^2 \tilde \sigma_{i}^2 - 2 \sum_{i=1}^n \sum_{\ell \neq i} \tilde C_{i\ell} \widehat{\sigma_i^2 \sigma_\ell^2}.
\end{align}
\begin{enumerate}
\item If the conditions of \ensuremath{^{\textrm{th}}}\ref{thm3} hold and
$\abs{\mathcal{B}} = O(1)$, then
$\frac{\hat \theta - \theta}{\hat {\mathbb{V}}[\hat \theta]^{1/2}} \xrightarrow{d} \mathcal{N}(0,1).$
\item If the conditions of \ensuremath{^{\textrm{th}}}\ref{thm3} hold, then $\liminf_{n \rightarrow \infty}\mathbb{P}\left( \theta \in \left[ \hat \theta \pm z_{\alpha} \hat{\mathbb{V}}[\hat \theta]^{1/2} \right] \right) \ge 1-\alpha$ where $z_{\alpha}^2$ denotes the $(1-\alpha)$'th quantile of a central $\chi^2_{1}$ random variable.
\end{enumerate}
\end{lem}
In the formula for $\hat{\mathbb{V}}[\hat \theta]$, the first term can be seen as a plug-in estimator and standard results for quartic forms imply that the expectation of this term is ${\mathbb{V}}[\hat \theta] + 2 \sum_{i=1}^n \sum_{\ell \neq i} \tilde C_{i\ell} {\sigma_i^2 \sigma_\ell^2}$. Hence, the second term is a bias correction which completely removes the bias when $\mathcal{B} = \emptyset$ and leaves a positive bias otherwise. In Appendix \ref{sec:cons} we establish validity of an adjustment to $\hat{\mathbb{V}}[\hat \theta]$ that utilizes an upward biased unconditional variance estimator for observations where it is not possible to construct $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$.
\begin{rem}\ensuremath{^{\textrm{th}}}\label{rem:std}
The purpose of the condition $\abs{\mathcal{B}} = O(1)$ in the above lemma is to ensure that the bias of $\hat{\mathbb{V}}[\hat \theta]$ grows small with the sample size. Because the bias of $\hat{\mathbb{V}}[\hat \theta]$ is non-negative, inference based on $\hat{\mathbb{V}}[\hat \theta]$ remains valid even when this condition fails, as stated in the second part of \ensuremath{^{\textrm{th}}}\ref{lem:varCS}. In practice, it may be useful for researchers to calculate the fraction of pairs that belong to $\mathcal{B}$ to gauge the extent to which inference might be conservative. Similarly, it may be useful to compute the share of observations where it is not possible to construct $\{\widehat{x_i'\beta}_{-i,s}\}_{s=1,2}$ to investigate whether upward bias in the standard error could lead to power concerns.
\end{rem}
\section{Weakly Identified Quadratic Forms of Growing Rank}\label{sec:weak}
In some settings where $r$ grows with the sample size, condition (ii) of \ensuremath{^{\textrm{th}}}\ref{thm3} may not apply. For example in two-way fixed effects models, it is possible that ``bottlenecks'' arise in the mobility network that lead the largest eigenvalues to dominate the others.
This section provides a theorem which covers the case where some of the squared eigenvalues $\lambda_1^2,\dots,\lambda_r^2$ are large relative to their sum $\sum_{\ell=1}^r \lambda_\ell^2$. To motivate this assumption, note that each eigenvalue of $\tilde A$ measures how strongly $\theta$ depends on a particular linear combination of the elements of $\beta$ relative to the difficulty of estimating that combination (as summarized by $S_{xx}^{-1}$). From \ensuremath{^{\textrm{th}}}\ref{lem:cons}, $\text{trace}(\tilde A^2) = \sum_{\ell=1}^r \lambda_\ell^2$ governs the total variability in $\hat \theta$. Therefore, \ensuremath{^{\textrm{th}}}\ref{thm4} covers the case where $\theta$ depends strongly on a few linear combinations of $\beta$ that are imprecisely estimated relative to the overall sampling uncertainty in $\hat \theta$. The following assumption formalizes this setting.
\begin{assumption}\ensuremath{^{\textrm{th}}}\label{ass:eig}
There exist a $c >0$ and a known and fixed $q \in \{1,\dots,r-1\}$ such that
\begin{align}
\frac{ \lambda_{q+1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} = o(1) \quad \text{and} \quad \frac{\lambda_{q}^2}{\sum_{\ell=1}^r \lambda_\ell^2} \ge c \quad \text{for all } n.
\end{align}
\end{assumption}
\ensuremath{^{\textrm{th}}}\ref{ass:eig} defines $q$ as the number of squared eigenvalues that are large relative to their sum. Equivalently, $q$ indexes the number of nuisance parameters in $b$ that are \emph{weakly identified} relative to their influence on $\theta$ and the uncertainty in $\hat \theta$. The assumption that $q$ is known is motivated by our discussion of Examples \ref{ex:R2}--\ref{ex:AKM} in Section \ref{sec:verify} and the theoretical literature on weak identification, which typically makes an ex-ante distinction between strongly and weakly identified parameters \citep[e.g.,][]{andrews2012estimation}. In Section \ref{sec:getq} we offer some guidance on choosing $q$ in settings where it is unknown.
\subsection{Limit Distribution}
Given knowledge of $q$, we can split $\hat \theta$ into a known function of $\hat {\mathsf{b}}_q= (\hat b_1,\dots,\hat b_{q})'$ and $\hat \theta_q$ where $\hat b_1,\dots,\hat b_{q}$ are OLS estimators of the weakly identified nuisance parameters:
\begin{align}
\hat {\mathsf{b}}_q & = \sum_{i=1}^n \mathsf{w}_{iq} y_i, & \mathsf{w}_{iq} &= ( w_{i1},\dots, w_{iq})', \\
\hat \theta_q &= \hat \theta - \sum_{\ell=1}^{q} \lambda_\ell (\hat b_\ell^2 - \hat{\mathbb{V}}[\hat b_\ell] ), & \hat{\mathbb{V}}[\hat b] &= \sum_{i=1}^n w_{i} w_{i}' \hat \sigma_{i}^2.
\end{align}
The main difficulty in proving the following \ensuremath{^{\textrm{th}}}\nameref{thm4} is to show that the joint distribution of $(\hat{\mathsf{b}}_q',\hat \theta_q)'$ is normal, which we do using the same variation of Stein's method that was employed for \ensuremath{^{\textrm{th}}}\ref{thm3}. The high-level conditions involve $\tilde x_{iq}$ and $\check x_{iq}$ which are the parts of $\tilde x_{i}$ and $\check x_{i}$ that pertain to $\hat \theta_q$ and are defined in the proof of \ensuremath{^{\textrm{th}}}\ref{thm4}.
\begin{thm}\ensuremath{^{\textrm{th}}}\label{thm4}
If $\max_i \mathsf{w}_{iq}'\mathsf{w}_{iq} = o(1)$, $\mathbb{V}[\hat \theta_q]^{-1} \max_i \left( (\tilde x_{iq}'\beta)^2 + (\check x_{iq}'\beta)^2 \right) = o(1)$, and \ensuremath{^{\textrm{th}}}\ref{ass:reg,ass:eig} hold, then
\begin{enumerate}
\item $\mathbb{V}[(\hat{\mathsf{b}}_q',\hat \theta_q)']^{-1/2}
\left(
(\hat{\mathsf{b}}_q',\hat \theta_q)'
- \mathbb{E}[(\hat{\mathsf{b}}_q',\hat \theta_q)']
\right)
\xrightarrow{d} \mathcal{N}\left(0, I_{q+1} \right)$
\item \label{eq:thm1} $\hat \theta = \sum_{\ell = 1}^{q} \lambda_\ell\left(\hat b_\ell^2-\mathbb{V}[\hat b_{\ell}]\right) + \hat \theta_q + o_p(\mathbb{V}[\hat \theta]^{1/2})$
\end{enumerate}
\end{thm}
\ensuremath{^{\textrm{th}}}\ref{thm4} provides an approximation to $\hat \theta$ in terms of a quadratic function of $q$ asymptotically normal random variables and a linear function of one asymptotically normal random variable. Here, the non-centralities $\mathbb{E}[\hat {\mathsf{b}}_q] = (b_1,\dots,b_q)'$ serve as nuisance parameters that influence both $\theta$ and the shape of the limiting distribution of $\hat \theta - \theta$. The next section proposes an approach to dealing with these nuisance parameters that provides asymptotically valid inference on $\theta$ for any value of $q$.
\subsection{Variance Estimation}
\label{sec:varq}
In \ensuremath{^{\textrm{th}}}\ref{thm4} the relevant variance is $\varSigma_q := \mathbb{V}[(\hat{\mathsf{b}}_q',\hat \theta_q)']$,
\begin{align}
\varSigma_q &= \sum_{i=1}^n \begin{bmatrix}
\mathsf{w}_{iq} \mathsf{w}_{iq}' \sigma_i^2 & 2\mathsf{w}_{iq}\left(\sum_{\ell \neq i} C_{i\ell q} x_\ell'\beta \right) \sigma_i^2 \\
2 \mathsf{w}_{iq}'\left(\sum_{\ell \neq i} C_{i\ell q} x_\ell'\beta \right) \sigma_i^2 & 4 \left(\sum_{\ell \neq i} C_{i\ell q} x_\ell'\beta \right)^2 \sigma_i^2 + 2 \sum_{\ell \neq i} C_{i\ell q}^2 \sigma_i^2 \sigma_\ell^2
\end{bmatrix},
\end{align}
where $C_{i\ell q}$ is defined in the proof of \ensuremath{^{\textrm{th}}}\ref{thm4}. Our estimator of this variance reuses the split sample estimators introduced for \ensuremath{^{\textrm{th}}}\ref{thm3}:
\begin{align}
\hat\varSigma_q &= \sum_{i=1}^n \begin{bmatrix}
\mathsf{w}_{iq} \mathsf{w}_{iq}' \hat \sigma_i^2 & 2 \mathsf{w}_{iq}\left(\sum_{\ell \neq i} C_{i\ell q} y_\ell \right) \tilde \sigma_i^2 \\
2 \mathsf{w}_{iq}'\left(\sum_{\ell \neq i} C_{i\ell q} y_\ell \right) \tilde \sigma_i^2 & 4 \left(\sum_{\ell \neq i} C_{i\ell q} y_\ell \right)^2 \tilde \sigma_i^2 - 2 \sum_{\ell \neq i} \tilde C_{i\ell q}^2 \widetilde{\sigma_i^2 \sigma_\ell^2}
\end{bmatrix}
\end{align}
where $\tilde C_{i\ell q}$ and $\widetilde{\sigma_i^2 \sigma_\ell^2}$ are defined in the proof of the next lemma which shows consistency of this variance estimator.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:4}
For $s =1,2$, suppose that $\widehat{x_i'\beta}_{-i,s}$ satisfies $\sum_{\ell \neq i}^n P_{i\ell,s} x_\ell'\beta=x_i'\beta$, $P_{i\ell,1} P_{i\ell,2} =0$ for all $\ell$, and $\lambda_{\max}(P_sP_s') = O(1)$. If the conditions of \ensuremath{^{\textrm{th}}}\ref{thm4} hold and $\abs{\mathcal{B}} = O(1)$, then $\varSigma_q^{-1}
\hat\varSigma_q \xrightarrow{p} I_{q+1}.$
\end{lem}
\begin{rem}
As in the case of variance estimation for \ensuremath{^{\textrm{th}}}\ref{thm3}, it may be that the design does not allow for construction of the predictions $\widehat{x_i'\beta}_{-i,1}$ and $\widehat{x_i'\beta}_{-i,2}$ used in $\hat\varSigma_q$. For such cases, Appendix \ref{sec:cons} proposes an adjustment to $\hat\varSigma_q$ which has a positive definite bias and therefore leads to valid (but conservative) inference when coupled with the inference method discussed in the next section.
\end{rem}
\section{Inference with Nuisance Parameters}\label{sec:variance}
In this section, we develop a two-sided confidence interval for $\theta$ that delivers asymptotic size control conditional on a choice of $q$. Our proposal involves inverting a minimum distance statistic in $\hat{\mathsf{b}}_q$ and $\hat \theta_q$, which \ensuremath{^{\textrm{th}}}\ref{thm4} implies are jointly normally distributed. To avoid the conservatism associated with standard projection methods \cite[e.g.,][]{dufour2001finite}, we seek to adjust the critical value downwards to deliver size control on $\theta$ rather than $\mathbb{E}[(\hat{\mathsf{b}}_q', \hat \theta_q)']$. However, unlike in standard projection problems (e.g., the problem of subvector inference), $\theta$ is a nonlinear function of $\mathbb{E}[\hat{\mathsf{b}}_q]$. To accommodate this complication, we use a critical value proposed by \cite{andrews2016geometric} that depends on the curvature of the problem.
\subsection{Inference With Known $q$}
The confidence interval we consider is based on inversion of a minimum-distance statistic for $(\hat{\mathsf{b}}_q',\hat \theta_q)'$ using the critical value proposed in \cite{andrews2016geometric}. For a specified level of confidence, $1-\alpha$, we consider the interval
\begin{align}
\hat C_{\alpha,q}^\theta &= \left[ \min_{(\dot b_1,\dots,\dot b_q,\dot \theta_q)'\in \hat{\mathsf{E}}_{\alpha,q}} \sum_{\ell=1}^{q} \lambda_\ell \dot b_\ell^2 + \dot \theta_q , \max_{(\dot b_1,\dots,\dot b_q,\dot \theta_q)'\in \hat{\mathsf{E}}_{\alpha,q}} \sum_{\ell=1}^{q} \lambda_\ell \dot b_\ell^2 + \dot \theta_q \right]
\shortintertext{where}
\hat{\mathsf{E}}_{\alpha,q} &= \left\{ (\mathsf{b}_q',\theta_q)' \in \mathbb{R}^{q+1} : \begin{pmatrix}
\hat{\mathsf{b}}_q - \mathsf{b}_q \\ \hat \theta_q - \theta_q
\end{pmatrix}'\hat \varSigma_q ^{-1} \begin{pmatrix}
\hat{\mathsf{b}}_q - \mathsf{b}_q \\ \hat \theta_q - \theta_q
\end{pmatrix} \le z_{\alpha,\hat \kappa_q}^2\right\}.
\end{align}
The critical value function, $z_{\alpha,\kappa}$, depends on the maximal curvature, $\kappa$, of a certain manifold (exact definitions of $z_{\alpha,\kappa}$ and $\kappa$ are given in Appendix \ref{app:conf}). Heuristically, $\kappa$ can be thought of as summarizing the influence of the nuisance parameter $\mathbb{E}[\hat {\mathsf{b}}_q]$ on the shape of $\hat \theta$'s limiting distribution. Accordingly, $z_{\alpha}^2 := z_{\alpha,0}^2$ is equal to the $(1-\alpha)$'th quantile of a central $\chi^2_{1}$ random variable. As $\kappa \rightarrow \infty$, $z_{\alpha,\kappa}^2$ approaches the $(1-\alpha)$'th quantile of a central $\chi^2_{q+1}$ random variable. This upper limit on $z_{\alpha,\kappa}$ is used in the projection method in its classical form as popularized in econometrics by \cite{dufour2001finite}, while the lower limit $z_{\alpha}$ would yield size control if $\theta$ were linear in $\mathbb{E}[\hat{\mathsf{b}}_q]$.
When $q=0$, the maximal curvature is zero and $\hat C_{0}^\theta$ simplifies to $[\hat \theta \pm z_{\alpha} \hat{\mathbb{V}}[\hat \theta]^{1/2} ]$.
When $q=1$, the maximal curvature is $\hat \kappa_1 = \frac{2\abs{\lambda_1}\hat{\mathbb{V}}[\hat b_{1}]}{\hat{\mathbb{V}}[\hat \theta_1]^{1/2}(1-\hat \rho^2)^{1/2}}$ where $\hat \rho$ is the estimated correlation between $\hat b_{1}$ and $\hat \theta_1$. This curvature measure is intimately related to eigenvalue ratios previously introduced, as $\hat \kappa_1^2$ is approximately equal to $\frac{2\lambda_1^2}{\sum_{\ell=2}^{r} \lambda_\ell^2}$ when the error terms are homoscedastic and $\beta=0$. A closed form expression for the $q=1$ confidence interval is provided in Appendix \ref{app:conf}. When $q>1$, inference relies on solving two quadratic optimization problems that involve $q+1$ unknowns, which can be achieved reliably using standard quadratic programming routines.
The following lemma shows that a consistent variance estimator as proposed in \ensuremath{^{\textrm{th}}}\ref{lem:4} suffices for asymptotic validity under the conditions of \ensuremath{^{\textrm{th}}}\ref{thm4} and Appendix \ref{sec:cons} establishes validity when only a conservative variance estimator is available.
\begin{lem}\ensuremath{^{\textrm{th}}}\label{lem:inf}
If $\varSigma_q ^{-1} \hat \varSigma_q \xrightarrow{p} I_{q+1}$ and the conditions of \ensuremath{^{\textrm{th}}}\ref{thm4} hold, then
\begin{align}
\liminf_{n \rightarrow \infty} \mathbb{P}\left( \theta \in \hat C_{\alpha,q}^\theta \right) \ge 1-\alpha.
\end{align}
\end{lem}
The confidence interval studied in \ensuremath{^{\textrm{th}}}\ref{lem:inf} constructs a $q+1$ dimensional ellipsoid $\hat{\mathsf{E}}_{\alpha,q}$ and maps it through the quadratic function $(\dot b_1,\dots,\dot b_q,\dot \theta_q) \mapsto \sum_{\ell=1}^{q} \lambda_\ell \dot b_\ell^2 + \dot \theta_q$. This approach ensures uniform coverage over any possible values of the nuisance parameters $b_1,\dots,b_q$ which are imprecisely estimated relative to overall sampling uncertainty in $\hat \theta$.
\begin{rem}
An alternative to \ensuremath{^{\textrm{th}}}\ref{lem:inf} is to conduct inference using a first-order Taylor expansion of $\sum_{\ell=1}^{q} \lambda_\ell \hat b_\ell^2 + \hat \theta_q$. This so-called ``Delta method'' approach is asymptotically equivalent to using the confidence interval $[\hat \theta \pm z_{\alpha} \hat{\mathbb{V}}[\hat \theta]^{1/2} ]$ studied in Section \ref{sec:dist}. However, the Delta method is not uniformly valid in the presence of nuisance parameters as approximate linearity can fail when $\min_{\ell \le q} b_\ell^2 = O(1)$. Section \ref{sec:verify} introduces a stochastic block model with $q=1$ and characterizes $b_1^2$ as the squared difference in average firm effects across two blocks multiplied by the number of between block movers. Thus the Delta method will potentially undercover unless there are strong systematic differences between the two blocks.
\end{rem}
\subsection{Choosing $q$}\label{sec:getq}
The preceding discussion of inference considered a setting where the number of weakly identified parameters was known in advance. In some applications, it may not be clear ex ante what value $q$ takes. In such situations researchers may wish to report confidence intervals for two consecutive values of $q$ (or their union).
This heuristic serves to minimize the influence of the specific value of $q$ picked, and both our simulations and empirical application suggest that $\hat C_{\alpha,q}^\theta$ barely varies with $q$ when $\frac{ \lambda_{ q+1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} < \frac{1}{10}$. Consequently, little power is sacrificed by taking the union.
This observation also suggests a heuristic threshold for choosing $q$; namely, to let $q$ be such that $\frac{ \lambda_{q}^2}{\sum_{\ell=1}^r \lambda_\ell^2} \ge \frac{1}{10}$ and $\frac{ \lambda_{ q+1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} < \frac{1}{10}$, with $q=0$ when $\frac{ \lambda_{1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} < \frac{1}{10}$.
A similar threshold rule can be motivated under a slight strengthening of \thref{ass:eig} which allows one to learn $q$ from the data.
\begin{customass}{2$^\prime$}\ensuremath{^{\textrm{th}}}\label{ass:eig2}
There exist a $c >0$, an $\epsilon >0$, and a fixed $q \in \{1,\dots,r-1\}$ such that
\begin{align}
\frac{ \lambda_{q+1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} = O(r^{-\varepsilon}) \quad \text{and} \quad \frac{\lambda_{q}^2}{\sum_{\ell=1}^r \lambda_\ell^2} \ge c \quad \text{for all } n.
\end{align}
\end{customass}
A threshold based choice of $q$ is the unique $\hat q$ for which
\begin{align}
\frac{ \lambda_{\hat q+1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} < c_r
\quad \text{and} \quad
\frac{\lambda_{\hat q}^2}{\sum_{\ell=1}^r \lambda_\ell^2} \ge c_r
\quad \text{for some } c_r \rightarrow 0,
\end{align}
with $\hat q=0$ when $\frac{ \lambda_{1}^2}{\sum_{\ell=1}^r \lambda_\ell^2} < c_r$.
Under Assumption 2$^\prime$, $\hat q = q$ in sufficiently large samples provided that $c_r$ is chosen so that $c_r r^{\varepsilon} \rightarrow \infty$. This condition is satisfied when $c_r$ shrinks slowly to zero, e.g., when $c_r \propto 1/\log(r)$.
\section{Verifying Conditions}\label{sec:verify}
We now revisit the examples of Section 2 and verify the conditions required to apply our theoretical results. Appendix \ref{app:sufficiency} provides further details on these calculations.
\begin{customthm}{1}(Coefficient of determination, continued)
Recall that $\theta = \sigma_{X\beta}^2 = \beta'A\beta$ where $A = \frac{1}{n}\sum_{i=1}^{n} (x_i-\bar x)(x_i-\bar x)'$ and $\tilde A = \frac{1}{n}(I_k - n S_{xx}^{-1/2} \bar x \bar x' S_{xx}^{-1/2} )$. Suppose \ensuremath{^{\textrm{th}}}\ref{ass:reg} holds.
\noindent \textbf{\emph{Consistency}} Consistency follows from \ensuremath{^{\textrm{th}}}\ref{lem:cons} since $\lambda_\ell = \frac{1}{n}$ for $\ell = 1,\dots,r$ where $r = \text{dim}(x_i)-1$. Thus $\text{trace}(\tilde A^2)=r/n^2 \le 1/n= o(1)$.
\noindent \textbf{\emph{Limit Distribution}} If $\text{dim}(x_i)$ is fixed, then $w_i'w_i = P_{ii} - \frac{1}{n}$ and \ensuremath{^{\textrm{th}}}\ref{thm2} applies under the standard ``textbook'' condition that $\max_{i} P_{ii} = o(1)$. If $\text{dim}(x_i) \rightarrow \infty$, then \ensuremath{^{\textrm{th}}}\ref{thm3} applies if $\mathbb{V}[\hat \theta]^{-1} \max_i (\check x_i'\beta )^2 = o(1)$ which follows if, e.g., $\max_{i} \frac{1}{\sqrt{r}}\sum_{\ell=1}^n \abs{M_{i\ell}} = o(1)$ where $M_{i\ell} = \mathbf{1}_{\{i = \ell\}} - x_i' S_{xx}^{-1} x_\ell$ (this condition holds in the next two examples). Equality among all eigenvalues excludes the weak identification setting of \ensuremath{^{\textrm{th}}}\ref{thm4}.
\noindent \textbf{\emph{Unbounded Mean Function}} Inspection of the proofs reveal that \ensuremath{^{\textrm{th}}}\ref{ass:reg}(iii), $\max_i (x_i'\beta)^2 = O(1)$, can be dropped if the above conditions are strengthened to $\max_{i,\ell} P_{ii} (x_\ell'\beta)^2 = o(1)$ when $\text{dim}(x_i)$ is fixed or $\max_{i,j} \frac{\abs{x_j'\beta}\left( 1+\sum_{\ell=1}^n \abs{M_{i\ell}}\right)}{\sqrt{r}} = o(1)$ when $\text{dim}(x_i) \rightarrow \infty$.
\end{customthm}
\begin{customthm}{2}(Analysis of covariance, continued)
Recall that $\theta = \sigma_{\alpha}^{2}=\frac{1}{n}\sum_{{g}=1}^{N} T_{g} \left(\alpha_{g}-\bar{\alpha}\right)^{2}$ where $y_{gt}= \alpha_g + x_{gt}'\delta + \varepsilon_{gt}$, $g$ index the $N$ groups, and $T_g$ is group size.
\noindent \textbf{\emph{No Common Regressors}} This is a special case of the previous example with $r=N-1$, $P_{ii} = T_{g(i)}^{-1}$ and $\check x_i =0$. Assumption 1(ii),(iii) requires $T_g \ge 2$ and $\max_g \alpha_g^2 = O(1)$. \ensuremath{^{\textrm{th}}}\ref{thm2} applies if the number of groups is fixed and $\min_g T_g \rightarrow \infty$, while \ensuremath{^{\textrm{th}}}\ref{thm3} applies if the number of groups is large. \ensuremath{^{\textrm{th}}}\ref{thm4} cannot apply as all eigenvalues are equal to $\frac{1}{n}$.
\noindent \textbf{\emph{Common Regressors}} To accommodate common regressors of fixed dimension, assume $\norm{\delta}^2 + \max_{g,t} \norm{x_{gt}}^2 = O(1)$ and that $\frac{1}{n}\sum_{g=1}^N \sum_{t=1}^{T_g} (x_{gt} - \bar x_g)(x_{gt} - \bar x_g)'$ converges to a positive definite limit. This is a standard assumption in basic panel data models \cite[see, e.g.,][Chapter 10]{wooldridge2010econometric}. Allowing such common regressors does not alter the previous conclusions: \ensuremath{^{\textrm{th}}}\ref{thm2} applies if $N$ is fixed and $\min_g T_g \rightarrow \infty$ since $w_i'w_i \le P_{ii} = T_{g(i)}^{-1} + O(n^{-1})$, \ensuremath{^{\textrm{th}}}\ref{thm3} applies if $N \rightarrow \infty$ since $\sum_{\ell=1}^n \abs{M_{i\ell}} = O(1)$, and \ensuremath{^{\textrm{th}}}\ref{thm4} cannot apply since $n \lambda_\ell \in [c_1,c_2]$ for $\ell =1,\dots,r$ and some $c_2 \ge c_1 > 0$ not depending on $n$.
\noindent \textbf{\emph{Unbounded Mean Function}} All conclusions continue to hold if $\max_{g,t} \alpha_g^2 + \norm{x_{gt}}^2 = O(1)$ is replaced with $\frac{\max_{g,t} \alpha_g^2 + \norm{x_{gt}}^2}{ \max \{N,\min_g T_g\}} = o(1)$ and $\sigma_\alpha^2 + \frac{1}{n}\sum_{g=1}^N \sum_{t=1}^{T_g} \norm{x_{gt}}^2 = O(1)$.
\end{customthm}
\begin{customthm}{3}(Random coefficients, continued)
For simplicity, consider the \textit{uncentered} second moment $\theta = \frac{1}{n} \sum_{g=1}^N T_g \gamma_g^2$ where $y_{gt}= \alpha_g + z_{gt}'\gamma_g + \varepsilon_{gt}$. Suppose \ensuremath{^{\textrm{th}}}\ref{ass:reg} holds and assume that $\max_{g,t} \alpha_g + \gamma_g^2 + z_{gt}^2 =O(1)$ and $\min_g S_{zz,g} \ge c > 0$ where $S_{zz,g}= \sum_{t=1}^{T_g}(z_{gt} - \bar z_{g})^2$. Note that $\min_g S_{zz,g} > 0$ is equivalent to full rank of $S_{xx}$ and $S_{zz,g}$ indexes how precisely $\gamma_g$ can be estimated.
\noindent \textbf{\emph{Consistency}} The $N$ eigenvalues of $\tilde A$ are $\lambda_g = \frac{T_g}{n}S_{zz,g}^{-1}$ for $g=1,\dots,N$ where the group indexes are ordered so that $\lambda_1 \ge \dots \ge \lambda_N$. Consistency follows from \ensuremath{^{\textrm{th}}}\ref{lem:cons} if $\lambda_1^{-1} = n\frac{S_{zz,1}}{T_1} \rightarrow \infty$. This is automatically satisfied with many groups of bounded size.
\noindent \textbf{\emph{Limit Distribution}} If $N$ is fixed and $\min_g S_{zz,g} \rightarrow \infty$, then \ensuremath{^{\textrm{th}}}\ref{thm2} applies. If $\frac{\sqrt{N}}{T_1} S_{zz,1} \rightarrow \infty$, then \ensuremath{^{\textrm{th}}}\ref{thm3} applies. If $\frac{\sqrt{N}}{T_2} S_{zz,2} \rightarrow \infty$, $ \frac{\sqrt{N}}{T_{1}} S_{zz,{1}} = O(1)$, and $S_{zz,1} \rightarrow \infty$, then \ensuremath{^{\textrm{th}}}\ref{thm4} applies with $q=1$. In this case, $\gamma_1$ is weakly identified relative to its influence on $\theta$ and the overall variability of $\hat \theta$. This is expressed through the condition $ \frac{\sqrt{N}}{T_{1}} S_{zz,{1}} = O(1)$ where $S_{zz,{1}}$ is the identification strength of $\gamma_1$, $T_1$ provides the influence of $\gamma_1$ on $\theta$ and $1/\sqrt{N}$ indexes the variability of $\hat \theta$.
\end{customthm}
\begin{customthm}{4}(Two-way fixed effects, continued)
In this final example, we restrict attention to the first-differenced setting $\Delta y_{{g}} = \Delta f_{{g}}' \psi + \Delta \varepsilon_{{g}}$ with $T_g=2$ and a large number of firms, $J \rightarrow \infty$. Our target parameter is the variance of firm effects $\theta = \sigma_{\psi}^{2} = \frac{1}{n}\sum_{{g}=1}^{N}\sum_{t=1}^{T_{g}} \left(\psi_{j\left({g},t\right)}-\bar{\psi}\right)^{2}$ and we consider \ensuremath{^{\textrm{th}}}\ref{ass:reg} satisfied; in particular, $\max_j \abs{\psi_j} = O(1)$.
\noindent \textbf{\emph{Leverages}} The leverage $P_{gg}$ of observation $g$ is less than one if the origin and destination firms of worker $g$ are connected by a path not involving $g$. Letting $n_g$ denote the number of edges in the shortest such path, one can show that $P_{gg} \le \frac{n_g}{1+n_g}$. Therefore, if $\max_g n_g < 100$ then \ensuremath{^{\textrm{th}}}\ref{ass:reg}(ii) is satisfied with $\max P_{gg} \le .99$. In our application we find $\max_g n_g =12$, leading to a somewhat smaller bound on the maximal leverage. The same consideration implies a bound on the model in levels since $P_{i(g,t)i(g,t)} = \frac{1}{2}(1+P_{gg})$.
\noindent \textbf{\emph{Eigenvalues}} The eigenvalues of $\tilde A$ satisfy the equality
\begin{align}
\lambda_\ell =\frac{1}{n \dot \lambda_{J+1-\ell}} \qquad \text{for} \quad \ell = 1,\dots,J
\end{align}
where $\dot \lambda_1 \ge \dots \ge \dot \lambda_J$ are the non-zero eigenvalues of the matrix $E^{1/2} \mathcal{L} E^{1/2}$. $\mathcal{L}$ is the normalized Laplacian of the employer mobility network and connectedness of the network is equivalent to full rank of $S_{xx}$ (see Appendix \ref{app:sufficiency} for definitions). $E$ is a diagonal matrix of employer specific ``churn rates'', i.e., the number of moves in and out of a firm divided by the total number of employees in the firm. $E$ and $\mathcal{L}$ interact in determining the eigenvalues of $\tilde A$. In \ensuremath{^{\textrm{th}}}\ref{ex:RC}, the quantities $\{T_\ell^{-1} S_{zz,\ell}\}_{\ell=1}^N$ played a role directly analogous to the churn rates in $E$, so in this example we focus on the role of $\mathcal{L}$ by assuming that the diagonal entries of $E$ are all equal to one.
\noindent \textbf{\emph{Strongly Connected Network}} The employer mobility network is \emph{strongly connected} if $\sqrt{J} \mathcal{C} \rightarrow \infty$ where $\mathcal{C} \in (0,1]$ is Cheeger's constant for the mobility network \cite[see, e.g.,][]{mohar1989isoperimetric,jochmans2016fixed}. Intuitively, $\mathcal{C}$ measures the most severe ``bottleneck'' in the network, where a bottleneck is a set of movers that upon removal from the data splits the mobility network into two disjoint blocks. The severity of the bottleneck is governed by the number of movers removed divided by the smallest number of movers in either of the two disjoint blocks. The inequalities $\dot \lambda_J \ge 1 - \sqrt{1-\mathcal{C}^2}$ \cite[][Theorem 2.3]{chung1997spectral} and ${\lambda_1^2}/{\sum_{\ell=1}^J \lambda_\ell^2} \le 4(\sqrt{J} \dot \lambda_J)^{-2}$ imply that a strongly connected network yields $q=0$, which rules out application of \ensuremath{^{\textrm{th}}}\ref{thm4}. Furthermore, a strongly connected network is sufficient (but not necessary) for consistency of $\hat \theta$ as $\sum_{\ell=1}^J \lambda_\ell^2 \le \frac{J}{n} ( \sqrt{n} \dot \lambda_J)^{-2}$.
\noindent \textbf{\emph{Weakly Connected Network}} When $\sqrt{J} \mathcal{C}$ is bounded, the network is \emph{weakly connected} and can contain a sufficiently severe bottleneck that a linear combination of the elements of $\psi$ is estimated imprecisely relative to its influence on $\theta$ and the total uncertainty in $\hat \theta$. The weakly identified linear combination in this case is a difference in average firm effects across the two blocks on either side of the bottleneck, which contributes a $\chi^2$ term to the asymptotic distribution. Below we use a stochastic block model to further illustrate this phenomenon. Our empirical application demonstrates that weakly connected networks can appear in practically relevant settings.
\noindent \textbf{\textit{Stochastic Block Model}} Consider a stochastic block model of network formation where firms belong to one of two blocks and a set of workers switch firms, possibly by moving between blocks. Workers' mobility decisions are independent: with probability $p_b$ a worker moves between blocks and with probability $1-p_b$ she moves within block. For simplicity, we further assume that the two blocks contain equally many firms and consider a semi-sparse network where $\frac{J\log(J)}{n} + \frac{\log(J)}{np_b} \rightarrow 0$.\footnote{The semi-sparse stochastic block model is routinely employed in the statistical literature on spectral clustering, see, e.g, \cite{sarkar2015role}.} In this model the asymptotic behavior of $\hat \theta$ is governed by $p_b$: the most severe bottleneck is between the two blocks and has a Cheeger's constant proportional to $p_b$. In Appendix \ref{app:sufficiency}, we use this model to verify the high-level conditions leading to \ensuremath{^{\textrm{th}}}\ref{thm3,thm4} and show that \ensuremath{^{\textrm{th}}}\ref{thm3} applies when $\sqrt{J} p_b \rightarrow \infty$, while \ensuremath{^{\textrm{th}}}\ref{thm4} applies with $q=1$ otherwise. The argument extends to any finite number of blocks, in which case $q$ is the number of blocks minus one. Finally, we show that $\hat \theta$ is consistent even when the network is weakly connected. To establish consistency we only impose $\frac{\log(J)}{np_b} \rightarrow 0$, which requires that the number of movers across the two blocks is large.
\end{customthm}
\section{Application}\label{sec:application}
Consider again the problem of estimating variance components in two-way fixed effect models of wage determination. \cite{card2016firms} note that plug-in wage decompositions of the sort introduced by AKM typically attribute $15\%$--$25\%$ of overall wage variance to variability in firm fixed effects. Given the bias and potential sampling
variability associated with plug-in estimates, however, it has been difficult to infer whether firm effects play a differentially important role in certain markets or among particular demographic groups.
In this section, we use Italian social security records to compute leave-out estimates of
the AKM wage decomposition and contrast them with estimates based upon the plug-in estimator of \citet{abowd1999high} and the homoscedasticity-corrected estimator of \cite{andrews2008high}. We then investigate whether
the variance components that comprise the AKM decomposition differ across age groups. While it is well known that wage inequality increases with age \citep{mincer1974schooling, lemieux2006increasing}, less is known about the extent to which firm pay premia mediate this phenomenon. Standard wage posting models \citep[e.g., ][]{burdett1998wage} suggest older workers have had more time to climb (and fall off) the job ladder and to receive outside offers \citep{bagger2014tenure}, which may result in more dispersed firm wage premia. But older workers have also had more time to develop professional reputations revealing their relative productivity, which should generate a large increase in the variance of person effects \citep{gibbons1992does,gibbons2005comparative}. The tools developed in this paper allow us to formally study these hypotheses.
\subsection{Sample Construction}
The data used in our analysis come from the Veneto Worker History (VWH) file, which provides the annual earnings
and days worked associated
with each covered employment spell taking place in the Veneto region of Northeast Italy
over the years 1984-2001. The VWH data have been used in a number of recent studies \citep{card2014rent,bartolucci2018identifying,serafinelli2019good,devicienti2019collective} and are well suited to the analysis of age differences because they provide precise information on dates of birth. These data are also notable for being publicly available, making the costs of replicating our analysis unusually low.\footnote{See \url{http://www.frdb.org/page/data/scheda/inps-data-veneto-workers-histories-vwh/doc_pk/11145} for information on obtaining the VHW.}
Our baseline sample consists of workers with employment spells taking place in the years 1999 and 2001, which provides us with a three year horizon over which to measure job mobility. In Section \ref{sec:bigT_AKM} we analyze a longer unbalanced sample spanning the years 1996--2001 and find that it yields similar results. For each worker-year pair, we retain the unique employment spell yielding the highest earnings in that year. Wages in each year are defined as earnings in the selected spell divided by the spell length in days. Workers are divided into two groups of roughly equal size according to their year of birth: ``younger'' workers born in the years 1965-1983 (aged 18-34 in 1999) and ``older'' workers born in the years 1937-1964 (aged 35-64 in 1999). Further details on our processing of the VWH records is provided in Appendix \ref{app:data}.
Table \ref{table1} reports the number of person-year observations available among
workers employed by firms in the region's largest connected set, along with the largest connected
set for each age group. Workers are classified as ``movers'' if they switch firms between 1999 and 2001. Comparing the number of movers to half the number of person-year observations reveals that roughly $21\%$ of all workers are movers. The movers share rises to $26\%$ among younger workers while only $16\%$ of older workers are movers, reflecting the tendency of mobility rates to decline with age. The average number of movers per connected
firm ranges from nearly 3 in the pooled sample to roughly 2 in the thinner age-specific samples, suggesting that many firms are
associated with only a single mover.
Our leave-out estimation strategy requires that each firm effect remain estimable
after removing any single observation. The second panel of Table \ref{table1} enforces this requirement by
restricting to firms that remain connected when any
mover is dropped (see Appendix \ref{sec:pruning} for computational details).
Pruning the sample in this way drops roughly half of the firms but less than a third of the movers and eliminates roughly 30\% of all workers regardless of their mobility status. These additional restrictions raise mean wages by roughly 5\% and lower the variance of wages by 5--10\% depending on the sample.
To assess the potential influence of these sample restrictions on our estimands of interest, we construct a third sample that further requires the firm effects to remain estimable after removing any two observations.\footnote{We thank an anonymous referee for this suggestion.} This ``leave-two-out connected set'' is also of theoretical interest because it provides a setting where the requirements for consistency of the variance estimator of Lemma \ref{lem:varCS} appear to be satisfied. On average, the leave-two-out connected sets have roughly half as many firms and 20\% fewer movers than the corresponding leave-one-out sets, and the average number of movers per firm
ranges from approximately 5.6 in the sample of older workers to 4.3 in the sample of younger workers. Restricting the sample in this way further raises mean wages by 3--4\% but yields negligible changes in variance, except among the sample of older workers, which experiences a nearly 7\% \emph{increase} in variance. We investigate below the extent to which these changes in unconditional variances reflect changes in the variance of underlying firm wage effects.
\subsection{AKM Model and Design Diagnostics}
\label{sec:estimates_AKM}
Consider the following simplified version of the AKM model:
\begin{align}
\label{AKM_appli}
y_{{g} t} = \alpha_{{g}} + \psi_{j({g},t)} + \varepsilon_{{g} t}. && ({g}=1,\dots,N, \ t=1,2)
\end{align}
We fit models of this sort to the VHW data after having pre-adjusted log wages for year effects in a first step. This adjustment is obtained by estimating an augmented version of the above model by OLS that includes a dummy control for the year 2001. Hence, $y_{{g} t}$ gives the log wage in year $t$ minus a year 2001 dummy times its estimated coefficient. This two-step approach simplifies computation without compromising consistency because the year effect is estimated at a $\sqrt{N}$ rate.
The bottom of Table \ref{table1} reports for each sample the maximum leverage $(\max_i P_{ii})$ of any person-year observation (Appendix \ref{sec:Lambda} discusses the computation of these leverages). While our pruning procedure ensures $\max_i P_{ii}<1$, it is noteworthy that $\max_i P_{ii}$ is still quite close to one, indicating that certain person-year observations remain influential on the parameter estimates. This finding highlights the inadequacy of asymptotic approximations that require the dimensionality of regressors to grow slower than the sample size, which would lead the maximum leverage to tend to zero.
The asymptotic results of Section \ref{sec:weak} emphasize the importance of not only the maximal leverage, but the number and severity of any bottlenecks in the mobility network. Figure \ref{fig:net} illustrates the leave-two-out connected set for older workers. Each firm is depicted as a dot, with the size of the dot proportional to the total number of workers employed at the firm over the years 1999 and 2001. Dots are connected when a worker moves between the corresponding pair of firms. The figure highlights the two most severe bottlenecks in this network, which divide the firms into three distinct blocks. Each block's firms have been shaded a distinct color. The blue block consists of only five firms, four of which are quite small, which limits its influence on the asymptotic behavior of our estimator. However, the green block has 51 firms with a non-negligible employment share of $9.5\%$. \ensuremath{^{\textrm{th}}}\ref{thm4} and the discussion in Section \ref{sec:verify} therefore suggest that the bottleneck between the green and the larger red block will generate weak identification and asymptotic non-normality, predictions we explore in detail below.
\subsection{Variance Decompositions}
Table \ref{table2} reports the results of applying to our samples three estimators of the AKM variance decomposition: the naive plug-in (PI) estimator $\hat \theta_{\text{PI}}$ originally proposed by AKM, the homoscedasticity-only (HO) estimator $\hat \theta_{\text{HO}}$ of \cite{andrews2008high}, and the leave-out (KSS) estimator $\hat \theta$. The PI estimator finds that the variance of firm effects in the pooled leave-one-out connected set accounts for roughly 20\% of the total variance of wages, while among younger workers firm effect variability is found to account for $31\%$ of overall wage variance. Among older workers, variability in firm effects is estimated to account for only $16\%$ of the variance of wages in the leave-one-out connected set.
Are these age differences driven by biases attributable to estimation error? Applying the HO estimator of \cite{andrews2008high} reduces the estimated variances of firm effects by roughly $18\%$ in the age-pooled sample, $27\%$ in the sample of younger workers, and $16\%$ in the sample of older workers. However, the KSS estimator yields further, comparably sized, reductions in the estimated firm effect variance relative to the HO estimator, indicating the presence of substantial heteroscedasticity in these samples. For instance, in the pooled leave-one-out sample, the KSS estimator finds a variance of firm effects that accounts for only 13\% of the overall variance of wages, while the HO estimator finds that firm effects account for 16\% of wage variance.
Moreover, while the plug-in estimates suggested that the firm effect variance was greater among older than younger workers, the KSS estimator finds the opposite pattern. The KSS estimator also finds that the pooled variance of firm effects exceeds the corresponding variance in either age-specific sample, a sign that mean firm effects differ by age. We explore this between age group component of firm variability in greater depth below.
A potential concern with analyzing the leave-one-out connected set is that worker and firm behavior in this sample may be non-representative of the broader (just-)connected set. To assess this possibility, we also report estimates for the leave-two-out connected set. Remarkably, the KSS estimator finds negligible differences in the variance of firm effects between the leave-one-out and leave-two-out samples for both the pooled sample and the sample of younger workers. Among older workers the estimated firm effect variance falls by about 11\% in the leave-two-out sample, though we show below that this difference may be attributable to sampling variation. The broad similarity between leave-one-out and leave-two-out KSS estimates is likely attributable to the fact that trimmed firms tend to be small and therefore contribute little to the person-year weighted variance of firm effects that has been the focus of the literature.
PI estimates of person effect variances are much larger than the corresponding estimates of firm effect variance, accounting for $66\%$--$88\%$ of the total variance of wages depending on the sample. The PI estimator also finds that person effects are much more dispersed among older than younger workers, which is in accord with standard models of human capital accumulation and employer learning. The estimated ratio of older to younger person effect variances in the leave-one-out sample is roughly 2.6. Applying the HO estimator reduces the magnitude of the person effect variance among all age groups, but boosts the ratio of older to younger person effect variances to 3.2. The KSS estimator yields further downward corrections to estimated person effect variances, leading the contribution of person effect variability to range from only $50\%$ to $80\%$ of total wage variance. Proportionally, however, the variability of older workers remains stable at 3.2 times that of younger workers.
PI estimates of the covariance between worker and firm effects are negative in both age-restricted samples, though not in the pooled sample. When converted to correlations, these figures suggest there is mild negative assortative matching of workers to firms. Applying the HO estimator leads the covariances to change sign in both age-specific samples, while generating a mild increase in the estimated covariance of the pooled sample. In all three samples, however, the HO estimates indicate very small correlations between worker and firm effects. By contrast, the KSS estimator finds a rather strong positive correlation of 0.21 among younger workers, 0.27 among older workers, and 0.28 in the pooled leave-one-out sample, indicating the presence of non-trivial positive assortative matching between workers and firms. While the patterns in the leave-two-out sample are broadly similar, the KSS correlation estimate among older workers is substantially smaller in the leave-two-out than the leave-one-out sample (0.18 vs 0.27).
Finally, we examine the overall fit of the two-way fixed effects model using the coefficient of determination.
The PI estimator of $R^2$ suggests the two-way fixed effects model explains more than $95\%$ of wage variation in the pooled sample, $91\%$ in the sample of younger workers, and 97\% in the sample of older workers. The HO estimator of $R^2$ is equivalent to the adjusted $R^2$ measure of \cite{theil1961economic}. The adjusted $R^2$ indicates that the two-way fixed effects model explains roughly $90\%$ of the variance of wages in the pooled sample, which is quite close to the figures reported in \cite{card2013workplace} for the German labor market. Applying the KSS estimator yields very minor changes in estimated explanatory power relative to the HO estimates. Interestingly, a sample size weighted average of the age group specific KSS $R^2$ estimates lies slightly below the pooled KSS estimate of $R^2$, which suggests allowing firm effects to differ by age group fails to appreciably improve the model's fit. We examine this hypothesis more carefully in Section \ref{sec:sorting}.
\subsection{Multiple Time Periods and Serial Correlation}
\label{sec:bigT_AKM}
Thus far, our analysis has relied upon panels with only two time periods. Table \ref{table3} reports KSS estimates of the variance of firm effects in an unbalanced panel spanning the years 1996--2001. To analyze this longer panel, we expand our set of time varying covariates to include unrestricted year effects and a third order polynomial in age normalized to have slope zero at age 40 as discussed in \cite{card2016firms}.\footnote{Pre-adjusting for age has negligible effects on the variance decompositions reported in Table \ref{table2} but is quantitatively more important in this longer panel. Age adjustments are particularly pronounced among younger workers who generally exhibit greater wage growth and tend to move rapidly to higher paying firms.} Allowing up to six wage observations per worker yields a substantially larger estimation sample with roughly three times more person-year observations in the age-pooled leave-one-out connected set than was found in Table \ref{table1}. For older workers, who have especially low mobility rates, allowing more time periods raises the number of person-year observations in the leave-one-out connected set by a factor of roughly 5.7 and more than triples the number of firms.
While these additional observations will tend to reduce the bias in the plug-in estimator, using longer panels may present two distinct sets of complications. First, the equivalence discussed in \ensuremath{^{\textrm{th}}}\ref{rem:nocluster} no longer holds, which implies that leaving a single person-year observation out is unlikely to remove the bias in estimates of the variance of firm effects when the errors are serially correlated. Second, pooling many years of data may change the target parameter if firm or person effects ``drift'' with time. The bottom rows of Table \ref{table3} probe for the importance of serial correlation by leaving out ``clusters'' of observations -- as described in \ensuremath{^{\textrm{th}}}\ref{rem:cluster} -- defined successively as all observations within the same worker-firm ``match'' and all observations belonging to the same worker; see Appendix \ref{app:cluster} for computational details. Because worker $g$'s person effect is not estimable when leaving that worker's entire wage history out, we estimate a within-transformed specification that eliminates the person effects in a first step.
Leaving out the match yields an important reduction in the variance of firm effects relative to leaving out a single person-year observation, indicating the presence of substantial serial correlation within match. By contrast, leaving out the worker turns out to have negligible effects on the estimated variance of firm effects, suggesting that serial correlation across-matches is negligible. As expected, pooling several years of data reduces the bias of the PI estimator: the magnitude of the difference between the PI estimates of the variance of firm effects and the leave-worker-out estimates tends to be smaller than the corresponding difference between the PI and KSS estimates of the variance of firm effects reported in Table \ref{table2}.
Remarkably, the firm effect variance estimates that result from leaving out either the match or worker are nearly identical to the KSS estimates reported in Table \ref{table2} for both the age-pooled samples and the samples of younger workers, suggesting the firm effects are relatively stable over this longer horizon. Among older workers, the leave-cluster-out estimates of the variance of firm effects are higher than those reported in Table \ref{table2}, which is unsurprising given that the number of firms under consideration more than tripled in this longer panel. Reassuringly, however, Table \ref{table3} reveals that the KSS estimates of the variance of firm effects among older workers in the leave-one-out and leave-two-out connected sets are very close to one another. The general stability of the KSS estimates of firm effect variances to alternate panel lengths may be attributable to the relatively placid macroeconomic conditions present in Veneto over this period, see the discussion in \cite{devicienti2019collective}.
Our leave-cluster-out exercises suggest researchers seeking to analyze longer panels may be able to avoid biases stemming from serial correlation by simply collapsing the data to match means in a first step and then analyzing these means using the leave-one-observation-out estimator. This two-step approach should substantially reduce computational time while generating only mild efficiency losses due to equal weighting of matches. In what follows, we revert to our baseline sample with exactly two observations per worker.
\subsection{Sorting and Wage Structure}\label{sec:sorting}
The KSS estimates reported in Table \ref{table2} indicate that older workers exhibit somewhat less variable firm effects and a stronger correlation between person and firm effects than younger workers. These findings might reflect lifecycle differences in the sorting of workers to firms or differences in the structure of firm wage effects across the two age groups.
Table \ref{table4} explores the sorting channel by projecting the pooled firm effects from the leave-one-out sample onto a constant, an indicator for being an older worker, the log of firm size, and the interaction of the indicator with log firm size. Because these projection coefficients are linear combinations of the estimated firm effects, we use the KSS standard errors proposed in equation \eqref{lincom} and analyzed in \ensuremath{^{\textrm{th}}}\ref{thm2}. For comparison, we also report a naive standard error that treats the firm effect estimates as independent observations and computes the usual Eicker-White ``robust'' standard errors. In all cases, the KSS standard error is at least twice the corresponding naive standard error and in one case roughly 24 times larger. In light of the consistency results of \ensuremath{^{\textrm{th}}}\ref{thm2}, this finding suggests the standard practice of regressing firm effect estimates on observables in a second step without adjusting the standard errors for correlation across firm effects can yield highly misleading inferences.
The first column of Table \ref{table4} shows that older workers tend to work at firms with higher average firm effects. Evidently older workers do occupy the upper rungs of the job ladder. The second column shows that this sorting relationship is largely mediated by firm size. An older worker at a firm with a single employee is estimated to have a mean firm wage effect 0.16 log points lower than a younger worker at a firm of the same size, an economically insignificant difference that is also revealed to be statistically insignificant when using the KSS standard error. As firm size grows, older workers begin to enjoy somewhat larger firm wage premia. Evaluated at the median firm size of 12 workers, the predicted gap between older and younger workers rises to 0.54 log points, a gap that we can distinguish from zero at the 5\% level using the KSS standard error but is still quite modest. We conclude that the tendency of older workers to be employed at larger firms is a quantitatively important driver of the firm wage premia they enjoy.
Figure \ref{fig:reg} investigates to what extent the firm wage effects differ between age groups. Using the age-restricted leave-one-out connected sets, we obtain a pair of age group specific firm effect estimates $\{\hat \psi_j^Y, \hat \psi_j^O\}_{j\in\mathcal{J}}$ for the set $\mathcal{J}$ of 8,578 firms present in both samples (see Appendix \ref{app:testOY} for details). Figure \ref{fig:reg} plots the person-year weighted averages of $\hat \psi_j^Y$ and $\hat \psi_j^O$ within each centile bin of $\hat \psi_j^O$. A person-year weighted projection of $\hat \psi_j^Y$ onto $\hat \psi_j^O$ yields a slope of only 0.501. To correct this plug-in slope estimate for attenuation bias, we multiply the unadjusted slope by the ratio of the PI estimate of the person-year weighted variance of $\psi_{j}^O$ to the corresponding KSS estimate of this quantity. Remarkably, this exercise yields a projection slope of 0.987, suggesting that, were it not for the estimation error in $\hat \psi_j^O$, the conditional averages depicted in Figure \ref{fig:reg} would be centered around the dashed 45 degree line. Converting this slope into a correlation using the KSS estimate of the person-year weighted variance of $\psi_{j}^Y$ yields a person-year weighted correlation between the two sets of firm effects of 0.89, which indicates the underlying $(\psi_{j}^Y,\psi_{j}^O)$ pairs are tightly clustered around this 45 degree line.
Theorem \ref{thm3} allows us to formally test the joint null hypothesis that the two sets of firm effects are actually identical, i.e., that both the slope and $R^2$ from a projection of $\psi_j^Y$ onto $\psi_j^O$ are one. We can state this hypothesis as $ H_0: \psi_j^O=\psi_j^Y$ for all $j\in \mathcal{J}.$ Using the test suggested in \ensuremath{^{\textrm{th}}}\ref{rem:testbig} we obtain a realized test statistic of $3.95$ which, when compared to the right tail of a standard normal distribution, yields a p-value on $H_0$ of less than $0.1\%$.
Hence, we can decisively reject the null hypothesis that older and younger workers face exactly the same vectors of firm effects. However, our earlier correlation results suggest that $H_0$ nonetheless provides a fairly accurate approximation to the structure of firm effects, at least among those firms that employ movers of both age groups.
\subsection{Inference}
We now study more carefully the problem of inference on the variance of firm effects. For convenience, the top row of Table \ref{table5} reprints our earlier KSS estimates of the variance of firm effects in each sample. Below each estimate of firm effect variance is a corresponding standard error estimate, computed according to the approach described in \ensuremath{^{\textrm{th}}}\ref{lem:varCS}. As noted in \ensuremath{^{\textrm{th}}}\ref{rem:std}, these standard errors will be somewhat conservative when there is a large share of observations for which no split sample predictions can be created. In the leave-one-out samples this share varies between 15\% and 22\%, indicating that the standard errors are likely upward biased. In the leave-two-out samples, however, this source of bias is not present as the split sample predictions always exist. The standard errors will also tend to be conservative when there is a large share of observation pairs in the set $\mathcal{B}$, for which there is upward bias in the estimator of the error variance product. However, for both the leave-one-out and leave-two-out samples, this share varies between only 0.03\% and 0.46\%, suggesting only a small degree of upward bias stems from this source.\fxnote{please check}
The next panel of Table \ref{table5} reports the $95\%$ confidence intervals that arise from setting $q=0$, $q=1$, or $q=2$. While the first interval employs a normal approximation, the latter two allow for weak identification by employing non-standard limiting distributions involving linear combinations of normal and $\chi^2$ random variables. We also report estimates of the curvature parameters $(\kappa_1,\kappa_2)$ used to construct the weak identification robust intervals. In the pooled samples both curvature parameters are estimated to be quite small, indicating that a normal approximation is likely to be accurate. Accordingly, setting $q>0$ has little discernible effect on the resulting confidence intervals in these samples. However, among older workers, particularly in the leave-two-out sample, we find stronger curvature coefficients suggesting weak identification may be empirically relevant. Setting $q>0$ in this sample widens the confidence interval somewhat and also changes its shape: mildly shortening the lower tail of the interval but lengthening the upper tail.
Treating the samples of younger and older workers as independent, the fact that the confidence intervals for the two age group samples overlap implies we cannot reject the null hypothesis that the firm effect variances are identical at the $(1-0.95^2)\times100 = 9.75\%$ level. The significance of the 0.23 log point difference between the leave-one-out and leave-two-out estimates of firm effect variance in the sample of older workers turns out to more difficult to assess. By the Cauchy-Schwartz inequality, the covariance between the leave-one-out and leave-two-out estimators is at most $(0.0026)^2(0.0014)^2=3.64\times10^{-6}$. Hence the standard error on the difference between the two estimators is at least 0.0012, which implies a maximal t-statistic of 1.92. Therefore, even when using a normal approximation, we find rather weak evidence against the null that the leave-one-out and leave-two-out estimands are equal. However, because the leave-one-out standard error estimator is likely upward biased, this finding is somewhat less conclusive than would typically be the case.
Theorem \ref{thm4} suggests two important diagnostics for the asymptotic behavior of our estimator are the Lindeberg statistics $\{ \max_i \mathsf{w}_{is}^2\}_{s=1,2} $ and the top eigenvalue shares $\{\lambda_{s}^2/\sum_{\ell=1}^r \lambda_\ell^2 \}_{s = 1,2,3}$.
The bottom panel of Table \ref{table5} reports these statistics for each sample. The top eigenvalue shares are fairly small in the pooled sample and among younger workers. A small top eigenvalue share indicates the estimator does not depend strongly on any particular linear combination of firm effects and hence that a normal distribution should provide a suitable approximation to the estimator's asymptotic behavior (i.e. that $q=0$). Accordingly, we find that the confidence intervals are virtually identical for all values of $q$ in both the pooled samples and the two samples of younger workers.
Among older workers the top eigenvalue share is 31\% in the leave-one-out sample and 58\% in the leave-two-out sample. The next largest eigenvalue share is, in both cases, less than 5\%, which suggests this is a setting where $q=1$. In line with this view, confidence intervals based upon the $q=1$ and $q=2$ approximations are nearly identical in both samples of older workers. The accuracy of these weak-identification robust confidence intervals hinges on the Lindeberg condition of Theorem \ref{thm4} being satisfied. One can think of the Lindeberg statistic $ \max_i \mathsf{w}_{is}^2 $ as giving an inverse measure of effective sample size available for estimating the linear combination of firm effects associated with the $s$'th largest eigenvalue. The fact that these statistics are all less than or equal to 0.05 implies an effective sample size of at least 20. We study in the Monte Carlo exercises below whether this effective sample size is sufficient to provide accurate coverage. Reassuringly, the sum of squared eigenvalues is quite small in all six samples considered, indicating that the leave out estimator is consistent also in our weakly identified settings.
\subsection{Monte Carlo Experiments}
\label{sec:MC}
We turn now to studying the finite sample behavior of the leave-out estimator of firm effect variance and its associated confidence intervals under a particular data generating process (DGP). Data were generated from the following first differenced model based upon equation (\ref{fd_model}):
\begin{align}
\Delta y_{g} = \Delta f_{g}' \hat \psi^{scale} + \Delta \varepsilon_{{g}}, && ({g}=1,\dots,N).
\end{align}
Here $\hat \psi^{scale}$ gives the $J\times1$ vector of OLS firm effect estimates found in the pooled leave-one-out sample, rescaled to match the KSS estimate of firm effect variance for that sample.
The errors ${\Delta \varepsilon_{{g}}}$ were drawn independently from a normal distribution with variances given by the following model of heteroscedasticity:
\begin{align}
\mathbb{V}[\Delta \varepsilon_{{g}}] = \exp(a_0 + a_1 B_{gg} + a_2 P_{gg} + a_3 \ln L_{g2} + a_4 \ln L_{g1}),
\end{align}
where $L_{gt}$ gives the size of the firm employing worker $g$ in period $t$. To choose the coefficients of this model, we estimated a nonlinear least squares fit to the ${\hat \sigma_g^2}$ in the pooled leave-one-out sample, which yielded the following estimates:
\begin{align}
\hat a_0=-3.3441, \quad \hat a_1=1.3951, \quad \hat a_2=-0.0037, \quad \hat a_3=-0.0012, \quad \hat a_4=-0.0086.
\end{align}
For each sample, we drew from the above DGP 1,000 times while holding firm assignments fixed at their sample values.
Table \ref{table6} reports the results of this Monte Carlo experiment. In accord with theory, the KSS estimator of firm effect variances is unbiased while the PI and HO estimators are biased upwards. As expected, the KSS standard error estimator exhibits a modest upward bias in the leave-one-out samples ranging from 15\% in the sample of older workers to 44\% among younger workers. In the leave-two-out sample, however, the standard error estimator exhibits biases of only 6\% or less. Unsurprisingly then, the $q=0$ confidence interval over-covers in both the pooled leave-one-out sample and the leave-one-out sample of younger workers. In the corresponding leave-two-out samples, however, coverage is very near its nominal level, both for the normal based $(q=0)$ and the weak identification robust $(q=1)$ intervals.
In the samples of older workers, the normal distribution provides a poor approximation to the shape of the estimator's sampling distribution, which is to be expected given the large top eigenvalues found in these designs. This non-normality generates substantial under-coverage by the $q=0$ confidence interval in the leave-two-out sample. Applying the weak identification robust interval in the leave-two-out sample of older workers yields coverage very close to nominal levels despite the fact that the effective sample size available for the top eigenvector is only about $20$.
In sum, the Monte Carlo experiments demonstrate that confidence intervals predicated on the assumption that $q=1$ can provide accurate size control in leave-two-out samples when the realized mobility network exhibits a severe bottleneck. We also achieved size control in leave-one-out samples, albeit at the cost of moderate over-coverage. Hence, in applications where statistical power is a first-order consideration, it may be attractive to restrict attention to leave-two-out samples, which tend to yield estimates of variance components very close to those found in leave-one-out samples but with substantially less biased standard errors.
\section{Conclusion}
We propose a new estimator of quadratic forms with applications to several areas of economics. The estimator is finite sample unbiased in the presence of unrestricted heteroscedasticity and can be accurately approximated in very large datasets via random projection methods. Consistency is established under verifiable design requirements in an environment where the number of regressors may grow in proportion to the sample size. The estimator enables tests of linear restrictions of varying dimension under weaker conditions than have been explored in previous work. A new distributional theory highlights the potential for the proposed estimator to exhibit deviations from normality when some linear combinations of coefficients are imprecisely estimated relative to others.
In an application to Italian worker-firm data, we showed that ignoring heteroscedasticity can substantially bias conclusions about the relative contribution of workers, firms, and worker-firm sorting to wage inequality. Accounting for serial correlation within a worker-firm match was found to be empirically important, while across match correlation appears to be negligible. Consequently, those studying longer panels may wish to collapse their data down to match level means and then apply the leave-observation-out estimator. Alternately, researchers can simply extract and analyze separately balanced panels of length two, which also facilitates analysis of the temporal stability of the firm and person effect variances.
Leave-out standard error estimates for the coefficients of a linear projection of firm effects onto worker and firm observables were found to be several times larger than standard errors that naively treat the estimated firm effects as independent. These results strongly suggest that researchers seeking to identify the observable correlates of high-dimensional fixed effects should consider employing the proposed standard errors, including when studying settings falling outside the traditional worker-firm setup \citep[e.g., ][]{finkelstein2016sources,chetty2018impacts}. Stratifying our analysis by birth cohort, we formally rejected the null hypothesis that older and younger workers face identical vectors of firm effects but found that the two sets of firm effects were highly correlated. Corresponding techniques can be used to study multivariate models.
A Monte Carlo analysis demonstrated that bottlenecks in the worker-firm mobility network can generate quantitatively important deviations from normality. The proposed inference procedure captured these deviations accurately with a weak identification robust confidence interval. In cases where the mobility network was strongly connected, accurate inferences were obtained with a normal approximation. Our results suggest that in typical worker-firm applications, the normal approximation is likely to suffice. However, when studying small areas, or sub-populations with limited mobility, accounting for weak identification can be quantitatively important.
\bibliographystyle{chicago}
\newpage
\bibliography{paper}
\newpage
\begin{figure}[H]
\begin{center}
\caption{Realized Mobility Network: Older workers}
\label{fig:net}
\includegraphics[width=\textwidth]{figure1.pdf}
\end{center}
{\footnotesize \underline{Note:} This figure provides a visualization of the design matrix $S_{xx}$ for the leave-two-out sample of older workers (see Table \ref{table1} for reference). The graph is plotted in the statistical software R using the \textit{igraph} package and the large-scale graph layout (DrL) using the option to concentrate firms from the same blocks. High weight mobility refers to observations that have $w_{i1}^2$ or $w_{i2}^2$ above $1/500$ and these observations form the bottlenecks between the three blocks.
}
\end{figure}
\addcontentsline{toc}{section}{Tables and Figures}
\newpage
\begin{figure}[H]
\begin{center}
\caption{Do Firm Effects Differ Across Age Groups?}
\label{fig:reg}
{\includegraphics[width = \textwidth]{figure_2.pdf}} \\
\end{center}
{\footnotesize \underline{Note:} This figure plots the mean of the estimated firm effects for younger workers ($\hat \psi_j^Y$) by centiles of the estimated firm effects for older workers ($\hat \psi_j^O$) in the sample of 8,578 firms for which both sets of effects are leave-one-out identified. Both sets of firm effects are demeaned within this estimation sample. ``PI slope'' gives the coefficient from a person-year weighted projection of $\hat \psi_j^Y$ onto $\hat \psi_j^O$. ``KSS slope'' adjusts for attenuation bias by multiplying the PI slope by the ratio of the plug-in estimate of the person-year weighted variance of $ \psi_j^O$ to the KSS adjusted estimate of the same quantity. ``PI correlation" gives the person-year weighted sample correlation between $\hat \psi_j^O$ and $\hat \psi_j^Y$ while ``KSS correlation" adjusts this correlation for sampling error in both $\hat \psi_j^O$ and $\hat \psi_j^Y$ using leave out estimates of the relevant variances. ``Test statistic'' refers to the realization of ${\hat \theta_{H_0}}/{\sqrt{\hat{\mathbb{V}}[\hat \theta_{H_0}]}}$ where $\hat \theta_{H_0}$ is the quadratic form associated with the null hypothesis that the firm effects are equal across age groups, see \ensuremath{^{\textrm{th}}}\ref{rem:testbig} and Appendix \ref{app:testOY} for details. From \ensuremath{^{\textrm{th}}}\ref{thm3}, ${\hat \theta_{H_0}}/{\sqrt{{\mathbb{V}}[\hat \theta_{H_0}]}}$ converges to a $\mathcal{N}(0,1)$ under the null hypothesis that $\psi_j^O=\psi_j^Y$ for all 8,578 firms.
}
\end{figure}
\includepdf[scale = 1, addtolist = {1,table,Table 1: Summary Statistics,table1}]{paper_revision_KSS_table1.pdf}
\includepdf[scale = 1, landscape=true, addtolist = {1,table,Table 2: Variance Decomposition,table2}]{paper_revision_KSS_table2.pdf}
\includepdf[scale = 1, landscape=true, addtolist = {1,table,Table 3: Variance Decomposition under different leave-out strategies,table3}]{paper_revision_KSS_table3.pdf}
\includepdf[scale = 1, landscape=true, addtolist = {1,table,Table 4: Projecting Firm Effects on Covariates,table4}]{paper_revision_KSS_table4.pdf}
\includepdf[scale = 1, landscape=true, addtolist = {1,table,Table 5: Inference on the Variance of Firm Effects,table5}]{paper_revision_KSS_table5.pdf}
\includepdf[scale = 1, landscape=true, addtolist = {1,table,Table 6: Montecarlo Results,table6}]{paper_revision_KSS_table6.pdf}