Optimal Shrinkage Estimation of Fixed Effects in Linear Panel Data Models

Optimal Shrinkage Estimation of Fixed Effects in Linear Panel Data Models

Introduction

Applied economists are often interested in unit-specific effects in linear panel data models.\footnote{Readers are referred to 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 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 (abaluck2020MortalityEffectsChoicea), time-drifting college counselor effects (mulhern2023beyond) and teacher effect on multiple outcomes (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 xie2012SUREEstimatesHeteroscedastica.

Consider the canonical EB approach, which assumes:

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 xie2012SUREEstimatesHeteroscedastica and chen2022gaussian, with the latter developing a novel nonparametric EB method that relaxes it. The canonical EB model ((ref)) 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 parametric/normal EB assumptions.\footnote{I refer to such assumptions simply as 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 hyperparameter $\lambda$ with an estimator $\hat{\lambda}$ based on the marginal distribution of $y_{j}$ implied by ((ref)). 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 chetty2014MeasuringImpactsTeachers.

A simple illustration of the main method begins with a multivariate version of ((ref)):

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)), as in the EB approach, I follow 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 URE estimator selects hyperparameters by minimizing the URE.\footnote{Both the strategy and terminology follow xie2012SUREEstimatesHeteroscedastica. Related work using this approach includes xie2016OptimalShrinkageEstimation and 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)), the optimality results require only a few bounded-moment conditions and no distributional assumptions.\footnote{That said, the method does rely on ((ref)), 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))—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 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)).

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.

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.

URE estimators

Fixed effects and the normal means model

I consider the following linear panel data model,

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.

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, 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 :

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

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)). The variance matrix $\Sigma_{j} \in \mathcal{S}_{T}^{+}$ is assumed to be known, or consistently estimable.

Class of shrinkage estimators

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:

where the location vector $\mu_{j} \in\mathbf{R}^{T} $ and the variance matrix $\Lambda \in \mathcal{S}_{T}^{+}$ are unknown. Under (ref) and (ref), the conditional expectation of $\theta_{j}$ given $y_{j}$\footnote{This is the posterior mean of $\theta_{j}$ if we interpret ((ref)) and (ref) as a Bayesian model, except for the fact that $\mu_{j}$ and $\Lambda$ are treated as unknown.} is

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}$.

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)), 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)) 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 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)), 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.}

Risk estimate and URE estimators

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 li1985SteinUnbiasedRisk, and has been introduced to this setting by xie2012SUREEstimatesHeteroscedastica.} The risk estimate is defined as $\mathrm{URE}(\mu, \Lambda) = \frac{1}{J}\sum_{j}\mathrm{URE}_{j}(\mu_{j}, \Lambda)$, where

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)) and ((ref)) 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.

Optimality of the URE estimators

I now show that the URE estimators defined in Section (ref) 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), 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

then the URE estimator has asymptotic risk as good as the oracle.

Therefore, establishing uniform consistency in ((ref)) 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) 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

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}$.

In the case where $\Sigma_{j}$ is diagonal for all $j$, Assumption (ref) (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) is enough to ensure uniform convergence of the first term of (ref).

To show convergence of the second term of ((ref)), 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

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) (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.

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). We now state the optimality of $\hat{\theta}^{\mathrm{URE}, \mathrm{m}}$ and $\hat{\theta}^{\mathrm{URE}, \mathrm{g}}$.

The proof of this theorem, provided in Appendix (ref), 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) 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), 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 robbins1964EmpiricalBayesApproach\footnote{The estimator obtains the Bayes risk of the model ((ref)) and ((ref)).} 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 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., jiang2009GeneralMaximumLikelihood, 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 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.

Forecasting $\theta_{T+1}$

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)) with $\mu_{j} = 0$. Write $\Lambda$ and $\Sigma_{j}$ in the following block matrix form:

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

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,

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

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

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}$.

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}$,

By the same argument made by Lemma (ref), if

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) that accounts for the fact that both $\theta_{j}$ and $\Sigma_{j}$ are now considered to be random draws.

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.

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 (ref) 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.

An application to teacher value-added

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.

Baseline model and data

I use a standard TVA model specified as the linear panel data model introduced in (ref), 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 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 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.

Estimation results and policy exercise

Figure (ref) 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) 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, 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.