Distribution Regression with Censored Selection

Distribution Regression with Censored Selection$^*$

\pagenumbering{arabic}

Introduction

Empirical researchers frequently encounter non-random sample selection in their data. A prominent example is the estimation of wage equations when the wages of non-working individuals are not observed. Sample selection bias arises when there are unobservables affecting both sample selection and the outcome of interest, such as employment status and offered wage. While the most widely used method for addressing sample selection bias is the classical Heckman selection model (HSM) introduced by heckman1974shadow, this model relies on parametric Gaussian error structures and assumes homogeneous effects of exogenous covariates. Recently, chernozhukov2023distribution (CFL) developed a distribution regression (DR) model that accounts for endogenous sample selection in an effort to generalize the HSM. They provide a method to analyze the conditional distribution, extending beyond the mean or median, while taking sample selection into account. Furthermore, they depart from Gaussian error structures and allow for heterogeneous effects of covariates on both the outcome distribution and the sample selection process.

This paper develops a distribution regression model with a censored selection rule, building on the distribution regression model with binary selection rule proposed by CFL. Implementing a censored selection rule requires more detailed information about the sample selection mechanism compared to using a binary selection rule. In the context of wage equations, this means modeling the employment selection equation with a censored work hours equation, instead of a binary employment status equation. While employing a censored selection rule does not aid in identification - since our model relies on the same identification assumptions as CFL - the main advantage is that it allows for the analysis of richer patterns of selection heterogeneity. First, this approach enables modeling the distribution of the selection variable, such as the distribution of work hours, which is not possible in sample selection models that rely on a binary selection rule. Second, it allows for the estimation of selection sorting (e.g., sign of selection) according to the level of the selection variable. For instance, while binary sample selection models enable estimating positive or negative selection into employment, our approach can examine selection sorting according to work hours, such as selection into full-time or overtime work.

An example where heterogeneous selection patterns, which can be revealed by a censored selection model, are of interest, is the study of the gender wage gap. blau2017gender identify differences in work hours and selection sorting across genders as part of contributing factors to this gap. In terms of work hours, they highlight a significant gender disparity, noting that women comprise a larger share of the part-time workforce. This higher prevalence of part-time work among women compared to men may widen the overall gender wage gap, as part-time workers typically earn lower hourly wages than full-time workers (blank1990parttime, hirsch2005part). On the other hand, there is literature examining how long work hours affect the gender wage gap among highly skilled workers (noonan2005pay, bertrand2010dynamics, goldin2014grand). These studies show that in business and law professions, the return to working long, inflexible hours is substantial, and this mechanism contributes more to the gender wage gap than underlying skill differences. Employing a censored selection rule and estimating the conditional distribution of work hours can help better understand the characteristics of workers with both short and long hours.

There is also a vast literature documenting different selection sorting patterns into employment by gender and their impacts on the gender wage gap. For example, mulligan2008selection find weaker evidence of convergence in the gender wage gap after accounting for sample selection. Similarly, arellano2017quantile, maasoumi2019gender, and chernozhukov2023distribution address sample selection bias in the wage distribution and highlight that examining wage inequality among the employed may present a distorted picture of overall labor market disparities. As there are significant gender differences in work hours, studying selection sorting patterns by work hours or worker types could offer valuable insights into the mechanisms behind the gender wage gap.

Given the advantages of a censored selection model, we highlight that implementing the model is often possible using available data to implement a binary selection model. As noted in fernandez2021nonseparable, the binary selection model is more widely adopted in empirical analysis compared to the censored selection model, which may be due to the widespread popularity of the HSM. \ifnum\value{page}=1 \origfootnote{It should be noted, however, that the HSM was initially developed for the case of censored selection in the context of estimation of wage equations when work hours are observed.} \else {\arabic{footnote}} \origfootnote{It should be noted, however, that the HSM was initially developed for the case of censored selection in the context of estimation of wage equations when work hours are observed.} \fi However, researchers often dichotomize censored selection variables to employ the binary selection model. Examples include the duration of unemployment, the length of a training program, and the magnitude or duration of welfare benefit receipt. Therefore, for researchers interested in analyzing the distribution of selection variables and heterogeneous selection sorting patterns, implementing a censored selection model is often feasible with existing data.

We analyze identification, estimation, and inference methods in a distribution regression model that corrects for sample selection bias using a censored selection rule. The parameters of interest are features of the joint distribution of the latent outcome and selection variable. A key distinction from CFL is that the observed selection variable is discrete or continuous, rather than binary. As a result, the local correlation parameter of the joint distribution, which measures the strength of selection or sorting, depends not only on the outcome but also on the selection variable. For identification, we rely on exclusion restrictions. In particular, we require the existence of a binary instrumental variable that affects the marginal distribution of the selection variable but does not affect the marginal distribution of latent outcome or the relationship (sorting) between the latent outcome and selection variables. The relevance condition for the marginal distribution and exclusion restriction for the sorting are local in that they only need to hold at a single level of the selection variable such as the censoring point. We propose a tractable estimation algorithm consisting of three steps involving (bivariate) probit regressions with sample selection corrections. Lastly, we provide functional central limit theorems for our estimators and uniform inference methods based on multiplier bootstrap.

We apply our method to wage and work hours data in the UK. The goal is to study selection sorting patterns across work hours (e.g., selection into full-time and overtime work) and conduct wage decompositions between gender by worker types (work hours). First, we find heterogeneous selection patterns by gender, marital status, year, and worker types. Specifically, single men exhibit positive selection into full-time work across all wage quantiles, while married women do so only at lower wage quantiles. Since married men and single women do not show significant selection effects, we interpret this as indicating that single men and married women have greater flexibility in selecting their work hours. Both genders display negative selection into overtime work, with the effect being stronger at the top of the wage distribution.

