Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
52,594 characters · 11 sections · 51 citation commands
INFERENCE ON A DISTRIBUTION FROM NOISY DRAWS
\def\spacingset#1 \spacingset{1}
\if11 \fi
\if01 {
} \fi
{\bf JEL Classification:} C14, C23
{\bf Keywords:} bias correction, estimation noise, nonparametric inference, measurement error, panel data, regression to the mean, shrinkage.
\spacingset{1.45}
\setcounter{equation}{0}
Let $\theta_1,\ldots,\theta_n$ be a random sample from a distribution $F$ that is of interest. Suppose that we only observe noisy measurements of these variables, say $\vartheta_1,\ldots,\vartheta_n$. A popular approach is to do inference on $F$ and its functionals using the empirical distribution of $\vartheta_1,\ldots,\vartheta_n$. This is common practice when analyzing panel data with heterogenous coefficients. In the literature on student achievement, for example, $\theta_i$ is a teacher effect, $\vartheta_i$ is an estimator of it obtained from data on student test scores, and we care about the distribution of teacher value-added (see, e.g., JacksonRockoffStaiger2014 for an overview). In the same vein, Guvenen2009, BrowningEjrnaesAlvarez2010, and MagnacRoux2019 estimate heterogenous earning profiles, while AhnChoiGaleKariv2014 find substantial heterogeneity in ambiguity aversion In a nonlinear fixed-effect model, the marginal effect is heterogenous across units and interest lies in the distribution of these effects as well as its functionals (Chamberlain1984, HahnNewey2004). Although the plug-in approach is popular, using $\vartheta_1,\ldots,\vartheta_n$ rather than $\theta_1,\ldots,\theta_n$ introduces bias that is almost entirely ignored in practice. BarrasGagliardiniScaillet2018, who are interested in the distribution of the skill of fund managers, find that not accounting for bias leads to substantial overestimation of tail mass and misses to pick up the substantial asymmetry in the skill distribution.
We analyze the properties of the plug-in estimator of $F$ in a location-scale setting where $$ \vartheta_i = \theta_i + \frac{\sigma_i}{\sqrt{m}} \, \varepsilon_i, \qquad \varepsilon_i \, \vert \, (\theta_i,\sigma_i^2) \sim \mathrm{i.i.d.}~(0,1), $$ where $m$ is a parameter that grows with $n$. As the variance of the (heteroskedastic) noise is $\sigma_i^2/m$, this device shrinks the noise as the sample size grows. This is a very natural asymptotic embedding in settings where $\vartheta_i$ is an estimator of $\theta_i$ obtained from a sample of size $m$, as in a panel data setting or meta-analysis Vivalt2015. It is related to, yet different from, an approach based on small measurement-error approximations as in Chesher1991,Chesher2017,\footnote{ Chesher1991 provides expansions for densities, while we focus on distribution and quantile functions. Chesher2017 discusses the impact of noise in the explanatory variables in a quantile-regression model; this is a different setup than the one considered here. EvdokimovZeleneev2020 use our device of measurement error that shrinks with the sample size to correct inference in generalized method-of-moment problems. } and has precedent in the analysis of fixed-effect models for panel data, although for different purposes, as discussed in more detail below (see, e.g., HahnKuersteiner2002 and AlvarezArellano2003).
Efron2011 essentially entertains the homoskedastic setting with normal noise, where $$ \vartheta_i \vert \, \theta_i \sim N(\theta_i,\sigma^2/m), $$ and defines selection bias as the tendency of the $\vartheta_i$'s associated with the (in magnitude) largest $\theta_i$'s to be larger than their corresponding $\theta_i$. He proposes to deal with selection bias by using the well-known Empirical Bayes estimator of Robbins1956, which here is equal to $$ \vartheta_i + \frac{\sigma^2}{m} \, \nabla^1 \log p(\vartheta_i), $$ where $p$ is the marginal density of the $\vartheta_i$ and $\nabla^1$ denotes the first-derivative operator. For example, when $\theta_i\sim N(0,\psi^2)$ this expression then yields the (infeasible) shrinkage estimator $$ \left(1-\frac{\sigma^2/m}{\sigma^2/m+\psi^2} \right) \, \vartheta_i, $$ a parametric plug-in estimator of which would be the JamesStein1961 estimator. More generally, non-parametric implementation would also require estimation of $p$ and its first derivative. Shrinkage to the overall mean (in this case zero) is intuitive, as selection bias essentially manifests itself through the tails of the empirical distribution of the $\vartheta_i$ being too thick.\footnote{The same shrinkage factor is applied to each $\vartheta_i$, a consequence of the noise being homoskedastic. How to deal with heteroskedastic noise in an Empirical Bayes framework is not obvious. Discussion and a recent contribution can be found in and WeinsteinMaBrownZhang2018.} Shrinkage is commonly-applied in empirical work (see, e.g., Rockoff2004; ChettyFriedmanRockoff2014). It should be stressed, though, that, while shrinkage improves on $\vartheta_1,\ldots,\vartheta_n$ in terms of estimation risk, it does not lead to preferable estimators of the distribution $F$ or its moments.
The approach taken here is different from Efron2011. Without making parametric assumptions on $F$, we calculate the (leading) bias of the naive plug-in estimator of the distribution, $$ \hat{F}(\theta) := n^{-1} \sum_{i=1}^n \mathrm{1}\lbrace \vartheta_i \leq \theta \rbrace. $$ This calculation allows to construct estimators that correct for the bias directly. In the James-Stein problem, where $\theta_i \sim N(\eta,\psi^2)$, for example, the bias under homoskedastic noise equals $$ - \frac{\theta-\eta}{2} \, \frac{\sigma^2/\psi^2}{m} \phi\left(\frac{\theta-\eta}{\psi}\right) + o(m^{-1}). $$ Thus, the empirical distribution is indeed upward biased in the left tail and downward biased in the right tail. A bias order of $m^{-1}$ implies incorrect coverage of confidence intervals unless $n/m^2\rightarrow 0$. We present plug-in and jackknife estimators of the leading bias and show that the bias-corrected estimators are asymptotically normal with zero mean and variance $F(\theta)\, (1-F(\theta))$ as long as $n/m^3\rightarrow0$. So, bias correction is preferable to the naive plug-in approach for typical data sizes encountered in practice, where $m$ tends to be quite small relative to $n$. We also provide corresponding bias-corrected estimators of the quantile function of $F$.
If the distribution of $\sigma_i\, \varepsilon_i$ is fully known, recovering $F$ is a (generalized) deconvolution problem that can be solved for fixed $m$. Deconvolution-based estimators are well studied (see, e.g., CarrollHall1988 and DelaigleMeister2008). However, they have a very slow rate of convergence and it is well documented that they can behave quite poorly in small samples.\footnote{There are also solutions to the measurement-error problem based on repeated measurements (or instrumental variables), coupled with suitable independent restrictions (see, for example, HorowitzMarkatou1996, LiVuong1998, Hu2008, HuSchennach2008, and BonhommeJochmansRobin2016b,BonhommeJochmansRobin2016a). These can be useful alternatives in static models for panel data, where the object of interest is the distribution of the random intercept, as in the work of HorowitzMarkatou1996, for example.} In response to this, Efron2016 has recently argued for a return to a more parametric approach. Our approach delivers intuitive estimators that enjoy the usual parametric convergence rate and are numerically well behaved. Although it does not deliver a fixed-$m$ consistent estimator, bias correction further ensures that size-correct inference can be performed, provided that $n/m^3$ is small. It is not clear how to conduct inference based on deconvolution estimators.
Working out the statistical properties of $\hat{F}$ (and of its quantile function) is non-trivial because $\hat{F}$ is a non-smooth function of the data $\vartheta_1,\ldots,\vartheta_n$. As such, the approach taken here is different from, and complementary to, recent work on estimating average marginal effects in panel data models, which only looks at smooth functionals such as the mean and variance (see, e.g., FernandezValLee2013; OkuiYanagi2016). The impact of noise on smooth transformations of the $\vartheta_i$ can be handled using conventional methods based on Taylor-series expansions. We contrast such an approach with our derivations below. How to perform inference on the quantiles of marginal effects in nonlinear panel models is a long-standing open question DhaeneJochmans2015, and the current work can be seen as a first step in that direction.
In work contemporaneous to our own, OkuiYanagi2018 derive the bias of a kernel-smoothed estimator of $F$ and its derivative. Such smoothing greatly facilitates the calculation of the bias, making it amenable to conventional analysis. However, it also introduces additional bias terms that require much stronger moment conditions as well as further restrictions on the relative growth rates of $n$, $m$, and the bandwidth that governs the smoothing. Nevertheless, the (leading) bias term obtained in OkuiYanagi2018 coincides with ours in Proposition (ref) below. Additional discussion on and comparison between the two different approaches is given in OkuiYanagi2018.
Let $F$ be a univariate distribution on the real line. We are interested in estimation of and inference on $F$ and its quantile function $ q(\tau) : = \inf_{\theta} \lbrace \theta: F(\theta) \geq \tau \rbrace. $ If a random sample $\theta_1,\ldots,\theta_n$ from $F$ would be available this would be a standard problem. We instead consider the situation where $\theta_1,\ldots,\theta_n$ themselves are unobserved and we observe noisy measurements $\vartheta_1,\ldots,\vartheta_n$, with variances $\sigma_1^2/m,\ldots,\sigma_n^2/m$ for a positive real number $m$ which, in our asymptotic analysis below, will be required to grow with $n$. We assume the following.
Our setup reflects a situation where the noisy measurements $\vartheta_1,\ldots,\vartheta_n$ converge in squared mean to $\theta_1,\ldots,\theta_n$ at the rate $m^{-1}$. A leading case is the situation where $\vartheta_i$ is an estimator of $\theta_i$ obtained from a sample of size $m$ that converges at the parametric rate.\footnote{ Everything to follow can be readily modified to different convergence rates as well as to the case where $$ \mathrm{var}(\vartheta_i \vert \, \theta_i,\sigma_i^2)= \sigma_i^2/m_i , $$ with $m_i := p_i m$ for a random variable $p_i\in (0,1]$. It suffices to redefine $\sigma_i^2$ as $\sigma_i^2 / p_i$. When the $\vartheta_i$ represent estimators this device allows for the sample size to vary with $i$. For example, in a panel data setting, it would cover unbalanced panels under a missing-at-random assumption. Further, the requirement that $\vartheta_i$ is unbiased can be relaxed to allow for standard non-linearity bias of order $m^{-1}$. We do not do this here as it is possible quite generally to reduce the bias down to $O(m^{-3/2})$, for example via a jackknife or bootstrap correction, making it negligible in our analysis below. Furthermore, the split-sample jackknife approach to bias correction that we discuss below would automatically take care of this additional $m^{-1}$ bias without modification. } We allow $\theta_i$ and $\sigma_i^2$ to be correlated, implying that the noise $\vartheta_i - \theta_i$ is not independent of $\theta_i$. Hence, we allow for measurement error to be non-classical. Recovering the distribution of $\theta_i$ from a sample of $(\vartheta_i,\sigma_i^2)$ is, therefore, not a standard deconvolution problem.
It is common to estimate $F(\theta)$ by $$ \hat{F}(\theta) : = n^{-1} \sum_{i=1}^n \mathrm{1}\lbrace \vartheta_i \leq \theta \rbrace, $$ the empirical distribution of the $\vartheta_i$ at $\theta$. As we will show below, under suitable regularity conditions, such plug-in estimators are consistent and asymptotically normal as $n\rightarrow\infty$ provided that $m$ grows with $n$ so that $n/m^2$ converges to a finite constant. The use of $\vartheta_1,\ldots,\vartheta_n$ rather than $\theta_1,\ldots,\theta_n$ introduces bias of the order $m^{-1}$, in general. This bias implies that test statistics are size distorted and the coverage of confidence sets is incorrect unless $n/m^2$ converges to zero.
The bias problem is easy to see (and fix) when interest lies in smooth functionals of $F$, $$ \mu:=E(\varphi(\theta_i)), $$ for a (multiple-times) differentiable function $\varphi$. An (infeasible) plug-in estimator based on $\theta_1,\ldots,\theta_n$ would be $$ \tilde{\mu} : = n^{-1} \sum_{i=1}^n \varphi(\theta_i). $$ Clearly, this estimator is unbiased and satisfies $ \tilde{\mu}\overset{a}{\sim} N(\mu,\sigma_\mu^2/n) $ as soon as $\sigma^2_\mu :=\mathrm{var}(\varphi(\theta_i))$ exists. For the feasible plug-in estimator of $\mu$, $$ \hat{\mu} := n^{-1} \sum_{i=1}^n \varphi(\vartheta_i), $$ under standard regularity conditions, a Taylor-series expansion of $\varphi_i(\vartheta_i)$ around $\theta_i$ yields $$ E(\hat{\mu}-\mu) = \frac{b_\mu}{m} + O(m^{-3/2}), \qquad b_\mu : = \frac{E(\nabla^2\varphi(\theta_i) \, \sigma_i^2)}{2}, $$ and $$ \mathrm{var}(\hat{\mu}) = \frac{\sigma^2_\mu}{n} + O\left(n^{-1}m^{-1} \right). $$ Hence, letting $z\sim N(0,1)$, we have $$ \frac{\hat{\mu}-\mu}{\sigma_\mu/\sqrt{n}} \overset{a}{\sim} z + \sqrt{\frac{n}{m^2}} \frac{b_\mu}{\sigma_\mu} \sim N(c\, b_\mu/\sigma_\mu,\sigma_\mu^2) , $$ as $n/m^2\rightarrow c^2<\infty$ when $n,m\rightarrow\infty$. The noise in $\vartheta_1,\ldots,\vartheta_n$ introduces bias unless $\varphi$ is linear. It can be corrected for by subtracting a plug-in estimator of $b_\mu/m$ from $\hat{\mu}$. Doing so, again under regularity conditions, delivers and estimator that is asymptotically unbiased as long as $n/m^3\rightarrow 0$.
The machinery from above cannot be applied to deduce the bias of $\hat{F}$ as it is a step function and, hence, non-differentiable. We will derive its leading bias under the following conditions. To state them, we let $$ \varepsilon_i := \frac{\vartheta_i-\theta_i}{\sigma_i/\sqrt{m}} $$ and write $f$ for the density function of $F$.
Assumption (ref) imposes smoothess on certain densities and conditional expectations but not on the estimator of $F$.
Define the function $$ \beta(\theta) : = \frac{E(\sigma_i^2 \vert \theta_i = \theta) \, f(\theta)}{2}, $$ which is well-behaved under Assumption (ref), and let $$ b_F(\theta) : = \beta^\prime(\theta) $$ be its derivative. We also introduce the covariance function $$ \sigma_F(\theta,\theta^\prime) : = F(\theta \wedge \theta^\prime ) - F(\theta) \, F(\theta^\prime), $$ where we use $\theta \wedge \theta^\prime $ to denote $\min\lbrace\theta,\theta^\prime\rbrace$. Proposition (ref) summarizes the large-sample properties of $\hat{F}$.
To illustrate the result suppose that $\sigma_i^2$ is independent of $\theta_i$ and that $\theta_i$ has density function $$ f(\theta) = \frac{1}{\psi} \phi\left( \frac{\theta-\eta}{\psi} \right), $$ as in the JamesStein1961 problem. Letting $\sigma^2$ denote the mean of the $\sigma_i^2$ an application of Proposition (ref) yields $$ b_F(\theta) = - \frac{\theta-\eta}{2} \, \frac{\sigma^2}{\psi^2} \, \phi\left( \frac{\theta-\eta}{\psi} \right). $$ Thus, $\hat{F}(\theta)$ is upward biased when $\theta<\eta$ and is downward biased when $\theta>\eta$. This finding is a manifestation of the phenomenon of regression to the mean (or selection bias, or the winner's curse; see Efron2011). It implies that the empirical distribution tends to be too disperse.
The bias in $\hat{F}$ translates to bias in estimators of the quantile function. A natural estimator of the quantile function is the left-inverse of $\hat{F}$. With this definition, the plug-in estimator of the $\tau$th-quantile is $$ \hat{q}(\tau) := \vartheta_{(\lceil \tau n \rceil)} , $$ where $\vartheta_{(\lceil \tau n \rceil)}$ is the $\lceil \tau n \rceil$th order statistic of our sample, where $\lceil a \rceil $ delivers the smallest integer at least as large as $a$.
To calculate the leading bias in $\hat{q}(\tau)$ observe that it is an (approximate) solution to the empirical moment condition $$ \hat{F}(q) - \tau = 0 $$ (with respect to $q$). From Proposition (ref) we know that $$ E(\hat{F}(q(\tau))) - \tau = \frac{b_F(q(\tau))}{m} + O(m^{-3/2}) , $$ uniformly in $\tau$, so the moment condition that defines the estimator $\hat{q}(\tau)$ is biased. Letting $$ b_q(\tau) : = -\frac{b_F(q(\tau))}{f(q(\tau))}, \qquad \sigma_q^2(\tau) : = \frac{\tau(1-\tau)}{f(q(\tau))^2}, $$ we obtain the following result.
As an example, when $\theta_i\sim N(\eta,\psi^2)$, independent of $\sigma_i^2$, we have $$ b_q(\tau) = \frac{\sigma^2/\psi^2}{2} \, (q(\tau)-\eta), $$ which, in line with our discussion on regression to the mean above, is positive for all quantiles below the median and negative for all quantiles above the median. The median itself is, in this particular case, estimated without plug-in bias of order $m^{-1}$. It will, of course, still be subject to the usual $n^{-1}$ bias arising from the nonlinear nature of the estimating equation.
Propositions (ref) and (ref) complement the existing results on the bias in smooth functionals (FernandezValLee2013; OkuiYanagi2016) of the distribution of heterogenous parameters in panel data models. Our calculations confirm that the order of the bias in the empirical distribution and in the quantile function is of the same order as in the smooth case, $m^{-1}$.
Importantly, our results validate a traditional jackknife approach to bias correction as in HahnNewey2004 and DhaeneJochmans2015. Such an approach exploits the fact that the bias is proportional to $m^{-1}$ and is based on re-estimating $\theta_1,\ldots,\theta_n$ from subsamples. The simplicity of such a method makes it very useful in panel data applications, for example.
To illustrate how the jackknife would work here, consider a stationary (balanced) $n\times m$ panel. Let $\vartheta_{i,m_1}$ be an estimator of $\theta_i$ constructed from the $n\times m_1$ subpanel consisting of the first $m_1$ cross sections only. Then $$ \hat{F}_{m_1}(\theta):= n^{-1} \sum_{i=1}^n \mathrm{1}\lbrace \vartheta_{i,m_1}\leq \theta \rbrace $$ is the plug-in estimator of $F(\theta)$ based on this subpanel alone. From Proposition (ref) it follows that $$ E(\hat{F}_{m_1}(\theta)) = F(\theta) + \frac{b_F(\theta)}{m_1} + O(m_1^{-3/2}). $$ Using the remaining $m_2 := m - m_1$ cross sections from the full panel we can equally calculate estimators $\vartheta_{i,m_2}$ and subsequently construct $$ \hat{F}_{m_2}(\theta):= n^{-1} \sum_{i=1}^n \mathrm{1}\lbrace \vartheta_{i,m_2}\leq \theta \rbrace, $$ for which $$ E(\hat{F}_{m_2}(\theta)) = F(\theta) + \frac{b_F(\theta)}{m_2} + O(m_2^{-3/2}) $$ follows in the same way. Consequently, $$ \tilde{b}_F(\theta) := {m_1} \hat{F}_{m_1}(\theta) + {m_2} \hat{F}_{m_2}(\theta) -m \hat{F}(\theta) $$ is a split-panel jackknife estimator of the leading bias term $b_F(\theta)$. Hence, $$ \tilde{F}(\theta):= \hat{F}(\theta) - \frac{\tilde{b}_F(\theta)}{m}. $$ is a nonparametric bias-corrected estimator.
A jackknife estimator of the quantile function can be defined in the same way. Moreover, let $\vartheta_{(\lceil \tau n \rceil), m_1}$ and $\vartheta_{(\lceil \tau n \rceil), m_2}$ be the $\lceil \tau n \rceil$ order statistic of the re-estimated quantities in the first and second subsample, respectively. Recall that $\vartheta_{(\lceil \tau n \rceil), m_1}$ is the (approximate) solution to $ \hat{F}_{m_1}(q) - \tau = 0, $ and so is our estimator of $q(\tau)$ as obtained from the information in the $n\times m_1$ subpanel only. As before, $$ \tilde{b}_q(\tau) := {m_1} \vartheta_{(\lceil \tau n \rceil), m_1} + {m_2} \vartheta_{(\lceil \tau n \rceil), m_2} -m \vartheta_{(\lceil \tau n \rceil)} $$ is a nonparametric estimator of $b_q(\tau)$ that gives rise to a jackknife bias-corrected estimator of the quantile function.
The large-sample behavior of these jackknife estimators is the same as for the analytic corrections in Propositions (ref) and (ref) below. The split-sample jackknife is simple to implement but requires access to the original data from which $\vartheta_1,\ldots,\vartheta_n$ were computed. This can be infeasible in meta-analysis problems, where each of the $\vartheta_i$ is an estimator constructed from a different data set that need not all be accessible. It can also be complicated in structural econometric models, where $\vartheta_i$ may be the solution to a cumbersome optimization programme that can be time-consuming to solve. We discuss an alternative bias-correction estimator next.
We will formulate regularity conditions for a plug-in estimator of the bias to be consistent under the maintained assumption that the $\sigma_1^2,\ldots, \sigma_m^2$ are known. We conjecture that, under suitable conditions, the results below will continue to go through when the $\sigma_i^2$ are replaced by estimators.
A bias-corrected estimator based on Proposition (ref) takes the form $$ \check{F}(\theta) : = \hat{F}(\theta) - \frac{\hat{b}_F(\theta)}{m}, \qquad \hat{b}_F(\theta) : = - \frac{(nh^2)^{-1} \sum_{i=1}^n \sigma_i^2 \, k^\prime\left(\frac{\vartheta_i-\theta}{h}\right)}{2}, $$ where $k^\prime$ is the derivative of kernel function $k$ and $h$ is a non-negative bandwidth parameter. Thus, we estimate the bias using standard kernel methods. For simplicity, we will use a Gaussian kernel throughout, so $k^\prime(\eta):=-\eta \, \phi(\eta)$.
We establish the asymptotic behavior of $\check{F}$ under the following conditions.
Assumption (ref) contains simple smoothness and boundedness requirements on the conditional density of $\theta_i$ given $\sigma_i^2$, as well as tail conditions on the marginal density of the $\theta_i$ and on the bias function $b_F(\theta)$.
We have the following result.
The implications of Proposition (ref) are qualitatively similar to those for smooth functionals discussed above. Indeed, for any fixed $\theta$, it implies that $$ \check{F}(\theta) \overset{a}{\sim} N(F(\theta), F(\theta)(1-F(\theta))/n) $$ as $n\rightarrow\infty$ and $m\rightarrow\infty$ with $n/m^4\rightarrow 0$. Thus, the leading bias is removed from $\hat{F}$ without incurring any cost in terms of (asymptotic) precision. Given the correction term, the sample variance of $$ \mathrm{1}\lbrace \vartheta_i\leq \theta \rbrace + \frac{1}{2}\frac{1}{mh^2} \sigma_i^2 \kappa^\prime\left(\frac{\vartheta_i-\theta}{h}\right) $$ is a more natural basis for inference in small samples than is that of $\mathrm{1}\lbrace \vartheta_i\leq \theta \rbrace$.
A data-driven way of choosing $h$ is by cross validation. A plug-in estimator of the integrated squared error $\int_{-\infty}^{+\infty} (\check{F}(\theta)-F(\theta))^2 \, d\theta$ (up to multiplicative and additive constants) is
where we use the shorthand $$ \underline{\phi}^\prime(\vartheta_i,\vartheta_j;h) : = \frac{1}{4} \frac{1}{\sqrt{2}h} \phi\left(\frac{\vartheta_i-\vartheta_j}{\sqrt{2}h} \right) \left(\frac{1}{2}-\frac{(\vartheta_i+\vartheta_j)^2}{4h^2} + \frac{\vartheta_i \vartheta_j}{h^2} \right). $$ See Appendix C for details on the derivation. The cross-validated bandwidth then is $\check{h}:=\arg\min_h v(h)$ on the interval $(0,+\infty)$.
Now turn the bias-corrected estimation of the quantile function. Proposition (ref) readily suggests a bias-corrected estimator of the form $$ \hat{q}(\tau) - \frac{\hat{b}_q(\tau)}{m}, \qquad \hat{b}_q(\tau) : = - \frac{\hat{b}_F(\hat{q}(\tau))}{\hat{f}(\hat{q}(\tau))}, $$ using obvious notation. While (under suitable regularity conditions) such an estimator successfully reduces bias it has the unattractive property that it requires a non-parametric estimator of the density $f$, which further shows up in the denominator.
An alternative estimator that avoids this issue is $$ \check{q}(\tau):= \vartheta_{(\lceil \hat{\tau}^* n \rceil)}, \qquad \hat{\tau}^* : = \tau + \frac{\hat{b}_F(\hat{q}(\tau))}{m} , $$ The justification for this estimator comes from the fact that $E(\hat{F}(q(\tau))) - \tau^* = O(m^{-2})$, where $\tau^* = \tau + {b_F(q(\tau))}/ {m}$, and its interpretation is intuitive. Given the noise in the $\vartheta_i$ relative to the $\theta_i$, the empirical distribution of the former is too heavy-tailed relative to the latter, and so $\hat{q}(\tau)$ estimates a quantile that is too extreme, on average. Changing the quantile of interest from $\tau$ to $\tau^*$ adjusts the naive estimator and corrects for regression to the mean.
The corrected estimator has the same asymptotic variance as the uncorrected estimator. It is well-known that plug-in estimators of $\sigma_q^2$ can perform quite poorly in small samples (MaritzJarrett1978). Typically, researchers rely on the bootstrap, and we suggest doing so here. Moreover, draw (many) random samples of size $n$ from the original sample $\vartheta_1,\ldots,\vartheta_n$ and re-estimate $q(\tau)$ by the bias-corrected estimator for each such sample. Then construct confidence intervals for $q(\tau)$ using the percentiles of the empirical distribution of these estimates. Note that, again, this bootstrap procedure does not involve re-estimation of the individual $\theta_i$.
To support our theory we provide simulation results for a JamesStein1961 problem where $\theta_i\sim N(0,\psi^2)$ and we have access to an $n\times m$ panel on independent realizations of the random variable $$ x_{it} \vert \, \theta_i\sim N(\theta_i,\sigma^2). $$ This setup is a simple random-coefficient model. It is similar to the classic many normal means problem of NeymanScott1948. While their focus was on consistent estimation of the within-group variance, $\sigma^2$, for fixed $m$, our focus is on between-group characteristics and the distribution of the $\theta_i$ as a whole. We estimate $\theta_i$ by the fixed-effect estimator, i.e., $$ \vartheta_i = m^{-1}\sum_{t=1}^m x_{it}. $$ The sampling variance of $\vartheta_i \vert \theta_i$ is $\sigma^2/m$. Rather than assuming this variance to be known we implement our analytical bias correction using the estimator $$ s_i^2 := (m-1)^{-1} \sum_{t=1}^m (x_{it}-\vartheta_i)^2 . $$ We do not make use of the fact that the $\vartheta_i$ are homoskedastic in estimating the noise or in constructing the bias correction. Moreover, the implementation of our procedure is non-parametric in the noise distribution.
A deconvolution argument implies that $$ \vartheta_i \sim N(0, \psi^2 + \sigma^2/m). $$ Thus, indeed, the empirical distribution of the fixed-effect estimator is too fat-tailed. In particular, the sample variance of $\vartheta_1,\ldots,\vartheta_n$, $$ \hat{\psi}^2 : = \frac{1}{n-1} \sum_{i=1}^n (\vartheta_i-\overline{\vartheta})^2, \qquad \overline{\vartheta} := n^{-1} \sum_{i=1}^n \vartheta_i, $$ is a biased estimator of $\psi^2$. To illustrate how this invalidates inference in typically-sized data sets we simulated data for $\psi^2=1$ (so $F$ is standard normal) and $\sigma^2=5$. The panel dimensions $(n,m)$ reported on are $(50,3)$, $(100,4)$, and $(200,5)$. Table (ref) shows the bias and standard deviation of $\hat{\psi}^2$ as well as the empirical rejection frequency of the usual two-sided $t$-test for the null that $\psi^2=1$. The nominal size is set to $5\%$. In practice, however, the test rejects in virtually all of the $10,000$ replications. The table provides the same summary statistics for the bias-corrected estimator $$ \check{\psi}^2 : = {\frac{1}{n-1} \sum_{i=1}^n\left( (\vartheta_i-\overline{\vartheta})^2-\frac{s_i^2}{m}\right)}. $$ The adjustment reduces the estimator's bias relative to its standard error and brings down the empirical rejection frequencies to just over their nominal value for the sample sizes considered.
A popular approach in empirical work to deal with noise in $\vartheta_1,\ldots,\vartheta_n$ is shrinkage estimation (see, e.g., ChettyFriedmanRockoff2014). This procedure is not designed to improve estimation and inference of $F$ or its moments, however. In the current setting, the (infeasible, parametric) shrinkage estimator is simply $$ \left(1 - \frac{\sigma^2/m}{\sigma^2/m+\psi^2} \right) \, \vartheta_i. $$ Its exact sampling variance is $$ \left(\frac{\psi^2}{\sigma^2/m+\psi^2} \right) \, \psi^2 = \psi^2 - \frac{\sigma^2/\psi^2}{m} + o(m^{-1}). $$ It follows that the sample variance of the shrunken $\vartheta_1,\ldots,\vartheta_n$ has a bias that is of the same order as that in the sample variance of $\vartheta_1,\ldots,\vartheta_n$. Interestingly, note that, here, this estimator overcorrects for the presence of noise, and so will be underestimating the true variance, $\psi^2$, on average.
The upper two plots in Figures (ref), (ref), and (ref) provide simulation results for the distribution function $F$ for the same Monte Carlo designs. The figures deal with the sample sizes $(50,3)$, $(100,4)$, $(200,5)$, respectively. The left plots contain (the average over the Monte Carlo replications of) the analytically bias-corrected estimator (solid blue line), with the bandwidth chosen according to a cross-validation procedure, together with $95\%$ confidence bands placed around in. Each of the plots also provide the average of the naive plug-in estimator (dashed red line), the empirical distribution of the Empirical-Bayes point estimates (dashed-dotted purple line), and the actual standard-normal distribution that is being estimated (solid black line).\footnote{Empirical Bayes was implemented non-parametrically (and correctly assuming homoskedasticity) based on the formula stated in the introduction using a kernel estimator and the optimal bandwidth that assumes knowledge of the normality of the target distribution.} The upper right plots in Figures (ref), (ref), and (ref) have the same structure, only now the bias-corrected estimator being plotted is the split-sample jackknife.
The simulations clearly show the substantial bias in the naive estimator. This bias becomes more pronounced relative to its standard error as the sample size grows and, indeed, $\hat{F}$ starts falling outside of the confidence bands (of the bias corrected estimator) as the sample size increases. The Empirical-Bayes estimator is less biased than $\hat{F}$. However, its bias is of the same order and so, as the sample size grows it does not move toward $F$ but, rather, towards $\hat{F}$.\footnote{Recall that the Empirical-Bayes estimator is not designed for inference on $F$ but, in stead, aims to minimize risk in estimating $\theta_1,\ldots,\theta_n$. In terms of RMSE it dominates $\vartheta_1,\ldots,\vartheta_n$. For the three sample sizes considered here, the RMSEs are $1.667$, $1.246$, and $1.000$ for the plug-in estimators and $1.233$, $1.018$, $.874$ for Empirical Bayes.} The confidence bands of $\check{F}$ and $\tilde{F}$ settle around $F$ as the sample grows. The results also show near identical performance of the split-sample approach and the analytical approach based on our bias formula. Indeed, the curves in the left and right plots are virtually indistinguishable.
The reduction in bias in our estimators of $F$ is again sufficient to bring the empirical size of tests in line with their nominal size. To see this Table (ref) provides empirical rejection frequencies of two-sided tests at the $5\%$ level for $F$ at each of its deciles using both $\hat{F}$ and $\check{F}$. The rejection frequencies based on the naive estimator are much too high for all sample sizes and deciles and get worse as the sample gets larger. Empirical size is much closer to nominal size after adjusting for noise, and this improvement is observed at all deciles of the distribution.
The lower two plots in Figures (ref), (ref), and (ref) provide corresponding simulation results for estimators of the deciles of $F$. The presentation is constructed around a QQ-plot of the standard normal, pictured as the black dashed-dotted line in each plot. Along the QQ-plot, the average (over the Monte Carlo replications) of the naive estimator (red), Empirical Bayes (purple), and the bias-corrected quantiles (blue) are shown by $*$ symbols. Again, the left plots deal with the analytical correction and the right plots show results for the split-sample approach. Confidence intervals around the corrected estimators (in blue,-o) are also again provided. Like the naive estimator, the Empirical Bayes estimators reported are the appropriate order statistics of $\vartheta_1,\ldots,\vartheta_n$, after shrinkage has been applied to each. Visual inspection reveals that the results are in line with those obtained for the distribution function. As the sample size grows, only $\check{q}$ successfully adjusts for bias arising from estimation noise in $\vartheta_1,\ldots,\vartheta_n$. Here, the split-sample correction is slightly more effective than our analytical approach.
We use quarterly panel data on a set of 48 consumer price index items in 52 US cities. The data span the period 1990--2007, yielding 72 time series observations. They were used by ParsleyWei1996, CruciniShintaniTsuruga2015, and OkuiYanagi2016,OkuiYanagi2018 to investigate the cross-sectional heterogeneity in deviations from the law of one price. Let $p_{cit}$ be the price of item $i$ in city $c$ at time $t$ and define the random variable $$ x_{cit} = \log \left( \frac{p_{cit}}{p_{1it}} \right) = \log (p_{cit}) - \log (p_{1it}) $$ for all $(52-1)\times 48=2448$ city/item combinations apart from the reference city (which here is Albuquerque, New Mexico). For each city/item combination we estimate the mean, standard deviation and first-order autocorrelation of $x_{cit}$ non-parametrically from the time dimension of our panel. Our interest lies in the distribution functions of their population counterparts. We estimate these three distributions by the empirical distributions of the cross-sectional estimates, and then correct for plug-in bias via the split-sample jackknife procedure. Our results complement the analysis of OkuiYanagi2018, which gives corresponding estimates of the associated density functions.
The results are collected in Figure (ref). The plots contain the empirical distribution functions (dashed red line) together with $95\%$ confidence bands based on the split-sample jackknife (shaded blue region). The correction for regression to the mean to the empirical distribution is clearly visible for the mean (left plot). It is also statistically significant, with the tails of the empirical distribution falling out of the confidence region. The sample standard deviation and autocorrelation obtained from the time series are biased estimators and so the empirical distribution function for these parameters (middle and right plot, respectively) suffer from an additional bias that is of the same order of magnitude as is the bias due to estimation noise (see the discussion on Footnote 4). The split-sample jackknife corrects for both these sources of bias automatically. Here, the bias adjustment leads to a pronounced shift of the empirical distribution; the corrected distribution functions all but stochastically dominate the naive plug-in estimators. The differences between the corrected and uncorrected functions are quantitatively large and, given the small standard error, they are also statistically significant.
In this paper, we have considered inference on the distribution of latent variables from noisy measurements. In an asymptotic embedding where the variance of the noise shrinks with the sample size, we have derived the leading bias in the empirical distribution function of the noisy measurements and suggested both an analytical and a jackknife correction. They provide a simple and numerically stable (approximate) solution to a generalized deconvolution problem that, in addition, yields valid inference procedures. The split-sample jackknife is particularly straightforward to implement and we recommend its use whenever possible.