Inference on High Dimensional Selective Labeling Models

Inference on High Dimensional Selective Labeling Models

\affil[1]{Department of Economics, Boston College} \affil[2]{Deptartment of Economics, Harvard University} \affil[3]{Department of Economics, Louisiana State University}

\vskip -.5in { \singlespace

The paper reconsiders the problem of inference on parameters in binary outcome models when these outcomes are subject to possibly endogenous censoring. Recently, these models have gained increasing interest in the computer science and machine learning literatures where the issue of endogenous sample selection is referred to as the {selective labels } problem. Such models are relevant in diverse empirical settings, including criminal justice, healthcare, and insurance. Notable recent studies in this area include lakkarujuetal2017, kleinbergetal, and \citet*{asheshsl}, which examine judicial bail decisions—where the outcome of whether a defendant fails to appear in court is observed only if the judge grants bail. Inference on such model parameters can be computationally challenging for two reasons. One is the nonconcavity of the bivariate likelihood function, and the other is the large number of covariates in each equation. Despite these hurdles, we propose a novel distribution free estimation procedure that is computationally friendly especially in the many covariates settings. The new method combines the semiparametric batched gradient descent algorithm introduced in \citet*{khanetal2021} with a novel sorting algorithm incorporated to control for selection bias. Asymptotic properties of the new procedure are established under increasing dimension conditions in both equations, and its finite sample properties are explored through a simulation study and an application using Stanford Open Policing Project data set pierson2020large. Extensions to models with endogenous treatment are also proposed.

{\bf Key Words:} Selective Label Models, Semiparametric Batched Gradient Descent, Selection Bias.

{\bf JEL Codes:} C14, C31, C35, C63.

\setcounter{page}{1}

{5pt} {5pt} {5pt} {5pt}

\setcounter{equation}{0}

Introduction

This paper addresses the challenge of inference on large-dimensional selective labeling models with binary outcomes. These models emerge in numerous domains where the observed binary outcomes result from the choices made by one of the agents within the system. Recently, they have garnered significant attention in the fields of computer science and machine learning, where this issue is known as the “selective labels problem”, or as an endogenous sample selection problem. Applications of these models span a wide range of areas, including criminal justice, healthcare, and insurance. For important recent work in this area, see for example lakkarujuetal2017, kleinbergetal and asheshsl. These authors focus on judicial bail decisions, where one observes the outcome of whether a defendant filed to return for their court appearance only if the judge in the case decides to release the defendant on bail. Letting $D_i$ denote the binary decision to grant bail, and $Y_i$ denote the binary outcome of the defendant returning for court appearance, they consider a model of the form

This process and the ensuing model can be best explained with the diagram below. The top node indicates the decision made by the agent which corresponds to a {\em yes} ($D_i=1$) or {\em no} ($D_i=0$) on individual $i$. The other observed dependent variable, corresponding to the two nodes beneath the top one, is denoted by $Y_i$, where $Y_i\in \{ 0,1,NA\}$ and denotes the resulting outcome (return to court in our example). The selective labels problem occurs because the observation of outcome $Y_i$ is constrained by the decision $D_i$ made by the agent:

Of course controlling for selection bias has a rich history in the econometrics literature, but usually for models where the outcome variable after selection is continuous. Seminal work include gronau, heckman1974 and in the semiparametric literature see, e.g. ahn1993 and neweyej\footnote{For work on a nonparametric sample selection model, see neweyetal}, who identify and estimate the model parameters as a semilinear or partially linear model for selected observations.\footnote{The partially linear model itself includes an expansive literature in the econometrics and statistics literature. Seminal work includes robinson1988, speckman1988. Recent developments for estimating this class of models includes, for large dimensional models, bellonietal2014, chernozhukhovetal2018, hsiaozhou2024, and for data combination, xaviermaurel2025.For Bayesian analysis of a wide class of simultaneous equation models, some of which can be nonlinear like in the Heckman model, see chib1, chib2. In these models unobserved components are assumed to be normally distributed, though finite mixtures of thereof can be considered. Nonparametric functions are considered there as is the case in robinson1988 but the dimension of the regressors both parametric and nonparmatric components of their semilinear model is restricted to be finite.}

More formally, the econometric model we consider:

The above system of equations is of a similar structure to that used in the classical sample selection model introduced in heckman1974. $I(\cdot)$ is the indicator function, $(z_{0,i}, \boldsymbol{Z}_i^{\mathrm{T}})^{\mathrm{T}}$ and $(x_{0,i}, \boldsymbol{X}_i^{\mathrm{T}})^{\mathrm{T}}$ denote vectors of observed regressors in selection and outcome equations, respectively. $D_i$ and $Y_i$ denote observed {\em binary} outcomes while $Y^*_i$ is only observed if $D_i=1.$ Also, $U_i$ and $V_i$ are unobserved random shocks while $\boldsymbol{\delta}_0$ and $\boldsymbol{\beta}_0$ are the parameters of interest. As we have mentioned above, this model differs from the classical selection model mainly in that we focus on the case where the outcome equation ((ref)) is also binary (and here we allow for large dimensional regressors).

The structure of the rest of the paper is organized as follows. In the next section we define our new (algorithmic) estimation procedures for the unknown regression coefficients in ((ref)) and ((ref)), which are designed to be computationally efficient, and hence, suitable to implement for models where the dimensions of $\boldsymbol{Z}_i$ and/or $\boldsymbol{X}_i$ are large. Section (ref) explores the asymptotic properties of the new methods, establishing their limiting distribution theories for models of increasing dimensions, which are gaining widespread and growing interest in the big data and machine learning literature, but have yet to be studied for this selective labeling model. Section (ref) explores the finite sample properties of our procedures by means of a simulation study and Section (ref) provides some empirical applications of our method. In Section (ref), we investigate extensions of the selective labeling models by allowing for endogenous treatment. Section (ref) concludes by summarizing our results and suggesting areas for future work. An Appendix collects tabular results from the simulation study and all the proofs of the main theorems.

Algorithmic Estimation Procedures

\setcounter{equation}{0}

This section introduces algorithmic estimation procedures for models ((ref)) and ((ref)), where $\boldsymbol{Z}_{e,i} \equiv (z_{0,i}, \boldsymbol{Z}_i^{\mathrm{T}})^{\mathrm{T}}\in \mathscr{Z}_e \subseteq R^{p_{Z}+1}$, $\boldsymbol{\delta}_0\in\mathscr{D}\subseteq R^{p_{Z}}$, $\boldsymbol{X}_{e,i} \equiv(x_{0,i}, \boldsymbol{X}_i^{\mathrm{T}})^{\mathrm{T}}\in \mathscr{X}_e \subseteq R^{p_{X}+1}$, and $\boldsymbol{\beta}_0\in\mathscr{B}\subseteq R^{p_{X} }$. $\boldsymbol{Z}_{e,i}$ and $\boldsymbol{X}_{e,i}$ are observed vectors of regressors in the selection and outcome equations, whose dimensions $p_Z$ and $p_X$ may increase with sample size $n$ but satisfy $\max\{p_Z, p_X\} \leq n$. $\boldsymbol{\delta}_0$ and $\boldsymbol{\beta}_0$ are unknown parameter vectors. $U_i, V_i$ are unobserved random shocks with joint distribution function $F(u,v)$, whose marginal distributions are given by $F_U(u)$ and $F_V(v)$.

We impose the following condition over the data set we observe.

The remainder of this section proposes two novel computationally efficient algorithms for estimating $\boldsymbol{\beta}_0$. The proposed methods are to first estimate the parameter vector $\boldsymbol{\delta}_0$ in the selection equation, and with that, use matching as in ahn1993 or series expansion as in neweyetal and neweyej to estimate the selection correction function, and finally estimate $\boldsymbol{\beta}_0$. We will not use rank estimation in either step because the dimension of $\boldsymbol{Z}_i$ and $\boldsymbol{X}_i$ potentially can be large and the resulting computational burden could be extremely heavy. Instead, we adopt iteration-based methods which feature simple implementation and fast computation speed.

The First-Step Estimator

We first introduce the algorithm for estimating $\boldsymbol{\delta}_0$. Define $\boldsymbol{\phi}_{\boldsymbol{\delta},q_{\boldsymbol{\delta}}}(\cdot) = \left(\phi_0(\cdot), \cdots, \phi_{q_{\boldsymbol{\delta}}}(\cdot)\right)^{\mathrm{T}}$, where $\phi_0(\cdot), \phi_1(\cdot), \cdots$ are a sequence of basis functions, and $q_{\boldsymbol{\delta}}$ is the order of sieve approximation used in the first-step estimation. For the choice of basis functions, see chen2007large and khanetal2021. The algorithm is described as follows.

\fbox{

}

Above algorithm is the sieve-based gradient descent estimator (SBGD) proposed by khanetal2021. Under some regularity conditions\footnote{The regularity conditions ensured point identification of $\boldsymbol{\delta}_0$. Recent work in khantamerwei2025 consider models where $\boldsymbol{\delta}_0$ is only partially identified and propose a two step procedure which converges to the identified set.}, khanetal2021 show that for $k$ sufficiently large, $\widehat{ \boldsymbol{\delta}}^k$ is consistent and asymptotically normally distributed under increasing dimensions.

The Second-Step Estimator

Denote the first-step estimator as $\widehat{\boldsymbol{\delta}}$. With $\widehat{\boldsymbol{\delta}}$ in hand, we now consider estimating $\boldsymbol{\beta}_0$. As in this SBGD, we will need to control for selection bias by explicitly estimating the selection correction function. To provide some intuition, suppose that we know the joint CDF of $U_i$ and $ V_i$, then the probability of $V_i<v$ conditioned on $U_i<u$ is given by \[ P(V_i<v|U_i<u) = \frac{F(u, v)}{F_U(u)} \equiv G(u, v). \] Define $Z_{0,i} = z_{0,i}+\boldsymbol{Z}_i^{\mathrm{T}}\boldsymbol{\delta}_0$ and $X_{0,i} = x_{0,i}+\boldsymbol{X}_i^{\mathrm{T}}\boldsymbol{\beta}_0$, we have that \[E(Y_i|\boldsymbol{Z}_{e,i}, \boldsymbol{X}_{e,i}, D_i = 1) = G(Z_{0,i},X_{0,i}).\] The key observation here is that $G(\cdot, \cdot)$ is increasing with respect to its second argument. Suppose further that we also know $\boldsymbol{\delta}_0$. Define $\widehat{X}_{i,k} = x_{0,i}+\boldsymbol{X}_i^{\mathrm{T}}\widehat {\boldsymbol{\beta}}^k$, then batch gradient descent algorithm based on the loss function in khanetal2021 immediately leads to the following iterative algorithm for estimating $\boldsymbol{\beta}_0$ using only observations with $D_i = 1$,

where $S_n=\sum_{i=1}^n D_i$ is the number of observations whose first-step outcome is 1. Obviously, update ((ref)) takes the index in the selection equation into consideration, so effectively controls for the selection bias.

However, since both $F(u,v)$ and $\boldsymbol{\delta}_0$ are unknown, the above algorithm is indeed infeasible. Note that the second issue can be easily resolved by plugging in our first-step estimator $\widehat{\boldsymbol{\delta}}$, while the first remains unsolved. An intuitive solution to such issue is to obtain an estimator for the conditional expectation $G(u, v)$ and then plug such estimator into update ((ref)). In this section we propose two methods to estimate such conditional expectation, one being local in nature and the other global.

The first local estimator uses matching to control for selection bias similar to ahn1993. To provide some intuition, suppose that the first-step estimator $\widehat {\boldsymbol{\delta}}$ is consistent and $\widehat{\boldsymbol{\beta}}^k$, the starting point in the $k$-th iteration, is close to $\boldsymbol{\beta}_0$. Define $\widehat{Z}_{i} = z_{0,i}+\boldsymbol{Z}_i^{\mathrm{T}}\widehat{\boldsymbol{\delta}}$, so long as $G$ is smooth enough, we have that \[ E\left(\left.Y_i\right|\boldsymbol{Z}_{e,i}, \boldsymbol{X}_{e,i}, D_i=1\right) = G(Z_{0,i}, X_{0,i}) \approx G(\widehat Z_{i}, \widehat X_{i,k}). \] The above result implies that for arbitrary $j\neq i$ such that $(\widehat Z_{j}, \widehat X_{j,k})$ is close enough to $(\widehat Z_{i}, \widehat X_{i,k})$, $Y_j$ can be used as a (noisy) replacement for $E\left(\left.Y_i\right|\boldsymbol{Z}_{e,i}, \boldsymbol{X}_{e,i}, D_i=1\right)$. This implies that $E\left(\left.Y_i\right|\boldsymbol{Z}_{e,i}, \boldsymbol{X}_{e,i}, D_i=1\right)$ can be estimated by a weighted combination of $Y_j$'s, where decreasing weight is assigned to each $Y_j$ as the distance between $(\widehat Z_{j}, \widehat X_{j,k})$ and $(\widehat Z_{i}, \widehat X_{i,k})$ increases. Such idea is similar to that in ahn1993.

To improve computational efficiency of the algorithm, in this paper we consider a nearest neighbor-type\footnote{Kernel-based weights are also easy to construct, which can be similarly done as in khanetal2021. However, constructing the weights involves $O(n^2)$ computational burdens in each round, which may cause heavy computation burdens, see Yao (2024).} weighting scheme. Define the Euclidean distance between $(\widehat Z_{j}, \widehat X_{j,k})$ and $(\widehat Z_{i}, \widehat X_{i,k})$ as \[d_{ij}^k = \left\Vert (\widehat Z_{j}, \widehat X_{j,k}) - (\widehat Z_{i}, \widehat X_{i,k})\right\Vert=\sqrt{\left(\widehat Z_{j} - \widehat Z_{i}\right)^2+ \left(\widehat{X}_{j,k} - \widehat{X}_{i,k}\right)^2}.\] For any $i$ with $D_i = 1$, rearrange the indices of $Y_j$ with $j\neq i$ and $D_j=1$ as $\varrho^k(i,1), \cdots, \varrho^k(i,S_n -1)$ such that $ d_{i,\varrho^k(i,1)}^k \leq \cdots \leq d_{i, \varrho^k(i,S_n -1)}^k$\footnote{If there is a tie then follow the original order of the indices.}. Then the weights based on $m$-nearest neighbor is given by

Given the above nearest-neighbor weighting scheme, the algorithm for estimating $\boldsymbol{\beta}_0$ based on the idea of matching is provided as follows.

\fbox{

}

We next propose an algorithmic approach which controls for selection bias by nonparametrically estimating the selection correction function globally using method of sieves. This was done for the standard selection model with linear outcome equation in \citet*{neweyetal} and neweyej. Let $\boldsymbol{\Phi}_q(\cdot, \cdot)$ be a $(q+1)^2$-dimensional vector of basis functions. Note that we can decompose $D_{i}Y_{i}$ as follows

where $\boldsymbol{\Pi}_q$ is the unknown pseudo true sieve parameter vector\footnote{Note that for any sequence of sieve functions $\{\Phi_{st}(u,v)\}_{s,t=0}^{\infty}$ that is complete in $C(R^2)$ space and any function $G(u,v)\in C(R^2)$, there exists a sequence of sieve coefficients $\{\pi_{st}\}_{s,t=0}^{\infty}$ such that $G(u,v)=\sum_{s,t=0}^{\infty}\pi_{s,t}\Phi_{st}(u,v)$. Then $\boldsymbol{\Pi}_q$ is the vector of first $(q+1)^2$ sieve coefficients. See chen2007large for more detailed discussion.}, and $\mathcal{E}_{i,q,k}$ is the error that can be decomposed as follows

When $G(u, v)$ is smooth enough, the first term on the right side of the above equation will be small as long as $q$ is large. Moreover, suppose again that the first-step estimator $\widehat {\boldsymbol{\delta}}$ is consistent and $\widehat{\boldsymbol{\beta}}^k$ is close to $\boldsymbol{\beta}_0$, then both second and third terms are small. Finally, the expectation of the last term conditioned on $\boldsymbol{Z}_{e,i}, \boldsymbol{X}_{e,i}$ and $D_i =1$ is zero. This naturally leads to an OLS-type estimator for $\boldsymbol{\Pi}_q$ given as follows

and the unknown conditional expectation function $G(u,v)$ can be estimated by $\widehat G^k(u,v)=\boldsymbol{\Phi}_q(u,v)^{\mathrm{T}}\widehat{\boldsymbol{\Pi}}_q^k$. Given the estimator of $G(u,v)$, we can plug it back to ((ref)), and conduct the update. The algorithm is formally detailed as follows.

\fbox{

}

\setcounter{equation}{0}

Statistical Properties

This section formally studies the statistical properties of the proposed iteration-based estimators. We start with introducing the conditions on the first-step estimator. In particular, we assume that $\widehat{\boldsymbol{\delta}}$ has the following asymptotic linear representation.

The next condition regulates the data generating process.

Statistical Properties of Matching-Based Estimator

This subsection studies the statistical properties of the matching-based BGD algorithm. To ease exposition, define $Z(\boldsymbol{\delta}) = z_0 + \boldsymbol{Z}^{\mathrm{T}}\boldsymbol{\delta}$, $Z_i(\boldsymbol{\delta}) = z_{0,i} + \boldsymbol{Z}_i^{\mathrm{T}}\boldsymbol{\delta}$, $X(\boldsymbol{\beta}) = x_0 + \boldsymbol{X}^{\mathrm{T}}\boldsymbol{\beta}$, and $X_i(\boldsymbol{\beta}) = x_{0,i} + \boldsymbol{X}_i^{\mathrm{T}}\boldsymbol{\beta}$. For any $\boldsymbol{\delta}\in\mathscr{D}$ and $\boldsymbol{\beta}\in\mathscr{B}$, define $\mathbb{Z}_{\boldsymbol{\delta}} = \{Z(\boldsymbol{\delta}): \boldsymbol{Z}_e\in\mathscr{Z}_e\}$ and $\mathbb{X}_{\boldsymbol{\beta}} = \{X(\boldsymbol{\beta}): \boldsymbol{X}_e\in\mathscr{X}_e\}$. We next introduce some additional technical conditions.

We finally introduce a condition that can be used as a sufficient condition for contraction mapping.

Given the above conditions, we now state our theorem regrading the convergence rate of the matching-based estimator, whose proof is provided in Appendix, Section (ref).

(ref) states that the matching-based estimator will be consistent when the number of nearest-neighbor points $m$ and the learning rate $\gamma_k$ are properly chosen. For the rate of convergence, the first term corresponds to the bias and the second term variance. It's intuitive that the bias increases with $m$, while the variance decreases with $m$. The optimal rate of $m$ is $\sqrt{\mathfrak{P}n\alpha_n}$, and the corresponding convergence rate of the matching-based estimator is $\sqrt{\log(n)\left(\frac{\mathfrak{P}^2}{n\alpha_n} + \frac{\mathfrak{P}\sqrt{\mathfrak{P}}}{\sqrt{n\alpha_n}}\right)}$. Note that due to the sizable bias, the resulting estimator is not guaranteed to be $1/\sqrt{n}$-consistent. This issue will be resolved for the sieve-based estimator. Nevertheless, the matching based estimator is easy to implement so can be used as a computationally efficient first-step estimator.

Statistical Properties of Sieve-Based Estimator

This section studies the statistical properties of the sieve-based BGD algorithm proposed in the previous section. We further introduce some technical conditions.

We finally introduce a condition that is similar to (ref).

Define \[\Xi_{1,n} = \sqrt{p_{X}}\mathscr{R}(q) + q^2C(0,q)^2\mathscr{R}_X(q)+\frac{\sqrt{p_{X}}q^4C(0,q)^3(p_{Z}C(1,q) + C(0,q)\sqrt{\log(n)p_{X}}}{\sqrt{n}}.\] Under the above conditions, we have the following result, whose proofs are provided in the Appendix, Section (ref).

(ref) provides asymptotic consistency of the sieve-based estimator for mildly increasing dimensionality as long as $\Xi_{1,n}\rightarrow 0$. Based on (ref), we can further establish the asymptotic linear representation for the sieve-based estimator. Define $\Xi_{2,n}$ as \[ \Xi_{2,n} = p_X\sqrt{p_{X}}q^4C(0,q)^{3}\Xi_{1,n}^{2}\left(q^{2}C(0,q)C(1,q)^2+C(2,q)\right). \] We have the following result.

(ref) establishes the asymptotic linear representation for our sieve-based estimator. Such a result is useful in that inference for the unknown parameter can be subsequently conducted. In particular, we have the following two corollaries.

The above corollary states that the linear combinations of arguments of $\boldsymbol{\widehat{\beta}}$ is $1/\sqrt{n}$-consistent and $\sqrt{n}\mathcal{W}^{\mathrm{T}}(\boldsymbol{\widehat{\beta}} - \boldsymbol{\beta}_0)$ is asymptotically normally distributed with asymptotic covariance matrix $\Sigma_{\mathcal{W},1} + \Sigma_{\mathcal{W},2}$. To estimate $\Sigma_{\mathcal{W},1} + \Sigma_{\mathcal{W},2}$, define $\widehat{X}_{i} = x_{0,i}+\boldsymbol{X}_i^{\mathrm{T}}\widehat{\boldsymbol{\beta}}$, \[\widehat {\boldsymbol{\pi}} = \left(\sum_{i=1}^n \boldsymbol{\phi}_{\delta,q_{\delta}}(\widehat{Z}_{i})\boldsymbol{\phi}_{\delta,q_{\delta}}(\widehat{Z}_{i})^{\mathrm{T}}\right)^{-1}\left(\sum_{i=1}^n \boldsymbol{\phi}_{\delta,q_{\delta}}(\widehat{Z}_{i}) D_i\right),\] \[\widehat{F}_{U}(u) = \boldsymbol{\phi}_{\delta,q_{\delta}}(u)^{\mathrm{T}}\widehat {\boldsymbol{\pi}}, \ \ \widehat{\nabla F}_{U}(u) = \nabla\boldsymbol{\phi}_{\delta,q_{\delta}}(u)^{\mathrm{T}}\widehat {\boldsymbol{\pi}},\] \[\widehat{\boldsymbol{\psi}}_{\delta}\left(\boldsymbol{Z}_e\right)=\boldsymbol{Z} -\frac{1}{n}\sum_{i=1}^{n}\boldsymbol{Z}_{i}\boldsymbol{\phi}_{\delta,q_{\delta}}\left(\widehat{Z}_{i} \right)^{\mathrm{T}}\left(\frac{1}{n}\sum_{i=1}^n\boldsymbol{\phi}_{\delta,q_{\delta}}\left(\widehat{Z}_{i} \right)\boldsymbol{\phi}_{\delta,q_{\delta}}\left(\widehat{Z}_{i} \right)^{\mathrm{T}}\right)^{-1}\boldsymbol{\phi}_{\delta,q_{\delta}}\left(z_0 + \boldsymbol{Z}^{\mathrm{T}}\widehat{\boldsymbol{\delta}}\right),\] \[\widehat{\boldsymbol{\Psi}}_{\delta}=\frac{1}{n}\sum_{i=1}^n\left(\widehat{\nabla F}_{U}\left(\widehat{Z}_{i}\right)\widehat{\boldsymbol{\psi}}_{\delta}\left(\boldsymbol{Z}_{e,i}\right)\boldsymbol{Z}_i^{\mathrm{T}}\right),\] \[\widehat{\mathcal{X}}_n\left(\nu_{Z},\nu_{X},\widehat{\boldsymbol{\beta}}\right)=\frac{1}{S_n}\sum_{i=1}^{n}D_{i}\boldsymbol{X}_{i}\boldsymbol{\Phi}_{q}\left(\widehat{Z}_{i},\widehat{X}_i\right)^{\mathrm{T}}\widehat{\boldsymbol{\Gamma}}_{n,q}^{-1}\left(\widehat{\boldsymbol{\delta}},\widehat{\boldsymbol{\beta}}\right)\boldsymbol{\Phi}_{q}\left(\nu_{Z},\nu_{X}\right),\] \[\widehat{\boldsymbol{\Pi}} _q = \left[\sum_{i=1}^n D_i\boldsymbol{\Phi}_q\left(\widehat{Z}_i, \widehat{X}_i\right)\boldsymbol{\Phi}_q\left(\widehat{Z}_i, \widehat{X}_i\right)^{\mathrm{T}}\right]^{-1}\times \left[\sum_{i=1}^n D_iY_i\boldsymbol{\Phi}_q\left(\widehat{Z}_i, \widehat{X}_i\right)\right],\] \[\widehat{G}(u, v) = \boldsymbol{\Phi}_q(u,v)^{\mathrm{T}}\widehat{\boldsymbol{\Pi}}_q, \widehat{\nabla_uG}(u, v) = \nabla_u\boldsymbol{\Phi}_q(u,v)^{\mathrm{T}}\widehat{\boldsymbol{\Pi}}_q, \widehat{\nabla_vG}(u, v) = \nabla_v\boldsymbol{\Phi}_q(u,v)^{\mathrm{T}}\widehat{\boldsymbol{\Pi}}_{q},\] \[\widehat{\boldsymbol{\Psi}}_S = \frac{1}{S_n}\sum_{i=1}^nD_i \widehat{\nabla_{v}G}\left(\widehat{Z}_{i},\widehat{X}_{i}\right)\left(\boldsymbol{X}_{i}-\widehat{\mathcal{X}}_n\left(\widehat{Z}_i,\widehat{X}_i,\widehat{\boldsymbol{\beta}}\right)\right)\boldsymbol{X}_{i}^{\mathrm{T}},\] \[\widehat{\Sigma}_{X,Z}=\frac{1}{S_n}\sum_{i=1}^nD_i\widehat{\nabla_{u}G}\left(\widehat{Z}_{i},\widehat{X}_{i}\right)\left(\boldsymbol{X}_{i}-\widehat{\mathcal{X}}_n\left(\widehat{Z}_i,\widehat{X}_i,\widehat{\boldsymbol{\beta}}\right)\right)\boldsymbol{Z}_{i}^{\mathrm{T}}.\] For any fixed $p_X\times 1$ vector $\mathcal{W}$, further define

Then $\Sigma_{\mathcal{W},1} + \Sigma_{\mathcal{W},2}$ is estimated by $\widehat{\Sigma}_{\mathcal{W},1}+\widehat{\Sigma}_{\mathcal{W},2}$.

Monte Carlo Simulations

\setcounter{equation}{0} This section conducts some simulation experiments to evaluate the performance of the proposed estimators. The data generating process we consider is ((ref)) and ((ref)). We generate $\widetilde z_{0,i},x_{0,i}\sim_{\text{i.i.d.}}N(0,1), \widetilde z_{1,i},x_{1,i}\sim_{\text{i.i.d.}}\text{Bernoulli}(0.5), \widetilde z_{2,i},x_{2,i}\sim_{\text{i.i.d.}}\text{Poisson}(2), z_{j,i}, x_{j,i}\sim_{\text{i.i.d.}} \frac{1}{\sqrt{2}}(\chi^2(1) - 1)$ for $j\geq 3$, and $\widetilde z_{j,i}$ is independent of $x_{l,i}$ for any $i$, $0\leq j\leq p_Z$ and $0\leq l\leq p_X$. Then, we set $z_{j,i} = \widetilde z_{j,i}$ for $0\leq j\leq 2$, and $z_{j, i} = \widetilde{z}_{j,i} + \sum_{l=3}^{p_X}\varsigma_{l}x_{l,i}$, where $\varsigma_l$ is randomly drawn from $0.4\cdot (U(0,1) - 0.5)$ in each round of . $(U_i, V_i)$ is i.i.d. across $i$ and is independent of $\boldsymbol{Z}_{e,i}$ and $\boldsymbol{X}_{e,i}$. We set \[\boldsymbol{\delta}_0 = \left(-1, 1, 0.5, -0.5, 1.2, -1.4, -2.2, 1.8, 0.05\cdot \boldsymbol{1}_{10}, -0.05\cdot \boldsymbol{1}_{10}\right)^{\mathrm{T}}\] and \[\boldsymbol{\beta}_0 = \left(1, -1, 0.5, -0.5, 0.8, -1.0, -1.8, 2.1, 0.05\cdot\boldsymbol{1}_{10}, -0.05\cdot\boldsymbol{1}_{10}\right)^{\mathrm{T}},\] where $\boldsymbol{1}_{10}$ is a $10\times 1$ row vector whose arguments are all 1. In the following we report simulation results under different sample sizes $n$ and setups of joint distributions of $U_i$ and $V_i$.

We consider four competing methods, two being parametric and two semiparametric. The first method is parametric estimation using nonlinear least squares (two-step NLS), which is similar to Heckman's two-step estimation but accounts for the binary response in the second stage. To apply such method, we assume that $U_i$ and $V_i$ in ((ref)) and ((ref)) have zero mean and unit variance, and are jointly normally distributed with covariance $\rho$. Then we can jointly estimate $\boldsymbol{\delta}_0$, $\boldsymbol{\beta}_0$ together with $\rho$. In particular, let $F_1(\cdot)$ and $F_2(\cdot, \cdot, \rho)$ denote the CDF's of univariate standard normal distribution and bivariate normal distribution with covariance matrix $[\sigma_{ij}]_{2\times2}$, where $\sigma_{11} = \sigma_{22} = 1$, and $\sigma_{12} = \sigma_{21} = \rho$. Also define $\overline{\boldsymbol{Z}}_{e,i} = \left(1, z_{0,i}, \boldsymbol{Z}_i^{\mathrm{T}}\right)^{\mathrm{T}}$, $\overline{\boldsymbol{X}}_{e,i} = \left(1, x_{0,i}, \boldsymbol{X}_i^{\mathrm{T}}\right)^{\mathrm{T}}$, $\overline{\boldsymbol{\delta}} = \left(c_{\delta, 0}, c_{\delta, 1}, \boldsymbol{\delta}^{\mathrm{T}}\right)^{\mathrm{T}}$, $\overline{\boldsymbol{\beta}} = \left(c_{\beta, 0}, c_{\beta, 1}, \boldsymbol{\beta}^{\mathrm{T}}\right)^{\mathrm{T}}$. In the first step, we minimize the following loss function \[ L_{1,n}(\overline{\boldsymbol{ \delta}}) = \frac{1}{n}\sum_{i=1}^n (D_i - F_1(\overline{\boldsymbol{ Z}}_{e,i}^{\mathrm{T}}\overline{\boldsymbol{ \delta}}))^2 \] and obtain the minimizer $\widehat{\overline{\boldsymbol{\delta}}}$. Then in the second step, we minimize the following loss function \[ L_{2,n}(\overline{\boldsymbol{\beta}}, \rho) = \frac{1}{S_n}\sum_{i=1}^nD_i\left(Y_i - \frac{F_2(\overline{\boldsymbol{Z}}_{e,i}^{\mathrm{T}}\widehat{\overline{\boldsymbol{\delta}}}, \overline{\boldsymbol{X}}_{e,i}^{\mathrm{T}}\overline{\boldsymbol{\beta}}, \rho)}{F_1(\overline{\boldsymbol{Z}}_{e,i}^{\mathrm{T}}\widehat{\overline{\boldsymbol{\delta}}})}\right)^2 \] and obtain the minimizer $\widehat{\overline{\boldsymbol{\beta}}}$. Then the two-step NLS estimators for $\boldsymbol{\delta}_0$ and $\boldsymbol{\beta}_0$ are given by $\widehat c_{\delta, 1}^{-1}\widehat{\boldsymbol{\delta}}$ and $\widehat c_{\beta, 1}^{-1}\widehat {\boldsymbol{\beta}}$.

The second method is parametric maximum likelihood estimation (MLE). Like in the first method, we also assume that $U_i$ and $V_i$ in ((ref)) and ((ref)) are jointly normally distributed, then the log-likelihood function is then given by

Suppose the MLE estimators are given by $\widehat{\overline{\boldsymbol{\delta}} }$ and $\widehat{\overline{\boldsymbol{\beta}}}$, then the MLE estimators for $\boldsymbol{\delta}_0$ and $\boldsymbol{\beta}_0$ are given by $\widehat c_{\delta, 1}^{-1}\widehat{\boldsymbol{\delta}}$ and $\widehat c_{\beta, 1}^{-1}\widehat{\boldsymbol{\beta}}$.

The third method is semiparametric estimation based on matching. In particular, we first obtain the estimator of $\delta_0$ in the first step, then we conduct Algorithm 1. To improve the computational efficiency, in the first and second step estimation, we use nearest neighbor matching with $m=[(\log(n))^{1.1}]$ and $m=[(\log(S_n))^{1.1}]$, where recall that $S_n = \sum_{i=1}^nD_i$. We update 2000 times for both first- and second-step estimation.

The fourth method is semiparametric estimation based on series approximation. For the sieve functions, we consider the Legendre polynomials used in khanetal2021, and use tensor products of one-variate sieve functions as the sieve functions for bivariate functions. The order of sieves is chosen to be 31 for the first-step estimation and 15 for the second-step estimation. Finally, the stopping rule is the same as in that of khanetal2021 with the tolerance being $10^{-6}$.

We report the bias and the root mean squared error (RMSE) of the second-step estimator for all methods. Let the simulation be repeated for $R$ times, and in the $r$-th round of repetition the estimator of $\beta_{j, 0}$ be $\widehat \beta_j^r$. Then the bias of the estimator is given by $B_j = |\frac{1}{R}\sum_{r=1}^R\widehat \beta_{j}^r - \beta_{j,0}|$, and the RMSE is given by $RMSE_j = \sqrt{\frac{1}{R}\sum_{r=1}^R(\widehat \beta_{j}^r - \beta_{j,0})^2}$. We report the bias and RMSE for $\beta_{1,0}$ through $\beta_{8,0}$. We also report the total bias and RMSE, which are defined as $B_{\text{all}} = \sum_{j=1}^{p_X}B_j$ and $RMSE_{\text{all}} = \sqrt{\sum_{j=1}^{p_X} RMSE_j^2}$. We choose $R=100$. Results are reported in (ref) to (ref).

(ref) through (ref) correspond to fours setups of the error terms. The first three setups feature model misspecification, while the last one features correct specification. Moreover, in the first setup both $U_i$ and $V_i$ are symmetric and are positively correlated. For the second step, $U_i$ and $V_i$ are also positive correlated but are both heavily skewed. The third setup corresponds to the case in which $U_i$ and $V_i$ are both bounded and are strongly (negatively) correlated with each other. Several insights can be drawn from the simulation results. First of all, when the model is misspecified, that is, the error terms $U_i,V_i$ are not jointly normally distributed, parametric estimators always have nonvanishing sizable bias, while the semiparametric sieve-based estimator constantly has minimal bias which is close to zero. This highlights the fact that the key advantage of semiparametric estimation lies in the robustness to model misspecifications. Further comparisons of the RMSE between different methods reveal that even though parametric estimators have $O(1/\sqrt{n})$ standard deviation, the sizable bias significantly contaminates the RMSE so that the sieve-based estimator always outperforms two parametric estimators in terms of RMSE. Notably, when the model is correctly specified ((ref)), both two-step NLS and joint MLE estimators have trivial bias, and joint MLE estimators have smallest RMSE among all methods. However, even in this case, the sieve-based semiparametric estimator is competitive compared with two-step NLS in terms of RMSE. We finally point out that the performance of matching-based semiparametric estimator is less promising compared with the sieve-based on in terms of both bias and RMSE. This may be explained by limited number of iterations (2000 iterations for both steps) or small neighborhood size $k$ ($k= \log(n)^{1.1}$ for first step and $k = \log(S)^{1.1}$ for second step). Nevertheless, as we mentioned above, the matching-based estimator can be used as a decent initial point for sieve-based estimator.

Empirical Application: Stanford Policing Data

In this section we apply our new algorithm to analyze the Stanford Open Policing Project data set pierson2020large. This large-scale data set records traffic stops made by police officers and various features related to the stops so can be used to provide insights into the potential racial disparities during police stops, which has been widely studied in the existing literature goel2016precinct.

In this paper we use the data set from Nashville, which contains a total of 3092351 raw observations and 2608109 observations after data clearing. To analyze this data set, we define the first-stage binary outcome as, after a stop has been made, whether the police officer further decides to conduct a search for the person or vehicle. Conditioned on that a search has been made, the second stage outcome is whether the police officer found any contraband items such as drugs or weapons. For the first stage equation, we consider a series of regressors including police precinct, the reason for the stop, as well as the gender, race, and age of the stopped subject. For the second stage equation, we consider all the regressors in the first stage but the precinct so that the precinct is used as the exclusion restriction. We provide detailed description of the variables in (ref).

When conducting SBGD estimation, we set the coefficient of Reason1 to be 1 for both selection and outcome equations, whose reason is detailed as follows. When police officer stops individuals to conduct investigation, it's highly likely that the police officer notice something unusual and hence a search is more likely to follow and contraband items are more likely to be found. This implies that the coefficients of Reason1 in both selection and outcome equations should be positive. We also note that the ratio of maximum and minimum eigenvalues of the covariate matrices for both selection and outcome equation is extremely large, which may lead to poor numerical performance during iterations. For example, in this case small learning rate $\gamma_k$ has to be chosen during the estimation of both stages to guarantee convergence of the algorithms, which leads to extensive rounds of iterations. To solve this issue, we use Gram-Schimidt method to orthogonalize the covariate matrix so that the resulting covariate matrix has unit column variance and zero column correlation, and the first column of the orthogonalized covariate matrix is equal to Reason1. See (ref) for more details of the empirical setups. The estimation results are reported in (ref).

Several insights can be drawn from the results. First of all, age is negatively related to the conditional probability of being searched after a stop is made; it is also negatively associated with the probability that illegal items are found. Second, compared with male individuals, female individuals are less likely to be searched after the stop, and conditioned on that the stopped individual is searched, female individuals are less likely to be found carrying illegal items. Finally, for racial disparities during stops, we find that the coefficient of Black in the selection equation is significantly positive, while that of white is negative. This implies that when other conditions are held constant, individuals are more likely to be searched (compared with Hispanic individuals) if they are black, while white people are less likely to be searched. We further notice that for the outcome equation, the coefficients of Black and White almost coincide with each other, implying that conditioned on being searched, black people are almost as likely as the white people to be found carrying contraband items. Such a result implies that even though black community may confront certain disparities in terms of a search decision, the hit rate (of finding illegal items) does not differ across races.

Selective Labeling with Endogenous Treatment

Multinomial Choice Models

We illustrate here how our proposed method can be used to estimate the standard multinomial response model where the dependent variable takes one of $J+1$ mutually exclusive and exhaustive alternatives numbered from $0$ to $J$. Following the notation similar to that used in KOT2021, for individual $i$, alternative $j$ is assumed to have an unobservable indirect utility $Y_{ij}^*$. The alternative with the highest indirect utility is assumed to be chosen. Thus the observed choice $Y_{ij}$ can be defined as \[ Y_{ij} = I\left(Y_{ij}^* > Y_{ij}^*, \forall k \neq j\right) \] with the convention that $Y_{ij} = 0$ indicates that the choice of alternative $j$ is not made by individual $i$. As is standard in the literature, an assumption of joint continuity of the indirect utilities rules out ties with probability one. In addition, we maintain the familiar linear form for indirect utilities\footnote{Our method can be applied to more general models with indirect utilities $y_{i\jmath}^{*}=u_{\jmath}(x_{i\jmath}'\beta_{0},-\epsilon_{i\jmath})$, $\jmath=1,2$, where $u_{\jmath}(\cdot,\cdot)$'s are unknown (to econometrician) $\mathbb{R}^{2}\mapsto\mathbb{R}$ functions strictly increasing in each of their arguments. It will be clear that our rank procedure does not rely on the additive separability of the regressors and error terms.}

where $\boldsymbol{\beta}_0$ is a $p$-dimensional vector of unknown preference parameters of interest whose first component is normalized to have absolute value 1 (scale normalization). Note that for alternative $j= 0$, the standard (location) normalization $Y_{i0}^* = 0$ is imposed. The vector $\boldsymbol{U}_i\equiv(U_{i1},...,U_{iJ})^{\mathrm{T}}$ of unobserved error terms, attained by stacking all the scalar idiosyncratic errors $U_{ij}$, is assumed to be jointly continuously distributed and independent of the $p\times J$-dimensional vector of regressors $\boldsymbol{X}_i\equiv(\boldsymbol{X}_{i1}^{\mathrm{T}},...,\boldsymbol{X}_{iJ}^{\mathrm{T}})^{\mathrm{T}}$\footnote{We impose the independence restriction here to simplify exposition. As explained in KOT2021, this matching-based approach allows $\boldsymbol{U}_i$ to be correlated with individual-specific regressors.}. We stress that model ((ref)) is rather general. By properly re-organizing $\boldsymbol{X}_{ij}$'s and $\boldsymbol{\beta}_0$, ((ref)) can accommodate both alternative-specific and individual-specific covariates\footnote{Note the identification of models with both alternative-specific and individual-specific regressors will need to take two steps, of which the first step only identifies the coefficients on alternative-specific regressors.}

Previous semiparametric contributions to estimating this model include lee95, who proposes a profile likelihood approach, extending the results in kleinspady for the binary response model. ahnetal2018 propose a two-step estimator that requires nonparametric methods but show the second step is of closed-form. shietal2018 also propose a two-step estimator in panel setups exploiting a cyclic monotonicity condition, which also requires a high dimensional nonparametric first stage, but whose second stage is not closed-form as ahnetal2018 is.

KOT2021 proposed a local rank procedure. We define it here for the special case where $J=2$, in which case the model is

One way to estimate $\boldsymbol{\beta}_0$ for this model proposed in KOT2021 was a weighted rank type estimator:\footnote{Here the weights correspond to binary, “exact" matches of each component of the vector $x_2$. For continuously distributed regressors they were replaced with kernel weights in KOT2021.} rank correlation estimator, analogous to the maximum rank correlation (MRC) estimator proposed in han1987non, defined as the maximizer, over the parameter space ${\cal B}$, of the following objective function

where above we denote pairs of individuals by $i,\ell$ and recall the second subscript denotes the choice. The motivations for the above estimator were robustness properties, notably when many of the regressors were discrete. Note that the matching of regressors is analogous to how we controlled for selection bias in the earlier part of the paper for the selective labeling model. We also note that the nonsmoothness and nonconvexity of the objective function make the implementation difficult, especially when there are a moderately large number of regressors.

For the model at hand we propose the following algorithm, which is analogous to the matching algorithm earlier in the paper (Algorithm 1) and keep notation as close as possible to that used there\footnote{ A sieve based algorithm could also be considered but we omit that here.}.

This section extends the previous estimation method to selective labeling models with endogenous treatments. In his seminal work, leebounds considers partial identification of the treatment parameters in treatment effect models with attrition but does not allow for explanatory variables. SEMENOVA2025106055 extends leebounds's bound and studies the asymptotic properties of the proposed bounds under both fixed and high dimensionality. For more discussion on leebounds's model and its recent development, see SEMENOVA2025106055 and references therein.

Extending our selective labeling model to allow for endogenous treatment status which is denoted by $T_i$ can be expressed as

where $\boldsymbol{R}_{e,i} = (r_{0,i}, \boldsymbol{R}_{i}^{\mathrm{T}})^{\mathrm{T}}\in R^{p_R+1}$ and $\boldsymbol{\varphi}_{0}\in R^{p_R}$. In the above system of equations, the observed binary variable $D_i$ indicates whether or not the $i$-th agent outcome variable is observed in the sample, $Y_i$ denotes the observed outcome variable for the selected sample, with selection governed by $D_i$, and $T_i$ denotes an observed binary variable indicating treatment status, whose coefficient in the outcome equation, $\tau_{2,0}$, is the main parameter of interest. Such system of equations generalizes leebounds's original model by explicitly modeling the determination of the treatment status. Moreover, since the unobserved random errors $(W_i, U_i,V_i)$ are potentially mutually correlated, the treatment status is allowed to be endogenous and correlated with both the selection status $D_i$ and the outcome $Y_i$. We also note that $\boldsymbol{R}_{e,i}, \boldsymbol{Z}_{e,i},\boldsymbol{X}_{e,i}$ can be large dimensional, which, similar to the selective labeling models, makes the estimation of the above system computationally intensive.

The above system of selective labeling model with endogenous treatment nests many models in important work in the literature. In a model where for $D_i=1$, $Y_i$ was linear and there were no regressors besides $T_i$, leebounds considered {\em partial} identification of $\tau_{2,0}$. In a model where there was no selection/attrition issue so $D_i$ was identical to 1, identification and estimation were considered in vytlacilyildiz, abhausmankhan, shaikhvytlacil. vytlacilyildiz attained point identification under a monotonicity condition as well as support conditions on exogenous covariates effecting $Y_i$. See also khanmaurelzhang and chenkhantang for point identification results in similar models. However, none of the methods proposed in the above papers are applicable to the model above where there are endogenous treatment, attrition, and selective labeling on the same time, even for low dimensional models. In contrast, the methods we introduced in Section (ref) can be extended to estimate $\boldsymbol{\varphi}_0 , \boldsymbol{\delta}_0, \boldsymbol{\beta}_0$, and $\tau_{1,0},\tau_{2,0}$ even when the covariates in the system are all large dimensional.

Before we demonstrate our main algorithm, we first show point identification for the parameters in the above system under a set of conditions. Denote the marginal CDF of $W_i$, the joint CDF of $W_i$ and $U_i$, and the joint CDF of $W_i$, $U_i$, and $V_i$ as $F_W(w)$, $F_{W,U}(w,u)$, and $F_{W,U,V}(w,u,v)$, respectively. Moreover, rearrange the order of the regressors such that $\boldsymbol{R} = (\boldsymbol{R}_c^{\mathrm{T}}, \boldsymbol{R}_d^{\mathrm{T}})^{\mathrm{T}}$, where $\boldsymbol{R}_c$ is continuous and $\boldsymbol{R}_d$ is discrete. Similarly, we write $\boldsymbol{Z} = (\boldsymbol{Z}_c^{\mathrm{T}}, \boldsymbol{Z}_d^{\mathrm{T}})^{\mathrm{T}}$ and $\boldsymbol{X} = (\boldsymbol{X}_c^{\mathrm{T}}, \boldsymbol{X}_d^{\mathrm{T}})^{\mathrm{T}}$. We impose the following conditions.

Based on the above conditions, we have the following identification results.

Given the identification results, now we can illustrate our estimation method in detail, which involves four steps. In particular, we first sequentially estimate $\boldsymbol{\varphi}_0$, $\boldsymbol{\delta}_0$, and $\boldsymbol{\beta}_0$ in the first three steps. Then with the estimators in hand, we finally estimate $\tau_{1,0}$ and $\tau_{2,0}$. First of all, note that the parameter $\boldsymbol{\varphi}_0$ in equation ((ref)) indicating the treatment status can be readily estimated using the SBGD algorithm proposed in khanetal2021 because it is a binary choice process. Denote the first-step estimator as $\widehat{\boldsymbol{\varphi}}$.

In the following, we use $\boldsymbol{\Phi}_q$ to denote generic vectors of sieves functions, whose length depends on $q$ and arguments may have different dimensions depending on the specific functions we would like to approximate. Correspondingly, we use $\boldsymbol{\Pi}_q$ to denote the pseudo true sieve parameters. Now we proceed to the second step. Note that

where $G_1(\cdot, \cdot)$ and $G_2(\cdot, \cdot)$ are both monotonically increasing in their second argument. Following the intuition of Algorithm 2 in Section (ref), we have the following estimation procedure for $\boldsymbol{\delta}_0$.

\fbox{

}

In the third step, we estimate $\boldsymbol{\beta}_0$ based on the first-step estimator $\widehat{\boldsymbol{\varphi}}$ and second-step estimator $\widehat{\boldsymbol{\delta}}$. The estimation procedure is also motivated by the following observation

where, similar to the previous step, $G_3(\cdot, \cdot, \cdot)$ and $G_4(\cdot, \cdot, \cdot)$ are both monotonically increasing in the third argument. Then we get the following estimation procedure for $\boldsymbol{\beta}_0$.

\fbox{

}

Denote the third-step estimator as $\widehat{\boldsymbol{\beta}}$, now we proceed to the final step for the estimation of the treatment effects $\tau_{1,0}$ and $\tau_{2,0}$. We first describe the estimator for $\tau_{1,0}$, whose idea can be similarly extended to the estimation of $\tau_{2,0}$. Note that under (ref), we can identify the following two functions

when $ (R_{0,i}, \boldsymbol{Z}_{e,i}^{\mathrm{T}})^{\mathrm{T}}$ is interior point of $\overline{\mathcal{Z}}_e$, then $\nabla_w F_{W,U}(w, u+\tau_{1,0})$ and $\nabla_w F_{W,U}(w, u) = -\nabla_w\left(F_U(u) - F_{W,U}(w, u)\right)$ can also be identified. So we can match the values of the derivatives of the CDF functions and their arguments to estimate $\tau_{1,0}$. Typically, under (ref), $\tau_{1,0}$ uniquely minimizes the following loss function

where $\omega_1(w,u)\geq 0$ is a well-chosen weight function. The above observation motivates the following algorithm for estimating $\tau_{1,0}$.

\fbox{

}

We finally describe the algorithm for estimating $\tau_{2,0}$. Under (ref), the following two probabilities can be identified

when $(R_{0,i}, Z_{0,i}, \boldsymbol{X}_i^{\mathrm{T}})^{\mathrm{T}}$ is interior point of $\overline{\mathcal{X}}_e$. So the cross partial derivative $\nabla_{wu}F_{W,U,V}(w, u+\tau_{1,0}, v+\tau_{2,0})$ and $\nabla_{wu}F_{W,U,V}(w, u, v)$ can be identified. Then since $\nabla_{wu}F_{W,U,V}(w, u, v)$ is non-decreasing with respect to $v$, under (ref), $\tau_{2,0}$ uniquely minimizes the following loss function \[ \int \left(\nabla_{wu}F_{W,U,V}(w, u+\tau_{1,0}, v+\tau_{2,0}) - \nabla_{wu}F_{W,U,V}(w, u+\tau_{1,0}, v+\tau_{2})\right)^2\omega_2(w,u,v)dwdudv, \] where $\omega_2(w,u,v)$ is some well-chosen weight function. This immediately leads to the following algorithm for estimating $\tau_{2,0}$.

\fbox{

}

\ \ \

Conclusions

This paper considers estimation and inference for large dimensional semiparametric selective labeling models. Statistically these models have a similar structure to sample selection models with binary, as opposed to linear outcome equations in the second stage. It is this binary/binary structure which makes computation of the model particularly difficult when compared to the standard selection model, especially for large dimensional (i.e many regressor) models.

To address this problem we propose novel algorithmic procedures which are computationally fast, and derive their asymptotic properties even for the case where the dimension increases with the sample size. We demonstrate the finite sample properties of our proposed procedures by a simulation study.

Our work here motivates areas for future research. For example to further ease implementation, a bivariate penalization scheme would be useful for model selection in this settings, and its asymptotic validity would need to be proven. Furthermore, the usefulness of our methods in other empirical settings in economics, biostatistics and medicine would be worthy of exploration.