Posterior Average Effects

TL;DR

This paper introduces posterior average effects (PAE), a method for estimating population averages of unobservables using sample-conditioned posterior means, justified by minimax error properties and applied to economic data.

Contribution

It develops the theoretical foundation for PAE, demonstrating their minimax optimality under misspecification and introducing a measure of informativeness for posterior conditioning.

Findings

PAE have minimum worst-case specification error under misspecification.

A measure of informativeness quantifies the robustness of PAE.

Empirical applications include neighborhood effects and income dynamics.

Abstract

Economists are often interested in estimating averages with respect to distributions of unobservables, such as moments of individual fixed-effects, or average partial effects in discrete choice models. For such quantities, we propose and study posterior average effects (PAE), where the average is computed conditional on the sample, in the spirit of empirical Bayes and shrinkage methods. While the usefulness of shrinkage for prediction is well-understood, a justification of posterior conditioning to estimate population averages is currently lacking. We show that PAE have minimum worst-case specification error under various forms of misspecification of the parametric distribution of unobservables. In addition, we introduce a measure of informativeness of the posterior conditioning, which quantifies the worst-case specification error of PAE relative to parametric model-based estimators. As illustrations, we report PAE estimates of distributions of neighborhood effects in the US, and of permanent and transitory components in a model of income dynamics.

Equations

$$Y_{ij}={\alpha }_i+{\epsilon }_{ij},i=1,...,n,j=1,...,J.$$

$$F_{\alpha }(a)=\mathbb{E}\left[1\{\alpha \le a\}\right],$$

$${\hat{F}}_{\alpha }^{\mathrm{F}E}(a)=\frac{1}{n}\sum _{i=1}^{n}1\{{\overline{Y}}_i\le a\},$$

$$\mathbb{E}[\alpha |Y=Y_i]={\hat{\mu }}_{\alpha }+\hat{\rho }({\overline{Y}}_i-{\hat{\mu }}_{\alpha }),i=1,...,n,$$

$${\hat{F}}_{\alpha }^{\mathrm{P}M}(a)=\frac{1}{n}\sum _{i=1}^{n}1\left\{{\hat{\mu }}_{\alpha }+\hat{\rho }({\overline{Y}}_i-{\hat{\mu }}_{\alpha })\le a\right\},$$

$${\hat{F}}_{\alpha }^{\mathrm{P}}(a)=\frac{1}{n}\sum _{i=1}^n\mathbb{E}\left[1\{\alpha \le a\}|Y=Y_i\right],$$

$${\hat{F}}_{\alpha }^{\mathrm{P}}(a)=\frac{1}{n}\sum _{i=1}^n\Phi \left(\frac{a-{\hat{\mu }}_{\alpha }-\hat{\rho }({\overline{Y}}_i-{\hat{\mu }}_{\alpha })}{{\hat{s}}_{\alpha }\sqrt{1-\hat{\rho }}}\right),$$

$${\hat{F}}_{\alpha }^{\mathrm{M}}(a)=\Phi \left(\frac{a-{\hat{\mu }}_{\alpha }}{{\hat{s}}_{\alpha }}\right),$$

$$\overline{\delta }={\left(\mathbb{E}[W{W}^{\prime }]\right)}^{-1}\mathbb{E}[W\alpha ].$$