Second, we decompose the difference between the male and female wage distributions into composition, wage structure, hours structure, and sorting effects for full-time and over-time workers. This is a distributional generalization of the Oaxaca-Blinder decomposition (kitagawa1955components; oaxaca1973male; blinder1973wage) which accounts for sample selection. The decomposition for full-time workers reveals that aggregated hours and sorting effects contribute to a smaller observed gender wage gap at the bottom and a larger gap at the top of the distribution. Also, while the observed gender wage gap among overtime workers is smaller than among employed or full-time workers, our decomposition suggests that this may be driven by different hours structure and sorting behaviors across genders. The decomposition of work hours distribution corroborates this implication in that the structure effect explains most of the difference in work hours, compared to the composition effect. These findings underscore the importance of considering the selection effect and incorporating flexible heterogeneity when studying the gender wage gap, consistent with the literature.

\paragraph{Literature Review.} The distribution regression model is a semiparametric generalization of the tobit type-3 model and is a variant of the heckman1979sample selection model. The tobit type-3 model was initially examined in a fully parametric setting, imposing additivity and normality, and estimated by maximum likelihood amemiya1978estimation,amemiya1979estimation. vella1993simple provided a two-step estimator based on estimating the generalized residual from the selection equation and including it as a control function in the outcome equation. honore1997estimation, chen1997semiparametric and lee2006semi provide semi-parametric estimators for this model. fernandez2021nonseparable showed nonparametric identification of local average and distributional effects in the selected population conditional on a control variable, the conditional distribution of the selection variable, even without exclusion restrictions. They relied on exclusion restrictions for an instrumental variable with large support to achieve identification of unconditional effects over the entire population. We focus on unconditional objects and do not require rich variation on the instrumental variable, but impose a restriction on the local dependence between the outcome and selection variable. chen2024estimation adopted a quantile regression model with a censored selection rule, building on the binary selection model of arellano2017quantile. They employ semiparametric quantile regression models for both the outcome and the selection equations. They do not rely on exclusion restrictions, but impose parametric structure on the sorting between the latent outcome and selection variables. In contract, our exclusion restrictions impose semiparametric structure on the sorting. Thus, while our sorting parameter is infinite dimensional, the sorting parameter of chen2024estimation is finite dimensional. Moreover, the quantile regression approach relies on continuity of the outcomes, whereas our distribution regression approach can accommodate continuous, discrete, and mixed outcomes.

\paragraph{Outline.} The rest of the paper is structured as follows. Section 2 analyzes identification under censored sample selection using a local Gaussian representation of the joint distribution of the latent outcome and selection variables. Section 3 presents the distribution regression model with a censored selection rule and provides estimation and inference results. Section 4 reports the results of the empirical application and Section 5 concludes.

Identification under Censored Selection

We first recall the local Gaussian representation (LGR) of CFL. We make use of this representation to analyze identification of the joint distribution of the latent outcome and selection variables under a censored selection rule.

Local Gaussian Representation

Let $(S^{*},Y^{*})$ be two random variables with joint cumulative distribution function (CDF) $F_{S^{*},Y^{*}}$ and marginal CDFs $F_{S^{*}}$ and $F_{Y^{*}}$. In the sample selection model, $Y^*$ will be the latent outcome of interest and $S^*$ the latent selection variable that determines when $Y^*$ is observed. The following lemma from CFL shows that $F_{S^{*},Y^{*}}$ has a representation at each point as a bivariate standard Gaussian distribution with local correlation parameter.

Lemma 1 establishes that $F_{S^{*},Y^{*}}$ can be represented by a sequence of standard bivariate normal distributions with correlation parameters that depend on the evaluation point $(s,y)$. The LGR is a convenient tool to analyze identification of the sample selection model because it connects naturally with the bivariate Gaussian distribution, which is the basis of the HSM. Thus, the local correlation parameter of the bivariate Gaussian distribution is constant and equal to the correlation coefficient, that is, $\rho(s,y) = \rho$ for all $(s,y)$. More generally, the local correlation parameter $\rho(s,y)$ is a measure of local dependence. In particular, $\rho(s,y)=0$ if and only if $F_{S^{*},Y^{*}}(s,y)= F_{S^{*}}(s) F_{Y^{*}}(y)$, i.e., $S^*$ and $Y^*$ are locally independent at $(s,y)$. Furthermore, $\rho(s,y)$ is positive if and only if the correlation between $1\left(S^{*}\leq s\right)$ and $1\left(Y^{*}\leq y\right)$ is positive. In the sample selection model, we will refer to $\rho(s,y)$ as the sorting parameter because it determines the local sign and extent of the selection.

It is convenient to introduce the following notation for the parameters of the LGR associated with the marginal distributions: \[ \mu(s) := \Phi^{-1}(F_{S^{*}}(s)),\;\nu(y) := \Phi^{-1}(F_{Y^{*}}(y)), \] where $\mu(s)\in\bar{\mathbb{R}}$, $\nu(y)\in\bar{\mathbb{R}}$ and $\bar{\mathbb{R}}=\mathbb{R}\cup\left\{ -\infty,\infty\right\} $ is the extended real number line. Thus, we can express the LGR of $F_{S^{*},Y^{*}}$ at point $(s,y)$ as \[ F_{S^{*},Y^{*}}(s,y)\equiv\Phi_{2}(\mu(s),\nu(y);\rho(s,y)). \] We note that the representation can be easily extended to CDFs conditional on covariates by allowing all parameters to depend on the value of covariates, which we shall use later in the identification analysis.

Now, we introduce the sample selection problem with a censored selection rule. The observed variables $(S,Y)$ can be defined in terms of the latent variables $(S^{*},Y^{*})$. We only observe $S^{*}$ censored from below at $0$ and the value of $Y^{*}$ whenever $S^{*}$ is larger than $0$. That is, we observe the random variables $S$ and $Y$ such that

