text:
\section{introduction}\label{sec:intro}
Sample selection is ubiquitous in empirical economics. For example, it arises naturally in the estimation of wage equations because we do not observe wages of individuals who do not work \cite{gronau74,heckman74}, and product demands because we do not observed quantities purchased by consumers who do not have access to the product.  Sample selection biases the estimation of causal or predictive effects when the reasons for not observing the data are related to the outcome of interest. For example, there is sample selection bias in the estimation of a wage equation whenever the employment status  and offered wage depend on common unobserved variables such as ability, motivation or skills. The most popular solution to the sample selection bias is the Heckman selection model (HSM) introduced in \citeasnoun{heckman74}. This classical model offers a convenient and parsimonious way to account  for sample selection by making parametric assumptions about the outcome and selection processes. Our development is motivated by the observation that, in addition to the parametric structure, the classical model imposes strong homogeneity assumptions on how covariates affect the outcome and selection processes and how the selection process operates itself. We develop a generalization of the HSM that relaxes all these three homogeneity restrictions. The resulting model is a semiparametric model, where key parameters are function-valued, thereby considerably generalizing the classical selection model.

Following the literature, we model sample selection using two latent variables for the selection and
outcome processes and relate the distribution of these variables with the distribution of the corresponding observed variables. Here we find convenient to work with a local Gaussian representation (LGR) of the joint distribution of the latent variables, which we introduce in the paper.
This representation is unique for any joint distribution  and might be of independent interest in other settings. The identification analysis is very transparent with the LGR. Thus, we show that the parameters of the LGR are partially identified in the presence of endogenous sample selection because there are only two free probabilities to identify three parameters. We rely on exclusion restrictions to point-identify the three parameters nonparametrically. These conditions require of a binary covariate that  does not affect the distribution of the latent outcome and dependence between the latent selection and outcome variables.

Once we have established nonparametric identification with the exclusion restrictions, we introduce a  flexible semiparametric distribution regression (DR) model with covariates for the LGR. This model generalizes the HSM by adding multiple sources of heterogeneity to the selection and outcome processes. Thus, it allows for observed and unobserved heterogeneity in selection sorting, together with unobserved heterogeneity in the effect of the covariates on the selection sorting and outcome.  In the case of the wage equation, the model can capture the presence of heterogenous returns to schooling across the wage distribution, or positive sorting at the top of the wage distribution and negative sorting at the bottom.  The model is semiparametric because its parameters are function-valued and can be applied without modification to continuous, discrete and continuous-discrete outcomes. We show how to construct interesting functionals of the model parameters such as actual and counterfactual distributions of latent and observed outcomes, which can be applied to policy evaluation, treatment effects, wage decompositions and discrimination analysis accounting for sample selection. In the case of wage decompositions, we show how to identify two new effects: a selection sorting effect and a selection structure effect. Selection sorting is determined by whether the employed individuals have higher or lower offered or latent wages than unemployed individuals with the same characteristics. Selection structure is determined by the proportion of employed individuals and how they are selected based on observed characteristics.

We develop a two-step estimator for the model parameters. The first step consists of a probit regression for the selection equation, which is identical to the first step in the Heckman two-step method \cite{heckman79}. The second step estimates multiple DRs with sample selection correction. The difference between these DRs and the standard DRs without sample selection is that we run  bivariate probits instead of univariate probits \cite{foresiperacchi95,chernozhukov+13inference}. We estimate functionals of the parameters using the plug-in method. We derive functional central limit theorems for all the estimators and show how to use these results  to perform uniform inference on function-valued parameters. This type of inference is useful to construct confidence bands and test hypotheses such as whether a coefficient or effect is uniformly zero, constant or positive. We implement the inference methods using Kolmogorov-Smirnov type statistics where the critical values are obtained via multiplier bootstrap \cite{gz84} applied to the estimated influence functions (as in \citeasnoun{ch06}, \citeasnoun{ks12} and \citeasnoun{cck13}). This bootstrap scheme is convenient in our setting because it avoids repeated computation of estimators in constructing the bootstrap draws of the statistic. We prove the validity of multiplier bootstrap by deriving bootstrap functional central limit theorems for all the estimators.

We apply our methods to study the relationship between wage and employment in the U.K. using updated data from 1978 to 2013. To this end we estimate wage equations for men and women  and carry out several wage decompositions accounting for endogenous selection into employment.  Here, we uncover  positive sorting among single men and negative sorting among married women. This difference in selection sorting  is consistent with assortative matching in the marriage market. It also explains a substantive proportion (but not all) of the gender wage gap at the top of the distribution, which is consistent with recent explanations based on glass ceiling theory.  We still find the strong evidence of the gender wage gap -- most of the gender wage gap in offered wages (as well as in observed wages) is accounted by differences in the wage structure that are often associated with gender discrimination in the labor market. The  effect of education is positive and increases along the distribution. All the heterogeneity that we find is inconsistent with the restrictions of the classical selection model, lending support for the use of the generalized HSM.

\paragraph{ \textbf{Literature review}}   The sample selection problem has a long history in statistics and econometrics. Classical references can be found in \citeasnoun{gronau74}, \citeasnoun{heckman74}, \citeasnoun{lee82}, \citeasnoun{goldberger83}, \citeasnoun[Section 10.7]{amemiya85}, \citeasnoun[Section 9.4]{maddala86}, \citeasnoun{manski89}, \citeasnoun{manski94},   and  \citeasnoun{vella98}. A popular solution to the problem is the HSM developed by Heckman in a sequence of papers \cite{heckman74,heckman76,heckman79,heckman90}.  This model has been extended in several dimensions.  \citeasnoun{lee83}, \citeasnoun{prieger02} and \citeasnoun{smith03} replaced the bivariate standard normal copula with other parametric copulas, and \citeasnoun{marchenko12} replaced the bivariate normal  by a bivariate t-distribution  to apply the HSM to heavy tailed data. \citeasnoun{ahn93}, \citeasnoun{powell94}, \citeasnoun{andrews98}, and \citeasnoun{newey99}  developed semiparametric versions of the HSM and \citeasnoun{das03} a nonparametric version, all focusing on location effect versions with homogeneous effects. None of the models considered in these extensions accommodates all the sources of heterogeneity allowed by our model.

\citeasnoun{ab17} proposed another extension of the HSM, which like our model allows for multiple sources of heterogeneity.\footnote{See \citeasnoun{arellano17} for a recent survey on sample selection in quantile models.}  Their method relies on  quantile regression to model the marginal distribution of the latent outcome coupled with a parametric model for the copula of the latent selection and outcome variables.  They estimate the model parameters using a three-step method where the first step is the same as in our method, but the second and third steps involve an iterative procedure that alternates between quantile regressions to estimate the outcome equation and nonlinear GMM to estimate the parameters of the copula. They also rely on numerical simulation to estimate functionals of the parameters such as actual and counterfactual distributions of the latent and observed outcomes. Compared to our method, they model the covariates effects as direct on the conditional quantile of the latent distribution, whereas  we model the covariate effects as direct
on the conditional latent distribution -- hence in our framework covariates affect the conditional quantiles indirectly. Further, their modeling approach imposes homogeneity on the copula function, which rules
out forms of copula heterogeneity across the distribution of the latent outcome, which are permitted in our approach. Moreover, their quantile regression model requires the latent outcome to be continuous, whereas  our distribution regression model can deal with any type of outcome and is therefore more widely applicable. Our method is computationally simpler as it does not involve any iteration between methods in the second step, nor does rely on numerical simulation to estimate functionals. The identification assumptions are also different and not nested: we impose more structure on the dependence between the outcome and selection processes, whereas they require more variation on the excluded covariates. We provide a more detailed comparison of the identifying assumptions in Appendix \ref{app:ab17}.\footnote{Finally, from a technical point of view, \citeasnoun{ab17} only derived pointwise limit theory for the estimators of the model parameters, whereas we derive functional limit theory for the estimators of the model parameters and related functionals.}

\paragraph{\textbf{Outline}} Section \ref{sec:id}  examines the identification problem under sample selection using a new representation of a joint distribution. Section \ref{sec:model} introduces the DR model with selection and associated functionals,  estimators of the model parameters and functionals, and a multiplier bootstrap method to perform functional inference. Section \ref{sec:theory} provides asymptotic theory for the estimation and inference methods, and Section \ref{sec:empirics}
 reports the results of the empirical application. Appendices \ref{app:id}--\ref{app:scores} give the proofs of the main results and additional technical results. The online Supplemental Material (SM) contains additional empirical results.

\section{Another View of the Sample Selection Problem}\label{sec:id}
\subsection{Local Gaussian Representation of a Joint Distribution}
We start by characterizing a local Gaussian representation (LGR) of the joint distribution of two random variables that is convenient to provide a new view of the  identification problem with sample selection and motivate our modeling choices later.

Let $Y^*$ and $D^*$ be two random variables with joint cumulative distribution function (CDF) $F_{Y^*,D^*}$ and marginal CDFs $F_{Y^*}$ and $F_{D^*}$. We label these variables with  asterisks because they will be latent variables when we introduce sample selection.  Our first result shows that $F_{Y^*,D^*}$ can be represented via a standard bivariate normal distribution at a point and with a correlation parameter that depend on the evaluation point $(y,d)$.

\begin{lemma}[LGR Result]\label{lemma:lgr} Let $F_{Y^*,D^*}$ be a joint CDF, then, for any $(y,d) \in \mathbb{R}^2$,
$$
F_{Y^*,D^*}(y,d) = \Phi_2(\mu(y), \nu(d); \rho(y,d)),
$$
where $\mu(y) \in \mathbb{R}$, $\nu(d) \in \mathbb{R}$, $\rho(y,d) \in [-1,1]$, and $\Phi_2(\cdot, \cdot; \rho)$ is the joint CDF of a standard bivariate normal random variable with parameter $\rho$. Moreover, the values of $\mu(y)$, $\nu(d)$ and  $\rho(y,d)$ are uniquely determined by   $\mu(y) = \Phi^{-1}(F_{Y^*}(y))$, $\nu(d) = \Phi^{-1}(F_{D^*}(d))$, and the solution in $\rho$ of
$$
F_{Y^*,D^*}(y,d) = \Phi_2(\Phi^{-1}(F_{Y^*}(y)), \Phi^{-1}(F_{D^*}(d)); \rho),
$$
where $\Phi$ is the standard normal CDF.  Hence, the representation is unique.
\end{lemma}

Lemma \ref{lemma:lgr} establishes that any joint CDF admits a unique representation as a sequence of standard bivariate normal distributions.  This result is stronger than the comprehensive property of the Gaussian copula that establishes that this copula includes the two Frechet bounds and independent copula by suitable choice of the correlation parameter,  e.g., \citeasnoun{smith03}.   Lemma \ref{lemma:lgr}  easily extends to CDFs conditional on covariates $Z$ by making all the parameters dependent on the value of $Z$.

The parameter $\rho(y,d)$ can be interpreted as a measure of local dependence.\footnote{See \citeasnoun{survey18} for a recent survey on measures of local dependence.} Thus, when $\rho(y,d) = 0$, the distribution $F_{Y^*,D^*}$ factorizes at $(y,d)$:
$$
F_{Y^*,D^*}(y,d) = \Phi_2(\Phi^{-1}(F_{Y^*}(y)), \Phi^{-1}(F_{D^*}(d)); 0) = F_{Y^*}(y) F_{D^*}(d),
$$
that is, the events $\{Y^*\leq y\}$ and $\{D^*\leq d\}$ are independent. Hence we can say that $Y^*$ and $D^*$ are ``locally independent" at $(y,d)$.\footnote{This concept is different from the local independence of \citeasnoun{doksum94}, which is local in only one of the variables. Thus, for example, our concept is symmetric in $Y^*$ and $D^*$ whereas the concept in \citeasnoun{doksum94} is not.}  In general, the discrepancy
$$
|\Phi_2(\Phi^{-1}(F_{Y^*}(y)), \Phi^{-1}(F_{D^*}(d)); \rho(y,d))- \Phi_2(\Phi^{-1}(F_{Y^*}(y)), \Phi^{-1}(F_{D^*}(d)); 0)|
$$
measures deviation away from independent factorization, thereby giving meaning to $\rho(y,d)$ as local dependence parameter.

In the LGR, the marginal CDFs of $Y^*$ and $D^*$ are represented by local Gaussian links
$$
F_{Y^*}(y) = \Phi(\mu(y)), \quad F_{D^*}(d) = \Phi(\nu(d)),
$$
and the copula of $Y^*$ and $D^*$ is represented by a local Gaussian copula
\begin{equation}\label{LGC}
\begin{array}{l}
 C_{Y^*,D^*}(u,v) =   \Phi_2(\Phi^{-1}(u), \Phi^{-1}(v); \rho(y_u,d_v)), \\
 \forall (u,v) \in [0,1] : \exists y_u \in \mathbb{R}: F_{Y^*}(y_u) = u, \quad \exists d_v  \in \mathbb{R}: F_{D^*}(d_v) = v.