$${\hat{\delta }}^{\mathrm{F}E}={\left(\sum _{i=1}^n{W}_i{W}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{W}_i{\overline{Y}}_i.$$

$${\hat{\delta }}^{\mathrm{P}}={\left(\sum _{i=1}^n{W}_i{W}_i^{\prime }\right)}^{-1}\sum _{i=1}^n{W}_i\left({\hat{\mu }}_{\alpha }+\hat{\rho }({\overline{Y}}_i-{\hat{\mu }}_{\alpha })\right),$$

$$Y_i=g_{\beta }(U_i,X_i),$$

$$\overline{\delta }={\mathbb{E}}_{f_0}\left[{\delta }_{\beta }(U,X)\right],$$

$$U={\left(\frac{\alpha -{\mu }_{\alpha }}{s_{\alpha }},\frac{{\epsilon }_1}{s_{\epsilon }},...,\frac{{\epsilon }_J}{s_{\epsilon }}\right)}^{\prime }.$$

$$\widehat{\delta}^{\mathrm{M}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{f_{\widehat{\sigma}}}\left[\delta_{\widehat{\beta}}(U,X)\,\big{|}\,X=X_{i}\right],$$

$$p_{\beta ,\sigma }(u|y,x)=\frac{f_{\sigma }(u|x)1\{u\in \mathcal{U}(y,x,\beta )\}}{\int f_{\sigma }(v|x)1\{v\in \mathcal{U}(y,x,\beta )\}dv}.$$

$$\widehat{\delta}^{\mathrm{P}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{p_{\widehat{\beta},\widehat{\sigma}}}\left[\delta_{\widehat{\beta}}(U,X)\,\Big{|}\,Y=Y_{i},X=X_{i}\right],$$

$${\Gamma }_{ϵ}=\left\{f_0:d(f_0,f_{{\sigma }_{\ast }})\le ϵ,{\mathbb{E}}_{P(\beta ,f_0)}[{\psi }_{\beta ,{\sigma }_{\ast }}(Y,X)]=0\right\}.$$

$$d(f_0,f_{\sigma })={\int }_{\mathcal{X}}{\int }_{\mathcal{U}}\varphi \left(\frac{f_0(u|x)}{f_{\sigma }(u|x)}\right)f_{\sigma }(u|x)f_X(x)dudx,$$

$${\hat{\delta }}_{\gamma }=\frac{1}{n}\sum _{i=1}^n{\gamma }_{\hat{\beta },\hat{\sigma }}(Y_i,X_i)+o_{P(\beta ,f_0)}(1).$$

$$b_{ϵ}(\gamma )={\sup }_{f_0\in {\Gamma }_{ϵ}}\left|{\mathbb{E}}_{P(\beta ,f_0)}[{\gamma }_{\beta ,{\sigma }_{\ast }}(Y,X)]-{\mathbb{E}}_{f_0}[{\delta }_{\beta }(U,X)]\right|.$$

$$b_{ϵ}(\gamma )=\left|{\mathbb{E}}_{\ast }[\gamma (Y,X)-\delta (U,X)]\right|$$

$$\displaystyle\,\,\,+\epsilon^{\frac{1}{2}}\bigg{\{}\frac{2}{\phi^{\prime\prime}(1)}{\mathrm{V}ar}_{*}\Big{(}\gamma(Y,X)-\delta(U,X)-\mathbb{E}_{*}\left[\gamma(Y,X)-\delta(U,X)\,|\,X\right]-\lambda^{\prime}\widetilde{\psi}(Y,X)\Big{)}\bigg{\}}^{\frac{1}{2}}+{\cal O}(\epsilon),$$

$$b_{ϵ}({\gamma }^{\mathrm{P}})={ϵ}^{\frac{1}{2}}{\left\{\frac{2}{{\varphi }^{\prime \prime }(1)}{\mathrm{Var}}_{\ast }\left(\delta (U,X)-{\mathbb{E}}_{\ast }[\delta (U,X)|Y,X]\right)\right\}}^{\frac{1}{2}}+\mathcal{O}(ϵ).$$

$$b_{ϵ}({\gamma }^{\mathrm{P}})={ϵ}^{\frac{1}{2}}{\left\{\frac{4}{{\varphi }^{\prime \prime }(1)}T\left(\frac{a-{\mu }_{\alpha }}{{\sigma }_{\alpha }},\sqrt{\frac{1-\rho }{1+\rho }}\right)\right\}}^{\frac{1}{2}}+\mathcal{O}(ϵ),$$

$${\gamma }^{\mathrm{P}}(y,x)={\mathbb{E}}_{\ast }[\delta (U,X)|Y=y,X=x].$$

$$b_{ϵ}({\gamma }^{\mathrm{P}})\le b_{ϵ}(\gamma )+\mathcal{O}(ϵ).$$

$$\overline{ϵ}$$

$$b_{ϵ}({\gamma }^{\mathrm{P}})\le b_{ϵ}(\gamma ).$$

$$b_{ϵ}({\gamma }^{\mathrm{P}})\le 2{\inf }_{\gamma }{b}_{ϵ}(\gamma ).$$

$$\underset{ϵ\to 0}{\lim }\frac{b_{ϵ}({\gamma }^{\mathrm{P}})}{b_{ϵ}({\gamma }^{\mathrm{M}})}=\frac{{\left\{{\mathrm{Var}}_{\ast }\left(v(U,X)-{\mathbb{E}}_{\ast }[v(U,X)|Y,X]\right)\right\}}^{\frac{1}{2}}}{{\left\{{\mathrm{Var}}_{\ast }\left(v(U,X)\right)\right\}}^{\frac{1}{2}}},$$

Full text

JEL codes:**  C13, C23.**

Keywords:**  model misspecification, robustness, sensitivity analysis, empirical Bayes, posterior conditioning, latent variables.**

1 Introduction

In many settings, applied researchers wish to estimate population averages with respect to a distribution of unobservables. This includes moments of individual fixed-effects in panel data, and average partial effects in discrete choice models, which are expectations with respect to some distribution of shocks or heterogeneity. The standard approach in applied work is to assume a parametric form for the distribution of unobservables, and to compute the average effect under that assumption. For example, in binary choice, researchers often assume normality of the error term, and compute average partial effects under normality. This “model-based” estimation of average effects is justified under the assumption that the parametric model is correctly specified.

In this paper, we consider a different approach, where the average effect is computed conditional on the observation sample. We refer to such estimators as “posterior average effects” (PAE). Posterior averaging is appealing for prediction purposes, and it plays a central role in Bayesian and empirical Bayes approaches (e.g., Berger, 1980, Morris, 1983). Here we focus instead on the estimation of population expectations. Our goal is twofold: to propose a novel class of estimators, and to provide a frequentist framework to understand when and why posterior conditioning may be useful in estimation. Our main result will show that PAE have robustness properties when the parametric model is misspecified.

PAE are closely related to empirical Bayes (EB) estimators, which are increasingly popular in applied economics. Consider a fixed-effects model of teacher quality, which is our main example. When the number of observations per teacher is small, the dispersion of teacher fixed-effects is likely to overstate that of true teacher quality, since teacher effects are estimated with noise. An alternative approach is to postulate a prior distribution for teacher quality — typically, a normal — and report posterior estimates, holding fixed the values of the mean and variance parameters. The hope is that such EB estimates, which are shrunk toward the prior, are less affected by noise than the teacher fixed-effects (e.g., Kane and Staiger, 2008, Chetty et al., 2014, Angrist et al., 2017). However, while EB estimates are well-justified predictors of the quality of individual teachers, it is not obvious how to aggregate them across teachers when the goal is to estimate a population average such as a moment or a distribution function.

As an example, suppose we wish to estimate the distribution function of teacher quality evaluated at a point. Since this quantity is an average of indicator functions, the PAE is simply an average of posterior means — that is, of EB estimates — of the indicator functions. This estimator is available in closed form. However, the PAE differs from the empirical distribution of the EB estimates of teacher effects. In particular, while the variance of EB estimates is too small relative to that of latent teacher quality, the PAE has the correct variance. Related applications of PAE include settings involving neighborhood/place effects (Chetty and Hendren, 2017, Finkelstein et al., 2017) or hospital quality (Hull, 2018).

Importantly, although posterior averages have desirable properties for predicting individual parameters, their usefulness for estimating population average quantities is not evident. For example, suppose that teacher quality is normally distributed. In this case, a model-based normal estimator of the distribution of teacher quality is consistent. Moreover, it is asymptotically efficient when means and variances are estimated by maximum likelihood. Hence, in the correctly specified case, there is no reason to deviate from the standard model-based approach and compute posterior estimators. The main insight of this paper is that, under misspecification — e.g., when teacher quality is not normally distributed — conditioning on the data using PAE can be beneficial.

To study estimators under misspecification, we focus on specification error, which is the population discrepancy between the probability limit of an estimator and the true parameter value. In our main results, we show that PAE have minimum worst-case specification error, where the worst case is computed in a nonparametric neighborhood of the reference parametric distribution (e.g., a normal). Specifically, we show that, when neighborhoods are defined in terms of the Pearson chi-squared divergence, PAE have minimum worst-case specification error within a large class of estimators, for any neighborhood size smaller than a threshold value that we characterize. In addition, when broadening the class of neighborhoods to $\phi$-divergences, we show that, while PAE do not have minimum worst-case specification error in general in fixed-size neighborhoods, they achieve minimum worst-case specification error under local misspecification, i.e. when the size of the neighborhood tends to zero.

In our examples and illustrations, we find that the information contained in the posterior conditioning is setting-specific. This is intuitive, since although PAE have minimum worst-case specification error under our conditions, the specification error is not zero in general and it varies between applications. PAE tend to behave better when the realizations of outcome variables (such as test scores) are more informative about the values of the unobservables (such as the quality of a teacher). Consistently with this intuition, our local result suggests quantifying the “informativeness” of the posterior conditioning using an easily computable $R^{2}$ coefficient.

While our theoretical results focus on population specification error, in practice PAE are also affected by sampling error, due to the fact that the sample size — e.g., the number of teachers — is not infinite. A common approach to account for both sampling variability and specification error is to focus on mean squared error. In general, PAE do not have minimum mean squared error: indeed, in finite samples, model-based estimators can have smaller mean squared error than PAE. In Bonhomme and Weidner (2018), we show how to construct estimators that minimize mean squared error under local asymptotic misspecification. However, such estimators depend on the neighborhood size. In contrast, PAE do not require taking a stand on the degree of misspecification through the size of the neighborhood, and they are simple to implement and do not depend on tuning parameters. To complement the theory, we report the results of a Monte Carlo simulation, where we compare the performance of the PAE to those of a model-based estimator and a nonparametric deconvolution-based estimator. We find that, while the model-based estimator tends to perform best under correct specification, the performance of the PAE appears less sensitive to misspecification than those of the model-based and nonparametric estimators.

To illustrate the scope of PAE for applications, we then consider two empirical settings. In the first one, we study the estimation of neighborhood/place effects in the US. Chetty and Hendren (2017) report estimates of the variance of neighborhood effects, as well as EB estimates of those effects. Our goal is to estimate the distribution of effects across neighborhoods. We find that, when using a normal prior as in Chetty and Hendren (2017), our posterior estimator of the distribution function of neighborhood effects across commuting zones is not normal. However, we also show through simulations and computation of our posterior informativeness measure that the signal-to-noise ratio in the data is not high enough to be confident about the exact shape of the distribution. Hence, in this setting, PAE inform our knowledge of the distribution of neighborhood effects, and motivate future analyses using more flexible model specifications and individual-level data.

In the second empirical illustration, our goal is to estimate the distributions of latent components in a permanent-transitory model of income dynamics (e.g., Hall and Mishkin, 1982, Blundell et al., 2008), where log-income is the sum of a random-walk component and a component that is independent over time. Researchers often estimate the covariance structure of the latent components in a first step. Then, in order to document distributions or to use the income process in a consumption-saving model, they often assume Gaussianity. However, there is increasing evidence that income components are not Gaussian (e.g., Geweke and Keane, 2000, Hirano, 2002, Bonhomme and Robin, 2010, Guvenen et al., 2016). We estimate posterior distribution functions of permanent and transitory income components using recent waves from the Panel Study of Income Dynamics (PSID). Our PAE estimates suggest some departure from Gaussianity, especially for the transitory income component.

We analyze several extensions. First, we describe the form of PAE in several models, including binary choice and censored regression. Second, we discuss how to construct confidence intervals and specification tests based on PAE. Lastly, we revisit the question of optimality of EB estimates for predicting individual parameters. By extending our misspecification analysis from worst-case specification error of sample averages to worst-case mean squared prediction error, we show that EB estimators remain optimal, up to small-order terms, under local deviations from normality.

Related literature and outline.

PAE are closely related to parametric EB estimators (Efron and Morris, 1973, Morris, 1983). For recent econometric applications of shrinkage methods (James and Stein, 1961, Efron, 2012), see Hansen (2016), Fessler and Kasy (2018), and Abadie and Kasy (2018). Recent contributions to nonparametric EB methods are Koenker and Mizera (2014) and Ignatiadis and Wager (2019).

Our analysis is also related to deconvolution and other nonparametric approaches. However, in our framework we allow for forms of misspecification under which the quantity of interest is not consistently estimable, and we search for estimators that have the smallest specification error.

In panel data settings, Arellano and Bonhomme (2009) study the asymptotic properties of random-effects estimators of averages of functions of covariates and individual effects. They show that, when the distribution of individual effects is misspecified whereas the other features of the model are correctly specified, PAE are consistent as $n$ and $T$ tend to infinity. By contrast, in our setup, only $n$ tends to infinity, and misspecification may affect the entire joint distribution of unobservables.

Our analysis also connects to the literature on robustness to model misspecification (e.g., Huber and Ronchetti, 2009, Kitamura et al., 2013, Andrews et al., 2017, 2020, Armstrong and Kolesár, 2018, Bonhomme and Weidner, 2018, Christensen and Connault, 2019). Here our aim is to propose and justify a class of simple, practical estimators.

The plan of the paper is as follows. In Section 2 we motivate the analysis by considering a fixed-effects model of teacher quality. In Section 3 we present our framework and derive our main theoretical results. In Section 4 we illustrate the use of PAE in two empirical settings. In Section 5 we describe several extensions. Finally, we conclude in Section 6. Replication codes are available as online material.

2 Motivating example: a fixed-effects model

To motivate the analysis, we start by considering the following model

[TABLE]

To fix ideas, we will think of $Y_{ij}$ as an average test score of teacher $i$ in classroom $j$, $\alpha_{i}$ as the quality of teacher $i$, and $\varepsilon_{ij}$ as a classroom-specific shock. There are $n$ teachers and $J$ observations per teacher. For simplicity, we abstract away from covariates (such as students’ past test scores), but those will be present in the framework we will introduce in the next section. Although here we focus on teacher effects, this model is of interest in other settings, such as the study of neighborhood effects, school effectiveness, or hospital quality, for example.

Suppose we wish to estimate a feature of the distribution of teacher quality $\alpha$. As an example, here we consider the distribution function of $\alpha$ at a particular point $a$,

[TABLE]

which is the percentage of teachers whose quality is below $a$. A first estimator is the empirical distribution of the fixed-effects estimates $\widehat{\alpha}_{i}=\overline{Y}_{i}=\frac{1}{J}\sum_{j=1}^{J}Y_{ij}$, for all teachers $i=1,...,n$; that is,

[TABLE]

where FE stands for “fixed-effects”. An obvious issue with this estimator is that $\overline{Y}_{i}=\alpha_{i}+\overline{\varepsilon}_{i}$ is a noisy estimate of $\alpha_{i}$, where $\overline{\varepsilon}_{i}=\frac{1}{J}\sum_{j=1}^{J}\varepsilon_{ij}$. Indeed, due to the presence of noise, for fixed $J$ and $n$ tends to infinity the distribution $\widehat{F}^{\mathrm{F}E}_{\alpha}$ tends to be too dispersed relative to $F_{\alpha}$ (although one can show that $\widehat{F}^{\mathrm{F}E}_{\alpha}(a)$ is consistent for $F_{\alpha}(a)$ as $J$ tends to infinity jointly with $n$ under mild conditions, see Jochmans and Weidner, 2018).

A different strategy is to model the joint distribution of $\alpha,\varepsilon_{1},...,\varepsilon_{J}$. A simple specification is a multivariate normal distribution with means $\mu_{\alpha}$ and $\mu_{\varepsilon}=0$, and variances $s_{\alpha}^{2}$ and $s_{\varepsilon}^{2}$. This specification can easily be made more flexible by allowing for different $s_{\varepsilon_{j}}^{2}$’s across $j$, for correlation between the different $\varepsilon_{j}$’s, or for means and variances being functions of covariates, for example. Under the assumption that all components are uncorrelated, $\mu_{\alpha}$, $s_{\alpha}^{2}$ and $s^{2}_{\varepsilon}$ can be consistently estimated for fixed $J$ as $n$ tends to infinity, using quasi-maximum likelihood or minimum distance based on mean and covariance restrictions.

Given estimates $\widehat{\mu}_{\alpha}$, $\widehat{s}_{\alpha}^{2}$, $\widehat{s}_{\varepsilon}^{2}$, we can compute empirical Bayes (EB) estimates (Morris, 1983) of the $\alpha_{i}$ as

[TABLE]

where the expectation is taken with respect to the posterior distribution of $\alpha$ given $Y=Y_{i}$ for $\widehat{\mu}_{\alpha}$, $\widehat{s}_{\alpha}^{2}$, $\widehat{s}_{\varepsilon}^{2}$ fixed, and $\widehat{\rho}=\frac{\widehat{s}_{\alpha}^{2}}{\widehat{s}_{\alpha}^{2}+\widehat{s}_{\varepsilon}^{2}/J}$ is a shrinkage factor. Here, $Y_{i}$ are vectors containing all $Y_{ij}$, $j=1,...,J$. The EB estimates in (3) are well-justified as predictors of the $\alpha_{i}$, since (when treating $\widehat{\mu}_{\alpha}$, $\widehat{s}_{\alpha}^{2}$, $\widehat{s}_{\varepsilon}^{2}$ as fixed) $\widehat{\mu}_{\alpha}+\widehat{\rho}(\overline{Y}_{i}-\widehat{\mu}_{\alpha})$ is the minimum mean squared error predictor of $\alpha_{i}$ under normality.

Given their rationale for prediction purposes, it is appealing to try and aggregate the EB estimates in order to estimate our target quantity $F_{\alpha}(a)$. A possible estimator is

[TABLE]

where PM stands for “posterior means”. For fixed $J$ as $n$ tends to infinity, the EB estimates tend to be less dispersed than the true $\alpha_{i}$, and $\widehat{F}^{\mathrm{P}M}_{\alpha}(a)$ is inconsistent in general. Indeed, while in large samples the variance of the fixed-effects estimates is $\rho^{-1}s_{\alpha}^{2}>s_{\alpha}^{2}$, the variance of the EB estimates is $\rho s_{\alpha}^{2}<s_{\alpha}^{2}$, where $\rho=\frac{s_{\alpha}^{2}}{s_{\alpha}^{2}+s_{\varepsilon}^{2}/J}$.

Instead of computing the distribution of EB estimates as in (4), a related idea is to compute the posterior distribution estimator

[TABLE]

where P stands for “posterior”. Using the normality assumption, we obtain

[TABLE]

where $\Phi$ denotes the distribution function of the standard normal. $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ is an example of a posterior average effect (PAE). One can check that it is consistent for fixed $J$ as $n$ tends to infinity, when the distribution of $\alpha,\varepsilon_{1},...,\varepsilon_{J}$ is normal. Under non-normality, $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ is generally inconsistent for fixed $J$ as $n$ tends to infinity. Moreover, the mean and variance of $\widehat{F}^{\mathrm{P}}_{\alpha}$ are $(1-\widehat{\rho})\widehat{\mu}_{\alpha}+\widehat{\rho}\frac{1}{n}\sum_{i=1}^{n}\overline{Y}_{i}$ and $(1-\widehat{\rho})\widehat{s}_{\alpha}^{2}+\widehat{\rho}^{2}\left[\frac{1}{n}\sum_{i=1}^{n}\overline{Y}_{i}^{2}-(\frac{1}{n}\sum_{i=1}^{n}\overline{Y}_{i})^{2}\right]$, respectively, which are consistent for $\mu_{\alpha}$ and $s_{\alpha}^{2}$ for fixed $J$ as $n$ tends to infinity.

The last estimator we consider here is directly based on the normal specification for $\alpha$,

[TABLE]

where M stands for “model”. This estimator enjoys attractive properties when the distribution of $\alpha,\varepsilon_{1},...,\varepsilon_{J}$ is indeed normal. In this case, $\widehat{F}^{\mathrm{M}}_{\alpha}(a)$ is consistent for fixed $J$ as $n$ tends to infinity, and it is efficient when $\widehat{\mu}_{\alpha}$ and $\widehat{s}_{\alpha}^{2}$ are maximum likelihood estimates. Moreover, the mean and variance of $\widehat{F}^{\mathrm{M}}_{\alpha}$ are $\widehat{\mu}_{\alpha}$ and $\widehat{s}_{\alpha}^{2}$, which are consistent irrespective of normality. However, when $\alpha,\varepsilon_{1},...,\varepsilon_{J}$ is not normally distributed, $\widehat{F}^{\mathrm{M}}_{\alpha}(a)$ is generally inconsistent for fixed $J$ as $n$ tends to infinity. Moreover, $\widehat{F}^{\mathrm{M}}_{\alpha}(a)$ only depends on the data through the mean $\widehat{\mu}_{\alpha}$ and the variance $\widehat{s}^{2}_{\alpha}$. In particular, $\widehat{F}^{\mathrm{M}}_{\alpha}$ is always normal, even when the data show clear evidence of non-normality.

Which one of these estimators should one use? The answer is not obvious, since they are all inconsistent as $n$ tends to infinity for fixed $J$ in general. In a framework that allows for misspecification of the normal distribution of $\alpha,\varepsilon_{1},...,\varepsilon_{J}$, we will show that the PAE $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ has minimum worst-case specification error in certain neighborhoods around the normal reference distribution. To our knowledge, unlike the other three estimators above, posterior estimators of distributions are novel to practitioners. They are easy to implement, and do not depend on additional tuning parameters. Our characterization provides a rationale for reporting them in applications, alongside other parametric and semi-parametric estimators.

Note that one may wish to relax normality by making the specification of $\alpha$, and possibly $\varepsilon_{j}$, more flexible. Deconvolution and nonparametric maximum likelihood estimators are often used for this purpose (e.g., Delaigle et al., 2008, Bonhomme and Robin, 2010, Koenker and Mizera, 2014). While these estimators may be consistent even when $\alpha$ is not normal, consistency relies on additional restrictions on the model. For example, the assumptions in Kotlarski (1967) require that $\alpha$, $\varepsilon_{1}$, …, $\varepsilon_{J}$ be mutually independent. By contrast, we do not impose any such additional conditions in our framework. In Section 3, we will show that asymptotically linear estimators have larger specification error than PAE under the form of misspecification that we consider.

To illustrate that an independence assumption among $\alpha$, $\varepsilon_{1}$, …, $\varepsilon_{J}$ can be restrictive, consider a situation where the researcher is concerned that the variance of $\varepsilon_{j}$ depends on $\alpha$. For instance, the variance of classroom-level shocks may depend on teacher quality. The presence of such conditional heteroskedasticity would invalidate conventional nonparametric deconvolution estimators. By contrast, we will show that $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ has minimum specification error in neighborhoods of distributions that allow for conditional heteroskedasticity. In Section 4 and the appendix, we will compare the finite-sample behavior of the parametric model-based estimator, the PAE, and a nonparametric deconvolution estimator, in data simulated from various specifications of model (1).

In model (1), the researcher may be interested in estimating other quantities. As an example, consider the coefficient in the population regression of teacher quality $\alpha$ on a vector of covariates $W$; that is,

[TABLE]

In applications, it is common to regress fixed-effects estimates on covariates to help interpret them (as in Dobbie and Fryer, 2013, among many others), and to compute

[TABLE]

Alternatively, one may regress the EB estimates of $\alpha_{i}$, as given by (3), on covariates (as in Angrist et al., 2017, and Hull, 2018, for example), and compute

[TABLE]

which is a PAE based on a normal reference specification for $\alpha$. We will see that, in our framework, the rationale for reporting $\widehat{\delta}^{\mathrm{P}}$ or $\widehat{\delta}^{\mathrm{F}E}$ depends on the form of misspecification that the researcher is concerned about.

The framework we describe next applies to the estimation of different quantities in a variety of settings. In Section 4 we apply PAE to model (1) and estimate the distribution of neighborhood/place effects in the US (Chetty and Hendren, 2017). In addition, we show that the permanent-transitory model of income dynamics (e.g., Hall and Mishkin, 1982) has a structure similar to model (1), and we report PAE estimates in this context. Lastly, in other models — such as static or dynamic discrete choice models and models with censored outcomes — our results motivate the use of PAE as complements to other estimators that researchers commonly report, and we provide examples in Section 5 and analyze them in the appendix.

3 Framework and main results

In this section we describe our framework to study PAE, and present our main results.

3.1 Model-based estimators and PAE

We consider the following class of models,

[TABLE]

where outcomes $Y_{i}$ and covariates $X_{i}$ are observed by the researcher, and $U_{i}$ are unobserved. The function $g_{\beta}$ is known up to the finite-dimensional parameter $\beta$. Our aim is to estimate an average effect of the form

[TABLE]

where $\delta_{\beta}$ is scalar, and known given $\beta$. Here $f_{0}$ denotes the true density of $U\,|\,X$. The expectation is taken with respect to the product $f_{0}f_{X}$, where $f_{X}$ is the marginal density of $X$. For conciseness we leave the dependence on $f_{X}$ implicit. While we focus on a scalar $\delta_{\beta}$, our results continue to hold in the vector-valued case, as we show at the end of this section. In Appendix S5, we discuss how to estimate quantities that depend on $f_{0}$ nonlinearly.

While the researcher does not know the true $f_{0}$, she has a reference parametric density $f_{\sigma}$ for $U\,|\,X$, which depends on a finite-dimensional parameter $\sigma$. We will allow $f_{\sigma}$ to be misspecified, in the sense that $f_{0}$ may not belong to $\{f_{\sigma}\}$. However, we will always assume that $g_{\beta}$ is correctly specified. In other words, misspecification will only affect the distribution of $U$ and its dependence on $X$, not the structural link between $(U,X)$ and outcomes.

To estimate $\overline{\delta}$ in (11), we assume that the researcher has an estimator $\widehat{\beta}$ that remains consistent for $\beta$ under misspecification of $f_{\sigma}$. More precisely, we will only consider potential true densities $f_{0}$ such that $\widehat{\beta}$ tends in probability to the true value $\beta$ under $f_{0}$. For example, in the fixed-effects model (1), consistent estimates of means and variances can be obtained in the absence of normality.

To map model (1) to the general notation of this section, note that in this case there are no covariates $X$, and the vector of unobservables $U$ is

[TABLE]

The vector $\beta$ is $\beta=(\mu_{\alpha},s_{\alpha}^{2},s_{\varepsilon}^{2})^{\prime}$. The reference distribution for $U$ is a standard multivariate normal, so the reference density $f_{\sigma}$ is known in this case — in other words, the parameter $\sigma$ in $f_{\sigma}$ can be omitted. We assume that the researcher has computed an estimator $\widehat{\beta}$, for example by quasi-maximum likelihood or minimum distance, which remains consistent for $\beta$ when $U$ is not normally distributed.

In certain applications, the reference density depends on some parameters $\sigma$ that cannot be consistently estimated absent parametric assumptions. In Appendix S6, we describe discrete choice and censored regression models that have this structure. In such settings, we assume that the researcher has an estimator $\widehat{\sigma}$ that tends in probability to some $\sigma_{*}$ under $f_{0}$. Unlike $\beta$, the parameter $\sigma_{*}$ is a model-specific “pseudo-true value” that is not assumed to have generated the data. However, in our leading example of model (1), as well as in the model’s generalizations that we study in our empirical illustrations in Section 4, the references to $\widehat{\sigma}$ and $\sigma_{*}$ can be omitted from all subsequent statements and derivations.

Given $\widehat{\beta}$, $\widehat{\sigma}$, a sample $\{Y_{i},X_{i},\,i=1,...,n\}$ from $(Y,X)$, and the parametric density $f_{\sigma}$, a model-based estimator of $\overline{\delta}$ is

[TABLE]

where, with some abuse of notation, the expectation with respect to $f_{\widehat{\sigma}}$ is computed only over $U$. When not available in closed form, this estimator can be computed by numerical integration or simulation under the parametric density $f_{\widehat{\sigma}}$. It is easy to see that, under standard conditions, $\widehat{\delta}^{\mathrm{M}}$ is consistent for $\overline{\delta}$ under correct specification; that is, when $f_{\sigma_{*}}$ is the true density of $U\,|\,X$.

To construct a posterior estimator, consider the posterior density $p_{\beta,\sigma}$ of $U\,|\,Y,X$. This posterior density is computed using Bayes rule, based on the prior $f_{\sigma}$ on $U\,|\,X$ and the likelihood of $Y\,|\,U,X$ implied by $g_{\beta}$. Formally, let ${\cal{U}}(y,x,\beta)=\{u\,:\,y=g_{\beta}(u,x)\}$. We define, whenever the denominator is non-zero,

[TABLE]

We will compute $p_{\beta,\sigma}$ analytically in our examples. In Appendix S5 we describe a simulation-based computational approach when an analytical expression is not available. We define the posterior average effect (PAE) as the posterior estimator

[TABLE]

where, again, the expectation is only taken over $U$. Under standard regularity conditions, it is easy to see that, like $\widehat{\delta}^{\mathrm{M}}$, the PAE $\widehat{\delta}^{\mathrm{P}}$ is consistent for $\overline{\delta}$ under correct specification.

From a Bayesian perspective, $\widehat{\delta}^{\mathrm{P}}$ is a natural estimator to consider when $\beta$ and $\sigma$ are known. Indeed, $\widehat{\delta}^{\mathrm{P}}$ is then the posterior mean of $\frac{1}{n}\sum_{i=1}^{n}\delta_{\beta}(U_{i},X_{i})$, where the prior on $U_{i}$ is $f_{\sigma}$, independent across $i$. An alternative Bayesian interpretation is obtained by specifying a nonparametric prior on $f_{0}$, and computing the posterior mean of $\overline{\delta}$ under this prior, as we discuss in Appendix S5 in the case where $U$ has finite support. However, a frequentist justification for $\widehat{\delta}^{\mathrm{P}}$ appears to be lacking in the literature. Indeed, under correct specification of $f_{\sigma}$, both estimators $\widehat{\delta}^{\mathrm{P}}$ and $\widehat{\delta}^{\mathrm{M}}$ are consistent, and, as we pointed out in the previous section, $\widehat{\delta}^{\mathrm{P}}$ may have a higher variance than $\widehat{\delta}^{\mathrm{M}}$. The key difference between model-based and posterior estimators is that $\widehat{\delta}^{\mathrm{P}}$ is conditional on the observation sample. An intuitive rationale for the conditioning is the recognition that realizations $Y_{i}$ may be informative about the values of the unknown $U_{i}$’s. We next formalize this intuition in a framework that accounts for specification error.

3.2 Neighborhoods, estimators, and worst-case specification error

Let $P(\beta,f_{0})$ denote the true density of $(Y,U,X)$, where as before we omit the reference to the marginal density of $X$ for conciseness. We assume that, under $P(\beta,f_{0})$, $\widehat{\beta}$ is consistent for the true $\beta$, and $\widehat{\sigma}$ is consistent for a model-specific “pseudo-true” value $\sigma_{*}$, where $\mathbb{E}_{P(\beta,f_{0})}[\psi_{\beta,\sigma_{*}}(Y,X)]=0$ for some moment function $\psi$. For example, $\widehat{\beta}$ and $\widehat{\sigma}$ may be the method-of-moments estimators that solve $\sum_{i=1}^{n}\psi_{\widehat{\beta},\widehat{\sigma}}(Y_{i},X_{i})=0$. In models with no $\sigma$ parameters, such as model (1) and its generalizations, we only assume that $\widehat{\beta}$ is consistent for $\beta$, and that $\mathbb{E}_{P(\beta,f_{0})}[\psi_{\beta}(Y,X)]=0$ for some $\psi$. Throughout, we take the estimators $\widehat{\beta}$ (and possibly $\widehat{\sigma}$), and the moment function $\psi$, as given. In particular, we do not address the question of optimal estimation of $\beta$ under misspecification.

Given a distance measure $d$ and a scalar $\epsilon\geq 0$, we define the following neighborhood of the reference density $f_{\sigma}$:

[TABLE]

This neighborhood consists of densities of $U\,|\,X$ that are at most $\epsilon$ away from $f_{\sigma_{*}}$, and under which $\widehat{\beta}$ and $\widehat{\sigma}$ converge asymptotically to $\beta$ and $\sigma_{*}$, respectively. The case $\epsilon=0$ corresponds to correct specification of the reference density, whereas $\epsilon>0$ corresponds to misspecification.

For ease of notation we omit the dependence of $\Gamma_{\epsilon}$ on $\beta$, $\sigma_{*}$, and $\psi$, all of which we consider fixed and given in this section. Indeed, we assume that the researcher has chosen an estimator $\widehat{\beta}$, and, depending on the setting, an estimator $\widehat{\sigma}$ — our theory is silent about where these choices come from — and that she has already observed their realized values in a large sample. The moment function $\psi$ is determined by this choice of estimators. Moreover, in large samples, the population values $\beta$ and $\sigma_{*}$ are arbitrarily close to the observed values $\widehat{\beta}$ and $\widehat{\sigma}$. In our setup, we only consider densities of unobservables $f_{0}$ that are consistent with those values, in the sense that the moment restriction $\mathbb{E}_{P(\beta,f_{0})}[\psi_{\beta,\sigma_{*}}(Y,X)]=0$ holds. This large-sample logic is consistent with our focus on specification error; see (16) below.

Note that the same logic might suggest imposing that other features of the joint population distribution of the data $(Y,X)$, such as means, covariances, higher-order moments, or even the entire distribution, be kept constant for all $f_{0}\in\Gamma_{\epsilon}$. Restricting neighborhoods in this way does not affect the results in this section, because those are valid for all possible $\psi$, and one could thus impose additional moment restrictions on $f_{0}$.

Let us denote the supports of $X$ and $U$ as ${\cal{X}}$ and ${\cal{U}}$, respectively. We assume that $d$ is a $\phi$-divergence of the form

[TABLE]

where $\phi$ is a convex function that satisfies $\phi(1)=0$ and $\phi^{\prime\prime}(1)>0$. This family contains as special cases the $\chi^{2}$ divergence (averaged over $X$), the Kullback-Leibler divergence, the Hellinger distance, and more generally the members of the Cressie-Read family of divergences (Cressie and Read, 1984). It is commonly used to measure misspecification, see Andrews et al. (2020) and Christensen and Connault (2019) for recent examples.

We focus on asymptotically linear estimators of $\overline{\delta}$ that satisfy, for a scalar non-stochastic function $\gamma$ and as $n$ tends to infinity,

[TABLE]

Note that $\widehat{\delta}_{\gamma}$ depends on $\widehat{\beta},\widehat{\sigma}$, but for conciseness we leave the dependence implicit in the notation. Many estimators can be written in this form (see, e.g., Bickel et al., 1993). Given an estimator $\widehat{\delta}_{\gamma}$, we define its $\epsilon$-worst-case specification error as

[TABLE]

We will take the worst-case specification error $b_{\epsilon}(\gamma)$ to be our measure of how well an estimator $\widehat{\delta}_{\gamma}$ performs under misspecification. It quantifies the maximum discrepancy, under any possible $f_{0}$ in the neighborhood $\Gamma_{\epsilon}$, between the probability limit of the estimator and the true parameter value. Under suitable regularity conditions, $\mathbb{E}_{P(\beta,f_{0})}[\gamma_{{\beta},{\sigma}_{*}}(Y,X)]-\mathbb{E}_{f_{0}}[\delta_{\beta}(U,X)]$ in (16) is the asymptotic bias of $\widehat{\delta}_{\gamma}$ under $P(\beta,f_{0})$.

By focusing on the worst-case specification error $b_{\epsilon}(\gamma)$, we abstract from other sources of estimation error. Importantly, we do not account for sampling variability. In Bonhomme and Weidner (2018), we study an alternative approach that consists in minimizing worst-case mean squared error under a local asymptotic — i.e., as $\epsilon$ tends to zero, $n$ tends to infinity, and $\epsilon n$ tends to a positive constant. Applying this approach to the present case gives estimators that have a smaller worst-case mean squared error than PAE in general. However, unlike PAE, minimum-MSE estimators depend on $\epsilon$, as we will discuss Subsection 3.5 below. Relative to such estimators, PAE do not require the researcher to take a stand on the degree of misspecification $\epsilon$, and they are easy to implement.

3.3 Result under small-$\epsilon$ misspecification

Before stating our first main result, we first characterize the worst-case specification error $b_{\epsilon}(\gamma)$ of estimators $\widehat{\delta}_{\gamma}$ for small $\epsilon$. For conciseness, in the remainder of this section we suppress the reference to $\beta,\sigma_{*}$ from the notation, and we denote as $\mathbb{E}_{*}$ and $\mathop{\mathrm{V}ar}_{*}$ expectations and variances that are taken under the reference model $P(\beta,f_{\sigma_{*}})$. All proofs are in Appendix S1.

Lemma 1**.**

Let \widetilde{\psi}(y,x)=\psi(y,x)-\mathbb{E}_{*}\left[\psi(Y,X)\big{|}X=x\right]. Suppose that one of the following conditions holds:

• (i)

$\phi(1)=0$, $\phi(r)$ is four times continuously differentiable with $\phi^{\prime\prime}(r)>0$ for all $r>0$, $\mathbb{E}_{*}[\psi(Y,X)]=0$, \mathbb{E}_{*}\big{[}\widetilde{\psi}(Y,X)\,\widetilde{\psi}(Y,X)^{\prime}\big{]}>0, and $\left|\gamma(y,x)\right|$, $\left|\delta(u,x)\right|$, $\left|\psi(y,x)\right|$ are bounded over the domain of $Y$, $U$, $X$.

• (ii)

Condition (ii) of Lemma S1 in Appendix S1 holds (this alternative condition allows for unbounded $\gamma$, $\delta$, $\psi$, but at the cost of stronger assumptions on $\phi(r)$).

Then, as $\epsilon$ tends to zero we have

[TABLE]

where \lambda=\left\{\mathbb{E}_{*}\big{[}\widetilde{\psi}(Y,X)\,\widetilde{\psi}(Y,X)^{\prime}\big{]}\right\}^{-1}\mathbb{E}_{*}\left[\left({\gamma}(Y,X)-{\delta}(U,X)\right)\widetilde{\psi}(Y,X)\right].

To derive the formula for the worst-case specification error in Lemma 1, we maximize the specification error with respect to $f_{0}$ subject to three contraints: $f_{0}$ belongs to an $\epsilon$-neighborhood of $f_{*}$, it is such that the moment condition is satisfied at $(\beta,\sigma_{*})$, and it is a density. In part $(i)$ we focus on the case where $\gamma$, $\delta$ and $\psi$ are bounded. This is satisfied, for example, if those functions and $g(u,x)$ are all continuous, and the domain of $U$ and $X$ is bounded. To accommodate situations where supports are unbounded, such as the example of Section 2, in part $(ii)$ we allow for unbounded functions $\gamma$, $\delta$ and $\psi$, which only requires existence of third moments under the reference distribution. To guarantee that $b_{\epsilon}(\gamma)$ is well-defined in the unbounded case, we require a regularization of the function $\phi(r)$ for large values of $r$.

Lemma 1 implies that the small-$\epsilon$ specification error of the PAE is, up to smaller-order terms, proportional to the within-$(Y,X)$ standard deviation of $\delta(U,X)$ under the reference model:

[TABLE]

In the fixed-effects model (1) of teacher quality, the worst-case specification error of the PAE $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ is

[TABLE]

where $T(a,b)=\varphi(a)\int_{0}^{b}\frac{\varphi(az)}{1+z^{2}}dz$ is Owen’s T function (Owen, 1956), and $\varphi$ is the standard normal density. The specification error decreases as the number $J$ of observations per teacher increases, and tends to zero as $J$ tends to infinity and the shrinkage factor $\rho$ tends to one.

The next theorem, which holds for all functions $\gamma(y,x)$, subject to regularity conditions, shows that the PAE has minimum worst-case specification error locally.

Theorem 1**.**

Suppose that the conditions of Lemma 1 hold, and let

[TABLE]

Then, as $\epsilon$ tends to zero we have

[TABLE]

3.4 Result under fixed-$\epsilon$ misspecification

To show our second main result, let us now focus on the case $\phi(t)=\frac{1}{2}(t-1)^{2}$; that is, we choose the distance measure $d(f_{0},f_{\sigma})$ to be the Pearson $\chi^{2}$ divergence. For this quadratic distance measure, we show that PAE satisfy a fixed-$\epsilon$ optimality result, which is valid for all values of $\epsilon$ that are smaller than

[TABLE]

where $\gamma^{\mathrm{P}}(y,x)$ is given by (17).

Theorem 2**.**

Assume that $\mathbb{E}_{*}[\psi(Y,X)]=0$, $\phi(t)=\frac{1}{2}(t-1)^{2}$, and that $\gamma(Y,X)$ and $\delta(U,X)$ have finite second moments under the reference model. Then, for $0<\epsilon\leq\overline{\epsilon}$, we have

[TABLE]

In Theorem 2 we show that $\gamma^{\mathrm{P}}$ is an exact minimizer of the function $b_{\epsilon}(\gamma)$. This is in contrast with Theorem 1, where we relied on a small-$\epsilon$ approximation. The condition $\epsilon\leq\overline{\epsilon}$ guarantees that, for $\gamma=\gamma^{\mathrm{P}}$, the constraint $f_{0}(u\,|\,x)\geq 0$ is non-binding in the optimization problem over $f_{0}$ in (16), implying that the problem has a simple analytic solution. Although, in many settings such as model (1), the parameter of interest $\overline{\delta}$ is not consistently estimable under our assumptions, Theorem 2 shows that PAE achieve the smallest possible worst-case specification error when the true distribution $f_{0}$ lies sufficiently close to the reference distribution $f_{\sigma_{*}}$, as measured according to the $\chi^{2}$ divergence.

If the distance measure $d(f_{0},f_{\sigma})$ is not a $\chi^{2}$-divergence, or if $\epsilon>\overline{\epsilon}$, then $\gamma^{\mathrm{P}}$ is not the exact minimizer of worst-case specification error $b_{\epsilon}(\gamma)$. Moreover, in such cases the estimator with minimum worst-case specification error depends on $\epsilon$ in general. However, one can still establish a fixed-$\epsilon$ bound on worst-case specification error, as the next result shows.

Theorem 3**.**

Let $\gamma^{\mathrm{P}}$ be as in (17), and assume that $\phi(r)$ is convex with $\phi(1)=0$. Then, for all $\epsilon>0$,

[TABLE]

In Theorem 3 we establish a fixed-$\epsilon$ bound on the worst-case specification error of PAE, which holds for all $\epsilon>0$ and all $\phi$-divergences such that $\phi$ is convex with $\phi(1)=0$. The infimum is taken over all possible functions $\gamma(y,x)$, subject to measurability conditions, which we implicitly assume throughout the paper. Although $\widehat{\delta}^{\mathrm{P}}$ may not minimize worst-case specification error for finite $\epsilon$, Theorem 3 shows that its worst-case specification error is never larger than twice the minimum worst-case specification error. In addition, the factor two in Theorem 3 cannot be improved upon in general, as we show in Appendix S5 in the context of a simple binary choice model.

3.5 Discussion

In this subsection, we discuss several features and implications of our main results given by Theorems 1 and 2.

Uniqueness.

In the absence of covariates and for known parameters $\beta$, $\sigma_{*}$, the proof of Theorem 1 shows that $\gamma^{\mathrm{P}}$ is the unique minimizer of the first-order worst-case specification error. Likewise, $\gamma^{\mathrm{P}}$ is also unique in Theorem 2. More generally, if covariates are present and the parameters $\beta$, $\sigma_{*}$ are estimated, then the leading order contribution of $b_{\epsilon}(\gamma)$ is minimized if and only if $\gamma(Y,X)=\gamma^{\mathrm{P}}(Y,X)+\omega(X)+\lambda^{\prime}\psi(Y,X)+o_{P_{*}}(1)$, for some $\lambda$ and $\omega$ such that $\mathbb{E}_{f_{X}}[\omega(X)]=0$ — see part (ii) of Theorem S1 in Appendix S1 for a formal statement. Hence, while the PAE is not the unique minimizer of the local worst-case specification error in this case, any minimizer differs from the PAE by a zero-mean function of $X$ and a linear combination of the moment function $\psi$. In addition, $\widehat{\delta}^{\mathrm{P}}$ has smallest variance within the class of minimum worst-case specification error estimators.

Form of misspecification.

Theorems 1 and 2 rely on specific distance measures, $\chi^{2}$ divergence for the latter and any member of the $\phi$-divergence family for the former. Under other distance measures, the PAE will not have minimum worst-case specification error in general.

Given a distance measure, the theorems are based on nonparametric neighborhoods that consist of unrestricted distributions of $U\,|\,X$, except for the moment conditions that pin down $\beta$ and $\sigma_{*}$. However, if one is willing to make additional assumptions on $f_{0}$ that further restrict the neighborhood, then one can construct estimators that are more robust than $\widehat{\delta}^{\mathrm{P}}$ within a particular class. As an example, consider the fixed-effects model (1). Suppose that, in addition to assuming that $\alpha$, $\varepsilon_{1}$, …, $\varepsilon_{J}$ are mutually uncorrelated, the researcher is willing to assume that they are fully independent. In that case, the distribution of $\alpha$ can be consistently estimated under suitable regularity conditions, provided $J\geq 2$ (Kotlarski, 1967, Li and Vuong, 1998). However, the PAE in (5) is inconsistent for fixed $J$ as $n$ tends to infinity. As a consequence, the PAE does not minimize worst-case specification error in a semi-parametric neighborhood that consists of distributions with independent marginals.

To elaborate further on this point, consider the coefficient $\overline{\delta}$ in the population regression of $\alpha$ on a covariates vector $W$, see (7). A possible estimator is the coefficient $\widehat{\delta}^{\mathrm{F}E}$ in the regression of the fixed-effects estimates $\overline{Y}_{i}$ on $W_{i}$, see (8). Under correct specification of the reference model, $\widehat{\delta}^{\mathrm{F}E}$ is consistent for $\overline{\delta}$. However, $\widehat{\delta}^{\mathrm{F}E}$ may be inconsistent under the type of misspecification that we allow for, since $\varepsilon_{j}$ and $W$ may be correlated under $f_{0}$. For example, $W$ (e.g., teacher absenteeism) may be influenced by $\alpha$ and factors that correlate with $\varepsilon_{j}$. Theorem 1 shows that, under such misspecification, the PAE $\widehat{\delta}^{\mathrm{P}}$ in (9) has minimum worst-case specification error locally. Nevertheless, if the researcher is confident that $W$ should not enter the outcome equation, and that it is independent of $\varepsilon_{j}$, then it is natural to report the consistent estimator $\widehat{\delta}^{\mathrm{F}E}$.

Posterior informativeness.

Our small-$\epsilon$ calculations can be used to compare the worst-case specification errors of the PAE $\widehat{\delta}^{\mathrm{P}}$ to that of the model-based estimator $\widehat{\delta}^{\mathrm{M}}$. To see this, let $\gamma^{\mathrm{M}}_{\beta,\sigma}(x)=\mathbb{E}_{f_{\sigma}}[\delta_{\beta}(U,X)\,|\,X=x]$. Using Lemma 1, the ratio of the two worst-case specification errors satisfies

[TABLE]

where $v(U,X)$ is the population residual of $(\delta(U,X)-{\gamma}^{\mathrm{M}}(X))$ on $\widetilde{\psi}(Y,X)$, under the parametric reference model; that is, $v(u,x)={\delta}(u,x)-{\gamma}^{\mathrm{M}}(x)+\lambda^{\prime}\widetilde{\psi}(g(u,x),x)$, where all functions are evaluated at $\beta,\sigma_{*}$, and $\lambda$ is as defined in Lemma 1 for the case $\gamma=\gamma^{\mathrm{M}}$. Intuitively, the robustness of $\widehat{\delta}^{\mathrm{P}}$ relative to $\widehat{\delta}^{\mathrm{M}}$ depends on how informative the outcome values $Y_{i}$ are for the latent individual parameters $\delta(U_{i},X_{i})$.

In practice, we will report an empirical counterpart to the small-$\epsilon$ limit of $1-\frac{b_{\epsilon}^{2}(\gamma^{\mathrm{P}})}{b_{\epsilon}^{2}(\gamma^{\mathrm{M}})}$. This quantity can be simply expressed as the $R^{2}$ in the population nonparametric regression of $v(U,X)$ on $Y,X$ under the reference model; that is,

[TABLE]

where with some abuse of notation here $v(U,X)$ denotes the sample residual of $({\delta}_{\widehat{\beta}}(U,X)-{\gamma}_{\widehat{\beta},\widehat{\sigma}}^{\mathrm{M}}(X))$ on $\widetilde{\psi}_{\widehat{\beta},\widehat{\sigma}}(Y,X)$, and expectations and variances are taken with respect to $P(\widehat{\beta},f_{\widehat{\sigma}})$. Using a term from Andrews et al. (2020) — albeit in a different setting — we refer to $R^{2}$ in (20) as a measure of the “informativeness” of the posterior conditioning, and we will report it in our illustrations. As an example, for $\widehat{F}^{\mathrm{P}}_{\alpha}(a)$ in model (1), the informativeness of the posterior conditioning is

[TABLE]

In this case the $R^{2}$ increases with the number $J$ of observations per teacher, and it tends to one as $J$ tends to infinity.

Multi-dimensional PAE.

For simplicity, in this section we have focused on the case where the target parameter $\overline{\delta}$ in (11) is scalar. However, our results can be extended to multi-dimensional parameters. The definition of worst-case specification error in (16) is then modified to

[TABLE]

where $\|\cdot\|$ is a norm over the vector space in which $\gamma(Y,X)$ and $\delta(U,X)$ take values.

If $\|\cdot\|_{*}$ denotes the corresponding dual norm, then we can rewrite $b_{\epsilon}(\gamma)=\sup_{\|v\|_{*}=1}b_{\epsilon}(\gamma,v)$, where b_{\epsilon}(\gamma,v)=\mathop{\mathrm{s}up}_{f_{0}\in\Gamma_{\epsilon}}\,\big{|}\mathbb{E}_{P(\beta,f_{0})}[v^{\prime}\gamma(Y,X)]-\mathbb{E}_{f_{0}}[v^{\prime}\delta(U,X)]\big{|}. Our minimum worst-case specification error results for PAE for scalar $\overline{\delta}$ then apply to $b_{\epsilon}(\gamma,v)$ for every given vector $v$, and the minimum-specification error properties are preserved after taking the supremum over the set of vectors $v$ with $\|v\|_{*}=1$. Thus, in the multi-dimensional case, PAE minimize worst-case specification error for small $\epsilon$ in the sense of Theorem 1, and for fixed $\epsilon$ under the conditions of Theorem 2. In our leading example of Section 2, suppose we are interested in the entire distribution function $F_{\alpha}$. In this case, the average effect is a function indexed by $a$. Taking the supremum norm $\|\cdot\|_{\infty}$ over distribution functions, we obtain that, as an estimator of $F_{\alpha}$, the PAE minimizes worst-case specification error under suitable conditions.

Mean squared error.

While we have shown that PAE minimize worst-case specification error locally under the conditions of Theorem 1, and for fixed $\epsilon$ under the conditions of Theorem 2, PAE generally do not have minimum mean squared error (MSE). To see this, let us assume that $\beta$ and $\sigma_{*}$ are known. In a local asymptotic framework where $n$ tends to infinity, $\epsilon$ tends to zero, and $n\epsilon$ tends to a positive constant, and under suitable regularity conditions, we show in Appendix S5 that the estimator with minimum worst-case MSE is given by

[TABLE]

which is a linear combination between the model-based estimator and the PAE. The model-based estimator $\widehat{\delta}^{\mathrm{M}}$, which has the smallest asymptotic variance, will be preferred when $\epsilon$ is small relative to $1/n$, while the PAE, which has smallest specification error, will be preferred when $\epsilon$ is large relative to $1/n$. However, in order to implement such estimators $\widehat{\delta}^{\mathrm{M}MSE}$ that minimize worst-case MSE, knowledge of $\epsilon$ is required. See Bonhomme and Weidner (2018) for an approach to minimum-MSE estimation.

4 Simulations and empirical illustrations

In this section, we study two empirical applications: we estimate the distribution of income neighborhood effects in the US, and the distributions of permanent and transitory earnings components in the PSID. We start the section by summarizing the results of a Monte Carlo simulation exercise, in samples generated from various specifications of model (1).

4.1 Monte Carlo simulation: summary of results

While Theorems 1 and 2 show that PAE minimize worst-case specification error under small-$\epsilon$ and fixed-$\epsilon$ misspecification, respectively, they are silent about other forms of estimation error. In Appendix S4 we report the results of a Monte Carlo simulation exercise, where we compare the performance of PAE and other estimators in finite sample in the fixed-effects model (1), for various specifications. Here we briefly summarize the results from the simulation exercise.

We compare the performance of four estimators: the fixed-effects estimator given by (2), the PAE given by (5), the model-based estimator given by (6), and a nonparametric kernel deconvolution estimator with normal errors (Stefanski and Carroll, 1990). We analyze two sets of data generating processes. When the reference normal distribution for $\alpha_{i}$ is correctly specified, the model-based estimator performs best, as expected. We find that, while the PAE has both larger bias and variance than the model-based estimator in this case, it is less biased and less variable than both the nonparametric deconvolution estimator and the fixed-effects estimator, especially when the number of measurements $J$ is small (see Appendix Figure S1).

We next turn to data generating processes where $\alpha_{i}$ is not normal, drawn from a skewed Beta distribution. We find that the model-based estimator is substantially biased in this case. The nonparametric deconvolution estimator has smallest bias when errors are normally distributed, but it is heavily biased when errors are non-normal. By contrast, although it has no consistency guarantees in these settings, the PAE tends to perform comparatively well in all situations, for bias and variance (see Appendix Figure S2).

Overall, the simulations complement our theory by highlighting that, beyond specification error, other sources of estimation error matter in practice. Under correct specification of the reference distribution, the model-based estimator should be preferred. At the same time, our results suggest that, at least in the particular settings we focus on, the performance of the PAE appears less sensitive to misspecification than those of the model-based and nonparametric deconvolution estimators. Moreover, we find that the robustness gains provided by the PAE depend on the signal-to-noise ratio and the informativeness of the posterior conditioning. We provide details on the simulations in Appendix S4.

4.2 Neighborhood effects

In this subsection and the next, we revisit two applications of models with latent variables. In our first illustration, we focus on a model of neighborhood effects following Chetty and Hendren (2017), using data for the US that these authors made public. In our second illustration, we study a permanent-transitory model of income dynamics (Hall and Mishkin, 1982, Blundell et al., 2008) using the PSID. In both cases, we rely on a normal reference specification and assess how and by how much the posterior conditioning informs the estimates of the parameters of interest.

Here we start with estimates of neighborhood (or “place”) effects reported in Chetty and Hendren (2017, CH hereafter). Those were obtained using individuals who moved between different commuting zones at different ages. The outcome variable that we focus on is the causal estimate of the income rank at age 26 of a child whose parents are at the 25 percentile of the income distribution. This is CH’s preferred measure of place effect.

CH report an estimate of the variance of neighborhood effects, corrected for noise. In addition, they report individual predictors. Here we are interested in documenting the entire distribution of place effects. To do so, we consider the model $\widehat{\mu}_{c}=\mu_{c}+\overline{\varepsilon}_{c}$, for each commuting zone $c$, where $\widehat{\mu}_{c}$ is a neighborhood-specific fixed-effects reported by CH, $\mu_{c}$ is the true effect of neighborhood $c$, and $\overline{\varepsilon}_{c}$ is additive estimation noise. CH also report estimates $\widehat{s}_{c}^{2}$ of the variances of $\overline{\varepsilon}_{c}$ for every $c$. When weighted by population, the fixed-effects estimates $\widehat{\mu}_{c}$ have mean zero. We treat neighborhoods as independent observations. The statistics we use for calculations are available at: https://opportunityinsights.org/paper/neighborhoodsii/. Given the aggregate data at hand, we necessarily need to assume that estimates $\widehat{\mu}_{c}$ are independent across neighborhoods $c$, although this might be restrictive in this setting.

We first estimate the variance of place effects $\mu_{c}$, following CH. We trim the top 1% percentile of $\widehat{s}_{c}^{2}$, and weigh all results by population weights. While this differs slightly from CH’s approach, which is based on $1/\widehat{s}_{c}^{2}$ precision weights and no trimming, we replicated the analysis using precision weights in the un-trimmed sample and found similar results. We have information about place effects in $C=590$ commuting zones $c$ in our sample, compared to 595 in the sample without trimming. We estimate a sizable variance of neighborhood fixed-effects: $\mathop{\mathrm{V}ar}(\widehat{\mu}_{c})=.077$. In turn, the mean of $\widehat{s}_{c}^{2}$ weighted by population is $\widehat{s}_{\overline{\varepsilon}}^{2}=.047$. Given those, we estimate the variance of place effects as $\widehat{s}_{\mu}^{2}=\mathop{\mathrm{V}ar}(\widehat{\mu}_{c})-\widehat{s}_{\overline{\varepsilon}}^{2}=.030$. In this setting, the shrinkage factor $\widehat{\rho}_{c}=\widehat{s}_{\mu}^{2}/(\widehat{s}_{\mu}^{2}+\widehat{s}_{c}^{2})$ exhibits substantial heterogeneity across commuting zones. Indeed, the mean of $\widehat{\rho}_{c}$ is .62, and its 10% and 90% percentiles are .21 and .93, respectively.

We use a normal with zero mean and variance $\widehat{s}_{\mu}^{2}$ as a prior for $\mu_{c}$. Then, we estimate the distribution function of neighborhood effects ${\mu}_{c}$ using the PAE given by (5); that is,

[TABLE]

where $\pi_{c}$ are population weights. In addition, in order to ease the visualization of the results, we will also report estimates of densities, which are the derivatives of the PAE of distribution functions. Note that the density of $\mu$ at $a$ can be approximated for arbitrarily small $h>0$ by the expectation of $\boldsymbol{1}\{|\mu-a|/h\}/2h$. Taking the limit of the corresponding PAE as $h$ tends to zero gives the derivative of $\widehat{F}^{\mathrm{P}}_{\mu}$ at $a$. We thus expect derivatives of PAE of distribution functions to enjoy similar minimum-worst-case specification error properties as PAE, but we do not formalize the required assumptions here.

In the top panel of Figure 1, we report several estimates of distribution functions. In the bottom panel, we report the corresponding density estimates. In the left graphs, we show nonparametric kernel estimates of the distribution function (respectively, density) of the fixed-effects $\widehat{\mu}_{c}$, weighted by population (in solid), together with the best-fitting normal (in dashed). The graphs show substantial non-normality of the fixed-effects estimates. In particular, the large variance appears to be driven by some large positive and negative estimates $\widehat{\mu}_{c}$. In the right graphs, we report the PAE $\widehat{F}^{\mathrm{P}}_{\mu}$ of the distribution function of true place effects $\mu_{c}$, with the associated density (in solid). In addition, we show the normal prior, with zero mean and variance $\widehat{s}_{\mu}^{2}$ (in dashed). The posterior distribution of neighborhood effects differs from the normal prior, although the two estimators have the same variance by construction. In comparison, neighborhood-specific empirical Bayes estimates have a substantially lower dispersion. In Appendix Figure S5 we report an estimate of their distribution function $\widehat{F}^{\mathrm{P}M}_{\mu}$ and associated density. While $\widehat{s}_{\mu}^{2}=.030$ and the variance associated with $\widehat{F}^{\mathrm{P}}_{\mu}$ is $.030$, the variance of the empirical Bayes estimates is only $.010$. In addition, a specification test that compares model-based estimator and PAE, which we describe in Appendix S5, suggests that these differences are statistically significant. Indeed, assuming independence across commuting zones, we obtain p-values below .01 at all deciles except the bottom two.

To assess how likely it is that the posterior estimator approximates the shape of the distribution of true neighborhood effects, we next perform two different exercises, based on a simulation and on numerical calculations motivated by our theory. We start with a Monte Carlo simulation, where $\mu_{c}$, for $c=1,...,C_{\mathrm{s}im}$, are log-normally distributed with zero mean and variance $\widehat{s}_{\mu}^{2}$, and $\overline{\varepsilon}_{c}$ are normally distributed independent of $\mu_{c}$ with zero mean. We consider three scenarios for the noise variances $\widehat{s}_{c}^{2}$: the estimates from CH, one-third of those values, and one-tenth of those values. In this exercise we again weigh by population. We show the results for $C_{\mathrm{s}im}=100,000$ simulated neighborhoods. In the left graphs of Figure 2 we see that, when the noise variances are the ones from the data, the posterior density is more skewed than the normal, yet the posterior shape is quite different from the true log-normal distribution of $\mu_{c}$. When reducing the noise variances in the middle and right graphs, the posterior distribution function and density estimates get closer to the log-normal ones. In the right graphs, where the shrinkage factor is .90 on average (as opposed to .62 in the data), the posterior distribution function and density approximate the highly non-normal shape of the true distribution of neighborhood effects very well.

We next turn to our posterior informativeness measure, which is given by equation (21). Note the $R^{2}$ coefficient varies along the distribution. We find that the weighted average $R^{2}$ across values of $a$ is 28%, where we weigh across cutoff values $a$ by the reference distribution for $\alpha$. This value is consistent with the message of Figure 2, since it suggests that, while the posterior conditioning informs the shape of the distribution of neighborhood effects, the signal-to-noise ratio is not high enough to be confident about the exact shape.

Lastly, we perform two additional exercises as robustness checks. Firstly, we incorporate the mean income $\overline{y}_{c}$ of permanent residents in county $c$ at the 25% percentile as a covariate. CH rely on information on permanent residents’ income to improve the accuracy of individual predictions. Here we use it to refine the reference distribution and to improve the estimation of the distribution of neighborhood effects. Specifically, our reference model for $\mu_{c}$ is then a correlated random-effects specification, where the mean depends on $\overline{y}_{c}$ linearly. Appendix Figure S6 shows small differences with our baseline estimates. Secondly, we re-do our main analysis at the county level, instead of the commuting zone level. In that case, the signal-to-noise ratio is lower, our posterior informativeness $R^{2}$ measure is 17% on average, and Appendix Figure S7 shows that the normal prior and the posterior distributions are closer to each other than in the case of commuting zones.

4.3 Income dynamics

In this subsection, we consider the following permanent-transitory model of household log-income,

[TABLE]

where $\varepsilon_{it}$ and $V_{it}$ are independent at all lags and leads, and independent of $\eta_{i0}$. This process is commonly used as an input for life-cycle consumption/savings models. Researchers often estimate covariances in a first step using minimum distance, and then impose a normality assumption for further analysis. However, there is increasing evidence that income components are not normally distributed. Instead of using a more flexible model — as has been done by Carlton and Hall (1978) and a large subsequent literature — here we compute posterior average effects. The advantages of this approach are that no additional assumptions are needed, and that implementation is straightforward.

We focus on six recent waves of the PSID 1999-2009 (every other year), see Blundell et al. (2016) for a description of the data. We use the same sample selection as in Arellano et al. (2017), and work with a balanced panel of $n=792$ households over $T=6$ periods. $Y_{it}$ are residuals of log total pre-tax household labor earnings on a set of demographics, which include cohort interacted with education categories for both household members, race, state, and large-city dummies, a family size indicator, number of kids, a dummy for income recipient other than husband and wife, and a dummy for kids out of the household. Our aim is to estimate the distributions of $\eta_{it}$ and $\varepsilon_{it}$. To do so, we compare normal model-based estimates with posterior estimates, by plotting distribution functions as well as the implied densities. The model’s structure is similar to that of the fixed-effects model (1), and analytical expressions for posterior estimators are easy to derive.

In the left graphs of Figure 3, we show the distribution of the permanent component $\eta_{it}$. In the right graphs, we show the distribution of the transitory component $\varepsilon_{it}$. We show PAE in solid, and model-based estimators in dashed. In the top panel we report estimates of distribution functions, and in the bottom panel we report the implied density estimates. The estimates show mild deviation from Gaussianity for the permanent component, and stronger evidence of non-Gaussianity for the transitory component. In particular, the latter shows excess kurtosis (i.e., “peakedness”) relative to the normal.

Several papers have already documented the presence of excess kurtosis in income components, particularly in transitory innovations, using parametric or semi-parametric methods. The estimates in Figure 3 share some qualitative similarities with recent findings in the literature. For example, the estimates of a flexible non-normal and non-linear model in Arellano et al. (2017, Figure 3) are quite similar to the PAE estimates in Figure 3 for permanent components. At the same time, their estimates of the distribution of transitory components show substantially more pronounced non-Gaussianity and excess kurtosis relative to PAE. This finding is in agreement with our posterior informativeness measure R2, which is 12% on average along the distribution for the permanent component, and 8% on average for the transitory component. This degree of informativeness suggests that posterior estimates may suffer from substantial specification error when the reference distribution is misspecified.

Overall, these empirical illustrations give two examples where, starting from a normal prior, the posterior conditioning is informative about the true unknown distributions. In both settings, PAE are not normal. Yet, as indicated by the $R^{2}$ values we report, the signal-to-noise ratios are not high enough to be certain about the exact shapes of the distributions of interest, thus motivating further analyses using non-normal specifications. PAE should be useful in other environments where model (1) and its extensions are widely used, for example in teacher value-added applications, where the signal-to-noise ratio is driven by the number of observations per teacher. Moreover, PAE are also applicable to other — nonlinear — econometric models, as we describe in the next section.

5 Complements and extensions

In this section, we outline several complements and extensions that we analyze in detail in the appendix.

5.1 PAE in other models

PAE are applicable to a variety of settings. In many econometric models, semi-parametric estimators — i.e., robust to distributional assumptions on unobservables — of $\beta$ parameters are available; see Powell (1994) for examples. In such models, PAE provide estimators of average effects that enjoy robustness properties when parametric assumptions are violated. In Appendix S6 we study static binary and ordered choice models, censored regression models, and panel data binary choice models. We also show how the White (1980) formula for robust standard errors in linear regression can be interpreted as a PAE.

5.2 Confidence intervals and specification test

Under correct specification of the reference model, it is easy to derive the asymptotic distributions of $\widehat{\delta}^{\mathrm{M}}$ and $\widehat{\delta}^{\mathrm{P}}$ using standard arguments. Moreover, under local misspecification, confidence intervals that account for both model uncertainty and sampling uncertainty can be constructed following Armstrong and Kolesár (2018) and Bonhomme and Weidner (2018). However, such confidence intervals require the researcher to set a value for the degree of misspecification $\epsilon$. In Appendix S5, we provide details on confidence intervals calculations. In addition, we explain how to construct a specification test of the reference model based on the difference $\widehat{\delta}^{\mathrm{P}}-\widehat{\delta}^{\mathrm{M}}$.

5.3 Robustness in prediction

In applications such as the fixed-effects model (1) of teacher quality, researchers are often interested in predicting the quality $\alpha_{i}$ of teacher $i$. Although our focus in this paper is on the estimation of population averages, it is interesting to see how different predictors perform under misspecification of the reference distribution. It is well-known that EB estimators minimize mean squared prediction error when the normal reference model is correctly specified. However, when normality fails, the best predictor is a different posterior mean, which does not generally coincide with the EB estimate based on a normal prior. Intuitively, conditioning on nonlinear functions of the data may improve prediction accuracy.

In Appendix S3 we use our framework — applied to worst-case mean squared prediction error instead of worst-case specification error of a sample average — to provide results on the robustness of EB estimators in the presence of misspecification. We show that EB estimators have minimum worst-case mean squared prediction error, up to smaller-order terms, under local deviations from normality. In addition, we derive a fixed-$\epsilon$, non-local risk bound in the spirit of Theorem 3.

6 Conclusion

Posterior averages are commonly used to predict individual parameters, such as teacher quality or neighborhood effects, and they play a central role in Bayesian and empirical Bayes approaches. In this paper, we have provided a frequentist justification for posterior conditioning when the goal of the researcher is to estimate a population average quantity. We have shown that posterior average effects (PAE) have minimum worst-case specification error under various forms of misspecification of parametric assumptions. PAE are simple to implement, and our analysis provides a rationale for reporting them in applications alongside other parametric and semi-parametric estimators, as well as a simple way to assess the informativeness of the posterior conditioning. As an example, Arnold et al. (2020) recently reported PAE to document judge heterogeneity in the context of bail decisions. While we have used a linear fixed-effects model as a running example due to its popularity, there are other possible applications, some of which we discuss in the appendix.

Acknowledgments.

We would like to thank two anonymous referees, Manuel Arellano, Tim Armstrong, Raj Chetty, Tim Christensen, Nathan Hendren, Peter Hull, Max Kasy, Derek Neal, Jesse Shapiro, Xiaoxia Shi, Danny Yagan, and audiences at various places for comments. Bonhomme acknowledges support from the NSF, Grant SES-1658920. Weidner acknowledges support from the Economic and Social Research Council through the ESRC Centre for Microdata Methods and Practice grant RES-589-28-0001 and from the European Research Council grants ERC-2014-CoG-646917-ROMIA and ERC-2018-CoG-819086-PANEDA.

Appendix S1 Proofs of Lemma 1, and Theorem 1 and 2

The following is an extended version of Lemma 1 and Theorem 1 in the main text, which also covers the case of unbounded functions $\gamma_{{\beta},{\sigma}_{*}}(y,x)$, $\delta_{\beta}(u,x)$ and $\psi_{\beta,\sigma_{*}}(y,x)$. In addition, we make explicit again the dependence on $\beta$ and $\sigma_{*}$, which we suppressed in the main text.

Lemma S1**.**

In addition to defining \widetilde{\psi}(y,x)=\psi(y,x)-\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\psi(Y,X)\big{|}X=x\right], let \widetilde{\gamma}(y,x)=\gamma(y,x)-\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\gamma(Y,X)\big{|}X=x\right] and \widetilde{\delta}(u,x)=\delta(u,x)-\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\delta(U,X)\big{|}X=x\right]. Suppose that $\phi(r)=\overline{\phi}(r)+\nu\,(r-1)^{2}$, with $\nu\geq 0$, and a function $\overline{\phi}(r)$ that is four times continuously differentiable with $\overline{\phi}(1)=0$ and $\overline{\phi}^{\prime\prime}(r)>0$, for all $r\in(0,\infty)$. Assume $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\psi_{{\beta},{\sigma}_{*}}(Y,X)=0$ and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)\,\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)^{\prime}\right]>0$. Furthermore, assume that one of the following holds:

• (i)

$\nu=0$, and the functions $\left|\gamma_{{\beta},{\sigma}_{*}}(y,x)\right|$, $\left|\delta_{\beta}(u,x)\right|$ and $\left|\psi_{{\beta},{\sigma}_{*}}(y,x)\right|$ are bounded over the domain of $Y$, $U$, $X$.

• (ii)

$\nu>0$, and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left|\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right|^{3}<\infty$, and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left|\psi_{{\beta},{\sigma}_{*}}(Y,X)\right|^{3}<\infty$.

Then, as $\epsilon\rightarrow 0$ we have

[TABLE]

where

[TABLE]

Theorem S1**.**

Suppose that the conditions of Lemma S1 hold, and let

[TABLE]

Then the following results hold as $\epsilon$ tends to zero.

• (i)

We have

[TABLE]

• (ii)

If we have $b_{\epsilon}(\gamma)=b_{\epsilon}(\gamma_{\beta,\sigma_{*}}^{\mathrm{P}})+o(\epsilon^{1/2}),$ then there exist $\lambda\in\mathbb{R}^{\dim\psi}$ and a function $\omega:{\cal X}\rightarrow\mathbb{R}$ with $\mathbb{E}_{f_{X}}[\omega(X)]=0$ such that

[TABLE]

