Large-Scale Estimation under Unknown Heteroskedasticity

-2cm Large-Scale Estimation under Unknown Heteroskedasticity

KEY\ WORDS: Empirical Bayes, Shrinkage Estimation, Unknown Heteroskedasticity, Nonparametric Estimation.

Introduction

The central task of this paper is the large-scale, optimal estimation of unit-specific parameters under unknown heteroskedasticity. Motivating examples\footnote{Examples of the following three applications are Chetty2018, Gilraine2020, and Liu2020 respectively.} of this include estimation of intergenerational mobilities of individual neighborhoods that are used to inform families on the areas with highest potentials for upward mobility, or estimation of teacher value-added that is necessary as inputs in high-stakes labor decisions. We work within the following model in which, or variations of it, the estimation problem of these applications are framed\footnote{We focus on (ref) for common $J_i=J$ and discuss variations such as heterogeneous $J_i$ subsequently in Section (ref). This is because the central feature of our analysis is the unknown heteroskedasticity.}:

for $j=1,\dots,J$ and the unobserved $(\mu_i,\sigma_i) \sim_{\text{iid}} G_0$ for $i=1,\dots,n$. Subsequently we also discuss settings where the point estimation is nested into more complex problems, such as optimal forecasting for a large collection of short time series, and demonstrate our results' relevance there. For now, let us consider the teacher value-added literature as a concrete example: the outcome $y_{ij}$ represents the test score of the $j^{th}$ student in teacher $i$'s class, often after partialling out student-specific controls such that the independence of $\epsilon_{ij}$ is reasonable. A common objective is to estimate $\left\{\mu_i\right\}_{i\in[n]}$, interpretable as how the average student (i.e., $\epsilon_{ij}=0$) will perform under each teacher and commonly termed as the teachers' mean value-added.

Relative to the existing large-scale estimation methods that impose homoskedasticity, we develop optimal estimators that allow for unknown heteroskedasticity. This is important in at least two respects. First, $\sigma_i$ is a key component of the optimal estimator of $\mu_i$, and restricting $\sigma_i=\sigma$ may lead to estimators with large risk and in some cases systematic under- or over-estimation of certain subsets of $\{\mu_i\}_{i\in[n]}$. Furthermore, $\sigma_i$ is important for the estimation of unit-specific quantiles

that we shall argue are of policy-relevance, and a generalization of the commonly targeted mean value-added $\mu_i$. For example, in the context of the teacher value-added literature, this means that we estimate and compare how the bottom (say) 10% student will perform under each teacher instead of focusing only on the average student. This will be a more relevant object than the means if policymakers are concerned about the lower tails of student performances under each teacher.

In many applications $J \ll n$, which carries two implications for the large-scale estimation. First, the small number of observations per unit implies usage of only unit $i$'s observations to estimate $\mu_i$ or $q_{\alpha,i}$ leads to imprecise estimates and therefore high risk. Second, the large $n$ allows us to borrow information across units to sharpen individual estimates which can reduce the average risk. More precisely, it is known that the estimator of $\mu_i$ or $q_{\alpha,i}$ minimizing the average (across $i$) risk under quadratic loss is its posterior mean that takes $G_0$, the distribution of unit-specific parameters $(\mu_i,\sigma_i)$, as the prior. This forms the oracle estimator\footnote{See the end of this Section for notational details.} $\mathbb{E}_{G_0}[\, \cdot \mid \mathcal{Y}_i ]$, where $\mathcal{Y}_i$ is the set of observations associated with unit $i$.

The oracle is infeasible because $G_0$ is unknown in practice, and feasible empirical Bayes versions that exploit the cross-sectional information to mimic the oracle $\mathbb{E}_{G_0}[\, \cdot \mid \mathcal{Y}_i ]$ have to be proposed. However, nonparametric estimation of $G_0$ using the cross-section is known to have slow rates of convergence because we are essentially estimating a density of unobservables. This suggests that emulating the oracle will be a difficult task. In response to this, this paper proves an extension of the well-known Tweedie's formula\footnote{We review this in Section (ref).} under homoskedasticity to the setting of unknown heteroskedasticity, where we show that the oracle depends on $G_0$ only through the density of some sufficient statistics under model (ref). Thus the problem of emulating the oracle is dramatically simpler, since the performance of feasible versions now only hinge on their ability to estimate a density of observables. We further exploit this insight to provide regret bounds for the proposed estimators that are of the order $(\log n)^\kappa /n$ for some $\kappa>0$. The upsot is that there is little cost to practitioners adopting a nonparametric empirical Bayes approach, whereas the benefits in terms of risk reduction is often substantial relative to the commonly used parametric modeling of $G_0$.