Note that the sample selection is modeled through a censored selection rule with a maximum function rather than using a binary selection model with an indicator function. \ifnum\value{page}=1 \origfootnote{The censored selection rule can be generalized to censoring points other than $0$ and to lower censoring. For example, if the censoring point is $\underline{s} \neq 0$, we can create $S^{**} := S^*-\underline{s}$ and $\tilde S = S - \underline{s}$ such that $\tilde S = \max(S^{**},0)$. If $S$ is lower censored, that is $S = \min(S^*,0)$, we can create $S^{**} := - S^*$ and $\tilde S = -S$, such that $\tilde S = \max(S^{**},0)$.} \else {\arabic{footnote}} \origfootnote{The censored selection rule can be generalized to censoring points other than $0$ and to lower censoring. For example, if the censoring point is $\underline{s} \neq 0$, we can create $S^{**} := S^*-\underline{s}$ and $\tilde S = S - \underline{s}$ such that $\tilde S = \max(S^{**},0)$. If $S$ is lower censored, that is $S = \min(S^*,0)$, we can create $S^{**} := - S^*$ and $\tilde S = -S$, such that $\tilde S = \max(S^{**},0)$.} \fi For example, if $S^*$ is desired work hours and $Y^*$ is offered wage, the observed work hours are censored at $0$ and we only observe the offered wage for those who work a positive number of hours.

Identification via Exclusion Restrictions

The goal is to determine if we can identify the parameters of the LGR of the latent variables $(S^*,Y^*)$, $(\mu(s),\nu(y),\rho(s,y))$, from the joint distribution of the observed variables $(S,Y)$. We first show that $\nu(y)$ and $\rho(s,y)$ are partially identified without further assumptions. Then we establish point identification using exclusion restrictions.

For any $s_{1}>s_{0}\geq0$, we can write the distribution of observed variables as

The selection probabilities pin down $(\mu(s_{0}),\mu(s_{1}))$ in (ref) as $\mu(s_{i})=\Phi^{-1}(P(S\leq s_{i}))$. However, in identifying the three parameters $(\nu(y),\rho(s_{0},y),\rho(s_{1},y))$, there are only two free probabilities (ref) and (ref), which generally leads to partial identification. Using additional values of $S$ and $Y$ does not help for identification as they introduce at least as many unknowns as equations.

We show that identification can be achieved using exclusion restrictions. Let $Z$ be a potential instrumental variable and $F_{S^{*},Y^{*}|Z}$ be the joint CDF of $(S^{*},Y^{*})$ conditional on $Z$. Then, $F_{S^{*},Y^{*}|Z}$ admits the LGR \[ F_{S^{*},Y^{*} \mid Z}(s,y\mid z)\equiv \Phi_{2}(\mu_{z}(s),\nu_{z}(y);\rho_{z}(s,y)), \] where $\mu_{z}(s) := \Phi^{-1}(F_{S^{*} \mid Z}(s \mid z))$ and $\nu_{z}(y) := \Phi^{-1}(F_{Y^{*} \mid Z}(y \mid z))$. We make the following assumption:

Assumption (ref) is the same exclusion restriction as in CFL when $s_{0}=0$. The case where $Z$ is binary is the most challenging for identification as it precludes identification coming from rich variation of $Z$ . When $Z$ is not binary, we only need Assumption (ref) to be satisfied for two values of $Z$. If it holds for more than two values, the model becomes overidentified, making the exclusion restrictions testable. If $s_{0}=0$, non-degeneracy states that sample selection occurs and that $Z$ varies within the selected population, and relevance requires that $Z$ affects the probability of selection. If $s_{0}>0$, non-degeneracy and relevance can be interpreted analogously in terms of the subpopulation with $S > s_0$ instead of the selected population. Outcome exclusion is a standard exclusion restriction, and it holds when $Y^{*}$ is independent of $Z$. Sorting exclusion is a local version of the copula invariance restriction from CFL and holds if selection sorting function at $s=s_{0}$ is independent of $Z$. We refer to CFL for more discussion on the interpretation of this condition. Note that we require the sorting exclusion condition only for $s=s_{0}$, so that $\rho_{z}(s,y)$ can still depend on $z$ when $s\neq s_{0}$. This implies that we can achieve identification of additional parameters from the censored selection model such as $(\mu_z(s_{1}),\rho_z(s_{1},y))$, relying on the same identification condition as in the binary selection model. The most common case is when $s_0=0$, but exclusion might also hold for other values. For example, in the labor supply, we might have variables that determine the decision of working part-time or overtime that affect neither offered wages nor sorting.

We now demonstrate how exclusion restrictions help identify all the parameters. Let Assumption (ref) hold. After imposing the exclusion restriction on the LGR, the target parameters are $(\mu_{0}(s_0),\mu_{1}(s_0),\mu_{0}(s_1),\mu_{1}(s_1),\nu(y),\rho(s_0,y),\rho_0(s_1,y),\rho_1(s_1,y))$. First, the parameters $\mu_{0}(s_0)$, $\mu_{1}(s_0)$, $\mu_{0}(s_1)$, and $\mu_{1}(s_1)$ are identified by \[ \mu_{z}(s)=\Phi^{-1}(\Pr(S\leq s\mid Z=z)),\ s\geq0,\; \ z\in\lbrace0,1\rbrace. \] At $s=s_0$, we have a system of two nonlinear equations and two unknowns \[ \Pr(S>s_0,Y\leq y\mid Z=z)=\Phi_{2}(\Phi^{-1}(\Pr(S>s_0\mid Z=z)),\nu(y);-\rho(s_0,y)), \ \ z\in\lbrace0,1\rbrace, \] where we have used the symmetry properties of $\Phi_2$. Using the same arguments as in CFL, it can be shown that this nonlinear system has a unique solution for $(\nu(y),\rho(s_0,y))$ by the global univariance result of gale1965jacobian. Lastly, $\rho_{z}(s,y)$ for $0 \leq s\neq s_0$ and $z\in\lbrace0,1\rbrace$ is identified from

where $\mu_{z}(s_0)$, $\mu_{z}(s)$, $\nu(y)$ and $\rho(s_0,y)$ are identified from the argument given above. Equation (ref) has unique solution because $\rho\to\Phi_{2}(\mu,\nu;\rho)$ is strictly increasing. In the theorem below, we formulate the identification of $\rho_z(s,y)$, as the identification results for other parameters are stated in CFL.

Distribution Regression with Censored Selection

Model

We consider a semi-parametric version of the LGR with covariates and impose the exclusion restrictions:

where $z'\nu(y) = x'\nu(y)$ and $g(z^{\prime}\rho(s_0,y)) = g(x^{\prime}\rho(s_0,y))$ are the exclusion restrictions, $S^{*}$ is the latent variable that determines selection, $Y^{*}$ is the latent outcome, $X$ is a vector of covariates related to $Y^*$ and $S^*$, $Z=(X,Z_{1})$, $Z_{1}$ is a vector of instrumental variables, the covariates that satify the exclusion restrictions, and $s_0$ is the level of selection variable where the exclusion restriction on $\rho_z(s,y)$ holds. In the empirical application in the next Section, $Y^{*}$ is offered wage, $S^{*}$ is desired work hours, $X$ includes labor market characteristics such as education or experience, $Z_{1}$ includes potential out-of-work income, and $s_0=0$ is the employment participation level.

The bivariate distribution regression (BDR) model (ref) is semiparametric in that the parameters $s\mapsto\mu(s)$, $y\mapsto\nu(y)$ and $(s,y)\mapsto\rho(s,y)$ are function-valued. $u\mapsto g(u)$ is a known link function with range $[-1,1]$ such as the Fisher trasnformation $g(u)=\tanh(u)$. We refer to $-z^{\prime}\mu(s)$ as the selection (hours) equation, $-x^{\prime}\nu(y)$ as the outcome (wage) equation, and $g(z^{\prime}\rho(s,y))$ as the selection sorting (hours-wage sorting) equation. Importantly, note that the covariates of the selection sorting equations are different for $s=s_0$ and $s \neq s_0$ as the exclusion restriction only applies when $s=s_0$.

The parameters $\left(\mu(s),\nu(y)\right)$ can be interpreted similarly to those in probit models. For example, $\nu(y)$ reveals the sign of partial effects of the corresponding covariates on the conditional distribution of the latent outcome, $$ F_{Y^{*} \mid Z}(y \mid z) = \Phi(-x^{\prime}\nu(y)), $$ such as the conditional probability that the latent wage $Y^{*}$ is less than or equal to $y$. Specifically, if $X$ is continuous, \[ \dfrac{\partial F_{Y^{*} \mid Z}(y \mid z)}{\partial x}=-\nu(y)\phi(-x^{\prime}\nu(y)) \] where $\phi$ is the standard normal PDF. The parameter $\mu(s)$ plays an equivalent role in the partial effects on the conditional distribution of the latent selection variable, $$ F_{S^{*} \mid Z}(s \mid z) = \Phi(-z^{\prime}\mu(s)), $$ such as the probability that work hours are less than or equal to $s$. The parameter $\rho(s,y)$ indicates the sign of the effect of the covariates on the sorting, since \[ \dfrac{\partial g(z^{\prime}\rho(s,y))}{\partial x}=\rho(s,y)\dot{g}(z^{\prime}\rho(s,y)) \] where $\dot{g}(u)=\mathrm{d} g(u)/\mathrm{d} \mathrm{u} =1-g(u)^{2}>0$ if $g(u)=\tanh(u)$.

The main parameter of interest in our empirical application is the selection sorting function $(z,s,y) \mapsto g(z'\rho(s,y))$. For illustration, assume that we include only a constant in selection sorting equation, i.e. $g(z'\rho(s,y)) = g(\rho(s,y)) =: \tilde \rho(s,y)$. By definition, $\tilde \rho(s,y)$ determines the conditional local dependence between $S^{*}$ and $Y^{*}$ at the point $(s,y)$. In the wage example, $y \mapsto \tilde \rho(0,y)$ can be interpreted as selection sorting into employment at the wage level $y$. Similarly, setting $s$ to the appropriate level corresponding to full-time or overtime work hours, $y \mapsto \tilde \rho(s,y)$ can be interpreted as sorting into full-time or overtime work at the wage level $y$. Therefore, as mentioned in the introduction, we can study heterogeneous patterns of sorting effects according to selection variable $S^{*}$ and outcome variable $Y^{*}$.

Now, we compare the HSM with censored selection to our model. HSM with outcome exclusion restriction assumes

where $(U,V)$ is standard bivariate Gaussian with correlation $\rho$ conditional on $Z$. The joint distribution of $(S^{*},Y^{*})$ conditional on $Z$ is therefore modeled as \[ F_{S^{*},Y^{*}\mid Z}(s,y\mid z)=\Phi_{2}\left(\frac{s-z^{\prime}\mu}{\sigma_{V}},\frac{y-x^{\prime}\nu}{\sigma_{U}};\rho\right). \] This is a special case of (ref) with \[ -z^{\prime}\mu(s)= (s-z^{\prime}\mu)/\sigma_{V},\ -z^{\prime}\nu(y)=(y-x^{\prime}\nu)/\sigma_{U},\ g(z^{\prime}\rho(s,y))=\rho. \]

We want to emphasize three key points. First, our model is a generalization of HSM with a censored selection rule. Second, HSM only allows for location effects of covariates in the selection and outcome equations, whereas our model permits heterogeneous effects of all covariates with respect to $(s,y)$. Finally, HSM sets the selection sorting parameter as a constant, while we allow for the effect of covariates and coefficients to change with $(s,y)$, representing a significant relaxation of the selection sorting mechanism. Furthermore, in terms of identification, the selection sorting exclusion restriction is trivially satisfied if we set the selection sorting parameter as a constant. The widespread use of the classical model in empirical research, despite its stricter assumptions, suggests that our more flexible selection sorting exclusion restriction can be similarly validated.

Functionals

We can construct functionals of the BDR's parameters and marginal distributions of the covariates, which may be of interest. For example, the marginal distributions of $Y^{*}$ and $S^{*}$ are

where $F_{X}$ and $F_{Z}$ are the CDFs of $X$ and $Z$. We can also construct counterfactual distributions by replacing $(\nu(y),\mu(s))$ or $(F_{X},F_{Z})$ by coefficients and distributions from different population groups. These distributions are useful to decompose the distribution of offered wages or desired work hours between different genders, races or time periods.

