The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
113,584 characters
Treatment Evaluation at the Intensive and Extensive Margins
\begin{titlepage}
\title{Treatment Evaluation at the Intensive and Extensive Margins \thanks{ \scriptsize We would like to thank Isaiah Andrews, Andres Aradillas-Lopez, Kevin Chen, Juan Carlos Escanciano, Patrik Guggenberger, Mark Henry, Keisuke Hirano, Edward Kennedy, Giovanni Mellace, Peter Schochet, Vira Semenova, Neil Shephard, Mikkel S\o lvsten, Elie Tamer, Davide Viviano, Chris Walker, as well as the participants of the Aarhus, Harvard, and PennState Econometrics Seminars for their comments and suggestions that helped to improve the work. All remaining errors are ours.}}
\author{Phillip Heiler\thanks{{\scriptsize Aarhus University. Department of Economics and Business Economics, TrygFonden's Centre for Child Research, Fuglesangs All\'e 4, 8210 Aarhus V, Denmark, email: [email removed]. Part of this research was conducted during main affiliation with Harvard University, Department of Economics. 02138 Cambridge MA, United States.}} \quad
Asbj\o rn Kaufmann\thanks{{\scriptsize Aarhus University. Department of Economics and Business Economics, Fuglesangs All\'e 4, 8210 Aarhus V, Denmark.}} \quad
Bezirgen Veliyev\thanks{{\scriptsize Aarhus University. Department of Economics and Business Economics, Fuglesangs All\'e 4, 8210 Aarhus V, Denmark.}}\\
}
\date{}
\maketitle
\thispagestyle{empty}
\begin{abstract} \singlespacing \small
\vspace{-4em}
This paper provides a solution to the evaluation of treatment effects in selective samples when neither instruments nor parametric assumptions are available. We provide sharp bounds for average treatment effects under a conditional monotonicity assumption for all principal strata, i.e.~units characterizing the complete intensive and extensive margins. Most importantly, we allow for a large share of units whose selection is indifferent to treatment, e.g.~due to non-compliance.
The existence of such a population is crucially tied to the regularity of sharp population bounds and thus conventional asymptotic inference for methods such as Lee bounds can be misleading. It can be solved using smoothed outer identification regions for inference. We provide semiparametrically efficient debiased machine learning estimators for both regular and smooth bounds that can accommodate high-dimensional covariates and flexible functional forms.
Our study of active labor market policy reveals the empirical prevalence of the aforementioned indifference population and supports results from previous impact analysis under much weaker assumptions.
\end{abstract}
\vspace{1em}
\noindent \textbf{Keywords:} Double/debiased machine learning; Lee bounds; Partial identification; Principal strata; Sample selection \\
\textbf{JEL classification:} C13, C14, C21
\end{titlepage}
\setcounter{page}{1}
\newpage
\section{Introduction}
This paper deals with the evaluation of causal effects of a binary treatment when the outcome is only selectively observed and no instruments are available. Such missing outcome data or sample selection is ubiquitous and a threat to internal validity \citep{heckman1974shadow,heckman1979sample}. Typical examples include outcome data such as wages that are only observed if units are employed \citep{lee2009training}, missing survey responses \citep{bernhardt2024howdoes}, or attrition \citep{zhang2003estimation}. In the context of impact evaluation, a treatment can often change the sample selection status for some units. For example, if the treatment is an active labor market policy such as job training, it likely affects the probability of employment for some participants, e.g.~due to accumulation of human capital (positive) or opportunity costs to job search (negative). For some units, however, the policy might be irrelevant with regards to employment or sample selection, but still relevant in terms of their earnings or other outcomes. Hence, there are potentially heterogeneous units at the intensive and extensive margins. These units can be classified within principal strata defined by their potential selection status.
We consider set identification, estimation, and inference for all principal strata, encompassing the complete intensive and extensive margins, under a weak conditional monotonicity assumption. Weak monotonicity allows for the simultaneous presence of always-takers, compliers, defiers, and never-takers.\footnote{In contrast to the instrumental variables literature, principal strata here refer to potential selection statuses when treatment is exogenously set to zero and one respectively. See Section \ref{sec_model}.}
The main challenge is the presence of units whose selection probabilities are unaffected by the treatment. The existence of such units is most apparent when the treatment is subject to additional non-compliance after assignment and intent-to-treat effects are analyzed. In this case, units who have information or preferences that lead them to not participate after treatment assignment will likely also not change their selection behavior, e.g.~selection into employment, as a consequence. Further, there are policies that do not boost or hinder employability for some units, e.g.~due to a lack of signal for their potential employers in the short run, but increase human capital boosting productivity and earnings. Similarly, in field experiments where responses cannot be enforced, an intervention might only affect response probabilities to follow-up surveys for some units.
From a discrete choice perspective, our setup admits units who are indifferent to the treatment at any conceivable level of their private information.
From a model selection perspective, unaffected units are equivalent to a sparsity constraint in the selection equation. Modern model selection such as $L_1$-regularization, subset selection, deep neural networks, boosted trees, and random forests are able to recover such potentially sparse structures in the selection equation.
Figure \ref{fig_p0x_INTRO1} contains an example from the job training program analyzed in this paper (Job Corps). It provides the histogram of the estimated relative selection probabilities with and without assignment to training using boosted trees. A value of $1.00$ indicates that assignment does not predict employment status for a given unit.
In this case we have that 2616 out of 9415 units stacking up at an exact[!] numerical one. This pattern is striking and persistent across evaluation periods and estimation methods, see Sections \ref{sec_model} and \ref{sec_empirical1}. In this study, there is significant non-compliance. Thus a large presence of estimated unaffected units seems plausible.
\begin{figure}[!h] \centering\caption{Job Corps Relative Sample Selection Probabilities}
\includegraphics[width = 0.8\textwidth, trim = 0 120 0 100, clip]{figures/JC/INTROv3_plot_p0x_m_3_t_180}
\label{fig_p0x_INTRO1}
\footnotesize \begin{justify}
Histogram of the estimated relative conditional employment selection probabilities under treatment and control. Time period $t=180$, method XGBoost, $n = 9415$. There are 2616 units at exact numerical $1.00$ ($P(\mathcal{X}^0) = 28.61\%$). For more details consider Section \ref{sec_empirical1}.
\end{justify}
\end{figure}
In general, sample selection with unaffected units yield \textit{mixture} distributions for relative conditional selection probabilities obtained in a first stage with point mass at $1.00$. While without harm from an identification perspective, it directly translates into a fundamental problem for estimation and inference. In particular, sharp population bounds are non-smooth functionals for which no uniformly unbiased estimator exists, i.e.~the problem is irregular by construction \citep{hirano2012impossibility}. More concretely, we show that not having a mixture is in fact necessary and sufficient for the population bound to be pathwise differentiable in the sense of \cite{bickel1993efficient}. This means that standard errors and statistical inference obtained from methods such as Lee bounds with ``naive'' asymptotics \citep{lee2009training} are invalid in the presence of unaffected units.
We propose a solution to the problem involving repeated smoothing steps that effectively act as covers for the sharp identified sets. The resulting smooth outer identification regions can be efficiently estimated at the parametric rate and always yield valid confidence regions for the true parameter of interest as the direction of any bias is known. Thus, the degree of smoothing can be chosen to maximize power in finite samples.
Under regularity, the smooth bounds can converge to the sharp identified set. On a high level, our results reveal that there is an identification-precision trade-off: Aiming for an outer identification region with high precision over a sharp region with low precision can significantly improve inference in finite samples.
For both the regular and the smoothed irregular case, we provide efficient influence functions and corresponding estimators using debiased machine learning. Under simple high level conditions, the estimators are asymptotically normal and reach their respective semiparametric efficiency bounds. The influence function for the intensive margin is equal to the one implied by the moment estimator in \cite{semenova2023generalized} under \textit{unknown} propensity scores and no moment selection. To our knowledge, this is the first paper to provide efficient influence functions that can be used for semiparametric estimation of compliers, defiers, and the combined extensive margin bounds in the sample selection model without exclusion.
In contrast to any regular bounds, the smooth counterparts require weaker convergence assumptions for the involved nuisance parameters and, even in the irregular case, no margin condition that controls the distribution of selection effects.
We also explore the role of the propensity score for efficiency, extending the results in \cite{hahn1998role} to sample selection without exclusion. The efficiency bound for the principal strata treatment effect bounds do not depend on the knowledge of the propensity score.
We quantify the efficiency gap when using the moment functions under knowledge of the propensity score as suggested by \cite{semenova2023generalized} and \cite{heiler2024heterogeneous}. These estimators are generally inefficient except in knife-edge cases. However, there is a trade-off from a practical perspective: Analogously to inverse probability weighting (IPW) estimation of average treatment effects, these estimators do not require conditional (truncated) outcome means as additional nuisance inputs. Thus, depending on the difficulty in estimating the latter, using these simpler moment functions might be preferable if propensity scores are known. All aforementioned properties regarding (ir)regularity and smoothing also directly apply to these inefficient moment functions and their smoothed counterparts.
Monte Carlo simulations suggest that, in irregular designs, inference using the smooth methods compares favorably to the alternatives from the literature that either ignore the irregularity of the problem or rely on trimming \citep{heiler2024heterogeneous}, or moment selection/switching methods \citep{andrews2010inference,semenova2023generalized}.
We extend the theory to other quantile-trimmed bounds obtained under additional stochastic dominance assumptions \citep{zhang2003estimation,huber2015sharp}. The components of our influence functions can also be used to construct heterogeneous treatment effect bounds in the sense of \cite{heiler2024heterogeneous} for any principal strata.
The empirical study is a comprehensive re-evaluation of the National Job Corps Study \citep{burghardt1999national,schochet2008does}. We demonstrate the empirical relevance of flexible estimation via machine learning methods, in particular with regards to sparsity in the selection equation. Aggregate results of the new semiparametrically efficient smoothing bounds with optimized learners rule out moderate to large negative impact of Job Corps on hourly earnings at the end of the evaluation period. The results over various time periods are surprisingly close to the more restrictive impact analysis by \cite{lee2009training} that relies on a stronger monotonicity assumption rejected by the data \citep{semenova2023generalized}.
The paper is organized as follows. Section \ref{sec_literature} discusses the relevant literature. Section \ref{sec_model} introduces the nonparametric sample selection model and provides an in-depth discussion of conditional monotonicity. Section \ref{sec_identification} presents the sharp and smooth outer identification regions for all principal strata. Section \ref{sec_regularity} provides the results on pathwise differentiability. Section \ref{sec_estimationinference} contains the assumptions for debiased machine learning estimation and asymptotic inference as well as the efficiency gap under known propensity scores. Section \ref{sec_extensions} extends the methodology to stochastic dominance bounds and heterogeneous effect bounds.
Section \ref{sec_empirical1} presents the empirical study. Section \ref{sec_conclusion} concludes.
\section{Literature} \label{sec_literature}
There is a large literature on sample selection models without exclusion under restrictive parametric assumptions \citep{heckman1979sample,staub2014tobit}. Bounds for causal effects under varying weaker assumptions also have a rich tradition, see \cite{molinari2020microeconometrics} for a comprehensive survey. We focus on research that does not rely on the use of additional exclusion restrictions or instruments and is connected to the monotonicity assumption. We follow a principal stratification approach that considers identification of causal effects for all latent subgroups characterized by their potential selection behavior as a function of a binary treatment \citep{frangakis2002principal}. This circumvents comparing systematically different groups when conditioning on realized selection behavior.
Bounds on causal effects for the intensive margin or always-takers under monotonicity have first been introduced by \cite{zhang2003estimation}, see also \cite{lee2009training} for an extension to a form of conditional monotonicity that excludes units whose selection behavior is unaffected by the treatment. Such ``Lee'' or ``Zhang-Rubin-Lee'' bounds \citep{andersen2023guide}
are defined by conditional expectations trimmed at conditional quantiles. \cite{imai2008sharp} shows the sharpness of such trimming bounds. \cite{semenova2023generalized} introduces orthogonal moments for estimation of intensive margin bounds that can be used in high-dimensional setups, see also \cite{olma2021nonparametric} for nonparametric estimation of truncated conditional expectations. \cite{huber2015sharp} derive bounds for compliers under strong monotonicity. We extend their ideas allowing for compliers and defiers simultaneously.
\cite{honore2020selection,honore2022sample} also discuss bounds in sample selection models without exclusion restrictions under additional parametric assumptions. Our bounds do not require any parametric structure.
\cite{bartalotti2021identifying} provide bounds using monotonicity and stochastic dominance within a marginal treatment effect (MTE) framework. They focus on strong monotonicity only and do not provide estimation or inference theory.
\cite{heiler2024heterogeneous} considers heterogeneous intensive margin bounds and misspecification robust inference under a strong margin assumption effectively ruling out subgroups for which selection probabilities do not depend on the treatment. Our new moment functions can also be used in conjunction with the approach by \cite{heiler2024heterogeneous} to provide heterogeneous effect bounds for any margin.
\cite{okamoto2023bounds} discusses sensitivity analysis of intensive margin and MTE bounds under partial violation of the conventional monotonicity assumption. His stochastic monotonicity assumption is distinct from our approach as it assumes strong monotonicity for a limited share of the population instead of a weaker form of monotonicity.
The use of intensive margin bounds has also been advocated by \cite{chen2023logs} when encountering dependent variables with (many) zeroes. All these papers consider either strong monotonicity or the restrictive version of weak monotonicity \citep{lee2009training}, with \cite{semenova2023generalized} being a notable exception.
Monotonicity also plays a crucial role for instrumental variable (IV) based methods \citep{imbens1994identification,angrist1995two,angrist2000fish,abadie2003semiparametric,heckman2007econometric2,sloczynski2020should,heiler2022efficient}. Similar to the sample selection setup, it is an assumption on latent subtypes in the population to react towards a change in the instrument weakly in a particular direction in terms of their treatment selection. It rationalizes the same choices as a structural model for treatment selection that is additive separable in observables and unobservables both in the population \citep{vytlacil2002independence} as well as numerically if estimated nonparametrically \citep{kline2019heckits}. This also applies to our setup. For nonparametric identification of local average treatment effects or other MTE-type parameters, the instrument has to be able to move certain units in terms of their treatment choice (first stage). Combining compliers and defiers for IV using a weak monotonicity assumption has been considered by \cite{kolesar2013estimation} and \cite{sloczynski2020should}.
The key distinction to the IV setup is that we are interested in the causal effect of the treatment accounting for endogenous selection. Thus, our treatment takes on the role of the instrument while the selection is fully endogenous as the treatment in the IV case.
The weak monotonicity assumption posited in this paper allows for the presence of all latent subtypes in the population but limits their presence within ex-ante unknown partitions of the covariate space. This can also by motivated by a structural selection model with mixed indices within these partitions. The IV framework does not allow for nonparametric identification in any population without movers (no first stage). We explicitly want to allow for units whose selection is not affected by the treatment. This can then generates the mixture distribution for the ``first stage'' relative conditional selection probabilities with point mass at one.
Such mixtures imply that the (generalized) ZR-Lee bounds are no longer smooth in the underlying data distribution and thus there exist no estimator sequence that is locally asymptotically unbiased in the sense of \cite{hirano2012impossibility}.
We construct a particular cover for the identification region that circumvents the non-regularity of the underlying functional of interest. Our identified set has non-empty interior and a zero duality gap as discussed in \cite{kaido2014asymptotically}. Thus, its semiparametric efficiency bound can be characterized and $\sqrt{n}$-consistent regular estimators of the support function exist in a uniform sense. As a by-product of the strict convexity of the smooth identified set, confidence intervals for the effects (not bounds) can be made more precise using modified critical values \citep{imbens2004confidence}. Under a strong margin assumption, our bounds converge converge to the standard sharp ZR-Lee bounds at a known rate.
Our estimands are ratios of a weighted conditional expectations normalized by their share. For the intensive margin, estimating the share is related to estimating an expected conditional outcome under unconfoundedness when setting the treatment status to the best conditional mean. For the latter in isolation, \cite{luedtke2016statistical} show that for the estimand to be pathwise differentiable, it is necessary and sufficient that units which are indifferent between treatments have an \textit{exceptional} law with zero conditional variance under both regimes. We extend their approach to a class of ratio estimands and show that, for the principal strata bounds, the point mass phenomenon is both necessary and sufficient for irregularity, i.e.~there are no \textit{exceptional} laws as the latter would imply a violation of the overlap assumption.
\cite{levis2023covariateassisted} also suggest smoothing for obtaining regular estimators of \cite{balke1997bounds} bounds in the presence of covariates. They propose confidence intervals with a worst-case bias correction arising from the particular approximation. Unlike their approach, our smoothing will always widen the identified set and thus one can navigate the identification-precision trade-off without the need for a worst-case bias component.
\cite{semenova2023generalized} uses the same weak monotonicity assumption as this paper but focuses exclusively on always-takers and does not address regularity and efficiency. To address the problem of misclassifying units, she uses a weak margin assumption and a shrinkage/moment-switching approach. The latter selects between different moment functions based on the estimated selection probabilities with a vanishing sequence of shrinkage parameters that obey a rate condition. There is no guidance on how to choose this parameter in finite samples.
Moreover, while asymptotically valid in handling potential misclassification, this moment selection effectively narrows the identified set leading to potential undercoverage in finite samples. Our approach guarantees at least nominal coverage and can be optimized with respect to power.
For always-taker bounds, \cite{heiler2024heterogeneous} suggests to round trimming threshold towards the closest value on a grid. While this approach similarly yields an outer identification region in the context of relative selection probabilities very close to one, it is unclear how to use in the case of a point mass for units whose selection is unaffected by the treatment, i.e.~an exact one, as rounding could be done in two directions leading to different impact of misclassification errors.
In other contexts, \cite{lee2021bounding} and \cite{pakel2023bounds} also note that there can be a trade-off between identification and precision. They argue that confidence sets obtained from outer bounds can well be tighter compared to using sharp bounds. A similar intuition applies for our smooth outer identification region. For a data combination problem, \cite{dHault2022partially} also consider outer bounds whose conservativeness arise from regularization.
We characterize the efficient influence functions and semiparametric efficiency bounds with and without knowledge of the propensity score. Our results generalize \cite{hahn1998role} to sample selection under weak conditional monotonicity and either (i) regular sharp bounds or (ii) regular smooth bounds that cover the irregular sharp bounds. We also contribute and make use of the expanding literature on the use of debiased machine learning methods for estimation of causal effects in economics \citep{chernozhukov2018double,chernozhukov2022locally} and its combination with partial identification \citep{heiler2024heterogeneous,semenova2023set,semenova2023generalized}. More specifically, we show that the particularities of the relative selection probabilities and the rate at which they can be learned from data crucially interact with the (ir)regularity of the parameters of interest.
\section{Model and Monotonicity} \label{sec_model}
Assume we observe iid data $W = (YS,S,D,X)'$ where $S \in \{0,1\}$ is a selection variable indicating whether an outcome is observed or not. $Y \in \mathcal{Y}$ is the partially unobserved outcome of interest. $D \in \{0,1\}$ is a binary treatment of interest. $X \in \mathcal{X}$ is a vector of predetermined covariates. We are interested in the causal effect of the treatment on the outcome defined in terms of potential outcomes and potential selection indicators. The realized but partially unobserved outcome is defined as \begin{align}
Y = DY(1) + (1-D)Y(0).
\end{align}
Selection is connected to potential selection indicators as \begin{align}
S = DS(1) + (1-D)S(0).
\end{align}
\begin{figure}[!h] \centering \caption{Nonparametric Sample Selection Model} \label{fig_dagSS1}
\begin{tikzpicture}[scale=1]
\tikzset{line width=1.5pt}
\node[ellipse,draw,line width = 1.2pt,dashed, drop shadow, fill = white] (-1) at (4.8,0) {\scriptsize$U$};
\node[ellipse,draw,line width = 1.2pt, drop shadow, fill = white] (0) at (2.4,1.2) {\scriptsize$X$};
\node[ellipse,draw,line width = 1.2pt, drop shadow, fill = white] (1) at (0,0) {\scriptsize$D$};
\node[ellipse,draw,line width = 1.2pt, drop shadow, fill = white] (2) at (3.6,-1.2) {\scriptsize$S$};
\node[ellipse,draw,line width = 1.2pt, drop shadow, fill = lightgray] (3) at (1.2,-1.2) {\scriptsize$Y$};
\node[ellipse,draw,line width = 1.2pt, drop shadow, fill = white] (4) at (2.4,-2.4) {\scriptsize$Y\times S$};
\path (1) edge (3);
\path (2) edge (4);
\path (3) edge (4);
\path (0) edge[out=west,in=north west] (1);
\path (0) edge (2);
\path (0) edge (3);
\path[<->] (0) edge[out=east,in=north west, dashed] (-1);
\path (-1) edge[dashed] (2);
\path (-1) edge[dashed] (3);
\path (1) edge (2);
\end{tikzpicture}\footnotesize \begin{justify}
A graphical representation of a nonparametric sample selection model. Nodes denote variables and edges are structural relationships. A missing arrow from one node to another is an exclusion restriction. Unobserved independent components at each node are omitted. Unobserved variables and its edges are dashed. Gray nodes are partially unobserved.
\end{justify} \vspace{-12pt}
\end{figure}
Figure \ref{fig_dagSS1} depicts a prototypical sample selection model with the relevant exclusion restrictions. Such a model implies the following conditional independence assumption for the treatment which we maintain throughout:
\begin{ass}[Conditional Independence] \label{ass_CIA}
$$D \perp\!\!\!\perp Y(d),S(d)~|~X=x \textit{ for all } x \in \mathcal{X} \textit{ and } d=0,1,$$
\end{ass}
where $\perp\!\!\!\perp$ denotes statistical independence. Thus, treatment is assumed to be exogenous conditional on $X$. However, selection is left essentially unrestricted, i.e.~it can be fully endogenous or exogenous.
In particular, there are no exclusion restrictions available for $S$, i.e.~we cannot use instruments or similar to account for endogenous selection.
The potential selection indicators define principal strata which represent how the treatment affects selection behavior. Table \ref{tab_principalStrata1} contains all strata using the nomenclature from the instrumental variables literature.
\begin{table}[!h]
\centering
\caption{Principal Strata}
\label{tab_principalStrata1}
\begin{tabular}{llcc} \hline \hline
Type & &$S(0)$ & $S(1)$ \\ \hline \\[-1.8ex]
(AT) &Always-takers & 1 & 1 \\
(C) &Compliers & 0 & 1 \\
(D) &Defiers & 1 & 0 \\
(NT) &Never-takers & 0 & 0 \\ [0.5ex] \hline \hline
\end{tabular}
\end{table}
To proceed, we postulate the existence of an unknown but identified partitioning of the covariate space $\mathcal{X}$ on which a weak monotonicity assumption applies to these principal strata\footnote{Note that this definition differs slightly from the conditional monotonicity assumption of \cite{semenova2023generalized} which does not yield a unique partitioning.}.
\begin{ass}[Weak/Conditional Monotonicity]\label{ass_monotone}
Assume there are subsets of the covariate space $\tilde{\mathcal{X}}^+,\tilde{\mathcal{X}}^-\subseteq \mathcal{X}$ such that \begin{align*}
S(1) \geq S(0) &\textit{ if } x \in \tilde{\mathcal{X}}^+, \\
S(1) \leq S(0) &\textit{ if } x \in \tilde{\mathcal{X}}^-.
\end{align*} \label{ass_weakMon1}
\end{ass}
\vspace{-2em}
Assumption \ref{ass_weakMon1} yields a distinct partitioning of $\mathcal{X} = {\mathcal{X}}^+ \cup {\mathcal{X}}^- \cup {\mathcal{X}}^0$ where \begin{align*}
\mathcal{X}^0 &= \tilde{\mathcal{X}}^+ \cap \tilde{\mathcal{X}}^-, \\
\mathcal{X}^+ &=\tilde{\mathcal{X}}^+ \backslash \mathcal{X}^0, \\
\mathcal{X}^- &=\tilde{\mathcal{X}}^- \backslash \mathcal{X}^0.
\end{align*}
The presence of a potentially non-empty $\mathcal{X}^0$, i.e. $P(\mathcal{X}^0) > 0$, will be crucial in what follows.
Substantially, weak monotonicity allows for the presence of all principal strata including defiers. However, it is a partial restriction of types within partitions. Table \ref{tab_partitions1} contains the mixture of types in the different partitions as well the population.
\begin{table}[!h]
\centering \footnotesize
\caption{Partitions, Observed Strata, and Admitted Latent Types}
\label{tab_partitions1}
\begin{minipage}{\linewidth}
\begin{subtable}{\linewidth}
\centering
\vspace{0.8em}
\begin{tabular}{l|cc}
$S~\backslash~ D$ & 0 & 1 \\[0.5ex] \hline \\[-1.5ex]
0 & NT/C & NT/D \\[1ex]
1 & AT/D & AT/C \\[0.5ex] \hline
\end{tabular}
\caption{$x\in \mathcal{X}$}
\vspace{1em}
\end{subtable}
\end{minipage} \\%
\begin{minipage}{0.3\linewidth}
\centering
\begin{subtable}{\linewidth}
\centering
\begin{tabular}{l|cc}
$S~\backslash~ D$ & 0 & 1 \\[0.5ex] \hline \\[-1.5ex]
0 & NT/C & NT \\[1ex]
1 & AT & AT/C \\[0.5ex] \hline
\end{tabular}
\caption{$x\in \mathcal{X}^+$}
\end{subtable}
\end{minipage}
\begin{minipage}{0.3\linewidth}
\centering
\begin{subtable}{\linewidth}
\centering
\begin{tabular}{l|cc}
$S~\backslash~ D$ & 0 & 1 \\[0.5ex] \hline \\[-1.5ex]
0 & ~~NT~~ & ~~NT~~ \\[1ex]
1 & AT & AT \\[0.5ex] \hline
\end{tabular}
\caption{$x\in \mathcal{X}^0$}
\end{subtable}
\end{minipage}
\begin{minipage}{0.3\linewidth}
\centering
\begin{subtable}{\linewidth}
\centering
\begin{tabular}{l|cc}
$S~\backslash~ D$ & 0 & 1 \\[0.5ex] \hline \\[-1.5ex]
0 & NT & NT/D \\[1ex]
1 & AT/D & AT \\[0.5ex] \hline
\end{tabular}
\caption{$x\in \mathcal{X}^-$}
\end{subtable}
\end{minipage}
\end{table}
For example, in partition $\mathcal{X}^0$, the population consists only of always-takers and never-takers while on $\mathcal{X}^+$ and $\mathcal{X}^-$ only defiers or compliers are ruled out, respectively. This assumption seems most relevant if we have some discrete covariates and/or continuous variable that affect selection behavior discontinuously e.g.~by threshold effects.
Assumption \ref{ass_weakMon1} is also consistent with a single-index structure in the selection equation that allows for sparsity with respect to the treatment on some parts of the covariate space.
\paragraph{\textbf{Example (Single Index Model):}} \label{sec_singleIndex1}
Consider a simple semiparametric selection model with treatment and selection that are additive separable in observables and unobservables
\begin{align}
D &= \mathbbm{1}(f(X) + U_D \geq 0), \notag \\
Y &= \mu(X,D) + U_Y, \notag \\
S &= \mathbbm{1}(g(D,X) + U_S \geq 0),
\end{align}
where $U_D \perp\!\!\!\perp U_Y,U_S$ conditional on $X=x$ then implies conditional independence as in Assumption \ref{ass_CIA}.
Weak monotonicity as in Assumption \ref{ass_monotone} here is implied by an unrestricted $g(D,X)$. This can generate a mixed index structure \begin{align}
g(d,x) = \begin{cases} g_+(d,x) &\textit{ if } x \in \mathcal{X^+},\\
g_-(d,x) &\textit{ if } x \in \mathcal{X^-},\\
g_0(x) &\textit{ if } x \in \mathcal{X}^0.\end{cases}
\end{align}
where $g_+(1,x) > g_+(0,x)$ and $g_-(0,x) > g_-(1,x)$. The fact that $g_0(x)$ does not depend on $d$ can be seen as a sparsity assumption within a partition of the covariate space. For example, consider the case of a simple parametric index in the selection equation with $X = (X_1,X_2)'$ where $X_1$ is continuous and $X_2$ discrete with support points $\{-1,0,1\}$
\begin{align}
S = \mathbbm{1}\big(\gamma_0 + \gamma_1X_1 + \gamma_2D\mathbbm{1}(X_2 = -1) + \gamma_3D\mathbbm{1}(X_2 = 1) + U_S \geq 0\big).
\end{align}
Here $X_2$ can separate the monotonicity types. In particular, $\gamma_2 < 0 < \gamma_3$ or $\gamma_2 > 0 > \gamma_3$ imply conditional monotonicity. If the signs of $\gamma_2$ and $\gamma_3$ are identical, then either $\mathcal{X}^+ = \{ \emptyset\} $ or $\mathcal{X}^- = \{ \emptyset\}$, implying the conventional \textit{unconditional} or \textit{strong} monotonicity $P(S(1) \geq S(0)) = 1$.
Whenever $P(\mathcal{X}^0) > 0$, conditional monotonicity and independence generate a \textit{mixture} distribution for the potential or relative conditional selection probabilities
\begin{align}
p_0(x) = \frac{P(S=1|D=0,X=x)}{P(S=1|D=1,X=x)} = \frac{P(S(0)=1|X=x)}{P(S(1)=1|X=x)}.
\end{align}
In particular, for $\mathcal{X}^+$ and $\mathcal{X}^-$ we have that $p_0(x) < 1$ and $p_0(x) > 1$, respectively, while there is potentially a point mass for which $p_0(x) = 1$ corresponding to $P(\mathcal{X}^0)>0$. Thus, conditional monotonicity allows for mixture distributions which are not necessarily smooth in $p_0(x)$. We seek to conduct inference on causal effects for all principal strata under such mixtures.
\begin{figure}[!h]
\centering
\caption{Example Distributions for $\hat{p}_0(x)$ at $t = 208$}
\label{fig:hist_all_208}
\begin{subfigure}{0.45\textwidth}
\caption{\footnotesize Logit}
\begin{flushright}
\includegraphics[width=0.66667\textwidth, trim =0 10 0 0 , clip]{figures/JC/v3_plot_p0x_m_1_t_208.pdf} \end{flushright}
\end{subfigure}
\begin{subfigure}{0.45\textwidth}
\caption{ \footnotesize Logit Interacted}
\begin{flushleft} \includegraphics[width=0.66667\textwidth, clip]{figures/JC/v3_plot_p0x_m_2_t_208.pdf} \end{flushleft}
\end{subfigure} \\
\begin{subfigure}{0.3\textwidth}
\caption{\footnotesize XGBoost}
\begin{flushright} \includegraphics[width=\textwidth, clip]{figures/JC/v3_plot_p0x_m_3_t_208.pdf} \end{flushright}
\end{subfigure}
\begin{subfigure}{0.3\textwidth}
\caption{\footnotesize Neural Network}
\begin{flushleft} \includegraphics[width=\textwidth, clip]{figures/JC/v3_plot_p0x_m_4_t_208.pdf} \end{flushleft}
\end{subfigure}
\begin{subfigure}{0.3\textwidth}
\caption{\footnotesize Random Forest}
\begin{flushleft} \includegraphics[width=\textwidth, clip]{figures/JC/v3_plot_p0x_m_5_t_208.pdf} \end{flushleft}
\end{subfigure}
\begin{justify} \footnotesize
The figure contains histograms for the estimated $p_0(x)$ using five different estimation methods on the research sample at week $t=208$. (a) Logit is a logistic regression with linear additive index. (b) Logit Interacted is a fully treatment interacted logistic model. (c) XGBoost are gradient boosted trees with cross-validated boosting steps using \verb|xgboost|. (d) Neural Network is a batch-trained, fully connected artificial feed-forward neural network with three ReLu hidden layers, drop-out regularization, and a final sigmoid layer using \verb|keras|. (e) Random Forest is an honest probability forest using \verb|grf| with default parameters.
\end{justify}
\end{figure}
Modern model selection methods such as $L_1$-regularization, subset selection, neural networks, boosting, and random forests are able to recover and produce potentially sparse structures in the selection equation. We present an example from relative employment probabilities and a job training treatment in Figure \ref{fig:hist_all_208}. It contains the predicted relative selection probabilities from Section \ref{sec_empirical1} using five different estimation methods: (a) logistic regression, (b) logistic regression with fully interacted treatment, (c) gradient boosted trees (XGBoost), (d) neural network, and (e) random forest, see Section \ref{sec_empirical1} and Appendix \ref{app_JC} for more details.
The methods capable of sparsity yield a clear mixture pattern for the estimated distribution of $p_0(x)$. In particular, (c), (d), and (e) produce an exact share of numerical ones with up to $45.16\%$ of the sample estimates. (a) imposes strong monotonicity and cannot produce exact ones in finite samples by construction. We suspect that (a) is most likely to be misspecified. (b) has relevant density around one, but also generates overly extreme values for the relative selection probabilities. The sparse methods agree on around 45\% -- 60\% of the close or equal to one classifications. This qualitative pattern can be observed along a long sequence of post-treatment periods. Over time there is an average distribution shift suggesting more positive employment effects. The sparsity pattern, however, is pervasive. We take this as a signal that the mixture distributions are likely reflecting of or a good approximation to an underlying structure and are not only spurious by-products of particular estimation algorithms motivating Assumption \ref{ass_weakMon1} with $P(\mathcal{X}^0) > 0$.
\paragraph{Remark:}
On a cautionary note, empirical distributions of individual predictions such as in Figure \ref{fig:hist_all_208} should not be over-interpreted. On one hand, we can have clear theoretical reasons why true $p_0(x)$ ought to be one for many units, e.g.~units that did not intend to join the program are likely to be unaffected by a non-binding assignment to treatment. On the other hand, these are just finite-sample estimates from models that are all likely misspecified to a smaller or larger extend. Extreme values that suggest economically unreasonable unit selection effects are likely driven by a high-variance model.
We can take a more agnostic perspective with regards to the evidence provided by the empirical distribution of $\hat{p}_0(x)$. In particular, in high dimensions or with complicated functional forms, convergence to true relative selection probabilities is hard to achieve for any model. In finite samples, producing a sparse representation by removing the impact of the treatment in parts of the covariate space could just be a method's solution to optimize fit. Methods for predicting conditional selection probabilities are not necessarily designed to best predict the implied relative $p_0(x)$. However, if a learner performs best among conventional performance metrics for probabilities and does so by producing sparsity with respect the treatment, we treat it as a potentially better approximation and need to deal with the point mass regardless of whether this is a true or approximate model.
\section{Identification of Sharp and Smooth Bounds} \label{sec_identification}
\subsection{Sharp Bounds}
\subsubsection{The General Case}
We denote $a = O(b)$ and $a=O_p(b)$ as $a\lesssim b$ and $a\lesssim_P b$, respectively. We first present identification of conditional causal effect bounds for any principal stratum, i.e.~parameters of the form \begin{align}
\beta(x,s_0,s_1) = E[Y(1)-Y(0)|S(0)=s_0,S(1)=s_1,X=x].
\end{align}
This covers all relevant components for any of the principal strata and margin effects.\footnote{Trivially, for a principal stratum effect to be well defined, we require that the corresponding population share is non-zero. We assume this throughout.}
Let $q_d(u,x)$ be the $u$-quantile of $Y$ conditional on $S=1, D=d, X=x$. For $d\in\{0,1\}$, we define \begin{align}
\beta_{1,d}(x,u) = E[Y|S=1,D=d,X=x,Y\leq q_d(u,x)], \\
\beta_{0,d}(x,u) = E[Y|S=1,D=d,X=x,Y\geq q_d(u,x)],
\end{align}
which yields
\begin{align}
\beta_{1,d}(x,1) = \beta_{0,d}(x,0) = E[Y|S=1,D=d,X=x].
\end{align}
In line with the literature on monotonicity bounds, we assume a continuous outcome.
\begin{ass}[Continuity]\label{ass_continuity}
For all $x\in \mathcal{X}$ and $d\in \{0,1\}$, the conditional outcome distribution $P(Y \leq y|S=1,D=d,X=x)$ is continuous.
\end{ass}
Moreover, there must be comparable units in terms of covariates in both selected treatment groups. \begin{ass}[Multiple Overlap]\label{ass_overlap} Let $m(x) = P(D=1|X=x)$ and $s(d,x) = P(S=1|D=d,X=x)$. For all $x\in \mathcal{X}$ and any $d \in \{0,1\}$, we have that $0 < m(x) < 1$ and $0 < s(d,x)< 1$.
\end{ass}
Assumptions \ref{ass_CIA} and \ref{ass_monotone} imply the specific bounds in terms of truncated means of observed conditional distributions. Assumptions \ref{ass_continuity} and \ref{ass_overlap} ensure that the corresponding conditioning sets are non-empty and the relevant truncation quantiles are unique. In particular, we obtain the following identification result for the sharp upper and lower bounds for all principal strata.
\begin{prop}[Identification] \label{prop_identification1}
Let $\mathcal{Y}_d(x)$ be the support of $Y(d)$ conditional on $X=x$ with $\overline{y}_d(x)$ and $\underline{y}_d(x)$ its respective infimum and supremum for both $d\in\{0,1\}$. Under Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_continuity}, and \ref{ass_overlap}, the sharp lower bounds for the principal strata are given by \begin{align*}
\beta_L(x,s_0,s_1) = \begin{cases}
\beta_{1,1}(x,\min\{p_0(x),1\}) - \beta_{0,0}(x,1-\min\{1/p_0(x),1\}) &\text{ if } s_0 = 1,~s_1 = 1 \\
\beta_{1,1}(x,1-\min\{p_0(x),1\}) - \beta_{0,0}(x,\min\{1/p_0(x),1\}) &\text{ if } s_0 \neq s_1 \\
\underline{{y}}_1(x) - \overline{{y}}_0(x) &\text{ if } s_0 = 0,~s_1 = 0,
\end{cases}
\end{align*}
and the sharp upper bounds by \begin{align*}
\beta_U(x,s_0,s_1) = \begin{cases}
\beta_{0,1}(x,1-\min\{p_0(x),1\}) - \beta_{1,0}(x,\min\{1/p_0(x),1\}) &\text{ if } s_0 = 1,~s_1 = 1 \\
\beta_{0,1}(x,\min\{p_0(x),1\}) - \beta_{1,0}(x,1-\min\{1/p_0(x),1\}) &\text{ if } s_0 \neq s_1 \\
\overline{{y}}_1(x) - \underline{{y}}_0(x) &\text{ if } s_0 = 0,~s_1 = 0.
\end{cases}
\end{align*}
\end{prop}
For never-takers, all lower and upper bounds are generally uninformative due to never observing them as part of a selected group. For defiers and compliers, the bounds are only informative when at least part of the conditional support is bounded, while for always-takers bounds are also informative without any support condition.
Unconditional bounds for any principal stratum are obtained by integrating the conditional bounds with respect to stratum-specific covariate distributions \begin{align}\label{def_betaB_general}
\beta_B(s_0,s_1) &= \int\beta_B(x,s_0,s_1)dP(x|S(0)=s_0,S(1)=s_1)
\end{align}
for $B\in\{L,U\}$, where \begin{align*}
dP(x|S(0)=s_0,S(1)=s_1) &= \begin{cases}
\frac{\min\{s(0,x),s(1,x)\}}{E[\min\{s(0,X),s(1,X)\}]}dP(x) &\text{ if } s_0 = 1,~ s_1 = 1 \\
\frac{\max\{0, s(1,x)-s(0,x)\}}{E[\max\{0, s(1,X)-s(0,X)\}]}dP(x) &\text{ if } s_0 = 0,~s_1 = 1 \\
\frac{\max\{0, s(0,x)-s(1,x)\}}{E[\max\{0, s(0,X)-s(1,X)\}]}dP(x) &\text{ if } s_1 = 1,~s_1 = 0 \\
\frac{1-\max\{s(0,x),s(1,x)\}}
{E[1-\max\{s(0,X),s(1,X)\}]}dP(x) &\text{ if } s_1 = 0,~s_1 = 0.
\end{cases}
\end{align*}
\subsubsection{Example I: Intensive Margin Bounds}\label{ex_extensive}
The upper and lower conditional intensive margin or always-taker bounds correspond to the ones given by
\cite{lee2009training} and \cite{semenova2023generalized}. Namely, we have that
\begin{equation}
\beta_L(x,1,1)=\beta_{L,1}(x,1,1)-\beta_{L,0}(x,1,1)
\end{equation}
where
\begin{align*}
\beta_{L,1}(x,1,1) &= \begin{cases}
E[Y|S=1,D=1,X=x,Y\leq q_1(p_0(x),x)] &\quad~\quad~\text{ if } x \in \mathcal{X}^+ ,\\
E[Y|S=1,D=1,X=x] &\quad~\quad~\text{ if } x \in \mathcal{X}^- \cup \mathcal{X}^0,
\end{cases} \\
\beta_{L,0}(x,1,1) &= \begin{cases}
E[Y|S=1,D=0,X=x] &\text{ if } x \in \mathcal{X}^+ \cup \mathcal{X}^0, \\
E[Y|S=1,D=0,X=x,Y\geq q_0(1-1/p_0(x),x)] &\text{ if } x \in \mathcal{X}^-. \\
\end{cases}
\end{align*}
Similarly, the conditional upper bounds are given by
\begin{equation}
\beta_U(x,1,1)=\beta_{U,1}(x,1,1)-\beta_{U,0}(x,1,1),
\end{equation}
where
\begin{align*}
\beta_{U,1}(x,1,1) &= \begin{cases}
E[Y|S=1,D=1,X=x,Y\geq q_1(1-p_0(x),x)] &\text{ if } x \in \mathcal{X}^+ \\
E[Y|S=1,D=1,X=x] &\text{ if } x \in \mathcal{X}^- \cup \mathcal{X}^0, \\
\end{cases} \\
\beta_{U,0}(x,1,1) &= \begin{cases}
E[Y|S=1,D=0,X=x] &\text{ if } x \in \mathcal{X}^+ \cup \mathcal{X}^0, \\
E[Y|S=1,D=0,X=x,Y\leq q_0(1/p_0(x),x)] &\text{ if } x \in \mathcal{X}^-.
\end{cases}
\end{align*}
The unconditional bounds for the effect at the intensive margin can then be recovered as \begin{align}\label{def_betaB}
\beta_B(1,1) &= \int\beta_B(x,1,1)dP(x|S(0)=S(1)=1) \notag \\
&= \frac{\int \beta_B(x,1,1)\min\{s(0,x),s(1,x)\}dP(x)}{\int\min\{s(0,x),s(1,x)\}dP(x)}.
\end{align}
\subsubsection{Example II: Extensive Margin Bounds}\label{ex_intensive}
The extensive margin is the combination of compliers and defiers where $s_0 \neq s_1$. Thus, for $B\in\{L,U\},$ the extensive margin effect bounds can be written as {\begin{align}
\beta_B(em)=\int \beta_B(x,em) dP(x|S(0) \neq S(1)).
\end{align}}
where \begin{align*}
dP&(x|S(0) \neq S(1)) \\
&= \frac{P(S(0) \neq S(1)|X=x) dP(x)}{P(S(0) \neq S(1))}\\
&=\frac{\left[P(S(0)=0,S(1)=1|X=x) + P(S(0)=1,S(1)=0|X=x) \right] dP(x)}{\int
\left[P(S(0)=0,S(1)=1|X=x) + P(S(0)=1,S(1)=0|X=x) \right] dP(x)} \\
&= \frac{|s(1,x) - s(0,x)| dP(x)}{\int
|s(1,x) - s(0,x)|dP(x)}.
\end{align*}
Combining compliers and defiers from Proposition \ref{prop_identification1} then yields the conditional lower bound
\begin{align*}
\beta_{L,1}(x,em) = \begin{cases}
E[Y|S=1,D=1,X=x,Y\leq q_1(1-p_0(x),x)] &\text{ if } x \in \mathcal{X}^+ ,\\
\underline{y}_1(x) &\text{ if } x \in \mathcal{X}^- \cup \mathcal{X}^0,
\end{cases}
\end{align*}
and
\begin{align*}
\beta_{L,0}(x,em) = \begin{cases}
\overline{y}_0(x) &\text{ if } x \in \mathcal{X}^+ \cup \mathcal{X}^0, \\
E[Y|S=1,D=0,X=x,Y\geq q_0(1/p_0(x),x)] &\text{ if } x \in \mathcal{X}^-. \\
\end{cases}
\end{align*}
The conditional upper bound $\beta_U(x,em)$ is found analogously from from Proposition \ref{prop_identification1}.
As a result, the unconditional bounds at the extensive margin are given by
\begin{align}
\beta_B(em)=\frac{\int \beta_B(x,em) |s(1,x)-s(0,x)|dP(x)}{\int |s(1,x)-s(0,x)|dP(x)}
\end{align}
\subsection{Smooth Bounds}\label{sec_smooth}
The bounds in Proposition \ref{prop_identification1} are sharp. However, they, as well as the strata densities, contain a variety of non-differentiable components. This poses a challenge to statistical inference. We demonstrate the specific regularity problem and its consequences for inference in Section \ref{sec_regularity}.
We now introduce the smooth bounds that provide an outer identification region to circumvent these problems. In particular, the identification region is constructed to eliminate the impact of misclassifying different monotonicity types while still remaining close to the sharp identified set. We focus on the intensive margin bounds for simplicity but the method works equivalently for all principal strata.
More specifically, the form of the bounds in Proposition \ref{prop_identification1} reveal two sources where classification of the partition will enter, (i) conditional effects bounds $\beta_B(x,1,1)$, and (ii) conditional always-taker share $\min\{s(0,x),s(1,x)\}$. We suggest to apply repeated, asymmetric smoothing to both components. For (i), consider first the conditional bound $\beta_L(x,1,1)$. Note that, for any $x \in \mathcal{X}$, the following parameter is a valid lower bound for $\beta_L(x,1,1)$
\begin{align}\label{smoot-bound-lower-term}
{\beta}_{L,h}(x,1,1) &= E[Y|S=1,D=1,X=x,Y\leq q_1(g_{1,h}(p_0(x)),x )] \notag \\
&\quad - E[Y|S=1,D=0,X=x,Y\geq q_0(1-g_{1,h}(1/p_0(x)),x)] \notag \\
&= \beta_{1,1}(x,g_{1,h}(p_0(x))) - \beta_{0,0}(x,1-g_{1,h}(1/p_0(x)))
\end{align}
where $g_{1,h}(z)$ is a smooth monotonic function indexed by a smoothing parameter $h \in \mathcal{H}_n \subset \mathbb{R}^+$ such that $0 \leq g_{1,h}(z) \leq \min\{z,1\}.$
This ensures that
\begin{equation}
{\beta}_{L,h}(x,1,1) \leq \beta_L(x,1,1).
\end{equation}
Next, consider that the smoothing of the conditional always-taker share is given by
$\min\{s(0,x),s(1,x)\} = s(1,x) \min\{p_0(x),1\}$.
We obtain lower and upper bounds using smooth functions
$g_{1,h}(z)$ and $g_{3,h}(z)$ such that $g_{1,h}(z) \leq \min\{z,1\} \leq g_{3,h}(z).$
In which direction to smooth, depends on the sign of ${\beta}_{L,h}(x,1,1)$. We decompose this into the difference of two non-negative terms via $z=\max\{z,0\}-\max\{-z,0\}.$
We approximate these two functions with a smoothing function $g_{4,h}(z)\leq \max\{z,0\},$ and the earlier $g_{2,h}(-z) \geq \max\{-z,0\}.$
The unconditional smooth lower effect bound is then given by
\begin{align}\label{def_smoothbeta}
\beta_{L,h}(1,1) &=
\frac{E\left[g_{4,h}({\beta}_{L,h}(x,1,1))g_{1,h}(p_0(X))s(1,X)\right]}{E[g_{3,h}(p_0(X))s(1,X)]} \notag \\
&\quad -\frac{E[g_{2,h}(-{\beta}_{L,h}(x,1,1)) g_{3,h}(p_0(X))s(1,X)]}
{E[g_{1,h}(p_0(X))s(1,X)]}.
\end{align}
The aforementioned relations together yield
\begin{equation}
\beta_{L,h}(1,1) \leq \beta_L(1,1).
\end{equation}
Analogously, the conditional smooth upper bound is given by
\begin{align}\label{smoot-bound-upper-term}
\beta_{U,h}(x,1,1)
&= \beta_{0,1}(x,1-g_{1,h}(p_0(x))) - \beta_{1,0}(x,g_{1,h}(1/p_0(x)))
\end{align}
which obeys the inequality $\beta_{U,h}(x,1,1) \geq \beta_U(x,1,1)$, and the unconditional smooth upper bound is defined as
\begin{align}\label{def_smoothbeta_upper}
\beta_{U,h}(1,1) &=
\frac{E\left[g_{2,h}({\beta}_{U,h}(x,1,1))g_{3,h}(p_0(X))s(1,X)\right]}{E[g_{1,h}(p_0(X))s(1,X)]} \notag \\
&\quad -\frac{E[g_{4,h}(-{\beta}_{U,h}(x,1,1)) g_{1,h}(p_0(X))s(1,X)]}
{E[g_{3,h}(p_0(X))s(1,X)]},
\end{align}
which analogously satisfies the reversed inequality
\begin{equation}
\beta_{U,h}(1,1) \geq \beta_U(1,1).
\end{equation}
Next, we propose a class of functions that satisfy the imposed conditions and show that the distance between smooth and sharp unconditional bounds is negligible as $h \to 0$ with bounded support.
\begin{ass}[Outcome Regularity I]\label{ass_bddmoment} For $d\in\{0,1\}$, the outcome distribution $P(Y \leq y| S=1, D=d, X=x)$ has uniformly bounded support
, i.e. $\sup_{x\in \mathcal{X}} (|\underline{y}_d(x)|+|\overline{y}_d(x)|)<\infty.$
\end{ass}
We also impose the following assumption about approximation functions.
\begin{ass}[Approximation Functions]\label{ass_errors-of-g}
Approximation functions $g_{i,h}, i=1, \ldots, 4$ satisfy for all $z \in \mathbb{R}$
\begin{align*}
g_{1,h}(z) \leq \min\{z,1\} \leq g_{3,h}(z) \quad \mbox{ and } \quad g_{4,h}(z)\leq \max\{z,0\} \leq g_{2,h}(z),
\end{align*}
and for each $i=1, \ldots, 4$
\begin{equation*}
\sup_{z \in \mathbb{R}} |g_{i,h}(z)-g_i(z)|
\lesssim h,
\end{equation*}
where $g_1(z)=g_3(z)=\min\{z,1\}$ and $g_2(z)=g_4(z)=\max\{z,0\}.$
In addition, $g_{1,h}(z) \geq 0$ for all $z \in \mathbb{R}$ and $h \in \mathcal{H}_n.$\footnote{All quantile trimming thresholds are well defined as long as $\inf_{x \in \mathcal{X}} g_{1,h}(p_0(x)) \geq 0$ and $\inf_{x \in \mathcal{X}} g_{1,h}(1/p_0(x)) \geq 0.$ Given the strong multiple overlap assumption, one can always find an upper bound for $h$ such that this applies. These large-$h$ values are not relevant from a practical perspective and we consider $h$ to be below such a ceiling in what follows.}
\end{ass}
There exist many functions with these properties. Figure \ref{fig_smoothapprox1} depicts an example using a LogSumExp function for smoothing via $g_{1,h}(z)$.\begin{figure}[!h]
\centering
\caption{Smooth Approximation: $g_{1,h}(z)$ Example}
\label{fig_smoothapprox1}
\includegraphics[trim = 0 130 0 110, width = 0.6\textwidth]{figures/plot_approx2.pdf}
\footnotesize \begin{justify}
The figure contains $\min\{p_0(x),1\}$ as well as two smooth approximations using $g_{1,h}(z) = 1 - h\log(1 + \exp(-(z - 1)/h))$ for $h=0.30$, $h=0.15$, and $h=0.0015$ (Asy). The size of the $h$-cover at $p_0(x) = 1$ determines how close the bounds are to the point identified effect for $x \in \mathcal{X}_0$. The $h$-cover for units with $p_0(x) \neq 1$ introduces an additional distance compared to the sharp trimming bounds.
\end{justify}
\end{figure}
We obtain the following Theorem.
\begin{theorem}\label{thm:convergence}
Let $\mathcal{P}$ be the set of probability measures satisfying Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_continuity}, \ref{ass_overlap}, \ref{ass_bddmoment}, and \ref{ass_errors-of-g}. Denote $\underline{g}_n = \inf_x g_{1,h}(p_0(x))$. Then, for $B \in \{L, U\},$ and $(s_0, s_1) \neq (0,0),$ we have that
\begin{equation*}
\sup_{P \in \mathcal{P}} |\beta_{B,h}(s_0,s_1) -\beta_B(s_0, s_1)|\lesssim \frac{h}{\underline{g}_n^2}.
\end{equation*}
\end{theorem}
Theorem \ref{thm:convergence} gives the approximation error of the smooth bounds. Under the strong overlap condition, $\underline{g}_n$ is bounded away from zero and thus the difference decays linearly in $h$. Without strong overlap, there is an implicit restriction on the rate at which the smooth relative probabilities $g_{1,h}(p_0(x))$ can approach zero.
\section{Regularity and Semiparametric Efficiency Bounds} \label{sec_regularity}
In this section, we present our main results on regularity and semiparametric efficiency. Again, we focus on the lower bound for the intensive margin for simplicity.
Table \ref{tab_EIFs1} provides the influence function for the always-taker lower bound and its smoothed counterpart. All other influence functions are in Appendix \ref{sec_proofs}. We obtain the following theorem.
\begin{ass}[Outcome Regularity II]\label{ass_density} For all $x\in \mathcal{X}$ and $d \in \{0,1\},$ the conditional outcome distribution $P(Y \leq y| S=1, D=d, X=x)$ has finite support on $[\underline{y}_d, \overline{y}_d]$ and a continuous density $f(y|X=x, D=d, S=1)$ bounded from below and above.
\end{ass}
\begin{ass}[Strong Multiple Overlap]\label{ass_strong_overlap} There exist constants $\underline{m},\underline{s} \in (0,1/2)$ such that \begin{align*}
\underline{m} < \inf_{x\in\mathcal{X}} m(x) \leq \sup_{x\in\mathcal{X}} m(x) < 1- \underline{m},
\end{align*}\begin{align*}
\underline{s} < \inf_{x\in\mathcal{X},d\in{\{0,1\}}} s(d,x) \leq \sup_{x\in\mathcal{X},d\in{\{0,1\}}} s(d,x) < 1- \underline{s}.
\end{align*} \label{ass_STRONGoverlap}
\end{ass}
Assumptions \ref{ass_density} and \ref{ass_strong_overlap} are stronger versions of continuity and overlap in Assumptions \ref{ass_continuity} and \ref{ass_overlap}, respectively.
Without strong overlap, there is an additional irregularity problem for the population bounds equivalently to irregular identification of average treatment effects under unconfoundedness with many extreme propensities scores \citep{khan2010irregular,HEILER2021valid}. We abstract from such issues in this paper to focus on the irregularity obtained from $P(\mathcal{X}^0) > 0$.
\begin{theorem}\label{thm_differentiable}
Suppose that Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_density}, and \ref{ass_strong_overlap} hold. For $B \in \{L, U \},$ $\beta_B(1,1)$, $\beta_B(0,1)$, $\beta_B(0,1),$ $\beta_B(\mbox{em})$ are pathwise differentiable if and only if $P(\mathcal{X}^0)=0.$
\end{theorem}
Theorem \ref{thm_differentiable} characterizes that the necessary and sufficient condition for the unconditional lower bound introduced in \eqref{def_betaB} to be pathwise differentiable is the absence of the point mass on the set $\mathcal{X}^0$, i.e.~$P(\mathcal{X}^0)=0$.\footnote{
When $P(\mathcal{X}^0)>0$, the result in Theorem \ref{thm_differentiable} that $\beta_L(1,1)$ is not pathwise differentiable is related but different to \cite{luedtke2016statistical}. Their target is essentially the denominator of $\beta_L(1,1)$ in \eqref{def_betaB}. However, divergence of the derivative of denominator and/or divergence of the numerator separately does not necessarily imply divergence of the ratio. Moreover, the numerator includes the conditional treatment effect bounds $\beta_L(1,1,x)$ and consequently conditional quantiles $q_d(u,x)$ which changes the analysis. } We also present the semiparametric efficiency bound for $\beta_L(1,1)$. It does not depend on the knowledge of the propensity score.
\begin{corr}\label{corr_var_bounds}
Let $1_{\mathcal{X}^{+}} = \mathbbm{1}(X \in \mathcal{X}^+)$ and $1_{\mathcal{X}^{-}} = \mathbbm{1}(X \in \mathcal{X}^-)$. Suppose that Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_density} and \ref{ass_strong_overlap} hold, and $P(\mathcal{X}^0)=0.$ The semiparametric efficiency bound $Var_a(\beta_L(1,1))$ for $\beta_L(1,1)$ is defined by
{\footnotesize\begin{align*}
&E[\min(s(0,X), s(1,X)]^2Var_a(\beta_L(1,1)) =
E\left[\frac{s(1,X) \sigma_1^2 (X) }{m(X)}+\frac{s(0,X) \sigma_0^2 (X) }{1-m(X)} \right] \\
&+E\left[1_{\mathcal{X}^{+}}(\beta_L(X,1,1)-\beta_L(1,1))^2 \frac{s(0,X)(1-s(0,X) m(X))}{1-m(X)} \right] \\
&+E\left[1_{\mathcal{X}^{-}}(\beta_L(X,1,1)-\beta_L(1,1))^2 \frac{s(1,X)(1-s(1,X)+ s(1,X) m(X))}{m(X)} \right] \\
&+E\left[ 1_{\mathcal{X}^{+}} \frac{s(1,X) q_1(p_0(X),X)^2 p_0(X) (1-p_0(X))}{m(X)} \right] \\
&+E\left[ 1_{\mathcal{X}^{-}} \frac{s(0,X) q_0(1-1/p_0(X),X)^2 p_0(X)^{-1} (1-p_0(X)^{-1})}{1-m(X)} \right]\\
&+E\left[1_{\mathcal{X}^{+}} (q_1(p_0(X),X)-\beta_{1,1}(X,p_0(X)))^2\left(\frac{ s(0,X)(1-s(0,X))}{1-m(X)}+
\frac{p_0(X)^2 s(1,X)(1-s(1,X))}{m(X)} \right) \right] \\
&+E\bigg[1_{\mathcal{X}^{-}}(q_0(1-1/p_0(X),X)-\beta_{0,0}(X,1-1/p_0(X)))^2 \\
&\quad \times \left(\frac{ p_0(X)^{-2} s(0,X)(1-s(0,X))}{1-m(X)}+
\frac{s(1,X)(1-s(1,X))}{m(X)} \right) \bigg]\\
&- 2E\left[ 1_{\mathcal{X}^{+}}
\frac{q_1(p_0(X),X) \beta_{1,1}(X,p_0(X)) s(1, X) p_0(X)(1-p_0(X))}{m(X)} \right]\\
&-2E\left[ 1_{\mathcal{X}^{-}}
\frac{q_0(1-1/p_0(X),X) \beta_{0,0}(X,1-1/p_0(X)) s(0, X) p_0(X)^{-1}(1-p_0(X)^{-1})}{1-m(X)} \right]\\
&+2E\left[ 1_{\mathcal{X}^{+}}
\frac{(\beta_L(X,1,1)-\beta_L(1,1))(q_1(p_0(X),X)-\beta_{1,1}(X,p_0(X))) s(0, X) (1-s(0,X))}{1-m(X)} \right] \\
&-2E\left[ 1_{\mathcal{X}^{-}}
\frac{(\beta_L(X,1,1)-\beta_L(1,1))(q_0(1-1/p_0(X),X)-\beta_{0,0}(X,1-1/p_0(X))) s(1, X) (1-s(1,X))}{m(X)} \right],
\end{align*} }
where
{\footnotesize \begin{align*}
\sigma_{1}^2(x) &=
Var[Y1_{\{ Y \leq q_1(\min\{p_0(x),1\},x) \}}|S=1,D=1,X=x], \\
\sigma_{0}^2(x) &=
Var[Y 1_{\{ Y \geq q_0(\max\{1-1/p_0(x),0\},x) \}}|S=1,D=0,X=x].
\end{align*} }
\end{corr}
The setting with full selection, $S=1$ almost surely, renders any truncation redundant and thus $\beta_L(x,1,1) = \beta_U(x,1,1) = \beta(x,1,1)$ are equal to the conditional average treatment effect. Consequently, the efficiency bound in Corollary \ref{corr_var_bounds} becomes
\begin{equation}
E\left[ \frac{\sigma_1^2 (X) }{m(X)}+\frac{\sigma_0^2 (X) }{1-m(X)}+(\beta(X,1,1)-\beta(1,1))^2 \right],
\end{equation} with
\begin{align*}
\sigma_{d}^2(x) &=
Var[Y|S=1,D=d,X=x],
\end{align*}
which is exactly the efficiency bound for the ATE obtained in Theorem 2 of \cite{hahn1998role}, i.e.~Corollary \ref{corr_var_bounds} nests this as a special case.
The smooth bounds \eqref{def_smoothbeta} circumvent the irregularity problem by relaxing the width of the identified set using smooth approximators, i.e.~the $g$-functions. If the latter are sufficiently smooth, then the resulting outer bounds are pathwise differentiable. In particular, we obtain the following Theorem:
\begin{theorem}\label{thm_differentiable_smooth}
Suppose that Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_errors-of-g}, \ref{ass_density}, \ref{ass_strong_overlap} hold. If $\sup_{z,j}g_{i,h}'(z) \lesssim 1$ for some $h>0$, then $\beta_{L,h}(1,1)$ is pathwise differentiable.
\end{theorem}
The corresponding efficiency bound is in Appendix \ref{thm_differentiable_smooth}. The additional smoothness assumption applies to the LogSumExp example shown in Section \ref{sec_smooth}.
\section{Estimation and Inference} \label{sec_estimationinference}
\subsection{Estimators}
We propose to construct estimators by using empirical analogues to the efficient influence functions in Table \ref{tab_EIFs1} with estimated nuisances.
Let $E_n[X] = n^{-1}\sum_i^nX_i$. For given nuisance estimates $\hat{\eta}$, any estimator $\hat{\beta}_B = \hat{\beta}_B(s_0,s_1)$ solves \begin{align}
E_n[\psi_{\beta_B}(W,\hat{\eta},\hat{\beta}_B)] = 0
\end{align}
Note that all principal strata and margin influence functions have the following linear structure
\begin{align} \label{eq_IF_linear1}
\psi_{\beta_B}(W,\eta,\beta_B) = \psi^{[S]}_{\beta_B}(W,\eta) - \psi^{[B]}_{\beta_B}(W,\eta)\beta_B.
\end{align}
Thus, the influence function-based estimators have a ratio structure of the form \begin{align} \label{eq_esthat1}
\hat{\beta}_{B} &= \frac{E_n[\psi^{[B]}_{\beta_B}(W,\hat{\eta})]}{E_n[\psi^{[S]}_{\beta_B}(W,\hat{\eta})] }.
\end{align}
For the smooth bounds, parameters are given by a sum
\begin{align}
\beta_{B,h} = \beta_{B,+,h} + \beta_{B,-,h},
\end{align}
whose components also have linear influence functions with structure
\begin{align} \label{eq_IF_linear_smooth}
\psi_{\beta_{B,+,h}}(W,\eta,\beta_B) &= \psi^{[S]}_{\beta_{B,+,h}}(W,\eta) - \psi^{[B]}_{\beta_{B,+,h}}(W,\eta)\beta_{B,+,h}, \\
\psi_{\beta_{B,-,h}}(W,\eta,\beta_{B,-,h}) &= \psi^{[S]}_{\beta_{B,-,h}}(W,\eta) - \psi^{[B]}_{\beta_{B,-,h}}(W,\eta)\beta_{B,-,h},
\end{align}
which combined with estimated nuisances yield moment estimators
\begin{align}
\hat{\beta}_{B,h}
&= \hat{\beta}_{B,+,h} + \hat{\beta}_{B,-,h} \notag \\
&= \frac{E_n[\psi^{[B]}_{\beta_{B,+,h}}(W,\hat{\eta})]}{E_n[\psi^{[S]}_{\beta_{B,+,h}}(W,\hat{\eta})]} + \frac{E_n[\psi^{[B]}_{\beta_{B,-,h}}(W,\hat{\eta})]}{E_n[\psi^{[S]}_{\beta_{B,-,h}}(W,\hat{\eta})]}.
\end{align}
For the intensive margin, we also consider the moment functions suggested by \cite{semenova2023generalized} and \cite{heiler2024heterogeneous} under known propensity scores that have influence function \begin{align}
\tilde{\psi}_{\beta_{B}}(W,\eta,\beta_B) = \psi_{\beta_B}(W,\eta,\beta_B) - \delta(D,X,\eta),
\end{align}
with mean-zero term $\delta(D,X,\eta)$ characterized by {\begin{align}
E&[\min\{s(0,X),s(1,X)\}]\delta(D,X,\eta) \notag \\ &= \mathbbm{1}_{(X \in \mathcal{X}^+)}s(0,X)\bigg[\beta_{1,1}(X,p_0(X))\bigg(1-\frac{D}{m(X)}\bigg) \notag \\ &\quad -\beta_{0,0}(X,0)\bigg(1-\frac{1-D}{1-m(X)}\bigg)\bigg] \notag \\
&+ \mathbbm{1}_{(X \in \mathcal{X}^-)}s(1,X)\bigg[\beta_{1,1}(X,1)\bigg(1-\frac{D}{m(X)}\bigg) \notag \\ &\quad - \beta_{0,0}(X,1-1/p_0(X))\bigg(1-\frac{1-D}{1-m(X)}\bigg)\bigg].
\end{align}}
Solving for the parameter at the sample average then yields the alternative estimator \begin{align} \label{eq_est_btilde}
\tilde{\beta}_{B} &= \frac{E_n[\tilde{\psi}^{[B]}_{\beta_{B}}(W,\hat{\eta})]}{E_n[\psi^{[S]}_{\beta_{B}}(W,\hat{\eta})] }
\end{align}
\begin{table}[!h] \caption{Influence Functions: Regular Bound $\beta = \beta_L(1,1)$} \label{tab_EIFs1}
\centering
{\footnotesize \begin{tabular}{lc|c} \hline \hline && \\[-1.5ex]
Partition & Component & {Non-centered Influence Function} \\ \hline
& & \\
$\mathcal{X}^+$& $\psi_{\beta}^{[L^+]}$ & $
\frac{SD}{m(X)} Y 1\{Y \le q_1(p_0(X),X)\}
-\frac{S(1-D)}{1-m(X)}
Y $ \\ &&
$- \frac{SD}{m(X)}q_1(p_0(X),X)
\left[1\{Y \le q_1(p_0(X),X)\}-p_0(X) \right]
$ \\ &&
$+q_1(p_0(X),X)\left[\frac{1-D}{1-m(X)}(S-s(0,X))- p_0(X)
\frac{D}{m(X)}(S-s(1, X)) \right]$ \\ &&
$+ s(0,X)\left[\beta_{1,1}(X,p_0(X))\left(1-\frac{D}{m(X)}\right) - \beta_{0,0}(X,0)\left(1-\frac{1-D}{1-m(X)}\right)\right]$
\\ && \\
~ & $\psi_{\beta}^{[S^+]}$ & $ s(0,X) + \frac{(1-D)(S-s(0,X))}{1-m(X)}$ \\
&& \\
$\mathcal{X}^-$& $\psi_{\beta}^{[L^-]}$ & $
\frac{SD}{m(X)} Y
-\frac{S(1-D)}{1-m(X)}Y1\{Y \geq q_0(1-1/p_0(X),X)\}
$ \\ &&
$- \frac{S(1-D)}{1-m(X)}q_0(1-1/p_0(X),X)
\left[1/p_0(X) - 1\{Y \geq q_1(p_0(X),X)\} \right]
$ \\ &&
$-q_1(p_0(X),X)\left[\frac{D}{m(X)}(S-s(1,X))- p_0(X)^{-1}
\frac{1-D}{1-m(X)}(S-s(0,X)) \right]$ \\ &&
$+ s(1,X)\left[\beta_{1,1}(X,1)\left(1-\frac{D}{m(X)}\right) - \beta_{0,0}(X,1-1/p_0(X))\left(1-\frac{1-D}{1-m(X)}\right)\right]$
\\ && \\
~ & $\psi_{\beta}^{[S^-]}$ & $ s(1,X) + \frac{D(S-s(1,X))}{m(X)}$ \\ \hline \\[-0.5ex]
\multicolumn{3}{c}{$ \psi_{\beta} = \big[\mathbbm{1}{(X \in \mathcal{X}^+)}\psi_{\beta}^{L^+} + \mathbbm{1}{(X \in \mathcal{X}^-)}\psi_{\beta}^{L^-}\big] - \big[\mathbbm{1}{(X \in \mathcal{X}^+)}\psi_{\beta}^{S^+} + \mathbbm{1}{(X \in \mathcal{X}^-)}\psi_{\beta}^{S^-}\big] \beta$ } \\[1ex] \hline && \\[-1.5ex] \multicolumn{1}{l}{Nuisance} && Definition \\ \hline && \\[-0.5ex]
$m(x)$ && $P(D=1|X=x)$ \\
$s(d,x)$ && $P(S=1|D=d,X=x)$ \\
$q_d(u,x)$ && $\inf\{q \in \mathcal{Y}: u \leq P(Y\leq q|S=1,D=d,X=x)\}$ \\
$\beta_{1,d}(x,u)$ && $E[Y|S=1,D=d,X=x,Y\leq q_d(u,x)]$ \\
$\beta_{0,d}(x,u)$ && $E[Y|S=1,D=d,X=x,Y\geq q_d(u,x)]$ \\ && \\[-0.5ex]
\hline \hline
\end{tabular} }
\end{table}
\begin{table}[!h] \caption{Influence Functions: Smooth Bound $\beta_{L,h}(1,1) = \beta_{L,+,h}(1,1) + \beta_{L,-,h}(1,1)$} \label{tab_EIFs_smooth}
\centering
{\footnotesize \begin{tabular}{l|c} \hline \hline \\[-1.5ex]
Function & {Definition} \\ \hline & \\[-0.5ex]
$f^{[L]}(j,k)$ & $\frac{SD}{m(X)}\frac{g'_{k,h}(\beta_{L,h}(X))g_{j,h}(p_0(X))}{g_{1,h}(p_0(X))}\left[Y\mathbbm{1}(Y \leq q_1(g_{1,h}(p_0(X)),X) - g_{1,h}(p_0(X))\beta_{1,1}(x,g_{1,h}(p_0(x))\right] $ \\
& $-\frac{S(1-D)}{1-m(X)}\frac{g'_{k,h}(\beta_{L,h}(X))g_{j,h}(p_0(X))}{p_0(X)g_{1,h}(1/p_0(X))}$\\&$\times \left[Y\mathbbm{1}(Y \geq q_0(1-g_{1,h}(1/p_0(X)),X)) - g_{1,h}(1/p_0(X))\beta_{0,0}(x,1-g_{1,h}(p_0(x)))\right]$ \\
& $+ g_{k,h}(\beta_{L,h}(X))\left[g_{j,h}(p_0(X))s(1,X) + g'_{j,h}(p_0(X))\frac{1-D}{1-m(X)}(S-s(0,X))\right]$\\
& $+ g_{k,h}(\beta_{L,h}(X))\left[g_{j,h}(p_0(X)) - p_0(X) g'_{j,h}(p_0(X))\frac{D}{m(X)}(S-s(1,X))\right]$ \\
& $-\frac{SD}{m(X)}\frac{g'_{k,h}(\beta_{L,h}(X))g_{j,h}(p_0(X))}{g_{1,h}(p_0(X))}q_1(g_{1,h}(p_0(X)),X)\left[\mathbbm{1}(Y\leq q_1(g_{1,h}(p_0(X)),X)-g_{1,h}(p_0(X))\right]$ \\
& $+ \frac{S(1-D)}{1-m(X)}\frac{g'_{k,h}(\beta_{L,h}(X))g_{j,h}(p_0(X))}{p_0(X)g_{1,h}(1/p_0(X))}q_0(1-g_{1,h}(1/p_0(X)),X)$ \\ & $\times \left[\mathbbm{1}(Y \geq q_0(1-g_{1,h}(1/p_0(X)),X) - g_{1,h}(1/p_0(X))\right]$ \\
& $+ g'_{k,h}(\beta_{L,h}(X))\left[\frac{1-D}{1-m(X)}(S-s(0,X)) - p_0(X)\frac{D}{m(X)}(S-s(1,X)) \right]$\\
& $\times \bigg[\frac{g'_{1,h}(p_0(X))g_{j,h}(p_0(X))}{g_{1,h}(p_0(X))}(q_1(g_{1,h}(p_0(X)),X) - \beta_{1,1}(x,g_{1,h}(p_0(x))) $ \\
& $+ \frac{g'_{1,h}(1/p_0(X))g_{j,h}(p_0(X))}{p^2_0(X)g_{1,h}(1/p_0(X))}(q_0(1-g_{1,h}(1/p_0(X)),X) - \beta_{0,0}(x,1-g_{1,h}(p_0(x))))\bigg]$ \\
& \\
$f^{[S]}(j)$ & $\left[g_{j,h}(p_0(X))s(1,X) + g'_{j,h}(p_0(X))\frac{1-D}{1-m(X)}(S-s(0,X))\right]$ \\
& $\left[g_{j,h}(p_0(X)) + p_0(X)g'_{j,h}(p_0(X))\frac{D}{m(X)}(S-s(1,X))\right]$ \\
& \\
$\psi^{[L]}_{\beta_{L,+,h}}$ & $f^{[L]}(1,4)$ \\
$\psi^{[S]}_{\beta_{L,+,h}}$ & $f^{[S]}(3)$ \\
$\psi^{[L]}_{\beta_{L,-,h}}$ & $f^{[L]}(3,5)$ \\
$\psi^{[S]}_{\beta_{L,-,h}}$ & $f^{[S]}(1)$ \\
$\psi_{\beta_{L,+,h}}$ & $\psi^{[S]}_{\beta_{L,+,h}} - \psi^{[L]}_{\beta_{L,+,h}}\beta_{L,+,h}$\\
$\psi_{\beta_{L,-,h}}$ & $\psi^{[S]}_{\beta_{L,-,h}} - \psi^{[L]}_{\beta_{L,-,h}}\beta_{L,-,h}$ \\[1ex] \hline \\[-1.5ex] Nuisance &Definition \\ \hline & \\[-0.5ex]
$m(x)$ & $P(D=1|X=x)$ \\
$s(d,x)$ & $P(S=1|D=d,X=x)$ \\
$q_d(u,x)$ & $\inf\{q \in \mathcal{Y}: u \leq P(Y\leq q|S=1,D=d,X=x)\}$ \\
$\beta_{L,h}(x)$ & $\beta_{1,1}(x,g_{1,h}(p_0(x)) - \beta_{0,0}(x,1-g_{1,h}(p_0(x)))$ \\
$\beta_{1,d}(x,u)$ & $E[Y|S=1,D=d,X=x,Y\leq q_d(u,x)]$ \\
$\beta_{0,d}(x,u)$ & $E[Y|S=1,D=d,X=x,Y\geq q_d(u,x)]$ \\[1ex] \hline \\[-1.5ex] & \multicolumn{1}{c}{$g$-functions} \\ \hline & \\[-0.5ex]
&\multicolumn{1}{c}{$g_{1,h}(z) \leq \min\{z,1\} \leq g_{3,h}(z)$} \\
&\multicolumn{1}{c}{$g_{4,h}(z) \leq \max\{z,0\} \leq g_{2,h}(z)$} \\
&\multicolumn{1}{c}{$g_{5,h}(z) = -g_{2,h}(-z)$} \\
\hline \hline
\end{tabular} }
\end{table}
\subsection{Definitions and Assumptions}
In this section, we provide and discuss the assumptions and results for large sample inference for the regular and irregular case.
Denote $||\cdot||_p$ as the $L_p$ norm. Let $\eta \in T$ where $T$ is a convex subset of some normed vector space. We assume all nuisance functions are cross-fitted as in \cite{chernozhukov2018double}, Definition 3.1 and use the LogSumExp function for smoothing, see Section \ref{sec_smooth}. Denote nuisance realization set $\mathcal{T}_n = (\mathcal{S}_{0,n}\times\mathcal{S}_{1,n}\times \mathcal{Q}_{0,n} \times \mathcal{Q}_{1,n} \times \mathcal{M}_n \times \mathcal{B}_{0,0,n}\times \mathcal{B}_{0,1,n}\times \mathcal{B}_{1,0,n}\times \mathcal{B}_{1,0,n}) \subset T$ as the set that with high probability contains estimators $\hat{\eta} = \{\hat{s}(0,X), \hat{s}(1,X),\hat{q}_0(u,X),\hat{q}_1(u,X),\hat{m}(u,X),\hat{\beta}_{0,0}(u,X),$ $\hat{\beta}_{0,1}(u,X),\hat{\beta}_{1,0}(u,X),\hat{\beta}_{1,1}(u,X)\}$ for nuisance quantities $\eta = \{s(0,X),s(1,X),q_0(u,X), $ $q_1(u,X),m(X),{\beta}_{0,0}(u,X), {\beta}_{0,1}(u,X),$ ${\beta}_{1,0}(u,X),{\beta}_{1,1}(u,X)\}$. Let their corresponding $L_p$ error rates be
\begin{align*}
\lambda_{s,n,p} &= \sup_{d\in\{0,1\}}\sup_{\hat{s}(d) \in \mathcal{S}_{d,n}} E[|\hat{s}(d,X) - s(d,X)|^p]^{1/p}, \\
\lambda_{q,n,p} &= \sup_{d\in\{0,1\}}\sup_{u \in \tilde{U}}\sup_{\hat{q}_d(u) \in \mathcal{Q}_{d,n}} E[|\hat{q}_d(u,X) - q_d(u,X)|^p]^{1/p}, \\
\lambda_{m,n,p} &= \sup_{\hat{m} \in \mathcal{M}_{n}} E[|\hat{m}(X) - m(X)|^p]^{1/p}, \\
\lambda_{b,n,p} &= \sup_{j,d\in\{0,1\}}\sup_{u \in \tilde{U}}\sup_{\hat{\beta}_{j,d}(u) \in \mathcal{B}_{j,d,n}} E[|\hat{\beta}_{j,d}(u,X) - \beta_{j,d}(u,X)|^p]^{1/p},
\end{align*}
where $\tilde{U}$ is a compact subset of $(0,1)$ containing the relevant quantile trimming threshold support unions.\footnote{For the regular case, this will be $([\mbox{supp}(p_0(X)) \cup \mbox{supp}(1-p_0(X))] \cap \mathcal{X}^{+}) \cup ([\mbox{supp}(1/p_0(X)) \cup \mbox{supp}(1-1/p_0(X))] \cap \mathcal{X}^{-})$ while for the smoothed bounds, the subset will depend on $h$ and is defined equivalently with all bounds replaced by their $g_{1,h}$-smoothed counterparts.}
\begin{ass}[Outcome Regularity III]
The outcome has
a continuous conditional density $f(y|X=x,D=d,S=s)$ that is uniformly bounded from above and away from zero with bounded first derivative for any $x\in \mathcal{X}$ and $s,d\in \{0,1\}$. \label{ass_regularEst}
\end{ass}
Assumption \ref{ass_regularEst} here strengthens the smoothness and moment assumptions. This is required to bound (higher order) terms in some of the expansions and a central limit theorem for transformations of the outcome with bounded weights.
\subsection{Regular Bounds}
We first consider the regular case that sets $P(\mathcal{X}^0) = 0$. For this, we require the density around of $s(1,x) - s(0,x)$ in a neighborhood around zero to be well-behaved. This neighborhood is allowed shrink at a rate proportional to the error over the realization set for the selection probabilities. In particular,
\begin{ass}[Margin]
There exist a finite constant $C > 0$ and $\alpha > 0$ such that \begin{align*}
P(|s(1,x) - s(0,x)| \leq \delta) \leq C\delta^{\alpha} \quad \text{for any }\ 0 \leq \delta \leq c_n
\end{align*}
and $c_n \lesssim \lambda_{s,n,2}^{1/(2+\alpha)}$. \label{ass_margin1}
\end{ass}
Note that this implies that $P(\mathcal{X}^0) = 0$ which we know is necessary for the existence of a first-order uniformly unbiased estimator by Theorem \ref{thm_differentiable}. Any continuous bounded density for $s(1,x) - s(0,x)$ is sufficient yielding $\alpha = \infty$. However, a density may also diverge around zero at a controlled rate.
We further require the following learning rates.
\begin{ass}[Machine Learning Bias I] Let $e_n = o(1)$. For all folds, the nuisance parameters obtained via cross-fitting belong to a shrinking neighborhood $\mathcal{T}_n$ around $\eta$ with probability of at least $1-e_n$, such that
\begin{align*}
\lambda_{s,n,1} + \lambda_{s,n,2}^{\frac{\alpha}{2 + \alpha} } + \lambda_{q,n,1} + \lambda_{q,n,2} + \lambda_{m,n,2} + \lambda_{b,n,2} = o(1)
\end{align*}
and \begin{align*}
\sqrt{n}(\lambda_{s,n,2}^{\frac{2\alpha}{2 + \alpha}} + \lambda_{q,n,2}^2 + \lambda_{m,n,2}^2 + \lambda_{b,n,2}^2) = o(1).
\end{align*} \label{ass_MLbias1}
\end{ass}
One can see that we require $L_1$ and $L_2$ consistency for quantiles and conditional selection probabilities. The $L_1$ requirement is due to the variance of the trimming indicator for the outcome. Moreover, to control the machine learning bias, the squared $L_2$ rates have to be sufficiently fast. For most nuisances, these are equivalent to the usual $o(n^{-1/4})$ root mean squared error (RMSE) requirement in the debiased machine learning literature \citep{chernozhukov2018double}. The more demanding rate requirement for the selection probabilities is due to difficulties of classifying positive and negative monotonicity types around the boundary. Given these rates, Assumption \ref{ass_margin1} is only really plausible for $\alpha > 2$ as, even in the parametric case where $\lambda_{s,n,2} \sim n^{-1/2}$, $\sqrt{n}\lambda_{s,n,2}^{2\alpha/(2+\alpha)} \nrightarrow 0$ when $\alpha \leq 2$. For more demanding selection probabilities in high dimensions or with complicated functional forms, the density needs to be increasingly well-behaved around the margin. As $\alpha \rightarrow \infty$, the requirement reduces to the typical $o(n^{-1/4})$ rate as in \cite{heiler2024heterogeneous}.
We obtain the following Theorem.
\begin{theorem} \label{thm_asyN_regular}
Under Assumptions \ref{ass_CIA}, \ref{ass_weakMon1}, \ref{ass_STRONGoverlap}, \ref{ass_regularEst}, \ref{ass_margin1}, and \ref{ass_MLbias1}, the regular estimator is asymptotically normal and semiparametrically efficient, i.e.\begin{align*}
\sqrt{n}(\hat{\beta}_L(1,1) - \beta_L(1,1)) \overset{d}{\rightarrow} \mathcal{N}\big(0,Var_a(\beta_L(1,1))\big),
\end{align*}
where $Var_a(\beta_L(1,1))$ is the semiparametric efficiency bound in Theorem \ref{thm_differentiable}.
\end{theorem}
\subsection{Smooth Bounds}
We now consider the large-sample behavior of the smooth bounds. We impose the following learning requirements.
\begin{ass}[Machine Learning Bias II] Let $e_n = o(1)$. For all folds, the nuisance parameters obtained via cross-fitting belong to a shrinking neighborhood $\mathcal{T}_n$ around $\eta$ with probability of at least $1-e_n$, such that
\begin{align*}
\lambda_{s,n,1} + \lambda_{s,n,2} + \lambda_{q,n,1} + \lambda_{q,n,2} + \lambda_{m,n,2} + \lambda_{b,n,2} = o(1)
\end{align*}
and \begin{align*}
\sqrt{n}(\lambda_{s,n,2}^2 + \lambda_{q,n,2}^2 + \lambda_{m,n,2}^2 + \lambda_{b,n,2}^2 ) = o(1).
\end{align*} \label{ass_MLbias2}
\end{ass}
We obtain the following Theorem.
\begin{theorem} \label{thm_asyN_smooth}
Under Assumptions \ref{ass_CIA}, \ref{ass_weakMon1}, \ref{ass_STRONGoverlap}, \ref{ass_regularEst}, and \ref{ass_MLbias2}, and $h > 0$ such that $\sup_{z,j}g_{j,h}'(z) \lesssim 1$, the smooth bound estimator is asymptotically normal and semiparametrically efficient, i.e.\begin{align*}
\sqrt{n}(\hat{\beta}_{L,h}(1,1) - \beta_{L,h}(1,1)) \overset{d}{\rightarrow} \mathcal{N}\big(0,Var_a(\beta_{L,h}(1,1))\big),
\end{align*}
where $Var_a(\beta_{L,h})$ is the semiparametric efficiency bound for $\beta_{L,h}$.
\end{theorem}
We can see that, without restricting the distribution of $s(1,x) - s(0,x)$, i.e.~allowing for point mass $P(\mathcal{X}^0) > 0$ or arbitrary density around zero, the smooth bounds are asymptotically normal under weaker assumptions than the non-smooth bounds with restricted distributions. In particular, they only require $L_1$ and $L_2$ consistency as well as root mean squared error rates for all nuisance quantities of order $o(n^{-1/4})$ as in standard debiased machine learning \citep{chernozhukov2018double}.
This is fundamentally due to the smoothing turning the target parameter into a pathwise differentiable object.
Our smooth estimator differs from \cite{semenova2023generalized} in the following way: \cite{semenova2023generalized} uses a moment shifting or shrinkage procedure based on the selection probabilities. These shifting methods are known to be quite sensitive to the choice of tuning/selection parameter \citep{andrews2010inference}.
Moreover, even when $P(\mathcal{X}^0) = 0$, the estimator by \cite{semenova2023generalized} in finite samples should tend to narrow estimated identified sets as observations close to $p_0(x) = 1$ use a moment for the standard difference in conditional outcome means without truncation. Thus, while asymptotically valid due to correct classification in the limit, it will tend to be over-reject a correct null hypothesis in finite samples. The smoothed estimators, on the other hand, are constructed to be always wider than the estimators using no smoothing. They therefore estimate an outer identification region and the estimated identified will always be at least as wide as the naive generalized Lee bounds. This suggests better finite-sample size control at the expense of some power, see Appendix \ref{sec_montecarlo} for Monte Carlo simulations.
\subsection{The Efficiency Gap and Known Propensity Scores}
If propensity scores are known, the moment functions from \cite{semenova2023generalized} or \cite{heiler2024heterogeneous} without correction can be used as well. The corresponding estimator \eqref{eq_est_btilde} does not require estimation of $\beta_{j,d}(u,X)$ for any $j,d\in\{0,1\}$ and thus all $\lambda_{b,n,p}$ terms that enter Assumption \ref{ass_MLbias1} and \ref{ass_MLbias2} can be omitted. However, there is an information loss from ignoring the conditional variation in the truncated means. In particular, we obtain the following Theorem.
\begin{theorem} \label{thm_efficiency_gap}
Under Assumptions \ref{ass_CIA}, \ref{ass_weakMon1}, \ref{ass_STRONGoverlap}, \ref{ass_regularEst}, \ref{ass_margin1}, and \ref{ass_MLbias1} with $\lambda_{b,n,p} = 0$ known, $\tilde{\beta}_{L}(1,1)$ is asymptotically normal \begin{align*}
\sqrt{n}(\tilde{\beta}_L(1,1) - \beta_L(1,1)) \overset{d}{\rightarrow} \mathcal{N}\big(0,Var_{b}(\beta_L(1,1))\big),
\end{align*}
with efficiency gap
\begin{align*}
&E[\min\{s(0,X),s(1,X)\}]^2(Var_{b}(\beta_L(1,1)) - Var_a(\beta_L(1,1))) \\
&={E\bigg[\mathbbm{1}_{\{X \in \mathcal{X}^+\}}s(0,X)^2 \bigg( \beta_{1,1}(X,p_0(X))\sqrt{\frac{1-m(X)}{m(X)}} - \beta_{0,0}(X,0)\sqrt{\frac{m(X)}{1-m(X)}}\bigg)^2 \bigg]} \\
&\quad + {E\bigg[\mathbbm{1}_{\{X \in \mathcal{X}^-\}}s(1,X)^2 \bigg( \beta_{1,1}(X,1)\sqrt{\frac{1-m(X)}{m(X)}} - \beta_{0,0}(X,1-1/p_0(X))\sqrt{\frac{m(X)}{1-m(X)}}\bigg)^2 \bigg]}.
\end{align*}
\end{theorem}
Theorem \ref{thm_efficiency_gap} demonstrates that, under known propensity scores, there is a trade-off between efficiency and learning rates for the truncated conditional means $\beta_{j,d}(u,X)$. A lack of efficiency from not including information about the conditional truncated potential outcome means in the respective influence function arises as the variance of the truncated variables are bounded from below by their conditional analogues. This is equivalent to the comparison between the variance of a dependent variable and its residual variance in standard regression analysis. Theorem \ref{thm_efficiency_gap} also nests the case without sample selection. In particular, when $S=1$ almost surely, the ATE is point identified $\beta_L(1,1) = \beta_U(1,1)$ and all truncated means simplify to their untruncated counterparts $\beta_{d,1}(x,1) = \beta_{d,0}(x,0) = E[Y|D=d,X=x]$.
In this case all correction terms vanish as well and the alternative estimator collapses to \begin{align*}
\tilde{\beta}_{L}(1,1) = {E_n\bigg[\frac{YD}{m(X)}\bigg]} - {E_n\bigg[\frac{Y(1-D)}{1-m(X)}\bigg]}.
\end{align*}
This is the standard Horvitz-Thompson-type IPW estimator for the ATE using true propensities which is known to not reach the semiparametric efficiency bound \citep{hahn1998role}. More precisely, Theorem \ref{thm_efficiency_gap} yields efficiency gap \begin{align*}
Var_{b}(\beta_L(1,1)) - Var_a(\beta_L(1,1)) = E\bigg[ \bigg( \beta_{1,1}(X,1)\sqrt{\frac{1-m(X)}{m(X)}} - \beta_{0,0}(X,0)\sqrt{\frac{m(X)}{1-m(X)}}\bigg)^2 \bigg].
\end{align*}
This comparison to the point-identified case suggests that, in the regular regime, one could potentially obtain an efficient estimator for the always-taker ATE bounds without modelling the conditional truncated quantiles analogously to \cite{hirano2003efficient} using nonparametrically estimated treatment propensities with sufficiently rich basis functions that asymptotically encode all information about the truncated conditional outcome means. We leave the development of such an alternative estimation approach for future research.
\section{Extensions} \label{sec_extensions}
\subsection{Mean Dominance}
The bounds in Proposition \ref{prop_identification1} can further be refined using a mean dominance assumptions. Mean dominance to bound causal effects in sample selection models has been previously suggested by \cite{zhang2003estimation} for always-takers and by \cite{huber2015sharp} for compliers under strong monotonicity without covariates.
\begin{ass}[Mean Dominance] \label{ass_dominance}
The conditional potential outcomes at the intensive margin are larger than at the extensive margin, i.e. \begin{align}
E[Y(d)|S(0)=1,S(1)=1,X=x] \geq E[Y(d)|S(0)\neq S(1), X=x]
\end{align}
for any $d=0,1$.
\end{ass}
This yields the following modified sharp bounds\footnote{
If only either complier or defier bounds are of interest, Assumption \ref{ass_dominance} can be relaxed to hold for the respective group only.}.
\begin{prop}[Identification with Mean Dominance] \label{prop_identification_dominance1}
Let $\mathcal{Y}_d(x)$ be the support of $Y(d)$ conditional on $X=x$ and $\overline{y}_d(x)$ and $\underline{y}_d(x)$ its respective infimum and supremum for $d\in\{0,1\}$. Under Assumptions \ref{ass_CIA}, \ref{ass_monotone}, \ref{ass_continuity}, \ref{ass_overlap}, and \ref{ass_dominance}, the lower bounds for the principal strata are given by \begin{align*}
\beta_L(x,s_0,s_1) = \begin{cases}
\beta_{1,1}(x,1) - \beta_{0,0}(x,1-\min\{1/p_0(x),1\}) &\text{ if } s_0 = 1,~s_1 = 1 \\
\beta_{1,1}(x,1-\min\{p_0(x),1\}) - \beta_{0,0}(x,0) &\text{ if } s_0 \neq s_1 \\
\underline{{y}}_1(x) - \overline{{y}}_0(x) &\text{ if } s_0 = 0,~s_1 = 0,
\end{cases}
\end{align*}
and the upper bounds by \begin{align*}
\beta_U(x,s_0,s_1) = \begin{cases}
\beta_{0,1}(x,1-\min\{p_0(x),1\}) - \beta_{1,0}(x,1) &\text{ if } s_0 = 1,~s_1 = 1 \\
\beta_{0,1}(x,0) - \beta_{1,0}(x,1-\min\{1/p_0(x),1\}) &\text{ if } s_0 \neq s_1 \\
\overline{{y}}_1(x) - \underline{{y}}_0(x) &\text{ if } s_0 = 0,~s_1 = 0.
\end{cases}
\end{align*}
\end{prop}
As the share of the principal strata are already identified via conditional monotonicity, aggregation to unconditional bounds uses the same conditional probability weights as in the case without mean dominance in \eqref{def_betaB_general} and thus the same moment function for the denominator or density part in \eqref{eq_esthat1}. For the numerator, all truncated quantities can be estimated as in the pure monotonicity case with the components in the influence function corresponding to the now untruncated means replaced by their standard untruncated augmented IPW counterparts for the conditional outcome means \citep{robins1994estimation}.
\subsection{Heterogeneous Bounds}
In many evaluation problems, objects of interest are heterogeneous effects, e.g.~the effect at a particular margin for an observable subgroup. \cite{heiler2024heterogeneous} argues that, in the nonparametric sample selection model and more broadly, there are two types of heterogeneity, 1) heterogeneous effects and 2) heterogeneity in the severity of the identification problem, i.e.~the width of the identified set. Exploiting these in combination can yield more precise inference and significant effect bounds even when unconditional aggregate bounds do not reject a null-effect. Our influence functions for any principal strata or margin can readily be used for the estimation of such heterogeneous effect bounds following the method of \cite{heiler2024heterogeneous}. In particular, let $f:\mathcal{X}\rightarrow \mathcal{Z}$ and define $Z = f(X)$ a low-dimensional subgroup or mapping from the covariate space. Recall that all presented influence functions have the following linear structure
\begin{align}
\psi_{\beta_B}(W,\eta,\beta_B) = \psi_{\beta_B}^{[B]}(W,\eta) - \psi^{[S]}_{\beta_B}(W,\eta)\beta_B,
\end{align}
where \begin{align}
E[\psi^{[B]}_{\beta_B}(W,\eta)|X=x] &= \beta_L(x,s_0,s_1)P(S(0)=s_0,S(1)=s_1|X=x), \notag \\
E[\psi^{[S]}_{\beta_B}(W,\eta)|X=x] &= P(S(0)=s_0,S(1)=s_1|X=x).
\end{align}
Thus, for any $z \in \mathcal{Z}$, \begin{align}
E[&Y(1)-Y(0)|S(0)=s_0,S(1)=s_1,Z=z] \notag\\
&\geq E[\beta_L(X,s_0,s_1)|S(0)=s_0,S(1)=s_1,Z=z] \notag\\
&= \frac{E[\beta_L(X,s_0,s_1)P(S(0)=s_0,S(1)=s_1|X)|Z=z]}{E[P(S(0)=s_0,S(1)=s_1|X)|Z=z]} \notag\\
&= \frac{E[\psi^{[B]}_{\beta_B}(W,\eta)|Z=z]}{E[\psi^{[S]}_{\beta_B}(W,\eta)|Z=z]},
\end{align}
by the law of iterated expectations. Hence we can obtain conditional effect bounds by aggregating the components of the presented influence functions within covariate partitions or by (nonparametric) projection onto (basis transformations of) $Z$. Regularity and/or the necessity for smoothing can be obtained under analogous conditions to the ones presented in this paper. For more details consider \cite{heiler2024heterogeneous}.
\section{Empirical Study: Labor Market Policy} \label{sec_empirical1}
\subsection{Job Corps and Data}
Job Corps (JC) is a US Department of Labor program aimed at empowering young individuals of economically or otherwise disadvantaged background with important labor market skills and qualifications. The intense program combines general education, vocational training, extracurricular activities, placement services, and more. Most participants are in residential slots in centers of varying size all over the US.
The National Job Corps Study, executed by Mathematica Policy Research in the mid-90s, analyzed the effects of the Job Corps program on labor market prospects. The experiment implemented stratified randomized assignment of applicants, incorporating over 15400 units between ages 16 to 24. Data on various outcomes such as income, job status, educational achievements, and criminal activities were gathered at different intervals.
The outcomes of the original Job Corps study were mixed. Initially, the program boosted the participants' educational achievements and income. However, most aggregate effects waned over time. There is evidence for increased earnings in the older participant population \citep{schochet2008does}. There is also well-document heterogeneity in terms of vocational and academic training returns as well a gender, see e.g.~\cite{flores2012estimating} or \cite{heiler2023effect}.
We employ the public use files which contain all survey data up to 208 weeks after the initial assignment \citep{burghardt1999national,schochet2008does}. Our main outcome measure is log hourly earnings as in \cite{lee2009training}. If workers are priced around their marginal productivity, this measure can be used as a proxy for increases in the latter caused by assignment to Job Corps.
We extract an extended set of covariates containing detailed information regarding demographics, employment, criminal history, education, health, expectations, regional characteristics, and other JC related information.\footnote{We collect all variables used in \cite{lee2009training} and additional information on individual background and behavior similar to \cite{flores2012estimating} and \cite{heiler2023effect}, as well as all variables that were used for stratification in the original experiment. Some of the latter have previously been overlooked but are necessary for correct specification of the propensity score, see Appendix \ref{app_JC_data1} for more details.}
All impact analysis is conditional on having a non-missing earnings entry, including zero, at weeks $t\in\{1,\dots,208\}$ identical to \cite{lee2009training} leaving a sample of 9.415 units.\footnote{\cite{lee2009training} suggests that non-response at this stage is a second order issue and initial treatment-control balance is preserved, see his Remark 2. We reassess this claim with the extended data using both formal tests for average differences as well as standardized differences and find little evidence for imbalance, see Appendix \ref{app_JC_data1}. The only exception is significantly larger worries to attend Job Corps in the treatment group.
Our prior is that worries are more likely negatively related to any treatment effect. Thus, we expect any potential bias to be negative, i.e.~earnings and employment effects are likely to be at least as large.}
In the following, results for nuisance functions are presented for the research sample. For any of the impact analysis as well as general descriptive statistics, we use the nationally representative design weights as in \cite{schochet2008does}.
For the impact evaluation, we analyze the effect per eligible participant, i.e.~the causal effect of assignment to Job Corps. We focus on the average treatment effect at the intensive margin as studied by \citep{lee2009training,semenova2023generalized}.
As participation or outside alternatives after assignment are not controlled, all Job Corps effects here are relative to the control state of not being assigned, including partaking in any alternative program and/or working.
Figure \ref{fig_JC_desc1} contains the weekly average earnings and employment rate for the 208 weeks after treatment assignment. There are systematic trends and differences between treatment and control groups in both earnings and employment over time.
There is a clear upwards trend in earnings and employment over time. The control group tends to be ahead early after assignment but catches up later in time exceeding the trajectory of the control group for both measures.
However, any differences between treatment and control earnings effectively consists of a mixture of units (the working) that could be made up of varying proportions of always-takers, compliers, and defiers, effectively leading to selected comparisons. Raw differences in employment rates are also uninformative about heterogeneity and composition of the intensive and extensive margins. We take this into account when constructing the bounds in what follows.
\begin{figure}
\centering \caption{Job Corps: Employment and Earnings} \label{fig_JC_desc1}
\begin{subfigure}{0.32\textwidth}
\includegraphics[width = \textwidth, trim = 0 80 0 50, clip]{figures/JC/plot_JC_weekly_Y}
\caption{Weekly Earnings in USD}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
\includegraphics[width = \textwidth, trim = 0 80 0 50, clip]{figures/JC/plot_JC_weekly_H}
\caption{Weekly Hours Worked}
\end{subfigure}
\begin{subfigure}{0.32\textwidth}
\includegraphics[width = \textwidth, trim = 0 80 0 50, clip]{figures/JC/plot_JC_weekly_S}
\caption{Weekly Employment Rate}
\end{subfigure} \footnotesize
\begin{justify}
Weakly earnings, hours worked, and employment averages of control and treated units. The $n=9415$ sample is restricted to units that report hours and earnings for all periods $t\in\{1,\dots,208\}$. All calculations use nationally representative design weights.
\end{justify}
\end{figure}
\subsection{Models and Estimation}
We give a brief overview over the estimation methods used. See Appendix \ref{app_JC_estimation} for more details. Treatment propensities are known by the design of the stratified experiment, see Appendix \ref{app_JC_data1}. All other nuisance parameters are estimated using 5-fold cross-fitting with all available observations and pre-treatment predictors at each time period reported. No cross-time information was used.
For the conditional selection probabilities we compare five different main specifications: (a) Logistic regression, (b) logistic regression with interacted treatment, (c) gradient boosted trees (XGBoost), (d) artificial neural network, and (e) random forest. These methods differ with respect to their capabilities to (i) allow for weak monotonicity, (ii) produce sparsity, (iii) be robust to extreme observations, and (iv) fit and predict employment probabilities and status well. Table \ref{tab:models_selection} contains an aggregate comparison between the different methods with respect to these properties. By construction, all models except for the simple logit allow for weak monotonicity. XGBoost, neural network, and random forest can produce sparsity.\footnote{Both the interacted logit and neural network tend to produce many or some relatively extreme estimates for the relative selection probabilities, see also Figure \ref{app_JC_p0x_all} in Appendix \ref{app_JC_results}. It seems questionable whether e.g.~extreme relative selection probabilities that suggest e.g.~a quadrupling in relative employment chances due to Job Corps are credible or mostly a by-product of an overly high variance model.} XGBoost is the model that has, by far, the best out-of-sample accuracy and loss as measured by the negative likelihood function outperforming \textit{all} other methods throughout 202/206 out of 208 periods, see Appendix \ref{app_JC_estimation} for the detailed results. This is in line with many studies that demonstrate boosted trees as or among the best performing methods for tabular data \citep{shwartz2022tabular}.
\begin{table}[!h] \centering\caption{Selection Model Estimates: Comparison} \label{tab:models_selection}
\footnotesize
\begin{tabular}{lccccc} \hline \hline \\[-1.5ex]
Method & Weak Monotonicity & Sparsity & Extreme Predictions & Accuracy & Likelihood \\[-0.5ex]
&&&&(rank)&(rank) \\[0.5ex]
Logit & \ding{55} & \ding{55} & \ding{55} & 3 & 4 \\
Logit Interacted & \ding{51} & \ding{55} & \ding{51} & 5 & 5 \\
XGBoost & \ding{51} & \ding{51} & \ding{55} & 1 & 1 \\
Neural Network & \ding{51} & \ding{51} & \ding{51} & 2 & 3 \\
Random Forest & \ding{51} & \ding{51} & \ding{55} & 4 & 2 \\ \hline
\end{tabular}
\begin{justify}
This table contains the different estimation methods, their properties as well as an overview over the resulting $p_0(x)$ estimates. Weak Monotonicity refers to Assumption \ref{ass_monotone}. Logit has that either $P(\mathcal{X}^+)$ or $P(\mathcal{X}^-)$ equal zero. Sparsity is achieved whenever $P(\mathcal{X}^0)$ can be larger than zero. Extreme Predictions refers to an average of at least 10 selection estimates outside of the $(0.25,4)$ interval. Accuracy and Likelihood contain the rank in the respective performance metric based on the cross-fitting error averaged over all periods $t\in\{1,\dots,208\}$.
\end{justify}
\end{table}
For conditional means, conditional quantiles, and truncated means required for the bound estimators, we use off-the-shelf honest regression and quantile forests across all specifications \citep{athey2019generalized}.
This was chosen to isolate the contribution of differences in relative selection probabilities which is the main issue addressed in this paper.\footnote{Quantile forests are also computationally attractive in high dimensions and have been used for sample selection bounds in previous research \citep{heiler2024heterogeneous,andersen2023guide}.}
For the impact analysis, we consider the bounds that (i) remove observations with $p_0(x) = 1$ and use the corresponding influence function (``Trim''), (ii) the moment shifting or selection method by \cite{semenova2023generalized} (``Shift''), (iii) smooth bounds (``Smooth''), and (iv) conventional Lee bounds using strong monotonicity. For the shifting method, we use the tuning parameter as suggested by \cite{semenova2023generalized}, $\rho_n = n^{-1/4}/\log n$. For both trimming and shifting, we report results with and without using the propensity score correction that affects efficiency. For the smooth bounds, we report the analysis on a grid of smoothing parameters.\footnote{On any fixed grid, inference is valid for any $h$, so one could pick the $h$ that provides the smallest confidence intervals to optimize power.} For inference on the effect we use \cite{imbens2004confidence} confidence intervals for all methods.
\subsection{Relative Selection Probabilities $p_0(x)$}
In this section, we consider the predicted relative selection probabilities close or equal to one for the methods that can produce sparsity. Figure \ref{fig:share_ones_001} depicts the share of probabilities for XGBoost, neural network, and random forest. The average amount of close-to-one units averaged over all time periods are $51.91\%$ (XGBoost), $18.75\%$ (Neural Net), and $42.36\%$ (Random Forest). Thus, there is a significant amount of units that can affect the various bounds differently.
In particular, for the XGBoost, almost all these probabilities are exact numerical ones, i.e.~the underlying model is sparse in some parts of the covariate space with respect to the treatment. Figure \ref{fig:share_ones_xgb} depicts the close-to-one and exact-one unit shares over all time periods. One can see that the trajectory of exact-one closely follows the close-to-one shares with an average rate of 44.36\%.
\begin{figure}
\centering \caption{Estimated $p_0(x)$: Sparsity and Margin}
\begin{subfigure}{0.48\textwidth}
\includegraphics[width = \textwidth, trim = 0 110 0 100, clip]{figures/JC/v3_plot_p0x_001}
\caption{$p_0(x)$ Close to One} \label{fig:share_ones_001}
\end{subfigure}
\begin{subfigure}{0.48\textwidth}
\includegraphics[width = \textwidth, trim = 0 110 0 100, clip]{figures/JC/v3_plot_p0x_xgb_both}
\caption{XGBoost: $p_0(x)$ Close or Equal to One}\label{fig:share_ones_xgb}
\end{subfigure} \footnotesize
\begin{justify} (a) contains the share of units for which each method returns relative selection probabilities $p_0(x)$ estimated within a 0.01 neighborhood of one at each post-assignment period $t=1,\dots,208$. (b) contains the results for XGBoost with relative selection probabilities $p_0(x)$ estimated exactly equal or within the 0.01 neighborhood of one respectively.
\end{justify}
\end{figure}
\subsection{Employment Effects}
Figure \ref{fig_JC_impact_employment1} contains the estimated average assignment effect on employment using the efficient estimator for all probability models. Estimates are almost identical. We replicate the finding that, in the short run, assignment reduces the likelihood of employment which only bounces back after about 1.5 years (the average duration of JC participation was around 8 months at the time). After about 2.25 years, positive employment effects tend to stabilize. In particular, most estimates are significantly negative for weeks 2 to 76 and significantly positive for weeks 116 to 208 ($p < 0.05$).
\begin{figure}[!h]
\centering \caption{Job Corps: Average Employment Effect Estimates} \label{fig_JC_impact_employment1}
\includegraphics[width = 0.6\textwidth, trim = 0 120 0 100, clip]{figures/JC/v3_plot_JC_impact_S_all} \footnotesize
\begin{justify}
Estimated average treatment effect of assignment to Job Corps on employment $E[S(1) - S(0)]$ using the efficient influence function with different selection probability models for $t\in\{1,\dots,208\}$. All calculations use nationally representative design weights.
\end{justify}
\end{figure}
\subsection{Causal Effects Bounds}
We now present the impact results for assignment to JC on log hourly earnings. We first present the results for the main evaluation in $t=208$ using the new methods and contrast them with results based on the literature and replications with comparable sample and design weights.
We then present the results using the best performing nuisance specification for methods that allow for conditional monotonicity for periods $t\in\{45,90,135,180,208\}$ as in \cite{lee2009training}. Robustness checks and results for additional methods, time periods, and smoothing parameters can be found in Appendix \ref{app:JC_effect_bounds_all1}.
\subsubsection{Main Effects at $t=208$}
Table \ref{tab:JC_results_208_overview1} contains the impact analysis for various methods and specifications at period $t=208$. All entries except for Lee (2009) are based on own calculations and with identical outcome and design weights. \cite{schochet2008does} is a simple treatment control difference for the sample of reported hourly log earnings. It is of similar magnitude when compared with the equivalent mean difference in the restricted Lee (2009) sample in line with the balancing analysis.
Switch, regular, and smooth bounds all admit flexible estimation of nuisance parameters. Using the best performing XGBoost model for selection, there are 4252 observations with exact $\hat{p}_0(x) = 1$. Thus, trimming bounds effectively have to use a much reduced sample. Nevertheless, estimates are in a similar range but with larger standard errors. The moment switching approach by \cite{semenova2023generalized} with known propensity scores provides overly wide bounds. This seems to be due to a strong sensitivity when there are many estimates for $p_0(x)$ around one. With unknown scores, estimates are tighter, however with standard errors much larger compared to alternatives.
\afterpage{
\begin{landscape}
\begin{table}[!h]
\caption{Impact Results at $t=208$: Overview}
\label{tab:JC_results_208_overview1}
\begin{threeparttable}
\centering \scriptsize
\begin{tabular}{llccccccc} \hline \hline
& Method & Impact Estimate & 95\% CI & Sample & Selection & Selection & Covariates & Sample \\
& & & & Selection & Model & Key Assumptions & & \\ \hline \\[-0.5ex]
Schochet et al.~(2008) & Mean Difference & 0.059 & $(0.031,~0.086)$ & \ding{55} & -- & Independence$^1$ & 0 & 10602 \\
& & & & & & & & \\
& & & & & & & & \\
Lee (2009) & Mean Difference & 0.058 & $(0.027,~0.088)$ & \ding{55} & -- & Independence$^1$ & 0 & 9415 \\
& & & & & & & & \\
& Bounds & $[-0.019,~0.093]$ & $(-0.049,~0.114)$ &\ding{51} & -- & Monotonicity & 0 & 9415 \\
& + Covariates & $[-0.012,~0.089]$ & $(-0.037,~0.112)$ &\ding{51} & Nonparametric & Monotonicity & $28/1/5^2$ & 9415 \\
& & & & & & & & \\
& Heckman (1979) & 0.015 & $(-0.008,~0.038)$ &\ding{51} & Probit & Exclusion$^3$, Normality & 28 & 9415 \\
& Das et al.~(2003) & 0.014 & $(-0.010,~0.038)$ &\ding{51} & Linear$^4$ & Exclusion$^3$, Single Index & $28^4$ & 9415 \\
& & & & & & & & \\
& & & & & & & & \\
Switch$^{5}$ & Known Propensity & $[-0.646,~1.514]$ & $(-0.438,~1.650)$ &\ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& Unknown Propensity & $[0.138,~0.162]$ & $(-0.003,~0.303)$ & \ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& & & & & & & & \\
& & & & & & & & \\
Trim$^{5,6}$ & Known Propensity & $[-0.008,~0.054]$ & $(-0.063,~0.111)$ &\ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& Unknown Propensity & $[-0.036,~0.077]$ & $(-0.201,~0.245)$ & \ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& & & & & & & & \\
Smooth Bounds$^{5,7}$ & $h = 5$ & $[-0.057,~0.275]$ & $(-0.095,~0.312)$ &\ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& $h = 1$ & $[-0.005,~0.113]$ & $(-0.040,~0.148)$ & \ding{51} & Nonparametric & Weak Mononicity & 77 & 9415 \\
& & & & & & & & \\ \hline
\end{tabular}
Impact estimates using different methods and samples with key assumptions for validity of the impact estimate. All estimates and intervals except for Lee (2009) are based on own calculations using design weights. Confidence intervals for point-identified methods use standard critical values. Confidence intervals for bounds use the \cite{imbens2004confidence} refinement. All methods use known propensity scores of the stratified experiment if necessary.
$^1$Independence of potential selection from potential outcomes.
$^2$28 covariates are used to generate a single predictive score. Bounds are calculated within 5 strata of this score.
$^3$Months employed is used as excluded variable in the selection model.
$^4$Semiparametric model also includes transformations of the original variables, see Lee (2009).
$^5$Best performing nuisance parameter specifications, see Appendix \ref{app_JC_estimation} for details.
$^6$Trim are the bounds that remove all units for which $\hat{p}_0(x) = 1$ with unknown/efficient and known/inefficient propensity score and related moment function.
$^7$Smoothing parameters $h$ scaled by 100.
\end{threeparttable}
\end{table}
\end{landscape}
}
Our preferred specification is smooth bounds with $h=1$. The impact estimate of $[-0.005,0.113]$ with 95\% confidence interval of $(-0.040,0.148)$ rules out large negative earnings effects and is surprisingly close and of similar width compared to results obtained under strong monotonicity as in \cite{lee2009training}. As strong monotonicity is rejected by the data \citep{semenova2023generalized}, this provides more evidence on the robustness of this final impact evaluation.
We also include a larger smoothing specification for comparison; for more values see Appendix \ref{app:JC_effect_bounds_all1}. As expected, the identified set becomes larger. In this case, the standard errors are still comparable such that overall less smoothing seems preferable.
\subsubsection{Secondary Effects and Robustness}
We now analyze the effect for multiple periods using all methods. We present the results for the two best nuisance models (XGBoost and neural network). For the switching and the trimming bounds, we only consider the unknown propensity score estimates, as they are more reliable across specifications. For smooth bounds we present two smoothing parameter results as above. Results for the other methods and parameters are robust with the exception of the interacted logistic model, that failed to converge in some folds.\footnote{See Appendix \ref{app_JC_results} for some histograms and descriptive statistics.}
\begin{table}[!h]
\begin{threeparttable}
\centering
\caption{Job Corps Impact Results, by week}
\label{tab:JC_results_both1}
\scriptsize
\begin{tabular}{llccccc}
\hline
\hline ~ \\[-0.5ex]
\multicolumn{7}{c}{Simple Monotonicity Bounds} \\[1ex]
~ & ~ & $t = 45$ & $t = 90$ & $t = 135$ & $t = 180$ & $t = 208$\\[-0.5ex] ~ \\
Lee (2009) & ~ & $[-0.074,~0.125]$ & $[0.042,~0.043]$ & $[-0.016,~0.076]$ & $[-0.033,~0.087]$ & $[-0.012,~0.089]$ \\
~ & ~ & $(-0.097,~0.152)$ & $(-0.005,~0.092)$ & $(-0.051,~0.099)$ & $(-0.064,~0.108)$ & $(-0.037,~0.112)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
~ \\[-0.5ex]
\multicolumn{7}{c}{Results for XGBoost Selection Model} \\[1ex]
~ & ~ & $t = 45$ & $t = 90$ & $t = 135$ & $t = 180$ & $t = 208$\\[-0.5ex] ~ \\
Smooth & $h = 5$ & $[-0.109,~0.145]$ & $[-0.112,~0.212]$ & $[-0.093,~0.213]$ & $[-0.072,~0.215]$ & $[-0.057,~0.275]$ \\
~ & ~ & $(-0.172,~0.199)$ & $(-0.161,~0.256)$ & $(-0.139,~0.250)$ & $(-0.107,~0.251)$ & $(-0.095,~0.312)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
~ & $h=1$ & $[-0.000,~0.161]$ & $[0.008,~0.096]$ & $[-0.028,~0.090]$ & $[-0.029,~0.093]$ & $[-0.005,~0.113]$ \\
~ & ~ & $(-0.068,~0.209)$ & $(-0.036,~0.135)$ & $(-0.066,~0.125)$ & $(-0.064,~0.128)$ & $(-0.040,~0.148)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
Trim
~ & ~ & $[-0.257,~0.126]$ & $[-0.175,~0.004]$ & $[-0.021,~0.152]$ & $[-0.180,-0.011]$ & $[-0.036,~0.077]$ \\
~ & ~ & $(-0.403,~0.269)$ & $(-0.745,~0.573)$ & $(-0.197,~0.327)$ & $(-0.328,0.134)$ & $(-0.201,~0.245)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
Shift
~ & ~ & $[-0.324,-0.076]$ & $[-0.021,-0.011]$ & $[0.022,~0.067]$ & $[0.075,~0.122]$ & $[0.138,~0.162]$ \\
~ & ~ & $(-0.470,~0.067)$ & $(-0.194,~0.162)$ & $(-0.116,~0.204)$ & $(-0.055,~0.253)$ & $(-0.003,~0.303)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
$|\hat{p}_0(x)-1|$ & $\leq \frac{n^{-1/4}}{\log(n)}$ & 2592 & 8680 & 6350 & 4871 & 5688 \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
\hline
\end{tabular}
\begin{tabular}{llccccc}
~ \\[-0.5ex]
\multicolumn{7}{c}{Results for Neural Network Selection Model} \\[1ex]
~ & ~ & $t = 45$ & $t = 90$ & $t = 135$ & $t = 180$ & $t = 208$\\[-0.5ex] ~ \\
Smooth & $h = 5$ & $[-0.177,~0.089]$ & $[-0.137,~0.214]$ & $[-0.089,~0.187]$ & $[-0.092,~0.182]$ & $[-0.058,~0.228]$ \\
~ & ~ & $(-0.234,~0.137)$ & $(-0.183,~0.255)$ & $(-0.132,~0.220)$ & $(-0.125,~0.215)$ & $(-0.094,~0.263)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
~ & $h=1$ & $[-0.469,~0.777]$ & $[-0.005,~0.105]$ & $[-0.023,~0.107]$ & $[-0.045,~0.112]$ & $[-0.022,~0.127]$ \\
~ & ~ & $(-0.543,~0.855)$ & $(-0.049,~0.147)$ & $(-0.061,~0.144)$ & $(-0.080,~0.149)$ & $(-0.058,~0.165)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
Trim & ~ & $[-0.347,~0.192]$ & $[-0.031,~0.003]$ & $[-0.107,~0.046]$ & $[-0.141,~0.038]$ & $[-0.157,~0.072]$ \\
~ & ~ & $(-0.481,~0.319)$ & $(-0.179,~0.149)$ & $(-0.223,~0.164)$ & $(-0.266,~0.159)$ & $(-0.284,~0.203)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
Shift & ~ & $[-0.424,-0.122]$ & $[-0.054,~0.011]$ & $[0.049,~0.108]$ & $[-0.062,~0.076]$ & $[0.046,~0.121]$ \\
~ & ~ & $(-0.571,~0.019)$ & $(-0.205,~0.162)$ & $(-0.078,~0.235)$ & $(-0.187,~0.199)$ & $(-0.084,~0.252)$ \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
$|\hat{p}_0(x)-1|$ & $\leq \frac{n^{-1/4}}{\log(n)}$ & 737 & 4842 & 3220 & 2180 & 3669 \\
~ & ~ & ~ & ~ & ~ & ~ & ~ \\
\hline
\end{tabular}
\footnotesize
This table contains the estimated always-taker bounds of assignment to JC on log hourly earnings using smoothed, trimmed, shifted, and regular (generalized) Lee bounds using $n = 9415$ for all weeks. Identified sets are in brackets and 95\% \cite{imbens2004confidence} confidence intervals in parenthesis below. The threshold $\frac{n^{-1/4}}{\log(n)}$ is the switching threshold suggested by \cite{semenova2023generalized} to choose between moments. Smoothing parameters are scaled by 100. All calculations use design weights.
\end{threeparttable}
\end{table}
Table \ref{tab:JC_results_both1} contains the estimates and confidence intervals for XGBoost and neural network selection scores. There are a few general observations. The choice of the selection model matters. In particular, the larger dispersion of predicted $p_0(x)$ for the neural net leads to wide intervals at $t=45$ for most methods. Trimming and switching seem more sensitive in general, see also Appendix \ref{app:JC_effect_bounds_all1}. For all other periods, results are relatively consistent across selection models. Trimming and switching bounds, however, tend to have relatively large standard errors compared to the smooth bounds. Results for additional methods reveal an identification-precision trade-off in some, but not all cases, see Appendix \ref{app:JC_effect_bounds_all1}.
The intervals for the smooth bounds at $h=1$ are, again, surprisingly close to the original \cite{lee2009training} specification. Even a strictly positive identified set at $t=90$ is recovered with XGBoost, albeit with larger width which seems more credible given the estimates of the surrounding periods. The Lee unconditional monotonicity can be rejected from the data \citep{semenova2023generalized}. However, it seems that the general finding for the later periods, with log hourly earnings effects in the range of around $-0.03$ to $0.11$, is robust to conditional monotonicity without parametric assumptions.
\section{Conclusion} \label{sec_conclusion}
This paper demonstrates the importance of heterogeneity in selection behavior with respect to estimation and inference on causal effects at the intensive and extensive margins in selected samples.
Allowing for different types of monotonicity crucially affects the (ir)regularity of sharp effect bounds and thus the ability to provide precise inference for the associated causal effects. This paper provides a solution in the form of
outer identification regions that can be estimated efficiently from a semiparametric perspective. The approach can handle modern machine learning methods that are able to better approximate such important heterogeneity in sample selection. We discover an empirically relevant trade-off between identification strength versus precision that is likely to apply to a much larger set of models and parameters beyond the ones considered in this paper, as in \cite{lee2021bounding} and \cite{pakel2023bounds}.
\addcontentsline{toc}{section}{References}
{\setstretch{1}
\bibliography{HKV}
}