Point-Identifying Semiparametric Sample Selection Models with No Excluded Variable

Point-Identifying Semiparametric Sample Selection Models with No Excluded Variable

Introduction

Sample selection poses a fundamental challenge in empirical economics, threatening the validity of research findings. When workers self-select into employment or patients choose whether to seek medical care, the resulting datasets systematically exclude critical segments of the population, distorting our understanding of economic relationships and leading to misguided policy prescriptions. Even carefully designed experiments are vulnerable, as systematic attrition introduces selection bias that undermines randomization. While economists have proposed various solutions, they often rely on strong assumptions or yield uninformative bounds. In particular, nonparametric and semiparametric selection models have long been thought to require exclusion restrictions, limiting their practical applicability. This paper challenges that conventional wisdom by introducing a class of semiparametric selection models that achieve point identification without exclusion restrictions. We further develop computationally tractable two-step plug-in estimators that can be readily implemented using standard statistical software.

heckman1974shadow, heckman1979sample pioneered correction methods for selection bias using a parametric model that assumes linearity in both selection and outcome equations, along with joint normality of the error terms:

where $Y^*$ is the latent outcome, $X$ and $Z$ are row vectors of exogenous covariates, $V$ and $\varepsilon$ are mean-zero unobserved heterogeneity terms that are joint normally distributed and independent of $(X,Z)$, with $Var(\varepsilon)$ normalized to $1$. Conditional on $X = x$, $Z = z$, and $D = 1$, the mean of the observed outcome is given by: $$E[Y|x,z,D=1] = \alpha + x\beta + \sigma_{V\varepsilon}\phi(z\gamma)/\Phi(z\gamma),$$ where $\sigma_{V\varepsilon}= Cov(V, \varepsilon),$ and $\phi(\cdot)$ and $\Phi(\cdot)$ denote the standard normal probability density function (p.d.f.) and cumulative distribution function (c.d.f.) respectively.

Although Heckman’s model can identify $(\alpha, \beta)$ when $X=Z$ due to the known functional form of selection bias, it is generally recommended to include at least one variable in $Z$ that is excluded from $X$ to strengthen identification. Without such an exclusion restriction, the numerical performance of both the two-step and maximum likelihood estimators can be poor, as highlighted in the debates in duan1984choosing, manning1987monte, and hay1984let.

By relaxing Heckman’s joint normality assumption, econometricians developed semiparametric approaches chamberlain1986asymptotic, ahn1993semiparametric, newey2009two. Later das2003nonparametric explored fully nonparametric selection models:

In both semiparametric and nonparametric models, the exclusion restriction is widely regarded as essential for identification.\footnote{lee2009bounds stated that “standard parametric or semiparametric methods for correcting for sample selection require exclusion restrictions that have little justification in this case” (p. 1072). Similarly, honore2020selection claimed that an exclusion restriction is the key identifying assumption in semiparametric selection models.} However, finding an excluded variable is often infeasible in empirical applications. Motivated by this challenge, lee2009bounds proposed a nonparametric bounds approach that does not require exclusion restrictions. By exploiting the selection monotonicity assumption, under which individuals who are observed without treatment would also be observed with treatment, Lee introduced a trimming procedure to adjust for missing data due to selection. His approach (henceforth Lee bounds) is intuitive and easily implemented, making it widely used in empirical studies, particularly in experiments where subjects tend to drop out.

Lee bounds, however, are often too wide to yield meaningful economic insights and its ability to incorporate covariate information is limited.\footnote{He proposes that weakly tighter bounds than the unconditional bounds can be obtained by averaging the group-specific effects, weighted by covariate density. However, in practice, discretization is necessary for continuous covariates, and handling a large number of covariates is often infeasible. semenova2023generalized generalizes Lee’s approach to a high-dimensional setting.} honore2020selection (henceforth HH) later demonstrated that $\beta$ is partially identified in (ref) without distributional assumptions on $(V, \varepsilon)$ and in the absence of exclusion restrictions. The HH model, serves as a semiparametric alternative to Lee's, provides tighter bounds than Lee bounds, as it imposes additional structural assumptions. Despite the growing popularity of partially identifying models, estimating identified sets and conducting inference remain challenging, particularly when the identified set is characterized by a large number of moment inequalities or the distributions of unobserved heterogeneity lack parametric restrictions.

Motivated by these challenges, this paper investigates semiparametric selection models that achieve point identification of $\beta$ without any excluded variable. Specifically, we relax the linear selection assumption in HH’s model, providing a middle ground between Lee’s and HH’s approaches. Unlike Lee’s framework, our method does not impose selection monotonicity, meaning it is not nested within Lee's. We challenge the prevailing belief that an exclusion restriction is necessary for semiparametric selection models. Our approach establishes point identification of $\beta$ under minimal assumptions without requiring scale normalization or identification at infinity, when there exists at least one continuous variable in $X$.

Our identification strategy leverages the nonlinearity of the conditional selection probability $p_0(X):= E[D|X]$, which is nonparametrically identified and easily verifiable in practice. Consequently, $\beta$ can be estimated via a partial linear regression, plugging in estimates of $p_0(\cdot)$ in the nonparametric components approximated by sieves, without needing to estimate both $p_0(X)$ and $E[Y|X]$ nonparametrically. Unlike many existing methods, we do not assume unobserved heterogeneity is independent of regressors, allowing for heteroskedasticity. We demonstrate that our estimators for $\beta$ are $\sqrt{n}$-consistent, semiparametrically efficient, asymptotically normal, and computationally scalable. As we maintain the linearity of the outcome equation, incorporating a large set of covariates is straightforward. Our proposed two-step semiparametric estimator performs exceptionally well in simulations and an empirical application on gender and racial wage gaps in the United State (US).