\end{array}
\end{equation}
The LGR is convenient because it separates $\mu(y)$ and $\nu(d)$ as two parameters determining the marginals of $Y^*$ and $D^*$ from  $\rho(y,d)$ as a parameter determining the dependence between $Y^*$ and $D^*$.\footnote{Note that the marginals of $F_{Y^*,D^*}$ do not identify separately the mean and variances of the local Gaussian representation.}

\citeasnoun{kolev06} developed a closely related result to (\ref{LGC}) for the copula. They established that the copula of any bivariate distribution can be represented by the bivariate Gaussian copula with a local correlation parameter. The LGR is different from the local Gaussian approximation of \citeasnoun{tjostheim13}, which approximates the distribution of a continuous bivariate random variable in a neighborhood of the point of interest by a bivariate normal distribution with local vector of means and variance-covariance matrix, see also \citeasnoun{hjort96}. As \citeasnoun{tjostheim13} noticed, unlike the LGR, a local Gaussian approximation that intersects with the joint distribution at the point of interest is not unique.

\subsection{Identification of Sample Selection Model}
We consider now the sample selection problem where we observe two random variables $D$ and $Y$, which can be defined in terms of the latent variables  $D^*$ and $Y^*$ as
   \begin{eqnarray*}
    D &=& 1(D^* \leq 0),\\
    Y &=& Y^* \  \ \text{if} \ \ D=1,
    \end{eqnarray*}
i.e., $D$ is an indicator for $D^* \leq 0$ and $Y^*$ is only observed when $D = 1$.  The goal is to identify features of the joint distribution of the latent variables from the joint distribution of the observed variables.

The joint CDF of the observed variables can be written in terms of the LGR of $F_{Y^*,D^*}$ as
$$
F_{Y,D}(y,d) = \Phi_2(\mu(y), \nu(0); \rho(y,0))1(d \geq 1) + [1 - \Phi(\nu(0))]1(d \geq 0).
$$
As shown below, the parameters of the LGR are partially identified because $F_{Y^*,D^*}$ is only partially identified. We proceed by characterizing the identified set for these parameters  and provide exclusion restrictions to achieve point identification. Since there is a one-to-one relationship between $F_{Y^*,D^*}$ and its LGR, the identified set for the parameters of the LGR determine the identified set for $F_{Y^*,D^*}$.  In what follows,  we simplify the notation to
$$\nu := \nu(0), \quad \rho(y) := \rho(y,0).$$
We can only hope to identify $\nu(d)$  and $\rho(y,d)$ at $d=0$ since we only observe whether $D^* \leq 0$.

To understand the source of the partial identification in terms of the LGR, note that in the presence of sample selection there are two free probabilities,
$
\Pr(D = 1)$ and $\Pr(Y \leq y \mid D = 1)$,
to identify three parameters,  $\mu(y)$, $\nu$ and $\rho(y)$. The parameter $\nu$ is pinned down by the selection probability as
   $$
   \nu = \Phi^{-1}(\Pr(D = 1)).
   $$
The parameters $\mu(y)$ and $\rho(y)$ are partially identified as the solutions in $(\mu,\rho)$ to
   $$
   F_{Y,D}(y,1) - \Pr(D = 0) = \Phi_2(\mu,\Phi^{-1}(\Pr(D = 1)); \rho).
   $$
These solutions  form a one-dimensional manifold in $\mathbb{R} \times [-1,1]$ because $
\partial \Phi_2(\mu, \cdot; \rho)/\partial \mu  > 0
$, $
\partial \Phi_2(\cdot, \cdot; \rho)/\partial \rho > 0
$, and  $
\partial^2 \Phi_2(\cdot, \cdot; \rho)/\partial \mu \partial \rho > 0
$ \cite{spivak65,munkres91}.  The identified set of $(\mu(y),\rho(y))$ can be shrunk using additional information such as that $\rho(y)$ is known to lie in a subinterval of $[-1,1]$, e.g. $|\rho(y)| < 0.2$.

We use exclusion restrictions to achieve point identification of the parameters of the LGR.
To state these restrictions in terms of the LGR, we start by extending the LGR to CDFs conditional on covariates. Let $Z$ be a random variable and $F_{Y^*,D^* \mid Z}$ be the joint CDF of $Y^*$ and $D^*$ conditional on $Z$. Then,  $F_{Y^*,D^* \mid Z}$ admits the LGR:
  \begin{equation*}
  F_{Y^*,D^* \mid Z}(y, d \mid z) = \Phi_2(\mu(y \mid z), \nu(d \mid z); \rho(y,d \mid z)),
  \end{equation*}
 where $\mu(y \mid z) \in \mathbb{R}$, $\nu(d \mid z) \in \mathbb{R}$, and $\rho(y,d \mid z) \in [-1,1]$.
This representation can be characterized using the same argument as in Lemma \ref{lemma:lgr} after fixing the value of the covariate $Z$ and letting all the parameters of the LGR to depend on this value. The exclusion restrictions are:

\begin{assumption}[Exclusion Restrictions]\label{ass:er} There is a binary random variable $Z$  that satisfies:
   \begin{enumerate}
    \item Non-Degeneracy: $0 < \Pr(D = 1) < 1$ and $0 < \Pr(Z=1 \mid D = 1) <1$.
   \item Relevance: $
    \Pr(D = 1 \mid Z = 0) < \Pr(D = 1 \mid Z = 1) < 1.
   $
   \item Outcome exclusion: $\mu(y \mid z) = \mu(y)$.
   \item Selection exclusion: $\rho(y,0 \mid z) = \rho(y,0)$.
   \end{enumerate}
\end{assumption}
The condition that $Z$ is binary is made to emphasize that our identification strategy does not rely on large variation of $Z$. If $Z$ is not binary we only require that Assumption \ref{ass:er} be satisfied for two values of $Z$. Part (1) requires that there is sample selection and that $Z$ has variation in the selected population. It is used to guarantee that all the probabilities employed in the identification analysis are well-defined. Part (2) requires that $Z$ affects the probability of selection and rules out corner cases. The condition $\Pr(D = 1 \mid Z = 1) < 1$ precludes identification at infinity (see Remark \ref{remark:idi}). The sign of the first inequality can be reversed by relabelling the values of $Z$.   Part (3) is a standard exclusion restriction, which is not sufficient for point identification in the presence of sample selection \cite{manski94,manski03}. It holds when $Y^*$ is independent of $Z$.\footnote{\citeasnoun{kitagawa2010} developed a test for the outcome exclusion.} Part (4) is an exclusion restriction in the local dependence between $Y^*$ and $D^*$, which holds if selection sorting is independent of $Z$.
We explain this condition more in detail below with an example and compare it with the identification condition in \citeasnoun{ab17} in Remark \ref{remark:ab17}.

We can get some intuition about the outcome and selection exclusion restrictions with an example. Consider a heteroskedastic bivariate normal model for the latent variables, i.e.,
  $$
  (Y^*,D^*) \mid Z = z \sim \mathcal{N}_2\left( \left[\begin{array}{c} \mu_{Y^*}(z) \\ \mu_{D^*}(z)\end{array}\right], \left[\begin{array}{cc} \sigma_{Y^*}(z)^2 & \sigma_{Y^*,D^*}(z) \\ \sigma_{Y^*,D^*}(z) & \sigma_{D^*}(z)^2 \end{array}\right] \right).
  $$
The outcome exclusion imposes that
  $$\frac{y-\mu_{Y^*}(z)}{\sigma_{Y^*}(z)} = \mu(y),$$
  whereas the selection exclusion imposes that
  $$\frac{\sigma_{Y^*,D^*}(z)}{\sigma_{Y^*}(z)\sigma_{D^*}(z)} = \rho.$$
  If $Y^*$ is independent of $Z$, the outcome exclusion holds and the selection exclusion boils down to the condition that  $\sigma_{Y^*,D^*}(z)/\sigma_{D^*}(z)$ does not depend on $z$.  In other words, the slopes of the linear regressions of $Y^*$ on $D^*$ conditional on $Z$ are the same when $Z=0$ and $Z=1$.

We now show how the presence of an exclusion restriction helps identify the parameters of the conditional LGR. Under the exclusion restriction the conditional LGR at $d=0$ simplifies to
  \begin{equation}\label{eq:clgr}
  F_{Y^*,D^* \mid Z}(y, 0 \mid z) = \Phi_2(\mu(y), \nu(z); \rho(y)), \ \  z \in \{0,1\},
  \end{equation}
  where $\nu(z) := \nu(0 \mid z)$ and $\rho(y) := \rho(y,0)$.
The CDF of the observed variables conditional on $Z$ can be related to this conditional LGR  as
$$
   F_{Y,D \mid Z}(y,d \mid z) = \Phi_2(\mu(y), \nu(z); \rho(y)) 1(d \geq 1) +  [1 - \Phi(\nu(z))] 1(d \geq0) , \ \  z \in \{0,1\}.
$$
 As before, $\nu(z)$ is identified from the conditional selection probability as
  \begin{equation}\label{eq:nu}
  \nu(z) = \Phi^{-1}\left( \Pr(D=1 \mid Z = z)\right), \quad z \in \{0,1\}.
  \end{equation}
Moreover, $\mu(y)$ and $\rho(y)$ are identified as the solution in $(\mu, \rho)$  to
  \begin{equation}\label{eq:murho}
  F_{Y,D \mid Z}(y,1 \mid z) - \Pr(D=0 \mid Z=z) = \Phi_2(\mu, \Phi^{-1}\left( \Pr(D=1 \mid Z = z)\right); \rho), \ \ z \in \{0,1\}.
 \end{equation}
This is a nonlinear system of two equations in two unknowns that has unique solution because the Jacobian is a P-matrix for all $\mu \in \mathbb{R}$ and $\rho \in (-1,1)$ by Theorem 4 of \citeasnoun{gale65}.

The following theorem provides a detailed identification analysis of the parameters of  the conditional LGR in \eqref{eq:clgr}. It includes boundary cases where the parameters  $\mu(y)$ and $\rho(y)$ can be either point or partially identified.
\begin{theorem}[Identification of LGR with Sample Selection]\label{theorem:id} Assume that  Assumption \ref{ass:er} holds. Then, $\nu(z)$ is identified by \eqref{eq:nu} and there are the following cases for the identification of $\mu(y)$ and $\rho(y)$:
\begin{enumerate}
  \item  If $ F_{Y,D \mid Z} (y,1 \mid 1) - \Pr(D=0 \mid Z=1) = F_{Y,D \mid Z} (y,1 \mid 0) - \Pr(D=0 \mid Z=0)> 0$,
  $$
  \rho(y) = 1, \ \ \mu(y) = \Phi^{-1}\left(F_{Y,D \mid Z} (y,1 \mid 1) - \Pr(D=0 \mid Z=1) \right).
   $$
  \item  If $ F_{Y,D \mid Z} (y,1 \mid 1)  <  1$  and $F_{Y,D \mid Z} (y,1 \mid 0) = 1$,
  $$
   \rho(y) = 1, \ \ \mu(y) = \Phi^{-1}\left(F_{Y,D \mid Z} (y,1 \mid 1) - \Pr(D=0 \mid Z=1) \right).

  $$
  \item If $ F_{Y,D \mid Z} (y,1 \mid 1)  =  F_{Y,D \mid Z} (y,1 \mid 0) = 1$,
  $$
  \rho(y) = 1, \ \   \mu(y) \in [\Phi^{-1}(\Pr(D=1 \mid Z=1)), + \infty).
  $$
  \item  If $F_{Y,D \mid Z} (y,1 \mid 1) > \Pr(D=0 \mid Z=1)$ and $F_{Y,D \mid Z} (y,1 \mid 0) = \Pr(D=0 \mid Z=0)$,
$$
 \rho(y) = - 1, \  \mu(y) =  \Phi^{-1}\left(F_{Y,D \mid Z} (y,1 \mid 1) \right).
  $$
  \item  If  $F_{Y,D \mid Z} (y,1 \mid 1) =  F_{Y,D \mid Z} (y,1 \mid 0) < 1$,
\begin{eqnarray*}
& \rho(y) = - 1, \  \mu(y) = \Phi^{-1}\left(F_{Y,D \mid Z} (y,1 \mid 1) \right).
\end{eqnarray*}
  \item If $F_{Y,D \mid Z} (y,1 \mid z) = \Pr(D=0 \mid Z=z)$, $z \in \{0,1\}$,
  $$
\rho(y) = - 1, \   \mu(y) \in (-\infty, \Phi^{-1}\left(F_{Y,D \mid Z} (y,1 \mid 1)\right)].
  $$