Notice that Theorem S1 in the main text is a special case of part (i) of Theorem S1. Part (ii) of Theorem S1 is discussed in Subsection 3.5 of the main text. The proof of Theorem S1 provides explicit expressions for $\lambda$ and $\omega(X)$ that appear in part (ii), namely $\lambda$ is the same as in the last line of Lemma S1, and \omega(x)=\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\big{[}\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)-\lambda^{\prime}\,\psi_{\beta,\sigma_{*}}(Y,X)\big{|}X=x\big{]}-\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\big{[}\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)-\lambda^{\prime}\,\psi_{\beta,\sigma_{*}}(Y,X)\big{]}.

S1.1 Proof of Lemma S1 (containing Lemma 1 as a special case)

We first introduce some additional notation and establish some helpful intermediate results. We write ${\cal B}$ and ${\cal S}$ for the set of possible values of the parameters $\beta$ and $\sigma$, respectively. Lemma S1 is for given values $\beta\in{\cal B}$ and $\sigma_{*}\in{\cal S}$, and given functions $\gamma_{\beta,\sigma_{*}}(y,x)$, $\delta_{\beta}(u,x)$, $\psi_{\beta,\sigma_{*}}(y,x)$, and those values and functions are also taken as given in following two intermediate lemmas. Remember also that $\Gamma_{\epsilon}$ depends on the function $\phi:[0,\infty)\rightarrow\mathbb{R}\cup\{\infty\}$, which is assumed to be strictly convex in Lemma S1. We define the corresponding function $\rho:\mathbb{R}\rightarrow\mathbb{R}\cup\{\infty\}$ by

