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.
56,568 characters
Optimal Shrinkage Estimation of Fixed Effects in Linear Panel Data Models
\maketitle
\begin{abstract}
Shrinkage methods are frequently used to improve the precision of least
squares estimators of fixed effects. However, widely used shrinkage estimators
guarantee improved precision only under strong distributional assumptions. I
develop an estimator for the fixed effects that obtains the best possible mean
squared error within a class of shrinkage estimators. This class includes
conventional shrinkage estimators and the optimality does not require
distributional assumptions. The estimator has an intuitive form and is easy to
implement. Moreover, the fixed effects are allowed to vary with time and to be
serially correlated, in which case the shrinkage optimally incorporates the
underlying correlation structure. I also provide a method to forecast fixed
effects one period ahead in this setting.
\end{abstract}
\newpage
\section{Introduction}
\label{sec:introduction}
Applied economists are often interested in unit-specific effects in linear panel
data models.\footnote{Readers are referred to \cite{walters2024empirical} for an
overview.} Estimation in such settings typically requires a large number of
unit-specific fixed effects. However, relatively small sample sizes at the unit
level yield noisy estimates of these effects. Empirical Bayes (EB) methods,
which shrink the least squares estimates, are frequently used to improve their
precision.
However, commonly used EB approaches guarantee such precision improvements only
under stringent assumptions. Moreover, existing methods that relax these
assumptions are typically limited to settings in which the true effects are
independent across units.
This restriction reduces their applicability in contexts where disaggregated
effects or effects on multiple outcomes are of interest.\footnote{While there is
a conceptual distinction between disaggregated effects and effects on multiple
outcomes, they both lead to the same econometric setting with non-exchangeable
effects within each unit. For brevity, I use the term disaggregated effects to
refer to both.} For instance, it has been noted by
\cite{chetty2014MeasuringImpactsTeachers} that allowing for a time drift is
crucial in the context of teacher value-added (TVA). Yet, since value-added is
likely correlated across years within a single teacher, methods based on the
assumption of independent effects cannot be directly
applied.\footnote{Other examples where such disaggregated effects are of
interest include insurance company-level effects
(\citealp{abaluck2020MortalityEffectsChoicea}), time-drifting college
counselor effects (\citealp{mulhern2023beyond}) and teacher effect on multiple
outcomes (\citealp{rose2022effects}).} The goal of the paper is to provide an
estimation method that 1) guarantees precision improvement under weak
conditions and 2) applies to settings with disaggregated effects. The proposed
method is derived by extending the analysis by \cite{xie2012SUREEstimatesHeteroscedastica}.
Consider the canonical EB approach, which assumes:\vspace{-5pt}
\begin{equation}
\label{eq:eb_model}
\theta_{j} \overset{i.i.d.}{\sim} N(0,
\lambda) \,\,\, \text{and} \,\,\, y_{j} \vert \theta_{j} \overset{indep}{\sim} N(\theta_{j},
\sigma^{2}_{j}),\vspace{-5pt}
\end{equation}
where $\theta_{j}$ is the true effect and $y_{j}$ is an estimator of
$\theta_{j}$ with known variance $\sigma^{2}_{j}$. An implicit assumption is
that the true effect $\theta_{j}$ is independent of the variance
$\sigma^{2}_{j}$. The importance of this prior invariance assumption has been
discussed in \cite{xie2012SUREEstimatesHeteroscedastica} and
\cite{chen2022gaussian}, with the latter developing a novel nonparametric EB
method that relaxes it. The canonical EB model (\ref{eq:eb_model}) therefore
relies on three components: normality of $\theta_{j}$, normality of
$y_{j}\vert \theta_{j}$ and independence of $\theta_{j}$ and $\sigma_{j}$. I
refer to these collectively as the \textit{parametric/normal EB
assumptions}.\footnote{I refer to such assumptions simply as \textit{EB
assumptions} hereafter.} The posterior mean of $\theta_{j}$ is
$\ \hat{\theta}_{j}(\lambda) := E[\theta_{j} \vert y_{j}] =
\frac{\lambda}{\lambda + \sigma_{j}} y_{j}$, and the EB estimator
$\hat{\theta}_{j}(\hat{\lambda})$ is obtained by replacing the unknown
\textit{hyperparameter} $\lambda$ with an estimator $\hat{\lambda}$ based on the
marginal distribution of $y_{j}$ implied by (\ref{eq:eb_model}). The risk
properties of $\hat{\theta}_{j}(\hat{\lambda})$ are therefore inherently
sensitive to these assumptions.
I propose an alternative shrinkage estimator with optimality properties that do
not rely on the EB assumptions, while retaining a simple form. The approach also
accommodates settings in which each unit-level effect can be decomposed into
multiple, possibly correlated, disaggregated effects. For convenience, I refer
to this disaggregated dimension as “time.” The proposed shrinkage method
incorporates the underlying correlation structure, in contrast to the canonical
EB estimator, which treats all effects as exchangeable.
In the context of time-varying fixed effects, I develop an optimal forecasting
method to predict the fixed effect one period ahead. A special case of this
forecasting method coincides with the estimator of
\cite{chetty2014MeasuringImpactsTeachers}.
A simple illustration of the main method begins with a multivariate version of
(\ref{eq:eb_model}):
\begin{equation}
\label{eq:eb_model_main}
\theta_{j} \overset{i.i.d.}{\sim} N(0,
\Lambda) \,\,\, \text{and} \,\,\, y_{j} \vert \theta_{j} \overset{indep}{\sim} N(\theta_{j},
\Sigma_{j}),
\end{equation}
where $\theta_{j}$ and $y_{j}$ are now $T$-dimensional vectors, and $\Lambda$
and $\Sigma_{j}$ are $T \times T$ matrices, reflecting the presence of
disaggregated effects. As before, $\Sigma_{j}$ is assumed to be known. The
posterior mean of $\theta_{j}$ is
$\hat{\theta}_{j}(\Lambda) := \operatorname{\mathbf{E}}[\theta_{j} \vert y_{j}] = \Lambda(\Lambda +
\Sigma_{j})^{-1}y_{j}$. Rather than using the marginal likelihood of $y_{j}$
implied by (\ref{eq:eb_model_main}), as in the EB approach, I follow
\cite{xie2012SUREEstimatesHeteroscedastica} and tune the hyperparameter
$\Lambda$ to minimize an estimate of the mean squared error (MSE). I refer to
this risk estimate as the unbiased risk estimate (URE), and the \textit{URE
estimator} selects hyperparameters by minimizing the URE.\footnote{Both the
strategy and terminology follow
\cite{xie2012SUREEstimatesHeteroscedastica}. Related work using this approach
includes \cite{xie2016OptimalShrinkageEstimation} and
\cite{brown2018EmpiricalBayesEstimates}.} I show
that this leads to an estimator that has minimum asymptotic MSE within the
class of estimators considered, which in this simple case is
$\{(\hat{\theta}_{j}(\Lambda))_{j=1}^{J}: \Lambda \text{ is positive
semidefinite}\}$. The estimators in this class take the intuitive form of
linearly shrinking $y_{j}$ according to their precision, and the class includes
widely used estimators such as the conventional EB estimators.
While the class of estimators is motivated by the model in
(\ref{eq:eb_model_main}), the optimality results require only a few
bounded-moment conditions and no distributional assumptions.\footnote{That said,
the method does rely on (\ref{eq:eb_model_main}), or a variation of it
introduced in later sections, to restrict the class of estimators.} The key reason is that
the URE remains “close” to the true risk under weak conditions, rendering such
assumptions unnecessary. To establish this optimality, I derive new results in a
multivariate normal means setting, which is a frequentist version of
(\ref{eq:eb_model_main})—that is, the problem of estimating the fixed mean
vectors $\{\theta_{j}\}_{j=1}^{J}$ based on observations
$y_{j} \overset{\mathrm{indep}}{\sim} N(\theta_{j}, \Sigma_{j})$ with known
$\Sigma_{j}$ for $j = 1, \dots, J$. Under heteroskedasticity, no estimator has
been shown to be risk-optimal (in the frequentist sense) unless $T = 1$, a case
handled by \cite{xie2012SUREEstimatesHeteroscedastica}. Allowing for $T > 1$ and
general forms of $\Sigma_{j}$, I derive an estimator that achieves the best
possible MSE within a class of estimators that take the form of posterior means
under (\ref{eq:eb_model_main}).
I use the proposed method to estimate a TVA model of public
schools of New York City.
I revisit the policy exercise of releasing the bottom 5\% of teachers according
to the estimated fixed effects. I find that, relative to the conventional
methods, the composition of released teachers changes by 25\% when using the
proposed estimation method and by 58\% when using the proposed forecast
method. An out-of-sample exercise shows that the average value-added of the
teachers released under the forecast method is about 20\% lower compared to the
case where the conventional estimator is used. Estimates from forecast method
are shown to reduce the MSE by 35\%, indicating that the choice of estimator
makes a significant difference, and that it is crucial to allow for the
value-added to vary with time.
\bigskip
\noindent \textbf{Notation.} Let $\lVert \cdot \rVert$ denote the Euclidean norm
for both vectors and matrices (i.e., the Frobenius norm in the latter case). For
any matrix $A$, $(A)_{ij}$ denotes its $(i,j)$ entry and $s_{k}(A)$ its $k$th
largest singular value. The set of positive semidefinite $k \times k$ matrices
is denoted by $\mathcal{S}_{k}^{+}$, and the $k \times k$ identity matrix is
denoted by $I_{k}.$ For any $\{W_{ijt}\}$, let
$\overline{W}_{jt} = {n^{-1}_{jt}}\sum_{i=1}^{n_{jt}}W_{ijt}$ denote the
$(j,t)$-level sample average.
\section{URE estimators}
\label{sec:ure-estimators}
\subsection{Fixed effects and the normal means model}
\label{sec:fixed-effects-normal}
I consider the following linear panel data model,
\begin{equation}
\label{eq:model}
Y_{ijt} = X_{ijt}'\beta + \theta_{jt} + \varepsilon_{ijt},
\end{equation}
where $t = 1, \dots, T$, $j = 1,\dots, J$, and $i = 1, \dots, n_{jt}$. Here,
$\{(Y_{ijt}, X_{ijt}')\}$ denotes the observed data, $\varepsilon_{ijt}$ the
idiosyncratic shock, and $\theta_{jt}$ the time-varying fixed effect, which is
the object of interest.\footnote{I treat $\theta_{jt}$ as random to be
consistent with the hierarchical models used to motivate the proposed
estimators. However, since no restriction on its potential dependence
structure with the observed covariates is imposed, it is still considered a
fixed effect.} Typically, $i$ denotes the individual level, $j$ the group level,
and $t$ the time dimension. The time-varying fixed effect for $j$,
$\theta_{j}:=(\theta_{j1},\dots \theta_{jT})'$ is assumed to be exchangeable
across $j$. For the idiosyncratic error terms, assume
$\overline{\varepsilon}_{j} = (\overline{\varepsilon}_{j1}, \dots,
\overline{\varepsilon}_{jT})'$ is independent across $j$ and independent of
$\theta_{j}$, with known variance $\Sigma_{j}$. In practice, a consistent
estimator is plugged in, which does not affect the asymptotic properties under
suitable conditions.
\begin{example}[Teacher value-added]
In the TVA model, $j$ corresponds to teacher, $t$ to year, and
$i$ to a student assigned to teacher $j$ in year $t$. The outcome variable
$Y_{ijt}$ is a measure of student achievement and $X_{ijt}$ is a vector of
student characteristics. The fixed effect $\theta_{jt}$ is the value-added of
teacher $j$ in year $t$
. Assuming that $\varepsilon_{ijt}$ is i.i.d across all $i$, $j$, and
$t$ with variance $\sigma^{2}_{\varepsilon}$, the variance of the estimate is
given by
$\Sigma_{j} = \sigma^{2}_{\varepsilon} \mathrm{diag}(1/n_{j1}, \dots,
1/n_{jT})$.
\end{example}
I consider the asymptotic regime where $ J \to \infty$ with $T$ and $n_{jt}$
fixed. This captures the common situation where the number of effects to be
estimated is large, with observations for each fixed effect unit being
relatively small. I assume that a consistent estimator $\hat{\beta}$ of $\beta$
is readily available, which is easy to obtain under standard assumptions (see,
for example, \citealp{wooldridge2010EconometricAnalysisCross} for a
textbook-level discussion). For example, the within estimator is consistent in
the present setting under minor regularity conditions.
Let $y_{jt}$ denote the least squares estimator for the fixed
effects
:
\begin{equation*}
\label{eq:1}
y_{jt} := \overline{Y}_{jt} - \overline{X}_{jt}'\hat{\beta} =
\overline{X}_{jt}'(\beta - \hat{\beta}) + \theta_{jt} +
\overline{\varepsilon}_{jt} = \theta_{jt} + \overline{\varepsilon}_{jt} + O_{p}(J^{-1/2}).
\end{equation*}
To see the connection between this estimator and the normal means model, note
that $y_{j}:= (y_{j1}, \dots, y_{jT})' \to_{d}
\theta_{j}+\overline{\varepsilon}_{j}$ for each $j \leq J$.\footnote{Under a mild boundedness condition on $ {X}_{ijt}$
that ensures $\sup_{j}\lVert \overline{X}_{j} \rVert = O_{p}(1)$, this
convergence is uniform over $j$.
} Further assuming that $\overline{\varepsilon}_{j}$ follows a normal
distribution (with known variance matrix $\Sigma_{j}$), it follows that
$ \left( \theta_{j}+\overline{\varepsilon}_{j} \right) \vert \theta_{j} \sim
N(\theta_{j}, \Sigma_{j})$ so that
\begin{equation}
\label{eq:normal_means_model}
y_{j} \vert \theta_{j} \sim N(\theta_{j}, \Sigma_{j}),
\end{equation}
approximately. Note that even if $\varepsilon_{ijt}$ is homoskedastic, the
different cell sizes $n_{jt}$ lead to heteroskedasticity of the estimators,
$\operatorname{\mathrm{var}} (\overline{\varepsilon}_{jt}) = \operatorname{\mathrm{var}}(\varepsilon_{ijt})/n_{jt}$. Due to
this connection, I now abstract away from the panel data model and focus on the
problem of estimating $\theta = (\theta_{1}', \dots, \theta_{J}')'$ after
observing the data $\{y_{j}\}_{j=1}^{J}$ that is generated according to
(\ref{eq:normal_means_model}). The variance matrix
$\Sigma_{j} \in \mathcal{S}_{T}^{+}$ is assumed to be known, or consistently
estimable.
\subsection{Class of shrinkage estimators}
\label{sec:class-shrink-estim}
The URE estimators will be shown to be optimal within a class of shrinkage
estimators that can be viewed as Bayes estimators under a certain hierarchical
model. To motivate the class of estimators considered, suppose that the true
effects are drawn from a normal distribution:
\begin{equation}\label{eq:second_level}
\theta_{j} \overset{\mathrm{indep}}{\sim} N(\mu_{j},\Lambda),
\end{equation}
where the location vector $\mu_{j} \in\mathbf{R}^{T} $ and the variance matrix
$\Lambda \in \mathcal{S}_{T}^{+}$ are unknown. Under
\eqref{eq:normal_means_model} and \eqref{eq:second_level}, the conditional
expectation of $\theta_{j}$ given $y_{j}$\footnote{This is the posterior mean of
$\theta_{j}$ if we interpret (\ref{eq:normal_means_model}) and
\eqref{eq:second_level} as a Bayesian model, except for the fact that
$\mu_{j}$ and $\Lambda$ are treated as unknown.} is
\begin{equation}\label{eq:post_mean}
\hat{\theta}_{j}(\mu_{j}, \Lambda) := \operatorname{\mathbf{E}}[\theta_{j} \vert y_{j}] = \left( I_{T} - \Lambda(\Lambda + \Sigma_{j})^{-1} \right)\mu_{j} + \Lambda(\Lambda + \Sigma_{j})^{-1}y_{j}
\end{equation}
The unknown parameters $\mu_{j}$ and $\Lambda$ are later tuned to achieve
desirable risk properties. The restriction one imposes on $\mu_{j}$ and
$\Lambda$ determines the class of estimators. I denote by
$\mathcal{M} \subset \mathbf{R}^{J \times T}$ and
$\mathcal{L} \subset \mathcal{S}_{T}^{+}$ the (possibly random) sets that
reflect this restriction. For example, setting $\mathcal{M} = \{(0,\dots,0)'\}$
corresponds to the common practice of shrinking to the origin, where the degree
(and direction) of shrinkage is determined by $\Lambda$. Also, if one takes
$\mathcal{L}$ to be the set of positive semidefinite Toeplitz matrices, the
estimators place more weight on estimates from nearby time periods.\footnote{In
this case, the dimension of $\Lambda$ is reduced to $T$, compared to
$T(T+1)/2$ when $\Lambda$ is left unrestricted. Hence, the restrictions can be
imposed for computational considerations as well.}
Analogous to the univariate case, I refer to $\Lambda(\Lambda +
\Sigma_{j})^{-1}$ as the shrinkage matrix. Unlike the univariate case where $T =
1$, this estimator not only shrinks the magnitude of the deviations from
$\mu_{j}$ but also applies a rotation. That is, it differentially shrinks linear
combinations of the elements of $y_{j}$, rather than each element
individually. This reflects the fact that the optimal amount of shrinkage
depends not only on the variances but also on the covariances, either directly
or through the sample sizes $n_{jt}$.
\begin{example}[Independent effects]
If $\mathcal{L} = \{\lambda I_{T}: \lambda \geq 0\}$ and
$\Sigma_{j} = \mathrm{diag}(\sigma^2_{j1}, \dots, \sigma^2_{jT})$, then the
$t$th component of $\hat{\theta}_{j}(\mu_{j}, \Lambda)$ is given as
\begin{equation*}
\left( 1 - \frac{\lambda}{\lambda + \sigma^{2}_{jt}} \right)
\mu_{jt} + \frac{\lambda}{\lambda + \sigma^{2}_{jt}} y_{jt},
\end{equation*}
which is the form of shrinkage estimators\footnote{More precisely, the
estimators used in the literature take this form without the time-varying
component, and thus effects are aggregated at the $j$ level so that the
subscript $t$ disappears.} used in the literature with a specific choice of
$\lambda$ and $\mu$.
\end{example}
\begin{example}[Perfectly correlated effects]
Let $\mathbf{1} $ denote the $T$-vector with all elements equal to
$1$. Consider the case where $\Lambda = \lambda \mathbf{1}\mathbf{1}'$. Let
$\Sigma_{j} = \sigma^{2} \mathrm{diag}(1/n_{j1}, \dots, 1/n_{jT})$ and denote
the teacher-level sample size by $n_{j} = \sum_{t=1}^{T}n_{jt}$. Then,
we have
$\hat{\theta}(0, \Lambda)= \mathbf{1} \frac{ \lambda}{\sigma^{2}/n_{j} +
\lambda} \left( \frac{1}{n_{j}}\sum_{t=1}^{T} n_{jt}y_{jt} \right)$. The
term $\frac{1}{n_{j}}\sum_{t=1}^{T} n_{jt}y_{jt}$ is a weighted mean of the
least squares estimators of teacher $j$ in year $t$, and thus is equal to the least
squares estimator for the teacher-level fixed effect without time drift. This
is the estimator used in the majority of the TVA
literature with an appropriate choice of $\lambda$.
\end{example}
Now, writing $\mu = (\mu_{1}', \dots, \mu_{J}')'$ and
$\hat{\theta}(\mu, \Lambda) = (\hat{\theta}_{1}(\mu_{1}, \Lambda)', \dots,
\hat{\theta}_{J}(\mu_{J}, \Lambda)')$, the class of estimators considered is
$$\widehat\Theta(\mathcal{M},\mathcal{L}) :=
\{\hat{\theta}(\mu, \Lambda): \mu \in \mathcal{M} \text{ and }
\Lambda \in \mathcal{L} \}.$$ This is precisely the set of estimators that take
the form of ``posterior means'' given in (\ref{eq:post_mean}), under the
restriction that $\mu \in \mathcal{M}$ and $\Lambda \in \mathcal{L}$. I refer to
$\mu$ and $\Lambda$ as the hyperparameters and $\mathcal{M}$ and $\mathcal{L}$
as the hyperparameter spaces.
I consider three different specifications of $\mathcal{M}$, with increasing
flexibility. The first specification, which is the simplest, takes
$\mathcal{M} = \mathcal{M}_{\mathrm{m}}:= \{(\overline{y}'_{J}, \dots,
\overline{y}_{J}')'\}$, where
$\overline{y}_{J} := \frac{1}{J}\sum_{j=1}^{J}y_{j}$. This corresponds to
shrinking to the grand mean. If the covariates in the panel data model
(\ref{eq:model}) included time dummies, then we have $\overline{y}_{J} = 0$,
which corresponds to the common practice of shrinking the least squares estimates
toward the origin.
The second specification takes
$\mathcal{M} = \mathcal{M}_{\mathrm{g}}:= \{(\mu_{0}', \dots, \mu_{0}')':
\mu_{0} \in \mathcal{B}\}$, which shrinks the data toward a general location
$\mu_{0}$. The restriction that $\mu_{0}$ lies in $\mathcal{B}$ is a technical
condition that ensures a certain boundedness property.\footnote{The set
$\mathcal{B}$ is defined as
$\mathcal{B} := \{\mu \in \mathbf{R}: \lvert \mu_{t}\rvert \leq
q_{1-\tau}(\{\lvert y_{jt}\rvert \}_{j=1}^{J}) \,\, \mathrm{for} \,\, t= 1,
\dots, T \},$ where $q_{1-\tau}(\{ \lvert y_{jt} \rvert \}_{j=1}^{J})$ denotes
the $1-\tau$ sample quantile of $\{ \lvert y_{jt} \rvert \}_{j=1}^{J}$. A
similar idea was used by \cite{brown2018EmpiricalBayesEstimates}. I recommend
a small $\tau$, such as $\tau = .01$. This restricts $\mu_{0}$ to be smaller
than the 99th percentile of the data in terms of magnitude. This is a mild
restriction that ensures the data are not shrunk to somewhere with almost no
observations nearby.} This estimator shrinks the data $y_{j}$ to a common but
general location $\mu_{0}$, which will later be chosen in a data-driven way to
minimize the risk.
For the third specification, consider a setting where we have additional data,
$Z_{jt} \in \mathbf{R}^{k} $. In the panel data model (\ref{eq:model}), these
are exactly the covariates that cannot be included as regressors because of the
inclusion of $(j,t)$-level fixed effects. I consider the specification
$\mathcal{M} = \mathcal{M}_{\mathrm{cov}} := \{(Z_{1}'\gamma, \dots,
Z_{J}'\gamma): \gamma \in \Gamma\}$, where $Z_{j} = (Z_{j1}, \dots, Z_{jT})'$,
which corresponds to a class of estimators that shrinks the least square
estimates $y_{j}$ to a linear combination of the covariates
$Z_{j}'\gamma$.\footnote{While I only consider the simple case of shrinking
toward a linear combination of the covariates, the estimator can be extended
in a straightforward manner to shrink toward more general functions of the
covariates.} Unlike the other two specifications, each estimate $y_{j}$ is
shrunk to a different location. Again, the restriction that $\gamma \in \Gamma$
is a technical restriction that ensures a certain boundedness property holds for
optimality results.\footnote{The set $\Gamma$ is defined as
$\Gamma := \{\gamma \in \mathbf{R}^{k}: \lVert \gamma \rVert \leq B \lVert
\hat{\gamma}^{\mathrm{OLS}} \rVert \}$ where $B$ is a large constant and
$\hat{\gamma}^{\mathrm{OLS}}$ is the OLS estimator obtained by regressing
$y_{j}$ on $Z_{j}$. This ensures that the OLS coefficient is a feasible
choice along with other coefficients with much larger magnitude as well.}
\begin{example}[Teacher value-added]
In a TVA setting, teacher (or teacher-year) level covariates are
frequently available. Such covariates cannot be included in (\ref{eq:model}),
since they are absorbed into the teacher-year fixed effects. However, one can use such
covariates to improve the precision of the teacher fixed-effect
estimates by using $\hat{\theta}^{\mathrm{URE}, \mathrm{cov}}$. Frequently available teacher level covariates include, for example,
gender, tenure, and union status of a teacher. Class size, which is almost
always available, is also an example of such covariate. Asymptotically, the inclusion
of such covariates are guaranteed to improve the MSE.
\end{example}
\subsection{Risk estimate and URE estimators}
\label{sec:risk-estimate}
To evaluate the performance of different estimators, I condition on $\theta$ and
use the compound MSE as the performance criterion. To be specific, writing the
compound loss as
$\ell(\theta, \hat{\theta}) := \frac{1}{J}(\hat{\theta} - \theta)' (\hat{\theta}
- \theta),$ the (conditional) compound MSE is given as
$R(\theta, \hat{\theta}) = \operatorname{\mathbf{E}}_{\theta} \ell(\theta, \hat{\theta})$.
This is a frequentist risk criterion since we condition on the parameters. For
what follows, I treat $\theta$ as fixed with the understanding that all claims
are conditional on a fixed sequence $\{\theta_{j}\}_{j=1}^{\infty}$.\footnote{By considering MSE
conditional on $\theta$, the optimality results established later on will not
depend on normality of $\theta$ which was used to motivate the class of estimators.} The
subscript $\theta$ in $\operatorname{\mathbf{E}} _{\theta}$ has been used to make clear that the
expectation is conditional on this sequence of $\theta$, but I omit the
subscript unless ambiguous otherwise.
Given this performance criterion, an optimal choice of the hyperparameters is
$ (\tilde{\mu}^{\mathrm{OL}}_{\mathcal{M}, \mathcal{L}},
\tilde{\Lambda}^{\mathrm{OL}}_{{\mathcal{M}, \mathcal{L}}}) := \operatorname*{arg\,min}_{(\mu,
\Lambda) \in \mathcal{M} \times \mathcal{L}} \ell(\theta, \hat{\theta}(\mu,
\Lambda))$, which gives the oracle loss ``estimator''
$\tilde{\theta}^{\mathrm{OL}}(\mathcal{M}, \mathcal{L}) :=
\hat{\theta}(\tilde{\mu}^{\mathrm{OL}}_{\mathcal{M}, \mathcal{L}},
\tilde{\Lambda}^{\mathrm{OL}}_{{\mathcal{M}, \mathcal{L}}})$. Of course, this is
infeasible because the loss function depends on the unobserved true mean
vectors. However, it turns out that the risk is estimable via Stein's unbiased
risk estimate (SURE). The idea is to choose the hyperparameters by minimizing
this risk estimate.\footnote{The idea of minimizing SURE to choose tuning
parameters has been around since at least \cite{li1985SteinUnbiasedRisk}, and
has been introduced to this setting by
\cite{xie2012SUREEstimatesHeteroscedastica}.}
The risk estimate is defined as
$\mathrm{URE}(\mu, \Lambda) = \frac{1}{J}\sum_{j}\mathrm{URE}_{j}(\mu_{j},
\Lambda)$, where
\begin{equation*}
\mathrm{URE}_{j}(\mu_{j}, \Lambda) := \operatorname{\mathrm{tr}} ( \Sigma_{j}) -
2 \operatorname{\mathrm{tr}} ((\Lambda+\Sigma_{j})^{-1}\Sigma_{j}^{2}) + (y_{j}-\mu_{j})'[(\Lambda+\Sigma_{j})^{-1}\Sigma_{j}^{2}(\Lambda+\Sigma_{j})^{-1}](y_{j}-\mu_{j}).
\end{equation*}
It is easy to show that $URE(\mu, \Lambda)$ is indeed an unbiased estimator of
the true risk as in $R(\theta, \hat{\theta}(\mu, \Lambda))$. The aim is to show
that choosing hyperparameters to minimize $\mathrm{URE}(\mu, \Lambda)$ is as
good as, in terms of asymptotic risk, choosing them by minimizing the true
loss. The URE estimator is given as
$\hat{\theta}^{\mathrm{URE}}(\mathcal{M}, \mathcal{L}) :=
\hat{\theta}(\hat{\mu}^{\mathrm{URE}}_{\mathcal{M}, \mathcal{L}},
\hat{\Lambda}^{\mathrm{URE}}_{{\mathcal{M}, \mathcal{L}}})$, where the
hyperparameters are chosen to minimize $\mathrm{URE}(\mu, \Lambda)$:
$ (\hat{\mu}^{\mathrm{URE}}_{\mathcal{M}, \mathcal{L}},
\hat{\Lambda}^{\mathrm{URE}}_{{\mathcal{M}, \mathcal{L}}}) := \operatorname*{arg\,min}_{(\mu,
\Lambda) \in \mathcal{M} \times \mathcal{L}} \mathrm{URE}(\mu, \Lambda)$.
In the EB framework, the hierarchical model given by
(\ref{eq:normal_means_model}) and (\ref{eq:second_level}) is taken as the true
data generating process, and the hyperparameters are estimated by using the
marginal distribution of the data implied by this model,
$y_{j} \overset{\text{indep}}{\sim} N(\mu_{j}, \Lambda + \Sigma_{j})$, either by
maximum likelihood or the method of moments. I denote the EB maximum likelihood
estimator (EBMLE) by
$\hat{\theta}^{\mathrm{EBMLE}}(\mathcal{M}, \mathcal{L}) =
\hat{\theta}(\hat{\mu}^{\mathrm{EBMLE}}, \hat{\Lambda}^{\mathrm{EBMLE}})$ where
$(\hat{\mu}^{\mathrm{EBMLE}}, \hat{\Lambda}^{\mathrm{EBMLE}})$ maximizes the
marginal likelihood subject to $\mu \in \mathcal{M}$ and
$\Lambda \in \mathcal{L}$.
For each of the three specifications $ \mathcal{M}_{m}$, $\mathcal{M}_{g}$ and
$\mathcal{M}_{\mathrm{cov}}$, I define the corresponding URE estimators as
$\hat{\theta}^{\mathrm{URE}, \mathrm{m}} :=
\hat{\theta}^{\mathrm{URE}}(\mathcal{M}_{\mathrm{m}}, \mathcal{S}_{T}^{+})$,
$\hat{\theta}^{\mathrm{URE}, \mathrm{g}} :=
\hat{\theta}^{\mathrm{URE}}(\mathcal{M}_{\mathrm{g}}, \mathcal{S}_{T}^{+})$ and
$\hat{\theta}^{\mathrm{URE}, \mathrm{cov}} :=
\hat{\theta}^{\mathrm{URE}}(\mathcal{M}_{\mathrm{cov}}, \mathcal{S}_{T}^{+})$. The
corresponding oracle estimators $\tilde{\theta}^{\mathrm{OL}, \mathrm{m}}$,
$\tilde{\theta}^{\mathrm{OL}, \mathrm{g}}$ and
$\tilde{\theta}^{\mathrm{OL}, \mathrm{cov}}$ are defined analogously.
\section{Optimality of the URE estimators}
\label{sec:optim-results-shrink}
I now show that the URE estimators defined in Section \ref{sec:risk-estimate}
achieve the smallest possible asymptotic MSE among all estimators
in the corresponding classes. I provide these optimality results for
$\hat{\theta}^{\mathrm{URE}, \mathrm{m}}$ and
$\hat{\theta}^{\mathrm{URE}, \mathrm{g}}$. An analogous result for
$\hat{\theta}^{\mathrm{URE}, \mathrm{cov}}$ is deferred to Appendix \ref{appsec:opt-cov},
as the arguments are similar. The main step in establishing such optimality
is to show that the corresponding UREs are uniformly close to the true
risk. Since I use an unbiased estimate of risk, this essentially reduces to a
uniform (weak) law of large numbers (ULLN) argument. I first establish a simple
high-level result for a generic URE estimator and then verify that the
conditions for this result hold for each estimator under appropriate lower-level
conditions.
The following lemma shows that if $\mathrm{URE}(\mu, \Lambda)$ is uniformly close to the
true loss in $L^{1}$, in the sense that
\begin{equation}
\label{eq:URE_loss_l1}
\sup_{(\mu, \Lambda) \in
\mathcal{M} \times \mathcal{L}} \,\,
\lvert \mathrm{URE}(\mu, \Lambda) - \ell(\theta,
\hat{\theta}(\mu, \Lambda)) \rvert \overset{L^{1}}{\to} 0,
\end{equation}
then the URE estimator has asymptotic risk as good as the oracle.
\begin{lemma}\label{lem:L1_oracle_generic} Suppose (\ref{eq:URE_loss_l1}) holds. Then,
\begin{equation}\label{eq:optimal_risk}
\underset{J \to \infty}{\lim} \,\big( R(\theta,
\hat{\theta}^{\mathrm{URE}}(\mathcal{M}, \mathcal{L})) - R(\theta, \tilde{\theta}^{\mathrm{OL}}(\mathcal{M}, \mathcal{L})) \big) = 0.
\end{equation}
\end{lemma}
\begin{proof}
See Appendix \ref{appsec:proof-theorem-l1-gen}.
\end{proof}
\noindent
Therefore, establishing uniform consistency in (\ref{eq:URE_loss_l1}) for
each estimator class is a key step. Then, since
$ R(\theta, \tilde{\theta}^{\mathrm{OL}}(\mathcal{M}, \mathcal{L})) =
\min_{\hat{\theta} \in \widehat\Theta(\mathcal{M}, \mathcal{L})} R(\theta,
\hat{\theta})$, it follows from Lemma \ref{lem:L1_oracle_generic} that the URE
estimators are asymptotically optimal within the class of estimators
$\widehat\Theta(\mathcal{M}, \mathcal{L})$. The difference between the URE and
the true loss can be decomposed as
\begin{align}
& \,\, \mathrm{URE}(\mu, \Lambda) - \ell(\theta, \hat{\theta}(\mu, \Lambda)) \nonumber \\
= & \textstyle ( \mathrm{URE}(0, \Lambda) - \ell(\theta, \hat{\theta}(0, \Lambda))) - \frac{2}{J} \sum_{j=1}^{J}\left( \mu_{j}' (\Lambda + \Sigma_{j})^{-1} \Sigma_{j}(y_{j}-\theta_{j}) \right). \label{eq:URE_minus_loss_genmu}
\end{align}
I show that both terms in the last line converge to $0$ in $L_{1}$, uniformly
over each hyperparameter space I consider. Note that the first term of the last
line does not depend on $\mu$, and thus is common for all three estimators.
The following assumption states that $y_{j}$'s are independent with (uniformly)
bounded fourth moments, and that the smallest eigenvalue of the variance of
$y_{j}$ is bounded away from zero. I write $y_{j} \sim (\theta_{j}, \Sigma_{j})$
to denote that $y_{j}$ follows a distribution such that
$ \operatorname{\mathbf{E}} y_{j} = \theta_{j} $ and $ \operatorname{\mathrm{var}} (y_{j}) = \Sigma_{j}$. The supremum
$\sup_{j}$ is taken over all $j \geq 1$, and likewise for $\inf_{j}$. Hence, the
assumption imposes conditions on the sequences
$ \left\{ \operatorname{\mathbf{E}} \lVert y_{j} \rVert \right\}_{j=1}^{\infty}$ and
$ \left\{s_{T}(\Sigma_{j}) \right\}_{j=1}^{\infty}$.
\begin{assumption}[Independent sampling and boundedness]
\label{assum:bounded}
$ \mathrm{(i)} \, y_{j} \overset{\mathrm{indep}}{\sim}
(\theta_{j}, \Sigma_{j})$, \\$ \mathrm{(ii)} \, \sup_{j} \operatorname{\mathbf{E}} \lVert y_{j} \rVert^{4} <
\infty$ and $ \mathrm{(iii)}\, 0 < \textstyle\inf_{j}s_{T}(\Sigma_{j}).$
\end{assumption}
In the case where $\Sigma_{j}$ is diagonal for all
$j$, Assumption \ref{assum:bounded} (iii) boils down to assuming that $\operatorname{\mathrm{var}}
(y_{jt})$ is bounded away from zero over $j$ and
$t$. Also, in the case where $\Sigma_{j} = \Sigma$ for all
$j$, the assumption trivially holds as long as $\Sigma$ is
invertible.
It turns out that Assumption \ref{assum:bounded} is enough to ensure uniform
convergence of the first term of \eqref{eq:URE_minus_loss_genmu}.
To show convergence of the second term of (\ref{eq:URE_minus_loss_genmu}), note
that for both $\hat{\theta}^{\mathrm{URE}, \mathrm{m}}$ and
$\hat{\theta}^{\mathrm{URE}, \mathrm{g}}$, the centering term $\mu_{j}$ does not
depend on $j$, so that I can write $\mu_{0} = \mu_{j}$ for all $j$. Hence, the
required convergence result is
\begin{equation}
\label{eq:sec_term}
\operatorname{\mathbf{E}}\Big[\textstyle\sup_{\mu_{0} \in \mathcal{M}_{0}, \Lambda \in \mathcal{S}^{+}_{T}} \left\lvert \frac{1}{J} \sum_{j=1}^{J} \mu_{0}' (\Lambda + \Sigma_{j})^{-1}
\Sigma_{j}(y_{j}-\theta_{j}) \right\rvert \Big] \to 0
\end{equation}
where $\mathcal{M}_{0} = \{\overline{y}_{J}\}$ for
$\hat{\theta}^{\mathrm{URE}, m}$ and $\mathcal{M}_{0} = \mathcal{B}$ for
$\hat{\theta}^{\mathrm{URE}, g}$. It is clear that some form of boundedness
condition on $\mathcal{M}_{0}$ is necessary for such a convergence result to
hold. For $\hat{\theta}^{\mathrm{URE}, \mathrm{m}}$, Assumption
\ref{assum:bounded} (ii) ensures this. For $\hat{\theta}^{\mathrm{URE},
\mathrm{g}}$, we require conditions that guarantee boundedness of $\mathcal{B}$
in a suitable sense, and the following assumption provides such a guarantee.
\begin{assumption}[Bounded sample quantiles]\label{assu:bdd_quan}
$ \limsup_{J}\operatorname{\mathbf{E}} q_{1-\tau}(\{ y_{jt}^{2} \}_{j=1}^{J}) < \infty$.
\end{assumption}
This assumption states that the expectation of the sample quantile of
$\{ y_{jt}^{2} \}_{j=1}^{J}$ is uniformly bounded. This is a rather mild
condition that is satisfied, for example, if the tail probabilities of
$y_{jt}^{2}$ vanish uniformly to 0. A lower level condition that is stronger but
easy to interpret is provided in Appendix \ref{sec:suff-cond-assumpt}. We now
state the optimality of $\hat{\theta}^{\mathrm{URE}, \mathrm{m}}$ and
$\hat{\theta}^{\mathrm{URE}, \mathrm{g}}$.
\begin{theorem}[Optimality of URE estimators]\label{thm:URE_opt_m_g} Suppose Assumption
\ref{assum:bounded} holds. Then,\\[.5ex]
(i) $\sup_{\Lambda \in \mathcal{S}_{T}^{+}} \big|
\mathrm{URE}(\overline y_{J}, \Lambda) - \ell(\theta, \hat{\theta}(\overline y_{J},
\Lambda)) \big| \overset{L^{1}}{\to} 0,$ and \\[.5ex]
(ii) if Assumption \ref{assu:bdd_quan} holds as well, then $\sup_{\mu \in \mathcal{B}, \Lambda \in \mathcal{S}_{T}^{+}} \big|
\mathrm{URE}(\mu, \Lambda) - \ell(\theta, \hat{\theta}(\mu,
\Lambda)) \big| \overset{L^{1}}{\to} 0.$
\end{theorem}
\begin{proof}
See Appendix \ref{sec:proof-main-theorems}.\footnote{\label{fn:
diff_with_xie}The proof technique used in related papers (e.g.,
\citealp{xie2012SUREEstimatesHeteroscedastica,
xie2016OptimalShrinkageEstimation})
does not go through in this setting mainly because 1)
the matrix hyperparameter $\Lambda$ governs the direction of
shrinkage as well as the magnitude and 2) normality of
$y_{j}$ is not assumed.
}
\end{proof}
The proof of this theorem, provided in Appendix
\ref{appsec:proof-theorem-URE}, relies on establishing a ULLN over
independent but non-identically distributed sequences of data, followed by
verifying uniform integrability to strengthen the mode of convergence from
convergence in probability to convergence in $L^{1}$. Note that the uniform
convergence is shown over the largest possible hyperparameter space for
$\Lambda$, $\mathcal{S}_{T}^{+}$, and thus convergence over any subset
$\mathcal{L} \subset \mathcal{S}_{T}^{+}$ follows immediately. I note that
Assumption \ref{assum:bounded} is stronger than necessary. However, this
stronger assumption is not particularly restrictive and has the clear
advantages of simplifying the proofs and being easy to
interpret.
By Lemma \ref{lem:L1_oracle_generic}, it follows that
$\hat{\theta}^{\mathrm{URE},\mathrm{m}}$ and
$\hat{\theta}^{\mathrm{URE},\mathrm{g}}$ are asymptotically optimal within
$\widehat\Theta(\mathcal{M}_{m}, \mathcal{S}_{T}^{+})$ and
$\widehat\Theta(\mathcal{M}_{g}, \mathcal{S}_{T}^{+})$, respectively. Note
that the optimality of the URE estimators requires only mild conditions on
the moments of the data, which is in contrast with the EB estimators that
require stringent distributional assumptions. The EB estimators are optimal
in the sense of \cite{robbins1964EmpiricalBayesApproach}\footnote{The
estimator obtains the Bayes risk of the model
(\ref{eq:normal_means_model}) and (\ref{eq:second_level}).} only when 1)
the normality assumptions for both the least squares estimator and the true
fixed effect hold and 2) the true fixed effect and variance of the least
squares estimator are independent.
The normality assumption on the true fixed effect is typically difficult to
justify.\footnote{Some evidence on the violation of such assumption in the
context of teacher value-added is provided in
\cite{gilraine2020NewMethodEstimatinga}.} The optimality results here are
conditional on a sequence of true mean vectors that satisfy a mild
boundedness condition. The independence between the true fixed effect and
the variance of the least squares estimator can be easily violated in
empirical settings as well. Since the variance of the least squares
estimator is inversely proportional to $n_{jt}$, the assumption is violated
if the fixed effect depends on $n_{jt}$. If teachers with
higher value-added teach more students, or if the size of the class is
related to teaching effectiveness, then such independence is unlikely to
hold.
The nonparametric EB literature (e.g.,
\citealp{jiang2009GeneralMaximumLikelihood},
\citealp{brown2009NonparametricEmpiricalBayes}) provides an alternative
method to relax the normality assumption on the true effects. In this
setting, the distribution of the true fixed effect is left unspecified,
thereby broadening the class of estimators. This approach, however, is
complementary to—but does not dominate—the URE approach, for two main
reasons. First, the risk properties of currently available nonparametric EB
methods still rely on the independence of the true fixed effect and the
variance of the least squares estimator (or at least a structured
relationship between the two, as in \citealp{chen2022gaussian}), as well as
on a normality assumption for the least squares estimators. Second, existing
approaches address only the case $T=1$. Even when extended to $T>1$, the
associated computation is likely to be infeasible for even moderate values
of $T$.
The URE estimators can also be shown to dominate the unshrunk (unbiased)
estimator, $\hat{\theta}^{\mathrm{ub}} =y$, which corresponds to using the
least squares estimators without any shrinkage in the context of fixed
effects. Because there is no $\Lambda \in \mathcal{S}_{T}^{+}$ such that
$\hat{\theta}(\mu_{j},\Lambda) = y$, the estimator
$\hat{\theta}^{\mathrm{ub}}$ is not included in any of the classes of
estimators I consider. However, a simple approximation argument can be used
to establish that
${\limsup}_{J \to \infty}\big( R(\theta, \hat{\theta}^{\mathrm{URE},
\mathrm{m}}) - R(\theta, \hat{\theta}^{\mathrm{ub}}) \big) \leq
0$.
This shows that using
$\hat{\theta}^{\mathrm{URE},m}$ is at least as good as using the unshrunk
estimator, providing a strong justification for shrinkage when the goal is
to improve precision. This is a property that EB methods do not enjoy unless
the EB assumptions are satisfied.
\section{Forecasting $\theta_{T+1}$}
\label{sec:optim-shrink-pred}
In settings where $t$ represents time, forecasts of future effects are often of
interest. For example, in teacher evaluation and retention, policies based on
forecasts of future effects—rather than past effects—can improve future student
outcomes. To this end, I consider the problem of predicting
$\theta_{T+1} = (\theta_{1,T+1}, \dots, \theta_{J,T+1})'$. The approach is
analogous to the URE framework proposed earlier: I first derive a class of
predictors based on a hierarchical model, and then tune the hyperparameters by
minimizing an unbiased estimate of the prediction error (UPE). The resulting
forecasts are referred to as UPE forecasts.
For simplicity, I consider the case where the fixed effects are demeaned, i.e.,
(\ref{eq:second_level}) with $\mu_{j} = 0$. Write $\Lambda$ and $\Sigma_{j}$ in
the following block matrix form:
\begin{equation*}
\Lambda =
\begin{pmatrix}
\Lambda_{-T} & \Lambda_{T,-T} \\
\Lambda_{T,-T}' & \lambda_{T}
\end{pmatrix}, \Sigma_{j} =
\begin{pmatrix}
\Sigma_{j,-T} & \Sigma_{j,T,-T} \\
\Sigma_{j,T,-T}' & \Sigma_{j,T}
\end{pmatrix}
=
\begin{pmatrix}
\Sigma_{j,1} & \Sigma_{j,1,-1}' \\
\Sigma_{j,1,-1} & \Sigma_{j,-1}
\end{pmatrix}
\end{equation*}
where $\Lambda_{-T}$, $\Sigma_{j,-T}$ and $\Sigma_{j,-1}$ are
$(T-1)\times (T-1)$ matrices.
A recommended choice for the hyperparameter space $\mathcal{L} \subset
\mathcal{S}^{+}_{T}$ is
\begin{equation*}
\mathcal{L} := \Big\{ \Lambda \in \mathcal{S}_{T}^{+}: s_{1}(\Lambda) \leq K
s_{1}\textstyle\Big( \frac{1}{J}\sum_{j=1}^{J}y_{j}y_{j}' \Big) \Big\}
\end{equation*}
for some large number $K$ that does not depend on $J$. Under the hierarchical
model,
$\frac{1}{J}\sum_{j=1}^{J}\operatorname{\mathbf{E}} y_{j}y_{j}' = \Lambda +
\frac{1}{J}\sum_{j=1}^{J}(\theta_{j}\theta_{j}' + \Sigma_{j})$, and thus
$\frac{1}{J}\sum_{j=1}^{J}y_{j}y_{j}'$ gives a sense of the scale of
$\Lambda$. By scaling this up by $K$, the bound becomes less restrictive. This
makes $\mathcal{L}$ bounded in a certain sense, which is necessary for the
optimality
argument.
The aim is to tune the hyperparameter in a way that it minimizes prediction
error of predicting $\theta_{T+1} := (\theta_{1,T+1}, \dots, \theta_{J,
T+1})'$. However, the challenge is that an unbiased estimator of this
prediction error is unavailable because we do not observe data for period
$T+1$. The strategy is to tune the hyperparameters by considering the problem of
predicting $\theta_{T} = (\theta_{1T}, \dots, \theta_{JT})'$ using the first
$T-1$ periods of data. Under a suitable stationarity condition, this will lead
to optimal hyperparameter selection for predicting $\theta_{T+1}$ with
$y_{j,-1}$ as well.
Consider the problem of estimating $\theta_{T}$ using observations from the
first $T-1$ periods. Let
$y_{j,-t} = (y_{j1}, \dots, y_{j, t-1}, y_{j,t+1}, \dots y_{j,T})'$ and
$y_{-t} = (y_{1,-t}', \dots, y_{J,-t}')'$ denote the vectors $y_{j}$ and $y$
with period-$t$ observations removed, respectively. The
class of estimators I consider is the posterior mean implied by the hierarchical
model,\vspace*{-5pt}
\begin{equation*}
\operatorname{\mathbf{E}} [\theta_{jT}\vert y_{-T}] = \Lambda_{T,-T}' (\Lambda_{-T} + \Sigma_{j,-T})^{-1}y_{j,-T}.\vspace*{-5pt}
\end{equation*}
\noindent Define the multiplicative factor as
$B(\Lambda,\Sigma_{-T}): = (\Lambda_{-T} + \Sigma_{-T})^{-1} \Lambda_{T,-T}$. The
performance criterion is the mean prediction error,
$\operatorname{\mathbf{E}} \mathrm{PE}(\Lambda; T)$, where the prediction error is given as\vspace*{-5pt}
\begin{equation*}
\mathrm{PE}(\Lambda; T) := \textstyle \frac{1}{J} \sum _{j=1}^{J}(B(\Lambda,\Sigma_{j,-T})'y_{j,-T} - \theta_{jT})^{2}.\vspace*{-5pt}
\end{equation*}
Similar to the URE estimator, I derive an estimator of the prediction error and
choose $\Lambda$ by minimizing this. An unbiased estimator of the mean
prediction error is given as
\begin{equation*}
\textstyle \mathrm{UPE}(\Lambda) = \frac{1}{J}\sum_{j=1}^{J} \left( (B(\Lambda, \Sigma_{j,-T})'y_{j,-T} -y_{jT})^{2} - \Sigma_{jT} + 2B(\Lambda, \Sigma_{j,-T})'\Sigma_{j,T,-T} \right).
\end{equation*}
Writing $\hat{\Lambda}^{\mathrm{UPE}}$ as the $\Lambda$ that minimizes
$\mathrm{UPE}(\Lambda)$, the proposed estimator for $\theta_{j,T+1}$ is
$B(\hat{\Lambda}^{\mathrm{UPE}}, \Sigma_{j,-T})'y_{j,-1}$.
\begin{remark}[Estimator of \citealp{chetty2014MeasuringImpactsTeachers}] While
I focus on predicting $\theta_{T}$ with the observations from the first $T-1$
periods, one can also consider predicting $\theta_{t}$ with observations
excluding the period $t$ observation. If $\Sigma_{j} = \Sigma$ with $\Sigma$
being diagonal, the $\Lambda$ that minimizes $\mathrm{UPE}(\Lambda)$ implies
$B(\Lambda, \Sigma_{-t}) = \hat{\beta}^{\mathrm{\mathrm{OLS}},t}$, which is
the OLS estimator of regressing $y_{jt}$ on $y_{j,-t}$. This leads to the
estimator used in \cite{chetty2014MeasuringImpactsTeachers}.
\end{remark}
Since the goal is to forecast $\theta_{T+1}$ rather than $\theta_{T}$, it is
necessary to establish that $\mathrm{UPE}(\Lambda)$ is a good estimator of the
prediction error for the problem of predicting $\theta_{T+1}$, \vspace{-5pt}
\begin{equation*}
\mathrm{PE}(\Lambda;T+1) = \textstyle \frac{1}{J} \sum_{j=1}^{J}(B(\Lambda, \Sigma_{j,-1})'y_{j,-1} - \theta_{j,T+1})^{2}. \vspace{-5pt}
\end{equation*}
By the same argument made by Lemma \ref{lem:L1_oracle_generic}, if \vspace{-5pt}
\begin{equation}\label{eq:UPE_conv}
\sup_{\Lambda \in \mathcal{L}}\,\, \lvert \mathrm{UPE}(\Lambda) -
\mathrm{PE}(\Lambda; T+1) \rvert \overset{L^{1}}{\to} 0, \vspace{-5pt}
\end{equation}
then $B(\hat{\Lambda}^{\mathrm{UPE}}, \Sigma_{j,-T})'y_{j,-1}$ obtains the
oracle mean prediction error, which is the mean prediction error of
$B(\tilde{\Lambda}, \Sigma_{j,-T})'y_{j,-1}$ with $ \tilde{\Lambda} := \operatorname*{arg\,min}_{\Lambda}\mathrm{PE}(\Lambda; T+1)$.
Due to the extrapolative nature of the method, a suitable stationarity
assumption is necessary. To formalize this notion of stationarity, I assume that
the pairs $\{ ((\theta_{j}', \theta_{j,T+1})', \Sigma_{j}) \}_{j=1}^{\infty}$
are drawn randomly from a joint density $f_{(\theta', \theta_{T+1})',
\Sigma}$. Let $f_{\Sigma}$ denote the marginal density of $\Sigma_{j}$ and let
$\mathrm{supp}(f_{\Sigma})$ denote its support. The following assumption is a
modified version of Assumption \ref{assum:bounded} that accounts for the fact
that both $\theta_{j}$ and $\Sigma_{j}$ are now considered to be random draws. \vspace{-5pt}
\begin{assumption}[Assumption \ref{assum:bounded} with random parameters\vspace{-.17in}] \label{assum:bounded_rand}
\begin{align*}
\mathrm{(i)}\,\, & \,\, y_{j} \vert \theta_{j}, \Sigma_{j} \overset{\mathrm{indep}}{\sim} (\theta_{j}, \Sigma_{j}), \hspace{4.9in}\\
\mathrm{(ii)}\, & \,\, \textstyle \sup_{j} \operatorname{\mathbf{E}} [\lVert
y_{j} \rVert^{4} \vert \theta_{j}, \Sigma_{j}] <
\infty
\text{, and} \hspace{4.9in}\\
\mathrm{(iii)} & \,\, \mathrm{supp}(f_{\Sigma}) \subset \{\Sigma \in \mathcal{S}_{T}^{+} : s_{T}(\Sigma) > \underline{s}_{\Sigma} \} \text{ for some } \underline{s}_{\Sigma} > 0.
\end{align*}\vspace*{-35pt}
\end{assumption}
To state the stationarity assumption, let $f_{\theta,\Sigma_{-T}}$ and
$f_{(\theta_{-1}', \theta_{T+1})',\Sigma_{-1}}$ denote the marginal densities
that correspond to $(\theta_{j},\Sigma_{j,-T})$ and
$((\theta_{j,-1}' , \theta_{j, T+1})',\Sigma_{j,-1})$, respectively. The
following assumption states that the distributions of
$(\theta_{j},\Sigma_{j,-T})$ and
$((\theta_{j,-1}' , \theta_{j, T+1})',\Sigma_{j,-1})$ are the same.
\begin{assumption}[Stationarity] $f_{\theta,\Sigma_{-T}} = f_{(\theta_{-1}',\theta_{T+1})',\Sigma_{-1}}$.\label{assu:stationarity_main}
\end{assumption}
\noindent Note that this assumption does not imply that the estimates $y_{j}$
themselves are stationary, and no restriction is imposed on the joint
distribution of the mean and variance.
The following theorem shows that these assumptions ensure that
\eqref{eq:UPE_conv} holds almost surely, where the almost sure convergence is
with respect to the randomness of the sequence
$\{ ((\theta_{j}', \theta_{j,T+1})', \Sigma_{j}) \}_{j=1}^{\infty}$. An
implication is that $B(\hat{\Lambda}^{\mathrm{UPE}}, \Sigma_{j,-T})'y_{j,-1}$
obtains the oracle (asymptotic) prediction error.
\begin{theorem}\label{thm:forecasting} Under Assumptions
\ref{assum:bounded_rand} and \ref{assu:stationarity_main},
(\ref{eq:UPE_conv}) holds almost surely.
\end{theorem}
\section{An application to teacher value-added}
\label{sec:empirical-study}
I now apply the proposed methods to estimate TVA in the public
schools of New York City (NYC). I show that allowing value-added to vary over
time and using the URE estimators (and forecasts) lead to substantially
different empirical results compared to the conventional approach.
\subsection{Baseline model and data}
\label{sec:baseline-model-data}
I use a standard TVA model specified as the linear panel data
model introduced in \eqref{eq:model}, where $Y_{ijt}$ is the (standardized) test
score in either English Language Arts (ELA) or math, and $X_{ijt}$ is a vector
of student characteristics. The covariates include: the previous year’s test
score, gender, ethnicity, special education status (SWD), English language
learner status (ELL), and eligibility for free or reduced-price lunch
(FL).\footnote{The results are not sensitive to which covariates are included
and/or interacted, as long as the previous year’s test scores are included.}
The teacher fixed effect $\theta_{jt}$ is indexed by time, allowing it to
vary over time. The error term $\varepsilon_{ijt}$ is assumed to be
i.i.d. across $i$, $j$, and $t$.
I use administrative data on all NYC public schools from academic years
2012–2013 to 2018–2019. As in \cite{bitler2019TeacherEffectsStudenta}, I
restrict the sample to 4th and 5th grade students. The analysis is based on ELA
scores, though results using math scores are similar. I further restrict the
sample to students whose ELA teachers are observed in all six years. The final
dataset includes $J = 1{,}185$ teachers across $T = 6$ years and 174{,}239
student-year observations. The coefficients are estimated using OLS with fixed
effects.
The average number of students per teacher per year is approximately 24.5, with
a standard deviation of about 11.7. This substantial variation in class size
implies considerable heteroskedasticity in the least squares
estimates. Moreover, regressing $\hat{\theta}_{jt}$ on class size $n_{jt}$
reveals a significant positive relationship, suggesting a potential dependence
between the variance of the least squares estimator and the true fixed
effect. Also, recent work by \cite{gilraine2020NewMethodEstimatinga} has noted
that the true value-added is unlikely to follow a normal distribution. These
patterns indicate that conventional EB approaches are likely suboptimal in this
context.
\subsection{Estimation results and policy exercise}
\label{sec:estimation-results}
\begin{figure}[t]
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=\textwidth]{est_distn}
\caption{Density plots for shrinkage estimates}
\label{fig:est_distn}
\end{subfigure}
\begin{subfigure}{.5\textwidth}
\centering
\includegraphics[width=1\textwidth]{shr_pattern}
\caption{Shrinkage patterns
}
\label{fig:shr_pattern}
\end{subfigure}\vspace{-10pt} \label{fig:shr_res}
\caption{Shrinkage results. The left panel shows the distribution of value-added
estimates across four estimators. The right panel shows average trajectories for
a group of teachers.}
\end{figure}
Figure \ref{fig:est_distn} shows the distribution of TVA
estimates under four different estimators: the conventional estimator (EBMLE
assuming time-invariant value-added; green), the EBMLE and URE estimators under
time-varying value-added (red and black), and the optimal UPE forecast
(blue).\footnote{For the EBMLE and URE estimators, I use
$\mathcal{M} = \mathcal{M}_{\mathrm{g}}.$} For the time-varying estimators,
the average across time within each teacher is used for comparability. All
shrinkage-based estimators, except for the UPE forecast, yield distributions
more concentrated around the mode than the least squares estimator (black dashed
line), reflecting the effect of shrinkage. Notably, the density plots indicate
substantial differences between the conventional method and those allowing for
time variation.\footnote{The UPE forecasts are more dispersed than the other
estimators, as they do not involve any temporal averaging. Under a
time-invariant fixed effect model, forecasts and current estimates coincide,
but here the divergence in distributions suggest that it may be misleading to
use an average of past performance to predict future value-added.}
Figure \ref{fig:shr_pattern} illustrates the shrinkage pattern of the URE
estimator.\footnote{Due to confidentiality restrictions, estimates are averaged
across a group of teachers.} Although the estimators are allowed to shrink
toward a general location ($\mathcal{M} = \mathcal{M}_{\mathrm{g}}$), the
optimal location turns out to be close to zero, so the shrinkage can be
interpreted as toward a horizontal line at zero. The URE estimator does not
shrink each point independently, but rather shrinks a smoothed version of the
trajectory, reflecting the structure of $\hat{\Lambda}^{\mathrm{URE}}$, which
has positive off-diagonal entries—consistent with positive serial correlation in
true value-added. In contrast, ignoring this correlation leads to over-shrinkage
by applying independent shrinkage at each time period. This highlights the
importance of allowing for serial dependence in the shrinkage procedure.
A standard policy simulation in the literature involves removing teachers in the
bottom 5\% of the value-added distribution and replacing them with average
teachers. I replicate this exercise, focusing on how the composition of the
bottom 5\% changes depending on the choice of estimator. Specifically, I compare
the sets of 60 teachers flagged for removal under three estimators: the
conventional time-invariant EBMLE, the URE estimator with time-varying
value-added, and the UPE forecast. Switching from the conventional estimator to
the URE changes the list by about 24\% (14 teachers), indicating that
alternative estimation strategies can significantly alter policy
outcomes.\footnote{In contrast, \cite{gilraine2020NewMethodEstimatinga} find
that a nonparametric EB method (assuming time-invariant value-added) leads to
minimal changes. This suggests that time variation is likely the main driver
of these differences.}
Moreover, when the goal of the policy is to improve student outcomes in the
following year, forecasts for the next period's value-added are arguably more
informative than a summary of past performance. When the value-added is allowed
to vary with time, one can use the UPE forecasts in such context. On the other
hand, if one specifies value-added to be time-invariant, past and future
value-added are the same by definition, and thus will release the bottom 5\%
according to the conventional estimator. This consideration turns out have a
large effect, with only 25 teachers (approximately 42\%) being released under
both estimators.
I also conduct an out-of-sample policy exercise: teachers are ranked based on
estimates using data from the first five years, and the bottom 5\% are removed
under each estimator. I then evaluate performance by measuring the average
value-added in year six using the least squares estimator as a proxy for the
true value-added. Again, the composition of removed teachers differs
substantially, with only 60\% overlap between the conventional and forecast-based
methods. Importantly, the average value-added of the removed teachers is 20\%
lower when using the UPE forecasts. This suggests that the different composition
is in the correct direction, highlighting that the choice of estimator leads to
more effective policy. Consistent with this observation, the MSE of the UPE
forecasts is also 35\% lower than that of the conventional estimator, further
underscoring the importance of incorporating time variation and the optimality
of the proposed methods.
\bibliography{EB.bib}
\newpage