\item Otherwise, $\mu(y)$ and $\rho(y)$ are point identified as the solution in $(\mu, \rho)$  to \eqref{eq:murho}. This solution exists and is unique.
\end{enumerate}
\end{theorem}

The boundary cases in Theorem \ref{theorem:id} are easy to detect. In practice,  partial identification usually occurs at extreme values of $y$. For example, case (3) arises for values of $y$ such that $Y > y$ a.s., and case (6) for values of $y$ such that $Y < y$ a.s.

\begin{remark}[Identification at Infinity]\label{remark:idi} When $\Pr(D = 1 \mid Z = 1) = 1$ and $\mu(y \mid z) = \mu(y)$, the conditional LGR at $z=1$ gives $F_{Y,D \mid Z}(y,1 \mid 1) = \lim_{\nu \nearrow + \infty} \Phi_2(\mu(y), \nu; \rho(y \mid 1)) = \Phi(\mu(y))$, which  identifies $\mu(y)$ by
   $$
  \mu(y) = \Phi^{-1}(F_{Y,D \mid Z}(y,1 \mid 1)),
  $$
 without the selection exclusion restriction.
This result is analogous to the identification at infinity of  \citeasnoun{chamberlain86} where $Z$ is continuous with unbounded support and
  $$
  \lim_{z \nearrow +\infty} \Pr(D = 1 \mid Z = z) = 1.
  $$
  Note that $\rho(y \mid z)$ is not point identified without further restrictions.\qed
\end{remark}

\begin{remark}[Comparison with \citeasnoun{ab17}, AB17]\label{remark:ab17}  Assumption \ref{ass:er} is not nested with the conditions that  AB17 used to show nonparametric identification of their model. We impose stronger restrictions in the dependence of the latent selection and outcome variables, but require less variation in the excluded covariate $Z$. We provide a more detailed comparison in Appendix \ref{app:ab17} \qed
\end{remark}