[TABLE]

For $t=\phi^{\prime}(r)$ we have $\rho(t)=r$, that is, for those values of $t$ the function $\rho(t)$ is simply the inverse function of the first derivative $\phi^{\prime}$. For $t<\inf_{r>0}\phi^{\prime}(r)$ we have $\rho(t)=0$, and for $t>\sup_{r>0}\phi^{\prime}(r)$ the value of $\rho(t)$ is defined to be $\infty$. The following lemma provides a characterization of the $\epsilon$-worst-case specification error $b_{\epsilon}(\gamma)$ that was defined in (16).

Lemma S2**.**

Let $\epsilon>0$. Assume that $\phi(r)$ is strictly convex with $\phi(1)=0$. Suppose that for $s\in\{-1,1\}$ and $x\in{\cal X}$ there exists $\lambda^{(1)}_{\beta,\sigma_{*}}(s,x)\in\mathbb{R}$, $\lambda^{(2)}_{\beta,\sigma_{*}}(s)>0$, $\lambda^{(3)}_{\beta,\sigma_{*}}(s)\in\mathbb{R}^{\dim\psi}$ such that

[TABLE]

satisfies

[TABLE]

Then the maximizer ($s=+1$) and minimizer ($s=-1$) of $\mathbb{E}_{P(\beta,f_{0})}\left[\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right]$ over $f_{0}\in\Gamma_{\epsilon}$ are given by

[TABLE]

and for the worst-case absolute specification error we therefore have

[TABLE]

The proof of Lemma S2 is given in Section S2. Notice that for $\phi(r)=r[\log(r)-1]$, when $d(f_{0},f_{\sigma_{*}})$ is the Kullback-Leibler divergence, we have $\rho(t)=\exp(t)$, and the worst case densities $f^{(s)}_{0}(u|x)$ in Lemma S2 are exponentially tilted versions of the reference density $f_{\sigma_{*}}(u|x)$. Lemma S2 shows that, more generally, the required “tilting function” is given by $\rho(t)$.

We impose $\phi(1)=0$ throughout the paper to guarantee that $d(f_{0},f_{\sigma_{*}})\geq 0$ (by an application of Jensen’s inequality). In addition, we now impose the normalization $\phi^{\prime}(1)=0$. This is without loss of generality, because we can always redefine $\phi(r)\mapsto\phi(r)-(r-1)\,\phi^{\prime}(1)$, which has no effect on $d(f_{0},f_{\sigma_{*}})$ and guarantees $\phi^{\prime}(1)=0$ for the redefined function.

The goal of the following lemma is to establish Taylor expansions of $\rho(t)$ and $\phi(\rho(t))$ around $t=0$ of the form

[TABLE]

where the remainder terms are defined by

[TABLE]

Notice that the expansions (S6) are trivially true by definition of $R_{1}(t)$ and $R_{2}(t)$, but the following lemma provides bounds on $R_{1}(t)$ and $R_{2}(t)$, which are useful for the proof of Lemma S1 afterwards.

Lemma S3**.**

For all $r\geq 0$ let $\phi(r)=\overline{\phi}(r)+\nu\,(r-1)^{2}$, for $\nu\geq 0$, and a function $\overline{\phi}:[0,\infty)\rightarrow\mathbb{R}\cup\{\infty\}$ that is four times continuously differentiable with $\overline{\phi}(1)=\overline{\phi}^{\prime}(1)=0$ and $\overline{\phi}^{\prime\prime}(r)>0$, for all $r\in(0,\infty)$. The lemma has two parts:

• (i)

Assume in addition that $\nu=0$. Then, there exist constants $c_{1}>0$, $c_{2}>0$ and $\eta>0$ such that for all $t\in[-\eta,\eta]$ we have

[TABLE]

and the functions $R_{1}(t)$ and $R_{2}(t)$ are continuous within $[-\eta,\eta]$.

• (ii)

Assume in addition that $\nu>0$. Then, there exist constants $c_{1}>0$ and $c_{2}>0$ such that the two inequalities in (S7) hold for all $t\in\mathbb{R}$, and the functions $R_{1}(t)$ and $R_{2}(t)$ are everywhere continuous.

The proof of Lemma S3 is given in Section S2.

Comment:

Part (i) and part (ii) of Lemma S3 give the same approximations of $\rho(t)$ and $\phi(\rho(t))$, but the difference is that in part (i) the result only holds locally in a neighborhood of $t=0$, while in part (ii) the inequalities are established globally for all $t\in\mathbb{R}$. Notice that the result of part (ii) cannot hold under the assumptions of part (i) only, because $\rho(t)$ is equal to infinity for all $t>t_{\sup}$, where $t_{\sup}=\sup_{r\in(0,\infty)}\phi^{\prime}(r)$ can be finite. The regularization $\phi(r)=\overline{\phi}(r)+\nu\,(r-1)^{2}$, with $\nu>0$, guarantees that $\rho(t)$ is finite and well-defined for all $t\in\mathbb{R}$. This property of the regularized $\phi(r)$ is key whenever the moment functions $\gamma$, $\delta$, $\psi$ are unbounded (i.e., for case (ii) of the assumptions of Lemma S1).

Using the intermediate Lemmas S2 and S3 we can now show Lemma S1, which contains Lemma 1 as a special case. In the following proofs we again drop the arguments $\beta$ and $\sigma_{*}$ everywhere for ease notation, and we write $\mathbb{E}_{*}$ and ${\mathrm{V}ar}_{*}$ for expectations and variances under the reference density $P(\beta,f_{\sigma_{*}})$. We also continue to use the normalization $\phi^{\prime}(1)=0$, which is without loss of generality, as explained above.

Proof of Lemma S1.

Additional notation and definitions:

Let $\lambda\in\mathbb{R}^{\dim\psi}$ be as defined in the statement of the lemma, and furthermore define

[TABLE]

For $s\in\{-1,+1\}$ and $\epsilon>0$, let

[TABLE]

with

[TABLE]

Here, we are explicit about the leading order terms (of order $\epsilon^{1/2}$), but the higher order terms (of order $\epsilon$) contain the coefficients $\lambda^{(1)}_{\mathrm{r}em}(s)\in\mathbb{R}$, $\lambda^{(2)}_{\mathrm{r}em}(s)\in\mathbb{R}$, and $\lambda^{(3)}_{\mathrm{r}em}(s)\in\mathbb{R}^{\dim\psi}$, which will only be specified in (S10) below. We can rewrite

[TABLE]

with

[TABLE]

Here, $t(u,x|s)$, $\lambda^{(1)}(s,x)$, $\lambda^{(2)}(s)$, etc, also depend on $\epsilon$, but we do not make this dependence explicit in our notation. Our goal is to apply Lemma S2 with $t_{\beta,\sigma_{*}}(u,x|s)$ in the lemma equal to $t(u,x|s)$ as defined here. However, in order to apply that lemma we need to satisfy the conditions (S5), which in current notation read

[TABLE]

The definition of $t(u,x|s)$ above is already designed to satisfy (S9) to leading order in $\epsilon$, but we still need to find $\lambda_{\mathrm{r}em}^{(1)}(s,x)$, $\lambda_{\mathrm{r}em}^{(2)}(s)$, $\lambda_{\mathrm{r}em}^{(3)}(s)$ such that (S9) holds exactly. Plugging the expansions (S6) into (S9), using the definition of $t(u,x|s)$, as well as \mathbb{E}_{*}\left[t_{(0)}(U,X|s)\big{|}X=x\right]=0, $\mathbb{E}_{*}\left\{[t_{(0)}(U,X|s)]^{2}\right\}=2\,\phi^{\prime\prime}(1)$, and $\mathbb{E}_{*}\psi(Y,X)\,t_{(0)}(U,X|s)=0$, we obtain

[TABLE]

Those conditions can be rewritten as follows

[TABLE]

Thus, as $\epsilon\rightarrow 0$ we have

[TABLE]

Notice that $\lambda_{\mathrm{r}em}^{(1)}(s,x)$, $\lambda_{\mathrm{r}em}^{(2)}(s)$, $\lambda_{\mathrm{r}em}^{(3)}(s)$ also appear implicitly on the right-hand sides of the equations (S10), because $t_{\mathrm{r}em}(u,x|s)$ depends on those parameters, and (S10) is therefore a system of equations for $\lambda_{\mathrm{r}em}^{(1)}(s,x)$, $\lambda_{\mathrm{r}em}^{(2)}(s)$, $\lambda_{\mathrm{r}em}^{(3)}(s)$. Our assumptions guarantee that the system (S10) has a solution for sufficiently small $\epsilon$, as will be explained below for the two different cases distinguished in the lemma.

Proof for case (i):

The assumptions for this case guarantee that $t(u,x|s)$ is uniformly bounded over $u$ and $x$. Part (i) of Lemma S3 guarantees existence of $c_{1}>0$, $c_{2}>0$, $\eta>0$ such that for all $t\in[-\eta,\eta]$ we have $\left|R_{1}(t)\right|\leq c_{1}$ and $\left|R_{2}(t)\right|\leq c_{2}$. For sufficiently small $\epsilon$ we have $t(u,x|s)\in[-\eta,\eta]$ for all $u$ and $x$, implying that as $\epsilon\rightarrow 0$ there exists a solution of (S10) that satisfies (S11), which in particular implies

[TABLE]

and by construction the conditions (S9) are satisfied for that solution. Thus, for sufficiently small $\epsilon$ the $t(u,x|s)$ defined above satisfies the conditions of Lemma S2. Applying that lemma we thus obtain that, for sufficiently small $\epsilon$, we have

[TABLE]

Again applying the expansion for $\rho(t)$ in (S6), and part (i) of Lemma S3 we thus obtain that

[TABLE]

This is what we wanted to show.

Proof for case (ii):

In this case, according to part (ii) of Lemma S3 the functions $R_{1}(t)$ and $R_{2}(t)$ are continuous and bounded over all $t\in\mathbb{R}$. In addition, we have assumed that $\mathbb{E}_{*}\left|\gamma(Y,X)-\delta(U,X)\right|^{3}<\infty$, and $\mathbb{E}_{*}\left|\psi(Y,X)\right|^{3}<\infty$, which guarantees that all of the expectations in (S10) are finite. We therefore again conclude that for small $\epsilon$ the equations (S10) have a solution such that (S12) holds. The remainder of the proof is equivalent to the proof of part (i), that is, we again apply Lemma S2 and Lemma S3 to obtain (S13).

S1.2 Proof of Theorem S1 (containing Theorem 1 as a special case)

Part (i):

We first want to show that $b_{\epsilon}(\gamma_{\beta,\sigma_{*}}^{\mathrm{P}})\leq b_{\epsilon}(\gamma)+{\cal O}(\epsilon)$. By applying Lemma S1 to both $\gamma_{\beta,\sigma_{*}}(y,x)$ and $\gamma_{\beta,\sigma_{*}}^{\mathrm{P}}(y,x)=\mathbb{E}_{p_{\beta,\sigma_{*}}}[\delta_{\beta}(U,X)\,|\,Y=y,X=x]$ we obtain, as $\epsilon\rightarrow 0$,

[TABLE]

where

[TABLE]