We can also use the model to construct distributions for the observed wage by worker type as defined by work hours. For $0\leq\underline{s}<\overline{s}$, \[ F_{Y}(y\mid\underline{s}<S^{*}\leq\overline{s})=F_{Y^*}(y\mid\underline{s}<S^{*}\leq\overline{s})=\frac{F_{S^{*},Y^{*}}(\overline{s},y)-F_{S^{*},Y^{*}}(\underline{s},y)}{F_{S^{*}}(\overline{s})-F_{S^{*}}(\underline{s})}, \] where \[ F_{S^{*},Y^{*}}(s,y)=\int\Phi_{2}(-z^{\prime}\mu(s),-x^{\prime}\nu(y);g(z^{\prime}\rho(s,y)))\mathrm{d} F_{Z}(z). \] We can construct corresponding counterfactual distributions by replacing $\mu(s)$, $\nu(y)$, $\rho(s,y)$, or $F_{Z}$. In the wage application, we will decompose the differences in the observed wage distribution between genders into changes in the worker composition $F_{Z}$ , wage structure $\nu(y)$, hours structure $\mu(s)$, and hours sorting $\rho(s,y)$. The hours structure and sorting effects are new to this model, compared to the decomposition analysis in CFL. We can construct quantiles and other functionals of the distribution of latent and observed outcomes by applying the appropriate operator.

Estimation

To estimate the model parameters and functionals of interest, we assume that we have a random sample of size $n$ from $(S,Y,Z)$, $\{(S_{i},Y_{i},Z_{i})\}_{i=1}^{n}$. In this section we will also assume that $s_0=0$, which corresponds to our empirical application.

Let $(\mathcal{S},\mathcal{Y})$ be the regions of $(S,Y)$ of interest, and let $D_{i}=1(S_{i}>0)$. For example, $\mathcal{Y} = [q_0,q_1]$, where $q_0$ and $q_1$ are quantiles of $Y$. For each $s\in\mathcal{S}$ and $y\in\mathcal{Y}$, we construct the indicators $J_{i}^{s}=1(S_{i}\leq s)$, $\bar{J}_{i}^{s}=1-J_{i}^{s}$, $I_{i}^{y}=1(Y_{i}\leq y)$, and $\bar{I}_{i}^{y}=1-I_{i}^{y}$. Note that $D_{i}=\bar{J}_{i}^{0}$. We replace the arguments of $(\mu(s),\nu(y),\rho(s,y))$ by subscripts to lighten the notation, that is, $(\mu_{s},\nu_{y},\rho_{sy})=(\mu(s),\nu(y),\rho(s,y))$.

Estimation is based on forming the joint conditional likelihood of $J_{i}^{s}$ and $I_{i}^{y}$ given $Z_i$ taking into account that $I_{i}^{y}$ is only observed when $D_{i}=1$. Thus, let $\xi_{sy} := (\mu_0,\mu_s,\nu_y,\rho_{0y},\rho_{sy})$, the likelihood is

We break down the maximization of the likelihood in 3 steps, which are illustrated in Figure (ref). The first step is a probit regression to estimate $\mu_{s}$ for each $s \in \{0\} \cup\mathcal{S}$. When $s=0$, we are estimating the probability of employment, which is identical to the first step of the estimation algorithm in CFL. The second step is a probit regression with sample selection correction to estimate $(\nu_{y},\rho_{0y})$ for each $y \in \mathcal{Y}$, which is also the same as the second step in CFL. As mentioned previously, the likelihood function incorporates adjustment to the sample selection. The final step is new relative to CFL and runs a bivariate probit regression to estimate $\rho_{sy}$ for each $(s,y) \in \mathcal{S} \times \mathcal{Y}$. This step estimates the local correlation parameters of the joint distribution plugging-in the parameters of the marginal distributions of $S^{*}$ and $Y^{*}$ estimated in the first and second steps. The 3 steps are summarized in Algorithm 1.

Estimators of the functionals of interest in Section 3.2 are constructed from the estimators of the parameters using the plug-in method. For example, the marginal distributions of the latent outcome and selection variable $(Y^{*},S^{*})$ are estimated with \[ \hat{F}_{Y^{*}}(y)=\dfrac{1}{n}\sum_{i=1}^{n}\Phi(-X_{i}^{\prime}\hat{\nu}(y))\;\text{for }y\in\overline{\mathcal{Y}},\;\hat{F}_{S^{*}}(s)=\dfrac{1}{n}\sum_{i=1}^{n}\Phi(-Z_{i}^{\prime}\hat{\mu}(s))\;\text{for }s\in\overline{\mathcal{S}} \] and the distribution of observed wage by worker with type is \[ \hat{F}_{Y}(y\mid\underline{s}<S^{*}\leq\overline{s})=\frac{\hat{F}_{S^{*},Y^{*}}(\overline{s},y)-\hat{F}_{S^{*},Y^{*}}(\underline{s},y)}{\hat{F}_{S^{*}}(\overline{s})-\hat{F}_{S^{*}}(\underline{s})},\;0\leq\underline{s}<\bar{s},y\in\overline{\mathcal{Y}} \] where \[ \hat{F}_{S^{*},Y^{*}}(s,y)=\dfrac{1}{n}\sum_{i=1}^{n}\Phi_{2}(-Z_{i}^{\prime}\hat{\mu}(s),-X_{i}^{\prime}\hat{\nu}(y);g(Z_{i}^{\prime}\hat{\rho}(s,y))). \]

Uniform Inference

We present asymptotic theory for the estimator of $(s,y) \mapsto \rho_{sy}$, the new function-valued parameter of interest, and show how to construct confidence bands for this parameter using multiplier bootstrap. The proof in Appendix (ref) includes results for all parameters in the model, $(\eta_{sy},\rho_{sy})$, where $\eta_{sy}:=(\mu_0,\mu_{s},\nu_y,\rho_{0y})$. We only present the result for $\rho_{sy}$ in the main text as it is the most novel result relative to CFL.

Let $L_{1}(\mu_s)$, $L_{2}(\theta_y,\mu_0)$, $L_{3}(\rho_{sy},\eta_{sy})$, $H_{1s}$, $H_{2y}$, $H_{3sy}$ and $\Sigma_{\rho}(s,y)$ be the functions defined in Appendix (ref). We make the following assumption to derive the asymptotic distribution of $\hat \rho_{sy}$:

Let $\mathcal{SY}:= \mathcal{S} \times \mathcal{Y}$, $\ell^{\infty}(\mathcal{SY})^{d_{\rho_{sy}}}$ be the set of bounded functions on $\mathcal{SY}$, $\leadsto$ denote weak convergence (in distribution), and $\to_{\Pr}$ denote convergence in probability. The first main result of this section is a functional central limit theorem for $\hat{\rho}_{sy}$.

We show how to construct confidence bands for function-valued parameters 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 $\rho_{sy}$ , that is the function $(s,y)\mapsto c^{\prime}\rho_{sy}$, where $c\in\mathbb{R}^{d_{\rho_{sy}}}$ and $(s,y)\in\mathcal{SY}$. The set $\textrm{CB}_{p}(c^{\prime}\rho_{sy})$ is a uniform (asymptotic) $p$-confidence band for $c^{\prime}\rho_{sy}$ on $\mathcal{SY}_0 \subseteq \mathcal{SY}$ if \[ \Pr\left[c^{\prime}\rho_{sy}\in \textrm{CB}_{p}(c^{\prime}\rho_{sy}),\forall(s,y)\in\mathcal{SY}_0\right]\to p \text{ as } n \to \infty. \] Let $\hat{\Sigma}_{\rho}(s,y)$ be a uniformly consistent estimator of $\Sigma_{\rho}(s,y)$, such as the one given in (ref). We construct \[ \textrm{CB}_{p}(c^{\prime}\rho_{sy})=c^{\prime}\hat{\rho}_{sy}\pm cv(p)\sqrt{c'\hat{\Sigma}_{\rho}(s,y)c/n}, \] where $cv(p)$ is a consistent estimator of the $p$-quantile of the maximal t-statistic \[ t_{\mathcal{SY}_0}=\sup_{(s,y)\in\mathcal{SY}_0} \sqrt{n}\frac{|c^{\prime}\hat{\rho}_{sy}-c^{\prime}\rho_{sy}|}{\sqrt{c'\hat{\Sigma}_{\rho}(s,y)c}}. \] We can obtain the critical value from the limit distribution of the stochastic process in Theorem (ref).

In practice, it is convenient to estimate the critical value using resampling methods. The multiplier bootstrap method of gine1984some is computationally attractive in our setting because it does not require parameter re-estimation and therefore avoids the nonlinear optimization in our estimation algorithm. The multiplier bootstrap is implemented using the following algorithm. Confidence bands for other functionals of the model parameter can be constructed using a similar method.

We make the following assumption about the multipliers:

We establish a functional central limit theorem for the bootstrap for $\hat \rho_{sy}$. Here we use $\leadsto_{\Pr}$ to denote bootstrap consistency, i.e. weak convergence conditional on the data in probability, which is formally defined in van1996weak. \ifnum\value{page}=1 \origfootnote{We recall this definition in Appendix (ref).} \else {\arabic{footnote}} \origfootnote{We recall this definition in Appendix (ref).} \fi

Employment Sorting and Wage Decompositions in the U.K.

We analyze sorting into full-time and overtime work across gender, marital status, and time in the UK labor market. Furthermore, we conduct wage decompositions to compare the wage distributions of men and women for full-time and overtime works.

Data

We use the same dataset as CFL, which contains repeated cross-sectional observations spanning from 1978 to 2013 data. The sources of the data are the Family Expenditure Survey (FES) from 1978 to 2001, the Expenditure and Food Survey (EFS) from 2002 to 2007, and the Living Costs and Food Survey (LCFS) from 2008 to 2013. The sample is constructed similarly to the previous studies that used the FES, such as gosling2000changing, blundell2003interpreting, and blundell2007changes. It includes individuals aged 23 to 59 years, excluding full-time students, self-employed individuals, those married with an absent spouse, and those with missing education or employees with missing wage information. This results in a sample of 258,900 observations, with 139,504 women and 119,396 men.

The outcome variable of interest, $Y$, is the logarithm of the real hourly wage rate. This is calculated as the ratio of weekly usual gross nominal earnings to weekly usual work hours, adjusted for inflation using the UK quarterly retail price index. The selection variable, $S$, is the weekly usual work hours. The covariates, $X$, include indicators for age at which individuals left school, a quartic polynomial in age, an indicator for being married or cohabiting, the number of children categorized by age, survey year indicators, and region. Table 1 presents descriptive statistics for these variables.

The instrumental variable, $Z_{1}$, is the potential out-of-work income benefit interacted with the marital status indicator, as previously used in blundell2003interpreting, blundell2007changes and CFL. This variable is constructed using the Institute of Fiscal Studies (IFS) tax and welfare-benefit simulation model (TAXBEN) for each individual. TAXBEN calculates the income of a tax unit if the individual were out-of-work, factoring in eligible unemployment and housing benefits. The primary sources of variation are the household demographic composition, housing costs, and policy changes across regions over time. We assume the outcome and selection sorting exclusion restrictions at $s_0=0$, that is, we assume that, conditional on observed covariates, the offered wage and the dependence between offered wage and desired work hours at $s=0$ do not depend on the benefit level. While the plausibility of the exclusion restriction is controversial, we assume these conditions are satisfied and refer to blundell2003interpreting and blundell2007changes for a discussion on the validity of outcome exclusion restriction. For the selection sorting exclusion restriction, we motivate its use by noting the widespread application of the classical HSM in empirical research, which trivially satisfies the restriction.

Empirical specifications

We estimate the censored distribution regression model separately for men and women. The selection (work hours) equation includes all the covariates mentioned earlier, while the outcome (wage) equation is identical except for the excluded covariates. Given that estimating the parameter of the selection sorting function is notoriously more challenging than estimating the parameters of the selection and outcome equations, we consider three simplified specifications for this function. The covariates included in the index $z^{\prime}\rho(s,y)$ are