We are not the first to consider nonlinearity in the selection process as a means to identify the linear index parameters in the latent outcome equation. For instance, ahn1993semiparametric and newey1993efficiency speculated that nonlinearity can yield nonzero semiparametric efficiency bounds, as the nonlinear terms in the selection equation may act as excluded variables. In this paper, we formally establish the conditions which suffice to identify the model parameters. More recently, escanciano2016identification introduced a more general model, assuming $E[Y|X] = F_0(X\beta_0, p_0(X))$, and demonstrated that $\beta_0$ can be point identified up to scale if $X$ contains at least two continuous variables and $p_0(X)$ is nonlinear in $X$. However, they normalize one of the continuous variables' coefficients to $1$, implying that full identification is not achieved. Unlike our proposed estimator, their approach requires nonparametric estimation of both $p_0(X)$ and $E[Y|X]$ even with linearity in the outcome equation, adding to its complexity.

pan2022semiparametric propose an integrated nonlinear least squares estimator chen2010integrated, chen2010non, chen2011semiparametric, chen2012semiparametric for escanciano2016identification's model, eliminating the need for scale normalization. However, they do not provide an identification argument, as they assume $\beta_0$ is already identified. Moreover, their estimation procedure relies on multiple layers of kernel regression, each requiring the selection of multiple tuning parameters, along with numerical optimization of an integrated criterion function. This results in a computational burden that scales exponentially with the sample size.

This paper is organized as follows. Section (ref) introduces our semiparametric selection model and establishes identification. Section (ref) presents our estimators and derives their asymptotic properties. Section (ref) evaluates finite-sample performance of our estimator via simulations. Section (ref) applies our method to estimating gender and racial wage disparities in the US. Section (ref) concludes.

The semiparametric selection model

We consider a sample selection model where the selection procedure is left unspecified:

Let $p_0(x) := P[D=1|X=x]$ represent the conditional selection probability given $X=x$. We impose the following assumption to identify $\beta_0$ without invoking an exclusion restriction.

The above assumption requires a continuously distributed variable in $X$, the smoothness of $p_0(\cdot)$ and $\lambda_0(\cdot).$ It also restricts the selection bias to only depend on the selection probability in a nonparametric form, $\lambda_0(\cdot).$ Additionally, it rules out a flat region in $p_0(\cdot).$ Our model is implied by a special case where

Suppose $\varepsilon$ is continuously distributed with the c.d.f. $F_\varepsilon(\cdot).$ We can normalize the selection process to $D = \mathbbm{1}[p_0(X) \ge U]$, where $p_0(X) = F_\varepsilon(h_0(X))$ and $U = F_\varepsilon(\varepsilon) \sim Unif(0,1)$. Unlike this special case, we do not impose the stochastic independence between $X$ and $(V,\varepsilon).$ For some $\beta$ and $\lambda(\cdot)$, define $l(x) := (\beta_0 - \beta) x$ and $b(p) := \lambda_0(p) - \lambda(p),$ both of which are deviations of $\beta$ and $\lambda(p)$ from the truth. For any observationally equivalent $\beta$ and $\lambda$ such that $E[Y|X, D=1] = X\beta + \lambda(p_0(X)),$ $l(x) + b(p) = 0$ identically.

Under Assumption (ref), $\beta_0$ can be point identified. We cannot separately identify the intercept $\alpha_0$ from the selection bias.\footnote{For point identification of the intercept $\alpha_0$, see heckman1990varieties and andrews1998semiparametric. Unlike our paper, both papers relies on the“identification at infinity” argument and an exclusion restriction.} Firstly, we consider the simplest case in which $X$ consists of only one continuous variable.

This proposition shows that more than nonlinearity is required in the selection process, as it rules out monotonicity of $p_0(\cdot)$. Consider $h_0(\cdot)$ in the special case described above. If $h_0(\cdot) = X + 0.5 X^2,$ $\beta_0$ is not point identified. When $X$ is binary, $p_0(X) = \gamma_0 +\gamma_1 X$ is fully nonparametric and therefore the parameters in the latent outcome equation are not point identified as shown in HH.

Now we consider more general cases where $X$ is multidimensional. Suppose two elements of $X$, $X_k$ and $X_j$, are continuously distributed. The following proposition shows that the model is point identifying.

In general, the marginal effect of $x_k$ on $p$ is not proportional to that of $x_j.$ Identification fails in the special case (ref) when $p_0(X) = F_\varepsilon(X\gamma)$ because $\frac{\partial p_0/\partial x_k}{\partial p_0/\partial x_j} =\frac{ \gamma_kf_\varepsilon(X\gamma)}{\gamma_jf_\varepsilon(X\gamma)} = \gamma_k/\gamma_j$ where $F_\varepsilon(\cdot)$ and $f_\varepsilon(\cdot)$ are the c.d.f. and p.d.f. of $\varepsilon.$ Therefore, the nonlinearity of $p_0(\cdot)$ enables us to identify the model parameters.

Lastly, we consider the case where $X_k$ is the only continuous variable in $X$. Without loss of generality, let $X_j$ be a binary variable for all $j \neq k$, as any discrete variable can be equivalently expressed as a set of dummy variables.