Here, to simplify $b_{\epsilon}(\gamma^{\mathrm{P}})$ we used that by the law of iterated expectations we have that $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}[\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)]-\mathbb{E}_{f_{\sigma_{*}}}[\delta_{\beta}(U,X)]=0$ (that is, the first term in $b_{\epsilon}(\gamma)$ is not present in $b_{\epsilon}(\gamma^{\mathrm{P}})$) and also $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\left(\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right)\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)\right]=0$ (that is, the vector $\lambda$ is equal to zero for $\gamma^{\mathrm{P}}$). We also use that under the reference model $\widetilde{\gamma}_{{\beta},{\sigma}_{*}}(Y,X)-\widetilde{\delta}_{\beta}(U,X)-\lambda^{\prime}\,\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)$ has zero mean, implying that its variance equals its second moment.

For any $\gamma_{{\beta},{\sigma}_{*}}(y,x)$ with $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}[\gamma_{{\beta},{\sigma}_{*}}(Y,X)]-\mathbb{E}_{f_{\sigma_{*}}}[\delta_{\beta}(U,X)]\neq 0$ we have $b_{\epsilon}(\gamma^{\mathrm{P}})\leq b_{\epsilon}(\gamma)$ for sufficiently small $\epsilon$, and the statement of the theorem thus holds in that case. In the following we therefore consider the case that $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}[\gamma_{{\beta},{\sigma}_{*}}(Y,X)]-\mathbb{E}_{f_{\sigma_{*}}}[\delta_{\beta}(U,X)]=0$. The expression for $b_{\epsilon}(\gamma)$ then simplifies to

[TABLE]

Again applying the law of iterated expectations we find that

[TABLE]

Using this we obtain

[TABLE]

Since $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left\{\left[\widetilde{\gamma}_{{\beta},{\sigma}_{*}}(Y,X)-\widetilde{\delta}_{\beta}(U,X)-\lambda^{\prime}\,\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)\right]-\left[\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right]\right\}^{2}\geq 0$ we thus conclude that

[TABLE]

and therefore we obtain that

[TABLE]

This is the first statement of the theorem. This concludes the proof of part (i) of Theorem S1, of which Theorem 1 in the main text is a special case.

Part (ii):

Next, let $\gamma_{\beta,\sigma_{*}}(y,x)$ be such that

[TABLE]

Then, the specification error expansions in (S14) are still valid, and using those we conclude that we must have

[TABLE]

because otherwise that term dominates all other terms in (S16). We also conclude that we must have

[TABLE]

for (S16) to hold. Furthermore, the calculation in (S15) is still valid here, and the inequality in the last display can therefore equivalently be rewritten as

[TABLE]

where we write $=$ instead of $\leq$, because the left hand side expression is non-negative. Applying Markov’s inequality we thus find that

[TABLE]

Defining

[TABLE]

we therefore obtain

[TABLE]

where in the last step we have used (S17) and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}[\psi_{\beta,\sigma_{*}}(Y,X)]=0$. Finally, notice that by construction we have

[TABLE]

S1.3 Proof of Theorem 2

We are going to show Theorem 2, which we restate here.

Theorem**.**

Assume that $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\psi_{{\beta},{\sigma}_{*}}(Y,X)=0$, $\phi(t)=\frac{1}{2}(t-1)^{2}$, and that $\gamma_{{\beta},{\sigma}_{*}}(Y,X)$ and $\delta_{{\beta}}(U,X)$ have finite second moments under the reference model. Then, for $0<\epsilon\leq\overline{\epsilon}$, we have

[TABLE]

with $\gamma_{\beta,\sigma_{*}}^{\mathrm{P}}(y,x)$ given by (S1) and $\overline{\epsilon}$ given by (18).

Proof of Theorem 2.

In the following proof, we again omit the arguments $\beta$, $\sigma_{*}$. By defining $Q(\gamma,f_{0}):=\mathbb{E}_{P(\beta,f_{0})}[\gamma_{{\beta},{\sigma_{*}}}(Y,X)]-\mathbb{E}_{f_{0}}[\delta_{\beta}(U,X)]$ we can rewrite (16) as

[TABLE]

Using that $\phi(t)=\frac{1}{2}(t-1)^{2}$, we find that the two worst-case distributions that maximize ($s=+$) and minimize ($s=-$) the function $Q(\gamma^{\mathrm{P}},f_{0})$ over $f_{0}\in\Gamma_{\epsilon}$ are given by

[TABLE]

Notice that our condition $\epsilon\leq\overline{\epsilon}$ guarantees that $f^{(s)}_{0}(u|x)\geq 0$ for all $u$, $x$, $s$.

We refer to Lemma S2 above for a more general derivation of the worst-case distribution (S18). The results of Lemma S2 simplifies here, because we assume that $\phi(t)=\frac{1}{2}(t-1)^{2}$ in the current theorem, which implies that the function $\rho(t)$ defined in (S4) is now given by $\rho(t)=\max(0,1+t)$, and we are also only interested in the worst-case distribution $f^{(s)}_{0}(u|x)$ at $\gamma=\gamma^{\mathrm{P}}$. Using (S18) we find that

[TABLE]

For a given $\gamma(y,x)$ we set

[TABLE]

Using that $\mathbb{E}_{*}[\gamma^{\mathrm{P}}(Y,X)-\delta(U,X)\,|\,Y=y,X=x]=0$, we then find that

[TABLE]

We now calculate

[TABLE]

where the first step is the definition of $b_{\epsilon}(\gamma)$, the second step is a property of the supremum, the third step uses (S18), the fourth step uses (S20), the fifth step uses (S19), and the sixth step uses that $\mathbb{E}_{*}[\gamma(Y,X)-\delta(U,X)]$ and $s_{\gamma}$ have the same sign. The result in the last display is exactly the statement of the theorem.

S1.4 Proof of Theorem 3

We are going to show Theorem 3, which we restate here.

Theorem**.**

Let $\gamma_{\beta,\sigma_{*}}^{\mathrm{P}}$ as in (S1). Then, for all $\epsilon>0$,

[TABLE]

The following lemma is useful for the proof of this theorem (Theorem 3 in the main text).

Lemma S4**.**

Let $\epsilon\geq 0$, $\beta\in{\cal B}$, $\sigma_{*}\in{\cal S}$, and let $\zeta:{\cal U}\times{\cal X}\rightarrow\mathbb{R}$. Then we have

[TABLE]

The proof of this lemma is given in Section S2. Notice that both Theorem 3 and Lemma S4 require that $\phi(r)$ is convex with $\phi(1)=0$, but they do not require $\phi^{\prime\prime}(1)>0$. For example, $\phi(r)=|r-1|/2$ is allowed here, which gives the total variation distance for $d(f_{0},f_{\sigma_{*}})$.

Proof of Theorem 3.

By definition we have

[TABLE]

By writing $\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)=\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)-\left[\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)\right]$ we obtain

[TABLE]

where in the second-to-last step we have used Lemma S4 with $\zeta(u,x)=\gamma_{{\beta},{\sigma}_{*}}(g_{\beta}(u,x),x)-\delta_{\beta}(u,x)$. We have thus shown that $b_{\epsilon}(\gamma^{\mathrm{P}})\leq 2\,b_{\epsilon}(\gamma)$ holds for any function $\gamma_{\beta,\sigma_{*}}(y,x)$, which implies that

[TABLE]

Appendix S2 Proofs of Technical Lemmas

Proof of Lemma S2.

In the following we assume that $f_{\sigma_{*}}(u|x)f_{X}(x)>0$ for all $(u,x)$ in the joint domain of $(U,X)$. This is without loss of generality, because we can define the joint domain of $(U,X)$ such that this is the case. With a slight abuse of notation we continue to write ${\cal U}\times{\cal X}$ for the joint domain, even though this need not be a product set.

To account for the absolute value in the definition of $b_{\epsilon}(\gamma)$ in (16) we let

[TABLE]

for $s\in\{-1,1\}$. We then have $b_{\epsilon}(\gamma)=\max_{s\in\{-1,1\}}b_{\epsilon}(\gamma,s)$. In the following we drop the arguments $\beta$ and $\sigma_{*}$ everywhere, that is, we simply write $g(u,x)$, $\gamma(y,x)$, $\delta(u,x)$, $f_{*}(u|x)$, $\psi(y,x)$, $\lambda^{(1)}(s,x)$, $\lambda^{(2)}(s)$, $\lambda^{(3)}(s)$ instead of $g_{\beta}(u,x)$, $\gamma_{\beta,\sigma_{*}}(y,x)$, $\delta_{\beta}(u,x)$, $f_{\sigma_{*}}(u|x)$, $\psi_{\beta,\sigma_{*}}(y,x)$, $\lambda^{(1)}_{\beta,\sigma_{*}}(s)$, $\lambda^{(2)}_{\beta,\sigma_{*}}(s)$, $\lambda^{(3)}_{\beta,\sigma_{*}}(s)$. The optimal $f_{0}(u|x)$ in the definition of $b_{\epsilon}(\gamma,s)$ solves, for $u,x\in{\cal U}\times{\cal X}$ almost surely under the reference distribution,

[TABLE]

where $\mu_{1}(s,x)\in\mathbb{R}$, $\mu_{2}(s)>0$, $\mu_{3}(s)\in\mathbb{R}^{\dim\psi}$ are Lagrange multipliers, which we choose to reparameterize as follows

[TABLE]

Those (reparameterized) Lagrange multipliers need to be chosen such that the constraints

[TABLE]

are satisfied. We need $\lambda^{(2)}(s)>0$ because the second constraint here is actually an inequality constraint ($\leq\epsilon$). Our assumptions guarantee that $f_{*}(u|x)>0$ and $f_{X}(x)>0$. We can therefore rewrite (S21) as follows,

[TABLE]

where $r=f_{0}\,f_{*}(u|x)$, the objective function was multiplied with $f_{\sigma_{*}}(u|x)f_{X}(x)$ (which does not change the value of the $\operatorname*{argmax}$), and $t(u,x|s)=t_{\beta,\sigma_{*}}(u,x|s)$ is defined in the statement of the lemma. Comparing the last display with the definition of $\rho(t)$ in (S4) we find that if $\rho\left[t(u,x|s)\right]<\infty$, then

[TABLE]

The condition $\rho\left[t(u,x|s)\right]<\infty$ is implicitly imposed in the statement of the lemma, because otherwise we could not have $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\,\rho\left[t_{\beta,\sigma_{*}}(U,X|s)\right]=1$. Using the result in the last display we find that the constraints (S22) are exactly the conditions (S5) imposed in the lemma. Under the conditions of the lemma we therefore have

[TABLE]

and from $b_{\epsilon}(\gamma)=\max_{s\in\{-1,1\}}b_{\epsilon}(\gamma,s)$ we thus obtain the statement of the lemma.

Proof of Lemma S3.

Part (i):

For $\nu=0$ we have $\phi=\overline{\phi}$. Our assumptions imply that there exists $\tau>0$ such that $\phi^{\prime}(r)$, $\phi^{\prime\prime}(r)$, $\phi^{\prime\prime\prime}(r)$ and $\phi^{\prime\prime\prime\prime}(r)$ are all uniformly bounded over $r\in[1-\tau,1+\tau]$. We can choose $\eta>0$ such that $[\rho(-\eta),\rho(\eta)]\subset[1-\tau,1+\tau]$. The conjugate of the convex function $\phi:\mathbb{R}\rightarrow\mathbb{R}$ is given by

[TABLE]

We have $\rho(t)=\phi_{*}^{\prime}(t)$, which is the inverse function of $\phi^{\prime}(r)$; that is, $\phi^{\prime}(\rho(t))=t$. We can express all derivatives of $\phi_{*}$ in terms of derivatives of $\phi$, for example, $\phi_{*}^{\prime\prime}(t)=1/\phi^{\prime\prime}(\rho(t))$ and $\phi_{*}^{\prime\prime\prime}(t)=-\phi^{\prime\prime\prime}(\rho(t))/[\phi^{\prime\prime}(\rho(t))]^{3}$. A Taylor expansion of $\rho(t)=\phi_{*}^{\prime}(t)$ around $t=0=\phi^{\prime}(1)$ reads

[TABLE]

where by the mean-value formula for the remainder term we have

[TABLE]

Similarly, a Taylor expansion of $\phi(\rho(t))=t\,\rho(t)-\phi_{*}(t)$ around $t=0$ reads

[TABLE]

where again by the mean-value formula for the remainder we have

[TABLE]

Continuity of $R_{1}(t)$ and $R_{2}(t)$ in a neighborhood of $t=0$ is also guaranteed by $\phi^{\prime}(r)$ being four times continuously differentiable in neighborhood around $r=1$. This concludes the proof of part (i).

Part (ii):

For $\nu>0$ the function $\phi(r)=\overline{\phi}(r)+\nu\,(r-1)^{2}$ still satisfies all the assumptions of part (i) of the lemma, that is, we can apply part (i) to find that there exists $\widetilde{c}_{1}>0$, $\widetilde{c}_{2}>0$ and $\eta>0$ such that for all $t\in[-\eta,\eta]$ we have

[TABLE]

What is left to show here is that there exists constant $c_{1}>0$ and $c_{2}>0$ such that (S7) also holds for $t<-\eta$ and for $t>\eta$.

We have $\phi^{\prime}(r)=\overline{\phi}^{\prime}(r)+\nu(r-1)$. Plugging in $r=\rho(t)$ we have $\phi^{\prime}(\rho(t))=t$, and therefore $t=\overline{\phi}^{\prime}(\rho(t))+\nu[\rho(t)-1]$. Our assumptions imply that $\overline{\phi}^{\prime}(\rho(t))>0$ for $t>0$ and $\overline{\phi}^{\prime}(\rho(t))<0$ for $t<0$. We therefore find that

[TABLE]

Using (S24) and (S25), and choosing $c_{1}=\max\left\{\widetilde{c}_{1},\,[1/\nu+1/\phi^{\prime\prime}(1)]/\eta\right\}$, we obtain

[TABLE]

for all $t\in\mathbb{R}$. This is the first inequality that we wanted to show.

Using again the convex conjugate defined in (S23) we have

[TABLE]

where in the second to last step we used that $r=1$ is one possible choice for $r\geq 0$, and we have $\phi(1)=0$, and in the last step we used that ${\mathrm{s}ign}[\rho(t)-1]={\mathrm{s}ign}(t)$. Our assumptions imply that $\phi(r)\geq 0$, that is, $\left|\phi(r)\right|=\phi(r)$. The result in the last display together with (S25) therefore give

[TABLE]

for all $t\in\mathbb{R}$. Using this and (S24), and choosing $c_{2}=\max\left\{\widetilde{c}_{2},\,[1/\nu+1/\{2\phi^{\prime\prime}(1)\}]/\eta\right\}$, we thus obtain

[TABLE]

for all $t\in\mathbb{R}$, which is the second inequality that we wanted to show. Continuity of $R_{1}(t)$ and $R_{2}(t)$ in $\mathbb{R}$ is also guaranteed by $\phi^{\prime}(r)$ being four times continuously differentiable in $r\in(0,\infty)$. This concludes the proof of part (ii).

Proof of Lemma S4.

Let $f_{0}\in\Gamma_{\epsilon}$. Remember the definition of the posterior density $p_{\beta,\sigma_{*}}(u\,|\,y,x)$ in (13). Define

[TABLE]

Then, for any $x\in{\cal X}$ we have $\widetilde{f}_{0}(u|x)\geq 0$, for all $u\in{\cal U}$, and $\int_{\cal U}\widetilde{f}_{0}(u|x)du=1$; that is, $\widetilde{f}_{0}(u|x)$ is a probability density over ${\cal U}$. Furthermore, by construction we have

[TABLE]

We also find that

[TABLE]

Furthermore, we have

[TABLE]

where we defined

[TABLE]

Using the definition of $p_{\beta,\sigma_{*}}(u\,|\,y,x)$ one can verify that $K_{\beta,\sigma_{*}}(\tilde{u}|u,x)\geq 0$, for all $\tilde{u}\in{\cal U}$, and $\int_{\cal U}K_{\beta,\sigma_{*}}(\tilde{u}|u,x)d\tilde{u}=\frac{\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[p_{\beta,\sigma_{*}}(u\,|\,Y,x)\right]}{f_{\sigma_{*}}(u\,|\,x)}=1$, almost surely (under $P(\beta,f_{\sigma_{*}})$) for $u\in{\cal U}$ and $x\in{\cal X}$. Thus, $K_{\beta,\sigma_{*}}(\tilde{u}|u,x)$ is a probability density over $\tilde{u}\in{\cal U}$, for all $u,x$. Also using that $\phi(r)$ is convex, we can therefore apply Jensen’s inequality to obtain

[TABLE]

Because $\widetilde{f}_{0}$ satisfies (S27) and (S28) we thus have $\widetilde{f}_{0}\in\Gamma_{\epsilon}$. We have thus shown that for every $f_{0}\in\Gamma_{\epsilon}$ there exists $\widetilde{f}_{0}\in\Gamma_{\epsilon}$ such that (S26) holds. Let $\widetilde{\Gamma}_{\epsilon}$ be the set of all such $\widetilde{f}_{0}$ obtained for an $f_{0}\in\Gamma_{\epsilon}$. Since $\widetilde{\Gamma}_{\epsilon}\subset\Gamma_{\epsilon}$ we find that

[TABLE]

Appendix S3 Robustness in prediction

Under squared loss, we wish to find a predictor $\gamma_{\widehat{\beta},\widehat{\sigma}}(Y_{i},X_{i})$, for some function $\gamma$, such that the worst-case mean squared prediction error is minimum. That is, our goal is to minimize

[TABLE]

with respect to $\gamma$. Similarly to our measure of worst-case specification error, here the mean squared prediction error is asymptotic, hence well-suited for settings with a large cross-section (e.g., settings with many teachers).

We first state the following local result, which is a direct generalization of Lemma 1.

Lemma S5**.**