The remaining paper is organized as follows. Section (ref) reviews the related literature. Section (ref) defines the model, risk and oracle estimators. Section (ref) proposes feasible versions of the oracles while Section (ref) provides theoretical guarantees for the feasible versions. Section (ref) discusses extensions to the model, including heterogeneous $J_i$ as well as constructing optimal forecasts in large-$n$ short-$T$ panels. Section (ref) employs numerical experiments to evaluate the finite-sample performances of the proposed estimators vis-\`a-vis alternative estimators. Section (ref) demonstrates the estimation methodology and objective through an application to estimation of teacher quality using a matched student-teacher dataset. Section (ref) concludes.

\noindentNotations. The symbol $[n]$ refers to the set $\{1,\dots,n\}$. The symbol $\mathcal{Y}_i$ refers to the set of observations associated with unit $i$ (i.e., $\{y_{ij}\}_{j\in[J]}$) and $\mathcal{Y}$ is taken to refer to the full set of observations (i.e., $\{y_{ij}\}_{i\in[n],j\in[J]}$), and finally $\mathcal{Y}^{(i)}$ refers to all observations except for $i$'s: $\{y_{lj}\}_{l\neq i,j\in[J]}$. A subscript on an expectation operator denotes objects that are conditioned upon. For instance, both $\mathbb{E}_{G}[\mu \mid \mathcal{Y}_i]$ and $\mathbb{E}[\mu \mid G,\mathcal{Y}_i]$ mean the same thing: that the integration is performed conditional on $\mathcal{Y}_i$ and $G$. We may sometimes also use a superscript on an expectation operator to make explicit the objects that are integrated out. For instance, $\mathbb{E}_{G}^{\mu}[\mu \mid \mathcal{Y}_i]$ emphasizes that we are integrating out the unknown parameter $\mu$.

Related Work

Our paper is within the empirical Bayes methodology, which goes back at least to Efron1973's interpretation of the James-Stein estimator as the posterior mean under a normal-normal hierarchical model with the normal prior estimated via method-of-moments. The James-Stein estimator can therefore be thought of as a parametric empirical Bayes methodology, where a particular class of prior is proposed, and feasible versions are proposed that seek to emulate the optimal estimator within the class of Bayes estimates implied by the class of priors. Recent papers in the parametric empirical Bayes methodology include Xie2012, and Kwon2023b, and the methodology has seen applications in work such as Chetty2014, Finkelstein2016.

Parametric empirical Bayes is however restrictive in the sense that the choice of a prior induces sometimes undesirable behavior in the consequent estimator. For instance, subject to the same number of observations, the posterior mean under a normal-normal hierarchical model imposes equal shrinkage for all units. Gilraine2020 studies the shortcomings of the parametric approach within an application of estimating teacher value-added under homoskedasticity and proposes nonparametric empirical Bayes estimators, a direction that we also take in this paper. Nonparametric empirical Bayes is more ambitious in that it targets the posterior mean without imposing a parametric form on the prior, and affords the flexbility that the parametric approach lacks. However, the nonparametric approach raises finite-sample concerns because $\mathbb{E}_{G_0}[\mu\mid \mathcal{Y}_i]$ depends on $G_0$, which suggests solving a deconvolution problem that is known to have slow rates of convergence (as slow as $(\log n)^{-1}$, see Fan1991).

In the restricted case of homoskedasticity $\sigma_i=\sigma$ for all $i$, theoretical progress was made possible by the Tweedie's representation of the posterior mean\footnote{Note that in the homoskedastic case, $G_0$ is supposed to represent the univariate distribution of $\mu_i$.} (see Efron2011):

where $y_i$ is the sample mean and $f_{G_0}(y_i \mid \sigma)$ is the (mixture) density of $y_i$. The insight of Tweedie's formula is that $\mathbb{E}_{G_0}[\mu\mid \sigma,\mathcal{Y}_i]$ comprises of the MLE $y_i$ and a Bayes correction term, for which we only require $f_{G_0}(y_i \mid \sigma)$ as opposed to the prior $G_0$.

This insight was instrumental for Jiang2009, which uses a certain nonparametric plug-in estimate $\hat{G}_0$, to establish the fast regret\footnote{Relative to the oracle that knows $G_0$.} convergence rates of $\mathbb{E}_{\hat{G}_0}[\mu\mid \sigma,\mathcal{Y}_i]$ in estimating $\{\mu_i\}_{i\in[n]}$. Their results were subsequently extended by Jiang2020 to the case of known heteroskasticity within a class of prior where the heteroskedasticity is assumed independent of $\mu_i$. More recent work by Chen2024 allowed for dependence between $\mu_i$ and the known heteroskedasticity, a feature that was demonstrated to be economically relevant, but where the dependence operates through the first two moments.\footnote{Remark (ref) elaborates in more precise terms the setting and contributions of our paper relative to these two papers. }

Our paper also adopts the nonparametric empirical Bayes approach, but is situated in a more general framework that allows for unknown heteroskedasticity. This substantially complicates the estimation of $\mu_i$ since $\sigma_i$, which regulates the Bayes correction in (ref), is now unknown and has to be estimated. Furthermore, the aforementioned insight from (ref) is no longer applicable since we do not know the form of $\mathbb{E}_{G_0}^\sigma [\sigma^2\frac{\partial}{\partial y}\log f_{G_0}(y_i \mid \sigma) \mid \mathcal{Y}_i]$ and in particular if it depends on $G_0$ solely through a density of observables. In this regard, the papers that adopt a similar framework are Gu2017, Gu2017a and Banerjee2023. Gu2017 and Gu2017a study the properties of $\mathbb{E}_{\hat{G}}[\mu\mid \mathcal{Y}_i]$ and $\mathbb{E}_{\hat{G}}[\sigma^2\mid\mathcal{Y}_i]$ for a nonparametric estimate $\hat{G}$ using simulations that show promising results, but do not provide theoretical guarantees. In contrast, we provide a significant extension of the Tweedie's formula that accommodates unknown heteroskedasticity and use it as the basis to study the theoretical properties of $\mathbb{E}_{\hat{G}}[\mu | \mathcal{Y}_i]$ for a nonparametric $\hat{G}$. Crucially, we do not restrict the form of dependence within $(\mu_i,\sigma_i)$.

Banerjee2023 also studies the optimal estimation of $\left\{\mu_i\right\}_{i=1}^n$ under unknown heteroskedasticity, albeit under a precision-weighted loss function. The precision-weighted loss function is easier to handle in that, roughly speaking, it transforms the heteroskedastic problem back into a homoskedastic one. However, the usage of the precision-weighted loss function is less appropriate for the types of applications we have in mind, where the estimation loss from each unit should be taken as equally important.

Finally, in addition to the estimation of $\{\mu_i\}_{i\in[n]}$, we also study the theoretical properties of similarly defined estimators for alternative estimation objectives such as unit-specific quantiles $\{q_{\alpha,i}\}_{i\in[n]}$ for some $\alpha\in(0,1)$ that we argue are a policy-relevant generalization of $\{\mu_i\}_{i\in[n]}$, and unit-specific variances $\{ \sigma_i^2\}_{i\in[n]}$ or standard deviations $\{\sigma_i\}_{i\in[n]}$ that may be relevant in other applications.

Model and Risk

Model

We maintain the following model throughout the paper:

for $j=1,\dots,J$ and $(\mu_i,\sigma_i) \sim_{\text{iid}} G_0$ for $i=1,\dots,n$, and discuss extensions to the model in Section (ref). We now briefly discuss how matched datasets such as the matched student-teacher dataset of our empirical application fits into (ref). Suppose that there are $\bar{J}$ students in total, with each student $j$ represented by $\tilde{\epsilon}_j$ drawn from a common distribution, who are then matched to the $n$ teachers. Let $i(j): [\bar{J}] \mapsto [n]$ be a function that returns the teacher identity of student $j$ resulting from this matching process. The outcomes of the students are then realized after the matching process as $\tilde{y}_j = \mu_{i(j)} + \sigma_{i(j)} \tilde{\epsilon}_j$. This fits into our framework (ref) by defining $J := | \{j\in[\bar{J}]:i(j)=i\} | $, and $\{\epsilon_{ij}\}_{j\in[J]} := \{\tilde{\epsilon}_{j}\}_{i(j)=i}$ and similarly for the $\{y_{ij}\}_{j\in[J]}$. The conditional independence assumption in (ref) then implies that the matching process between students and teachers did not depend on $\{\tilde{\epsilon}_{j}\}_{j=1}^{\bar{J}}$. More realistically, the variable $\tilde{y}_{j}$ may be taken to have a set of relevant covariates partialled out from the raw test scores $\tilde{y}_j^*$:

For example, the variables $x_{j}$ in the teacher-value added literature typically include student-specific covariates such as student demographics or lagged test scores such that the independence assumption in (ref) is reasonable. We abstract from the estimation of the common vector $\beta_0$ and take it to be known as is common in the literature, see e.g., Gilraine2020 and Kwon2023b.

An essential difference between (ref) and models commonly studied in the empirical Bayes literature is the heterogeneous and unknown $\{\sigma_i^2\}_{i\in[n]}$. This is not only empirically more plausible, but also allows for meaningful comparisons of units beyond their mean parameters. We provide two areas where the unit-specific quantile parameters will be of relevance.

Therefore, in the subsequent theoretical analysis we also study the estimation problems involving $\{q_{\alpha,i}\}_{i\in[n]}$, and even $\{\sigma_i\}_{i\in[n]}$ or $\{\sigma_i^2\}_{i\in[n]}$ that may be of interest in other areas. Finally, in anticipation of the theoretical results, we provide several remarks relevant to the existing nonparametric empirical Bayes literature.

Risk

The discussion in this section will be framed in terms of the mean value-added $\{\mu_i\}_{i\in[n]}$, because these are the de-facto object of interest within the empirical Bayes literature. Nonetheless, our theoretical results equally apply to other objects of interests including $\{q_{\alpha,i}\}_{i\in[n]}$ as we show in Section (ref).

An estimator $\hat{\mu}$, which maps the observations $\mathcal{Y}$ to $\mathbb{R}^n$, is evaluated under the quadratic loss function

Define the (integrated) risk

The independence across $i$ together with the Rao-Blackwell Theorem\footnote{See, e.g., Theorem 7.8 of Casella1998.} imply the integrated risk is minimized, over all borel mappings of the data $\mathcal{Y}$ to $\mathbb{R}^n$, by

and

are the sufficient statistics for $(\mu_i,\sigma_i^2)$. Note that similar statements also hold for the decision problems involving $\{\sigma_i\}_{i\in[n]}$ or $\{\sigma_i^2\}_{i\in[n]}$, and we denote their corresponding optimal estimators as $\hat{\sigma}^{*}$ and $\hat{\sigma}^{2*}$.

Finally, even though our interest in $\hat{\mu}^*$ stems from the decision problem (ref), we note that there are other reasons why $\hat{\mu}^*$, and thus our proposed estimator, may be of interest. We provide one such setting where $\hat{\mu}^*$ is also essential. This is a decision problem where we wish to identify units with $\mu_i$ below a certain threshold $c$ under an absolute loss function.

Thus in the context of teacher value-added literature, $\hat{\mu}^*$ will be of interest when we wish to identify teachers with poor performance defined in terms of an absolute threshold, e.g., mean value-added below a certain fixed value, and where the cost of misidentification increases and is symmetric in terms of the absolute loss. This result is complementary to existing results, such as in Gu2023 that seek to identify units within a relative threshold (say bottom $10\%$ of the population), or Kline2024 that focus on assigning ranks to units in the population under alternative loss functions.

Estimation

Characterizing Optimal Estimators

As a first step to studying feasible estimators, we prove a Tweedie-type characterization of the oracles.

Theorem (ref) states that the posterior mean of the model parameters such as $\mu$ depend on $G_0$ only through the mixture density $f_{G_0}(y,s^2)$, which is a density of observables. The representation of $\mathbb{E}_{G_0} [\mu \mid y,s]$ looks similar to its counterpart under a model where $\sigma^2$ is known:

In particular Theorem (ref)(i) tells us that relative to the case of known heteroskedasticity, the posterior mean of $\mu$ under unknown heteroskedasticity involves replacing the unknown $\sigma$ with its posterior mean, and an additional Bayes correction term $\frac{1}{J}\frac{\partial}{\partial y} \mathbb{E}_{G_0}[\sigma^2 \mid y,s]$.

The characterization of $\mathbb{E}_{G_0} [\sigma \mid y,s]$ in Theorem (ref) will be relevant in the estimation of teacher-specific quantiles (ref), which in our model (ref) takes a particularly simple form:

where $\Phi(\cdot)$ is the standard normal CDF. In particular, the quantiles are linear functions of the mean $\mu$ and standard deviation $\sigma$. Theorem (ref) then provides us with a Tweedie-type characterization of the optimal estimator of $q_{\alpha}$,

where the minimization is over all borel mappings of the observations $\mathcal{Y}$ to $\mathbb{R}^n$.

Theorem (ref) and Lemma (ref) will be used in the following sections to both propose and study the theoretical properties of feasible counterparts to the oracles.

Proposed Estimators

We first define truncated forms of the optimal estimators\footnote{ We use a slightly different representation of $\hat{\mu}^*$ that is more convenient for theoretical analysis -- see second to last line of (ref) in Appendix (ref) for derivation. } given in Theorem (ref) and Lemma (ref):

Note that when $G=G_0$ and $\rho=0$, we recover $(\hat{\mu}^*,\hat{\sigma}^{2*},\hat{\sigma}^*,\hat{q}_\alpha^*)$ respectively.

In operationalizing the estimators, we replace the unknown $f_{G_0}$ with a leave-one-out (approximate) nonparametric maximum likelihood estimator (NPMLE). That is, our estimator for $\mu_i$ is $ \hat{\mu}\big( y_i,s_i \mid \hat{G}^{(i)},\rho \big) $ where $\hat{G}^{(i)}\in\mathcal{G}$, the set of all bivariate distributions on $\mathbb{R}\times\mathbb{R}_{++}$, and

where $\eta_n \asymp \tfrac{1}{n}$. It is well-known by Lindsay1983 that the exact NPMLE exists and takes the form of a discrete distribution with at most $n-1$ support points.

We now explain our choice of a leave-one-out NPMLE $f_{\hat{G}^{(i)}}$ as opposed to a common NPMLE $f_{\hat{G}}$ that uses all observations $\mathcal{Y}$ for all $i\in[n]$:

For one, in terms of finite-sample performances, we show using numerical experiments in Section (ref) that the leave-one-out $f_{\hat{G}^{(i)}}$ can often lead to a finite-sample improvement over $f_{\hat{G}}$, intuitively by reducing overfitting in small samples. In terms of the theoretical analysis, the leave-one-out estimator allows us to side-step the analysis of certain entropy bounds for the optimal estimators. This is convenient in analyzing the extended Tweedie's formula of Theorem (ref) that is substantially more complex than in the case of homoskedasticity as in (ref). This is especially so because we study the theoretical properties of estimators for multiple estimation objectives. Nonetheless, we have constructed the proofs in a modular fashion such that pieces of the proofs can be incorporated into a theoretical study of estimators using a common $f_{\hat{G}}$ if desired.

A potential shortcoming of the leave-one-out NPMLE is that the maximization problem has to be solved once for every $i\in[n]$, which may turn out to be computationally expensive for large $n$. We view this as a minor obstacle because, firstly, the leave-one-out NPMLE is directly parallelizable and this alone alleviates the associated time cost. Furthermore, as we do in the simulations and empirical application, one may first compute the full NPMLE $f_{\hat{G}}$ and use $\hat{G}$ as the initial guess downstream for the subsequent leave-one-out estimates $f_{G^{(i)}}$. Assuming that the maximum only changes slightly from leaving out each observation, as is observed in the simulations, the additional cost of the (parallelized) second step will be marginal.

Theoretical Results

The set-up for our theoretical analysis is a triangular array, where the subscript $n$ indicates the row of the array:

for $j=1,\dots,J$ and $(\mu_{n,i},\sigma_{n,i}) \sim_{iid} G_{0,n}$ for $i=1,\dots,n$.

In preparation for the main convergence result, we clarify certain notations and state the assumptions underlying the result. First, the notation $G(A \times B)$ for a bivariate distribution $G$ and sets $A,B\subset\mathbb{R}$ denotes the probability that $G$ places on the event $(\mu\in A,\sigma\in B)$. Second, $(\hat{\mu},\hat{\sigma},\hat{\sigma}^2,\hat{q}_{\alpha})$ for the remainder of this section are taken to denote the estimator forms of (ref) with $G$ replaced by the leave-one-out estimator $\hat{G}^{(i)}$ and $\rho$ by $\rho_n$, both more clearly explicated in Assumptions (ref) and (ref). Finally, the notation $X\lesssim Y$ means that there exists a constant $C$, possibly dependent on $(\gamma_1,\gamma_2,\underline{\sigma},L)$ defined below, such that $X \leq C Y $.

Assumption (ref) essentially requires that the support of $G_{0,n}$ grows at at most logarithmic rates, and can conceivably be relaxed to allow for sub-Gaussian tails, which would then accommodate distributions with potentially infinite support such as Gaussian mixtures for $\mu_{i,n}$ and Gamma mixtures (suitably truncated away from 0) for $\sigma_{i,n}^2$. What is crucial for our proofs however, is the lower-bound $\underline{\sigma}$ on the support of $\sigma$, which is used to guarantee sufficient smoothness of the mixture density $f_{G_{0,n}}$.

Assumption (ref) requires that there exist approximate NPMLEs with supports that lie within the boundaries of the data. Note that it also requires knowledge of a lower bound $\underline{\sigma}$ on $\{\sigma_{n,i}\}_{i\in[n]}$. The final assumption is on the truncation parameter $\rho_n$.

We now state the main theoretical result of the paper.

Theorem (ref) provides regret bounds for the estimation of $\{\mu_i\}_{i\in[n]}$, $\{\sigma_i\}_{i\in[n]}$ and $\{\sigma_i^2\}_{i\in[n]}$, which all converge at the parametric rate (up to logarithmic factors). This suggests that the cost of taking a nonparametric approach will be small in finite-samples. In contrast, the benefits in terms of risk reduction can be large relative to assuming a particular parametric form for $G_{0,n}$, as documented in Gilraine2020. Similar results are well-known in the cases of homoskedasticity or known heteroskedasticity, and Theorem (ref) can thus be seen as an extension to the setting where the heteroskedasticity is unknown and has to be estimated.

Theorem (ref) also extends ratio-optimality type results, originally established by Jiang2009 for homoskedastic models\footnote{ These results were similarly established by Brown2009 and Liu2020 for $f$-modeling estimators. }, to accommodate unknown heteroskedasticity and more general estimation objectives. For example, supposing

such that the oracle risk is not too small, then Theorem (ref) implies the following ratio-optimality result:

Note that the concept of ratio-optimality, which targets the relative risk, is a much harder concept than the usual empirical Bayes optimality which only require that

Finally recall that by Lemma (ref)$, \hat{q}_{\alpha,i}^*$ (i.e., the oracle estimator of $q_{\alpha,i}$) is a linear function of $\hat{\mu}_i^*$ and $\hat{\sigma}_i^*$. Theorem (ref) therefore also provides bounds on the regret of using $\{\hat{q}_{\alpha,i}\}_{i\in[n]}$ to estimate $\{q_{\alpha,i}\}_{i\in[n]}$.

Extensions

This section considers three extensions of the basic model (ref), including to the literature on optimal forecasting in large-$n$ short-$T$ panels, and we focus on the estimation of $\left\{\mu_i\right\}_{i\in[n]}$ to keep the discussion succinct.

Heterogeneous Sample Sizes. The first extension consists of settings where the number of observations per unit, $J_i\in \mathbb{N}$, is heterogeneous and furthermore informative about the underlying parameters $(\mu_i,\sigma_i)$. For instance, teachers have differing class sizes and a teacher with higher value-added $\mu_i$ may be systemically assigned smaller classes. Failure to take this into account will result in under- or over-estimation of $\left\{\mu_i\right\}_{i=1}^n$; see Chen2024 for fuller arguments in the setting of estimating neighborhood effects under homoskedasticity or known heteroskedasticity. In this case, we augment the model to be

for $j=1,\dots,J_i$ and $(\mu_i,\sigma_i,J_i) \sim_{\text{iid}} H_0$ for $i=1,\dots,n$, and modify the decision problem to

where $\mathcal{J}:=\{J_i\}_{i\in[n]}$ and the minimization is over all borel mappings of $(\mathcal{Y},\mathcal{J})$ to $\mathbb{R}^n$. The oracle estimator $\check{\mu}^*$ takes a similar form as in Theorem (ref) but with conditioning on $J_i$ everywhere:

Thus, when the heterogeneity in the discrete $J_i$ is limited such that there are sufficient observations for each value $J$, one may pool observations of teachers with identical (or similar) $J_i$ to estimate the density $f_{H_0}(y,s | J_i=J)$ for each $J$. The fast regret convergence rates provided in Theorem (ref) suggest that this can still work reasonably well in practice.

Informative Covariates. The discussion above applies equally to the empirically relevant setting where the researcher has access to covariate(s) $x_i\in\mathbb{R}^m$ that is potentially informative about the parameters $(\mu_i,\sigma_i)$ and that may help sharpen estimation. For instance in the teacher value-added literature, the years of experience of teacher $i$ is likely to be positively correlated with her value-added $\mu_i$. In such settings, by suitably modifying the model and decision problem as above, we may once again show that the oracles depend only on the distribution of $(\mu_i,\sigma_i,x_i)$ through the conditional density $f(y,s^2 \mid x)$.

Dynamic Panel Forecasting. Consider the problem of constructing optimal forecasts for a large collection of short time series, taking the basic dynamic panel model studied in Liu2020 as an example:

for $i=1,\dots,n$ and $t=1,\dots,T$. The motivating application is the stress testing of bank-holding companies as part of regulatory requirements, where the key step is being able to first reliably forecast bank balance-sheet variables under observed macroeconomic and finanacial conditions. Furthermore, as detailed in Liu2020, frequent mergers in the banking industry and changes in regulatory environments post-2008 makes a large-$n$ small-$T$ framework suitable to this task. In such a model, the optimal one-step-ahead forecast is $\mathbb{E}_{G_0}[y_{it} \mid \mathcal{Y}_{i,t-1}] = \rho y_{it-1} + \mathbb{E}_{G_0}[\mu_i \mid \mathcal{Y}_{i,t-1}]$, where $\mathcal{Y}_{i,t-1}$ is the vector of observations up to time period $t-1$ associated with unit $i$. The conditional mean $\mathbb{E}_{G_0}[\mu_i \mid \mathcal{Y}_{i,t-1}]$ takes the form (ref) with $f_{G_0}(y_i | \sigma)$ replaced by $f_{G_0}(y_i| \sigma,y_{i,0})$ where $y_i$ is the sufficient statistic for $\mu_i$, which reveals that the oracle forecasts' reliability hinges upon the homoskedasticity assumption, because $\sigma^2$ regulates the Bayes correction in estimating each and every $\mu_i$.

Within this context, our model adds an additional layer of heterogeneity through the unit-level varying $\sigma_i^2$:

for $i=1,\dots,n$ and $t=1,\dots,T$. Indeed it is reasonable to expect in such an application, that the heterogeneity in $\{\sigma_i^2\}_{i\in[n]}$ to be comparable or larger than is present within $\{\mu_i\}_{i\in[n]}$. As a result, the oracle forecasts that allow for unknown heteroskedasticity will likely have substantial improvements over those that impose homoskedasticity by allowing for unit-specific Bayes correction. In such a model, emulating the oracles again reduces to emulating $\mathbb{E}_{G_0}[\mu_i \mid \mathcal{Y}_{i,t-1}]$, which now takes the form in Theorem (ref) with $f_{G_0}(y_i,s_i^2)$ replaced by $f_{G_0}(y_i,s_i^2| y_{i,0})$ and $s_i^2$ is the sufficient statistic for $\sigma_i^2$. If one is willing to accept independence between $y_{i,0}$ and $(\mu_i,\sigma_i^2)$, then the theoretical results in Section (ref) directly delivers fast regret\footnote{Regret that is relative to the oracle forecasts, i.e., that knows the distribution $G_0$.} convergence rates of using feasible forecasts that plug in the leave-one-out NPMLE of the density of $(\mu_i,\sigma_i^2)$. In the case of unrestricted dependence between $y_{i,0}$ and $(\mu_i,\sigma_i^2)$, the proof and results developed here, leveraging on the expressions derived in Theorem (ref), may be extended under e.g., suitable smoothness conditions on the density of $y_i,s_i^2 \mid y_{i,0}$.

Simulations

We now utilize numerical experiments to evaluate the finite-sample performance of the proposed estimators and alternative sets of estimators.

Estimators

In implementing the proposed estimators of Section (ref), we adopt a leave-one-out $n^{1/2}$-order sieve MLE $f_{\hat{G}^{(i)}(n^{1/2})}$-- see Remark (ref). We call this set of estimators het.

As alternatives, we include an additional set of estimators using the (full) sieve-MLE $f_{\hat{G}(n^{1/2})}$; again, see Remark (ref). We call this $\textsc{het}_{\text{full}}$. This is included to show the finite-sample improvements from the leave-one-out strategy within the context of nonparametric empirical Bayes estimation. We also include a homoskedastic version of het, i.e., imposing homoskedasticity $\sigma_i^2 = \sigma^2$ and where $\sigma^2$ is estimated using $\tfrac{1}{n}\textstyle\sum_{i=1}^ns_i^2$. We call this $\textsc{hom}$. Finally the set of naive estimators, which uses $(y_i,s_i^2)$ to directly estimate $(\mu_i,\sigma_i^2)$ and $y_i + \Phi^{-1}(\alpha)\cdot s_i$ to estimate $q_{i,\alpha}$, is also included and referenced as naive.

Data Generating Process

The DGP for the simulations proceeds in two stages, where conditional on $\left\{(\mu_i,\sigma_i)\right\}_{i=1}^n$ the second stage generates observations $y_{ij}$ according to model (ref) with $J=15$.

In the first stage, we adopt a relatively simple DGP and generate $(\mu_i,\sigma_i)$ by coupling Gaussian and Gamma random variables\footnote{$\gamma(\kappa,\lambda)$ denotes the Gamma distribution with shape $k$ and scale $\lambda$, and $F_{\gamma(\kappa,\lambda)}^{-1}(\cdot)$ denotes the quantile function of $\gamma(\kappa,\lambda)$. $\Phi(\cdot)$ denotes the standard normal CDF.}:

Thus $\mu_i$ is a Gaussian variable with mean $\alpha$ and variance $\nu$, whereas $\sigma_i^2$ is a Gamma variable with shape $\kappa$ and scale $\lambda$. The parameter $\rho$ controls the correlation between $(\mu_i,\sigma_i^2)$. In the calibrated DGP, we set $(\alpha,\nu,\kappa,\lambda,\rho)$ such that the first five moments of $(\mu_i,\sigma_i^2)$ matches those estimated from the dataset of our teacher value-added application:

Figure (ref) displays the marginal and joint densities of $(\mu_i,\sigma_i^2)$ under this calibration. We also perturb this calibrated DGP along a single dimension (e.g., changing $\mathbb{V}[\mu_i]$ while keeping the other four moments constant) and observe the performances of the estimators under each perturbation.

Results

Relative Regrets I

Table (ref) displays\footnote{Appendix (ref) explores further numerical experiments, when the estimation objective is $\{\sigma_i\}_{i\in[n]}$ and $\{\sigma_i^2\}_{i\in[n]}$ instead.} the relative regrets involving the MSE loss $l(\cdot,\cdot)$ as defined in (ref):

for varying estimators $\hat{\mu}$ or $\hat{q}_{0.1}$ and DGPs, and where the expectation $\hat{\mathbb{E}}$ is taken across the simulation rounds. The 3\textsuperscript{rd} column displays results for the calibrated DGP, whereas the DGP for the 4\textsuperscript{th} column onward each perturb the calibrated DGP along one dimension. For example, the 4\textsuperscript{th} column displays results for the DGP with the same moments as the calibrated DGP, except that now $\mathbb{V}[\sigma_i^2] = 0.018$ so that $\sigma_i^2$ is relatively hetergeneous. The 5\textsuperscript{th} column displays results for the DGP with the same moments as calibrated DGP, except that now $\mathbb{V}[\mu_i] = 0.0023$ so that $\mu_i$ is relatively homogeneous. The top and bottom half of Table (ref) displays the relative-regret when the estimation target is $\{\mu_i\}_{i\in[n]}$ and $\{q_{0.1,i}\}_{i\in[n]}$ respectively.

We first observe that the calibrated DGP is one where $\{\sigma_i^2\}_{i\in[n]}$ is relatively homogeneous as compared to $\{\mu_i\}_{i\in[n]}$. As a result, the hom estimator applies a nearly optimal Bayes correction (see $(i)$ of Theorem (ref) and (ref)) for the estimation of $\{\mu_i\}_{i\in[n]}$ and is competitive with het in terms of relative regret. As we shall see later however, this similarity of hom and het in terms of average performance in MSE across units masks important differences between the two methodologies in terms of detecting units with low $\mu_i$, an exercise of wide policy-relevance. Nonetheless, for large $n$ in 3\textsuperscript{rd} column, the relative regret of het disappears whereas hom still suffers a non-zero relative regret because of the (limited) heterogeneity in $\{\sigma_i^2\}_{i\in[n]}$.

Furthermore, when the estimation target shifts to $\{q_{0.1,i}\}_{i\in[n]}$, the performance of hom deteriorates relative to het. This is because $\sigma_i^2$ not only matters for the Bayes correction in estimating $\mu_i$, but also directly as part of the estimation target: $q_{0.1,i}=\mu_i+\Phi^{-1}(0.1)\cdot\sigma_i$ where $\Phi^{-1}(0.1)=-1.28$. The suboptimality induced by hom's homoskedasticity assumption turns out to be large in this case.

The relative-regret differential between het and hom is also especially large in the case of heterogeneous $\{\sigma_i^2\}_{i\in[n]}$ (see 4\textsuperscript{th} column) or strong dependence between $\mu_i $ and $\sigma_i^2$ (see 6\textsuperscript{th} column). Recall that hom, which assumes homogeneous $\sigma_i=\sigma$, implicitly assumes no dependence between $(\mu_i,\sigma_i)$. When there is meaningful heterogeneity in $\{\sigma_i^2\}_{i\in[n]}$ and dependence between $(\mu_i,\sigma_i)$, hom is unable to take advantage of the information in $s_i^2$ regarding $\sigma_i^2$ to sharpen the estimation of $\mu_i$ and so incurs a larger relative regret. On the other hand, het is fully capable of exploiting this information for estimating $\{\mu_i\}_{i\in[n]}$ and its relative regret remains similar from 3\textsuperscript{rd} column to 4\textsuperscript{th} or 6\textsuperscript{th} columns.

The next observation is that as $\{\mu_i\}_{i\in[n]}$ becomes homogeneous (from 3\textsuperscript{rd} to 5\textsuperscript{th} column), the relative regrets of all estimators increase. This is because the oracle knows the distribution of the parameters and incurs a much smaller MSE in the case of homogeneous $\{\mu_i\}_{i\in[n]}$. The cost of estimating the distribution for the feasible estimators thus implies a higher regret relative to the oracle. Nonetheless, the fast convergence rates established in Section (ref) imply that het is able to achieve small relative regrets for realistic sample sizes.

Furthermore, it is within this setting (5\textsuperscript{th} column), where the relative regrets of the estimators are large, that the finite-sample gains from the leave-one-out strategy is most substantial. The het estimator improves on $\textsc{het}_{\text{full}}$ that does not apply the leave-one-out estimation by over 10 percentage points when the sample size is very small. The improvement fades away as $n$ becomes large, because then the complexity of the model relative to the number of observations (represented by the order of the sieve, $n^{1/2}$), decreases and overfitting no longer becomes an issue for $\textsc{het}_{\text{full}}$. We note however that if the order of the sieve were to be increased (e.g., to $n^{3/4}$ or even considering the actual NPMLE of sieve order $n$), the finite-sample improvements of het relative to $\textsc{het}_{\text{full}}$ will be larger for the same sample sizes as in Table (ref).

In sum, we conclude that het dominates hom and naive for both targets $\{\mu_i\}_{i\in[n]}$ and $\{q_{0.1,i}\}_{i\in[n]}$ across all specifications and in particular even for very small sample sizes. Thus, there is little finite-sample cost to adopting het in practice, and much to be gained in terms of the relative regrets especially when $\{\sigma_i^2\}_{i\in[n]}$ is hetergeneous or when there exists strong dependence between $\mu_i$ and $\sigma_i^2$.

Relative Regrets II

Figure (ref) displays the relative regrets involving the loss function $l_{\text{d}}(\cdot,\cdot,c)$ as defined in (ref):

for the left and the right panels of the figure respectively. Intuitively, the loss functions underlies a decision problem of identifying units with $\mu_i$ or $q_{0.1,i}$ below a fixed threshold $c$, where the loss from each unit increases symmetrically in the magnitude of misidentification. The DGP utilized is the calibrated DGP with $n=4000$, and the threshold value $c$ is varied on the $x$-axis of Figure (ref).

From the left panel, we observe that het dominates hom for almost all values of $c$. Furthermore, even though the heterogeneity in $\{\sigma_i^2\}_{i\in[n]}$ is limited under this DGP, the relative regret differential between het and hom widens considerably for small values of $c$. This is likely due to the negative dependence between $\mu_i$ and $\sigma_i^2$: units with the smallest $\mu_i$ tend to have the largest $\sigma_i^2$, and hom, which assumes a common $\sigma^2$ for all units, will suffer the largest regret at detecting units with small $\mu_i$ because it is precisely at these regions that the homoskedasticity assumption is most misplaced.

Interestingly, the story seems to be reversed for detecting units with small $q_{i,0.1}$, in that the regret differential decreases for smaller values of $c$. It turns out that there is a happy coincidence for hom in this case, specifically for small values of $\mu_i$: the (negative) bias in the estimation of $\mu_i$ arising from the homoskedasticity assumption cancels out the (negative) bias in the estimation of $\sigma_i^2$ through the form of the quantile: $q_{0.1,i}= \mu_i - 1.28\cdot \sigma_i$.\footnote{To elaborate in detail: in this DGP the units with small $\mu_i$ tend to have larger $\sigma_i$, and so the homoskedasticity assumption falsely imputes small $\sigma_i$s to these units; thus hom estimators of $\sigma_i$ for these units will be too small. As a result of this, the Bayes correction for the estimation of $\mu_i$ is smaller than optimal and the hom estimates of these $\mu_i$ will also be too small.}

Finally, the naive methodology, which simply plugs in $(y_i,s_i^2)$ for the unknown $(\mu_i,\sigma_i^2)$, performs poorly at all values of $c$. In fact, it has relative regret that is uniformly above $0.15$ in the right panel. This is because it identifies many units as having $\mu_i$ or $q_{0.1,i}$ below $c$, but the vast majority of these are likely due to noise in the sufficient statistics than a true reflection of the unit-specific parameters $(\mu_i,\sigma_i^2)$.

Empirical Application

We conclude with an application of estimating teacher quality using a matched student-teacher dataset from the North Carolina Education Research Data Center.

Data

The dataset is a panel dataset on elementary school students from grades 3--5 in North Carolina. Each observation is linked to a (student, year, grade, subject) cell with covariates including student demographics information such as ethnicity, gender, economic status, etc. We follow Gilraine2020 in restricting observations to students that are matched to teachers and have a lagged test score in the subject associated to that observation.

Estimation

To keep the analysis simple, we further restrict attention to observations from year 2019 and where the class size associated with each observation is between 14 to 22 students. This is to exclude special education classes and to abstract from potential effects of class size heterogeneity. Doing so leaves us with a total of 50,063 observations, each associated with a unique student, and a total of $n=2,841$ teachers in the sample. Table (ref) displays the summary statistics for some of the covariates present in the restricted dataset.

We take a standard model of teacher value-added:

where $y_{j}^*$ is student $j$'s 2019 reading test score, and $i(\cdot)$ is the matching function that returns student $j$'s teacher identity in 2019. The covariates $x_j$ include student $j$'s 2018 reading test score, ethnicity, gender, family economic status, whether the student is an english language learner, and age. The test scores are standardized to have mean 0 and standard deviation 1 within each (year,grade) cell. The identification assumption is that conditional on the covariates, the student effect $\epsilon_j$ is independent of her matched teacher's profile $\left(\mu_{i(j)},\sigma_{i(j)}\right)$. We follow the literature in first running a regression that purges the covariates' effects\footnote{That is, we estimate the common coefficients $\beta_0$ using within-teacher variation and move the consequent estimated covariates $x_j^\prime \hat{\beta}$ to the left-hand side, defining $y_j:=y_j^*-x_j^\prime\hat{\beta}$. The estimates $\hat{\beta}$ are consistent for $\beta_0$ because the relevant sample size is $n$, which is large.} which then leads to the model (ref) posited in Section (ref), written in terms of the matching function notation:

for $j=1,\dots,J_i$ and $i=1,\dots,n$, where $J_i:=\lvert\left\{j:i(j)=i\right\}\rvert$.

The difference with the standard model, such as in Gilraine2020, is that we allow for teacher-level heteroskedasticity. In a model where teachers differ in both $\mu_i$ and $\sigma_i$, there are multiple ways to define and measure teacher quality. For example, conditional on teachers having the same $\mu_i$, whether a high $\sigma_i$ is preferred depends on the policymaker's preferences. The most common way to measure a teacher's performance is by her mean value-added $\mu_i$, which measures how the average student ($\epsilon_j=0$) would have performed under each teacher. In what follows we also consider defining teacher performances by $q_{0.1,i}$, which measures how the 10\textsuperscript{th} percentile student would have performed under each teacher. This implies we penalize teachers with high $\sigma_i$, because $q_{0.1,i}=\mu_i-1.28\cdot\sigma_i$.

Estimators

With the model (ref) in hand, we construct the teacher-specific sufficient statistics

and estimate $(\mu_i,q_{0.1,i})$ using firstly the proposed set of estimators (called het), secondly a modification of het that assumes homoskedasticity (called hom), and thirdly naive that uses the teacher-specific sufficient statistics $(y_i,s_i^2)$ directly; see Section (ref) for more details. Note that we abstract from class size heterogeneity in the estimation which we justify by our sample restriction to a relatively narrow band of class sizes between 14 to 22 students.

Results

The main finding here is that teachers with small $\mu_i$ will also tend to have large $\sigma_i$, as is shown in the left panel of Figure (ref) that plots the het estimates of $\sigma_i$ against those of $\mu_i$. This implies differences in teacher quality is exacerbated, in that teachers who differ in how the average student performs under them (that is, $\mu_i$) tend to differ even more in how the lower-percentile students will perform under them (e.g., $q_{0.1,i}$). To see this, the right panel of Figure (ref) displays the test score distribution of (ref) for two representative teacher profiles $(\mu,\sigma)$ constructed from the het estimates; see the notes of Figure (ref) for construction details. The negative correlation between $\mu_i$ and $\sigma_i$ implies that while the two teachers' $\mu$ differ by 0.42 (i.e., 0.42 standard deviations of the raw test scores), their $q_{0.1}$ differ by $0.51$, a $22\%$ increase in difference.

We turn now to explore in detail the individual estimates for the two targets, $\{\mu_i\}_{i\in[n]}$ and $\{q_{0.1,i}\}_{i\in[n]}$, displayed in Figure (ref). To give an orientation to the panels of this figure: the top-left panel displays the scatter plot of the naive estimate of $\mu_i$ against the het estimate of the same $\mu_i$, for every $i$. Now with reference to the first row which compares het and naive estimates, we see that there is substantial shrinkage from using het for both targets. The shrinkage also appears to be non-linear, particularly for estimation of $\{\mu_i\}_{i\in[n]}$ which has less (more) shrinkage for smaller (larger) values of $y_i$. Turning to the second row of Figure (ref), it appears from first sight that the estimates for $\{\mu_i\}_{i\in[n]}$ are very similar across the het and \textsc{hom} methodologies. There are however still systematic differences between the two, particularly so when we consider the estimation of $\{q_{0.1,i}\}_{i\in[n]}$. The estimates from \textsc{het} are more spread out compared to those from \textsc{hom}, because \textsc{het} allows for heterogeneity in $\sigma_i$ whereas \textsc{hom} does not. Furthermore, it appears from the best-fit line that \textsc{hom} generally overpredicts and underpredicts $q_{0.1,i}$ relative to \textsc{het} near the left and right end points respectively. This may be explained by the fact that a small $\mu_i$ is associated with a large $\sigma_i$, which \textsc{hom} does not take into account and thus spuriously compresses the estimates of $\{q_{0.1,i}\}_{i\in[n]}$.

Finally, recall by Lemma (ref) that the oracle $\hat{\mu}^*$ (or $\hat{q}_{0.1}^*$) is also optimal for the problem of detecting teachers with $\mu_i$ (or $q_{0.1,i}$) below a fixed threshold. We thus explore the various methodologies in the context of this problem, again according to the two cases where teacher quality is measured by $\{\mu_i\}_{i\in[n]}$ and by $\{q_{0.1,i}\}_{i\in[n]}$. The left (and right) panel of Figure (ref) plots, by methodology, the number of teachers whose estimates of $\mu_i$ (and $q_{\alpha,i}$) are flagged as being less than the threshold value which is varied on the $x$-axis. Unsurprisingly, naive flags the greatest number of teachers for every threshold value. It is likely that most of these were flagged because of the large variability of the estimates $(y_i,s_i^2)$ rather than a true reflection of the value-added $\mu_i$.

Another notable feature from the left panel of Figure (ref) is that hom flags more teachers than het for almost all threshold values. For instance, it flags out 25 more teachers than het for the threshold value of $-0.24$, or an increase of 40% in relative terms. As we have seen from Figure (ref), teachers with smaller $\mu_i$ tend to have larger $\sigma_i$. The hom methodology, which assumes $\sigma_i=\sigma$ for all $i$, applies a less-than-optimal Bayes correction (relative to the oracle) for this set of teachers and the additional variability in the hom estimates results in more teachers being flagged. On the contrary, het allows for heterogeneous variances and dependence within $(\mu_i,\sigma_i)$ and therefore does not suffer from this problem. This difference may also be seen from the bottom-left panel of Figure (ref), where \textsc{hom} underpredicts relative to \textsc{het} in the left tail.

On the right panel where teacher quality is defined in terms of $\{q_{0.1,i}\}_{i\in[n]}$, hom now flags significantly fewer teachers than het for every threshold value. This is again due to the same combination of hom imposing a common $\sigma$ for all teachers and the negative correlation within $(\mu_i,\sigma_i)$, albeit through a different channel. Recall that $ q_{0.1,i} = \mu_i - 1.28\sigma_i $. As a result, hom understates how poorly the 10\textsuperscript{th}-percentile students will do under teachers with smaller $\mu_i$ whereas het, which incorporates the negative correlation in $(\mu_i,\sigma_i)$, does not. In absolute terms, the implications of both methodologies are very different. The hom methodology generally flags about 100 fewer teachers than does \textsc{het} as a result of its a-priori homoskedasticity restriction -- a restriction that Figure (ref) suggests to be misguided.

Conclusion

This paper contributes to the growing empirical Bayes (EB) literature by proposing nonparametric EB estimators with regret bounds of the order $\tfrac{1}{n}(\log n)^\kappa$. The upshot is that nonparametric EB methods come at little cost in typical sample sizes and in fact with much benefits in terms of risk reduction. Importantly, these results are developed within a framework that features unknown heteroskedasticity. Allowing for this feature is essential for two reasons: firstly, in the estimation of the commonly targeted mean parameters by applying the optimal level of unit-specific Bayes correction. This is important within many micreconometrics applications and extends to other areas such as the literature on optimal forecasting for a large collection of short time series. Unknown heteroskedasticity is also essential for comparing units along richer dimensions of heterogeneity, such as the performances of lower-percentile students under each teacher. These (i.e., unit-specific quantiles) we have argued are also relevant in other areas of literature such as in the estimation of neighborhood effects. We demonstrated how the estimation objective and methodology fit together in an empirical application using a matched dataset on North Carolina elementary school students and teachers. Our results suggest that the detection of low-performing teachers -- a typical exercise with high-stakes consequences -- is sensitive to the estimator choice regardless of how teacher quality is defined. In particular, using nonparametric EB estimators that allow for dependence among the teacher-specific parameters is crucial for delivering optimal estimates of teacher quality and detecting low-performing teachers.

\setcounter{page}{1}

Proofs for Sections (ref) and (ref)

Proof of Theorem (ref)

Proofs of Lemmas (ref), (ref), and (ref)

Proofs for Section (ref)

Notations. Let $\phi(\cdot)$ denote the standard normal density, and $\gamma(\cdot | k,\theta)$ denote the density of a gamma random variable with shape $k$ and scale $\theta$. We sometimes use $\{ A \}$ as an indicator function that takes on a value of 1 when the event $A$ occurs and 0 otherwise; for instance, $\{ x \geq y \}$ takes on a value of 1 when $x \geq y$ and 0 otherwise. Let $C>0$ denote a generic constant depending only on $\underline{\sigma}$ and $k$. Let $\bar{\mu}_{G}$ and $\bar{\sigma}_{G}$ denote the maximum supports of $\mu$ and $\sigma$ for a distribution $G$ respectively. Let $\bar{\mu}_{\mathcal{G}}$ and $\bar{\sigma}_{\mathcal{G}}$ denote the maximum supports for $\mu$ and $\sigma$ of the class of distributions $\mathcal{G}$ respectively. That is, $\bar{\mu}_{\mathcal{G}} := \sup_{G \in \mathcal{G}}\bar{\mu}_{G}$. We will further define a class of distributions in addition to $\mathcal{G}_n$ defined in Assumption (ref):

where $\tilde{\gamma}_1 := \gamma_1 + \gamma_2 + \tfrac{1}{2}$ and $\tilde{\gamma}_2 := \gamma_2 + \tfrac{1}{2}$. Finally, we will sometimes drop the subscript $n$ (which indicates the row of the triangular array) for notational compactness.

Proof of Theorem (ref)

In this section, we split the proof of Theorem (ref) into three parts, each corresponding to one of the three bounds. We provide the proof for the regret bounds of estimating $\{\sigma_i^2\}_{i\in[n]}$ in detail, and then specify how it is modified to complete the proofs for $\{\sigma_i\}_{i\in[n]}$ and $\{\mu_i\}_{i\in[n]}$ respectively.

Leave-one-out regret

The following reduces the regret of using leave-one-out NPMLEs to a form more amenable for analysis.

Support of NPMLE maximization

First recall the definition of $\tilde{\mathcal{G}}$ in (ref).

MSEs within region $\mathcal{R}$

This section provides adaptations of results from Jiang2009 that reduce the MSE from using a density estimate in place of the true density $f_{G_{0,n}}$ to the Hellinger distance between the density estimate and the true density. In this regard, Theorems (ref), (ref) and (ref) are adaptations of Proposition 3 of Jiang2020 to the setting of unknown heteroskedasticity for the estimation objectives of $\{\sigma_i^2\}_{i\in[n]}$, $\{\sigma_i\}_{i\in[n]}$ and $\{\mu_i\}_{i\in[n]}$ respectively. Lemmas (ref) is used in the proofs of all three theorems, whereas Lemmas (ref) and (ref) are specifically for handling the additional differentiation term required in the proof of Theorem (ref).

Recall that the indicator function $\mathcal{R}$ is defined as in (ref). We maintain the following assumption for this subsection:

The following lemma (and Lemma (ref)) is an extension of Lemma 1 of Jiang2009 from a univariate to a bivariate density of sufficient statistics, which is relevant for our framework of unknown heteroskedasticity.

MSEs over region $\mathcal{R}^c$

Recall that the indicator function $\mathcal{R}$ is defined as in (ref), with $\mathcal{R}^c$ being its complement.

Auxiliary Lemmas

NPMLE density estimation

This section provides adaptations of results in Ghosal2001 on maximum likelihood estimation of normal mixtures to accommodate our setting\footnote{That is, to accommodate bivariate distributions on $(\mu,\sigma)$ with possibly growing supports for both $\mu$ and $\sigma$ and also multiple observations per unit.}. In this regard, Lemmas (ref), (ref)+(ref) and (ref) are extensions of Lemma 3.4, Theorem 3.3 and Theorem 4.1 of Ghosal2001 respectively.

\noindentNotations. The notation $G(A \times B)$ for a bivariate distribution $G$ and sets $A,B\subset\mathbb{R}$ denotes the probability that $G$ places on the event $(\mu\in A,\sigma\in B)$. The notation $X\lesssim Y$ means that there exists a constant $C$, possibly dependent on $(\gamma_1,\gamma_2,\underline{\sigma},L)$ but independent of $\epsilon$, such that $X \leq C Y $. Furthermore, all (generic) constants $C$ defined are also possibly dependent on $(\gamma_1,\gamma_2,\underline{\sigma},L)$ but independent of $\epsilon$. Let $\vee_i x_i$ denote the maximum of $\{x_i\}_{i\in[k]}$. Let $|| \cdot ||_{\infty}$ denote the sup-norm and for two densities $f$ and $f'$, let $d(f,f')$ denotes their Hellinger distance. Let $N_{[\,]}(\epsilon,\mathcal{F},m)$ and $N(\epsilon,\mathcal{F},m)$ denote the $\epsilon$-bracketing and $\epsilon$-covering numbers respectively for the set $\mathcal{F}$ under metric $m$ (see, e.g., Ghosal2001 for precise definitions).

First introduce the following definition

where $(\gamma_1,\gamma_2,\underline{\sigma},L)>0$ are taken to be constants throughout our analysis.

Let $x\in\mathbb{R}^k$ and define

for $G\in \mathcal{G}(\epsilon)$ and $k$ is also a constant throughout our analysis. Finally let