Using these specifications, we estimate the selection sorting parameters $\rho(s,y)$ on $\mathcal{SY}_0 = \mathcal{S}_0 \times \mathcal{Y}_0 $ and plug-in these estimates to the selection sorting function $(s,y)\to g\left(z^{\prime}\rho(s,y)\right)$ for the population group of interest defined by the value of $z$, where $\mathcal{S}_0 = \{0,34,40\}$, the thresholds corresponding to part-time, full-time and overtime work, and $\mathcal{Y}_0$ is the interval between the sample quantiles of $Y$ with orders $.1$ and $.9$. For example, using specification 2, we estimate the selection sorting function for married men $(s,y)\to g\left(\rho_{M}(s,y)\right)$ using the sample of men. In Algorithms 1 and 2, we use a grid $\overline{\mathcal{SY}}_0 = \mathcal{S}_0 \times \overline{\mathcal{Y}}_0$, where $\overline{\mathcal{Y}}_0$ is a grid of sample quantiles of $Y$ with indices $\left\{ 0.1,0.11,...,0.9\right\} $. We present 95% confidence bands for the selection sorting function using Algorithm 2 with $B=500$ bootstrap repetitions.

Our specification of $\mathcal{S}_0$ involves the definition of part-time, full-time, and overtime work, with thresholds set at 34 and 40 hours. This means we classify individuals working 34 hours or less as part-time workers, those working from 35 to 40 hours as full-time workers, and those working 41 hours or more as overtime workers. Although we have not been able to find formal definitions of these worker types, we have adopted thresholds that are commonly accepted in the literature (mulligan2008selection, maasoumi2019gender, blundell2021wages). Moreover, as shown in Figure (ref), which presents histograms of work hours by gender, our data exhibits clear bunching at the thresholds of 34 and 40 hours. Similarly, garnero2014part used thresholds defined by the bunching observed in their histograms to categorize worker types. As provided in Table 1, 83 percent of men and 66 percent of women are employed. Among the employed, the proportions of part-time, full-time, overtime workers are 5, 50, 45 percent for men and 50, 41, 9 percent for women, respectively. This shows a significant difference in worker types by gender, with part-time work being more prevalent among women and overtime work more common among men. These findings motivate further analysis of selection patterns according to work hours by gender.

Main results

Employment sorting

Figure (ref) presents estimates of selection sorting functions using Specification 1. The results for selection sorting into employment in the first column replicate those reported in CFL. The estimated sorting functions indicate that selection into full-time and overtime work exhibit distinct patterns compared to selection into employment. Specifically, men show positive selection into employment and full-time work but negative selection into overtime work. In contrast, women exhibit negative selection into employment and overtime work but positive selection into full-time work. These findings showcase how the censored selection model is able to uncover heterogeneous selection sorting patterns based on work hours, which cannot be captured using a binary selection model. In the subsequent section, we discuss results using Specifications 2 to 4 to identify the primary drivers of the different sorting patterns observed in Specification 1, focusing on marital status and changes over time.

Figures (ref) and (ref) present the results for selection into full-time and overtime work using Specification 2. We first discuss the findings for selection into full-time work. Figure (ref) shows that men exhibit positive sorting across the wage distribution, whereas women show positive sorting only at the lower wage quantiles. Figure (ref) indicates that these sorting patterns are particularly pronounced among single men and married women. We interpret this as a result of single men and married women having greater flexibility in selecting their work hours. Traditionally, men have been the primary financial supporters of their families after marriage, which may allow single men and married women more freedom in choosing work hours, leading to significant positive sorting patterns. Thus, single men and married women are relatively more likely to work full time when they have relative comparative advantage based on unobservable characteristics, say ability.

To explain why married women exhibit positive selection only at the bottom of the distribution, we adopt the concepts of voluntary and involuntary part-time work. According to dunn2018chooses, part-time workers can be classified into two categories based on their reasons for working part-time. Involuntary part-time workers seek full-time work but work part-time because they can only find part-time jobs. In contrast, voluntary part-time workers prefer part-time work due to non-economic reasons, such as childcare, school, or medical issues. cam2012involuntary, using UK Labour Force Survey data, examined the determinants of involuntary part-time work. Their findings suggest that women working part-time with dependent children are more likely to be voluntary, especially when they are in a couple. Based on this, we posit that single male part-time workers are predominantly involuntary, meaning only those with high potential wages are likely to work full-time. On the other hand, married women in part-time work are a mix of voluntary and involuntary workers, which explains why we observe a selection effect only at the bottom of the distribution.

For selection into overtime work, Figure (ref) shows that both genders exhibit a negative selection effect, which is more pronounced at the top of the wage distribution. Figure (ref) further reveals that this negative selection effect is more pronounced for married men and women. One potential explanation for this pattern is that higher-wage jobs are more likely to be task-based in the sense that the work hours are determined by the time needed to complete the assigned tasks. Consequently, workers who require overtime at high wage levels are more likely to be those with relatively lower ability, which is reflected in the negative selection effect observed in our estimates.

Figure (ref) shows the selection sorting function using Specification 3, revealing that sorting has become more positive and homogeneous across wage quantiles over time. For women, this aligns with maasoumi2019gender, who showed that selection into full-time shifts from negative to positive over time. However, while they estimate a single selection parameter for each year, we estimate a selection sorting function across wage quantiles. Figures (ref) and (ref) in the Appendix (ref) indicate that these sorting trends are primarily driven by single men and married women, similar to the findings in Specifications 1 and 2.

Wage decomposition

We next use the DR model to decompose changes in the distribution of the observed wages between women and men for full-time and overtime workers. The goal is to investigate the factors contributing to the gender wage gap according to worker type. We do not include the result for part-time workers, as there are too few part-time male workers in our dataset, leading to unreliable results.

We identify four components corresponding to different inputs of the DR model: (1) wage structure $\nu(y)$; (2) selection sorting $\rho(s,y)$; (3) selection structure $\mu(s)$; and (4) composition $F_{Z}$. Let $Q_{Y\langle t,j,r,k\rangle}^{(\underline{s},\overline{s}]}$ represent the (counterfactual) quantiles of workers with work hours within the $(\underline{s},\overline{s}]$ interval, with $t$-wage structure, $j$- hours sorting, $r$-hours structure, and $k$-composition. The actual quantile in group $t$ corresponds to $Q_{Y\langle t,t,t,t\rangle}^{(\underline{s},\overline{s}]}$. We consider two groups, labeled 0 and 1, representing different demographic populations, such as women and men. The observed wage quantile difference between group 1 and 0 can be decomposed as:

In this decomposition, the first term on the right-hand side captures the wage structure effect, the second reflects the hours-wage sorting effect, the third represents the hours structure effect, and the fourth indicates the composition effect. This generalizes the classical Oaxaca-Blinder decomposition (kitagawa1955components; oaxaca1973male; blinder1973wage) to analyze distributional gap of the outcome and to take sample selection into account. Furthermore, since we are adopting censored selection rule, we can decompose wage distribution gaps according to worker type by work hours, in contrast to CFL decomposing wage distribution gap only of employed. While it is well-known that the order of component extraction in such decompositions can influence the results, we estimate the decomposition using different orderings of the components and find that the main conclusions remain consistent.

Let $F_{Y\langle t,j,r,k\rangle}^{(\underline{s},\overline{s}]}$ be the (counterfactual) distribution of wages, which can be expressed as

where $\nu_{t}$ is the coefficient of the wage equation in group $t$, $\rho_{j}$ is the coefficient of the sorting function in group $j$, $\mu_{r}$ is the coefficient of the selection equation in group $r$, and $F_{Z_{k}}$ is the distribution of characteristics in group $k$. We construct a plug-in estimator of $F_{Y\langle t,r,j,k\rangle}^{(\underline{s},\overline{s}]}(y)$ by appropriately using the estimates of model parameters and covariate distributions from the two groups. Lastly, we use the generalized quantile operator $$Q_{\tau}(F)= \int_{-\infty}^0 1(F(y) \leq \tau) \mathrm{d} y + \int_{0}^{\infty} 1(F(y) > \tau) \mathrm{d} y, $$ to estimate $Q_{Y\langle t,r,j,k\rangle}^{(\underline{s},\overline{s}]}(\tau)$. \ifnum\value{page}=1 \origfootnote{The operator $Q_\tau$ monotonically arrange the estimator of the distribution before applying the (left) inverse.} \else {\arabic{footnote}} \origfootnote{The operator $Q_\tau$ monotonically arrange the estimator of the distribution before applying the (left) inverse.} \fi

Figures (ref) and (ref) present the estimated quantile functions for observed wages among full-time and overtime working men and women, along with the relative contributions of each component to the wage distribution gap between men (group 1) and women (group 0) based on Specification 2. For both full-time and overtime workers, men exhibit higher wage quantiles than women. However, note that the wage gap at each quantile index is smaller for overtime workers compared to full-time workers starting from $\tau=0.14$ and onward. Specifically, the wage gap for full-time workers is up to 2.5 times larger than that for overtime workers at the same quantile index. In the decomposition analysis below, we explore the components contributing to the reduced gender wage gap among overtime workers. In the decomposition results, a positive contribution indicates that the component accounts for a proportion of the estimated gap, while a negative contribution suggests that the component is associated with narrowing the estimated gap.

For full-time workers, most of the gender wage gap is attributed to differences in the wage structure, that is, differences in returns to observed charactersitics. However, the hours structure and hours-wage sorting effects also explain an important portion of the gap. At lower wage quantiles, the hours structure has a negative contribution. This implies that women with similar observed characteristics to men are less likely to work full-time, leading to a reduced observed wage gap. At higher wage quantiles, the hours-wage sorting effect partially explains the gap, similar to the decomposition result for the employed in CFL.

For overtime workers, the percentages exceed $100$ or fall below $-100$ because the wage distribution gaps generated by the counterfactual distributions are larger than the observed gap, which is used as the denominator. The decomposition procedure shows that if we impose the women's wage structure onto the men's wage distribution, the resulting counterfactual distribution yields much lower quantiles than the observed quantiles for women. \ifnum\value{page}=1 \origfootnote{Note that $[Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]}]/[Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}] > 1$ if and only if $Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]} < Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}$.} \else {\arabic{footnote}} \origfootnote{Note that $[Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]}]/[Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}] > 1$ if and only if $Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]} < Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}$.} \fi The largely negative aggregate selection effect --- which includes both the hours structure and hours-wage sorting effects --- suggests that selection behavior is a key factor in reducing the observed wage gap. For example, men with similar observed characteristics to women are more likely to work overtime, which lowers the observed wage quantiles for men relative to women as both gender exhibit negative sorting into overtime work. In fact, the estimated probability of working overtime with male hours structure and male composition is 0.37, whereas the counterfactual probability with female hours structure and male composition is 0.07. This supports the idea that women with the same observed characteristics as men are less likely to work overtime.

Work hours decomposition

Lastly, we decompose the difference in the distribution functions of observed work hours between genders using counterfactual distributions. Let

where $\mu_{r}$ is the coefficient of the selection equation in group $r$, and $F_{Z_{k}}$ is the distribution of characteristics in group $k$. We decompose the difference in the distribution functions as

The first term is the structure effect, showing differences in how work hours are selected across genders, and the second term is the composition effect. \ifnum\value{page}=1 \origfootnote{Note that we are using distribution functions for work hours decomposition. Therefore, we subtract the women's distribution function from the men's, whereas the order of subtraction is reversed when using quantile functions.} \else {\arabic{footnote}} \origfootnote{Note that we are using distribution functions for work hours decomposition. Therefore, we subtract the women's distribution function from the men's, whereas the order of subtraction is reversed when using quantile functions.} \fi

The distribution of work hours for men first-order stochastically dominates the distribution for women. The majority of the difference is explained by differences in structure, whereas differences in composition have very little explanatory power. This implies that differences in how men and women choose work hours, such as according to education level, marital status, and number of children, lead to the gender gap in work hours rather than differences in characteristics. This aligns with the wage decomposition result in Section 4.3.2 in that composition effects have little explanatory power in gender wage gap.

Conclusion

We propose a distribution regression model with censored sample selection rule. Compared to the classical HSM, our model enables the analysis of the distribution, moves away from parametric Gaussian assumptions on the error structure, and accommodates rich patterns of heterogeneity in the effects of covariates on the outcome and selection process. Compared to the distribution regression model using a binary selection rule, our approach allows for the examination of selection effects based on both the outcome and selection variables. We apply the developed model to wage and work hours data from the U.K. to investigate selection into full-time and overtime work. Additionally, we conduct a decomposition of wage and work hours distributions by gender for these worker types.