In addition to defining \widetilde{\psi}(y,x)=\psi(y,x)-\mathbb{E}_{*}\left[\psi(Y,X)\big{|}X=x\right], let \widetilde{\gamma}(y,x)=\gamma(y,x)-\mathbb{E}_{*}\left[\gamma(Y,X)\big{|}X=x\right] and \widetilde{\delta}(u,x)=\delta(u,x)-\mathbb{E}_{*}\left[\delta(U,X)\big{|}X=x\right]. Suppose that $\phi(r)=\overline{\phi}(r)+\nu\,(r-1)^{2}$, with $\nu\geq 0$, and a function $\overline{\phi}(r)$ that is four times continuously differentiable with $\overline{\phi}(1)=0$ and $\overline{\phi}^{\prime\prime}(r)>0$, for all $r\in(0,\infty)$. Assume $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\psi_{{\beta},{\sigma}_{*}}(Y,X)=0$ and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left[\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)\widetilde{\psi}_{\beta,\sigma_{*}}(Y,X)^{\prime}\right]>0$. Furthermore, assume that one of the following holds:

• (i)

$\nu=0$, and the functions $\left|\gamma_{{\beta},{\sigma}_{*}}(y,x)\right|$, $\left|\delta_{\beta}(u,x)\right|$ and $\left|\psi_{{\beta},{\sigma}_{*}}(y,x)\right|$ are bounded over the domain of $Y$, $U$, $X$.

• (ii)

$\nu>0$, and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left|\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right|^{6}<\infty$, and $\mathbb{E}_{P(\beta,f_{\sigma_{*}})}\left|\psi_{{\beta},{\sigma}_{*}}(Y,X)\right|^{3}<\infty$.

Then, as $\epsilon\rightarrow 0$ we have

[TABLE]

where

[TABLE]

Let $\gamma^{\mathrm{P}}$ as in (17), so $\gamma^{\mathrm{P}}_{\widehat{\beta},\widehat{\sigma}}(Y_{i},X_{i})$ is the empirical Bayes estimate of $\delta_{\beta}(U_{i},X_{i})$. Under correct specification of the reference density $f_{\sigma}$, the posterior mean $\gamma^{\mathrm{P}}_{\beta,\sigma_{*}}(Y_{i},X_{i})$ is the minimum mean squared error predictor of $\delta_{\beta}(U_{i},X_{i})$ under squared loss. Under misspecification of $f_{\sigma}$, Lemma S5 implies that the leading term of the worst-case mean squared error is minimized at $\gamma=\gamma^{\mathrm{P}}$. Moreover, the lemma also implies the stronger result that the first-order term in the expansion of the worst-case mean squared prediction error (which is a multiple of $\epsilon^{\frac{1}{2}}$) is also minimized at $\gamma^{\mathrm{P}}$, provided the following condition holds almost surely:

[TABLE]

While (S29) is restrictive in general, it is satisfied in the fixed-effects model (1), when the researcher wishes to predict the quality $\alpha_{i}$ of teacher $i$. Indeed, in that case (S29) is equivalent to the posterior skewness of $\alpha_{i}$ being zero, when using the normal reference model as the prior. Since the normal distribution is symmetric, (S29) is satisfied, and the empirical Bayes estimator $\gamma^{\mathrm{P}}_{\widehat{\beta},\widehat{\sigma}}(Y_{i},X_{i})=\widehat{\mu}_{\alpha}+\widehat{\rho}(\overline{Y}_{i}-\widehat{\mu}_{\alpha})$ has minimum worst-case mean squared prediction error, up to second-order terms in $\epsilon^{\frac{1}{2}}$.

We also have a fixed-$\epsilon$ bound in the spirit of Theorem 3.

Theorem S2**.**

Let $\gamma_{\beta,\sigma_{*}}^{\mathrm{P}}$ as in (S1). Then, for all $\epsilon>0$,

[TABLE]

Theorem S2 shows that EB estimators are optimal, up to a factor of at most four, in terms of worst-case mean squared prediction error. In model (1), when ${\varepsilon}_{1},...,{\varepsilon}_{J}$ are normally distributed and $\alpha_{1},...,\alpha_{N}$ are parameters belonging to an $L^{2}$ ball, empirical Bayes James-Stein estimators are known to be optimal in terms of asymptotic minimax mean squared error since they achieve the Pinsker bound (see Wasserman, 2006, Chapter 7). Here, by contrast, we consider a worst case computed in a set of unrestricted, possibly non-normal joint distributions of $\alpha,\varepsilon_{1},...,\varepsilon_{J}$.

Proof of Lemma S5.

This statement of the lemma is obtained from Lemma S1 by replacing $(\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X))$ by $\left(\gamma_{{\beta},{\sigma}_{*}}(Y,X)-\delta_{\beta}(U,X)\right)^{2}$. The proof is obtained by the same replacement from the proof of Lemma S1.

Proof of Theorem S2.

By definition we have

[TABLE]

Using that $(a-b)^{2}\leq 2(a^{2}+b^{2})$ with $a=\delta_{\beta}(U,X)-\gamma_{{\beta},{\sigma}_{*}}(Y,X)$ and $b=\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(Y,X)-\gamma_{{\beta},{\sigma}_{*}}(Y,X)$ we obtain

[TABLE]

We furthermore have

[TABLE]

where in the first step we used the definition of $\gamma^{\mathrm{P}}_{{\beta},{\sigma}_{*}}(y,x)$, in the second step we applied the Cauchy-Schwarz inequality, and in the last line we used Lemma S4 and the definition of $e_{\epsilon}(\gamma)$. Combining the results of the last two displays we obtain that

[TABLE]

Appendix S4 Simulations

Here we provide details on the simulations summarized in Subsection 4.1. We consider four data generating processes (DGP) based on model (1). We use a standard normal as the (possibly misspecified) reference model for $\alpha_{i}$. In the first two DGP, we draw $\alpha_{i}$ from a standard normal distribution. Hence, the reference normal model for $\alpha_{i}$ is correctly specified in this case. We compare two specification for $\varepsilon_{ij}$. In DGP 1, $\varepsilon_{ij}$ are i.i.d. standard normal. In DGP 2, we model heteroskedastic errors as ${\cal{N}}(0,s(\alpha_{i})^{2})$, where $s(\alpha_{i})=\boldsymbol{1}\{\alpha_{i}>0\}\times.1+\boldsymbol{1}\{\alpha_{i}\leq 0\}\times 1.41$, so the variance of errors is the same as in DGP 1. In the next two DGP, we draw $\alpha_{i}$ from a Beta distribution with parameters $(11,1)$, shifted and rescaled such that $\alpha_{i}$ has mean 0 and variance 1. This distribution is skewed to the left, and the reference normal model for $\alpha_{i}$ is thus misspecified. In DGP 3, $\varepsilon_{ij}$ are i.i.d. standard normal. In DGP 4, we model heteroskedastic errors as ${\cal{N}}(0,s(\alpha_{i})^{2})$, where $s(\alpha_{i})=\boldsymbol{1}\{\alpha_{i}>0\}\times.1+\boldsymbol{1}\{\alpha_{i}\leq 0\}\times 1.61$, so the variance of errors is the same as in the other DGP.

In all DGP, we compare the performance of four estimators: the fixed-effects estimator given by (2), the PAE given by (5), the model-based estimator given by (6), and a nonparametric kernel deconvolution estimator with normal errors (Stefanski and Carroll, 1990). Unlike the other three estimators, the nonparametric kernel deconvolution estimator requires choosing a tuning parameter. We use the MISE-minimization approach of Delaigle and Gijbels (2004). To implement bandwidth selection and estimator, we use the codes available on Aurore Delaigle’s page: https://researchers.ms.unimelb.edu.au/$\sim$aurored/links.html#Code. To estimate the variance components $s_{\alpha}^{2}$ and $s_{\varepsilon}^{2}$, we use a minimum-distance estimator based on the empirical covariance matrix of the $Y_{ij}$. We set $\mu_{\alpha}=0$ and do not estimate it.

We run 1000 simulations for $n=1000$, and report estimates of $F_{\alpha}(a)=\mathbb{E}[\boldsymbol{1}\{\alpha\leq a\}]$ for a grid of $a$ values. We focus on bias (in absolute value), standard deviation, and root-MSE. We vary the number of measurements by taking $J\in\{2,20\}$. When $J=2$, the shrinkage factor is $\rho=.67$. In addition, the average value of the informativeness measure (21), over simulations and $a$ values, is $R^{2}=27\%$. When $J=20$, the shrinkage factor is $\rho=.95$. In addition, the average informativeness value is $R^{2}=69\%$.

In Figure S1, we show the results of the simulations in the cases where $\alpha_{i}$ is normally distributed. In this case, we expect the model-based estimator (in dashed lines) to be particularly well-behaved, since it is consistent and, for a suitable weighting of the minimum-distance estimator of $s_{\alpha}^{2}$ and $s_{\varepsilon}^{2}$, efficient as well. In DGP 1, $\varepsilon_{ij}$ are normal homoskedastic, so both the PAE (in solid) and the deconvolution estimator (in dotted) are consistent for fixed $J$ as $n$ tends to infinity. In the top two panels of Figure S1, we see that the bias of these two estimators is indeed small. Yet, the model-based estimator has the smallest bias, as well as the smallest variance and root-MSE.

In the bottom two panels of Figure S1, $\varepsilon_{ij}$ are heteroskedastic, so for fixed $J$ as $n$ tends to infinity neither the PAE nor the deconvolution and fixed-effects estimators are consistent, yet the model-based estimator remains consistent in this case. We see that the PAE and the model-based estimator behave comparably to the case of DGP 2, with rather small biases and root-MSE, with smallest bias and root-MSE achieved by the model-based estimator. However, the performance of the deconvolution and fixed-effects estimators worsens relative to the homoskedastic case, especially when $J=2$. This suggests that the PAE is less sensitive to this particular form of misspecification than the deconvolution and fixed-effects estimators.

We next turn to cases where the reference distribution of $\alpha_{i}$ is misspecified, which is the focus of our theory. In Figure S2, we show the results of the simulations when $\alpha_{i}$ is distributed as a shifted and rescaled Beta with parameters $(11,1)$. In this case, the model-based estimator is substantially biased, as shown by the left column in Figure S2, and the bias remains similar when varying $J$ and the distribution of $\varepsilon_{ij}$. In the top two panels of Figure S2, $\varepsilon_{ij}$ are normal homoskedastic, so the deconvolution estimator is consistent as $n$ tends to infinity for fixed $J$. When $J=2$, the deconvolution estimator is biased, which is likely to reflect the ill-posedness of the estimation problem. When $J=20$, the bias and root-MSE of the deconvolution estimator are small. Interestingly, although it has no consistency guarantees in this DGP, the PAE performs relatively well. Indeed, while the PAE is biased for $J=2$, its variance is smaller than the one of the deconvolution estimator. In addition, when $J=20$ and the posterior conditioning is more informative, the performance of the PAE and deconvolution estimator improves, the latter still dominating the former.

In the bottom two panels of Figure S2, $\varepsilon_{ij}$ are heteroskedastic. We see that this form of misspecification has large effects on the performance of the deconvolution estimator, especially when $J=2$ so the signal-to-noise ratio is lower. However, the performance of the PAE is very similar to the homoskedastic case: it is slightly less biased than the model-based estimator when $J=2$, and it is substantially less biased and has small root-MSE when $J=20$. In these designs the fixed-effects estimator and the PAE perform similarly.

Appendix S5 Extensions

In this section of the appendix we consider a number of issues in turn: how to compute PAE when they are not available in closed form, how to estimate quantities of interest that are nonlinear in $f_{0}$, whether the constant two appearing in Theorem 3 can be improved upon, how to construct confidence intervals, how to perform specification tests, how to derive the form of minimum-MSE estimators, and how to interpret PAE as Bayesian estimators in models where $U$ has finite support.

S5.1 Computation

$\widehat{\delta}^{\mathrm{P}}$ can be computed in closed form in simple models, such as all the examples in this paper. However, in complex models such as structural models, the likelihood function or posterior density may not be available in closed form. A simple approach in such cases is to proceed by simulation.

Specifically, for all $i=1,...,n$ we first draw $U_{i}^{(s)}$, $s=1,...,S$ according to $f_{\widehat{\sigma}}(\cdot\,|\,X_{i})$, and compute $Y_{i}^{(s)}=g_{\widehat{\beta}}(U_{i}^{(s)},X_{i})$. Then, we regress $\delta_{\widehat{\beta}}(U_{i}^{(s)},X_{i})$ on $Y_{i}^{(s)}$, for $s=1,...,S$. Any nonparametric/machine learning regression estimator can be used for this purpose. This procedure requires virtually no additional coding given simulation codes for outcomes and counterfactuals.

S5.2 Nonlinear effects

The researcher may be interested in a nonlinear function of $f_{0}$. Specifically, here we abstract from covariates $X$ and focus on $\overline{\delta}=\varphi_{\beta}(f_{0})$, for some functional $\varphi_{\beta}$. As an example, in the fixed-effects model (1), $\overline{\delta}$ may be the Gini coefficient of $\alpha$. The analysis in the linear case applies verbatim to this case, since under regularity conditions

[TABLE]

which is linear in $f_{0}$, up to smaller-order terms. Here $\nabla\varphi_{\beta}$ denotes the gradient of $\varphi_{\beta}(f)$ with respect to $f$. In Appendix S6 we report model-based and posterior estimates of Gini coefficients based on simulated data.

S5.3 The constant in Theorem 3

The binary choice model that we describe in Section S6 below is helpful to see that the global bound in Theorem 3, which depends on the constant two, cannot be improved upon in general. To see this, consider the binary choice model (S33) of Section S6 with three simplifications: $X$ consists of a single value, $b$ is known, and $\sigma_{*}=1$ is fixed. We assume that $x^{\prime}b>X^{\prime}b$.

In this example, for $\epsilon$ large enough the worst-case blue specification errors of $\widehat{\delta}^{\mathrm{M}}$ and $\widehat{\delta}^{\mathrm{P}}$ are

[TABLE]

and

[TABLE]

respectively.

From this, we first see that the specification error of the posterior estimator is smaller than twice that of the model-based estimator. In addition, taking $X^{\prime}b=0$ and $x^{\prime}b=\eta$, we have, for small $\eta$,

[TABLE]

This shows that two is indeed the smallest possible constant in Theorem 3.

S5.4 Confidence intervals

Consider first the correctly specified case. Suppose that $\widehat{\beta}$ and $\widehat{\sigma}$ are asymptotically linear in the sense that, for some mean-zero function $h$, we have

[TABLE]

Then, under standard conditions (e.g., Newey and McFadden, 1994), we have

[TABLE]

Here, $\Sigma_{11}=\mathop{\mathrm{V}ar}_{*}\left(G_{1}^{\prime}h(Y,X)+\mathbb{E}_{*}[\delta(U,X)\,|\,X]\right)$, \Sigma_{12}=\mathop{\mathrm{C}ov}_{*}\big{(}G_{1}^{\prime}h(Y,X)+\mathbb{E}_{*}[\delta(U,X)\,|\,X],\allowbreak G_{2}^{\prime}h(Y,X)+\mathbb{E}_{*}[\delta(U,X)\,|\,Y,X]\big{)}, $\Sigma_{21}=\Sigma_{12}$, and $\Sigma_{22}=\mathop{\mathrm{V}ar}_{*}\left(G_{2}^{\prime}h(Y,X)+\mathbb{E}_{*}[\delta(U,X)\,|\,Y,X]\right)$, for $G_{1}=\partial_{\beta,\sigma}\mathbb{E}_{\beta,\sigma_{*}}\,[\delta_{\beta}(U,X)]$ and $G_{2}=\mathbb{E}_{\beta,\sigma_{*}}\left\{\partial_{\beta,\sigma}\mathbb{E}_{p_{\beta,\sigma_{*}}}\,[\delta_{\beta}(U,X)\,|\,Y,X]\right\}$, where $\partial_{\theta}g(\theta_{1})$ denotes the gradient of $g(\theta)$ at $\theta=\theta_{1}$. Note that in (S31) we allow $\delta_{\beta}$ to be non-smooth in $\beta$ (e.g., an indicator function).

Consider next the locally misspecified case. A simple possibility to ensure uniform coverage within an $\epsilon$-neighborhood is to add $b_{\epsilon}(\gamma)$ on both sides of a standard confidence interval of $\overline{\delta}$. For example, one may construct the 95% interval

[TABLE]

for $\widehat{\Sigma}_{22}=\mathop{\mathrm{V}ar}_{*}\left(G_{2}^{\prime}h(Y,X)+\mathbb{E}_{*}[\delta(U,X)\,|\,Y,X]\right)$, where expectations and variances are taken with respect to $P(\widehat{\beta},f_{\widehat{\sigma}})$, and $\delta$, $G_{2}$, and $h$ are evaluated at $\widehat{\beta}$ and $\widehat{\sigma}$. Note that this confidence interval requires setting a value for $\epsilon$. Building on Hansen and Sargent (2008), Bonhomme and Weidner (2018) propose to interpret $\epsilon$ by relating it to the local power of a specification test.

S5.5 Specification test

Using the asymptotic distribution of $(\widehat{\delta}^{\mathrm{M}},\widehat{\delta}^{\mathrm{P}})$ under correct specification of $f_{\sigma}$, we obtain

[TABLE]

where $\widetilde{\Sigma}=\mathop{\mathrm{V}ar}_{*}\left(\mathbb{E}_{*}[\delta(U,X)\,|\,Y,X]-\mathbb{E}_{*}[\delta(U,X)\,|\,X]+(G_{2}-G_{1})^{\prime}h(Y,X)\right)$. Hence, under correct specification,

[TABLE]

Plugging-in a consistent empirical counterpart for $\widetilde{\Sigma}$ in this expression, we obtain a simple test of correct specification of the parametric density $f_{\sigma}$.