\section{Distribution Regression Model with Sample Selection}\label{sec:model}
\subsection{The Model}
We consider a semiparametric version of the LGR with covariates:
\begin{equation}\label{eq:drm}
    F_{Y^*,D^*}(y,0 \mid Z=z)   = \Phi_2(-x^\prime\beta(y), -z^\prime \pi;\rho(x'\delta(y))),
\end{equation}
where $Y^*$ is the latent outcome of interest, which can be continuous, discrete or mixed continuous-discrete;  $D^*$ is a latent variable that determines sample selection;   $X$ is a vector of covariates;  $Z=(Z_1,X)$; and $Z_1$ are excluded covariates, i.e., observed covariates that satisfy the exclusion restrictions. The excluded covariates avoid reliance on functional form assumptions to achieve identification. The model for the LGR consists of three indexes. We shall refer to $-x^\prime\beta(y)$ as the outcome equation, to $-z^\prime \pi$ as the selection equation, and to $\rho(x'\delta(y))$ as the selection sorting equation.  We observe the selection indicator $D = 1(D^* > 0)$ and the outcome $Y=Y^*$ when $D=1$.\footnote{The minus signs  in \eqref{eq:drm} are included to take into account that the selection is defined by $D^* > 0$ instead of $D^* \leq 0$. We use this definition to facilitate the interpretation of the parameters and the comparison with the classical Heckman selection model; see Example \ref{ex:hsm}.}   In the empirical application that we consider below, $Y^*$ is offered wage, $D^*$ is the difference between offered wage and reservation wage, $D$ is an employment indicator, $Y$ is the observed wage, $X$ includes labor market characteristics such as education, age, number of children and marital status, and $Z_1$ includes measures of out-of-work income. We shall discuss the validity of these measures as excluded covariates in Section \ref{sec:empirics}.

The model \eqref{eq:drm} is semiparametric because $y\mapsto \beta(y)$ and $y\mapsto \delta(y)$ are unknown functions, i.e.  infinite dimensional parameters in general. This flexibility allows the effect of $X$ on the outcome and selection sorting to vary across the distribution. For example, it allows the return to education to vary across the distribution, the selection sorting to be different for high and low educated individuals, or to have positive selection sorting  at the upper tail and negative at the bottom tail or vice versa. The function $u \mapsto \rho(u)$ is a known link  with range $[-1,1]$, e.g. the Fisher transformation \cite{fisher15}, $\rho(u) = \tanh(u)$. The corresponding distribution of $Y^*$ conditional on $Z$ is
   $$F_{Y^*}(y\mid Z=z)= \lim_{\nu \nearrow + \infty} F_{Y^*,D^*}(y,v \mid Z=z)=\Phi(-x^\prime \beta(y)), \ \ z = (x,z_1).$$
The selection bias arises because this distribution is different from the distribution of the observed outcome $Y$, i.e.
    $$F_{Y^*}(y\mid Z=z)\neq F_{Y}(y\mid Z=z,D=1) = \frac{\Phi_2(-x^\prime\beta(y), z^\prime \pi;-\rho(x'\delta(y)))}{\Phi(z^\prime \pi)}.$$

\begin{example}[HSM]\label{ex:hsm} Consider the \citeasnoun{heckman74} sample selection model (HSM):
\begin{eqnarray*}
D^* &=& Z^\prime\pi + V, \\
Y^* &=& X^\prime \beta+\sigma U,
\end{eqnarray*}
where $(U,V)$ is  independent of $Z$ and has standard bivariate normal distribution with parameter $\rho$, such that
    $$
        F_{Y^*,D^*}(y,0 \mid Z=z)
         = \Phi_2\left(\frac{y-x^\prime \beta}{\sigma}, -z^\prime \pi; \rho\right).
$$
This is a special case of model \eqref{eq:drm} with
$$\beta_1(y)=(\beta_1-y)/\sigma, \ \ \beta_{-1}(y)= \beta_{-1}/\sigma, \ \ \rho(x'\delta(y)) = \rho.$$
The HSM therefore imposes strong homogeneity restrictions in the selection process and effect of the covariates on the outcome and selection sorting.  Thus, only the intercept of $\beta(y)$ varies with $y$, and $\rho(x'\delta(y))$ is invariant to both $x$ and $y$.\qed
\end{example}

The model \eqref{eq:drm} has multiple data generating process representations as nonseparable systems. One  example is
   \begin{eqnarray*}
        D^* &=& Z^\prime\pi + V, \ \ V\mid Z \sim \mathcal{N}(0,1), \\
        0 &=& X^\prime \beta(Y^*) + \rho(X'\delta(Y^*))V+\sqrt{1-\rho(X'\delta(Y^*))^2} U, \ \  U\mid Z \sim \mathcal{N}(0,1),
        \end{eqnarray*}
where $U$ and $V$ are independent. For example, in the wage application $V$ can be interpreted as unobserved net benefit of working, and $U$ as unobserved skills or innate ability net of $V$. This representation is similar to the HSM in Example \ref{ex:hsm} with the difference that the equation for $Y^*$ is nonseparable.\footnote{Note that in Example \ref{ex:hsm} the equation for $Y^*$ can be written as $0= (X'\beta - Y^*)/\sigma + \rho V+\sqrt{1-\rho^2} \widetilde{U}$, where $\widetilde{U}$ is standard normally distributed and independent of $V$ and $Z$.}

\subsection{Functionals}
There are several functionals of the parameters of the model \eqref{eq:drm} that can be of interest. One is the marginal distribution of the latent outcome $Y^*$
   $$F_{Y^*}(y) = \int F_{Y^*}(y\mid Z=z)dF_Z(z)=\int \Phi(-x^\prime \beta(y))dF_X(x),$$
where $F_Z$ and $F_X$ are the marginal distributions of $Z$ and $X$, respectively.  In the case of the wage application, $F_{Y^*}$ corresponds to the distribution of the offered wage, which is a potential or latent outcome free of selection.   We can also construct counterfactual distributions by combining coefficients $\beta(y)$ and distributions $F_X$ from different populations or groups. These distributions are useful to decompose the distribution of offered wages between females and males or between blacks and whites, which can be used to uncover discrimination in the labor market.  Another functional is the probability of selection
$$
\Pr(Y = 1) = \int \Pr(Y=1 \mid Z=z) dF_Z(z) = \int \Phi(z^\prime \pi)dF_Z(z),
$$
which can be use to decompose differences in employment rates between employment structure effects, $\pi$, and composition effects, $F_Z$.

We can also use the model to construct distributions for the observed outcome using that
    \begin{align*}
    F_Y(y)

    &=
    \int \frac{\Phi_2\left( -x^\prime\beta(y), z^\prime \pi; -\rho(x'\delta(y)) \right)}{\Phi(z^\prime \pi)}dF_Z(z\mid D=1) \\
    &= \frac{\int \Phi_2\left( -x^\prime\beta(y), z^\prime \pi; -\rho(x'\delta(y)) \right) dF_Z(z)}{\int \Phi(z^\prime \pi) dF_Z(z)},
    \end{align*}
where the second equality follows from the Bayes rule. We can again construct counterfactual distributions by changing $\beta(y)$, $\pi$, $\delta(y)$ and $F_Z$. In the wage application, we will decompose the differences in the wage distribution between genders or across time into changes in the worker composition  $F_Z$, wage structure $\beta(y)$, selection structure $\pi$, and selection sorting $\delta(y)$. Both selection effects are new to this model.

\begin{remark}[Selection effects]\label{remark:se} To interpret the selection effects, it is useful to consider a simplified version of the model without covariates where
$    F_Y(y;\pi,\rho)
= \Phi_2\left( -\beta, \pi; -\rho \right)/\Phi (\pi).
$
Here we drop the dependence of $\beta$ and $\rho$ on $y$ to lighten the notation, and make explicit the dependence of $F_Y$ on the selection parameters $\pi$ and $\rho$ to carry out comparative statics with respect to them. Then, by the properties of the normal distribution
$$
    \frac{\partial  F_Y(y;\pi,\rho)}{\partial \rho} = - \frac{\phi_2(-\beta,\pi; -\rho)}{ \Phi(\pi)} < 0,
$$
and
$$
    \frac{\partial F_Y(y;\pi,\rho)}{\partial \pi} \propto \Phi\left(\frac{-\beta + \rho \pi}{\sqrt{1-\rho^2}}  \right) \Phi(\pi) - \int_{-\infty}^{\pi} \Phi\left(\frac{-\beta + \rho x}{\sqrt{1-\rho^2}}  \right) \phi(x) dx \left\{\begin{array}{c} < 0 \text{ if } \rho < 0, \\  = 0 \text{ if } \rho = 0, \\ > 0 \text{ if } \rho > 0, \end{array}\right.
$$
where $\Phi$ and $\phi$ are the standard normal CDF and probability density function (PDF), and  $\phi_2(\cdot, \cdot; \rho)$ be the joint PDF of a standard bivariate normal random variable with parameter $\rho$.\footnote{To obtain the  derivative we use that $\Phi_2\left( -\beta, \pi; -\rho \right) = \int_{-\infty}^{\pi} \Phi\left( \frac{-\beta + \rho x}{\sqrt{1-\rho^2}}\right) \phi(x) dx$.}
Increasing $\rho$ therefore shifts the distribution to the right (increases quantiles) because it makes selection sorting more positive while the size of the selected population is fixed. The effect of increasing $\pi$ is more nuanced and depends on the sign of $\rho$. Intuitively, $\pi$ affects the size of the selected population and the relative importance of observables and unobservables in the selection. For example, when selection sorting is negative, increasing the size of the selected population by increasing $\pi$ shifts the distribution of the right (increases quantiles) because the newly selected individuals have smaller (more negative) selection unobservables that correspond to larger (more positive) outcome unobservables. In other words, the newly selected individuals are relatively less adversely selected.

The sign of the selection effects might be different in the presence of covariates if the variation in the parameters changes the composition of the selected population. Consider the following extreme example with only one covariate based on the wage application. Let the covariate be an indicator for high skills. Assume that high-skilled workers are relatively more likely to participate than low-skilled workers, there is no selection sorting on unobservables, which corresponds to $\rho(x'\delta(y)) =0$ in the model, and the distribution of offered wages for high-skilled workers first-order stochastically dominates the same distribution for low-skilled workers. In this case increasing the probability of participation for high-skilled workers, which corresponds to increasing the component of $\pi$ associated to the high skill indicator in the model, both increases the overall probability of participation and shifts the distribution of observed wages to the right (increases quantiles), despite the lack of selection sorting. Intuitively, the distribution of observed wages is a mixture of the distribution of wages for employed high-skilled and low-skilled workers, and we are increasing the relative proportion of employed high-skilled workers. The opposite holds if the distribution of offered wages for high-skilled workers is first-order stochastically dominated by the same distribution for low-skilled workers. \qed
\end{remark}

Quantiles and other functionals of the distributions of latent and observed outcomes can be constructed by applying the appropriate operator. For example, the $\tau$-quantile of the latent outcome is  $Q_{Y^*}(\tau) = \mathbf{Q}_{\tau}(F_{Y^*})$, where $\mathbf{Q}_{\tau}(F) := \inf\{y \in \mathbb{R} : F(y) \geq \tau\}$ is the quantile or left-inverse operator.

\subsection{Estimation} To estimate the model parameters and functionals of interest, we assume that we have a random sample of size $n$ from $(D,DY,Z)$, $\{(D_i, D_i Y_i,Z_i)\}_{i=1}^n$,  where we use $D Y$ to indicate that we only observe $Y$ when $D=1$.

Before describing the estimators, it is convenient to introduce some notation.   Let  $\mathcal{Y}$ be the region of interest of $Y$, and denote $\theta_y := (\beta(y), \delta(y))$, where we replace the arguments in $y$ by subscripts to lighten the notation.\footnote{If the support of $Y$ is finite,  $\mathcal{Y}$ can be the entire support, otherwise $\mathcal{Y}$ should be a subset of the support excluding low density areas such as the tails.}

The estimation relies on the relationship between conditional distributions and binary regressions. Thus, the CDF of $Y$ at a point $y$ conditional on $X$ is the expectation that an indicator that $Y$ is less than $y$ conditional on $X$,
$$
F_{Y \mid X}(y \mid x) = {\mathrm{E}}[1(Y \leq y) \mid X = x].
$$
 To implement this idea,  we construct the set of indicators for the selected observations
$$
I_{yi} = 1(Y_i\leq y) \text{ if } D_i = 1,
$$
for each $y \in \mathcal{Y}$. In the presence of sample selection, we cannot just run a probit binary regression of $I_{yi}$ on $X_i$ to estimate the parameter $\beta(y)$ as in \citeasnoun{foresiperacchi95} and \citeasnoun{chernozhukov+13inference}. The problem is similar to running least squares in the HSM. Instead, we use that
 \begin{align*}
      \ell_i(\pi,\theta_y)=\left[1-\Phi(Z_i^\prime\pi)\right]^{1-D_i}&\times \Phi_2(-X_i^\prime\beta(y),Z_i^\prime\pi;-\rho(X_i'\delta(y)))^{D_iI_{yi}} \\ &\times\Phi_2(X_i^\prime\beta(y),Z_i^\prime\pi;\rho(X_i'\delta(y)))^{D_i(1-I_{yi})}
      \end{align*}
is the likelihood of $(D_i, I_{yi})$ conditional on $Z_i$. This likelihood is the same as the likelihood of a bivariate probit model or more precisely a probit model with sample selection \cite{zellner65,poirier80,vandeven81}.

We estimate the model parameters using a computationally attractive two-step method to maximize the average log-likelihood, similar to the Heckman two-step method. The first step is a probit regression for the probability of selection to estimate $\pi$, which is identical to the first step in the Heckman two-step method. The second step consists of multiple distribution regressions (DRs) with sample selection corrections to estimate $\beta(y)$ and $\delta(y)$ for each value of $y \in \mathcal{Y}$. These steps are summarized in the following algorithm:

\begin{algorithm}[Two-Step DR Method]\label{alg:tsdr}   (1) Run a probit for the selection equation to estimate $\pi$:
$$\hat{\pi} = \argmax_{c \in \mathbb{R}^{d_{\pi}} }L_1(c) = \frac{1}{n} \sum_{i=1}^n \left[D_{i}\log{\Phi(Z_i^\prime c)}+(1-D_{i})\log{\Phi(-Z_i^\prime c)}\right], \ \ d_{\pi} := \dim \pi.$$
(2) Run multiple DRs with sample selection correction to estimate $\theta_y$: for each $y \in \mathcal{Y}$
\begin{align*} \hat{\theta}_y =  \argmax_{t=(b,d) \in \Theta} L_2(t, \hat \pi) = \frac{1}{n} \sum_{i=1}^n D_{i}  \left[ \right.&I_{yi} \log{\Phi_2\left(-X_i^\prime b,Z_i^\prime \hat{\pi}; -\rho(X_i'd)\right)} \\
 + & \left.(1-I_{yi}) \log{\Phi_2\left(X_i^\prime b,Z_i^\prime \hat{\pi}; \rho(X_i'd) \right)} \right],
\end{align*}
where $\Theta \in \mathbb{R}^{d_{\theta}}$ is a compact parameter set, and
  $$
  d_{\theta} := \dim \theta_u, \ \ \rho(u) = \tanh(u) = \frac{e^u - e^{-u}}{e^u + e^{-u}} \in [-1,1], \quad \frac{\partial \rho(u)}{\partial u} > 0.
  $$
\end{algorithm}

In practice we replace the set  $\mathcal{Y}$ by a finite grid $\bar{\mathcal{Y}}$ if $\mathcal{Y}$ contains many values.

The estimators of the functionals of interested are constructed from the estimators of the parameters using the plug-in method. For example, the estimator of the distribution of the latent outcome is
\begin{equation}\label{eq:ldisthat}
\hat{F}_{Y^*}(y)=\frac{1}{n}\sum_{i=1}^n\Phi(-X_i^\prime \hat{\beta}(y)),
\end{equation}
the estimator of the probability of selection is
$$
\hat{\Pr}(Y = 1) =  \frac{1}{n}\sum_{i=1}^n \Phi(z^\prime \hat \pi),
$$
and the estimators of the counterfactual distributions of the observed outcome are constructed from
\begin{equation}\label{eq:odisthat}
\hat{F}_Y(y \mid D=1)=  \frac{\sum_{i=1}^n \Phi_2(-X_i^\prime \hat{\beta}(y),Z_i^\prime \hat{\pi};-\rho(X_i'\hat{\delta}(y)))}{\sum_{i=1}^n \Phi(Z_i^\prime \hat{\pi})},
\end{equation}
by choosing the estimators of $ \hat{\beta}(y)$, $\hat{\pi}$, and $\hat{\delta}(y)$ and the sample values of $Z$ appropriately.  Estimators of quantiles and other functionals of these distributions are obtained by applying the operators that define the functionals to the estimator of the distribution. For example, the estimator of the $\tau$-quantile of the latent outcome is $\hat Q_{Y^*}(\tau) = \mathbf{Q}_{\tau} (\hat F_{Y^*})$.\footnote{If the function $y \mapsto \hat F_{Y^*}(y)$ is not increasing, we monotonize it  before taking the left-inverse. See Section \ref{sec:theory} for details.}

\subsection{Uniform Inference}
The model parameters and functionals of interest are generally function-valued. We show how to construct confidence bands for them that can be used to test functional hypotheses such as the entire function be zero, non-negative or constant. To explain the construction consider the case where the functional of interest is a linear combination of the model parameter $\theta_y$, that is the function $y\mapsto c'\theta_y,$ $y\in\mathcal{Y}$, where $c \in \mathbb{R}^{d_{\theta}}$.  The set $CB_p(c'\theta_y)$ is an asymptotic $p$-confidence band for $c'\theta_y$ if it satisfies
  $$
  \Pr\left[c'\theta_y \in CB_p(c'\theta_y), \forall y \in\mathcal{Y}\right] \to p, \text{ as } n \to \infty.
  $$
We form  $CB_p(c'\theta_y)$  as
  $$CB_p(c'\theta_y)=c'\hat{\theta}_y\pm cv(p)SE(c'\hat{\theta}_y),$$
  where $\hat{\theta}_y$ is the estimator of $\theta_y$ defined in Algorithm \ref{alg:tsdr},  $SE(c'\hat{\theta}_y)$ is the standard error of $c'\hat{\theta}_y$, and $cv(p)$ is a critical value, i.e. a consistent estimator of the $p$-quantile of the statistic
  $$t_{\mathcal{Y}} = \sup_{y\in \mathcal{Y}} \frac{|c'\hat{\theta}_y-c'\theta_y|}{SE(c'\hat{\theta}_y)}.$$
We obtain the standard error and critical value from the limit distribution of $\hat{\theta}_y$  derived in Section \ref{sec:theory}.

In practice, it is convenient to estimate the critical value using resampling methods. Multiplier bootstrap is computationally attractive in our setting because it does not require parameter re-estimation and therefore avoids the nonlinear optimization in both steps of Algorithm \ref{alg:tsdr}. The multiplier bootstrap is implemented using the following algorithm:
\begin{algorithm}[Multiplier Bootstrap]\label{alg:mb} (i) For $b \in 1,\ldots, B$ and the finite grid $\bar{\mathcal{Y}} \subseteq \mathcal{Y}$, repeat the steps: (1) Draw the bootstrap multipliers $\{\omega_{i}^b : 1 \leq i \leq n \}$ independently from the data and normalized them to have zero mean,
   $$
    \omega_i^b = \tilde \omega_i^b - \sum_{i =1}^n \tilde  \omega_i^b/n, \ \ \tilde \omega_i^b \sim \text{ i.i.d. } \mathcal{N}(0,1).
    $$
(2) Obtain the bootstrap estimator of the model parameter
    $$
    \hat \theta^b_y = \hat \theta_y + n^{-1} \sum_{i = 1}^n \omega_i^b \ \widehat \psi_i(\hat \theta_y, \hat \pi),
    $$
    where $\widehat \psi_i(\hat \theta_y, \hat \pi)$ is an estimators of the influence function of $\hat{\theta}_y$ given in \eqref{eq:if}.
(2) Construct bootstrap realization of maximal t-statistic $t_{\mathcal{Y}}$ for the functional of interest,
    $$
    t_{\mathcal{Y}}^b = \max_{y \in \bar{\mathcal{Y}}} \frac{|c'\hat \theta^b_y - c'\widehat \theta_y|}{SE(c'\hat{\theta}_y)}, \ \ SE(c'\hat{\theta}_y) = \sqrt{c' \widehat \Sigma_{\theta_y\theta_y} c},
    $$
    where $\widehat \Sigma_{\theta_y\theta_y}$ is an estimator of the asymptotic variance-covariance matrix of $\hat{\theta}_y$ given in \eqref{eq:se}.
(ii) Compute the critical value $cv(p)$ as the simulation $p$-quantile of $t_{\bar{\mathcal{Y}}}^b$,
$$
cv(p) =  p-\text{quantile of }  \{t_{\mathcal{Y}}^b : 1 \leq b \leq B\}
$$
\end{algorithm}

 The centering of the multipliers in step (i1) of the algorithm is a finite sample adjustment. Confidence bands for other functionals of the model parameter can be constructed using a similar bootstrap method.

\section{Asymptotic Theory}\label{sec:theory}
 We derive asymptotic theory for the estimators of the model parameters and functionals of interest.

\subsection{Limit distributions} We first introduce some notation that is useful to state the assumptions that we make to derive the limit distribution of the estimators. Let  $\widetilde S_1 :=   \partial_{\pi} L_1(\pi)$ and  $\widetilde S_{2y} :=  \partial_{\theta_y} L_2(\theta_y, \pi)$ be the scores of the first and second steps   in Algorithm \ref{alg:tsdr} evaluated at the true parameter values, and  $H_1 := {\mathrm{E}} \left[ \partial_{\pi \pi'} L_1(\pi) \right]$ and $H_{2y} := {\mathrm{E}} \left[ \partial_{\theta_y \theta_y} L_2(\theta_y, \pi) \right]$ be the corresponding expected Hessians. Let
\begin{equation}\label{eq:cf}
\Sigma_{\theta_y\theta_{\tilde y}} := H_{2y}^{-1} \left\{n {\mathrm{E}} \left[ \widetilde S_{2y} \widetilde S_{2\tilde y}' \right] - J_{21y} H_1^{-1} J_{21\tilde y}' \right\} H_{2\tilde y}^{-1},
\end{equation}
where  $J_{21y} := {\mathrm{E}} \left[ \partial_{\theta_y \pi'} L_2(\theta_y, \pi) \right]$, $d_{\pi} := \dim \pi$,  and $d_{\theta} := \dim \theta_y$.

\begin{assumption}[DR Estimator with Sample Selection]\label{ass:fclt} (1) Random sampling: $\{(D^*_i, Y^*_i,Z_i)\}_{i=1}^n$ is a sequence of independent and identically distributed copies of $(D^*,Y^*,Z)$. We observe $D = 1(D^* > 0)$ and $Y = Y^*$ if $D=1$.
(2) Model: the distribution of $(D^*,Y^*)$ conditional on $Z$ follows the DR model \eqref{eq:drm}. (3) The support of $Z$, $\mathcal{Z}$,  is a  compact set.  (4) The support of $Y$ is either finite or a bounded interval. In the second case, the density function of $Y$ conditional on $X$ and $D=1$, $f_{Y \mid X,D}(y \mid x, 1)$, exists, is uniformly bounded above, and is uniformly continuous in $(y,x)$ on the support of $(Y,X)$ conditional on $D=1$
(5) Identification and non-degeneracy: the equations ${\mathrm{E}}[\partial_{\pi} L_1(\tilde \pi)] = 0$ and ${\mathrm{E}}[\partial_{\theta_y} L_2(\tilde \theta_y, \tilde \pi)] = 0$ posses a unique solution at $(\tilde \pi, \tilde \theta_y) = (\pi,\theta_y)$ that lies in the interior of a compact set $\Pi \times \Theta \subset \mathbb{R}^{d_{\pi} + d_{\theta}}$ for all $y \in \mathcal{Y}$; and the matrices $H_1$, $H_{2y}$  and $\Sigma_{\theta_y\theta_{y}}$ are nonsingular for each $y \in \mathcal{Y}$.
\end{assumption}

Part (1) is a standard condition about the sampling and selection process, which is designed for cross sectional  data.  Part (2) imposes the semiparametric DR model on the LGR of the conditional distribution of $(D^*,Y^*)$ at $d=0$.  Part (3) imposes some compactness conditions, which can be generalized at the cost of more complicated proofs.  Part (4) covers continuous, discrete and mixed continuous-discrete outcomes. Part (5) imposes directly identification and that the variance-covariance matrix of the first-step estimator and the covariance function of the second-step estimator are well-behaved. Note that $H_1$, $H_{2y}$ and $J_{21y}$ are finite by Part (3).
More primitive conditions for part (5) can be found in the conditional maximum likelihood literature, e.g.,  \citeasnoun{newey94}.

The main result of this section is a functional central limit theorem for $\hat{\theta}_y$.  Let $\ell^{\infty}(\mathcal{Y})$ be the set of bounded functions on $\mathcal{Y}$, and $\leadsto$ denote weak convergence (in distribution).
\begin{theorem}[FCLT for $\hat{\theta}_y$]\label{thm:fclt} Under Assumption \ref{ass:fclt},
$$
\sqrt{n} (\hat \pi - \pi) = -H_1^{-1} \widetilde{S}_1 + o_P(1) \leadsto Z_{\pi} \sim \mathcal{N}(0, -H_1^{-1}), \text{ in } \mathbb{R}^{d_{\pi}}
$$
and
$$
\sqrt{n} (\hat \theta_y - \theta_y)= - H_{2y}^{-1} \sqrt{n} \left( \widetilde S_{2y} - J_{21y} H_1^{-1} \widetilde S_1 \right) + o_P(1) \leadsto Z_{\theta_y} \text{ in } \ell^{\infty}(\mathcal{Y})^{d_{\theta}},
$$
where  $y \mapsto Z_{\theta_y}$ is a zero-mean Gaussian process with uniformly continuous sample paths and covariance function $
\Sigma_{\theta_y\theta_{\tilde y}}$,  $y,\tilde y \in \mathcal{Y}$, defined in \eqref{eq:cf}.
\end{theorem}
The first order term in the limit of $\sqrt{n} (\hat \theta_y - \theta_y)$ is the sample average of the influence function of $\hat \theta_y$. We  construct an estimator of the covariance function $
\Sigma_{\theta_y\theta_{\tilde y}}$ based on this function. Thus, we form
\begin{equation}\label{eq:se}
\widehat \Sigma_{\theta_y\theta_{\tilde y}} =  n^{-2} \sum_{i=1}^{n}  \widehat \psi_{i}(\hat \theta_y, \hat \pi) \widehat \psi_{i}(\hat \theta_{\tilde y}, \hat \pi)'.
\end{equation}
Here, $\widehat \psi_{i}$ is an estimator of the influence function of $\hat \theta_y$,
\begin{equation}\label{eq:if}
\widehat \psi_{i}(t,c) = - \widehat H_{2y}(t,c)^{-1} \left( S_{2yi}(t,c) - \widehat J_{21y}(t,c) \widehat H_1(c)^{-1} S_{1i}(c) \right),
\end{equation}
where $S_{1i}(c)$ and $S_{2iy}(t,c)$ are the individual scores of the first and second steps of Algorithm \ref{alg:tsdr},
\begin{eqnarray*}
S_{1i}(c) &:=&  \partial_{c} L_{1i}(c), \  L_{1i}(c) :=  D_{i}\log{\Phi(Z_i^\prime c)}+(1-D_{i})\log{\Phi(-Z_i^\prime c)}, \\ S_{2yi}(t,c) &:=& \partial_{t} L_{2yi}(t,c), \  t=(b,d)\\
L_{2yi}(t,c) &:= & D_{i}  \left[ I_{yi} \log{\Phi_2\left(-X_i^\prime b,Z_i^\prime c; -\rho(x'd)\right)}
 +  (1-I_{yi}) \log{\Phi_2\left(X_i^\prime b,Z_i^\prime c; \rho(x'd) \right)}\right],
\end{eqnarray*}
and
\begin{eqnarray*}
\widehat H_1(c) :=  \partial_{c c'} L_1(c), \ \widehat H_{2y}(t,c) :=  \partial_{t t'} L_2(t, c),  \  \widehat J_{21y}(t,c) :=  \partial_{t c'} L_2(t, c),
\end{eqnarray*}
are estimators of $H_1$, $H_{2y}$, and $J_{21y}$ when evaluated at $c = \hat \pi$ and $t = \hat \theta_y$.

We now establish a functional central limit theorem for the estimators of functionals of the model parameters. This result is based on expressing the functional as a suitable operator of the model parameters and using the functional delta method \cite[Chapter 3.9]{vdV-W}. To present the result in a concise manner, we consider a generic functional $$u \mapsto \Delta_u = \varphi_u(\pi, \theta_\cdot, F_Z),$$
where $u \in \mathcal{U}$, a totally bounded metric space, and $\varphi_u$ is an operator that maps $\mathbb{D}_{\Delta}$ to the set $\ell^\infty(\mathcal{U})$, where $\Delta_\cdot$ takes values. Here $\mathbb{D}_{\Delta}$ denotes the space for  the parameter tuple $(\pi, \theta_\cdot, F_Z)$; this space is not stated here explicitly, but is restricted by the regularity conditions of the previous section.  Here
we identify $F_Z$ with an integral operator $f \mapsto \int f(z) d F_Z(z)$ taking values in $\ell^\infty(\mathcal{F})$ that acts on a Donsker set of bounded measurable functions $\mathcal{F}$, which includes indicators of rectangular sets;  see \citeasnoun{chernozhukov+13inference} and examples below.
The parameter space $\mathbb{D}_{\Delta}$ is a subset of a normed space $\mathbb{D}:=\mathbb{R}^{d_{\pi}} \times \ell^{\infty}(\mathcal{Y})^{d_{\theta}} \times \ell^{\infty}(\mathcal{F}) $.  In this notation, the plug-in estimator of the functional $\Delta_u$ is $$\hat \Delta_u = \varphi_u(\hat \pi, \hat \theta_y, \hat F_Z),$$
where $\hat \pi$ and $\hat \theta_y$ are the estimators of the parameters defined in Algorithm \ref{alg:tsdr} and $\hat F_Z$ is the empirical distribution of $Z$.

We provide some examples.  The distribution of the latent outcome is given by:
$$
F_{Y^*}(y) = \varphi_y(\pi, \theta_y, F_Z) = \int \Phi(-x'\beta_y) dF_Z(z),
$$
$\mathcal{F}$  contains $\{\Phi(-\cdot'\beta_y) : y \in \mathcal{Y} \}$ as well as the indicators of all rectangles in $\overline{\mathbb{R}}^{d_z}$, $\overline{\mathbb{R}} := \mathbb{R} \cup \{-\infty, +\infty\}$, $d_z = \dim Z$,  and $\mathcal{U} = \mathcal{Y}$.
The quantile function of the latent outcome is $$Q_{Y^*}(\tau) = \varphi_{\tau}(\pi, \theta_y, F_Z) =   \mathbf{Q}_{\tau} \mathbf{R} F_{Y^*},$$ $\mathcal{F}$  is the same as for the distribution of the latent outcome, $\mathcal{U}$ is a closed subset of $(0,1)$ including the quantile indexes of interest, $\mathbf{R}$ is the non-decreasing rearrangement operator, and $\mathbf{Q}_{\tau}$ is the left-inverse (quantile) operator. The distribution of the observed outcome is given by:
$$
   F_Y(y \mid D=1) = \varphi_y(\pi, \theta_y, F_Z) = \frac{\int \Phi_2\left( -x^\prime\beta(y), z^\prime \pi; -\rho(x'\delta(y)) \right) dF_Z(z)}{\int \Phi(z^\prime \pi) dF_Z(z)},
$$
$\mathcal{F}$  contains $\{\Phi_2\left( -\cdot^\prime\beta(y), \cdot^\prime \pi; -\rho(\cdot'\delta(y)) \right) : y \in \mathcal{Y} \}$ as well as the indicators of all rectangles in $\overline{\mathbb{R}}^{d_z}$,  and $\mathcal{U} = \mathcal{Y}$.

The following result is a corollary of Theorem \ref{thm:fclt} by the functional delta method. Let $UC(\mathcal{Y},\xi)$ be the set of functions on $\mathcal{Y}$ that are uniformly continuous with respect to $\xi$, a standard metric on $\mathbb{R}$, and $UC(\mathcal{F},\lambda)$ be the set of functionals on $\mathcal{F}$ that are uniformly continuous with respect to $\lambda$, where $\lambda(f,\tilde f) = [\Pr(f - \tilde f)^2]^{1/2}$ for any $f, \tilde f \in \mathcal{F}$.

\begin{corollary}[FCLT for $\hat \Delta_u$]\label{cor:fclt} Suppose that Assumption \ref{ass:fclt} holds, and $(p, t_y, F) \mapsto \varphi_\cdot(p, t_y, F)$, from $\mathbb{D}_{\Delta} \subset \mathbb{D}$ to $ \ell^\infty(\mathcal{U})$  is Hadamard differentiable at $(\pi, \theta_y, F_Z)$, tangentially to $\mathbb{R}^{d_{\pi}} \times UC(\mathcal{Y},\xi)^{d_{\theta}} \times UC(\mathcal{F},\lambda)$ with derivative $(p, t_y, F) \mapsto \varphi_\cdot'(p, t_y, F)$ that is defined and continuous on $\mathbb{R}^{d_{\pi}} \times \ell^{\infty}(\mathcal{Y})^{d_{\theta}} \times \ell^{\infty}(\mathcal{F})$. Then,
$$
\sqrt{n}(\hat \Delta_u -  \Delta_u) \leadsto Z_{\Delta_u} := \varphi_u'(Z_{\pi}, Z_{\theta_y}, Z_{F}) \text{ in } \ell^{\infty}(\mathcal{U}),
$$
where $Z_{\pi}$  and $Z_{\theta_y}$ are the random limits in Theorem \ref{thm:fclt}, $Z_{F}$ is a tight $F_Z$-Brownian bridge, and $u \mapsto  Z_{\Delta_u}$ is a tight zero-mean Gaussian process.
\end{corollary}

\begin{remark}[Hadamard Differentiable Functionals]\label{remark:hd} The distributions of the latent and observed outcome together with counterfactual distributions constructed thereof are examples of Hadamard differentiable functions. In the case of the latent outcome, the result follows from the Hadamard differentiability of the counterfactual operator in \citeasnoun{chernozhukov+13inference}. In the case of the observed outcome, the result follows from the differentiability of the counterfactual operator and the composition rule for Hadamard derivatives applied to the ratio of two functions. Quantile (left-inverse) functionals of these distributions are Hadamard differentiable under additional conditions that guarantee that the quantile operator is Hadamard differentiable. These include that the outcome variable be continuous with density bounded above and away from zero \cite{cfg10}. Then the Hadamard differentiability of the quantile function follows from the composition rule for Hadamard derivatives.\qed
\end{remark}

\begin{remark}[Inference on Quantile Functions]\label{remark:qinfer} There are two alternatives to construct confidence bands for quantile functions. The first approach is the standard method based on characterizing the limit distribution of the estimator of the quantile function using the delta method, which relies on the Hadamard differentiability of the inverse operator. As we mention in Remark \ref{remark:hd}, this differentiability requires of additional conditions including that the outcome variable be continuous. The second approach applies to any type of outcome variable. It is based on the generic method of \citeasnoun{cfmw16} that inverts confidence bands for distribution functions into confidence bands of quantile function. This method does not rely on the delta method and is therefore more robust to modeling assumptions and widely applicable.  It  has the  shortcoming, however,  that the bands might not be centered at the point estimate of the quantile function. We apply the second method to obtain most of the results in the empirical application.\qed
\end{remark}

\subsection{Multiplier Bootstrap} We make the following assumption about the bootstrap multipliers of Algorithm \ref{alg:mb}:
\begin{assumption}[Multiplier Bootstrap]\label{ass:mb} The multipliers  $(\omega_{1},...,\omega_{n})$ are i.i.d. draws from a
random variable $\omega \sim \mathcal{N}(0,1)$, and are
independent of  $\{(D^*_i, Y^*_i,Z_i)\}_{i=1}^n$ for all $n$.
\end{assumption}

Let
$$
\hat{\theta}_y^b  = \hat \theta_y + n^{-1} \sum_{i = 1}^n \omega_i \ \widehat \psi_i(\hat \theta_y, \hat \pi)
$$
be the multiplier bootstrap version of  $\hat{\theta}_y$. We establish a functional central limit theorem for the bootstrap for  $\hat{\theta}_y$.
  Here we use $\leadsto_{\Pr}$ to denote
bootstrap consistency, i.e. weak convergence conditional on the data in
probability, which is formally defined in Appendix \ref{app:notation}.

\begin{theorem}[Bootstrap FCLT for $\hat{\theta}_y$]\label{thm:bfclt} Under the conditions of Theorem \ref{thm:fclt} and Assumption \ref{ass:mb},
$$
\sqrt{n} (\hat \theta_y^b - \hat \theta_y) \leadsto_{\Pr} Z_{\theta_y} \text{ in } \ell^{\infty}(\mathcal{Y})^{d_{\theta}},
$$
where  $y \mapsto Z_{\theta_y}$ is the same Gaussian process as in Theorem \ref{thm:fclt}.
\end{theorem}

The following result is a corollary of Theorem \ref{thm:bfclt} by the functional delta method for the bootstrap \cite[Chapter 3.9]{vdV-W}. Let $\hat \Delta^b_u = \varphi_u(\hat \pi^b, \hat \theta^b_y, \hat F^b_Z)$, be the multiplier bootstrap version of $\hat \Delta_u$ where $$\hat \pi^b = \hat \pi - n^{-1} \sum_{i=1}^n \omega_i  \ \hat{H}_1(\hat \pi)^{-1} S_{1i}(\hat \pi),$$ and $\hat F^b_Z$ is the weighted empirical distribution of $Z$  that uses  $(1+\omega_1, \ldots, 1+\omega_n)$ as sampling weights.

\begin{corollary}[Bootstrap FCLT for $\hat \Delta_u$] Suppose that the conditions of Corollary \ref{cor:fclt} and Assumption \ref{ass:mb} hold. Then,
$$
\sqrt{n}(\hat \Delta^b_u - \hat \Delta_u) \leadsto_{\Pr} Z_{\Delta_u} \text{ in } \ell^{\infty}(\mathcal{U}),
$$
where $Z_{\Delta_u}$ is the same process as in Corollary \ref{cor:fclt}.
\end{corollary}

\section{Wage Decompositions in the UK}\label{sec:empirics}  We apply the DR model with sample selection to carry out wage decompositions accounting for endogenous employment participation using data from the United Kingdom.

\subsection{Data} The data come from the U.K. Family Expenditure Survey (FES) for the years 1978 to 2001, Expenditure and Food Survey (EFS) for the years 2002 to 2007, and Living Costs and Food Survey (LCFS) for the years 2008 to 2013. Despite the differences in the name, these surveys contain comparable information. Indeed, the FES was combined to the National Food Survey to form the EFS, which was renamed LCFS when it became a module of the Integrated Household Survey. The data from the FES has been previously used by \citeasnoun{gosling00}, \citeasnoun{brs03}, \citeasnoun{bgim07} and \citeasnoun{ab17} to study wage equations in the U.K. labor market. We are not aware of any previous use of the data from the EFS and LCFS for this purpose.\footnote{See \citeasnoun{rv18} for another recent application of the data to the analysis of female labour force participation. }
The three surveys contain repeated cross-sectional observations for women and men.  The selection of the sample is similar to the previous work that used the FES. Thus, we keep individuals with ages between 23 to 59 years, and drop full-time students, self-employed workers, those married with spouse absent, and those with missing education or employees whose wages are missing. This leaves a sample of 258,900 observations, 139,504 of them correspond to women and 119,765 to men. The sample size per survey year and gender ranges from 2,197 to 4,545.

The outcome of interest, $Y$, is the logarithm of real hourly wage rate. We construct this variable as the ratio of the weekly usual gross main nominal earning to the weekly usual working hours, deflated by the U.K. quarterly retail price index. The selection variable, $D$, is an indicator for being employed.\footnote{For data before 1990, $D=0$ if the individual is in one of the following status: seeking work, sick but seeking work, sick but not seeking work, retired and unoccupied. For those in and after 1990, $D=0$ if the individual is seeking work and available, waiting to start work, sick or injured, retired or unoccupied.} The covariates, $X$, include 5  indicators for age when ceasing school ($\leq$15, 16, 17--18, 19--20, 21--22 and $\geq$ 23), a quartic polynomial in age, an indicator of being married or cohabiting, 6 variables with the number of kids by age categories (1, 2, 3--4, 5--10, 11--16,  and 17-18), 36 survey year indicators, and 11 region indicators (Northern 5.48\%, Yorkshire 9.56\%, North Western 10.20\%, East Midlands 7.36\%, West Midlands 9.13\%, East Anglia 5.31\%, Greater London 10.06\%, South Eastern 16.82\%, South Western 7.94\%, Wales 4.99\%, Scotland 8.92\%, and Northern Ireland 4.23\%).\footnote{In the rest of the paper we shall refer to an individual being married or cohabiting as married.}

The excluded covariate, $Z_1$,  is a potential out-of-work income benefit interacted with the marital status indicator used before in  \citeasnoun{brs03} and \citeasnoun{bgim07}. This benefit is constructed with the Institute for Fiscal Studies (IFS) tax and welfare-benefit model (TAXBEN). TAXBEN is a static tax and benefit micro-simulation model of taxes on personal incomes, local taxes, expenditure taxes, and entitlement to benefits and tax credits that operates on large-scale, representative, household surveys \cite{b09}. It is designed to calculate the income of a tax unit if the individual considered were out of work.\footnote{Our definition of the out-of-work benefit income is slightly different from  the definition of \citeasnoun{brs03} and \citeasnoun{bgim07}. They calculated it as the income of a tax unit if all the individuals within the tax unit were out of work. In our view our definition might better reflect the opportunity cost or outside value option of working that the individual faces.} It is composed of eligible unemployment and housing benefits, which are determined by the demographic composition of the tax unit and the housing costs that the tax unit faces. These costs vary by region and over time due to numerous policy changes that have occurred over time. There is no consensus in the literature about the validity of this variable as excluded covariate. In our case the outcome and exclusion restrictions imply that, conditional on the observed covariates, the offered wage and dependence between offered wage and net reservation wage do not depend on the level of the benefit.  We shall assume that the exclusion restrictions are satisfied  and refer to \citeasnoun{brs03} and \citeasnoun{bgim07} for a discussion on the plausibility of the outcome restriction.  In Appendix \ref{app:wage}, we provide an example of a model for offered and reservation wages  that satisfies the exclusion restrictions.

Table \ref{table:summstat} reports means and standard deviations of all the variables used in the analysis. We report these statistics for the entire sample, and by employment status and gender. The overall employment rate is 74\%. Women are 17\% less likely to be employed than men, and the unconditional gender wage gap is 33\%. Overall, women and men are similar in terms of covariates. Both working men and women are relatively more highly educated, younger, and more likely to be married than their non-working counterparts. Having young children and high out-of-work benefits is negatively associated with employment for women but not for men.

\begin{footnotesize}
\begin{table}
\caption{Summary Statistics}\label{table:summstat}
\def\sym#1{\ifmmode^{#1}\else\(^{#1}\)\fi}
\begin{tabular}{l*{6}{c}}
\hline\hline
                &\multicolumn{2}{c}{Full}         &\multicolumn{2}{c}{Male}         &\multicolumn{2}{c}{Female}    \\
                &     All        &  Employed         &     All        &   Employed         &   All         & Employed \\
\hline
Log Hourly Wage &                  &    2.38         &                  &    2.54         &                  &    2.21         \\
                &                  &  (0.54)         &                  &  (0.51)         &                  &  (0.52)         \\
Employed        &    0.74         &                  &    0.83         &                  &    0.66        &                  \\
                &  (0.44)         &                  &  (0.38)         &                  &  (0.47)         &                \\
\multicolumn{7}{l}{Ceased School at} \\
\hspace{3mm} $\leq$ 15       &    0.33         &    0.30         &    0.33         &    0.31         &    0.33         &    0.29         \\
                &  (0.47)         &  (0.46)         &  (0.47)         &  (0.46)         &  (0.47)         &  (0.45)         \\

\hspace{3mm}16       &    0.31         &    0.30         &    0.32         &    0.32         &    0.30         &    0.29         \\
                &  (0.46)         &  (0.46)         &  (0.47)         &  (0.47)         &  (0.46)         &  (0.45)         \\

\hspace{3mm}17-18     &    0.18         &    0.19         &    0.16         &    0.17         &    0.20         &    0.22         \\
                &  (0.38)         &  (0.39)         &  (0.37)         &  (0.37)         &  (0.40)         &  (0.41)         \\

\hspace{3mm}19-20      &   0.04         &   0.05         &   0.04         &   0.04         &   0.04         &   0.05         \\
                &  (0.20)         &  (0.21)         &  (0.20)         &  (0.20)         &  (0.21)         &  (0.22)         \\

\hspace{3mm}21-22      &   0.09         &    0.11         &   0.09         &    0.10         &   0.09         &    0.12         \\
                &  (0.29)         &  (0.31)         &  (0.29)         &  (0.30)         &  (0.29)         &  (0.32)         \\

\hspace{3mm}$\geq$23     &   0.05         &   0.05         &   0.06         &   0.06         &   0.04         &   0.04         \\
                &  (0.21)         &  (0.22)         &  (0.23)         &  (0.24)         &  (0.19)         &  (0.20)         \\
Age            &    40.13         &    39.84         &    40.22         &    39.76         &    40.06         &    39.92         \\
                &  (10.43)         &  (10.10)         &  (10.40)         &  (10.11)         &  (10.45)         &  (10.08)         \\

Married        &    0.76         &    0.79         &    0.78         &    0.81         &    0.75         &    0.76         \\
                &  (0.43)         &  (0.41)         &  (0.42)         &  (0.39)         &  (0.43)         &  (0.43)         \\
\multicolumn{7}{l}{Number of children with age} \\
\hspace{3mm} 0--1          &   0.06         &   0.05         &   0.06         &   0.06         &   0.06         &   0.03         \\
                &  (0.24)         &  (0.22)         &  (0.24)         &  (0.25)         &  (0.24)         &  (0.18)         \\

\hspace{3mm} 2          &   0.05         &   0.04         &   0.05         &   0.05         &   0.05         &   0.03         \\
                &  (0.23)         &  (0.21)         &  (0.23)         &  (0.23)         &  (0.23)         &  (0.18)         \\

\hspace{3mm} 3--4         &    0.10         &   0.09         &   0.10         &    0.10         &    0.11         &   0.07         \\
                &  (0.32)         &  (0.29)         &  (0.31)         &  (0.32)         &  (0.32)         &  (0.27)         \\

\hspace{3mm} 5--10        &    0.32         &    0.29         &    0.30         &    0.30         &    0.33         &    0.28         \\
                &  (0.64)         &  (0.61)         &  (0.62)         &  (0.62)         &  (0.65)         &  (0.59)         \\

\hspace{3mm} 11--16       &    0.30         &    0.30         &    0.28         &    0.29         &    0.32         &    0.32         \\
                &  (0.63)         &  (0.62)         &  (0.61)         &  (0.61)         &  (0.64)         &  (0.63)         \\

\hspace{3mm} 17--18       &   0.03         &   0.04         &   0.03         &   0.03         &   0.04         &   0.04         \\
                &  (0.19)         &  (0.19)         &  (0.18)         &  (0.18)         &  (0.19)         &  (0.20)         \\

Benefit Income  &    5.44         &    5.50         &    5.25         &    5.29         &    5.60         &    5.73         \\
                &  (0.74)         &  (0.78)         &  (0.70)         &  (0.72)         &  (0.73)         &  (0.78)         \\
\hline
Observations    &   258,900         &   190,765         &   119,396         &    98,764         &   139,504         &    92,001         \\
\hline\hline
\multicolumn{7}{l}{\footnotesize \emph{Notes}: all the entries are means with standard deviations in parentheses.}\\
\multicolumn{7}{l}{\footnotesize \emph{Source:} FES/EFS/LCFS Data.}\\
\end{tabular}

\end{table}

\end{footnotesize}

Figure \ref{fig:trends} provides some background on the U.K. labor market using our data. The left panel shows that  over 36 years the average wages of working men and women have continuously grown and the unconditional gender wage gap  has progressively narrowed from $46\%$ to $24\%$. The middle panel indicates that the growth of average wage has come together with an increase in wage inequality for both working men and women until 2000. The positive trend in wage inequality has continued for men after 2000, but not for women. The right panel shows opposite trends in the employment rate for men and women, where the gender employment gap has steadily and sharply reduced from 34\% to 8\%.

\begin{figure}    \includegraphics[width=0.32\textwidth,page=1]{Inequality.pdf}
    \includegraphics[width=0.32\textwidth,page=2]{Inequality.pdf}
    \includegraphics[width=0.32\textwidth,page=3]{Inequality.pdf}
\caption{Trends in U.K. labor market 1978-2013 by gender: left panel reports the average of the log wage rate, the middle panel reports the 90-10 percentile  spread of the log wage rate, and the right panel reports the employment rate}\label{fig:trends}
\end{figure}

\subsection{Empirical Specifications} We estimate the DR model for different samples and carry out several wage decompositions where we compare the distributions of men and women, or the distributions over time within genders. The specifications of the selection and outcome  equations include all the covariates  described above except for the excluded covariates in the outcome equation. The parameter of the selection sorting function is notoriously more difficult to estimate than the parameters of the selection and outcome equations. We consider four simplified specifications of the sorting function where  the covariates included in the index $X'\delta(y)$ are:
\begin{itemize}
\item Specification 1: a constant.
\item Specification 2: a constant and the marital status indicator.
\item Specification 3: a constant and a linear trend on the year of the survey.
\item Specification 4: a constant and a linear trend on the year of the survey interacted with the marital status indicator.
\end{itemize}
We also experimented with other specifications that include the education indicators, indicators of survey year, or age.  We do not report these results because they do not show any clear pattern mainly due to imprecision in the estimation of the parameter $\delta(y)$.\footnote{The main results on the coefficients of the wage equation and  wage decompositions presented below are not sensitive to the specification of the sorting equation.}

\subsection{Model Parameters}
We report point estimates and 95\% confidence bands for the coefficients of the education and marital status indicators in the outcome equation and  the correlation function in the selection sorting. Estimates and 95\% confidence bands for the coefficients of the selection equation, coefficients of the fertility indicators in the outcome equation and coefficients in the selection sorting function are given in the SM. The estimates are obtained with Algorithm \ref{alg:tsdr} replacing $\mathcal{Y}$ by a finite grid containing the sample quantiles of log real hourly wage with indexes $\{0.10, 0.11, \ldots, 0.90 \}$ in the pooled sample of men and women. We report all the estimates as a function of the quantile index. The confidence bands are constructed by Algorithm \ref{alg:mb} with $B=500$ bootstrap repetitions and the same finite grid  as for the estimates. We also report  estimates from the HSM of Example \ref{ex:hsm} with dash lines as a benchmark of comparison.\footnote{We report  estimates of $\beta/\sigma$ in the outcome equation of the HSM for comparability; see Example \ref{ex:hsm}.}

The estimates of the coefficients of the education and marital status indicators in the outcome equation are reported in Figure \ref{fig:coef_men_spec1} for men and Figure \ref{fig:coef_women_spec1} for women. These estimates correspond to specification 1. Estimates for specifications 2--4 are given in the SM. For all the specifications and genders, we find that the returns to education are heterogenous across the distribution and broadly increasing in the years of education (age leaving school). In contrast, the classical HSM model does not allow the heterogeneity and estimates constant coefficients. This restriction is rejected empirically, lending support for
the generalized HSM. The coefficient of the marital status indicator is uniformly positive for men, whereas is negative but mainly statistically not different from zero for women. We cannot reject that this coefficient is homogeneous across the distribution for both men and women.\footnote{We find more heterogeneity in the coefficient of the marital status indicator in the specifications 2 and 4 that include marital status in the selection sorting function.}

\begin{figure}
    \includegraphics[height=0.75\textheight, width=0.85\textwidth,page=1]{Estimates_11378_500.pdf}
    \caption{Estimates and 95\% confidence bands for coefficients of education and marital status in the outcome equation: specification 1 for men}\label{fig:coef_men_spec1}
\end{figure}

\begin{figure}
    \includegraphics[height=0.75\textheight, width=0.85\textwidth,page=3]{Estimates_11378_500.pdf}
    \caption{Estimates and 95\% confidence bands for coefficients of education and marital status in the outcome equation:  specification 1 for women}\label{fig:coef_women_spec1}
\end{figure}

Figures \ref{fig:rho_men_spec1}--\ref{fig:rho_men_spec4} display the estimates of the sorting effect functions for specifications 1--4, respectively. The estimates of the coefficients of these functions for specifications 3 and 4 are given in the SM. Figure \ref{fig:rho_men_spec1} shows  positive selection sorting for men and negative selection sorting for women.  In both cases we cannot reject that the sorting is constant across the distribution. This finding is refined in Figure \ref{fig:rho_men_spec2}, where we  uncover that the positive male sorting comes mainly from bachelors, whereas the negative female sorting comes from married women.  This pattern is consistent with a marriage market where there is assortative matching in offered wages given observable characteristics, where women with high potential wages are married to highly paid working men and decide not to work \cite{neal04}. Figure \ref{fig:rho_men_spec3} shows that the sorting homogeneity found in the pooled sample hides some heterogeneity across time. Thus, we find that the male sorting is heterogeneous in the early years,  negative at the bottom and positive  at the top  of the distribution, and progressively becomes homogenous. The female sorting is more homogenous over time, but also displays a positive trend, specially at the bottom of the distribution. Figure \ref{fig:rho_men_spec4} shows that the trends in sorting are driven by married individuals at the bottom of the distribution and single individuals at the top of the distribution.\footnote{We do not report confidence bands for specifications 3 and 4 to avoid cluttering. The confidence bands for the coefficients of the selection sorting function in the SM show that the results on the trends are statistically significant.}

\begin{figure}
    \includegraphics[height=0.5\textwidth, width=0.49\textwidth,page=2]{Estimates_11378_500.pdf}
    \includegraphics[height=0.5\textwidth, width=0.49\textwidth,page=4]{Estimates_11378_500.pdf}
    \caption{Estimates and 95\% confidence bands for the selection sorting function: specification 1}\label{fig:rho_men_spec1}
\end{figure}

\begin{figure}
    \includegraphics[height=.8\textwidth, width=0.4\textwidth,page=2]{Estimates_11378_couple_500.pdf}
    \includegraphics[height=.8\textwidth, width=0.4\textwidth,page=4]{Estimates_11378_couple_500.pdf}
    \caption{Estimates and 95\% confidence bands for the selection sorting function: specification 2}\label{fig:rho_men_spec2}
\end{figure}

\begin{figure}
    \includegraphics[height=0.5\textwidth, width=0.49\textwidth,page=2]{Estimates_11378_lintime_500.pdf}
    \includegraphics[height=0.5\textwidth, width=0.49\textwidth,page=5]{Estimates_11378_lintime_500.pdf}
    \caption{Estimates and 95\% confidence bands for the selection sorting function: specification 3}\label{fig:rho_men_spec3}
\end{figure}

\begin{figure}
    \includegraphics[height=.8\textwidth, width=0.4\textwidth,page=2]{Estimates_11378_lintimexcouple_500.pdf}
    \includegraphics[height=.8\textwidth, width=0.4\textwidth,page=5]{Estimates_11378_lintimexcouple_500.pdf}
    \caption{Estimates and 95\% confidence bands for the selection sorting function: specification 4}\label{fig:rho_men_spec4}
\end{figure}

\subsection{Distributions of Offered and Observed Wages, and Wage Decompositions}
Figure \ref{fig:latent} shows point estimates of the quantiles of offered and observed wages for men and women based on specification 1. Estimates for the other specifications and confidence bands for all the specifications are given in the SM.  The offered wage is a latent variable defined for all the individuals that is free of sample selection. As we showed in Section \ref{sec:model}, the distributions of both types of wages can be expressed as functionals of the model parameters, and estimated using the plug-in estimators \eqref{eq:ldisthat} and \eqref{eq:odisthat}.\footnote{The model-based estimator of the observed distribution in  \eqref{eq:odisthat} produces almost identical estimates to the empirical distribution of the observed wages.} We find opposite signs in the sample selection bias for men and women.  The quantiles of the observed wages are below the quantiles of latent wages for men, but the opposite holds for women.  This pattern is consistent with the sign of the estimates of the selection sorting function, where we found positive sorting for men and negative sorting for women.

\begin{figure}
    \includegraphics[height=0.49\textwidth, width=0.49\textwidth,page=1]{Comp_Dist.pdf}
    \includegraphics[height=0.49\textwidth, width=0.49\textwidth,page=2]{Comp_Dist.pdf}
    \caption{Estimates of the quantiles of observed and offered (latent) wages: specification 1}\label{fig:latent}
\end{figure}

Figure \ref{fig:latent2} compares the quantile function of offered wages  between men and women and carries out a gender discrimination analysis  based on specification 1. The estimates and 95\% confidence bands for the other specifications are reported in the SM. The gender discrimination analysis is based on the counterfactual distributions
$$
F_{Y^* \langle j, k \rangle}(y) = \int \Phi(-x'\beta_j(y)) dF_{X_k}(x),
$$
where $\beta_j(y)$ is the coefficient of the wage equation in group $j$, $F_{X_k}$ is the distribution of the characteristics in group $k$, and $j$ and $k$ are group indices for women and men.  $F_{Y^* \langle j, k \rangle}$ corresponds to the distribution of offered wages that we would observed when the wage structure is as in group $j$ and the distribution of  characteristics is as in group $k$. We decompose the  difference in the quantile functions of the latent wages between women (group 1) and men (group 0) using the counterfactual distributions as
$$
F_{Y^* \langle 1, 1 \rangle} -  F_{Y^* \langle 0, 0\rangle}  = [F_{Y^* \langle 1, 1 \rangle} -  F_{Y^* \langle 0, 1 \rangle}] + [F_{Y^* \langle 0, 1 \rangle} -  F_{Y^* \langle 0, 0 \rangle}],
$$
where the first term is the wage structure or discrimination effect and the second term is the composition effect. We obtain estimates of the counterfactual distributions and quantile functions using the plug-in estimator in \eqref{eq:ldisthat} and the operators $ \mathbf{Q}_{\tau}$ and $\mathbf{R}$.
We find that the wages offered to women are between $21$ and $40\%$ lower than the wages offered to men at the same quantile index. The majority of this difference  is explained by differences in the wage structure, $\beta(y)$, whereas differences in composition, $F_Z$, have very little explanatory power. This result can be interpreted as evidence of gender discrimination in the labor market.

\begin{figure}
    \includegraphics[height=0.49\textwidth, width=0.49\textwidth,page=17]{DFQF_11378_500.pdf}
    \includegraphics[height=0.49\textwidth, width=0.49\textwidth,page=21]{DFQF_11378_500.pdf}
    \caption{Estimates and 95\% confidence bands for the quantiles of offered (latent) wages and decomposition between women and men: specification 1}\label{fig:latent2}
\end{figure}

 We next use the DR model to decompose changes in the distribution of the observed wage between women and men, and between the first and second halves of the sample period for each gender. We extract four components that correspond to different inputs of the DR model:
\begin{enumerate}
\item Selection (employment) sorting: $\delta(y)$.
\item Selection (employment) structure: $\pi$.
\item Outcome (wage) structure: $\beta(y)$.
\item Composition: $F_Z$.
\end{enumerate}
To define the effects of these components, let $F_{Y\langle t, s, r, k\rangle}$ be the counterfactual distribution of wages when the sorting is as in group $t$, the employment structure is as in group $s$, the wage structure is as  in group $r$, and the composition of the population is as in group $k$. The actual distribution in group $t$ therefore corresponds to $F_{Y\langle t, t, t, t\rangle}$.  We assume that there are two groups indexed by $0$ and $1$ that correspond to demographic populations such as men and women, or time periods such as the first and second halves of the sample years. Then, we can decompose the distribution of observed wage between group $1$ and group $0$ as:
\begin{multline*}
F_{Y\langle 1, 1, 1, 1 \rangle} -  F_{Y\langle 0, 0, 0, 0 \rangle}  = [F_{Y\langle 1, 1, 1, 1 \rangle} -  F_{Y\langle 0, 1, 1, 1 \rangle}] + [F_{Y\langle 0, 1, 1, 1 \rangle} -  F_{Y\langle 0, 0, 1, 1 \rangle}] \\ + [F_{Y\langle 0, 0, 1, 1 \rangle} -  F_{Y\langle 0, 0, 0, 1 \rangle}] + [F_{Y\langle 0, 0, 0, 1 \rangle} -  F_{Y\langle 0, 0, 0, 0 \rangle}],
\end{multline*}
where the first term in square brackets of the right hand side is a sorting effect, the second an employment structure effect, the third a wage structure effect, and the forth a composition effect. This is a distributional version of the classical Oaxaca-Blinder decomposition that accounts for sample selection \cite{kitagawa55,oaxaca73,blinder73}. It is well-known that the order of extraction of the components in this type of decompositions might matter. As a robustness check, we estimate the decomposition changing the ordering of the components. In results not reported, we find that the main findings are not sensitive to the change of ordering.

In terms of the DR model, the counterfactual distribution can be expressed as the functional
$$
F_{Y\langle t, s, r, k\rangle}(y) = \frac{\int \Phi_2\left( -x^\prime\beta_r(y), z^\prime \pi_s; -\rho(x'\delta_t(y)) \right) dF_{Z_k}(z)}{\int \Phi(z^\prime \pi_s) dF_{Z_k}(z)},
$$
where $\delta_t$ is the coefficient of the sorting function in group $t$, $\pi_s$ is the coefficient of the employment equation in group $s$, $\beta_r$ is the coefficient of the wage equation in group $r$, and $F_{Z_k}$ is the distribution of characteristics in group $k$. Given random samples for groups $0$ and $1$, we construct a plug-in estimator of $F_{Y\langle t, s, r, k\rangle}$ by suitably combining the estimators of the model parameters and distribution of covariates from the two groups.

Figure \ref{fig:dec_spec1} reports estimates of the quantile functions of observed wages for men and women, together with the relative contributions of each component to the decomposition between men (group $0$) and women (group $1$) based on specification 1. The bands for the contributions are joint for all the components and rely on the delta method; see Remark \ref{remark:qinfer}. Estimates of the components of the decomposition and the analysis based on specifications 2--4 are given in the SM.  The distribution for men first order stochastically dominates the distribution for women.  Most of this gender wage gap is explained by differences in the wage structure, i.e. differences in the returns to observed characteristics that might be associated to gender discrimination. However,  differences in sorting and employment structure also account for an important percentage of the gap, specially at the top of the distribution. Thus, we uncover that the negative female sorting explains about 30--40\% of the gap at the top of the distribution. A possible explanation is that women with very high potential wages decide not to work  because there are no high-paid jobs available to them due to glass ceiling \cite{abv03}. The negative contribution of the employment structure can be explained by the order of the decomposition where we are applying the male employment structure to the female distribution with positive male sorting. In this case we are increasing the proportion of employed women, where the added women come from a pool with lower positive selection, and this negative effect is not reversed by a change in the composition of the working women; see Remark \ref{remark:se} for more details. The aggregate selection effect, defined as the sum of the selection sorting and selection structure effects, is positive and statistically significant at the top of the distribution; see Figure \ref{fig:adec_spec1} in the SM. Differences in the composition of the characteristics contribute very little to explain the gender gap.  Finally, the estimates from the HSM in dash lines pick up the average contributions of the components, but miss all the heterogeneity across the distribution.

\begin{figure}
    \includegraphics[height=.49\textwidth, width=.49\textwidth,page=14]{CondDFQF_11378_500.pdf}
        \includegraphics[height=.49\textwidth, width=.49\textwidth,page=43]{CondDFQF_11378_500.pdf}
    \caption{Estimates and 95\% confidence bands for the quantiles of observed wages and decomposition between men and women in specification 1}\label{fig:dec_spec1}
\end{figure}

Figures \ref{fig:dec2_women_spec1} and \ref{fig:dec2_men_spec1} report estimates of the quantile functions of observed wages for the first and second halves of the sample period, together with the relative contributions of each component to the decomposition between second half (group $0$) and first half (group $1$) based on specification 1 for women and men, respectively. Estimates of the components of the decompositions  are given in the SM. The distribution for the second half first order stochastically dominates the distribution for the first half in both cases. For women, the most important components are the wage structure and composition effects in this order. The importance of the wage structure is decreasing along the distribution, whereas the importance of the composition is increasing. Composition and wage structure  are also the most important components for men. The small contributions of the selection sorting component to the change in the distribution of wages between the two time period for both genders seem to contradict the linear time trends that we found in the coefficient of the sorting selection function. This might be explained by the inability of a coarse partition of the sample into two halves to capture the gradual increase in selection sorting, together with the changes in the composition.

\begin{figure}
    \includegraphics[height=.49\textwidth, width=.49\textwidth,page=14]{CondDFQF_fhalf_500.pdf}
        \includegraphics[height=.49\textwidth, width=.49\textwidth,page=44]{CondDFQF_fhalf_500.pdf}

    \caption{Estimates and 95\% confidence bands for the quantiles of observed wages and decomposition between first and second half of the sample period for women in specification 1}\label{fig:dec2_women_spec1}
\end{figure}

\begin{figure}
    \includegraphics[height=.49\textwidth, width=.49\textwidth,page=14]{CondDFQF_mhalf_500.pdf}
        \includegraphics[height=.49\textwidth, width=.49\textwidth,page=44]{CondDFQF_mhalf_500.pdf}

    \caption{Estimates and 95\% confidence bands for the quantiles of observed wages and decomposition between first and second half of the sample period for men in specification 1}\label{fig:dec2_men_spec1}
\end{figure}

\subsection{Discussion}
The main findings can be summarized as: (1)  heterogeneous positive effect of education and homogeneous effect of being married on offered wages for both genders; (2) positive sorting for men and negative sorting for women driven by single men and married women, which is consistent with assortative matching in the marriage market; (3) heterogeneity in selection sorting decreases gradually over time; (4) differences in returns to characteristics in the wage equation, which might  be associated to gender discrimination in the labor market,  account for most of the gender wage gap; (5) selection sorting on unobservables explains up to 39\% of the gender wage gap at the top of the distribution, which can be taken as evidence of glass ceiling; and (6) changes in the structure of the wage equation and composition of the characteristics account for most of the differences in the wage distribution between the two halves of the sample period within each gender.

We  compare and contrast these findings with previous results from the  literature that studied similar issues. These results were obtained from different data and/or using different methodology. \citeasnoun{bgim07} applied a bound approach that does not require of exclusion restrictions to study the evolution of wage inequality using the FES data for the period 1978--2000.  They assumed positive sorting for men and women in some of their estimates  to make the bounds more informative.  Interestingly, they mentioned the possibility that the assumption is violated for married women due to assortative matching in the marriage market.\footnote{In results not reported, we find that the negative sorting for married women is robust to the definition of the out-of-work benefit income variable. Thus, we find similar estimates using the income of a tax unit if all the individuals within the tax unit were out of work as the excluded covariate, as in \citeasnoun{bgim07}.} They also found evidence against the validity of out-of-work benefit income as a valid excluded covariate for men.  \citeasnoun{ab17} using the same data from the FES, also found positive sorting for men, stronger  for single than for married men, using an alternative methodology that combines quantile regression for the marginal distributions with a parametric model for the copula. Contrary to our findings, they also found positive selection for women, which is statistically significant only for married women.  \citeasnoun{mr08} estimated a HSM using data from the US-CPS for the periods 1975-1979 and 1995-1999. They found that the selection sorting for women shifted from negative to positive between the two periods. We also find for the UK that the sorting for most women has a  positive trend over time, but remains negative even in 2013 for most of the distribution.  \citeasnoun{mw18} applied the methodology of \citeasnoun{ab17} to data from the US-CPS for the period 1976--2014. They also found negative sorting for women at the beginning of the sample period that became positive during the 90s, and positive sorting for men throughout the entire period. \citeasnoun{bertrand17} pointed out multiple possible explanations for the glass ceiling based on the field of education, psychological attributes or preferences for job flexibility that are compatible with our finding on the importance of sorting on unobservables at the top of the distribution. None of the previous papers distinguished between the selection sorting and selection structure effects.

One limitation of our dataset is that it does not contain a direct measure of work experience. As a final robustness check, we find that the results are not sensitive to the exclusion of college graduates from the sample by redoing the analysis excluding all the individuals who cease school after age 18. This is a relevant exclusion because work experience is a more relevant determinant of wage for highly educated workers.\footnote{These results are available from the authors upon request.}

\section{Monte Carlo Simulation}

We conduct a Monte Carlo simulation calibrated to the empirical application to study the properties of the estimation and inference methods in small samples. The data generating process is the HSM of Example \ref{ex:hsm} with the values of the covariates and parameters calibrated to the data for women in the last ten years of the sample (2004--2013). We do not use the entire dataset to speed up computation. We generate 500 artificial datasets and estimate the DR-model with the same specifications for the selection and outcome equations as in the empirical application and specification 1 for the selection sorting function, i.e. $\rho(x'\delta(y)) = \rho(y)$.

Figures \ref{figure:Sim_col}, \ref{figure:Sim_married} and \ref{figure:Sim_rho} report the biases, standard deviations and root mean square errors for the estimators of the coefficients of the college (age when ceasing school 21--22) and marital status indicators in the outcome equation, and  $\rho(y)$ in the selection sorting function, as a function of the quantile indexes of the values of log real hourly wage in the data used in the calibration.\footnote{We find similar results for the other coefficients of the outcome equation. We do not report these reports for the sake of brevity.} Although these coefficients are constant in the HSM, we do not impose this condition in the estimation. The estimates are obtained with Algorithm \ref{alg:tsdr} replacing $\mathcal{Y}$ by a finite grid containing the sample quantiles of log real hourly wage with indexes $\{0.10, 0.11, \ldots, 0.90 \}$ in the original subsample of women in the last ten years of the sample.  All the results are in percentage of the true value of the parameter. As predicted by the asymptotic theory, the biases are all small relative to the standard deviations and root mean squared errors. The estimation error increases for all the coefficients as we move away from the median towards tail values of the outcome.

\begin{figure}[htbp]
    \includegraphics[width=0.32\textwidth,page=1]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=2]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=3]{Sim_f_113104_coef.pdf}
    \caption{Bias, SD and RMSE for the coefficient of the college indicator in the outcome equation}
    \label{figure:Sim_col}
\end{figure}

\begin{figure}[htbp]
    \includegraphics[width=0.32\textwidth,page=5]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=6]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=7]{Sim_f_113104_coef.pdf}
    \caption{Bias, SD and RMSE for the coefficient of the marital status indicator in the outcome equation}
        \label{figure:Sim_married}
\end{figure}

\begin{figure}[htbp]
    \includegraphics[width=0.32\textwidth,page=9]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=10]{Sim_f_113104_coef.pdf}
    \includegraphics[width=0.32\textwidth,page=11]{Sim_f_113104_coef.pdf}
    \caption{Bias, SD and RMSE for coefficient $\rho(y)$ in the selection sorting equation}
    \label{figure:Sim_rho}
\end{figure}

Table \ref{table:simcb} shows results on the finite sample properties of 95\% confidence bands for the coefficients of the indicators of college and marital status in the outcome equation and  $\rho(y)$ of the selection sorting function. The confidence bands are constructed by Algorithm \ref{alg:mb} with $B=200$ bootstrap repetitions and the same grid of values $\bar{\mathcal{Y}}$ as for the estimators.  We report the average length of the confidence bands integrated over threshold values, the average value of the estimated critical values, and the empirical coverages of the confidence bands. For comparison, we also report the coverage of pointwise confidence bands using the normal distribution, i.e. with critical value equal to 1.96. The last row computes the ratio of the standard error averaged across simulations to the simulation standard deviation, integrated over threshold values. We find that the bands have coverages close to the nominal level. As expected, pointwise bands severely undercover the entire functions. The standard errors based on the asymptotic distribution provide a fair approximation to the sampling variability of the estimator.

\begin{table}[h]
    \caption{Properties of 95\% Confidence Bands}
    \centering
    \label{table:simcb}
    \begin{tabular}{lccc}
        \hline\hline
        & College & Married & $\rho(y)$ \\
        \hline
        Average Length & 0.38 & 0.16 & 0.35 \\
        Average Critical Value & 2.91 & 2.89 & 2.88 \\
        Coverage uniform band (\%) & 96 & 98 & 96 \\
        Coverage pointwise band (\%) & 68 & 64 & 67 \\
        Average SE/SD & 1.04 & 1.05 & 1.07 \\
        \hline\hline
        \multicolumn{4}{l}{\begin{minipage}{9cm}\footnotesize Notes: Nominal level of critical values is 95\%. 500 simulations with 200 bootstrap draws.\end{minipage}}
    \end{tabular}
\end{table}

\section{Conclusion}

 We develop a distribution regression model with sample selection that accommodates rich patterns of heterogeneity in the effects of covariates on outcomes and selection.  The model is semiparametric in nature, as it has function-valued parameters, and is able to considerably generalize the classical selection model of \citeasnoun{heckman74}. Furthermore, the model allows for richer covariate effects than the previous semiparametric generalizations which allowed the location effects for covariates.  We propose to  estimate the model by a process of probit regressions, indexed by threshold-dependent parameters. We show that the resulting estimators of the function-valued parameters are approximately Gaussian and concentrate in a $1/\sqrt{n}$ neigborhood of the true values.  We present an extensive wage
decomposition analysis for the U.K. using new data, generating both new findings and demonstrating the power of the method.

\bibliographystyle{econometrica}