Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
86,929 characters · 18 sections · 14 citation commands
Distribution Regression with Sample Selection and UK Wage Decomposition
\pagenumbering{arabic}
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 gronau74,heckman74. Sample selection biases the estimation of causal and predictive effects when the reasons for not observing the data are related to the outcome of interest. In the wage example, there is sample selection bias whenever the employment status and offered wage depend on common unobserved variables such as ability, motivation or skills. The most widely used solution to the sample selection bias is the Heckman selection model (HSM) introduced in \citeasnoun{heckman74}. The classical HSM provides a parsimonious and convenient way to account for sample selection by assuming parametric Gaussian distributions for the outcome and selection processes and imposing strong homogeneity assumptions on the impact of exogenous covariates. In this paper, we present a generalization of the HSM that allows for expressly non-Gaussian structures and eliminates the strong homogeneity restrictions, resulting in a semi-parametric model with some key parameters represented as nonparametric functions. This generalization enables a tractable econometric approach, similar in ease to the classical HSM.
To illustrate the central concepts, we take labor supply as an example. Let $Y^*$ denote the latent offered wage (outcome), and $D^*$ be the net disutility from working (selection), e.g., difference between reservation and offered wage. We observe the employment status indicator equal to one when the offered wage is greater than the reservation wage, $D = 1(D^* \leqslant 0)$, and wage, $Y= Y^*$, only in the case of employment, $D=1$. Assume that the drivers of the outcome and selection obey the additive structure:\footnote{We can have other exogenous covariates $X$ affecting all model components. We interpret the discussion here as conditional on $X$ having taken a fixed value. We suppress the dependence on $X$ here for notational convenience.}
where $(\mu, \nu(Z))$ are the means of the outcome and selection, $Z$ is an instrumental variable that shifts the disutility of working but does not affect the offered wage, and $(U,V)$ are centered stochastic shocks independent of $Z$. The classical HSM results from restricting $(U,V)$ to follow a Gaussian distribution:
where $\Phi_2(\cdot,\cdot;\rho)$ denotes the standard bivariate normal distribution function with correlation $\rho$, and $\sigma_U$ and $\sigma_V$ the standard deviations of $U$ and $V$. This model allows one to identify the distribution of the offered wage from the distribution of observed wage and provides tractable inference. The classical HSM has been challenged for its reliance on parametric assumptions to obtain identification, and the data often reject these assumptions (e.g., due to clustering of offered and observed wages at the minimum wage and other levels). These challenges motivate the generalizations we consider in this study.
We generalize the classical HSM in several dimensions. We start by relaxing the Gaussian assumption on the marginal distributions of the stochastic shocks in ((ref)):
where $F_U$ and $F_V$ are non-parametric marginal distributions of the stochastic shocks $U$ and $V$, which can be continuous, discrete or mixed, and $\Phi^{-1}$ is the quantile function of the standard normal. The distribution of the latent offered wage $Y^*$ is characterized by the parameters $(\mu, F_U)$. Despite the still partly Gaussian nature, the model can astutely generate drastically non-Gaussian distributional shapes, as demonstrated by \citeasnoun{wasserman-nonparanormal} and illustrated in Figure 1. Our contribution is to demonstrate that the distribution of the latent offered wage $Y^*$ is still identified from the observed distribution of wages and employment if the instrumental variable $Z$ takes on at least two values, and to provide tractable estimation and inference methods. Thus, by using ((ref)) we break away from parametric assumptions and Gaussian behaviors of the classical sample selection model but maintain the same level of tractability.
\paragraph{The generalized HSM} The semi-parametric model ((ref)) presents an expressive generalization of the classical model, but for us it merely serves as a starting point to fix ideas. We can get even more flexibility and expressivity while retaining the tractability of the identification and inference. We further relax (ref)-(ref) by allowing the effect of $Z$ on $D^*$ to be non-additive and the dependence between $U$ and $V$ to be heterogeneous. The general form of our model specifies the joint distribution of latent outcome and selection drivers as
where the entire conditional distribution of $D^*$ can now depend on $Z$; and $\rho(y,d)$ is a local correlation parameter that measures the strength of selection or sorting. The model imposes exclusion restrictions -- $Z$ affects the marginal distribution of $D^*$, but does not affect the marginal distribution of $Y^*$ nor the strength of selection. We show that this model admits identification and inference as convenient as the classical model, requiring again only that $Z$ takes on at least two values. Moreover, the model is overidentified (hence testable) when $Z$ takes on more than two values.
Finally, we incorporate exogenous covariates yielding the model:
where the covariates $X$ can affect both the marginal conditional distributions of $Y^*$ and $D^*$ and the sorting mechanism. This model expressly allows for heterogenous effects of $X$'s on different parts of distribution. The data on U.K. wages lends strong empirical support for this property, with marital status and other variables heterogenously impacting different parts of the wage distribution. The heterogeneity in $X$ property also contrasts sharply with the classical model ((ref)) or its generalization ((ref)), which confine the covariates to shift the location of the distribution, but not other shape properties; see \citeasnoun{mr08} for a use of the Gaussian model where covariates shift both location and scale, but not other properties of the distribution. The data firmly reject the location shift restrictions.
Our econometric development proceeds by specifying a flexible pointwise approximation to ${\mathrm{P}}(Y^* \leqslant y, D^* \leqslant d \mid X, Z)$ based upon the model's logic, akin to those used in univariate distribution regression (DR) models foresiperacchi95,chernozhukov+13inference.\footnote{Flexible refers to the use of regressors constructed from original raw regressors by taking series transformations and interactions, for example.} This leads to a tractable Heckman-type approach to estimating the model parameter functions. The first step consists of a probit regression for the selection equation, as in the Heckman two-step method heckman79. The second step estimates multiple bivariate probit regression with sample selection correction. We derive functional central limit theorems for all the estimators and their bootstrapped versions and use these results to perform uniform inference on function-valued parameters, for example, the distribution of the offered wages as well as counterfactual distributions induced by various hypotheses. This type of inference provides simultaneous confidence bands and hypotheses tests about the functions of interest. Specifically, we rely on multiplier bootstrap gz84, applied to the estimated influence functions as in \citeasnoun{lewbel95consistent}, \citeasnoun{hansen1996inference}, \citeasnoun{ch06}, \citeasnoun{ks12} and \citeasnoun{cck13}.
We utilize the proposed methods to analyze data on wages and employment in the UK, covering the period from 1978 to 2013. By estimating the conditional wage distributions for men and women and conducting wage decompositions that account for the endogeneity of employment selection, we shed light on the relationship between wage and employment in the UK. The results highlight the existence of positive sorting among single men and negative sorting among married women, which is consistent with the assortative matching hypothesis in the marriage market. Furthermore, the findings show that this difference in selection sorting explains a substantial portion of the gender wage gap at the upper end of the distribution, aligning with recent theories centered on glass ceiling. We still find that most of the gender gap in offered and observed wages can be attributed to differences in the wage structure that are often associated with gender discrimination in the labor market. The effect of education on wages is positive and increases with the distribution. In summary, our results demonstrate the significance of considering selection and flexible heterogeneity in empirical analysis and provide support for the use of the generalized Heckman selection model ((ref)). Our results also align with those of \citeasnoun{mw18} which applied a quantile-regression based methodology of \citeasnoun{ab17} to U.S. data over a similar period.
This paper contributes to the extensive literature on sample selection in economics and statistics. The classical references in this field include works by \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 commonly utilized approach to sample selection is the Heckman selection model developed by Heckman heckman74,heckman76,heckman79,heckman90. Several extensions to the HSM have been made, including parametric extensions by \citeasnoun{lee83}, \citeasnoun{prieger02} and \citeasnoun{smith03} using different distributions, and a bivariate t-distribution extension by \citeasnoun{marchenko12} for heavy-tailed data. Semi-parametric versions have been developed by \citeasnoun{ahn93}, \citeasnoun{powell94}, \citeasnoun{andrews98}, and \citeasnoun{newey99}, while \citeasnoun{das03} introduced a nonparametric version with additive shocks, all of which focus on location effect models with homogeneous effects. None of these extensions can accommodate all sources of heterogeneity considered in our model. Previous studies on partial identification of unrestricted sample selection models include works by \citeasnoun{manski1990}, \citeasnoun{balke1994counterfactual}, \citeasnoun{manski94}, \citeasnoun{heckman2001instrumental}, and \citeasnoun{manski03}, among others. Our work contributes with new partial identification results to this development, with the main emphasis on achieving point identification. In the online Supplemental material, we consider relaxed forms of the exclusion restriction on the sorting that obtain a nested sequence of bounds, which starts from a single point and ends with the (agnostic) \citeasnoun{balke1994counterfactual} bounds. This sequence provides a form of sensitivity analysis, analogous in spirit to \citeasnoun{christensen:sense}'s sensitivity analysis for structural models.
\citeasnoun[AB17]{ab17} proposed another extension of the HSM, which like our model accounts for multiple sources of heterogeneity. Their approach is based on quantile regression, which models the marginal distribution of the latent outcome, and a parametric copula model that links the latent selection and outcome variables. However, compared to our approach, this model requires the latent variables $(Y^*,D^*)$ to be continuous, while our approach accommodates mixed discrete-continuous distributions, which is important given that offered wages are often constrained to be above a minimum wage level, and observed wages have point masses at minimum wage and other levels due to wage rounding. Additionally, AB17 models the effect of covariates on the conditional quantile of the latent distribution, while we model the effect on the conditional latent distribution. The identification assumptions of the two models are also distinct, as our approach imposes more structure on the dependence between the outcome and selection processes, whereas AB17 requires continuous variation in the instrument $Z$ and real analyticity of the copula function governing the dependence of the shocks in the outcome and selection structures. The analyticity condition implies extrapolability of the conditional distribution of outcome to the case where there is no selection, enabling identification at infinity. \ \ In this case, however, identification does not automatically imply consistent estimability, because the continuation of an analytical function is known to be an ill-posed problem. \ \ To overcome the latter problem, regularization must be used: AB17 impose parametric assumptions on the copula as a means of such regularization to carry-out estimation.
Finally, there are relevant connections to the work of \citeasnoun{wasserman-nonparanormal} on the estimation of graphical models using the semi-parametric family ((ref)). While our approach considers the more general model ((ref)), we contribute by offering a DR based approach, which can incorporate covariates $(X,Z)$ in a general yet tractable manner. Furthermore, our focus is on addressing the sample selection problem in the context of these models.
\paragraph{Outline} Section (ref) examines the identification problem under sample selection using a new representation of a joint distribution. Section (ref) 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) reports the results of the empirical application. The proofs of the results of Section (ref) are gathered in Appendix (ref). The online Supplemental Material (SM) contains deferred discussions of Sections (ref)--(ref), the asymptotic theory for our estimation and inference methods, additional empirical results, a Monte Carlo simulation calibrated to the empirical application and other technical results.
Our first result shows that any joint distribution of two random variables has a local Gaussian representation (LGR). This is a useful representation for econometric analysis, because it naturally nests the joint normal distribution and rich semi-parametric models that allow for nonparametric marginal distributions-- with much more room to spare. Indeed, we shall use the LGR 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.
Lemma (ref) establishes that any bivariate CDF admits a unique pointwise representation by standard bivariate normal distributions. This result is related to Sklar's copula representation of the joint distributions but is different, especially in using the localized correlation. Also 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}. We also note that Lemma (ref) easily extends to CDFs conditional on covariates by making all the parameters dependent on the value of the covariates.
In what follows, it is convenient to define the LGR parameters: $$\mu(y) := \Phi^{-1}(F_{Y^*}(y)), \ \ \nu(d) := \Phi^{-1}(F_{D^*}(d)),$$ noting that $\mu(y) \in \overline{\mathbb{R}}$ and $\nu(d) \in \overline{\mathbb{R}}$, where $\overline{\mathbb{R}} := \mathbb{R} \cup \{-\infty,+\infty\}$ is the extended real number line. By LGR, the mapping between $F_{Y^*,D^*}(y)$ and its LGR parameters $(\mu(y), \nu(d), \rho(d,y))$ is bijective. Thus LGR carries the same information as the joint CDF.
The implied correlation $\rho(y,d)$ is a key parameter that measures local dependence. Indeed, $\rho(y,d) = 0$ if and only if the distribution $F_{Y^*,D^*}$ factorizes at $(y,d)$ because $$ 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, $\rho(y,d) = 0$ if and only if the events $\{Y^*\leqslant y\}$ and $\{D^*\leqslant d\}$ are independent. Moreover, $\rho(y,d)$ is positive if and only if correlation of $1(Y^*\leqslant y)$ and $1(D^* \leqslant d)$ is positive. Thus, if $\rho(y,d)$ is positive everywhere then $Y^*$ and $D^*$ are positively quadrant dependent lehmann66. We prove these assertions and provide additional discussion in Appendix (ref) of the SM.\footnote{See also \citeasnoun{survey18} for a recent review of copula-based measures of local dependence.}
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
i.e., $D$ is an indicator for $D^* \leqslant 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.
We can write the distribution of the observed variables as $$ {\mathrm{P}}(D = 1) = \Phi(\nu) \text{ and } {\mathrm{P}}(Y \leqslant y, D=1) = \Phi_2(\mu(y), \nu; \rho(y)), $$ where $\mu(y)$, $\nu := \nu(0)$ and $\rho(y) := \rho(y,0)$ are the parameters of LGR for the latent distribution $F_{Y^*,D^*}$. The identified set for these parameters determine the identified set for $F_{Y^*}(y)$. Note that in particular, $$ F_{Y^*}(y) = \Phi(\mu(y)). $$ As shown below, $\mu(y)$ and $\rho(y)$ are partially identified. We proceed by characterizing the identified set for these parameters and provide exclusion restrictions to achieve point identification. We also examine the role of relaxed exclusion and other restrictions in reducing the size of the identified set in Appendix (ref) of the SM.
To understand the source of partial identification under sample selection, note that there are two free probabilities, $ {\mathrm{P}}(D = 1)$ and ${\mathrm{P}}(Y \leqslant y, D = 1)$, to identify three parameters, $\mu(y)$, $\nu$ and $\rho(y)$. The selection probability pins down $\nu$: $ \nu = \Phi^{-1}({\mathrm{P}}(D = 1)). $ The parameters $\mu(y)$ and $\rho(y)$ are partially identified by the set of solutions in $(\mu,\rho)$ to the equation $$ {\mathrm{P}}(Y \leqslant y, D = 1) = \Phi_2(\mu,\Phi^{-1}({\mathrm{P}}(D = 1)); \rho). $$ These solutions form a one-dimensional manifold in $\mathbb{R} \times (-1,1)$ spivak65,munkres91.\footnote{This is 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 $.}
To state the exclusion restrictions, let $Z$ be a candidate instrumental 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:
with parameters $\mu(y \mid z) \in \overline{\mathbb{R}}$, $\nu(d \mid z) \in \overline{\mathbb{R}}$, and $\rho(y,d \mid z) \in [-1,1]$. The exclusion restrictions are the following.
The condition that $Z$ is binary emphasizes that our identification strategy does not rely on large variation of $Z$. If $Z$ is not binary, we only require that Assumption (ref) be satisfied for two values of $Z$. If it is satisfied for more than two values of $Z$, then the model is overidentified and the exclusion restrictions become testable.\footnote{We leave the development of such specification test to future research.} Non-degeneracy states that there is sample selection and that $Z$ has variation in the selected population. Relevance requires that $Z$ affects the probability of selection. The condition ${\mathrm{P}}(D = 1 \mid Z = 1) < 1$ precludes identification at infinity, which we analyze separately in Remark (ref). The sign of the first inequality can be reversed by relabelling the values of $Z$. Outcome exclusion is a standard exclusion restriction, which is not sufficient for point identification in the presence of sample selection balke1994counterfactual,manski94,heckman2001instrumental,manski03. It holds when $Y^*$ is independent of $Z$.\footnote{\citeasnoun{kitagawa2010} developed a test for the outcome exclusion.} Selection sorting exclusion requires the local correlation function to be independent of $Z$.\footnote{ \citeasnoun{torgovitsky2010identification} previously used this type of restriction to analyze identification of nonseparable models with continuous endogenous explanatory variables.} AB17 consider an alternative condition to selection sorting exclusion based on real analyticity and continuous variation of $Z$. We refer to Appendices (ref) and (ref) in the SM for a comparison with the analyticity approach and another alternative approach based on imposing semi-parametric structure on the sorting mechanism when $Z$ takes on more than 2 values.
We can get some intuition about the outcome and selection exclusion restrictions with the parametric and semi-parametric models of labor supply given in (ref)--(ref). These models trivially satisfy the conditions stated above, because $Z$ only affects the latent disutility $D^*$ through the location, and the local correlation parameter is constant and does not depend on $Z$ (or $y$ for that matter). Yet the selection sorting exclusion does entail some loss of generality, which motivates forms of relaxed conditions on selection sorting that we consider in Appendix (ref) of the SM.
We now show how the presence of exclusion restrictions helps identify the parameters. Under Assumption (ref) the conditional LGR at $d=0$ simplifies to
where we simplify the notation $\nu(z) := \nu(0 \mid z) \text{ and } \rho(y) := \rho(y,0).$ We can relate this representation to the conditional distribution of observed variables as
As before, $\nu(z)$ is identified from the conditional selection probability:
Moreover, $\mu(y)$ and $\rho(y)$ are identified as the solution in $(\mu, \rho)$ to
This is a nonlinear system of two equations in two unknowns. The result below states that the solution exists and is unique.
The result follows by showing that the Jacobian of the equations in (ref) is a P-matrix for all $\mu \in \mathbb{R}$ and $\rho \in (-1,1)$, so uniqueness follows by the global univalence result in Theorem 4 of \citeasnoun{gale65}. We defer to Appendix (ref) the analysis of the case $\rho(y)^2=1$ due to its more technical nature. This case deals with boundary situations that can occur at extreme values of $y$ or other non-typical circumstances. In these situations, we can have either point or partial identification. They are easy to detect empirically.
We provide a further interpretation of the selection exclusion restriction in a formulation of the sample selection problem in terms of the propensity score. The sample selection process is often represented as $D = 1\{V \geqslant \nu(Z)\}$, where $\nu(z) := \Phi^{-1} (p(Z))$, with $p(Z):={\mathrm{P}}(D=1 \mid Z=z)$ being the propensity score, and $V \mid Z \sim N(0,1)$ is the unobserved selection normal score. Here we shall maintain Assumption (ref)(1)--(3). The following theorem shows the condition that selection sorting exclusion imposes on the relationship between $Y^*$ and $V$. Let $$F_{Y^*,V \mid Z}(y,v \mid z)= \Phi_2(\Phi^{-1}(F_{Y^*}(y)), v; \tilde{\rho}(y,v \mid z) )$$ be the LGR of the joint CDF of $(Y^*,V)$ conditional on $Z$, where we use that $Y^*$ is independent of $Z$ and $V \mid Z \sim N(0,1)$.
Theorem (ref) shows that selection exclusion holds whenever the implied correlation $\tilde \rho(y,v \mid z)$ between $Y^*$ and $V$ does not depend on $v$ and $z$. A sufficient condition is that $(Y^*,V)$ are jointly independent of $Z$ and $\tilde \rho(y,v) = \tilde \rho(y)$ for all $v$ in the support of $\nu(Z)$, where $\tilde \rho(y,v)$ is the local dependence parameter of the joint CDF of $(Y^*,V)$. In the context of labor supply, $Y^*$ is latent wage and $V$ is the ranking in the conditional distribution of utility/net benefit of employment. The classical HSM in ((ref)) and its semi-parametric generalization in ((ref)) trivially satisfy this condition, because they have constant local correlation, $\tilde \rho(y, v)= \rho$, that is, the strength of selection does not vary with the level of offered wage or the employment ranking. Relative to such traditional models, the condition $\tilde \rho(y,v) = \tilde \rho(y)$ leaves some room for additionally flexibility -- namely, the strength of selection can vary with the level of offered wage. For example, $y \mapsto \tilde \rho(y)$ increasing means that workers are more likely to select into the labor force if the offered wage's rank is higher (holding the employment ranking fixed).
The selection exclusion $\tilde \rho(y,v) = \tilde \rho(y)$ is essentially equivalent to the single index restriction:
where $y \mapsto a(y)$ and $y \mapsto b(y)$ are nonparametric functions linked to the LGR of $Y^*$ and $V$, $F_{Y^*,V}(y,v)= \Phi_2(\mu(y), v; \tilde{\rho}(y) )$, via $a(y) = \mu(y)/\sqrt{1-\tilde{\rho}(y)^2}$ and $b(y) = - \tilde \rho(y)/\sqrt{1-\tilde{\rho}(y)^2}$. This index form enhances interpretability and also connects to the econometric literature that utilizes the simplifying index restrictions; e.g., \citeasnoun{ichimura1993semiparametric}, \citeasnoun{klein1993efficient}, \citeasnoun{powell94} and \citeasnoun{vytlacil:equiv}. The proof of the equivalence in (ref) is given in Appendix (ref). Homogeneous models, such as the classical Heckman's labor supply model or the semi-parametric generalization ((ref)), impose the much stronger single index condition: $$ {\mathrm{P}}(Y^* \leqslant y \mid V=v) = \Phi \left(a(y) + b v\right), $$ for all $y$, for $b = -\rho/\sqrt{(1- \rho^2)}$ being a rescaled correlation coefficient; note that the Gaussian restriction on $Y^*$ of the HSM corresponds to the further linearity restriction: $a(y) = a y$. This connection highlights the generalization our model brings: in the general model $a(y) \neq a y$ allows for expressive departures from Gaussianity, and $b(y) \neq b$ allows for much richer patterns of dependence across the conditional distribution. The restriction (ref) is related to a linear model assumption on the expectation of $Y^*$ conditional on $V$ that was previously used in \citeasnoun{brinch2017beyond} and \citeasnoun{kowalski2016doing} to identify marginal treatment effects in treatment effects settings with endogenous treatment assignments and discrete instruments. One difference, however, is that our restriction is a single index restriction on the conditional distribution.\footnote{Single index restrictions on the conditional distribution like (ref) and linearity restrictions on the conditional expectation are not nested. This can be seen, for example, from the relationship $$ {\mathrm{E}}[Y^* \mid V = v] = \int_{-\infty}^{\infty} [1(y > 0) - {\mathrm{P}}(Y^* \leqslant y \mid V=v)] \mathrm{d} y. $$} Moreover, \citeasnoun{brinch2017beyond} and \citeasnoun{kowalski2016doing} analysis imposes that $(Y^*,V)$ are jointly independent of $Z$. In Appendix (ref), we provide an alternative sufficient condition for the single index restriction that does not rely on this joint independence.
We now elaborate on the value that the variation $\rho(y)$ with respect to $y$ brings to the table. What does this property mean in terms of the flexibility of dependence patterns between net offered wage $Y^*$ and employment ranking $V$? The standard models that have homogeneous $\rho(y) = \rho$ or equivalently homogeneous $b(y) = b$, are only able to generate limited forms of dependence. The left panel of Figure (ref) presents an example of the contour of the joint pdf of $(Y^*,V)$ when $\rho(y) = .14$. In contrast, our model is able to generate richer forms of dependence between $(Y^*, V)$ as illustrated in the middle and right panels. The middle panel has the same correlation between $Y^*$ and $V$ as the left panel, but we see that the dependence increases from small values as we move towards the upper-right tail corner. The right panel is another example, where the correlation between $Y^*$ and $V$ is zero, but $(Y^*, V)$ are positively dependent when $Y^*>0$ and negative dependent when $Y^*<0$. In all examples, the marginal distributions are fixed to be normal, but the joint distributions are not normal, except for the left panel.
The above discussions aim to support the idea of using sorting exclusion, but one may and should challenge it in applications. To this end, we emphasize three points: First, as we mentioned in the discussion of Assumption (ref), selection exclusion is testable when $Z$ takes on more than 2 values. Second, one can work with the relaxed forms of exclusion restrictions, namely the sign restrictions and r-relaxed exclusion restrictions and compute the bounds provided by Theorem (ref) in the SM. Third, if richer instruments are available, one can still obtain results for identification under analyticity or other regularity assumptions on the local correlation function without selection sorting exclusion; see Appendix (ref) in the SM.
We consider a semi-parametric version of the LGR with covariates:
where $Y^*$ is the latent outcome of interest, which can be either continuous, discrete or mixed; $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.\footnote{It is understood that the models are flexible in the following sense. Given $x$ and $z$, we can generate constructed regressors $t(x)$ and $b(z)$ as technical transformations of $x$ and $z$, for example, by taking powers or splines of the components and their interactions. We then can reassign notation $x \leftarrow t(x)$ and $z \leftarrow b(z)$. This convention entails no loss of generality provided that $t(x)$ and $b(z)$ contain the same information as $x$ and $z$. In this paper, we do not formally consider the case where the dimensions of $t$ and $b$ grow with the sample size, but this is possible along the lines of the series literature, with the inference being operationally equivalent to the case with the fixed number of series terms, provided that the number of series terms is small compared to the sample size and that the approximation errors are negligible; e.g., \citeasnoun{newey:series} and \citeasnoun{chen:Chapter}.} The excluded covariates avoid reliance on functional form assumptions to achieve identification.
The model (ref) is a semi-parametric version of (ref), where we replace the components $\Phi^{-1}(F_{Y^*}(y\mid X))$, $\Phi^{-1}(F_{D^*}(d \mid Z,X))$, and $\rho(y,d \mid X)$ by three indexes.\footnote{Here we include the covariates $X$ in $Z$ to lighten the notation and denote the covariates that satisfy the exclusion restriction by $Z_1$ instead of $Z$.} 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 (ref) are included to take into account that the selection is defined by $D^* > 0$ instead of $D^* \leqslant 0$. We use this definition to facilitate the interpretation of the parameters and the comparison with the classical HSM; see Example (ref).} In the empirical application that we consider below, $Y^*$ is offered wage, $D^*$ is the net utility from working, e.g., 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).
The model (ref) is semi-parametric 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 equations 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.\footnote{The parametric copula model of AB17 imposes that the sign of the sorting is the same across the latent wage distribution.} The function $u \mapsto \rho(u)$ is a known link with range $[-1,1]$, e.g. the Fisher transformation 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)$
The parameters $\beta(y)$ and $\pi$ have the same interpretation as in probit models. Their signs are informative about the signs of the partial effects of the corresponding covariates on the conditional distribution of the outcome or probability of selection, and ratios of their components yield ratios of partial effects. Indeed, if $X$ is continuous, $$ \frac{\partial F_{Y^*}(y\mid Z=z)}{\partial x} = - \beta(y) \phi(-x^\prime \beta(y)), $$ where $\phi$ is the standard normal PDF. The parameter $\delta(y)$ determines the sign of the effect of the covariates on the sorting because $$ \frac{\partial \rho(x'\delta(y))}{\partial x} = \delta(y) \dot{\rho}(x'\delta(y)), $$ and $\dot{\rho}(u) = \partial \rho(u)/ \partial u = 1 - \tanh(u)^2 > 0$. Appendix (ref) in the SM provides additional discussion on the interpretation of the model parameters.
Several key functionals of the model's parameters (ref) can be of interest. One is the marginal distribution of the latent outcome $Y^*$ $$F_{Y^*}(y) = F_{Y^*}(y; \beta, F_X) := \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. Using the formula above, we can also construct counterfactual distributions by replacing $\beta(y)$ and $F_X$ by coefficients and distributions from different populations or groups, $\bar \beta(y)$ and $\bar F_X$. These distributions are useful to decompose the distribution of offered wages between females and males or between blacks and whites, which can be the basis to uncover discrimination in the labor market. Another functional is the probability of selection $$ {\mathrm{P}}(D=1) = {\mathrm{P}}(D = 1; \pi, F_Z) := \int {\mathrm{P}}(D=1 \mid Z=z) dF_Z(z) = \int \Phi(z^\prime \pi)dF_Z(z), $$ which we can also use to define counterfactuals and employ them 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
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.
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) \geqslant \tau\}$ is the quantile or left-inverse operator.
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 \leqslant y) \mid X = x]. $$ To implement this idea, we construct the set of indicators for the selected observations $$ I_{yi} = 1(Y_i\leqslant 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
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 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:
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 interest are constructed from the estimators of the parameters using the plug-in method. For example, the estimator of the distribution of the latent outcome and the estimator of the probability of selection are
and the estimators of the counterfactual distributions of the observed outcome are constructed from
by choosing the estimators of $ \widehat{\beta}(y)$, $\widehat{\pi}$, and $\widehat{\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. See Section (ref) of SM for details.
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 being 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 $$ {\mathrm{P}}\left[c'\theta_y \in CB_p(c'\theta_y), \text{ for all } y \in\mathcal{Y}\right] \to p. $$ We form $CB_p(c'\theta_y)$ as $CB_p(c'\theta_y):=c'\widehat{\theta}_y\pm cv(p)SE(c'\widehat{\theta}_y),$ where $\widehat{\theta}_y$ is the estimator of $\theta_y$ defined in Algorithm (ref), $SE(c'\widehat{\theta}_y)$ is the standard error of $c'\widehat{\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'\widehat{\theta}_y-c'\theta_y|}{SE(c'\widehat{\theta}_y)}.$$
We obtain the standard error and critical value from the limit distribution of the stochastic process $y \mapsto \widehat{\theta}_y$ derived in Section (ref) of the SM. 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). The multiplier bootstrap is implemented using the following 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.
We apply the DR model with sample selection to carry out wage decompositions accounting for endogenous employment participation using data from the United Kingdom.
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 labor 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,396 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 ($\leqslant$15, 16, 17--18, 19--20, 21--22 and $\geqslant$ 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.} We provide descriptive statistics of the variables and some background on the U.K. labor market using our data in Section (ref) of the SM.
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 b09. It is designed to calculate the income of a tax unit if the individual was considered 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 an excluded covariate. In our case, the outcome and selection sorting exclusions 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 the Introduction, we stated a rich semi-parametric generalization of Heckman's labor supply model that trivially satisfies the selection sorting exclusion, motivating their use in our analysis.
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:
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 wage decompositions presented below are not sensitive to the specification of the sorting equation.}
We report point estimates and 95% confidence bands for the local correlation function $y \mapsto \rho(x'\delta(y))$ in the selection sorting equation. Estimates and 95% confidence bands for the coefficients of the selection, outcome and selection sorting equations are given in Section (ref) of the SM. The estimates are obtained with Algorithm (ref) 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) with $B=500$ bootstrap repetitions and the same finite grid as for the estimates. We also report estimates from the HSM of Example (ref) with dashed lines as a benchmark of comparison.
Figures (ref)--(ref) display the estimates of the sorting functions for specifications 1--3, respectively. Figure (ref) 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), 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 neal04. Figure (ref) 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, especially at the bottom of the distribution. Figure (ref) in the SM 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 $\delta(y)$ in the SM show that the results on the trends are statistically significant.}
Figure (ref) shows point estimates of the quantiles of offered and observed wages for men and women based on specification 4. 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), the distributions of both types of wages can be expressed as functionals of the model parameters, and estimated using the plug-in estimators (ref) and (ref).\footnote{The model-based estimator of the observed distribution in (ref) produces almost identical estimates to the empirical distribution of the observed wages.} We find different sample selection biases for men and women. Thus, the quantiles of the observed wages are below the quantiles of latent wages for men, consistent with the positive selection in the sorting equation, whereas they are similar for women. In this case, the bias coming from the negative sorting is almost exactly offset by the difference in the composition between working and non-working women.
Figure (ref) compares the quantile function of offered wages between men and women and carries out a gender wage gap analysis based on specification 4. The estimates and 95% confidence bands for the other specifications are reported in the SM. The gender wage gap analysis is based on the counterfactual distributions $$ F_{Y^* \langle j, k \rangle}(y) = F_{Y^*}(y; \beta_j, F_{X_k}) = \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 (ref) and the quantile and increasing rearrangement operators, see Section (ref) of the SM. 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_X$, have very little explanatory power. This result can be interpreted as evidence of gender discrimination in the labor market.
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: (1) selection (employment) sorting: $\delta(y)$; (2) selection (employment) structure: $\pi$; (3) outcome (wage) structure: $\beta(y)$; and (4) composition: $F_Z$. 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:
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 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) = F_Y(y; \beta_r, \pi_s, \delta_t, F_{Z_k}) = \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) 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 4. The bands for the contributions are joint for all the components and rely on the delta method; see Remark (ref) in the SM. Estimates of the components of the decomposition and the analysis based on specifications 1--3 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. However, differences in sorting and employment structure also account for an important percentage of the gap, especially 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 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.\footnote{While the sign of the employment structure contribution changes with the order of the decomposition, neither its importance nor the significance of the contributions of the other components are sensitive to this order.} 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) 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) 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 dotted lines pick up the average contributions of the components, but miss all the heterogeneity across the distribution.
Figures (ref) and (ref) 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 2 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.
The main findings can be summarized as: (1) 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; (2) heterogeneity in selection sorting decreases gradually over time; (3) differences in returns to characteristics in the wage equation, which might be associated with gender discrimination in the labor market, account for most of the gender wage gap; (4) 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 (5) 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. AB17 using the same data from the FES, also found positive sorting for men, stronger for single than for married men, but they employed 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 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.}
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 semi-parametric in nature, as it has function-valued parameters, and is able to considerably generalize the classical selection model of \citeasnoun{heckman74}. Furthermore, the model accounts for richer covariate effects than the previous semi-parametric generalizations which allowed only the location effects for covariates. We propose to estimate the model by a process of bivariate 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. Our identification approach is especially designed for the sample selection problem but can be applied to other settings. In work in progress, we show how the selection exclusion can be used to identify causal effects in treatment effects models with endogeneity in the presence of an instrumental variable that can be binary.