The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
79,790 characters
Distribution Regression with Censored Selection$^*$
\title[DR with Censored Selection]{Distribution Regression with Censored Selection$^*$}
\author[Fern\'andez-Val \and Hong]{Iv\'an Fern\'andez-Val \and Seoyun Hong$^\dag$}
\date{\today}
\thanks{$^*$ We thank Kevin Lang, Ruli Xiao, and seminar and conference participants at Boston University, Munich Econometrics Workshop 2023, Midwest Econometrics Group Conference 2024, and BC-BU Greenline Econometrics Workshop 2024 for helpful comments.}
\thanks{$^\dag$ Boston University}
\begin{abstract}
We develop a distribution regression model with a censored selection rule, offering a semi-parametric generalization of the Heckman selection model. Our approach applies to the entire distribution, extending beyond the mean or median, accommodates non-Gaussian error structures, and allows for heterogeneous effects of covariates on both the selection and outcome distributions. By employing a censored selection rule, our model can uncover richer selection patterns according to both outcome and selection variables, compared to the binary selection case. We analyze identification, estimation, and inference of model functionals such as sorting parameters and distributions purged of sample selection. An application to labor supply using data from the UK reveals different selection patterns into full-time and overtime work across gender, marital status, and time. Additionally, decompositions of wage distributions by gender show that selection effects contribute to a decrease in the observed gender wage gap at low quantiles and an increase in the gap at high quantiles for full-time workers. The observed gender wage gap among overtime workers is smaller, which may be driven by different selection behaviors into overtime work across genders.
\pagenumbering{gobble}
\end{abstract}
\maketitle
\newpage
\pagenumbering{arabic}
\section{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 \citet{heckman1974shadow},
this model relies on parametric Gaussian error structures and assumes
homogeneous effects of exogenous covariates. Recently, \citet{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. \citet{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 (\citet{blank1990parttime}, \citet{hirsch2005part}). On the other hand, there is literature examining how long work hours affect the gender wage gap among highly skilled workers (\citet{noonan2005pay}, \citet{bertrand2010dynamics}, \citet{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, \citet{mulligan2008selection} find weaker evidence of convergence in the gender wage gap after accounting for sample selection. Similarly, \citet{arellano2017quantile}, \citet{maasoumi2019gender}, and \citet{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 \citet{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
\renewcommand{\arabic{footnote}}{\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
(\citet{kitagawa1955components}; \citet{oaxaca1973male}; \citet{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{\textbf{Literature Review.}} The distribution regression model is a semiparametric generalization of the tobit
type-3 model and is a variant of the \citet{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 \citep[e.g.,][]{amemiya1978estimation,amemiya1979estimation}.
\citet{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. \citet{honore1997estimation}, \citet{chen1997semiparametric} and \citet{lee2006semi} provide semi-parametric estimators for this
model. \citet{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. \citet{chen2024estimation} adopted a quantile regression model with a censored selection rule, building on the binary selection model of \citet{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 \citet{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{\textbf{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.
\section{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.
\subsection{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.
\begin{lem}[LGR]\label{lemma:lgr}
$F_{S^{*},Y^{*}}$ can be represented as
\[
F_{S^{*},Y^{*}}(s,y)=P(S^{*}\leq s,Y^{*}\leq y)\equiv\Phi_{2}\left(\Phi^{-1}(F_{S^{*}}(s)),\Phi^{-1}(F_{Y^{*}}(y));\rho(s,y)\right), \ \ (s,y) \in \mathbb{R}^2,
\]
where $\Phi_{2}(\cdot,\cdot;\rho)$ is the joint CDF of a standard
bivariate Gaussian random variable with correlation parameter $\rho\in[-1,1]$. In the previous representation, $\rho(s,y)\in[-1,1]$ is a local correlation parameter that
varies with the evaluation point $(s,y)$.
\end{lem}
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
\begin{eqnarray*}
S & = & \max(S^{*},0),\\
Y & = & Y^{*}\ \text{ if }\ S>0.
\end{eqnarray*}
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
\renewcommand{\arabic{footnote}}{\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.
\subsection{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
\begin{eqnarray}
\Pr(S\leq s_{i}) & = & \Phi(\mu(s_{i})),\quad i\in\{0,1\}, \label{eq:sel-probs} \\
\Pr(S>s_{0},Y\leq y) & = & \Phi(\nu(y))-\Phi_{2}(\mu(s_{0}),\nu(y);\rho(s_{0},y)),\label{eq:prob1}\\
\Pr(s_{0}<S\leq s_{1},Y\leq y) & = & \Phi_{2}(\mu(s_{1}),\nu(y);\rho(s_{1},y))-\Phi_{2}(\mu(s_{0}),\nu(y);\rho(s_{0},y)) \label{eq:prob2}.
\end{eqnarray}
The selection probabilities pin down $(\mu(s_{0}),\mu(s_{1}))$ in
\eqref{eq:sel-probs} 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 \eqref{eq:prob1} and \eqref{eq:prob2}, 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:
\begin{assumption}[Exclusion Restrictions]\label{ass:excl} Assume that $Z$ is binary and for some $s_{0}\geq0$:
\begin{enumerate}
\item \textit{Non-degeneracy: $0<\Pr(S>s_{0})<1$ and $0<\Pr(Z=1\mid S>s_{0})<1$ }
\item \textit{Relevance: $\Pr(S>s_{0}\mid Z=0)<\Pr(S>s_{0}\mid Z=1)<1$ }
\item \textit{Outcome exclusion: $\nu_{z}(y)=\nu(y)$ for all $y\in\mathbb{R}$
and $z\in\{0,1\}$ }
\item \textit{Sorting exclusion: $\rho_{z}(s_{0},y)=\rho(s_{0},y)$ for
all $y\in\mathbb{R}$ and $z\in\{0,1\}$. }
\end{enumerate}
\end{assumption}
Assumption \ref{ass:excl} 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{ass:excl} 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{ass:excl} 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
\citet{gale1965jacobian}. Lastly, $\rho_{z}(s,y)$ for $0 \leq s\neq s_0$ and $z\in\lbrace0,1\rbrace$
is identified from
\begin{equation}\label{eq:rho_id}
\Pr(0 < S\leq s,Y\leq y\mid Z=z)=\Phi_{2}(\mu_{z}(s),\nu(y);\rho_{z}(s,y)) - \Phi_{2}(\mu_{z}(s_0),\nu(y);\rho(s_0,y)),
\end{equation}
where $\mu_{z}(s_0)$, $\mu_{z}(s)$, $\nu(y)$ and $\rho(s_0,y)$ are
identified from the argument given above. Equation \eqref{eq:rho_id} 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.
\begin{thm}[Identification of $\rho_z(s,y)$]\label{thm:id}
Suppose that Assumption 1 holds. Then, $\rho_z(s,y)$ is
point identified as the unique solution to \eqref{eq:rho_id}, for $z\in\{0,1\}$, $0 \leq s \neq s_0$ and $y \in \mathbb{R}$.
\end{thm}
\section{Distribution Regression with Censored Selection}
\subsection{Model }
We consider a semi-parametric version of the LGR with covariates and impose the exclusion restrictions:
\begin{equation}\label{eq:drmodel}
F_{S^{*},Y^{*}\mid Z}(s,y\mid z) = \Phi_{2}(\mu_{z}(s),\nu_{z}(y);\rho_{z}(s,y)) = \Phi_{2}(-z^{\prime}\mu(s),-z^{\prime}\nu(y);g(z^{\prime}\rho(s,y))),
\end{equation}
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 \eqref{eq:drmodel} 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
\begin{eqnarray*}
Y^{*} & = & X^{\prime}\nu+\sigma_{U}U\\
S^{*} & = & Z^{\prime}\mu+\sigma_{V}V
\end{eqnarray*}
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 \eqref{eq:drmodel}
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.
\subsection{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
\begin{align*}
F_{Y^{*}}(y) & =\int\Phi(-x^{\prime}\nu(y))\mathrm{d} F_{X}(x), \ \
F_{S^{*}}(s) =\int\Phi(-z^{\prime}\mu(s))\mathrm{d} F_{Z}(z),\quad s\geq0,
\end{align*}
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.
\subsection{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
\begin{align*}
f_i(\xi_{sy})=\Phi(-Z_i^\prime\mu_0)^{1-D_i}&\times \left[\Phi_{2}(Z_{i}^{\prime}\mu_{0},-X_{i}^{\prime}\nu_{y};-g(X_{i}^{\prime}\rho_{0y})) - \Phi_{2}(Z_{i}^{\prime}\mu_{s},-X_{i}^{\prime}\nu_{y};-g(Z_{i}^{\prime}\rho_{sy})) \right]^{I^y_{i} J^s_{i}}
\\
&\times\left[\Phi_{2}(Z_{i}^{\prime}\mu_{0},X_{i}^{\prime}\nu_{y};g(X_{i}^{\prime}\rho_{0y})) - \Phi_{2}(Z_{i}^{\prime}\mu_{s},X_{i}^{\prime}\nu_{y};g(Z_{i}^{\prime}\rho_{sy}))\right] ^{\bar I^y_{i} J^s_{i}}
\\
&\times \Phi_{2}(Z_{i}^{\prime}\mu_{s},-X_{i}^{\prime}\nu_{y};-g(Z_{i}^{\prime}\rho_{sy}))^{I^y_{i} \bar J^s_{i}} \\
& \times \Phi_{2}(Z_{i}^{\prime}\mu_{s},X_{i}^{\prime}\nu_{y};g(Z_{i}^{\prime}\rho_{sy}))^{\bar I^y_{i}\bar J^s_{i}}
\end{align*}
We break down the maximization of the likelihood in 3 steps, which are illustrated in Figure \ref{fig1}. 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.
\bigskip
\begin{figure} [ht!]
\centering
\begin{subfigure}{0.3\textwidth}
\includegraphics[width=\textwidth]{Step1_update.png}
\caption{Step 1}
\end{subfigure}
\hfill
\begin{subfigure}{0.3\textwidth}
\includegraphics[width=\textwidth]{Step2_updated.png}
\caption{Step 2}
\end{subfigure}
\hfill
\begin{subfigure}{0.3\textwidth}
\includegraphics[width=\textwidth]{Step3_updated.png}
\caption{Step 3}
\end{subfigure}
\caption{Estimation steps} \label{fig1}
\end{figure}
\newpage
\begin{algorithm}[Three-Step CDR Method]\label{alg:est}
Let $\overline{\mathcal{S}}$ and $\overline{\mathcal{Y}}$ be finite grids covering $\mathcal{S}$ and $\mathcal{Y}$, respectively.
\begin{enumerate}
\item For each $s\in \{0\} \cup \overline{\mathcal{S}}$, estimate $\mu_{s}$
by running probit regressions of $J_{i}^{s}$ on $Z_{i}$.
\begin{align*}
\hat{\mu}_{s} & =\operatorname*{arg\,max}_{\mu_{s}}\dfrac{1}{n}\sum_{i=1}^{n}\bar{J}_{i}^{s}\log\Pr(S_{i}>s \mid Z_i)+J_{i}^{s}\log\Pr(S_{i}\leq s \mid Z_i)\\
& =\operatorname*{arg\,max}_{\mu_{s}}\dfrac{1}{n}\sum_{i=1}^{n}\bar{J}_{i}^{s}\log\Phi(Z_{i}^{\prime}\mu_{s})+J_{i}^{s}\log\Phi(-Z_{i}^{\prime}\mu_{s}).
\end{align*}
\item For each $y \in \overline{\mathcal{Y}}$, estimate $\theta_{y}=(\nu_{y},\rho_{0y})$
by a probit regression with sample selection correction.
\begin{align*}
\hat{\theta}_{y} & =\operatorname*{arg\,max}_{\theta_{y}}\dfrac{1}{n}\sum_{i=1}^{n}D_{i}\left[\bar{I}_{i}^{y}\log P(S_{i}>0,Y_{i}>y \mid Z_i) +I_{i}^{y}\log P(S_{i}>0,Y_{i}\leq y \mid Z_i)\right]\\
& =\operatorname*{arg\,max}_{\theta_{y}}\dfrac{1}{n}\sum_{i=1}^{n}D_{i}\left[\bar{I}_{i}^{y}\log\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},X_{i}^{\prime}\nu_{y};g(X_{i}^{\prime}\rho_{0y}))+I_{i}^{y}\log\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},-X_{i}^{\prime}\nu_{y};-g(X_{i}^{\prime}\rho_{0y}))\right].
\end{align*}
\item For each $(s,y) \in \overline{\mathcal{S}} \times \overline{\mathcal{Y}}$, estimate
$\rho_{sy}$ by a bivariate probit regression,
\begin{align*}
\hat{\rho}_{sy}= & \operatorname*{arg\,max}_{\rho_{sy}}\dfrac{1}{n}\sum_{i=1}^{n}D_{i}\Big[\bar{J}_{i}^{s}\bar{I}_{i}^{y}\log P(S_{i}>s,Y_{i}>y \mid Z_i)+\bar{J}_{i}^{s}I_{i}^{y}\log P(S_{i}>s,Y_{i}\leq y \mid Z_i)\\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+J_{i}^{s}\bar{I}_{i}^{y}\log P(0 < S_{i}\leq s,Y_{i}>y \mid Z_i)+J_{i}^{s}I_{i}^{y}\log P(0 < S_{i}\leq s,Y_{i}\leq y \mid Z_i)\Big]\\
= & \operatorname*{arg\,max}_{\rho_{sy}}\dfrac{1}{n}\sum_{i=1}^{n}D_{i}\Big[\bar{J}_{i}^{s}\bar{I}_{i}^{y}\log\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},X_{i}^{\prime}\hat{\nu}_{y};g(Z_{i}^{\prime}\rho_{sy}))+\bar{J}_{i}^{s}I_{i}^{y}\log\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},-X_{i}^{\prime}\hat{\nu}_{y};-g(Z_{i}^{\prime}\rho_{sy}))\\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+J_{i}^{s}\bar{I}_{i}^{y}\log\left\{ \Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},X_{i}^{\prime}\hat{\nu}_{y};g(X_{i}^{\prime}\hat{\rho}_{0y}))-\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},X_{i}^{\prime}\hat{\nu}_{y};g(Z_{i}^{\prime}\rho_{sy}))\right\} \\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+J_{i}^{s}I_{i}^{y}\log\left\{ \Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},-X_{i}^{\prime}\hat{\nu}_{y};-g(X_{i}^{\prime}\hat{\rho}_{0y}))-\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},-X_{i}^{\prime}\hat{\nu}_{y};-g(Z_{i}^{\prime}\rho_{sy}))\right\} \Big].
\end{align*}
\end{enumerate}
\end{algorithm}
\bigskip
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))).
\]
\begin{remark}[Computation] When implementing Step 3 of Algorithm 1, we encountered the issue that some
probabilities were evaluated as negative in certain regions of the
parameter space. Specifically, the model predictions of $P(0 < S_{i}\leq s,Y_{i}>y \mid Z_i)$ or $P(0 < S_{i}\leq s,Y_{i}\leq y \mid Z_i)$
were negative when the prediction of $P(S_{i} > 0,Y_{i}>y \mid Z_i)$ was smaller
than the prediction of $P(S_{i}> s,Y_{i}>y \mid Z_i)$ or the prediction of $P(S_{i} > 0,Y_{i}\leq y \mid Z_i)$ was
smaller than the prediction of $P(S_{i}> s,Y_{i}\leq y \mid Z_i)$, respectively. Using simulations, we confirmed
that this region of negative predicted probabilities occurs when the parameters are far from their true
values, suggesting it is unlikely to affect estimation.
However, to avoid computational problems, we adopted a smooth transformation
to keep the estimated probabilities non-negative when they
fell below a certain threshold, say $\tau$.
Let $p$ represent the estimated probability, and let $\epsilon$
be a small number such that $0<\epsilon<\tau$, where we set $\epsilon=\tau/2$
for implementation. We use the smooth transformation:
\begin{align*}
f(p) & =1(p\geq\tau)\cdot p+1(p<\tau)\cdot\left\{ (\tau-\epsilon)\cdot\tanh\left\{ \dfrac{p-\tau}{\tau-\epsilon}\right\} +\tau\right\} \\
f^{\prime}(p) & =1(p\geq\tau)+1(p<\tau)\cdot\left\{ 1-\tanh^{2}\left\{ \dfrac{p-\tau}{\tau-\epsilon}\right\} \right\}
\end{align*}
We adopt this transformation because it has a continuous derivative
at $p=\tau$, which is used for gradient calculation. Therefore, the
Step 3 is implemented in practice as:
\begin{align*}
\hat{\rho}_{sy}= & \operatorname*{arg\,max}_{\rho_{sy}}\dfrac{1}{n}\sum_{i=1}^{n}D_{i}\Big[\bar{J}_{i}^{s}\bar{I}_{i}^{y}\log f\left(\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},X_{i}^{\prime}\hat{\nu}_{y};g(Z_{i}^{\prime}\rho_{sy}))\right)\\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; +\bar{J}_{i}^{s}I_{i}^{y}\log f\left(\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},-X_{i}^{\prime}\hat{\nu}_{y};-g(Z_{i}^{\prime}\rho_{sy}))\right)\\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+J_{i}^{s}\bar{I}_{i}^{y}\log f\left(\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},X_{i}^{\prime}\hat{\nu}_{y};g(X_{i}^{\prime}\hat{\rho}_{0y}))-\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},X_{i}^{\prime}\hat{\nu}_{y};g(Z_{i}^{\prime}\rho_{sy}))\right)\\
& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;+J_{i}^{s}I_{i}^{y}\log f\left(\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{0},-X_{i}^{\prime}\hat{\nu}_{y};-g(X_{i}^{\prime}\hat{\rho}_{0y}))-\Phi_{2}(Z_{i}^{\prime}\hat{\mu}_{s},-X_{i}^{\prime}\hat{\nu}_{y};-g(Z_{i}^{\prime}\rho_{sy}))\right)\Big].
\end{align*}
\end{remark}
\subsection{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{app:proof} 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{app:expressions}. We make the following assumption to derive the asymptotic distribution of $\hat \rho_{sy}$:
\begin{assumption}[Model and Sampling]\label{ass:fclt} (1) Random sampling: $\{(S^*_i, Y^*_i,Z_i)\}_{i=1}^n$ is a sequence of independent and identically distributed copies of $(S^*,Y^*,Z)$. We observe $S = \max(S^*,0)$ and $Y = Y^*$ if $S > 0$.
(2) Model: the joint distribution of $(S^*,Y^*)$ conditional on $Z$ follows the BDR model \eqref{eq:drmodel} with the exclusion restrictions. (3) The support of $Z$, $\mathcal{Z}$, is a compact set. (4) The sets $\mathcal{Y}$ and $\mathcal{S}$
are bounded intervals, which are compact subsets of the supports of $Y$ and $S$.
The conditional density function $f_{S^* \mid Z}(s \mid z)$ exists, is uniformly bounded above, and is uniformly continuous in $(s,z)$ on $\mathcal{S} \times \mathcal{Z}_1$, where $\mathcal{Z}_1$ is the support of $Z$ conditional on $S>0$, and the conditional density function $f_{Y^* \mid X}(y \mid x)$ exists, is uniformly bounded above, and is uniformly continuous in $(y,z)$ on $\mathcal{Y}\times \mathcal{X}_1$, where $\mathcal{X}_1$ is the support of $X$ conditional on $S>0$. (5) Identification and non-degeneracy: the equations $\mathrm{E}[\partial_{\mu_s} L_1(\tilde \mu_s)] = 0$, $\mathrm{E}[\partial_{\theta_y} L_2(\tilde \theta_y, \tilde \mu_0)] = 0$ and $\mathrm{E}[\partial_{\rho_{sy}} L_3(\tilde \rho_{sy}, \tilde \eta_{sy})] = 0$ posses a unique solution at $(\tilde \mu_s, \tilde \theta_y, \tilde \rho_{sy}) = (\mu_s, \theta_y, \rho_{sy})$ that lies in the interior of a compact set $\mathcal{M} \times \Theta \times \mathcal{R} \subset \mathbb{R}^{d_{\mu} + d_{\theta} + d_{\rho}}$ for all $y \in \mathcal{Y}$ and $s \in \mathcal{S}$; and the minimum eigenvalues of the matrices $H_{1s}$, $H_{2y}$, $H_{3sy}$ and $\Sigma_{\rho}(s,y)$ are bounded away from zero uniformly over $y \in \mathcal{Y}$ and $s \in \mathcal{S}$.
\end{assumption}
\begin{remark}[Assumption \ref{ass:fclt}(4)] The condition on the joint density function imposes smoothness on the coefficients of the BDR because
\ifnum\value{page}=1
\origfootnote{For a formal result, see Lemma \ref{lem:conz} in Appendix \ref{app:proof}.}
\else
\renewcommand{\arabic{footnote}}{\arabic{footnote}}
\origfootnote{For a formal result, see Lemma \ref{lem:conz} in Appendix \ref{app:proof}.}
\fi
$$
f_{S^*,Y^* \mid Z}(s,y \mid z) = \dfrac{\partial \Phi_{2}(-z^{\prime}\mu(s),-z^{\prime}\nu(y);g(z^{\prime}\rho(s,y)))}{\partial s \partial y}.
$$
We focus on the case where both $Y$ and $S$ are continuous because it is the most challenging, but it is easy to accommodate cases where $Y$ and/or $S$ are discrete by suitably modifying Assumption \ref{ass:fclt}(4). For example if $Y$ is continuous and $S$ discrete, Assumption \ref{ass:fclt}(4) becomes: $\mathcal{Y}$ is a bounded interval, $\mathcal{S}$
is finite, and the conditional density function $f_{Y^* \mid X}(y \mid x)$ exists, is uniformly bounded above, and is uniformly continuous in $(y,x)$ on $\mathcal{Y}\times \mathcal{X}_1$, where $\mathcal{X}_1$ is the support of $X$ conditional on $S > 0$. The restrictions that $\mathcal{Y}$ and $\mathcal{S}$ be compact subsets of the supports is imposed to avoid far tail estimation. They can be dispensed by imposing structure on the BDR's parameters at the tails.
\end{remark}
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}$.
\begin{thm}[FCLT for $\hat \rho_{sy}$]\label{thm:fclt} Under Assumption \ref{ass:fclt},
\[
\sqrt{n}(\hat{\rho}_{sy}-\rho_{sy})=-H_{3sy}^{-1}\sqrt{n}\left(S_{3sy}+J_{3sy}(\hat{\eta}_{sy}-\eta_{sy})\right)+o_{P}(1)\leadsto Z_{\rho_{sy}} \text{ in } \ell^{\infty}(\mathcal{SY})^{d_{\rho_{sy}}},
\]
where $S_{3sy} := \partial_{\rho_{sy}}L_{3}(\rho_{sy},\eta_{sy})$,
$H_{3sy} := \mathrm{E} \left[\partial_{\rho_{sy}\rho_{sy}^{\prime}}L_{3}(\rho_{sy},\eta_{sy})\right]$,
$J_{3sy} := \mathrm{E} \left[\partial_{\rho_{sy}\eta_{sy}^{\prime}}L_{3}(\rho_{sy},\eta_{sy})\right]$, and $(s,y)\mapsto Z_{\rho_{sy}}$ is a zero-mean Gaussian process with
uniformly continuous sample paths and covariance function
$\Sigma_{\rho_{sy}\rho_{s'y'}} = \mathrm{E}[Z_{\rho_{sy}} Z_{\rho_{s'y'}}],$
for $(s,y),(s',y')\in\mathcal{SY}$. Also, $\hat{\Sigma}_{\rho}(s,y)$, the estimator of $\Sigma_{\rho}(s,y) := \Sigma_{\rho_{sy}\rho_{sy}}$ defined in \eqref{eq:se} is uniformly consistent, that is,
$$
\sup_{(s,y) \in \mathcal{SY}} \|\hat{\Sigma}_{\rho}(s,y) - \Sigma_{\rho}(s,y) \| \to_{\Pr} 0.
$$
The analytical expressions of $S_{3sy}$, $H_{3sy}$, $J_{3sy}$ and $\Sigma_{\rho}(s,y)$ are given in Appendix \ref{app:expressions}.
\end{thm}
\medskip
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 \eqref{eq:se}.
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{thm:fclt}.
In practice, it is convenient
to estimate the critical value using resampling methods. The multiplier
bootstrap method of \citet{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.
\medskip
\begin{algorithm}[Multiplier Bootstrap]\label{alg:mboot}
For $b\in1,\ldots,B$ and a finite grid $\overline{\mathcal{SY}}_0$ covering $\mathcal{SY}_0$,
\begin{enumerate}
\item Draw bootstrap multipliers $\{\omega_{i}^{b}:1\leq i\leq n\}$
independently from the data and normalized them to have zero mean,
$\omega_{i}^{b}=\tilde{\omega}_{i}^{b}-\sum_{i=1}^{n}\tilde{\omega}_{i}^{b}/n,\ \ \tilde{\omega}_{i}^{b}\sim\text{ i.i.d. }\mathcal{N}(0,1)$.
\item Obtain a bootstrap draw of the estimator of $\rho_{sy}$ as
\[
\hat{\rho}_{sy}^{b}=\hat{\rho}_{sy}+n^{-1}\sum_{i=1}^{n}\omega_{i}^{b} \hat \psi_{3sy,\hat \xi}(W_i),
\]
where $\hat \psi_{3sy,\hat \xi}(W_i)$ is the estimator of the influence
function of $\hat{\rho}_{sy}$ given in \eqref{eq:influence}.
\item Construct bootstrap realization of maximal t-statistic
$$t_{\overline{\mathcal{SY}}_0}^{b}=\max_{(s,y)\in\overline{\mathcal{SY}}_0} \sqrt{n} \frac{|c'\hat{\rho}_{sy}^{b}-c'\widehat{\rho}_{sy}|}{\sqrt{c'\hat{\Sigma}_{\rho}(s,y)c}},$$
where $\hat{\Sigma}_{\rho}(s,y)$ be the estimator of $\Sigma_{\rho}(s,y)$ given in \eqref{eq:se}.
\end{enumerate}
Compute the critical value $cv(p)$ as the simulation $p$-quantile
of $t_{\overline{\mathcal{SY}}_0}^{b}$, that is, $$cv(p)=p-\text{quantile of }\left\{t_{\overline{\mathcal{SY}}_0}^{b}:1\leq b\leq B\right\}.$$
\end{algorithm}
We make the following assumption about the multipliers:
\begin{assumption}[Multiplier Bootstrap]\label{ass:mb} The multipliers $(\omega_{1},...,\omega_{n})$ are i.i.d. draws from a
random variable $\omega \sim \mathcal{N}(0,1)$, and are
independent of $\{(S^*_i, Y^*_i,Z_i)\}_{i=1}^n$ for all $n$.
\end{assumption}
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 \cite{van1996weak}.
\ifnum\value{page}=1
\origfootnote{We recall this definition in Appendix \ref{app:proof}.}
\else
\renewcommand{\arabic{footnote}}{\arabic{footnote}}
\origfootnote{We recall this definition in Appendix \ref{app:proof}.}
\fi
\begin{thm}[Multiplier Bootstrap FCLT for $\hat \rho_{sy}$]\label{thm:fclt_boot}
Under Assumptions \ref{ass:fclt} and \ref{ass:mb},
\[
\sqrt{n}(\hat{\rho}^b_{sy}-\hat \rho_{sy}) \leadsto_{\Pr} Z_{\rho_{sy}} \text{ in } \ell^{\infty}(\mathcal{SY})^{d_{\rho_{sy}}},
\]
where $Z_{\rho_{sy}}$ is the same Gaussian process as in Theorem \ref{thm:fclt}. In particular,
$$
\Pr\left[c^{\prime}\rho_{sy}\in \textrm{CB}_{p}(c^{\prime}\rho_{sy}),\forall(s,y)\in\overline{\mathcal{SY}}_0\right]\to p \text{ as } n \to \infty,
$$
where $\textrm{CB}_{p}(c^{\prime}\rho_{sy})$ is the confidence band obtained in Algorithm \ref{alg:mboot}.
\end{thm}
\begin{remark}[Estimation and Inference on Functionals] Estimation and inference on functionals is performed using the plug-in rule and multiplier bootstrap. FCLT and bootstrap multiplier FCLT can be derived using the functional delta method. We refer to Corollaries D.1 and D.2 of CFL for precise statements of these results.
\end{remark}
\section{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.
\subsection{Data}
We use the same dataset as CFL, which contains repeated cross-sectional
observations spanning from 1978 to 2013 \citep{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 \citet{gosling2000changing}, \citet{blundell2003interpreting},
and \citet{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.
\bigskip
\begin{table}[h]
\centering \caption{Descriptive Statistics of the Sample}
\begin{tabular}{l*{6}{c}}
\hline\hline &\multicolumn{2}{c}{\textbf{Full}} &\multicolumn{2}{c}{\textbf{Male}} &\multicolumn{2}{c}{\textbf{Female}} \\
& All & Employed & All & Employed & All & Employed \\
\hline Log Hourly Wage & & 2.38 & & 2.54 & & 2.21 \\
\hline
Employed & 0.74 & & 0.83 & & 0.66 & \\
Weekly work hours & 26.9 & 36.5 & 35.5 & 42.9 & 19.6 & 29.7 \\
Part-time & 0.19 & 0.26 & 0.02 & 0.05 & 0.18 & 0.50 \\
Full-time & 0.34 & 0.46 & 0.19 & 0.50 & 0.15 & 0.41 \\
Overtime & 0.20 & 0.28 & 0.17 & 0.45 & 0.03 & 0.09 \\
\hline
\multicolumn{7}{l}{Ceased School at} \\ \hspace{3mm} $\leq$ 15 & 0.33 & 0.30 & 0.33 & 0.31 & 0.33 & 0.29 \\ \hspace{3mm}16 & 0.31 & 0.30 & 0.32 & 0.32 & 0.30 & 0.29 \\ \hspace{3mm}17-18 & 0.18 & 0.19 & 0.16 & 0.17 & 0.20 & 0.22 \\ \hspace{3mm}19-20 & 0.04 & 0.05 & 0.04 & 0.04 & 0.04 & 0.05 \\ \hspace{3mm}21-22 & 0.09 & 0.11 & 0.09 & 0.10 & 0.09 & 0.12 \\ \hspace{3mm}$\geq$23 & 0.05 & 0.05 & 0.06 & 0.06 & 0.04 & 0.04 \\ Age & 40.13 & 39.84 & 40.22 & 39.76 & 40.06 & 39.92 \\ Married & 0.76 & 0.79 & 0.78 & 0.81 & 0.75 & 0.76 \\ Benefit Income & 5.44 & 5.50 & 5.25 & 5.29 & 5.60 & 5.73 \\ \hline Observations & 258,900 & 190,765 & 119,396 & 98,764 & 139,504 & 92,001 \\ \hline\hline
\end{tabular}
\end{table}
The instrumental variable, $Z_{1}$, is the potential out-of-work
income benefit interacted with the marital status indicator, as previously
used in \citet{blundell2003interpreting}, \citet{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 \citet{blundell2003interpreting} and \citet{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.
\subsection{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
\begin{itemize}
\item Specification 1: a constant $z^{\prime}\rho(s,y)=\rho(s,y)$.
\item Specification 2: a constant and the marital status indicator $z^{\prime}\rho(s,y)=\rho_{M}(s,y)m+\rho_{S}(s,y)(1-m)$.
where $m=1$ if the individual is married.
\item Specification 3: a constant and a linear trend on the year of the
survey $z^{\prime}\rho(s,y)=\rho_{0}(s,y)+\rho_{1}(s,y)t$ where $t=\text{year}-1978$.
\item Specification 4: a constant and a linear trend on the year of the
survey interacted with the marital status indicator $z^{\prime}\rho(s,y)=\rho_{0}(s,y)+\rho_{1}(s,y)t+\rho_{2}(s,y)m+\rho_{3}(s,y)t\times m.$
\end{itemize}
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.
\bigskip
\begin{figure} [ht!]
\centering
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Hist_male.pdf}
\end{subfigure}
\begin{subfigure}{0.49\textwidth}
\includegraphics[width=\textwidth]{Hist_female.pdf}
\end{subfigure}
\caption{Histograms of work hours by gender} \label{fig2}
\end{figure}
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 (\citet{mulligan2008selection}, \citet{maasoumi2019gender}, \citet{blundell2021wages}). Moreover, as shown in Figure \ref{fig2}, which presents histograms of work hours by gender, our data exhibits clear bunching at the thresholds of 34 and 40 hours. Similarly, \citet{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.
\subsection{Main results}
\subsubsection{Employment sorting}
Figure \ref{fig3} 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.
\begin{figure} [ht!]
\centering
\includegraphics[width=\textwidth]{figure3.pdf}
\caption{Estimates and 95\% confidence bands for the selection sorting function using Specification 1} \label{fig3}
\end{figure}
Figures \ref{fig4} and \ref{fig5} 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{fig3} shows that men exhibit positive sorting across the wage distribution, whereas women show positive sorting only at the lower wage quantiles. Figure \ref{fig4} 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.
\begin{figure} [H]
\centering
\includegraphics[width=.85\textwidth, height=9cm]{figure4.pdf}
\caption{Estimates and 95\% confidence bands for the selection sorting into full-time using Specification 2} \label{fig4}
\end{figure}
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 \citet{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. \citet{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{fig3} shows that both genders
exhibit a negative selection effect, which is more pronounced
at the top of the wage distribution. Figure \ref{fig5} 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.
\begin{figure} [H]
\centering
\includegraphics[width=.85\textwidth, height=9cm]{figure5.pdf}
\caption{Estimates and 95\% confidence bands for the selection sorting into overtime using Specification 2} \label{fig5}
\end{figure}
Figure \ref{fig6} 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 \citet{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{fig9} and \ref{fig10} in the Appendix \ref{app:spec4} indicate that these sorting
trends are primarily driven by single men and married women, similar
to the findings in Specifications 1 and 2.
\begin{figure} [H]
\centering
\includegraphics[width=.85\textwidth, height=9cm]{figure6.pdf}
\caption{Selection sorting function using Specification 3} \label{fig6}
\end{figure}
\subsubsection{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:
\begin{align*}
\underset{\text{Total}}{\underbrace{Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}}} & =\underset{\text{Wage Structure}}{\underbrace{[Q_{Y\langle1,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]}]}}+\underset{\text{Selection Sorting}}{\underbrace{[Q_{Y\langle0,1,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,1,1\rangle}^{(\underline{s},\overline{s}]}]}}\\
& +\underset{\text{Selection Structure}}{\underbrace{[Q_{Y\langle0,0,1,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,0,1\rangle}^{(\underline{s},\overline{s}]}]}}+\underset{\text{Composition}}{\underbrace{[Q_{Y\langle0,0,0,1\rangle}^{(\underline{s},\overline{s}]}-Q_{Y\langle0,0,0,0\rangle}^{(\underline{s},\overline{s}]}]}}.
\end{align*}
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 (\citet{kitagawa1955components};
\citet{oaxaca1973male}; \citet{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
\begin{footnotesize}
\[
F_{Y\langle t,j,r,k\rangle}^{(\underline{s},\overline{s}]}(y)=\frac{\int[\Phi_{2}(-z^{\prime}\mu_{r}(\overline{s}),-x^{\prime}\nu_{t}(y);g(z^{\prime}\rho_{j}(\overline{s},y)))-\Phi_{2}(-z^{\prime}\mu_{r}(\underline{s}),-x^{\prime}\nu_{t}(y);g(z^{\prime}\rho_{j}(\underline{s},y)))]dF_{Z_{k}}(z)}{\int[\Phi(-z^{\prime}\mu_{r}(\overline{s}))-\Phi(-z^{\prime}\mu_{r}(\underline{s}))]dF_{Z_{k}}(z)},
\]
\end{footnotesize}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
\renewcommand{\arabic{footnote}}{\arabic{footnote}}
\origfootnote{The operator $Q_\tau$ monotonically arrange the estimator of the distribution before applying the (left) inverse.}
\fi
Figures \ref{fig7} and \ref{fig8} 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.
\begin{figure} [H]
\centering
\includegraphics[width=.85\textwidth]{decomp_ft.pdf}
\caption{Estimates and 95\% confidence bands for the quantiles of observed wages and decomposition between full-time working men and women in Specification 2} \label{fig7}
\end{figure}
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
\renewcommand{\arabic{footnote}}{\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.
\begin{figure} [H]
\centering
\includegraphics[width=.85\textwidth]{decomp_ot.pdf}
\caption{Estimates and 95\% confidence bands for the quantiles of observed wages and decomposition between overtime working men and women in Specification 2} \label{fig8}
\end{figure}
\subsubsection{Work hours decomposition}
Lastly, we decompose the difference in the distribution functions of observed work hours between genders using counterfactual distributions. Let
\begin{align*}
F_{S<r,k>}(s) = \int \Phi(-z^{\prime}\mu_{r}(s)) \mathrm{d} F_{Z_{k}}(z), \quad s \geq 0,
\end{align*}
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
\begin{align*}
\underset{\text{Total}}{\underbrace{F_{S<0,0>}-F_{S<1,1>}}}=\underset{\text{Structure}}{\underbrace{[F_{S<0,0>}-F_{S<1,0>}]}}+\underset{\text{Composition}}{\underbrace{[F_{S<1,0>}-F_{S<1,1>}]}}.
\end{align*}
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
\renewcommand{\arabic{footnote}}{\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.
\begin{figure}[H]
\centering
\includegraphics[width=.85\textwidth]{hours_decomp_new.pdf}
\caption{Estimates and 95\% confidence bands for the distribution of work hours and decomposition between men and women} \label{fig9}
\end{figure}
\section{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.
\bibliographystyle{ecta}
\bibliography{Reference}