Unless the required change in $x_k$ to achieve $p''$ with $x_j = 0$ from $p'$ with $x_j = 1$ is constant for all values of $x_k',$ the model parameters are identified. This proposition rules out the case in which $p_0(X) = F_\varepsilon(X\gamma)$ because $p'' = F_\varepsilon(x''\gamma) = F_\varepsilon(x'\gamma - \gamma_j) = F_\varepsilon\left(x'\gamma + (x_k''' - x_k')\gamma_k\right)$ and hence $x_k'-x_k''' = \gamma_j/\gamma_k$ for all $x_k'.$

The identification results in this section show that our semiparametric selection model can point identify the model parameters without an excluded variable as long as there is at least one continuously distributed covariate. In applied economic studies, continuous variables such as age, income, and price are not uncommon. Hence, our semiparametric model can be quite generally applicable to many modern data sets. Furthermore, it is natural that the true selection process exhibit some degree of nonlinearity. As the conditional selection probability is nonparametrically identified, it is simple to check the nonlinearity in the selection probability. When there exists strong empirical evidence of selection nonlinearity, our model becomes a strong alternative to the HH's model, as our model is more robust and point-identifying.

The estimators

Given the identification results, we now propose a class of tractable two-step plug-in semiparametric estimators. We first begin by demonstrating that these estimators are consistent, semiparametrically efficient, and asymptotically normal when the individual selection probability, $p_{i}=p_{0}\left(X_{i}\right)=\mathbb{E}\left[D_{i}|X_{i}\right]$, is observed. Subsequently, we prove that replacing $p_i$ with $\hat{p}_i$, a consistent estimator of $p_i$, does not affect the asymptotic behavior of our estimators when $\hat{p}_i$ converges to $p_i$ at a sufficiently fast rate as $n \to \infty$. Then we describe how the estimator can be practically implemented.

When $p_i$ is known

Suppose we have an i.i.d. sample, $\{W_i\}_{i=1}^n$, where $W_i:=(Y_i, D_i, X_i, p_i).$ As $p_i$ is observed, the model parameters $\theta_{0}=\left(\beta_{0},\lambda_{0} \right)\in \Theta := \mathcal{B}\times \Lambda$ can be estimated by the least squares (LS) procedure: \[ \tilde{\theta}_{n}=\left(\tilde{\beta}_{n},\tilde{\lambda}_{n}\right) =\operatorname*{arg\,min}_{\left(\beta,\lambda\right)\in\mathcal{B}\times\Lambda} \frac{1}{n}\sum_{i=1}^{n} D_{i}\left(Y_{i}-X_{i}'\beta-\lambda\left(p_{i}\right)\right)^{2}. \] As $\lambda$ is infinite-dimensional, we approximate the unknown function $\lambda\in\Lambda$ by sieves, $\lambda_{n}\in\Lambda_{n},$ where $\Lambda_{n}$ is an approximating function space (such as polynomials, trigonometric polynomials, splines, and orthogonal wavelets) that becomes dense in $\Lambda$ as $n\rightarrow\infty$.

Let $\Lambda_{n}=\left\{\lambda_n\left(\cdot\right) =R_{K\left(n\right)}\left(\cdot\right)'\gamma:\gamma\in\mathbb{R}^{K\left(n \right)}\right\}$, where $R_{K\left(n\right)}=\left[r_{1}\left(\cdot\right), \ldots,r_{K\left(n\right)}\left(\cdot\right)\right]'$ denote a vector of basis functions. Let $\Theta_{n}=\mathcal{B}\times \Lambda_{n}$ be the sieve space for $\theta=\left(\beta,\lambda\left( \cdot\right)\right)$ and $\bar{K}\left(n\right)=\dim(\beta)+K\left(n\right)$. For notational simplicity, we omit the subscript $K\left(n \right)$ and write $R_{K\left(n\right)}\left(\cdot\right)=R\left(\cdot \right)$. Define $$D=diag\left(D_{1},\ldots,D_{n}\right),\quad X=\left[X_{1},\ldots,X_{n}\right]',\quad y=\left(Y_{1},\ldots,Y_{n}\right)', \quad p_0=(p_1, \cdots, p_n)',$$ For an arbitrary vector $\tilde{p}:= (\tilde{p}_1, \cdots, \tilde{p}_n)$ where $\tilde{p}_i \in [0,1],$ define $R\left(\tilde{p}\right)=D\left[R\left(\tilde{p}_{1}\right),\ldots,R\left(\tilde{p}_{n}\right)\right]'$, and $Q\left(\tilde{p}\right)=R\left(\tilde{p}\right)\left(R\left(\tilde{p}\right)'R\left(\tilde{p}\right)\right)^{-1}R\left(\tilde{p}\right)'$. Then the LS estimator of $\beta$ is written as: \[ \tilde{\beta}_{n} =\left(\left(DX\right)'\left(I-Q\left(p_{0}\right)\right)\left(DX\right)/n\right)^{-1} \left(DX\right)'\left(I-Q\left(p_{0}\right)\right)\left(Dy\right)/n. \]

We establish asymptotic properties of our sieve LS (SLS) estimator using the results in chen2007. Let $\Theta$ be equipped with a norm $\lVert\theta\rVert_{s}=\left|\beta\right|_{e}+\lVert \lambda\rVert_{\infty}$, where $\left|\cdot\right|_{e}$ denotes the Euclidean norm and $\lVert \lambda\rVert_{\infty}=\sup_{p\in\left[0,1\right]}\left|\lambda\left(p\right)\right|$ is the supremum norm. We introduce the H\"{o}lder class of functions. Let $\left[m\right]$ be the largest nonnegative integer such that $\left[m\right] <m$. A real-valued function $\lambda$ on $\left[0,1\right]$ is said to be in the H\"{o}lder space $\Lambda^{m}\left(\left[0,1\right]\right)$ if it is $\left[m\right]$ times continuously differentiable on $\left[0,1\right]$ and \[ \max_{\ell\leq\left[m\right]}\sup_{p} \left|\frac{\partial^{\ell}\lambda\left(p\right)}{\partial p^{\ell}}\right| +\sup_{p,p'} \left|\frac{\partial^{\left[m\right]}\lambda\left(p\right)}{\partial p^{\left[m\right]}} -\frac{\partial^{\left[m\right]}\lambda\left(p'\right)}{\partial p^{\left[m\right]}}\right| /\left|p-p'\right|^{m-\left[m\right]} \] is finite. We impose the following conditions to derive asymptotic properties of the SLS estimator.

Under these regularity conditions, the consistency of the estimator is obtained by Proposition 3.3 of chen2007 as $\lVert\tilde{\theta}_{n}-\theta_{0}\rVert= O_{p}\left(n^{-m/\left(2m+1\right)}\right).$

We now show that the parametric components of the SLS estimator, $\tilde{\beta}_{n},$ is asymptotically normal. Let $\tilde{X}_i=D_iX_i -\mathbb{E}\left[D_iX_i|p_i,D_i=1\right]$.

Assumption (ref)(i) is satisfied when $\beta_0$ and $\lambda_0$ are identified. Applying Proposition 4.5 of chen2007, we obtain the following asymptotic normality of $\tilde{\beta}_n$.

If the error term is homoskedastic i.e., $\sigma_{0}(X_i, D_i)$ is constant, the SLS estimator $\tilde{\beta}_n$ is semiparametrically efficient. When the error term exhibits heteroskedasticity, efficient estimation can be achieved through the sieve generalized least squares (SGLS) estimator. However, in applied economic studies, the standard practice is to report heteroskedasticity-robust standard errors rather than employing the GLS approach. Consequently, in our simulations and empirical applications, we proceed with robust standard errors calculated using the asymptotic variance formula reported in Proposition (ref).

When $p_i$ is replaced by $\hat{p}_i$

As $p_i$ is never observed in practice, $\tilde{\beta}_n$ is an infeasible estimator. Suppose $p_i$ is consistently estimated by an estimator $\hat{p}_i := \hat{p}_{n} \left(X_i\right)=p_{0}\left(X_i\right)+O_{p}(n^{-1/3})$. Define $\hat{p}= (\hat{p}_1, \cdots, \hat{p}_n)'.$ Replacing $p_0$ with $\hat{p}$ in $\tilde{\beta}_n,$ we yield the following feasible estimator: \[ \hat{\beta}_{n} =\left(\left(DX\right)'\left(I-Q\left(\hat{p}\right)\right) \left(DX\right)/n\right)^{-1}\left(DX\right)'\left(I-Q\left(\hat{p}\right)\right)\left(Dy\right)/n. \] We will show that $\hat{\beta}_n$ and $\tilde{\beta}_n$ have the same asymptotic distribution by leveraging the results in song2012, song2014. Define \[

\] Let $d_X$ denote $\dim(X)$ and define $H\left(a,b\right)=\left(b_{2}-a_{2}\right)^{-1}\left(b_{1}-a_{1}\right)$ where $a_{2}$ and $b_{2}$ are the right $d_X\times d_X$ subblocks, and $a_{1}$ and $b_{1}$ are the left $d_X\times1$ subblocks of $a$ and $b$. Then, we can write $\lVert \hat{\beta}_{n}-\tilde{\beta}_{n}\rVert$ as $\lVert H\left(\hat{a}_{n}\left(\hat{p}\right),\hat{b}_{n}\right) -H\left(\hat{a}_{n}\left(p_{0}\right),\hat{b}_{n}\right)\rVert$. From the continuously differentiability of $H$, we have \[ \lVert H\left(\hat{a}_{n}\left(\hat{p}\right),\hat{b}_{n}\right) -H\left(\hat{a}_{n}\left(p_{0}\right),\hat{b}_{n}\right)\rVert \leq C\lVert\hat{a}_{n}\left(\hat{p}\right)-\hat{a}_{n}\left(p_{0} \right)\rVert+o_{p}\left(\lVert\hat{a}_{n}\left(\hat{p}\right)- \hat{a}_{n}\left(p_{0}\right)\rVert\right). \] For the RHS of the above inequality, observe that

where $a_{0}\left(p\right):=\mathbb{E}\left[D_iX_i\mathbb{E}\left[Z_i|p_i,D_i=1\right]\right] $ with $Z_i=\left[Y_i;X_i\right]$.

As $\lVert\hat{p}-p_{0}\rVert=o_{p}\left(1\right)$, $A_{n}(\hat{p})$ and $A_{n}(p_0)$ in (ref) can be shown to be $O_{p}\left(n^{-1/2}\right).$ First, we observe that {

} We obtain the asymptotic linear representation of $\sqrt{n}A_n\left(\tilde{p}\right)$ as follows:

for some i.i.d. $\psi_{i}\left(\cdot\right)$ such that $\mathbb{E}[\psi_{i}\left(\tilde{p}\right)]=0,$ where $o_{p}\left(1\right)$ is uniform local around $p_{0}.$ Then, applying the maximal inequality yields the stochastic equicontinuity of $A_n(\cdot)$ as shown in andrews1994empirical. By Lemma B3 in song2014, we rewrite the first element of the RHS of (ref) to: \[ \frac{1}{n}\sum_{i=1}^{n}D_{i}\mathbb{E}\left[X_{i}|\tilde{p}_{i},D_{i}=1\right] \left(Z_{i}-\mathbb{E}\left[Z_{i}|\tilde{p}_{i},D_{i}=1\right]\right)+o_{p}\left(n^{-1/2}\right), \] uniformly over $\tilde{p}\in B\left(p_{0};c_{n}\right):=\left\{\tilde{p}:\lVert \tilde{p}-p_{0}\rVert<c_{n}\right\}$, which implies that $\sqrt{n}A_n\left(\tilde{p}\right)$ is, uniformly over $p\in B\left(p_{0};c_{n}\right)$, equal to \[

\]

From this uniform linear representation, we have $\sup_{\tilde{p}\in B\left(p_{0};c_{n}\right)}\left|\sqrt{n}A_n\left(\tilde{p}\right)\right|=O_{p}\left(1\right)$, which means that both $A_{n}(\hat{p})$ and $A_{n}(p_0)$ are $O_p(n^{-1/2}).$ Furthermore, we have

whose asymptotic variance, $\mathbb{E}\left[\left(\psi_{i}\left(\tilde{p}\right)-\psi_{i}\left(p_{0}\right) \right)^{2}\right],$ goes to 0 as $\tilde{p}\rightarrow p_{0}$ under minor regularity conditions for $\psi_{i}\left(\tilde{p}\right)$. This result also works when $\tilde{p} = \hat{p}$ as shown in Lemma (ref) provided in the appendix. Hence we conclude $A_n\left(\hat{p}\right)-A_n\left(p_0\right)=o_{p}\left(n^{-1/2}\right)$ as $\hat{p}\rightarrow p_{0}$.

Lastly, we show that the term $B_{n}$ in (ref) is also $o_p(n^{-1/2})$. We extend song2014's results, which are derived under the linear selection procedure, to the cases with nonlinear, possibly nonmonotone selection. Under regularity conditions, the function $a_0\left(\tilde{p}\right)$ is sufficiently smooth in $\tilde{p}$ around $p_{0}$ so that there exist constants $C>0$ and $\varepsilon\in\left(0,1/2\right]$ such that for each $\eta\in\left(0,\varepsilon \right]$,

The formal proof of (ref) is provided in Lemma (ref) in the appendix. Hence, $\lVert a_{0}\left(\hat{p}\right)-a_{0}\left(p_{0}\right)\rVert=O_{p}\left(\eta_{n}^{2}\right)$ if $\lVert\hat{p}-p_{0}\rVert\leq\eta_{n}$. As we consider $\hat{p}$ converging to $p_0$ at a cube-root rate, we obtain $\lVert a_{0}\left(\hat{p}\right)-a_{0}\left(p_{0}\right)\rVert =o_{p}\left(n^{-1/2}\right)$ by taking $\eta_{n}= n^{-1/3}\log n$. This concludes that $\lVert\hat{\beta}_{n}-\tilde{\beta}_{n}\rVert=o_{p}\left(n^{-1/2}\right),$ which implies that $\hat{\beta}_{n}$ and $\tilde{\beta}_{n}$ have the same asymptotic distribution.

Practical implementations

Given the asymptotic results, one can consider a wide range of estimators in the first stage estimation of $p_i =P[D_i=1|X_i].$ Any consistent nonparametric estimator $\hat{p}_n(\cdot)$ that converges to $p_0$ at a cube-root rate or faster can be employed. The convergence rate of the first stage estimation depends on the smoothness of $p_0$ and the number of continuous elements in $X_i,$ denoted as $d_c$. With a relatively low $d_c$, it is quite feasible that standard kernel or sieve estimators converge faster than the $n^{-1/3}$ rate. In a high-dimensional setting, a high level of smoothness for $p_0$ is necessary to ensure a sufficiently fast rate. Recall that $m$ denotes the H\"{o}lder smoothness of $p_0$. The rate condition is satisfied when $m > d_c$. Suppose that the true selection process is $D_i =\mathbbm{1}[g_k(X_i) \ge U_i]$ where $g_k(\cdot)$ is a k-th order polynomial of $X_i$ and $U_i$ is continuously distributed with the c.d.f. $F_U(\cdot)$ which belongs to a smooth parametric class. Consequently, $p_0(X_i) = F_U(g_k(X_i))$ and $p_0(\cdot)$ is infinitely continuously differentiable, implying $m = \infty.$ Therefore, if one would like to impose the selection procedure outlined above, standard kernel or sieve methods can be employed with a high $d_c$. Alternatively, additional structural restrictions can be imposed, such as additivity, i.e. $p_0(X_{1i}, X_{2i}) = p_{10}(X_{1i})+p_{20}(X_{2i}).$

In practical implementation, we propose the following simple two-step procedure.

In the simulations and empirical studies conducted in this paper, we employ the sieve maximum likelihood estimator in Step 1, using piecewise polynomial basis functions. In Step 2, we regress $Y_i$ on $X_i$ and piecewise polynomial transformations of $\hat{p}_i$ using ordinary least squares (OLS). The OLS standard errors of $\hat{\beta}$ provided in standard statistical programs such as Stata, R, and Matlab are asymptotically valid standard errors under homoskedasticity. In practice, researchers often desire heteroskedasticity-robust or cluster-robust standard errors. The same 'sandwich' formula can be used to compute robust standard errors. Therefore, the two-step procedure outlined here can be readily utilized in any statistical software without necessitating a new implementation package, which makes our proposal particularly attractive to applied researchers.

Simulations

In this section, we evaluate the finite sample performances of our semiparametric estimator using known data-generating processes (DGPs). For each DGP, we repeat 1,000 iterations in each of which we draw a Monte Carlo sample of size $n = 5,000.$ We first investigate the single-covariate case using the following DGP:

In this instance, $\beta_0$ is not separately identified from $E[V|X, D=1] = \lambda_0(p_0).$ The identification of $\beta_1$ is contingent upon the parameter values $\alpha = (\alpha_0, \alpha_1, \alpha_2, \alpha_3)$. This is because the conditional selection probability $p_0(X) = E[D|X]$ must not exhibit strict monotonicity. We employ the proposed two-step sieve-based approach to estimate the model. For alternative estimators, we consider the ordinary least squares (OLS) estimator conditional on selection ($D=1$) assuming random selection, which is commonly referred to as the two-part model (TPM), and the maximum likelihood estimator under the Heckman selection model (HSM), both of which are misspecified. We also compare our estimator to the oracle estimator, which incorporates the true functional form of $p_0(X)$ given the selection bias is expressed using the inverse Mills ratio.\footnote{In the oracle estimation, we first estimate $\alpha$ using probit regression of $D$ on $Z = (1, X_1, X_1^2, X_1^3)$. Subsequently, we employ $\lambda_0(\hat{p}_0(Z)) = \phi(Z\hat{\alpha})/\Phi(Z\hat{\alpha})$ to correct the selection bias.}

We consider two selection designs: (a) $\alpha = (0.6, 1.50, -0.5, -0.05)$ and (b) $\alpha = (0.4, 1.50, 0.2, 0.05)$. In both designs, the parameter values are chosen to ensure that the selection probability, $P(D=1),$ is approximately 60%. Let $h(X) := \alpha_0 + \alpha_1 X + \alpha_2 X^2 + \alpha_3 X^3$ be the selection index. The shape of $h$ is displayed in Figure (ref) and the performances of estimators are reported in Table (ref) and Figure (ref). Under design (a), $h(\cdot)$ is not monotone, so our two-step sieve estimator for $\beta_1$ is well centered around the true value. Its root-mean-squared error (RMSE) is close to that of the oracle estimator. Conversely, with design (b), $\beta_1$ remains unidentified so it suffers from a large RMSE as expected. In both designs, the OLS exhibits substantial misspecification bias. The Heckman’s MLE performs poorly in the non-monotone design due to misspecification but performs very well in the monotone design. This is because the selection index is close to linear in the effective support of $X$ and the error distribution is correctly specified in the monotone design.

We next generate Monte Carlo samples from the following DGP (referred to as DGP1 henceforth), where $X$ consists of two continuously distributed variables and the unobservables are joint normally distributed as (ref):

$X_1$ and $X_2$ are drawn from the standard normal distribution and independent of each other. The parameter values are set as: $$\alpha = (1.5, 0.5, -0.5, 0.2, 0.5, 1.0, -0.5), \quad \beta = (0.5, 0.5, 0.25).$$ The average selection probability across Monte Carlo samples is 52%.

In the latest design (DGP2 henceforth), we consider the scenario where $X$ consists of a continuously distributed variable, $X_1 \sim N(0,1)$, and a binary variable, $X_2 \sim Bernoulli(0.5)$. The remaining elements of DGP2 are otherwise identical to DGP1 except the selection process:

The parameter values are set as: $$\alpha = (0.2, -0.2, -0.5, 0.3, 0.1, 0.5, -0.3, 0.2), \quad \beta = (0.5, 0.5, 0.25).$$ The average selection probability is 66% under this DGP.

In both DGPs, we have at least one continuous covariate and the selection probability function $p_0(\cdot)$ exhibits sufficient nonlinearity, so our model point identifies $\beta_1$ and $\beta_2$. As showcased in Figures (ref)-(ref), the TPM and the Heckman selection model are misspecified and hence the OLS and MLE suffer from large bias for both DGPs. Table (ref) displays the RMSE and mean bias of each estimator. Heckman's MLE works particularly badly in DGP2. In contrast, our semiparametric estimator performs exceptionally well in DGP1 for both parameters as the oracle estimator outperforms our estimator by a very slight margin in terms of root-mean-squared errors (RMSE) and mean bias. In DGP2, it performs similarly to the oracle estimator for $\beta_1$, but shows a larger RMSE (0.130) than the oracle estimator (0.085) for $\beta_2$, possibly due to limited variations in $X_2$.

Finally, we evaluate the performance of lee2009bounds's and honore2020selection's bounds approaches using DGP2. We do not consider DGP1 because there is no binary treatment variable for which Lee’s bounds are applicable. We use a sample size of 100,000 instead of 5,000, which we use for point estimators. This is because the HH bounds are not reliably estimated with a moderate sample size, with which the bounds are often empty (in 93 iterations out of 1,000) or uninformative (including zero within the bounds in 615 iterations out of 1,000). With the 100,000 sample size, both bounds are reliably estimated. Figure (ref) displays the box plots of HH’s and Lee’s bounds. It is not surprising to observe that the Lee bounds consistently contain the true parameter value for $\beta_2$ because the bounds are very wide in this setup. The Lee bounds are never informative about the sign of the treatment effect as they include zero in every simulation under DGP2. In contrast, the HH bounds are significantly tighter than the Lee bounds. However, in most iterations, the HH bounds are not informative and never contain the true value because the model misspecifies the selection process.

These simulation exercises clearly demonstrate the practical usefulness of our semiparametric estimator. When at least one continuous covariate is present and the selection process exhibits nonlinearity, our estimator performs exceptionally well even with a modest sample size. The first-stage selection procedure is nonparametrically identified so assessing the nonlinearity in the selection equation is practically easy. In contrast to HH’s partially identifying linear selection model, our semiparametric model offers greater flexibility by not imposing linearity in the first stage while still point-identifying the parameters of interest. Consequently, our estimator can serve as a valuable alternative when the Lee bounds are excessively wide to provide meaningful insights. However, if there is no continuous variable or the selection process is genuinely linear, the HH bounds would be an excellent alternative to the Lee bounds.

Empirical application: gender and racial wage gaps in the US

We now demonstrate the empirical usefulness of our semiparametric model and its estimator using real-world data. We estimate the gender and racial wage disparities in the US. The reservation wage varies between different genders and ethnicities. Upon selection into employment, the distribution of unobserved factors can differ from that of the unemployed. Therefore, the effect of sample selection on observed wages should be taken into account to accurately calculate wage gaps. Following mora2008nonparametric and honore2020selection where they focus on racial wage gaps, we analyze Current Population Survey (CPS) data on wages from Arizona, California, New Mexico, and Texas. The data set covers the years 2003–2016 and includes 129,907 women. Among them, 26,698 are third-generation Mexican-Americans, while 103,209 are non-Hispanic whites. The remaining 118,418 men comprise 21,402 third-generation Mexican-Americans and 97,016 non-Hispanic whites. All individuals in the sample are aged between 25 and 62. In terms of employment, the percentage of women working is 64% for third-generation Mexican-Americans and 61% for non-Hispanic whites. The employment rates for men are 71% for Mexican-Americans and 67% for non-Hispanic whites, respectively.

The gender wage gap is estimated for Mexican-Americans and non-Hispanic whites separately to nonparametrically control for ethnicity. We use the log inflation-adjusted hourly wage as the outcome variable. In the latent outcome equation, we estimate the coefficient on the female dummy with age, age squared, experience, experience squared, education dummies (less than high school, some college, college, and advanced degree such as master's and doctorate, with high school as the base category), dummies for being a veteran and being married, state dummies (New Mexico as a base state), and year dummies as control variables. Age and experience serve as continuously distributed covariates in the selection equation for our semiparametric model. For the racial wage gap, we estimate the coefficient on the Mexican-American dummy with the same set of control variables separately for men and women.

We first estimate the Lee and HH bounds. As the Lee bounds are fully nonparametric, we compute the bounds conditional on the education level (high school and college) with no other covariates. For the HH bounds, we use the full set of control variables. We closely follow honore2020selection's implementation except that we employ probit regression in the first stage estimation of selection parameters in lieu of logit. The results are still very similar to the original results of HH with the logit first stage. Table (ref) presents the estimated bounds. For the racial wage disparities, the Lee bounds are not informative for college graduates, as they contain zero. For high school graduates, the bounds range from -25% to -7.5% for men and from -21% to -4.1% for women. In contrast, the HH bounds are highly informative and significantly narrower than the Lee bounds. The HH bounds range between -11.4% and -10.3% for men and between -8.9% and -6.6% for women. Regarding the gender wage gap, the Lee bounds suggest substantially lower wages for females, ceteris paribus. The bounds are wider for high school graduates (-32.5% to -10.5% for Mexican-Americans and -37.2% to -14.2% for whites). For college graduates, the bounds lie between -24.0% and -15.2% for Mexican-Americans, and between -28.4% and -11.9% for whites, indicating a potentially smaller gender wage gap among college graduates. The HH bounds for the gender wage gaps (-21.9% to -14.4% for Mexican-Americans and -22% to 17% for whites) are narrower than the Lee bounds but not as tight as for the racial gaps.

For point estimators, like in the simulation experiments, we consider the two-part model (“TPM”) using the OLS conditional on employment assuming random selection, the Heckman selection two-step estimator (“HS 2step”) and MLE (“HS MLE”), and our proposed semiparametric two -step estimator (“KL”). In the first stage estimation for the selection probability, we employ the sieve maximum likelihood estimator and predict $\hat{p}_0(\cdot),$ by including piecewise-polynomial (cubic b-spline) basis functions of age and experience with 5 knots, and their interactions with dummy variables. Most coefficients on sieve terms in the first stage estimation are highly significant across all the subsamples, indicating strong nonlinearity in the selection process. Given the prediction for the selection probability, $\hat{p}_i$, from the first stage, we estimate a partial linear model where the bias correction term $\lambda(\cdot)$ is approximated by cubic b-spline basis functions with 7 knots.

The estimation results are shown in Table (ref) for the racial wage gaps. It is surprising that the OLS assuming random selection and the Heckman selection approach using the MLE produce the same estimate of the racial wage gap for both men (-11.3%) and women (-7.8%). The Heckman two-step estimator, on the other hand, gives quite different results from the OLS and Heckit MLE with inflated standard errors. As it does not exploit the full information in the model, the two-step estimator tend to be less reliable. The OLS and Heckman MLE generally produce almost identical coefficient estimates for all covariates. In Heckman's approach, the null hypothesis of no correlation between the error terms cannot be rejected. Both OLS and Heckman MLE estimates are contained in the HH bounds, meaning that the linear selection models fail to capture any selection bias. As we can see in the first stage estimation, linearity of the selection process is strongly rejected, so the linear selection models are misspecified regardless of the assumption on the error terms. On the contrary, our semiparametric estimator shows a smaller magnitude of the racial wage disparity which is outside the HH bounds for both (-8.7%) men and women (-6.5%). Figure (ref)(a)-(b) compares our point estimates with the bounds estimates.

Our estimator also corrects selection bias in the coefficient estimates for the other covariates. The semiparametric estimator yields smaller wage premiums for higher education degrees (particularly for advanced degrees) for both men and women. Veteran status provides a higher wage premium for women than men, whose veteran premium is virtually negligible. Married men earn significantly higher wages than unmarried men, while married and unmarried women exhibit no significant difference in their hourly wage rates. The standard errors of our semiparametric estimates remain comparable to those obtained using Heckman MLE. These results effectively demonstrate the strong efficiency of our semiparametric estimator. There is a minimal difference in the estimated state fixed effects between the estimators. Both men and women are the highest-paid in California, followed by Arizona. The wage premium associated with residing in California and Arizona is higher for women compared to men by approximately 5% points relative to their counterparts in New Mexico.

The results on the gender wage gap also show interesting patterns as shown in Table (ref). For Mexican Americans, the OLS and Heckit MLE again produce the same estimates (around -19.5%). In contrast, our estimator indicates a smaller magnitude of the gender wage disparity (-18%). As the HH bounds in this case are quite wide, all the point estimates are contained in the bounds. For non-Hispanic whites, the patterns are quite the opposite. Heckit MLE seems to over-correct the selection bias, delivering a much smaller magnitude of the gender wage gap (-15.9%) than OLS (-20.9%). It also indicates much larger premiums on higher degrees (college and advanced degrees) compared to high school diploma than the OLS. These patterns are totally flipped in the semiparametric estimation. Our estimator produces a very similar estimate of the gender wage gap (-21.1%) to the OLS, while it produces lower wage premiums of higher degrees. Interestingly, the Heckman MLE estimate does not lie in the HH bounds, whereas our semiparametric estimate is still contained in the bounds as shown in Figure (ref)(c)-(d).

Finally, we incorporate heteroskedasticity of the error term in our semiparametric model and compute the heteroskedasticity-robust standard errors of the coefficients. Table (ref) presents the results. The robust standard errors are generally almost identical to the standard errors computed under the homoskedasticity assumption. The robust standard errors tend to be slightly larger than the non-robust errors, but occasionally slightly smaller.

This empirical application demonstrates that the widely used bounds approach proposed by lee2009bounds can yield uninformative bounds in analyzing crucial labor market outcomes, such as wages. The HH bounds offer a potential alternative, as they tend to provide tighter bounds. However, even the feasible non-sharp version of the HH bounds (as the sharp characterization relies on an uncountable infinity of moment inequalities) are computationally intensive. Moreover, the inference for these bounds hinges on resampling, which can be computationally demanding. In contrast, our semiparametric estimator is straightforward to implement in standard statistical packages like Stata and R, as it is a simple two-step plugin estimator. Any nonparametric estimator that satisfies the rate condition outlined in Section (ref) can be used for the first stage of estimation, which calculates the selection probability. The second stage can then be executed using standard partial linear regression. The asymptotically valid standard errors are computationally straightforward and incorporating heteroskedasticity is also very tractable. The estimator is efficient, as demonstrated in this application and simulations. Therefore, our estimator presents a valuable alternative that can be easily applied in cases where the bounds approaches fail to provide informative results, while it remains more robust than linear selection models. For researchers interested in correcting sample selection bias without resorting to unjustifiable parametric or distributional assumptions, we recommend reporting point estimates from our semiparametric selection model.

Concluding remarks

In this paper, we investigate point identification and efficient estimation of semiparametric selection models without imposing an exclusion restriction. We do not restrict the selection process to be linear, demonstrating identification of the model parameters when there is at least one continuous covariate and the linearity of the selection process is violated. The primary objective of our paper is to challenge the long-held belief that an exclusion restriction is necessary for semiparametric selection models. Bounds approaches for selection models are often motivated by this misconception. We present convenient and practical semiparametric estimators that accommodate non-monotone selection, heteroskedastic error, multiple control variables, and simple asymptotically valid inference. Our recommendation for applied researchers is to report point estimates using our semiparametric method when their preferred bounds are not sufficiently informative. The identifying conditions are readily verifiable in practice, as researchers simply need to ensure the presence of a continuous variable in the data and reject the linearity of the selection process.

In our simulations and empirical applications, we demonstrate that our method provides more robust estimates of parameters of interest compared to linear selection models such as the Heckman selection model and honore2020selection's model. While our semiparametric approach is not necessarily nested within Lee’s fully nonparametric model, it imposes more restrictive assumptions than Lee’s. Our model can permit parameter heterogeneity but it does not allow treatment effects to vary across different subpopulations. Extending our results to the case where the treated group and the untreated group have different treatment effects beyond the assumption made in honore2024sample would be an intriguing avenue for future research. Another promising research direction would be identification of semiparametric sample selection models with endogenous regressors.