S5.6 Minimum local worst-case MSE estimator

Here we explain why $\widehat{\delta}^{\mathrm{M}MSE}$ in (22) gives the estimator with minimum worst-case MSE, in a local asymptotic framework where $n$ tends to infinity, $\epsilon$ tends to zero, and $n\epsilon$ tends to a positive constant. We only consider the case where $\beta$ and $\sigma_{*}$ are known and not estimated; that is, we have $\psi(y,x)=0$. Then, finding $\gamma^{\mathrm{M}MSE}(y,x)$ such that $\widehat{\delta}^{\mathrm{M}MSE}$ minimizes worst-case MSE over $f_{0}\in\Gamma_{\epsilon}$ can, to leading order in $\epsilon$ and $n^{-1}$, be shown to be equivalent to minimizing

[TABLE]

See Bonhomme and Weidner (2018) for details.

Next, applying Lemma 1 and noting that $\mathbb{E}_{*}[\gamma(Y,X)-\delta(U,X)]=0$ is required for MSE minimization (since adding a constant to $\gamma(y,x)$ such that $\mathbb{E}_{*}[\gamma(Y,X)-\delta(U,X)]=0$ has no effect on the higher order bias terms in Lemma 1, nor on ${\mathrm{V}ar}_{*}[\gamma(Y,X)]$, it is optimal to eliminate the leading bias term $\mathbb{E}_{*}[\gamma(Y,X)-\delta(U,X)]$ in this way), we find that to leading order in $\epsilon$ and $n^{-1}$ the worst-case MSE reads

[TABLE]

This expression for the approximate worst-case MSE depends on the distribution of $X$, which is unknown. For the minimum local worst-case specification error result in Theorem 1, it does not matter that the distribution of $X$ is unknown, because that distribution is identified from the sample as $n\rightarrow\infty$. However, for the MSE result here we have to take a stand on how to deal with the randomness in the observed covariates. In the following we condition on the observed sample of covariates, and replace all population expectations over $X$ by sample averages over $X_{i}$, $i=1,\ldots,n$. We write $\widehat{\mathbb{E}}_{X}$ for those sample averages. The worst-case MSE objective function in the last display then reads

[TABLE]

By the law of total variance we have

[TABLE]

In the following we can ignore the term \mathbb{E}_{*}\left\{{\mathrm{V}ar}_{*}\left[\delta(U,X)\,\big{|}\,Y,X\right]\big{|}\,X\right\}, because it does not depend on $\gamma(y,x)$. Then, the leading approximation to the worst-case MSE is given by the sample average over $X$ of

[TABLE]

Clearly, if for any given $X=x$ we find $\gamma(y,x)$ that minimizes this objective function, then its expected value over the sample distribution of $X$ is also minimized. The corresponding first-order condition for $\gamma^{\mathrm{M}MSE}(Y,X)$ reads

[TABLE]

The solution to this first-order condition is

[TABLE]

where we have now written $\widehat{\mathbb{E}}_{X}$ as $\frac{1}{n}\sum_{i=1}^{n}$.

The corresponding minimum local MSE estimator for $\overline{\delta}=\mathbb{E}_{*}\left[\delta_{\beta}(U,X)\right]$ is then given by

[TABLE]

which is the result stated in equation (22) of the main text.

S5.7 Finite support

Here we consider the case where $U$ has finite support and takes the values $u_{1},u_{2},...,u_{K}$ with probability $\omega_{1}^{0},...,\omega_{K}^{0}$. Here we abstract away from $\beta$, $\sigma$, and covariates $X$.

Injective and non-injective models.

Let $\delta_{k}=\delta(u_{k})$, and denote $g_{k}=g(u_{k})$ where $Y=g(U)$. Let $\overline{g}_{1},...,\overline{g}_{L}$ denote the $L\leq K$ equivalence classes of $g_{1},...,g_{K}$. We will denote as $\ell(k)\in\{1,...,L\}$ the index corresponding to the equivalence class of $g_{k}$, for all $k$. In addition, let $n_{\ell}=\sum_{i=1}^{n}\boldsymbol{1}\{Y_{i}=\overline{g}_{\ell}\}$ for all $\ell$, and denote $\omega_{k}^{U}=f(u_{k})$ for all $k$.

It is useful to distinguish two cases. When $g$ is injective, $K=L$ and $\mathbb{E}_{p(f)}[\delta(U)\,|\,g(U)=g_{k}]=\delta_{k}$. So we have $\widehat{\delta}^{\mathrm{P}}=\frac{1}{n}\sum_{k=1}^{K}n_{k}\delta_{k}$. This estimator does not depend on the assumed $f$. Moreover, as $\mathop{\mathrm{m}in}_{k=1,...,K}\,n_{k}$ tends to infinity we have

[TABLE]

Hence $\widehat{\delta}^{\mathrm{P}}$ is consistent for $\overline{\delta}$, irrespective of the choice of the reference density $f$, provided $\omega_{k}^{U}>0$ for all $k$.

When $g$ is not injective, $K\neq L$ and we have

[TABLE]

Moreover,

[TABLE]

Hence,

[TABLE]

Through $\overline{\delta}_{\ell}^{U}$, $\widehat{\delta}^{\mathrm{P}}$ depends on the prior $\omega^{U}$ in general, even as $\mathop{\mathrm{m}in}_{\ell=1,...,L}\,n_{\ell}$ tends to infinity.

Bayesian interpretation.

From a Bayesian perspective, one may view $\omega^{0}$ as a parameter, and put a prior on it. A simple conjugate prior specification is a Dirichlet distribution $\omega\sim\text{Dir}(K,\alpha)$, where $\alpha_{k}>0$ for $k=1,...,K$. We will focus on the posterior mean

[TABLE]

for a Dirichlet prior with $\alpha_{k}=M\omega_{k}^{U}$ for all $k$, where $M>0$ is a constant.

For all $\ell$, let $\overline{\alpha}_{\ell}=\sum_{k=1}^{K}\boldsymbol{1}\{\ell(k)=\ell\}\alpha_{k}$, and $\overline{\omega}_{\ell}=\sum_{k=1}^{K}\boldsymbol{1}\{\ell(k)=\ell\}\omega_{k}$. $(\overline{\omega}_{1},...,\overline{\omega}_{L})$ follows the Dirichlet distribution $\text{Dir}(L,\overline{\alpha})$. Moreover, for all $k$, $\omega_{k}/\overline{\omega}_{\ell(k)}$ is a component of a Dirichlet distribution with mean $\alpha_{k}/\overline{\alpha}_{\ell(k)}$.

Unlike the $\overline{\omega}_{\ell}$’s, the $\omega_{k}/\overline{\omega}_{\ell(k)}$’s are not updated in light of the data since they do not enter the likelihood. Notice the link with the Bayesian analysis of partially identified models in Moon and Schorfheide (2012): here the $\overline{\omega}_{\ell}$’s are identified but the $\omega_{k}$’s are not, since for identical $g_{k}$’s the data provides no information to discriminate across $\omega_{k}$’s.

As a result, we have

[TABLE]

It thus follows that

[TABLE]

Hence, under a diffuse Dirichlet prior centered around $\omega^{U}$, the Bayesian posterior mean coincides with the PAE.

Appendix S6 Posterior average effects in various settings

In this section, we provide additional examples of models where PAE may be of interest, and we show illustrative simulations for two models.

S6.1 Models

Linear regression.

Consider the linear regression

[TABLE]

Suppose that $\mathbb{E}[XU]=0$, and that the OLS estimator $\widehat{b}$ is consistent for $b$. Suppose also that the researcher is interested in the average effect $\overline{\delta}=\mathbb{E}_{f_{0}}[U^{2}XX^{\prime}]$. In this example $\overline{\delta}$ is multi-dimensional; see Appendix S5.

In this context, a model-based approach consists in modeling $U\,|\,X$, say, as a normal with zero mean and variance $s^{2}$, and computing

[TABLE]

where $\widehat{s}^{2}=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-X_{i}^{\prime}\widehat{b})^{2}$ is the maximum likelihood estimator of $s^{2}$ under normality.

By contrast, a PAE is

[TABLE]

This is the central piece in the White (1980) variance formula. $\widehat{\delta}^{\mathrm{P}}$ remains consistent for $\overline{\delta}$ absent normality or homoskedasticity of $U$. In this very special case, $\widehat{\delta}^{\mathrm{P}}$ is thus fully robust to misspecification of $f_{s}$, since $U_{i}$ is a deterministic function of $Y_{i}$, $X_{i}$ and $b$.

Censored regression.

Consider next the censored regression model

[TABLE]

In this model, $b$ can be consistently estimated under weak conditions. For example, Powell’s (1986) symmetrically trimmed least-squares estimator is consistent for $b$ when $U\,|\,X$ is symmetric around zero, under suitable regularity conditions. In this setting, suppose that we are interested in a moment of the potential outcomes $Y_{i}^{*}$, such as $\overline{\delta}=\mathbb{E}_{f_{0}}[h(Y^{*})]$ for some function $h$. As an example, the researcher may wish to estimate a feature of the distribution of wages using a sample affected by top- or bottom-coding.

Following a model-based approach, let us assume that $U\,|\,X\sim{\cal{N}}(0,s^{2})$, and estimate ${s}^{2}$ using maximum likelihood. A model-based estimator is then $\widehat{\delta}^{\mathrm{M}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}_{f_{\widehat{s}}}[h(X_{i}^{\prime}\widehat{b}+U)]$. By contrast, a PAE is

[TABLE]

This estimator relies on actual $Y$’s for uncensored observations, and on imputed $Y$’s for censored ones.

The censored regression model illustrates an aspect related to the class of neighborhoods that our theoretical characterizations rely on. In model (S32), the researcher might want to impose that $U\,|\,X$ be symmetric around zero, which is the main assumption for consistency of the Powell (1986) estimator. It is possible to construct estimators that minimize local worst-case specification error in an $\epsilon$-neighborhood that only consists of symmetric densities $f_{0}$. However, PAE may no longer have minimum specification error in this class. More generally, the assumptions that justify the use of a particular estimator $\widehat{\beta}$ may suggest further restrictions on the neighborhood. Our worst-case specification error results are based on a class where such restrictions are not imposed. Indeed, the only additional restriction on $f_{0}$, beyond belonging to an $\epsilon$-neighborhood around $f_{\sigma_{*}}$, is that the population moment condition $\mathbb{E}_{P(\beta,f_{0})}[\psi_{\beta,\sigma_{*}}(Y,X)]=0$ is assumed to hold, and we do not impose further restrictions that might be natural in order to justify the validity of this moment condition.

Binary choice.

Consider now the binary choice model

[TABLE]

In this model, Manski (1975, 1985) shows that $b$ is identified up to scale as soon as the median of $U\,|\,X$ is zero, under sufficiently large support of $X$. In addition, he provides conditions for consistency of the maximum score estimator $\widehat{b}$, again up to scale. Manski’s conditions, however, are not sufficient to consistently estimate the average structural function (ASF, Blundell and Powell, 2004)

[TABLE]

Let us take as reference parametric distribution for $U\,|\,X$ a normal with zero mean and variance ${s}^{2}$, and let $\widehat{s}^{2}$ denote the maximum likelihood estimator of $s^{2}$ given $\widehat{b}$, based on normality. Specifically, $\widehat{s}$ maximizes the probit log-likelihood $\sum_{i=1}^{n}Y_{i}\log\Phi\left(\frac{X_{i}^{\prime}\widehat{b}}{s}\right)+(1-Y_{i})\log\Phi\left(-\frac{X_{i}^{\prime}\widehat{b}}{s}\right)$. A model-based estimator of the ASF is $\widehat{\delta}^{\mathrm{M}}(x)=\Phi\left(\frac{x^{\prime}\widehat{b}}{\widehat{s}}\right)$, and a posterior estimator is

[TABLE]

Unlike $\widehat{\delta}^{\mathrm{M}}(x)$, the posterior ASF estimator $\widehat{\delta}^{\mathrm{P}}(x)$ depends directly on the observations of the binary $Y_{i}$’s, in addition to the indirect data dependence through $\widehat{b}$ and $\widehat{s}^{2}$. In the next subsection we present simulations from an ordered choice model, which suggest that the informativeness of the posterior conditioning — and the robustness properties of posterior estimators compared to model-based estimators — depend crucially on the support of the dependent variable.

Panel data discrete choice.

Our last example is the panel data model

[TABLE]

When $\varepsilon_{it}$ are i.i.d. standard logistic, $b$ can be consistently estimated using the conditional logit estimator (Andersen, 1970, Chamberlain, 1984). However, additional assumptions are needed to consistently estimate average partial effects such as the effect of a discrete shift of $\Delta$ along the $k$-th component of $X$,

[TABLE]

where $e_{k}$ is a vector of zeros with a one in the $k$-th position.

The standard approach is to postulate a parametric random-effects specification for the conditional distribution of $\alpha$ given $X_{1},...,X_{T}$, and to compute an average effect $\widehat{\delta}^{\mathrm{M}}$ with respect to that distribution. By contrast, a posterior estimator is computed conditional on the observations $Y_{i1},...,Y_{iT}$, for every individual $i$. As $T$ tends to infinity, such estimators are robust to misspecification of $\alpha$, provided $\varepsilon_{t}$ is correctly specified (Arellano and Bonhomme, 2009). Our analysis shows that they also have robustness properties when $T$ is fixed and $n$ tends to infinity.

Aguirregabiria et al. (2018) show that conditional logit-like estimators can also be used to consistently estimate parameters in structural dynamic discrete choice settings. As an example, they study the Rust (1987) model of bus engine replacement in the presence of unobserved heterogeneity in maintenance and replacement costs. In such structural models, estimating average welfare effects of policies requires averaging with respect to the distribution of unobservables. PAE provide an alternative to the standard parametric model-based approach in this context.

S6.2 Simulations

Here we report the results of two simulation exercises, based on the fixed-effects model (1), and on an ordered choice model.

S6.2.1 Fixed-effects model

Skewness.

Let us consider the fixed-effects model (1). Suppose the parameter of interest is the skewness of $\alpha$

[TABLE]

For example, it is of interest to estimate the skewnesses of income components and how they evolve over time (Guvenen et al., 2014). Since the normal distribution is symmetric, the model-based normal estimator of skewness is simply $\widehat{\delta}^{\mathrm{M}}=0$, irrespective of the observations $Y_{ij}$. Hence, $\widehat{\delta}^{\mathrm{M}}$ is not informed by the data, even when the empirical distribution of the fixed-effects $\overline{Y}_{i}=\frac{1}{J}\sum_{j=1}^{J}Y_{ij}$ indicates strong asymmetry.

By contrast, a PAE based on a normal reference distribution is

[TABLE]

It can be verified that

[TABLE]

where $\widehat{\rho}=\frac{\widehat{s}_{\alpha}^{2}}{\widehat{s}_{\alpha}^{2}+\widehat{s}_{\varepsilon}^{2}/J}$. Under mild conditions, and in contrast with $\widehat{\delta}^{\mathrm{M}}$, the posterior estimator $\widehat{\delta}^{\mathrm{P}}$ is consistent for the true skewness of $\alpha$ as $J$ tends to infinity. However, $\widehat{\delta}^{\mathrm{P}}$ is biased for small $J$ in general.

To provide intuition about the magnitude of the bias, we simulate data where all latent components are independent, $\varepsilon_{j}$ are standard normal, and $\alpha$ follows a skew-normal distribution (e.g., Azzalini, 2013) with zero mean, variance 1, and skewness $\approx.47$ corresponding to the skew-normal parameter $\delta=.99$. We take $n=1000$, and run $100$ simulations varying $J$ from 1 to 30. We estimate means and variances using minimum-distance based on first and second moment restrictions.

In the left panel of Figure S3 we show the results. We see that the model-based estimator is equal to zero irrespective of the number $J$ of individual measurements. By contrast, the posterior estimator converges to the true skewness of $\alpha$ as $J$ increases, although it is biased for small $J$.

Gini coefficient.

We next focus on the Gini coefficient of $\alpha$:

[TABLE]

In this case, a model-based estimator is

[TABLE]

while a PAE is, following (S30),

[TABLE]

where

[TABLE]

In the right panel of Figure S3 we show the simulation results. We see that in this case also the model-based estimator is insensitive to $J$. The posterior estimator has a lower bias, especially for larger $J$.

S6.2.2 Ordered choice model

We next consider the ordered choice model

[TABLE]

for a sequence of known thresholds $-\infty=\mu_{0}<\mu_{1}<...<\mu_{J-1}<\mu_{J}=+\infty$. This model may be of interest to analyze data on wealth or income, say, where only a bracket containing the true observation is recorded. We focus on the average structural function

[TABLE]

We take as reference distribution $U\,|\,X\sim{\cal{N}}(0,s^{2})$. In the simulated data generating process, $U$ is independent of $X$, distributed as a re-centered $\chi^{2}$ with mean zero and variance one. We simulate a scalar standard normal $X$. We set $n=1000$, $b_{1}=.5$, $b_{0}=0$, $s=1$, and $\mu$ as uniformly distributed between $-2$ and $2$. We estimate $b$ up to scale using maximum score (Manski, 1985). Specifically, using maximum score we regress $\boldsymbol{1}\{Y_{i}\leq j\}$ on $X_{i}$ and a constant, for all $j$, imposing that the coefficient of $X_{i}$ is one. We then regress the $J$ estimates on a common constant and the $\mu_{j}$, and obtain the implied estimate for $b$ by rescaling. For computation of maximum score, we use the mixed integer linear programming algorithm of Florios and Skouras (2008).

In Figure S4 we report the results for $J=3$ (left) and $J=10$ (right). We see that, when $J=3$, model-based and posterior estimators are similarly biased. By contrast, when $J=10$, the posterior estimator aligns well with the true average structural function, even though the model-based estimator is substantially biased.

Appendix S7 Additional empirical results