Estimating Functionals of the Joint Distribution of Potential Outcomes with Optimal Transport

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.

92,316 characters

Estimating Functionals of the Joint Distribution of Potential Outcomes with Optimal Transport


	\maketitle


	\begin{abstract}
	Many causal parameters depend on a moment of the joint distribution of potential outcomes.
	Such parameters are especially relevant in policy evaluation settings, where noncompliance is common and accommodated through the model of \cite{imbens1994late}.
	This paper shows that the sharp identified set for these parameters is an interval with endpoints characterized by the value of optimal transport problems.
	Sample analogue estimators are proposed based on the dual problem of optimal transport.
	These estimators are $\sqrt{n}$-consistent and converge in distribution under mild assumptions.
	Inference procedures based on the bootstrap are straightforward and computationally convenient.
	The ideas and estimators are demonstrated in an application revisiting the National Supported Work Demonstration job training program.
	I find suggestive evidence that workers who would see below average earnings without treatment tend to see above average benefits from treatment.
\end{abstract}

\bigskip

\begin{singlespace}
	\begin{center}
		\textbf{Keywords:} potential outcomes, treatment effects, partial identification, bounds, \\optimal transport
	\end{center}
\end{singlespace}

	\pagenumbering{gobble}

	\clearpage

	\pagenumbering{arabic}

	\section{Introduction}
\label{Section: introduction}

Researchers studying the causal effects of a binary treatment see an observation's treated or untreated outcome, but never both.
As a result, the data identify the marginal distributions of each potential outcome, but not their joint distribution.
This ``fundamental problem of causal inference'' \citep{holland1986statistics} leaves parameters depending on the joint distribution partially identified.


In this paper I study a wide class of parameters that depend on a moment of the joint distribution of potential outcomes.
My setting is the canonical potential outcomes framework with binary treatment, a binary instrument satisfying a monotonicity restriction, and finitely supported covariates \citep{imbens1994late, abadie2003semiparametric}.
In this setting, I show the sharp identified set for such parameters is an interval with endpoints characterized by the value of optimal transport problems.
I propose sample analogue estimators based on the dual problem of optimal transport, which facilitates both computation and asymptotic analysis.
Through the functional delta method, I show these estimators converge in distribution allowing for straightforward inference procedures based on the bootstrap.


The proposed estimators are especially attractive due to their wide applicability and computational simplicity.
The class of parameters under study is broad, including the correlation between potential outcomes, the probability of benefitting from treatment, and many more examples discussed in section \ref{Section: setting and parameter class}.
As argued in \cite{heckman1997making}, such parameters are of particular interest to policymakers and economists carrying out econometric policy evaluation.
Noncompliance with the assigned treatment status is common in these settings.
Most studies accomodate noncompliance with the same framework adopted in this paper, and could make use of these estimators with no additional identifying assumptions.
Computing the estimator and constructing confidence sets entails nothing more challenging than solving linear programming problems, for which there are fast and efficient algorithms readily available.


This paper contributes to a large econometrics literature studying parameters of the joint distribution of potential outcomes.
Many papers in this literature focus on a subset of the parameters considered here, especially the cumulative distribution function (cdf) or quantiles of treatment effects \citep{manski1997monotone, heckman1997making, firpo2007efficient, fan2010sharp, fan2012confidence, firpo2019partial, callaway2021bounds, frandsen2021partial}.
This limited focus allows greater use of known analytical expressions when deriving sharp bounds, especially the famed Makarov bounds on the cdf and Fr\'echet-Hoeffding bounds on the joint distribution.
Several recent works develop methods applicable to broad parameters classes by employing procedures that do not require analytical expressions for the identified set.
\cite{russell2021sharp} studies continuous functionals of the joint distribution of discrete potential outcomes, through a computationally intensive (sometimes infeasible) search over all permissible distributions of model primitives.
\cite{fan2023partial} study parameters identified through moment conditions in several incomplete data settings -- including potential outcomes -- by searching over an infinite dimensional space of smooth copulas.
This paper occupies a middle ground: by focusing on parameters that depend on a scalar moment of the joint distribution and working with optimal transport, I obtain expressions for the bounds with tractable sample analogues.
This approach allows consideration of a wide variety of parameters while maintaining computational tractability.

This paper also contributes to a growing literature on applications of optimal transport to econometrics; see \cite{galichon2017survey} for a recent survey. Several recent working papers utilize optimal transport for issues related to casual inference, including inverse propensity weighting \citep{dunipace2021optimal}, matching on covariates \citep{gunsilius2021matching}, and obtaining counterfactual distributions \citep{torous2021optimal}. In concurrent and highly complementary work, \cite{ji2023model} consider a very similar class of parameters to the present paper and also propose inference based on the dual problem of optimal transport. Their focus, accomodating non-discrete covariates without resorting to parametric models, leads to theory based on cross fitting and high-level assumptions on first stage estimators. The goal of the present paper is to provide simple, low-level conditions and computationally convenient estimators in the common case where covariates are discrete. This leads to theory based on Hadamard directional differentiability and the functional delta method quite distinct from that of \cite{ji2023model}.

The remainder of this paper is organized as follows. Section \ref{Section: setting and parameter class} formalizes the setting and introduces the class of parameters under study. Optimal transport is introduced in section \ref{Section: optimal transport}, and used in identification in section \ref{Section: identification}. Section \ref{Section: estimators} proposes the estimators and contains the asymptotic results.
Section \ref{Section: application} contains the application, showing suggestive evidence that the the National Supported Work Demonstration job training program was especially beneficial for workers who would otherwise see below average incomes.
Section \ref{Section: extensions} discusses straightforward extensions, and section \ref{Section: conclusion} concludes.



	\section{Setting and parameter class}
\label{Section: setting and parameter class}

\subsection{Setting}
\label{Section: setting and parameter class, subsection setting}

Consider a potential outcomes framework with binary treatment, a binary instrument, and finitely supported covariates (\cite{imbens1994late}, \cite{abadie2003semiparametric}). Let $Y$ denote the scalar, real-valued outcome of interest and $D \in \{0,1\}$ indicate treatment status. Further let $Y_1$ denote the potential outcome when treated and $Y_0$ the potential outcome when untreated. The observed outcome $Y$ is given by
\begin{equation}
	Y = D Y_1 + (1-D) Y_0. \label{Display: observed and potential outcomes}
\end{equation}
The difference in potential outcomes, $Y_1 - Y_0$, is called the treatment effect.

The binary instrument is denoted $Z \in \{0,1\}$. Let $D_1$ denote the treatment status when $Z = 1$, and $D_0$ the treatment status when $Z = 0$. The observed treatment status $D$ is given by
\begin{equation}
	D = Z D_1 + (1-Z) D_0. \label{Display: observed and potential treatment status}
\end{equation}
It is assumed that the instrument itself does not affect the outcome.\footnote{\protect\linespread{1}\protect\selectfont
	One could hypothesize potential outcomes varying with the value of the instrument, i.e. $Y_{dz}$ for each $(d,z)$. The exposition here implicitly assumes \textit{instrument exclusion}, also known as the \textit{Stable Unit Treatment Value Assumption}: that $P(Y_{d1} = Y_{d0}) = 1$ for each $d$.

}
Units with $1 = D_1 > D_0 = 0$ are known as \textit{compliers}.

Assumption \ref{Assumption: setting} formalizes the setting.
\begin{restatable}[Setting]{assumption}{assumptionSetting}
	\label{Assumption: setting}
	\singlespacing

	$\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$ is an i.i.d. sample with $(Y, D, Z, X) \sim P$,
	\begin{align}
		&Y \in \mathcal{Y} \subseteq \mathbb{R}, &&D \in \{0,1\}, &&Z \in \{0,1\}, &&X \in \mathcal{X} = \{x_1, \ldots, x_M\} \subseteq \mathbb{R}^{d_x} \label{Assumption: setting, support}
	\end{align}
	where $Y$, $D$, and $Z$ are related to $(Y_1, Y_0, D_1, D_0)$ through equations \eqref{Display: observed and potential outcomes} and \eqref{Display: observed and potential treatment status}, and the random vector $(Y_1, Y_0, D_1, D_0, Z, X)$ satisfies
	\begin{enumerate}[label=(\roman*), noitemsep]
		\item Instrument independence: $(Y_1, Y_0, D_1, D_0) \perp Z \mid X$, \label{Assumption: setting, instrument independence}
		\item Monotonicity: $P(D_1 \geq D_0) = 1$, \label{Assumption: setting, monotonicity}
		\item Existence of compliers: $P(D_1 > D_0, X = x) > 0$ for each $x$, and \label{Assumption: setting, existence of compliers}
		\item $P(X = x, Z = z) > 0$ for each $(x,z)$. \label{Assumption: setting, common support}
	\end{enumerate}
\end{restatable}

Assumption \ref{Assumption: setting} is essentially equivalent to assumption 2.1 in \cite{abadie2003semiparametric}, with the addition that covariates are finitely supported.
Instrument independence is sometimes referred to as \textit{ignorability}, and satisfied in most randomized controlled trials, where $Z$ indicates being assigned to treatment. Monotonicity is typically a weak assumption in such settings.

It is worth emphasizing that this setting nests the case where treatment is exogenous. Specifically, when $D_1 = 1$ and $D_0 = 0$ (degenerately), every unit is a complier. In this case equation \eqref{Display: observed and potential treatment status} shows treatment status equals the instrument: $D = Z$. Instrument independence simplifies to $(Y_1, Y_0) \perp D \mid X$, and monotonicity is trivially satisfied.

\subsubsection{Distributions of compliers}
\label{Section: setting and parameter class, subsection setting, subsubsection distributions of compliers}

Interest focuses on the distribution of compliers. Such focus is especially policy relevant when ``the policy is the instrument'' i.e., the proposed change in policy is to assign $Z=1$ to all units. \cite{abadie2003semiparametric} shows that assumption \ref{Assumption: setting} suffices to identify the marginal distributions of $Y_1$ and $Y_0$ for the subpopulation of compliers.
\begin{restatable}[\cite{abadie2003semiparametric}]{lemma}{lemmaLATEIVMarginalDistributionIdentification}
	\label{Lemma: identification, LATE IV marginal distribution identification}
	\singlespacing

	Suppose assumption \ref{Assumption: setting} holds. Then the marginal distributions of $Y_d$ conditional on $D_1 > D_0$ and $X = x$, denoted $P_{d \mid x}$, are identified by
	\begin{align}
		E_{P_{d \mid x}}[f(Y_d)] &\equiv E[f(Y_d) \mid D_1 > D_0, X = x] \notag \\
		&= \frac{E[f(Y) \mathbbm{1}\{D = d\} \mid Z = d, X = x] - E[f(Y) \mathbbm{1}\{D = d\} \mid Z = 1-d, X = x]}{P(D = d \mid Z = d, X = x) - P(D = d \mid Z = 1-d, X = x)} \label{Display: lemma, identification, LATE IV marginal distribution identification, conditional distributions}
	\end{align}
	for any integrable function $f$. Furthermore, the distribution of $X$ conditional on $D_1 > D_0$ is identified by
	\begin{align}
		s_x &\equiv P(X = x \mid D_1 > D_0) \notag \\
		&= \frac{\left[P(D = 1 \mid Z = 1, X = x) - P(D = 1 \mid Z = 0, X = x) \right]P(X = x) }{\sum_{x'} \left[P(D = 1 \mid Z = 1, X = x') - P(D = 1 \mid Z = 0, X = x') \right]P(X = x')} \label{Display: lemma, identification, LATE IV marginal distribution identification, conditional probability X = x}
	\end{align}
\end{restatable}

The joint distribution of potential outcomes is not identified. This is a result of the fundamental problem of causal inference: there is no unit where both $Y_1$ and $Y_0$ are observed, and as a result the joint distribution of $(Y_1, Y_0)$ is not identified for any subpopulation. Let $P_{1,0}$ denote the joint distribution of $(Y_1,Y_0)$ conditional on compliance, and $P_{1, 0 \mid x}$ denote the joint distribution conditional on compliance and $X = x$. These are related through the law of iterated expectations; for any function $c(y_1,y_0)$ with values in $\mathbb{R}$,
\begin{equation*}
	E_{P_{1,0}}[c(Y_1,Y_0)] = E[E[c(Y_1,Y_0) \mid D_1 > D_0, X] \mid D_1 > D_0] = \sum_x s_x E_{P_{1,0 \mid x}}[c(Y_1,Y_0)].
\end{equation*}
This relation can also be expressed as $P_{1,0} = \sum_x s_x P_{1,0 \mid x}$.

A joint distribution with marginals $P_{1 \mid x}$ and $P_{0 \mid x}$ is called a \textit{coupling} of $P_{1 \mid x}$ and $P_{0 \mid x}$. $P_{1,0 \mid x}$ is such a coupling, and is otherwise unrestricted by assumption \ref{Assumption: setting}.
Thus the identified set for $P_{1,0 \mid x}$ is the set of distributions $\pi_{1,0 \mid x}$ for $(Y_1,Y_0)$ with marginals $\pi_{1 \mid x} = P_{1 \mid x}$ and $\pi_{0 \mid x} = P_{0 \mid x}$, denoted
\begin{equation}
	\Pi(P_{1 \mid x}, P_{0 \mid x}) = \left\{\pi_{1,0 \mid x} \; : \; \pi_{1 \mid x} = P_{1 \mid x}, \; \pi_{0 \mid x} = P_{0 \mid x}\right\}. \label{Display: identified set for P10x}
\end{equation}
Moreover, the identified set for $P_{1,0}$ is $\left\{\pi_{1,0} = \sum_x s_x \pi_{1,0 \mid x} \; : \; \pi_{1,0 \mid x} \in \Pi(P_{1 \mid x}, P_{0 \mid x})\right\}$.

\subsection{Parameter class}
\label{Section: setting and parameter class, subsection parameter class}

The idea at the core of this paper is to bound a moment of the joint distribution of potential outcomes by optimization. Accordingly, the focus is on scalar parameters of the form
\begin{equation}
	\gamma = g(\theta, \eta) \label{Display: gamma = g(theta, eta)}
\end{equation}
where $g$ is a known function and $\theta = E_{P_{1,0}}[c(Y_1, Y_0)] \in \mathbb{R}$ is a scalar moment of the joint distribution of $(Y_1,Y_0)$ conditional on compliance. The function $c$ is known, and referred to as a ``cost function'' in connection with the optimal transport literature. This class of parameters is broad, as illustrated by the examples given below. In each of these examples $\eta$ is a finite collection of moments of the marginal distributions conditional on compliers: $\eta = (E_{P_1}[\eta_1(Y_1)], E_{P_0}[\eta_0(Y_0)]) \in \mathbb{R}^{K_1 + K_0}$. The formal results focus on this case, but could be generalized to allow $\eta$ to be other point identified nuisance parameters.

The following conditions are stronger than necessary for identification of the sharp identified set of $\gamma$, but will be used when constructing and studying estimators. Assumption \ref{Assumption: cost function} places restrictions on the cost function to ensure optimal transport can be used characterize and estimate the sharp identified set for $\theta$.
\begin{restatable}[Cost function]{assumption}{assumptionCostFunction}
	\label{Assumption: cost function}
	\singlespacing

	Either
	\begin{enumerate}[label=(\roman*)]
		\item $c(y_1,y_0)$ is Lipschitz continuous and $\mathcal{Y}$ is compact, or  \label{Assumption: cost function, smooth costs}

		\item $c(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$ for a known $\delta \in \mathbb{R}$ and the cumulative distribution functions $F_{d \mid x}(y) = P(Y_d \leq y \mid D_1 > D_0, X = x)$ are continuous. \label{Assumption: cost function, CDF}
	\end{enumerate}
\end{restatable}

Assumption \ref{Assumption: cost function} covers every example listed below. Continuous cost functions $c$ are given a unified analysis, but for reasons discussed in section \ref{Section: optimal transport} discontinuous cost functions must be handled on a case-by-case basis. I focus on the leading case of interest in applications, $c(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$, corresponding to the cumulative distribution of treatment effects. The approach developed in this paper could likely be generalized to cover other discontinuous cost functions; for example, results in the appendix allow estimation of the sharp lower bound of $P((Y_1, Y_0) \in C)$ for any open, convex set $C \subseteq \mathbb{R}^2$.

Assumption \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} requires the cdfs $F_{d \mid x}$ be continuous. As discussed in section \ref{Section: identification}, this ensures the set being estimated is the sharp identified set for the parameter of interest. However, the estimation and inference results of section \ref{Section: estimators} hold \textit{regardless} of whether the cdfs are continuous or not; when the cdfs are not continuous, the estimand is a valid outer identified set.

Under assumptions \ref{Assumption: setting} and \ref{Assumption: cost function}, the sharp identified set for $\theta$ is an interval $[\theta^L, \theta^H]$. Assumption \ref{Assumption: parameter, function of moments} contains conditions on $g$ and $\eta$.

\begin{restatable}[Function of moments]{assumption}{assumptionParameterFunctionOfMoments}
	\label{Assumption: parameter, function of moments}
	\singlespacing

	The parameter is $\gamma = g(\theta, \eta) \in \mathbb{R}$, where
	\begin{align*}
		&\theta = E[c(Y_1, Y_0) \mid D_1 > D_0] \in \mathbb{R}, &&\eta = E
		\begin{bmatrix}
			\eta_1(Y_1), \eta_0(Y_0) \mid D_1 > D_0
		\end{bmatrix}
		\in \mathbb{R}^{K_1 + K_0}
	\end{align*}
	for known functions $g$, $c$, $\eta_1$ and $\eta_0$ such that
	\begin{enumerate}[label=(\roman*)]
		\item $E[\lVert \eta_d(Y)\rVert^2] < \infty$ for $d = 1,0$,  \label{Assumption: parameter, function of moments, nuisance moments have finite variance}
		\item $g(\cdot, \eta)$ is continuous, and \label{Assumption: parameter, function of moments, g is continuous}
		\item the functions
		\begin{align*}
			&g^L(t^L, t^H, e) = \min_{t \in [t^L, t^H]} g(t, e), &&g^H(t^L, t^H, e) = \max_{t \in [t^L, t^H]} g(t, e)
		\end{align*}
		are continuously differentiable at $(t^L, t^H, e) = (\theta^L, \theta^H, \eta)$. \label{Assumption: parameter, function of moments, sup and inf of g are differentiable}
	\end{enumerate}
\end{restatable}
Note that when $\theta$ itself is of interest, assumption \ref{Assumption: parameter, function of moments} is satisfied with $g(\theta, \eta) = \theta$. Assumption \ref{Assumption: parameter, function of moments} \ref{Assumption: parameter, function of moments, g is continuous} ensures the identified set for $\gamma$ is the interval $[\gamma^L, \gamma^H]$, and assumption \ref{Assumption: parameter, function of moments} \ref{Assumption: parameter, function of moments, sup and inf of g are differentiable} is used to apply the delta method. It is straightforward to show assumption \ref{Assumption: parameter, function of moments} \ref{Assumption: parameter, function of moments, sup and inf of g are differentiable} holds when $g$ is continuously differentiable in both arguments and $g(\cdot, \eta)$ is strictly increasing, as the latter condition implies $g^L(\theta^L, \theta^H, \eta) = g(\theta^L, \eta)$ and $g^H(\theta^L, \theta^H, \eta) = g(\theta^H, \eta)$ and the former condition implies they are continuously differentiable. This argument applies to every parameter listed below. When $g$ is differentiable but $g(\cdot, \eta)$ is not monotonic, it is often possible to use the implicit function theorem applied to first order conditions to derive sufficient conditions for the corresponding $\operatorname*{arg\,min}$ and $\operatorname*{arg\,max}$ to be differentiable, and thus for assumption \ref{Assumption: parameter, function of moments} \ref{Assumption: parameter, function of moments, sup and inf of g are differentiable} to hold.

\subsubsection{Examples}
\label{Section: setting and parameter class, subsection parameter class, subsubsection examples}

The following examples are intended both to fix ideas and illustrate the broad scope of the parameter class described above.

\begin{restatable}[Summary statistics]{example}{exampleSummaryStatistics}
	\singlespacing
	Many summary statistics can be rewritten in the form $\gamma = g(\theta,\eta)$. For example, suppose interest is in the variance of treatment effects for compliers: $\gamma = \text{Var}(Y_1 - Y_0 \mid D_1 > D_0)$. This parameter can be rewritten as
	\begin{equation*}
		\gamma = \text{Var}(Y_1 - Y_0 \mid D_1 > D_0) = E_{P_{1,0}}[(Y_1 - Y_0)^2] - (E_{P_1}[Y_1] - E_{P_0}[Y_0])^2,
	\end{equation*}
	This parameter fits the form $\gamma = g(\theta, \eta)$ required of display \eqref{Display: gamma = g(theta, eta)}, with $\theta = E_{P_{1,0}}[(Y_1 - Y_0)^2]$, $\eta = (\eta^{(1)}, \eta^{(2)}) = (E_{P_1}[Y_1], E_{P_0}[Y_0])$, and $g(\theta, \eta) = \theta - (\eta^{(1)} - \eta^{(2)})^2$. The cost function $c(y_1,y_0) = (y_1 - y_0)^2$ satisfies assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} when $\mathcal{Y}$, the support of the outcome $Y$, is bounded.

	Similarly, suppose the researcher is interested in the correlation between $Y_1$ and $Y_0$ for compliers. Set $\gamma = \text{Corr}(Y_1, Y_0 \mid D_1 > D_0)$, which can be rewritten as
	\begin{equation*}
		\gamma = \text{Corr}(Y_1, Y_0 \mid D_1 > D_0) = \frac{E_{P_{1,0}}[Y_1 Y_0] - E_{P_1}[Y_1]E_{P_0}[Y_0]}{\sqrt{E_{P_1}[Y_1^2] - (E_{P_1}[Y_1])^2}\sqrt{E_{P_0}[Y_0^2] - (E_{P_0}[Y_0])^2}}
	\end{equation*}
	This parameter also fits the form $\gamma = g(\theta, \eta)$ in display \eqref{Display: gamma = g(theta, eta)}, with $\theta = E_{P_{1,0}}[Y_1 Y_0]$, $\eta = (\eta^{(1)}, \eta^{(2)}, \eta^{(3)}, \eta^{(4)}) = (E_{P_1}[Y_1], E_{P_1}[Y_1^2], E_{P_0}[Y_0], E_{P_0}[Y_0^2])$, and $g(\theta,\eta) = \frac{\theta - \eta^{(1)} \times \eta^{(3)}}{\sqrt{\eta^{(2)} - (\eta^{(1)})^2}\sqrt{\eta^{(4)} - (\eta^{(3)})^2}}$. The cost function $c(y_1,y_0) = y_1 y_0$ satisfies assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} when $\mathcal{Y}$ is bounded.
\end{restatable}

\begin{restatable}[Expected percent change]{example}{exampleExpectedPercentChange}
	\singlespacing
	The expected percent change in the outcome can be written as $100 \times E\left[\frac{Y_1 - Y_0}{Y_0} \mid D_1 > D_0\right] \%$. This is a unit-invariant causal parameter that is a natural summary measure when $Y_0$ exhibits considerably variation. For example, a treatment effect of $Y_1 - Y_0 = 5$ is typically of greater economic significance when the untreated outcome is small, say $ Y_0 = 10$, than when $Y_0 = 100$.

	The expected percent change is proportional to
	\begin{equation*}
		\gamma = E\left[\frac{Y_1 - Y_0}{Y_0} \mid D_1 > D_0\right] = E_{P_{1,0}}\left[\frac{Y_1 - Y_0}{Y_0}\right],
	\end{equation*}
	which fits the form of display \eqref{Display: gamma = g(theta, eta)}, with $\gamma = \theta =  E_{P_{1,0}}\left[\frac{Y_1 - Y_0}{Y_0}\right]$. The cost function $c(y_1,y_0) = \frac{y_1 - y_0}{y_0}$ satisfies assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} when $\mathcal{Y}$ is bounded and bounded away from zero.
\end{restatable}

\begin{restatable}[Equitable policies]{example}{exampleCovarianceTreatmentEffectUntreatedOutcome}
	\label{Example: equitable policies}
	\singlespacing

	Policy makers are often interested in whether a policy is equitable -- that is, whether the benefits are concentrated among those who would have undesirable outcomes without treatment.

	One parameter that speaks to these concerns is the covariance between treatment effects and untreated outcomes among compliers: $\gamma = \text{Cov}(Y_1 - Y_0, Y_0\mid D_1 > D_0)$. Notice that $\gamma < 0$ implies those with below average $Y_0$ tend to see above average treatment effects. This parameter can be rewritten as
	\begin{equation*}
		\gamma = \text{Cov}(Y_1 - Y_0, Y_0\mid D_1 > D_0) = E_{P_{1,0}}[(Y_1 - Y_0) Y_0] - (E_{P_1}[Y_1] - E_{P_0}[Y_0])E_{P_0}[Y_0]
	\end{equation*}
	and fits the form $g(\theta, \eta)$ with $\theta = E_{P_{1,0}}[(Y_1 - Y_0)Y_0]$, $\eta = (E_{P_1}[Y_1], E_{P_0}[Y_0])$, and $g(\theta,\eta) = \theta - (\eta^{(1)} - \eta^{(2)})\eta^{(2)}$. The cost function $c(y_1,y_0) = (y_1 - y_0)y_0$ satisfies assumpion \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} when $\mathcal{Y}$ is bounded.

	Many related parameters share a sign with $\text{Cov}(Y_1 - Y_0, Y_0 \mid D_1 > D_0) $ and are also suitable for such an analysis. One such example is the OLS slope when regressing $Y_1 - Y_0$ on $Y_0$ and a constant: $\gamma = \frac{\text{Cov}(Y_1 - Y_0, Y_0 \mid D_1 > D_0)}{\text{Var}(Y_0 \mid D_1 > D_0)}$. This parameter can be rewritten as
	\begin{equation*}
		\gamma = \frac{\text{Cov}(Y_1 - Y_0, Y_0 \mid D_1 > D_0)}{\text{Var}(Y_0 \mid D_1 > D_0)} = \frac{E_{P_{1,0}}[(Y_1 - Y_0) Y_0] - (E_{P_1}[Y_1] - E_{P_0}[Y_0])E_{P_0}[Y_0]}{E_{P_0}[Y_0^2] - (E_{P_0}[Y_0])^2}
	\end{equation*}
	where $\theta = E_{P_{1,0}}[(Y_1 - Y_0)Y_0]$, $\eta = (E_{P_1}[Y_1], E_{P_0}[Y_0], E_{P_0}[Y_0^2])$, and $g(\theta,\eta) = \frac{\theta - (\eta^{(1)} - \eta^{(2)})\eta^{(2)}}{\eta^{(3)} - (\eta^{(2)})^2}$.
\end{restatable}

\begin{restatable}[Proportion that benefit]{example}{exampleShareBenefiting}
	\singlespacing
	The share of compliers benefiting from treatment, written
	\begin{equation*}
		\gamma = P(Y_1 > Y_0 \mid D_1 > D_0),
	\end{equation*}
	is naturally of interest in applications where theory gives little indication whether the treatment will have a positive or negative effect. For example, \cite{allcott2020welfare} study the effect of deactivating facebook on subjective well-being. The authors find significant positive average effects of deactivation, but find substantial heterogeneity in follow-up interviews.

	This parameter fits the form of display \eqref{Display: gamma = g(theta, eta)}, with $\gamma = \theta = E_{P_{1,0}}[\mathbbm{1}\{Y_1 - Y_0 \leq 0\}]$. The cost function $c(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 \leq 0\}$ satisfies assumption \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} if the cdfs $F_{d \mid x}(y)$ are continuous.

	The share benefiting from treatment is also of particular interest when the intervention comes at a financial cost and the outcome of interest is a pecuniary return. Common examples include job training programs intended to increase a worker's income (e.g. the National Supported Work Demonstration studied in \cite{couch1992long}) or management practices intended to raise a firm's accounting profit (e.g. the employee referral program studied in \cite{friebel2023employee}). To illustrate, suppose the researcher observes $\{R_i, C_i, D_i, Z_i\}_{i=1}^n$, where $R$ is observed revenue and $C$ is the observed cost. These are related to treatment status $D \in \{0,1\}$, potential revenues $(R_1, R_0)$, and potential costs $(C_1, C_0)$ by
	\begin{align*}
		&R = D R_1 + (1-D)R_0, &&C = D C_1 + (1-D) C_0
	\end{align*}
	The observed profit, $Y = R - C$, is related to treatment status by
	\begin{align*}
		Y = D \underbrace{(R_1 - C_1)}_{\coloneqq Y_1} + (1-D)\underbrace{(R_0 - C_0)}_{\coloneqq Y_0}
	\end{align*}
	The probability the change in revenue exceeds the change in cost is
	\begin{align*}
		P(R_1 - R_0 > C_1 - C_0 \mid D_1 > D_0) = P(Y_1 > Y_0 \mid D_1 > D_0)
	\end{align*}
\end{restatable}

\begin{restatable}[Quantiles]{example}{exampleQuantilesTreatmentEffects}
\singlespacing
\label{Example: pseudo quantiles}

Suppose the parameter of interest is any $q_\tau$ solving
\begin{equation}
	P(Y_1 - Y_0 \leq q_\tau) = \tau \label{Display: pseudo quantile definition}
\end{equation}
This parameter has a similar interpretation to the $\tau$-th quantile.\footnote{The $\tau$-th quantile is usually defined as the unique  value $\tilde{q}_\tau = \inf\{y \; ; \; P(Y_1 - Y_0 \leq y) \geq \tau\}$. When the $\tau$ level set of the cumulative distribution function $P(Y_1 - Y_0 \leq \cdot)$ is nonempty, the $\tau$-th quantile has the interpretation that $100 \times \tau \%$ of the population has treatment effect less than or equal to $\tilde{q}_\tau$. Every $q_\tau$ solving \eqref{Display: pseudo quantile definition} has the same interpretation.}
$q_\tau$ cannot be viewed as $\gamma = g(\theta,\eta)$. However, by viewing $\theta(\delta) = P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0) = E_{P_{1,0}}[\mathbbm{1}\{Y_1 - Y_0 \leq \delta\}]$ as a function of $\delta$, the results below can be adapted to construct a confidence set for the identified set of this parameter as described in section \ref{Section: extensions, subsection quantiles}.
\end{restatable}





	\section{Optimal Transport}
\label{Section: optimal transport}

This section defines and discusses optimal transport, which is used to characterize the identified set and construct estimators.

Given any marginal distributions $P_1$ and $P_0$ and a ``cost function'' $c(y_1, y_0)$, the Monge-Kantorovich formulation of \textbf{optimal transport} is the problem of choosing a coupling $\pi \in \Pi(P_1,P_0)$ to minimize $E_\pi[c(Y_1,Y_0)]$:
\begin{equation}
	OT_c(P_1, P_0) = \inf_{\pi \in \Pi(P_1, P_0)} E_\pi[c(Y_1,Y_0)]. \label{Defn: optimal transport primal problem, main text}
\end{equation}
This minimization problem in \eqref{Defn: optimal transport primal problem, main text} is referred to as the \textbf{primal problem}, and will be used to characterize the identified set of $\theta$.

The dual problem of optimal transport will be used to construct and analyze estimators. Let $\Phi_c$ denote the set of functions $\varphi(y_1)$ and $\psi(y_0)$ whose pointwise sum is less than $c(y_1,y_0)$:
\begin{equation}
	\Phi_c = \left\{(\varphi,\psi) \; ; \; \varphi(y_1) + \psi(y_0) \leq c(y_1, y_0)\right\}. \label{Defn: Phi_c, main text}
\end{equation}
The \textbf{dual problem} chooses a pair of functions in $\Phi_c$ to maximize the sum of the corresponding expectations:
\begin{equation}
	\sup_{(\varphi,\psi) \in \Phi_c} E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)]. \label{Display: optimal transport dual problem, main text}
\end{equation}
When the cost function is lower semicontinuous and bounded from below, the primal problem is attained and \textbf{strong duality} holds:
\begin{equation}
	OT_c(P_1, P_0) = \min_{\pi \in \Pi(P_1, P_0)} E_\pi[c(Y_1, Y_0)] = \sup_{(\varphi,\psi) \in \Phi_c} E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)]. \label{Display: strong duality, main text}
\end{equation}
The dual problem will be used to construct and analyze estimators. Indeed, the identification of $P_{d \mid x}$ in lemma \ref{Lemma: identification, LATE IV marginal distribution identification} suggests straightforward sample analogues estimating $E_{P_{d \mid x}}[f(Y_d)]$ for a given $f$, which makes it possible to form a sample analogue of the dual problem.

Although it is clear how to form a sample analogue of the dual problem, it is not immediately clear how to analyze the resulting estimator. Fortunately, the dual problem can often be simplified by restricting the maximization problem to a smaller set of functions. Estimators based on this restricted dual problem can then be studied with empirical process techniques.

The dual feasible set is restricted with the concept of $c$-concavity. Notice the dual problem's objective is monotonic, in the sense that $\varphi(y_1) \leq \tilde{\varphi}(y_1)$ for all $y_1$ implies
\begin{align*}
	E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)] \leq E_{P_1}[\tilde{\varphi}(Y_1)] + E_{P_0}[\psi(Y_0)].
\end{align*}
Increasing $\psi$ pointwise will also increase the dual objective. Speaking loosely, any function pair $(\varphi,\psi) \in \Phi_c$ for which the constraint $\varphi(y_1) + \psi(y_0) \leq c(y_1,y_0)$ is ``slack'' cannot be a solution to the dual problem and can therefore be ignored. This motivates the definition of the \textbf{$c$-transforms} of a function $\varphi$:
\begin{align*}
	&\varphi^c(y_0) = \inf_{y_1} \{c(y_1, y_0) - \varphi(y_1)\}, &&\varphi^{cc}(y_1) = \inf_{y_0} \{c(y_1,y_0) - \varphi^c(y_0)\}.
\end{align*}
For any pair of functions $(\varphi, \psi) \in \Phi_c$, these definitions imply $\psi(y_0) \leq \varphi^c(y_0)$, $\varphi(y_1) \leq \varphi^{cc}(y_1)$, and $\varphi^{cc}(y_1) + \varphi^c(y_0) \leq c(y_1,y_0)$. Further $c$-transformations are irrelevant because $(\varphi^{cc})^c = \varphi^c$, so a function $\varphi$ is called \textbf{$c$-concave} if $\varphi^{cc} = \varphi$. If the $c$-transforms are integrable, the dual problem can be restricted to $c$-concave conjugate pairs, $(\varphi^{cc}, \varphi^c)$. Furthermore, $c$-concave functions often ``inherit'' properties of the cost function $c$; for example, if $c$ is Lipschitz continuous then $\varphi^c$ and $\varphi^{cc}$ are Lipschitz continuous as well. These properties can be used to define sets of functions $\mathcal{F}_c$ and $\mathcal{F}_c^c$ (depending on the cost function $c$ but not on the distributions $P_1$, $P_0$) such that
\begin{equation}
	\sup_{(\varphi,\psi) \in \Phi_c} E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)] =  \sup_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)]. \label{Display: strong duality, with smaller feasible set}
\end{equation}

Two cases suffice for the parameters considered in this paper. When the cost function $c(y_1,y_0)$ is Lipschitz continuous and $\mathcal{Y}$ is compact, define
\begin{align}
	\mathcal{F}_c &= \left\{\varphi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \;  -\lVert c \rVert_\infty \leq \varphi(y_1) \leq \lVert c \rVert_\infty, \; \lvert \varphi(y_1) - \varphi(y_1') \rvert \leq L \lvert y_1 - y_1'\rvert \right\} \label{Defn: F_c for smooth costs} \\
	\mathcal{F}_c^c &= \left\{\psi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; -2\lVert c \rVert_\infty \leq \psi(y_0) \leq 0, \; \lvert \psi(y_0) - \psi(y_0') \rvert \leq L \lvert y_0 - y_0'\rvert \right\} \label{Defn: F_c^c for smooth costs}
\end{align}
where $\lVert c \rVert_\infty = \sup_{(y_1,y_0)} \lvert c(y_1, y_0) \rvert$ and $L$ is the Lipschitz constant of $c$. When $c(y_1,y_0) = \mathbbm{1}\{(y_1,y_0) \in C\}$ for an open, convex set $C$, let
\begin{align}
	\mathcal{F}_c &= \left\{\varphi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; \varphi(y_1) = \mathbbm{1}\{y_1 \in I\} \text{ for some interval } I\right\} \label{Defn: F_c for indicator costs of convex C} \\
	\mathcal{F}_c^c &= \left\{\psi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; \psi(y_0) = -\mathbbm{1}\{y_0 \in I^c\} \text{ for some  interval } I\right\} \label{Defn: F_c^c for indicator costs of convex C}
\end{align}

Equation \eqref{Display: strong duality, with smaller feasible set} shows the optimal transport functional $OT_c(P_1, P_0)$ depends only on the values of $E_{P_1}[\varphi(Y_1)]$ and $ E_{P_0}[\psi(Y_0)]$ for $(\varphi,\psi) \in \mathcal{F}_c \times \mathcal{F}_c^c$. For any set $A$, let $\ell^\infty(A)$ denote the space of real-valued bounded functions defined on $A$, equipped with the supremum norm: $\ell^\infty(A) = \left\{f : A \rightarrow \mathbb{R} \; ; \; \lVert f \rVert_\infty = \sup_{a \in A} \lvert f(a) \rvert < \infty \right\}$.
Optimal transport can be viewed as the map $OT_c : \ell^\infty(\mathcal{F}_c) \times \ell^\infty(\mathcal{F}_c^c) \rightarrow \mathbb{R}$ given by
\begin{equation}
	OT_c(P_1, P_0) = \sup_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} E_{P_1}[\varphi(Y_1)] + E_{P_0}[\psi(Y_0)]. \label{Defn: optimal transport dual formulation, main text optimal transport}
\end{equation}
This problem will be referred to as the \textbf{restricted dual problem}. Estimators formed with this map can be studied with empirical process techniques.

In summary, $OT_c(P_1, P_0)$ will be viewed as the functional in \eqref{Defn: optimal transport primal problem, main text} when considering identification, and as the functional given in \eqref{Defn: optimal transport dual formulation, main text optimal transport} when considering estimation. By ensuring $c$ is either Lipschitz continuous or the indicator of an open convex set, strong duality and $c$-concavity ensures these functionals agree on the space of probability distributions.



	\section{Identification}
\label{Section: identification}

Recall the parameter of interest is $\gamma = g(\theta, \eta)$, where $\eta$ is a point identified parameter, $\theta = E_{P_{1,0}}[c(Y_1,Y_0)] \in \mathbb{R}$, and $g$ and $c$ are known functions.

Begin by rewriting $\theta = E_{P_{1, 0}}[c(Y_1,Y_0)] = E[c(Y_1,Y_0) \mid D_1 > D_0]$ with the law of iterated expectations:
\begin{equation*}
	\theta = E[E[c(Y_1,Y_0) \mid D_1 > D_0, X] \mid D_1 > D_0] = E[\theta_X \mid D_1 > D_0] = \sum_x s_x \theta_x
\end{equation*}
where $s_x = P(X = x \mid D_1 > D_0)$ and $\theta_x = E[c(Y_1,Y_0) \mid D_1 > D_0, X = x] = E_{P_{1,0 \mid x}}[c(Y_1,Y_0)]$. As noted in section \ref{Section: setting and parameter class, subsection setting, subsubsection distributions of compliers}, the identified set for $P_{1,0 \mid x}$ is the set of couplings of $P_{1\mid x}$ and $P_{0 \mid x}$, denoted $\Pi(P_{1 \mid x},P_{0 \mid x})$. Thus the identified set for $\theta_x$ is $\Theta_{I,x} = \left\{t \in \mathbb{R} \; : \; t = E_\pi[c(Y_1,Y_0)] \text{ for some } \pi \in \Pi(P_{1 \mid x},P_{0 \mid x})\right\}$. $\Pi(P_{1 \mid x},P_{0 \mid x})$ is convex, implying that $\Theta_{I,x}$ is an interval. Let $\theta_x^L$ and $\theta_x^H$ denote its lower and upper endpoint respectively.

To ensure the restricted dual problem can be used for estimation, $\theta_x^L$ and $\theta_x^H$ are characterized through an optimal transport problem with a suitable cost function $c$. When assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} holds ($c(y_1,y_0)$ is Lipschitz continuous and $\mathcal{Y}$ is compact), define
\begin{align}
	&c_L(y_1, y_0) = c(y_1,y_0), &&c_H(y_1, y_0) = -c(y_1,y_0) \notag \\
	&\theta^L(P_{1 \mid x}, P_{0 \mid x}) = OT_{c_L}(P_{1 \mid x}, P_{0 \mid x}), &&\theta^H(P_{1 \mid x}, P_{0 \mid x}) = -OT_{c_H}(P_{1 \mid x}, P_{0 \mid x}). \label{Display: thetaL, thetaH when c is continuous}
\end{align}
Note that $\theta_x^L = \theta^L(P_{1 \mid x}, P_{0 \mid x})$ and $\theta_x^H = \theta^H(P_{1 \mid x}, P_{0 \mid x})$.

The cumulative distribution function of $Y_1 - Y_0$ corresponds to the cost function $c(y_1, y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$, which is not lower semicontinuous. This challenge is circumvented by a small change in the cost function. When assumption \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} holds (the cost function is $c(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$) define
\begin{align}
	&c_L(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 < \delta\}, &&c_H = \mathbbm{1}\{y_1 - y_0 > \delta\} \notag \\
	&\theta^L(P_{1 \mid x}, P_{0 \mid x}) = OT_{c_L}(P_{1 \mid x}, P_{0 \mid x}), &&\theta^H(P_{1 \mid x}, P_{0 \mid x}) = 1 - OT_{c_H}(P_{1 \mid x}, P_{0 \mid x}) \label{Display: thetaL, thetaH when c is for CDF}
\end{align}
It follows from definitions that $\theta_x^H = \theta^H(P_{1 \mid x}, P_{0 \mid x})$. Moreover, $c_L(y_1,y_0) \leq c(y_1,y_0)$ implies $\theta^L(P_{1 \mid x}, P_{0 \mid x})$ is a valid lower bound for $\theta_x$. It is sharp if $P_{1 \mid x}$, $P_{0\mid x}$ have continuous cumulative distribution functions, in which case $\theta_x^L =\theta^L(P_{1 \mid x}, P_{0 \mid x})$. It is worth emphasizing again that the estimation and inference results of section \ref{Section: estimators} hold \textit{regardless} of whether the cdfs are continuous or not; when the cdfs are not continuous, the estimand is a valid outer identified set.

Under assumptions \ref{Assumption: setting} and \ref{Assumption: cost function}, the identified set for $\theta = E_{P_{1,0}}[c(Y_1,Y_0)] = E[c(Y_1,Y_0) \mid D_1 > D_0]$ is the compact interval $[\theta^L, \theta^H]$ with endpoints
\begin{align*}
	&\theta^L = E[\theta_X^L \mid D_1 > D_0] = \sum_x s_x \theta_x^L, &&\theta^H = E[\theta_X^H \mid D_1 > D_0] = \sum_x s_x \theta_x^H
\end{align*}
Under assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments}, the identified set for $\gamma$ is $[\gamma^L, \gamma^H]$, with endpoints
\begin{align}
	&\gamma^L = g^L(\theta^L, \theta^H, \eta) = \inf_{t \in [\theta^L, \theta^H]} g(t, \eta), &&\gamma^H = g^H(\theta^L, \theta^H, \eta) = \sup_{t \in [\theta^L, \theta^H]} g(t, \eta) \label{Display: main idea, gamma bounds}
\end{align}
The following theorem summarizes the discussion above. Let $\theta^L(\cdot, \cdot)$ and $\theta^H(\cdot, \cdot)$ be given by \eqref{Display: thetaL, thetaH when c is continuous} or \eqref{Display: thetaL, thetaH when c is for CDF} depending on the cost function, and set
\begin{align}
	&\theta_x^L = \theta^L(P_{1 \mid x}, P_{0 \mid x}), &&\theta_x^H = \theta^H(P_{1 \mid x}, P_{0 \mid x}), \label{Defn: theta_x^L, theta_x^H formal definition} \\
	&\theta^L = \sum_x s_x \theta_x^L, &&\theta^H = \sum_x s_x \theta_x^H, \label{Defn: theta^L, theta^H formal definition} \\
	&\gamma^L = g^L(\theta^L, \theta^H, \eta), &&\gamma^H = g^H(\theta^L, \theta^H, \eta) \label{Defn: gamma^L, gamma^H formal definition}
\end{align}


\begin{restatable}[Identification of functions of moments]{theorem}{theoremIdentificationFunctionOfMoments}
	\label{Theorem: identification, function of moments}
	\singlespacing

	Suppose assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments} are satisfied. Then the sharp identified set for $\gamma$ is $[\gamma^L, \gamma^H]$.

\end{restatable}
All results are proven in the appendix.

It is worth pausing to consider the role of covariates. When covariates are available, ignoring them leads to wider bounds that are not sharp. Specifically, the marginal distributions $P_1$ and $P_0$ could be used to form a lower bound on $\theta$ with $\theta^L(P_1, P_0) = \inf_{\pi \in \Pi(P_1,P_0)} E_\pi[c_L(Y_1,Y_0)]$. This bound minimizes over the whole set $\Pi(P_1, P_0) = \left\{\pi_{1,0} \; ; \; \pi_1 = P_1, \pi_0 = P_0\right\}$, but the identified set for $P_{1,0}$ is the subset of $\Pi(P_1, P_0)$ given by $\left\{\pi_{1,0} = \sum_x s_x \pi_{1, 0 \mid x} \; ; \; \pi_{1,0 \mid x} \in \Pi(P_{1 \mid x}, P_{0 \mid x})\right\}$. The bound defined through equations \eqref{Defn: theta_x^L, theta_x^H formal definition} and \eqref{Defn: theta^L, theta^H formal definition} is found while enforcing the additional constraints that $\pi_{1,0 \mid x} \in \Pi(P_{1 \mid x}, P_{0 \mid x})$ for each $x$. These additional constraints imply $\theta^L(P_1, P_0) \leq \theta^L$,
\begin{comment}
\begin{align}
	\theta^L(P_1, P_0) &= \inf_{\pi_{1,0} \in \Pi(P_1, P_0)} E_{\pi_{1,0}}[c_L(Y_1,Y_0)] \notag \\
	&\leq \sum_x s_x \times \left(\inf_{\pi_{1,0 \mid x} \in \Pi(P_{1 \mid x}, P_{0 \mid x})} E_{\pi_{1,0 \mid x}}[c_L(Y_1,Y_0)]\right) = \theta^L \label{Display: covariates shrink the identified set}
\end{align}
\end{comment}
and similarly $\theta^H \leq \theta^H(P_1, P_0)$.

Extreme cases illustrate when covariates are informative. If $X$ is independent of $(Y_1,Y_0)$ conditional on $D_1 > D_0$, then $P_{d \mid x} = P_d$ for each $x$, $\Pi(P_{1 \mid x}, P_{0 \mid x}) = \Pi(P_1, P_0)$, and the inequalities above hold as equalities. On the other hand, if $P_{d \mid x}$ is degenerate for either $d = 1$ or $d = 0$, then there is only one possible coupling of $P_{1 \mid x}$ and $P_{0 \mid x}$. Since $\Pi(P_{1 \mid x}, P_{0 \mid x})$ is a singleton, $\theta_x^L = \theta_x^H$ and $\theta_x = E[c(Y_1,Y_0) \mid D_1 > D_0, X = x]$ is point identified. If this occurs for all $x \in \mathcal{X}$, $\theta$ and $\gamma$ are point identified.



\begin{remark}[Makarov bounds]
	\singlespacing
	\label{Remark: Makarov bounds}

	The proof of theorem \ref{Theorem: identification, function of moments} given in the appendix uses properties of optimal transport to argue that under assumptions \ref{Assumption: setting} and \ref{Assumption: cost function} \ref{Assumption: cost function, CDF}, $[\theta^L, \theta^H]$ is the sharp identified set for $P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0)$. Nonetheless, it is interesting to note that the proof shows
	\begin{align*}
		\theta_x^L &= OT_{c_L}(P_{1 \mid x}, P_{0 \mid x}) = \sup_y \{F_{1 \mid x}(y) - F_{0 \mid x}(y - \delta)\} \\
		\theta_x^H &= 1 - OT_{c_H}(P_{1 \mid x}, P_{0 \mid x}) = 1 - \sup_y\{F_{0 \mid x}(y - \delta) - F_{1 \mid x}(y)\} = 1 + \inf_y\{F_{1\mid x}(y) - F_{0 \mid x}(y - \delta)\}
	\end{align*}
	which are the Makarov bounds on $P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0, X = x)$ studied in \cite{fan2010sharp}.
\end{remark}

\begin{remark}[Pointwise vs. uniformly sharp CDF bounds]
	\label{Remark: Pointwise vs. uniformly sharp CDF bounds}
	\singlespacing

	Under assumptions \ref{Assumption: setting} and \ref{Assumption: cost function} \ref{Assumption: cost function, CDF}, $[\theta^L, \theta^H]$ is the sharp identified set for $P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0)$ at the \textit{point} $\delta$. Viewing these bounds as functions of $\delta$, $\theta^L(\delta)$ and $\theta^H(\delta)$ are not \textit{uniformly} sharp bounds for the cumulative distribution function $P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0)$, in the sense that not every CDF $F(\cdot)$ satisfying $\theta^L(\delta) \leq F(\delta) \leq \theta^H(\delta)$ for all $\delta$ could be the CDF of $Y_1 - Y_0$. See \cite{firpo2019partial} for a detailed discussion of this point.
\end{remark}


	\section{Estimators}
\label{Section: estimators}

Sample analogues of the expressions identifying $P_{1 \mid x}$, $P_{0 \mid x}$, and $s_x$ in lemma \ref{Lemma: identification, LATE IV marginal distribution identification} provide convenient plug-in estimators of $\gamma^L$ and $\gamma^H$.

The following notation simplifies expressions for the sample analogues. Let $P$ denote the distribution of an observation $(Y, D, Z, X)$, and $f$ be a real-valued function. Use $P(f)$ to mean $E_P[f(Y, D, Z, X)]$. Similarly, let $P_{d \mid x}(f) = E_{P_{d \mid x}}[f(Y_d)] = E[f(Y_d) \mid D_1 > D_0, X = x]$. Let $\mathbb{P}_n$ denote the empirical distribution formed from the sample $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$, and $\mathbb{P}_n(f) = \frac{1}{n}\sum_{i=1}^n f(Y_i, D_i, Z_i, X_i)$. The following indicator function notation also simplifies expressions:
\begin{gather*}
	\mathbbm{1}_{d,x,z}(D,X,Z) = \mathbbm{1}\{D = d, X = x, Z = z\}, \\
	\mathbbm{1}_{x,z}(X,Z) = \mathbbm{1}\{X = x, Z = z\}, \quad \quad \quad \mathbbm{1}_x(X) = \mathbbm{1}\{X = x\}
\end{gather*}
For example, $P(D = d, X = x, Z = z)$ shortens to $P(\mathbbm{1}_{d,x,z})$, and $\frac{1}{n}\sum_{i=1}^n \mathbbm{1}\{D_i = 1, X_i = x, Z_i = 0\}$ to $\mathbb{P}_n(\mathbbm{1}_{1,x,0})$.

The probabilities $p_{d,x,z} = P(\mathbbm{1}_{d,x,z})$, $p_{x,z} = P(\mathbbm{1}_{x,z})$, and $p_x = P(\mathbbm{1}_x)$ are estimated with empirical analogues:
\begin{align*}
	&\hat{p}_{d,x,z} = \mathbb{P}_n(\mathbbm{1}_{d,x,z}), &&\hat{p}_{x,z} = \mathbb{P}_n(\mathbbm{1}_{x,z}), &&\hat{p}_x = \mathbb{P}_n(\mathbbm{1}_x)
\end{align*}
In this notation, $s_x = P(X = x \mid D_1 > D_0)$ and its empirical analogue $\hat{s}_x$ are
\begin{align}
	&s_x = \frac{(p_{1,x,1}/p_{x,1} - p_{1,x,0}/p_{x,0})p_x}{\sum_{x'} (p_{1,x',1}/p_{x',1} - p_{1,x',0}/p_{x',0})p_x'}, &&\hat{s}_x = \frac{(\hat{p}_{1,x,1}/\hat{p}_{x,1} - \hat{p}_{1,x,0}/\hat{p}_{x,0})\hat{p}_x}{\sum_{x'} (\hat{p}_{1,x',1}/\hat{p}_{x',1} - \hat{p}_{1,x',0}/\hat{p}_{x',0})\hat{p}_{x'}} \label{Display: T_1, map to conditional distributions s_x, main text}
\end{align}
The maps $P_{d \mid x}$ and their empirical analogues are
\begin{align}
	P_{d \mid x}(f) &= \frac{P(\mathbbm{1}_{d,x,d} \times f)/p_{x,d} - P(\mathbbm{1}_{d,x,1-d} \times f)/p_{x,1-d}}{p_{d,x,d}/p_{x,d} - p_{d,x,1-d}/p_{x,1-d}} \notag \\
	\hat{P}_{d \mid x}(f) &= \frac{\mathbb{P}_n(\mathbbm{1}_{d,x,d} \times f)/\hat{p}_{x,d} - \mathbb{P}_n(\mathbbm{1}_{d,x,1-d} \times f)/\hat{p}_{x,1-d}}{\hat{p}_{d,x,d}/\hat{p}_{x,d} - \hat{p}_{d,x,1-d}/\hat{p}_{x,1-d}} \label{Display: T_1, map to conditional distributions, main text}
\end{align}
Under assumption \ref{Assumption: parameter, function of moments}, $\eta = (\eta_1, \eta_0) = (E_{P_1}[\eta_1(Y_1)], E_{P_0}[\eta_0(Y_0)])$. Each vector $\eta_d \in \mathbb{R}^{K_d}$ has coordinates $\eta_d^{(k)} = \sum_x s_x P_{d \mid x}(\eta_d^{(k)})$. Empirical analogues $\hat{\eta} = (\hat{\eta}_1, \hat{\eta}_0)$ are formed by $\hat{\eta}_d^{(k)} = \sum_x \hat{s}_x \hat{P}_{d \mid x}(\eta_d^{(k)})$.

Computing $\hat{P}_{d \mid x}(f)$ for a known $f$ is straightforward:
\begin{align*}
	\hat{P}_{d \mid x}(f) &= \frac{\frac{1}{\hat{p}_{x,d}}\frac{1}{n}\sum_{i=1}^n \mathbbm{1}_{d,x,d}(D_i, X_i, Z_i) f(Y_i) - \frac{1}{\hat{p}_{x,1-d}}\frac{1}{n}\sum_{i=1}^n \mathbbm{1}_{d,x,1-d}(D_i,X_i,Z_i)f(Y_i)}{\hat{p}_{d,x,d}/\hat{p}_{x,d} - \hat{p}_{d,x,1-d}/\hat{p}_{x,1-d}} \\
	&= \sum_{i=1}^n \omega_{d,x,i} \times f_i
\end{align*}
where $f_i = f(Y_i)$ and the weights $\omega_{d,x,i}$ can be computed directly from data:
\begin{equation}
	\omega_{d,x,i} = \frac{1}{n} \times \frac{\mathbbm{1}_{d,x,d}(D_i, X_i, Z_i)/\hat{p}_{x,d} - \mathbbm{1}_{d,x,1-d}(D_i, X_i, Z_i)/\hat{p}_{x,1-d}}{\hat{p}_{d,x,d}/\hat{p}_{x,d} - \hat{p}_{d,x,1-d}/\hat{p}_{x,1-d}} \label{Display: weights to compute Pdx(f)}
\end{equation}

Sample analogue estimators of $\gamma^L$ and $\gamma^H$ are based on equations \eqref{Display: thetaL, thetaH when c is continuous}, \eqref{Display: thetaL, thetaH when c is for CDF}, \eqref{Defn: theta_x^L, theta_x^H formal definition}, \eqref{Defn: theta^L, theta^H formal definition}, and \eqref{Defn: gamma^L, gamma^H formal definition}. These expressions involve the optimal transport functional $OT_c(P_{1 \mid x}, P_{0 \mid x})$. The sample analogue of the simplified dual problem discussed in section \ref{Section: optimal transport} is written
\begin{equation}
	OT_c(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) = \sup_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} \hat{P}_{1 \mid x}(\varphi) + \hat{P}_{0 \mid x}(\psi) \label{Defn: optimal transport dual formulation, main text estimators}
\end{equation}
Here $\mathcal{F}_c$, $\mathcal{F}_c^c$, and the functions $\theta^L(\cdot)$, $\theta^H(\cdot)$ are defined according to the cost function:
\begin{enumerate}[label=(\roman*)]
	\item When assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} holds (the cost function $c(y_1,y_0)$ is Lipschitz continuous and $\mathcal{Y}$ is compact), $\mathcal{F}_c$ and $\mathcal{F}_c^c$ are given by:
	\begin{align*}
		\mathcal{F}_c &= \left\{\varphi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \;  -\lVert c \rVert_\infty \leq \varphi(y_1) \leq \lVert c \rVert_\infty, \; \lvert \varphi(y_1) - \varphi(y_1') \rvert \leq L \lvert y_1 - y_1'\rvert \right\} \\
		\mathcal{F}_c^c &= \left\{\psi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; -2\lVert c \rVert_\infty \leq \psi(y_0) \leq 0, \; \lvert \psi(y_0) - \psi(y_0') \rvert \leq L \lvert y_0 - y_0'\rvert \right\}
	\end{align*}
	and $\theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$, $\theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$ are analogues of equation \eqref{Display: thetaL, thetaH when c is continuous}:
	\begin{align*}
		&c_L(y_1, y_0) = c(y_1,y_0), &&c_H(y_1, y_0) = -c(y_1,y_0) \\
		&\theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) = OT_{c_L}(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}), &&\theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) = -OT_{c_H}(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}).
	\end{align*}

	\item When assumption \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} holds (the cost function is $c(y_1, y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$), $\mathcal{F}_c$ and $\mathcal{F}_c^c$  are given by:
	\begin{align*}
		\mathcal{F}_c &= \left\{\varphi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; \varphi(y_1) = \mathbbm{1}\{y_1 \in I\} \text{ for some interval } I\right\} \\
		\mathcal{F}_c^c &= \left\{\psi : \mathcal{Y} \rightarrow \mathbb{R} \; ; \; \psi(y_0) = -\mathbbm{1}\{y_0 \in I^c\} \text{ for some  interval } I\right\}
	\end{align*}
	and $\theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$, $\theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$ are analogues of equation \eqref{Display: thetaL, thetaH when c is for CDF}:
	\begin{align*}
		&c_L(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 < \delta\}, &&c_H = \mathbbm{1}\{y_1 - y_0 > \delta\} \\
		&\theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) = OT_{c_L}(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}), &&\theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) = 1 - OT_{c_H}(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})
	\end{align*}
\end{enumerate}
The sample analogue estimators are given by
\begin{align}
	&\hat{\theta}_x^L = \theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}), &&\hat{\theta}_x^H = \theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}), \label{Defn: theta_x^L, theta_x^H estimators} \\
	&\hat{\theta}^L = \sum_x \hat{s}_x \hat{\theta}_x^L, &&\hat{\theta}^H = \sum_x \hat{s}_x \hat{\theta}_x^H, \label{Defn: theta^L, theta^H estimators} \\
	&\hat{\gamma}^L = g^L(\hat{\theta}^L, \hat{\theta}^H, \hat{\eta}), &&\hat{\gamma}^H = g^H(\hat{\theta}^L, \hat{\theta}^H, \hat{\eta}) \label{Defn: gamma^L, gamma^H estimators}
\end{align}

The optimization problems in $\theta^L(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$ and $\theta^H(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$ are especially straightforward when treatment is exogenous. Recall the claim of equation \eqref{Display: strong duality, with smaller feasible set}: the supremum of $P_{1 \mid x}(\varphi) + P_{0 \mid x}(\psi)$ over the larger set $\Phi_c$ is the same value when restricted to $\Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)$. The argument behind this claim uses monotonicity of the maps $P_{d \mid x}$. When treatment is exogenous, $\hat{P}_{d \mid x}$ corresponds to a probability distribution and is therefore also monotonic. Thus the claim holds replacing $P_{d \mid x}$ with $\hat{P}_{d \mid x}$, implying the function classes $\mathcal{F}_c$ and $\mathcal{F}_c^c$ can be ignored in computation:
\begin{align}
	OT_c(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) &= \sup_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} \hat{P}_{1 \mid x}(\varphi) + \hat{P}_{0 \mid x}(\psi) = \sup_{(\varphi,\psi) \in \Phi_c} \hat{P}_{1 \mid x}(\varphi) + \hat{P}_{0 \mid x}(\psi) \notag \\
	&= \sup_{\{\varphi_i, \psi_j\}_{i,j}} \sum_{i=1}^n \omega_{1, x,i} \varphi_i + \sum_{j=1}^n \omega_{0, x, j} \psi_j \label{Display: computation with exogenous treatment} \\
	&\hspace{1 cm} \text{s.t. } \varphi_i + \psi_j \leq c(Y_i, Y_j) \text{ for all } 1 \leq i, j \leq n \notag
\end{align}
the final problem in this display is a linear programming problem with $2n$ choice variables and $n^2$ constraints, and can be further simplified by removing choice variables (and the corresponding constraints) whose weights $\omega_{d,x,i}$ equal zero. Many weights do equal zero, as only observations with $X_i = x$ correspond to nonzero weights.

When there is noncompliance in the sample, $\hat{P}_{d \mid x}$ does not correspond to a probability distribution. This is easily seen by noting that for observations $i$ where $Z_i$ differs from $D_i$, the weight $\omega_{d,x,i}$ defined in \eqref{Display: weights to compute Pdx(f)} is negative. Nonetheless, it remains computationally tractable to search over $\Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)$. For example, when the cost function is continuous $OT_c(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x})$ remains a linear programming problem, with additional linear constraints enforcing $\lvert \varphi_i + \psi_j \rvert \leq L\lvert Y_i - Y_j\rvert$, $-\lVert c \rVert_\infty \leq \varphi_i \leq \lVert c \rVert_\infty$, and $-2\lVert c \rVert_\infty \leq \psi_j \leq 0$.

\subsection{Asymptotic analysis}

The estimators proposed above are especially attractive because they are a (Hadamard directionally) differentiable map of the empirical distribution. Specifically, there exists a collection of functions $\mathcal{F}$ and a map $T : \ell^\infty(\mathcal{F}) \rightarrow \mathbb{R}^2$ described by equations \eqref{Display: T_1, map to conditional distributions s_x, main text}, \eqref{Display: T_1, map to conditional distributions, main text}, \eqref{Defn: theta_x^L, theta_x^H estimators}, \eqref{Defn: theta^L, theta^H estimators}, and \eqref{Defn: gamma^L, gamma^H estimators} such that
\begin{align*}
	&(\hat{\gamma}^L, \hat{\gamma}^H) = T(\mathbb{P}_n), &&(\gamma^L, \gamma^H) = T(P)
\end{align*}
The set $\mathcal{F}$ consists of the functions in $\mathcal{F}_c$, $\mathcal{F}_c^c$, and the coordinate functions defining $\eta$, multiplied by various indicator functions. It is formally defined in appendix \ref{Appendix: weak convergence}. Under assumption \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments}, $\mathcal{F}$ is a Donsker set and $T(\cdot)$ is continuous at $P$, which implies the esimators are consistent:
\begin{equation}
	(\hat{\gamma}^L, \hat{\gamma}^H) = T(\mathbb{P}_n) \overset{p}{\rightarrow} T(P) = (\gamma^L, \gamma^H) \label{Display: consistency, main text}
\end{equation}


\subsubsection{Weak convergence}
\label{Section: estimators, subsection weak convergence}

The map $T(\cdot)$ is not only continuous under assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments}, but Hadamard directionally differentiable. An application of the functional delta method gives the conclusion $\sqrt{n}((\hat{\gamma}^L, \hat{\gamma}^H) - (\gamma^L, \gamma^H))$ converges in distribution, a result stated formally in theorem \ref{Theorem: weak convergence, weak convergence of estimators} below.

In order to build hypothesis tests or construct confidence intervals based on the asymptotic distribution of $\sqrt{n}((\hat{\gamma}^L, \hat{\gamma}^H) - (\gamma^L, \gamma^H))$, one must be able to estimate the asymptotic distribution. This is possible under assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments}, but involves a more complex procedure described in section \ref{Section: estimators, subsection inference, subsubsection consistent alternative}. Under an additional assumption, a straightforward bootstrap will do.

For each instance of the restricted dual problem used in defining $T(\cdot)$, the set of maximizers
\begin{equation}
	\Psi_c(P_{1 \mid x}, P_{0 \mid x}) = \operatorname*{arg\,max}_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} P_{1 \mid x}(\varphi) + P_{0 \mid x}(\psi)
\end{equation}
is nonempty. If the solutions are suitably unique for each instance, the map $T(\cdot)$ is fully Hadamard differentiable at $P$ and a straightforward bootstrap will consistently estimate the asymptotic distribution.

Assumption \ref{Assumption: full differentiability} states this high-level uniqueness condition, while the following lemma \ref{Lemma: weak convergence, simple sufficient conditions for full differentiability} gives low-level sufficient conditions for it to hold. Let $\mathcal{Y}_{d,x}$ be the support of $Y$ conditional on $D = d$ and $X = x$, and $\mathbbm{1}_{\mathcal{Y}_{d,x}}(y) = \mathbbm{1}\{y \in \mathcal{Y}_{d,x}\}$ be the indicator function for this set.
\begin{restatable}[]{assumption}{assumptionFullDifferentiability}
	\label{Assumption: full differentiability}
	\singlespacing

	For each $x \in \mathcal{X}$, each $c \in \{c_L, c_H\}$, and any $(\varphi_1,\psi_1), (\varphi_2,\psi_2) \in \Psi_c(P_{1 \mid x},P_{0 \mid x})$, there exists $s \in \mathbb{R}$ such that
	\begin{align*}
		&\mathbbm{1}_{\mathcal{Y}_{1,x}} \times \varphi_1 = \mathbbm{1}_{\mathcal{Y}_{1,x}} \times (\varphi_2 + s), \; P\text{-a.s.} &&\text{ and } &&\mathbbm{1}_{\mathcal{Y}_{0,x}} \times \psi_1 = \mathbbm{1}_{\mathcal{Y}_{0,x}} \times (\psi_2 - s), \; P\text{-a.s.}
	\end{align*}

\end{restatable}


\begin{restatable}[]{lemma}{lemmaSimpleConditionsForFullDifferentiability}
	\label{Lemma: weak convergence, simple sufficient conditions for full differentiability}
	\singlespacing

	Suppose that
	\begin{enumerate}[label=(\roman*)]
		\item assumption \ref{Assumption: cost function} \ref{Assumption: cost function, smooth costs} holds, with cost function $c(y_1,y_0)$ that is continuously differentiable, and \label{Assumption: full differentiability lemma, cost is continuously differentiable}
		\item for each $(d,x)$, the support of $P_{d \mid x}$ is $\mathcal{Y}_{d, x}$, which is a bounded interval. \label{Assumption: full differentiability lemma, support of compliers}
	\end{enumerate}
	Then assumption \ref{Assumption: full differentiability} holds.
\end{restatable}

When treatment is exogenous, condition \ref{Assumption: full differentiability lemma, support of compliers} of lemma \ref{Lemma: weak convergence, simple sufficient conditions for full differentiability} simplifies to the assumption that the distribution of $Y_d \mid X = x$ has bounded support $[y_{d,x}^\ell, y_{d,x}^u]$. In general, this condition requires the support of $Y_d$ for the subpopulation of compliers with covariate value $x$ is a bounded interval that contains the support of the relevant subpopulation of non-compliers. Specifically, the support of $Y_1$ for compliers is a bounded interval containing the support of $Y_1$ for always-takers, and the support of $Y_0$ for compliers is a bounded interval containing the support of $Y_0$ for never-takers.

Assumption \ref{Assumption: full differentiability} can hold even when the conditions of lemma \ref{Lemma: weak convergence, simple sufficient conditions for full differentiability} do not. For example, when interest is in the cumulative distribution function and assumption \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} is satisfied, the dual problem is essentially optimizing over the difference of CDFs (see remark \ref{Remark: Makarov bounds}). Although the cost functions are not continuously differentiable, it is still plausible for this optimization problem to have a unique solution in well-behaved cases. For further discussion of uniqueness of the dual solutions of optimal transport, see \cite{staudt2022uniqueness}.

The following theorem gives the main weak convergence result.
\begin{restatable}[]{theorem}{theoremWeakConvergenceOfEstimators}
	\label{Theorem: weak convergence, weak convergence of estimators}
	\singlespacing

	Suppose assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments} hold, and let $\mathbb{G}$ be the weak limit of $\sqrt{n}(\mathbb{P}_n - P)$ in $\ell^\infty(\mathcal{F})$. Then $T$ is Hadamard directionally differentiable at $P$ tangentially to the support of $\mathbb{G}$, and
	\begin{align*}
		\sqrt{n}((\hat{\gamma}^L, \hat{\gamma}^H) - (\gamma^L, \gamma^H)) = \sqrt{n}(T(\mathbb{P}_n) - T(P)) \overset{L}{\rightarrow} T_P'(\mathbb{G})
	\end{align*}
	If assumption \ref{Assumption: full differentiability} also holds, then $T_P'$ is linear on the support of $\mathbb{G}$ and $T_P'(\mathbb{G})$ is bivariate normal.
\end{restatable}

\subsection{Inference}
\label{Section: estimators, subsection inference}

To make use of the weak convergence result of theorem \ref{Theorem: weak convergence, weak convergence of estimators} for inference, this section develops methods of estimating the law of $T_P'(\mathbb{G})$ by utilizing the bootstrap. The ``exchangeable bootstrap'' procedures discussed in \cite{vaart1997weak} are computationally convenient for reasons discussed below. These procedures define a new map $\mathbb{P}_n^* \in \ell^\infty(\mathcal{F})$ pointwise with
\begin{equation}
	\mathbb{P}_n^*(f) = \frac{1}{n}\sum_{i=1}^n W_i f(Y_i, D_i, Z_i, X_i) \label{Defn: exchangeable bootstrap}
\end{equation}
for nonnegative random variables $\{W_i\}_{i=1}^n$ independent of the data $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$, and satisfying technical conditions omitted here. I focus on two notable examples, the nonparametric bootstrap of \cite{efron1979bootstrap} and the ``Bayesian'' bootstrap of \cite{rubin1981bayesian}. Either bootstrap can be used to estimate the asymptotic distribution. The Bayesian bootstrap may be preferable in small samples for reasons discussed below.
\begin{restatable}[Nonparametric bootstrap]{definition}{definitionNonparametricBootstrap}
	\label{Definition: exchangeable bootstrap, nonparametric bootstrap}
	\singlespacing

	Let $(W_1, \ldots, W_n) \sim \text{Multinomial}(n, (1/n, \ldots, 1/n))$ be independent of the data $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$. Define $\mathbb{P}_n^* \in \ell^\infty(\mathcal{F})$ pointwise with \eqref{Defn: exchangeable bootstrap}.
\end{restatable}

\begin{restatable}[Bayesian bootstrap]{definition}{definitionBayesianBootstrap}
	\label{Definition: exchangeable bootstrap, bayesian bootstrap}
	\singlespacing

	Let $\{\xi_i\}_{i=1}^n$ be i.i.d. exponentially distributed random variables with mean $1$, independent of the data $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$. Set $W_i = \xi_i/(n^{-1}\sum_{i=1}^n \xi_i)$, and define $\mathbb{P}_n^* \in \ell^\infty(\mathcal{F})$ pointwise with \eqref{Defn: exchangeable bootstrap}.
\end{restatable}


The map $\mathbb{P}_n^*$ in \eqref{Defn: exchangeable bootstrap} can be used to compute $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*}) = T(\mathbb{P}_n^*)$ in much the same way that $T(\mathbb{P}_n)$ is computed. Specifically, bootstrap analogues of $\hat{p}_{d,x,z}$, $\hat{p}_{x,z}$, and $\hat{p}_x$ are given by
\begin{align*}
	&\hat{p}_{d,x,z}^* = \frac{1}{n}\sum_{i=1}^n W_i \mathbbm{1}_{d,x,z}(D_i, X_i, Z_i), &&\hat{p}_{x,z}^* = \frac{1}{n}\sum_{i=1}^n W_i \mathbbm{1}_{x,z}(X_i, Z_i), &&\hat{p}_x^* = \frac{1}{n}\sum_{i=1}^n W_i \mathbbm{1}_x(X_i),
\end{align*}
and the bootstrap analogue of $\hat{s}_x$ is
\begin{equation*}
	\hat{s}_x^* = \frac{(\hat{p}_{1,x,1}^*/\hat{p}_{x,1}^* - \hat{p}_{1,x,0}^*/\hat{p}_{x,0}^*)\hat{p}_x^*}{\sum_{x'} (\hat{p}_{1,x',1}^*/\hat{p}_{x',1}^* - \hat{p}_{1,x',0}^*/\hat{p}_{x',0}^*)\hat{p}_{x'}^*}
\end{equation*}
The maps $\hat{P}_{d \mid x}$ have bootstrap analogues
\begin{equation*}
	\hat{P}_{d \mid x}^*(f) = \frac{\mathbb{P}_n^*(\mathbbm{1}_{d,x,d} \times f)/\hat{p}_{x,d}^* - \mathbb{P}_n^*(\mathbbm{1}_{d,x,1-d} \times f)/\hat{p}_{x,1-d}^*}{\hat{p}_{d,x,d}^*/\hat{p}_{x,d}^* - \hat{p}_{d,x,1-d}^*/\hat{p}_{x,1-d}^*} = \sum_{i=1}^n \omega_{d,x,i}^* f_i
\end{equation*}
where $f_i = f(Y_i)$ and $\omega_{d,x,i}^*$ are bootstrap versions of the weights in \eqref{Display: weights to compute Pdx(f)}:
\begin{equation}
	\omega_{d,x,i}^* = \frac{W_i}{n} \times \frac{\mathbbm{1}_{d,x,d}(D_i, X_i, Z_i)/\hat{p}_{x,d}^* - \mathbbm{1}_{d,x,1-d}(D_i, X_i, Z_i)/\hat{p}_{x,1-d}^*}{\hat{p}_{d,x,d}^*/\hat{p}_{x,d}^* - \hat{p}_{d,x,1-d}^*/\hat{p}_{x,1-d}^*} \label{Display: weights to compute Pdx*(f)}
\end{equation}
Finally, $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*})$ can be computed with
\begin{align}
	&\hat{\theta}_x^{L*} = \theta^L(\hat{P}_{1 \mid x}^*, \hat{P}_{0 \mid x}^*), &&\hat{\theta}_x^{H*} = \theta^H(\hat{P}_{1 \mid x}^*, \hat{P}_{0 \mid x}^*), \label{Defn: theta_x^L, theta_x^H boostrap} \\
	&\hat{\theta}^{L*} = \sum_x \hat{s}_x^* \hat{\theta}_x^{L*}, &&\hat{\theta}^{H*} = \sum_x \hat{s}_x^* \hat{\theta}_x^{H*}, \label{Defn: theta^L, theta^H boostrap} \\
	&\hat{\gamma}^{L*} = g^L(\hat{\theta}^{L*}, \hat{\theta}^{H*}, \hat{\eta}^*), &&\hat{\gamma}^{H*} = g^H(\hat{\theta}^{L*}, \hat{\theta}^{H*}, \hat{\eta}^*) \label{Defn: gamma^L, gamma^H boostrap}
\end{align}

\subsubsection{Simple bootstrap with full differentiability}
\label{Section: estimators, subsection inference, subsubsection simple bootstrap}

Under assumption \ref{Assumption: full differentiability}, estimating the distribution of $T_P'(\mathbb{G})$ is straightforward.

\begin{restatable}[]{theorem}{theoremBootstrapWorksWithFullDifferentiability}
	\label{Theorem: inference, bootstrap, bootstrap works with full differentiability}
	\singlespacing

	Suppose assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, \ref{Assumption: parameter, function of moments}, and \ref{Assumption: full differentiability} hold, and let $\mathbb{P}_n^*$ be given by definition \ref{Definition: exchangeable bootstrap, nonparametric bootstrap} or \ref{Definition: exchangeable bootstrap, bayesian bootstrap}. Then conditional on $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$,
	\begin{equation*}
		\sqrt{n}(T(\mathbb{P}_n^*) - T(\mathbb{P}_n)) \overset{L}{\rightarrow} T_P'(\mathbb{G})
	\end{equation*}
	in outer probability.
\end{restatable}

It is worth emphasizing the computationally convenience of the bootstrap $\mathbb{P}_n^*$ given in \eqref{Defn: exchangeable bootstrap} when treatment is exogenous. The weights given in display \eqref{Display: weights to compute Pdx*(f)} simplify to
\begin{equation}
	\omega_{d,x,i}^* = \frac{W_i}{n} \times \frac{\mathbbm{1}\{D_i = d, X_i = x\}}{\hat{p}_{x,d}^*} \label{Display: Pdx for bootstrap written as weighted sum}
\end{equation}
As these weights are nonnegative and sum to one, $\hat{P}_{d \mid x}^*$ is a probability distribution. Accordingly, $\theta^L(\hat{P}_{1 \mid x}^*, \hat{P}_{0 \mid x}^*)$ and $\theta^H(\hat{P}_{1 \mid x}^*, \hat{P}_{0 \mid x}^*)$ can be computed ignoring the function classes $\mathcal{F}_c$ and $\mathcal{F}_c^c$ for the same reasons discussed around display \eqref{Display: computation with exogenous treatment}:
\begin{align*}
	OT_c(\hat{P}_{1 \mid x}^*, \hat{P}_{0 \mid x}^*) &= \sup_{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c)} \hat{P}_{1 \mid x}^*(\varphi) + \hat{P}_{0 \mid x}^*(\psi) = \sup_{(\varphi,\psi) \in \Phi_c} \hat{P}_{1 \mid x}^*(\varphi) + \hat{P}_{0 \mid x}^*(\psi) \notag \\
	&= \sup_{\{\varphi_i, \psi_j\}_{i,j}} \sum_{i=1}^n \omega_{1, x,i}^* \varphi_i + \sum_{j=1}^n \omega_{0, x, j}^* \psi_j \\
	&\hspace{1 cm} \text{s.t. } \varphi_i + \psi_j \leq c(Y_i, Y_j) \text{ for all } 1 \leq i, j \leq n \notag
\end{align*}

A researcher utilizing the nonparametric bootstrap runs the risk of a boostrap draw including no observations with $\mathbbm{1}\{D_i = d, X_i = x\}$. As $\hat{p}_{x,d}^* = \frac{1}{n}\sum_{i=1}^n W_i \mathbbm{1}\{D_i = d, X_i = x\}$, this would result in the formula in \eqref{Display: Pdx for bootstrap written as weighted sum} attempting to divide by zero. This problem cannot arise when using the Bayesian bootstrap suggested in \ref{Definition: exchangeable bootstrap, bayesian bootstrap}; in this procedure $W_i > 0$ for each $i$, and thus $\hat{p}_{x,d}^* = \frac{1}{n}\sum_{i=1}^n W_i \mathbbm{1}\{D_i = d, X_i = x\} > 0$ as long as $\hat{p}_{d,x} > 0$.


\subsubsection{Alternative for directional differentiability}
\label{Section: estimators, subsection inference, subsubsection consistent alternative}

The solutions to optimal transport may not be unique as assumption \ref{Assumption: full differentiability} requires. As emphasized in the statement of theorem \ref{Theorem: weak convergence, weak convergence of estimators}, assumption \ref{Assumption: full differentiability} is not needed to obtain the asymptotic distribution of the estimators. However, without assumption \ref{Assumption: full differentiability} the procedure suggested by lemma \ref{Theorem: inference, bootstrap, bootstrap works with full differentiability} may not consistently estimate that limiting distribution. When in doubt, researchers can make use of an alternative procedure based on the results of \cite{fang2019inference} and described below.

Additional notation is needed to describe this alternative. Let $\eta_{d, x}^{(k)} = P_{d \mid x}(\eta_d^{(k)})$, and let $T_1(\cdot)$ denote the  ``first stage'' function computing $P_{1 \mid x}$, $P_{0 \mid x}$, $\eta_{1,x}$, $\eta_{0, x}$, and $s_x$ for each $x$:
\begin{align*}
	T_1(P) = \left(\left\{P_{1 \mid x}, P_{0 \mid x}, \eta_{1,x}, \eta_{0, x}, s_x\right\}_{x \in \mathcal{X}}\right)
\end{align*}
Here $\{a_x\}_{x \in \mathcal{X}} = (a_{x_1}, \ldots, a_{x_M})$. Let $\{\kappa_n\}_{n=1}^\infty$ be a sequence in $\mathbb{R}$ satisfying $\kappa_n \uparrow \infty$ and $\kappa_n / \sqrt{n} \rightarrow 0$. Define the set of empirical approximate maximizers:
\begin{equation*}
	\widehat{\Psi}_{c, x} = \left\{(\varphi,\psi) \in \Phi_c \cap (\mathcal{F}_c \times \mathcal{F}_c^c) \; ; \; OT_c(\hat{P}_{1 \mid x}, \hat{P}_{0 \mid x}) \leq \hat{P}_{1 \mid x}(\varphi) + \hat{P}_{0 \mid x}(\psi) + \frac{\kappa_n}{\sqrt{n}} \right\}
\end{equation*}
and the maps
\begin{equation*}
	\widehat{OT}_{c, x}'(H_1, H_0) = \sup_{(\varphi,\psi) \in \widehat{\Psi}_{c,x}} H_1(\varphi) + H_0(\psi),
\end{equation*}
and
\begin{align*}
	&\widehat{T}_{2,T_1(P)}'\left(\{H_{1,x}, H_{0,x}, h_{\eta_1, x}, h_{\eta_0, x}, h_{s, x}\}_{x \in \mathcal{X}}\right) \\
	&\hspace{1 cm} = \left(\left\{\widehat{OT}_{c_L, x}'(H_{1,x}, H_{0,x}), -\widehat{OT}_{c_H, x}'(H_{1,x}, H_{0,x}), h_{\eta_1, x}, h_{\eta_0, x}, h_{s, x}\right\}_{x \in \mathcal{X}}\right)
\end{align*}

The alternative procedure uses the conditional law of
\begin{equation*}
	\hat{D}_4 \hat{D}_3 \widehat{T}_{2, T_1(P)}'\left(\sqrt{n}(T_1(\mathbb{P}_n^*) - T_1(\mathbb{P}_n))\right)
\end{equation*}
given the data, where $\hat{D}_4$ and $\hat{D}_3$ are matrices given by
\begin{align*}
	&\hat{D}_3 = \underset{(2+d_\eta) \times M(3 + d_\eta)}{
		\begin{bmatrix}
			\hat{D}_{3, x_1} & \hat{D}_{s, x_2} & \ldots & \hat{D}_{s, x_M}
	\end{bmatrix}},
	&&\hat{D}_{3,x} = \underset{(2+d_\eta) \times (3 + d_\eta)}{
		\begin{bmatrix}
			\hat{s}_x & 0 & 0 & 0 & \hat{\theta}_x^L \\
			0 & \hat{s}_x & 0 & 0 & \hat{\theta}_x^H \\
			0 & 0 & \hat{s}_x I_{K_1} & 0 & \hat{\eta}_{1,x} \\
			0 & 0 & 0  & \hat{s}_x I_{K_0}& \hat{\eta}_{0,x} \\
	\end{bmatrix}}, \\
	&D_4 =
	\underset{2 \times (2+d_\eta)}{
		\begin{bmatrix}
			\nabla g^L(\hat{\theta}^L, \hat{\theta}^H, \hat{\eta})^\intercal \\
			\nabla g^H(\hat{\theta}^L, \hat{\theta}^H, \hat{\eta})^\intercal \\
	\end{bmatrix}},
\end{align*}

\begin{restatable}[]{theorem}{theoremFangAndSantosAlternativeWorks}
	\label{Theorem: inference, bootstrap, Fang and Santos alternative works}
	\singlespacing

	Suppose assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments} hold, let $\mathbb{P}_n^*$ be given by definition \ref{Definition: exchangeable bootstrap, nonparametric bootstrap} or \ref{Definition: exchangeable bootstrap, bayesian bootstrap}, and $\{\kappa_n\}_{n=1}^\infty \subseteq \mathbb{R}$ satisfy $\kappa_n \rightarrow \infty$ and $\kappa_n / \sqrt{n} \rightarrow 0$. Then conditional on $\{Y_i, D_i, Z_i, X_i\}_{i=1}^n$,
	\begin{equation*}
		\hat{D}_4 \hat{D}_3 \widehat{T}_{2,T_1(P)}(\sqrt{n}(T_1(\mathbb{P}_n^*) - T_1(\mathbb{P}_n))) \overset{L}{\rightarrow} T_P'(\mathbb{G})
	\end{equation*}
	in outer probability.
\end{restatable}

\subsubsection{Confidence sets}
\label{Section: estimators, subsection confidence sets}

Theorems \ref{Theorem: inference, bootstrap, bootstrap works with full differentiability} and \ref{Theorem: inference, bootstrap, Fang and Santos alternative works} make it straightforward to conduct inference. For example, a simple confidence set for the identified set $[\gamma^L, \gamma^H]$ is given by
\begin{equation*}
	\left[\hat{\gamma}^L - \hat{c}_{1-\alpha}/\sqrt{n}, \hat{\gamma}^H + \hat{c}_{1-\alpha}/\sqrt{n}\right]
\end{equation*}
where $\hat{c}_{1-\alpha}$ is a consistent estimator of the $1-\alpha$ quantile of $\max\{T_P'(\mathbb{G})^{(1)}, -T_P'(\mathbb{G})^{(2)}\}$. When assumptions \ref{Assumption: setting} through \ref{Assumption: full differentiability} hold, let $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*}) = T(\mathbb{P}_n^*)$. When assumptions \ref{Assumption: setting} through \ref{Assumption: parameter, function of moments} hold but assumption \ref{Assumption: full differentiability} is doubtful, let $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*}) = (\hat{\gamma}^L, \hat{\gamma}^H) + \frac{1}{\sqrt{n}}\hat{D}_4 \hat{D}_3 \widehat{T}_{2,T_1(P)}(\sqrt{n}(T_1(\mathbb{P}_n^*) - T_1(\mathbb{P}_n)))$. In either case, compute
\begin{equation*}
	\hat{c}_{1-\alpha} = \inf\left\{c \; ; \; P\left(\max\left\{\sqrt{n}(\hat{\gamma}^{L*} - \hat{\gamma}^L), -\sqrt{n}(\hat{\gamma}^{H*} - \hat{\gamma}^H)\right\} \leq c \mid \{Y_i, D_i, Z_i, X_i\}_{i=1}^n\right) \geq 1-\alpha\right\}
\end{equation*}
through simulation:
\begin{enumerate}
	\item Compute $(\hat{\gamma}^L, \hat{\gamma}^H) = T(\mathbb{P}_n)$ and, if necessary, $\hat{D}_4$, and $\hat{D}_3$.
	\item Generate $N$ boostrap samples, $\{W_{i,b}\}_{i=1}^n$ for each $b = 1, \ldots, N$ according to definition \ref{Definition: exchangeable bootstrap, nonparametric bootstrap} or \ref{Definition: exchangeable bootstrap, bayesian bootstrap}. For each bootstrap sample $b$, compute $(\hat{\gamma}_b^{L*}, \hat{\gamma}_b^{H*})$ as described above.
	\item Let $\hat{c}_{1-\alpha}$ be the $1-\alpha$ quantile of $\{\max\{\sqrt{n}(\hat{\gamma}_b^{L*} - \hat{\gamma}^L), -\sqrt{n}(\hat{\gamma}_b^{H*} - \hat{\gamma}^H)\}_{b=1}^N$.
\end{enumerate}
Under the further assumption that the cumulative distribution function of $\max\{T_P'(\mathbb{G})^{(1)}, -T_P'(\mathbb{G})^{(2)}\}$ is continuous and strictly increasing at its $1-\alpha$ quantile,
\begin{equation*}
	\lim_{n \rightarrow \infty} P\left([\gamma^L, \gamma^H] \subseteq \left[\hat{\gamma}^L - \hat{c}_{1-\alpha}/\sqrt{n}, \hat{\gamma}^H + \hat{c}_{1-\alpha}/\sqrt{n}\right] \right) = 1-\alpha
\end{equation*}
\begin{comment}
	The derivation: let $M : \mathbb{R}^2 \rightarrow \mathbb{R}$ be given by $M(x_1, x_2) = \max\{x_1,-x_2\}$. Then
	\begin{equation*}
		\hat{c}_{1-\alpha} = \inf\left\{c \; ; \; P\left(M(\sqrt{n}(T(\mathbb{P}_n^*) - T(\mathbb{P}_n)))\leq c \mid \{Y_i, D_i, Z_i, X_i\}_{i=1}^n\right) \geq 1-\alpha\right\}
	\end{equation*}
	The bootstrap ``works'' here because $M$ is Lipschitz continuous; thus if $f \in \text{BL}_1$, then $f\circ M$ is Lipschitz and bounded; see \cite{kosorok2008introduction} proposition 10.7. Then notice that
	\begin{align*}
		&P\left([\gamma^L, \gamma^H] \subseteq \left[\hat{\gamma}^L - \hat{c}_{1-\alpha}/\sqrt{n}, \hat{\gamma}^H + \hat{c}_{1-\alpha}/\sqrt{n}\right] \right) \\
		&\hspace{1 cm} = P\left(\hat{\gamma}^L - \hat{c}_{1-\alpha}/\sqrt{n} \leq \gamma^L, \gamma^H \leq \hat{\gamma}^H + \hat{c}_{1-\alpha}/\sqrt{n} \right) \\
		&\hspace{1 cm} = P\left(\sqrt{n}(\hat{\gamma}^L - \gamma^L) \leq \hat{c}_{1-\alpha}, -\sqrt{n}(\hat{\gamma}^H - \gamma^H) \leq \hat{c}_{1-\alpha}\right) \\
		&\hspace{1 cm} = P\left(M(\sqrt{n}(T(\mathbb{P}_n) - T(P))) \leq \hat{c}_{1-\alpha}\right) \\
		&\hspace{1 cm} \rightarrow 1-\alpha
	\end{align*}
	where the limit claim follows from $\hat{c}_{1-\alpha} \overset{p}{\rightarrow} c_{1-\alpha}$ (implied by the CDF being strictly increasing at $c_{1-\alpha}$, implying the quantile function is continuous at $c_{1-\alpha}$) and the CDF of $M(T_P'(\mathbb{G}))$ being continuous at $c_{1-\alpha}$ implies, through weak convergence of $M(\sqrt{n}(T(\mathbb{P}_n) - T(P))) - \hat{c}_{1-\alpha} \overset{L}{\rightarrow} M(T_P'(\mathbb{G})) - c_{1-\alpha}$, that $P\left(M(\sqrt{n}(T(\mathbb{P}_n) - T(P))) \leq \hat{c}_{1-\alpha}\right) \rightarrow P(M(T_P'(\mathbb{G})) \leq c_{1-\alpha}) = 1-\alpha$.
\end{comment}
Confidence sets for the parameter could be constructed following \cite{imbens2004confidence}.

	\section{Application: job training experiment}
\label{Section: application}

In this section I demonstrate the estimators in revisiting the famous National Supported Work Demonstration program (\cite{lalonde1986evaluating}). This program was implemented in the 1970s with the aim of helping socially and economically disadvantaged workers obtain job skills. Those randomly selected into the program were guaranteed a job lasting six to eighteen months, and frequently met with a counselor to discuss performance.

I make use of the ``LaLonde'' sample studied in \cite{diamond2013genetic}. This sample consists of male participants and includes 297 treated and 425 control observations. The outcome of interest is real earnings in 1978. Observed covariates include age, years of education, real earnings in months 13 to 24 prior to randomization, and indicators for whether a participant is a high school dropout, black, hispanic, or married. Averages and standard deviations of these covariates by treatment status are reported in table \ref{Table: balance table for NSW}:

\begin{table}[H]
	\begin{center}
				\caption{Balance table}
		\label{Table: balance table for NSW}
		\begin{tabular}{l||cccccccc}
			\toprule
			&       base inc. &   age &  yrs. educ &  HS dropout &  black &  hispanic &  married & $N$ \\ \hline \hline
			\multirow{2}{*}{control} & 3672.49 & 24.45 & 10.19 &      0.81 &   0.80 &  0.11 &     0.16 & \multirow{2}{*}{425} \\
			& (6521.53) &  (6.59) & (1.62) &      (0.39) &   (0.40) &  (0.32) &     (0.36) & \\ \hline
			\multirow{2}{*}{treated} & 3571.00 & 24.63 & 10.38 &      0.73 &   0.80 &  0.09 &     0.17 & \multirow{2}{*}{297} \\
			& (5773.13) &  (6.69) &  (1.82) &      (0.44) &   (0.40) &  (0.29) &     (0.37) & \\
			\bottomrule
		\end{tabular}
	\end{center}
	\vspace{-0.3 cm}
	\hspace{0.5 cm} {\scriptsize \textit{Note:} Standard deviations in parentheses.}
\end{table}

There is no reported noncompliance, so I interpret the setting as one of exogenous treatment. The parameter of interest is the OLS slope coefficient of regressing treatment effects on a constant and $Y_0$:
\begin{align*}
	\gamma = \frac{\text{Cov}(Y_1 - Y_0, Y_0)}{\text{Var}(Y_0)} = \frac{E_{P_{1,0}}[(Y_1 - Y_0) Y_0] - (E_{P_1}[Y_1] - E_{P_0}[Y_0])E_{P_0}[Y_0]}{E_{P_0}[Y_0^2] - (E_{P_0}[Y_0])^2}
\end{align*}
as described in example \ref{Example: equitable policies}, the sign of this parameter describes who receives larger benefits from treatment: $\gamma < 0$ implies those with below average untreated outcomes tend to see above average treatment effects.

Discretized versions of baseline income and age are found to be informative covariates. Baseline income is binned as: $[0,0]$ or $(0, \infty)$, while age is binned as $(16,20]$, $(20, 26]$, or $(26, \infty)$. $X$ is the cartesian product of bins. The resulting $(d,x)$ bins have a minimum of 31 observations per bin, and an average of 60.2 observations per bin.

The point estimates are $(\hat{\gamma}^L, \hat{\gamma}^H) = (-1.73, -0.004)$. The negative upper bound point estimates suggests that the treatment was especially beneficial for participants who would otherwise have incomes below average (for the eligible population). Covariates are found to be informative, especially for the upper bound. Ignoring covariates, the lower bound point estimate is $-1.78$ and the upper bound point estimate is $0.189$. The $95\%$ confidence set for the identified based on 500 bootstrap draws is $[-1.94, 0.20]$, suggesting $\gamma$ may still be zero or slightly positive once accounted for sample uncertainty.



	\section{Extensions}
\label{Section: extensions}

This section briefly describes simple extensions.

\subsection{Conditioning on $X \in A$}

In many applications parameters conditional on a covariate taking a particular value are of interest. For example, the share of compliers of a particular demographic benefiting from treatment is $P(Y_1 > Y_0 \mid D_1 > D_0, \text{demographic})$.

Such parameters can be written in the form
\begin{align*}
	\gamma_A = g(\theta_A, \eta_A)
\end{align*}
where for a known set $A \subseteq \mathcal{X}$,
\begin{align*}
	&\theta_A \equiv E[c(Y_1,Y_0) \mid D_1 > D_0, X \in A], &&\eta_A \equiv E[\eta_1(Y_1), \eta_0(Y_0) \mid D_1 > D_0, X \in A]
\end{align*}
The identified set for $\gamma_A$ is straightforward to characterize and estimate. First note that
\begin{align*}
	\theta_A = E[\theta_X \mid D_1 > D_0, X \in A] = \frac{1}{s_A} \sum_{x \in A} s_x \theta_x
\end{align*}
where $s_A = \sum_{x \in A} s_x$. The proof of theorem \ref{Theorem: identification, function of moments} shows that the sharp identified set for $(\theta_{x_1}, \ldots, \theta_{x_M})$ is in fact $[\theta_{x_1}^L, \theta_{x_1}^H] \times \ldots \times [\theta_{x_M}^L, \theta_{x_M}^H]$. It follows that the sharp identified set for $\theta_A$ is $[\theta_A^L, \theta_A^H]$, where
\begin{align*}
	&\theta_A^L = \frac{1}{s_A} \sum_{x \in A} s_x \theta_x^L, &&\theta_A^H = \frac{1}{s_A} \sum_{x \in A} s_x \theta_x^H
\end{align*}
and the sharp identified set for $\gamma_A$ is $[\gamma_A^L, \gamma_A^H]$ where
\begin{align*}
	&\gamma_A^L = \min_{t \in [\theta_A^L, \theta_A^H]} g(t, \eta_A), &&\gamma_A^H = \max_{t \in [\theta_A^L, \theta_A^H]} g(t, \eta_A),
\end{align*}

Let $\hat{s}_x$, $\hat{\theta}_x^L$, and $\hat{\theta}_x^H$ be as defined in section \ref{Section: estimators}. Let $\hat{s}_A = \sum_{x \in A} \hat{s}_x$ and
\begin{align*}
	&\hat{\theta}_A^L = \frac{1}{\hat{s}_A} \sum_{x \in A} \hat{s}_x \hat{\theta}_x^L, &&\hat{\theta}^H(A) = \frac{1}{\hat{s}_A} \sum_{x \in A} \hat{s}_x \hat{\theta}_x^H \\
	&\hat{\gamma}_A^L = \min_{t \in [\hat{\theta}_A^L, \hat{\theta}_A^H]} g(t, \hat{\eta}_A), &&\hat{\gamma}_A^H = \max_{t \in [\hat{\theta}_A^L, \hat{\theta}_A^H]} g(t, \hat{\eta}_A),
\end{align*}
Under assumptions \ref{Assumption: setting}, \ref{Assumption: cost function}, and \ref{Assumption: parameter, function of moments}, $\sqrt{n}((\hat{\gamma}_A^L, \hat{\gamma}_A^H) - (\gamma_A^L, \gamma_A^H)$ will converge weakly. With assumption \ref{Assumption: full differentiability} the straightforward bootstrap will consistently estimate its asymptotic distribution.

\subsection{Quantiles}
\label{Section: extensions, subsection quantiles}

Example \ref{Example: pseudo quantiles} considers the parameter $q_\tau$ solving
\begin{equation*}
	P(Y_1 - Y_0 \leq q_\tau \mid D_1 > D_0) = \tau
\end{equation*}
As noted in that example, the sharp identification results for $P(Y_1 - Y_0 \leq \delta \mid D_1 > D_0)$ can be adapted to characterize the sharp identified set for $q_\tau$. First view the bounds on the cumulative distribution function as functions of $\delta$:
\begin{align*}
	&c_{L,\delta}(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 < \delta\}, &&c_{H,\delta} (y_1,y_0) = \mathbbm{1}\{y_1 - y_0 > \delta\}, \\
	&\theta_x^L(\delta) = OT_{c_{L,\delta}}(P_{1 \mid x}, P_{0 \mid x}), &&\theta_x^H(\delta) = 1 - OT_{c_{H, \delta}}(P_{1 \mid x}, P_{0 \mid x}) \\
	&\theta^L(\delta) = \sum_x s_x \theta_x^L(\delta) &&\theta^H(\delta) = \sum_x s_x \theta_x^H(\delta)
\end{align*}
Let $Q_{I, \tau}$ denote the sharp identified set for $q_\tau$.

\begin{restatable}[Identification of $q_\tau$]{lemma}{lemmaIdentificationPseudoQuantile}
	\label{Lemma: identification, pseudo quantile}
	\singlespacing

	Suppose assumptions \ref{Assumption: setting} and \ref{Assumption: cost function} \ref{Assumption: cost function, CDF} hold. Then $q \in Q_{I, \tau}$ if and only if $\theta^L(q) \leq \tau \leq \theta^H(q)$.
\end{restatable}

Lemma \ref{Lemma: identification, pseudo quantile} implies that inverting a test of $H_0 : \theta^L(q) \leq \tau \leq \theta^H(q)$ against the alternative $H_1 : \tau < \theta^L(q) \text{ or } \theta^H(q) < \tau$ will lead to valid confidence sets for $q_\tau$.

\begin{remark}
	\label{Remark: pseudo quantile interpretation and definition}
	\singlespacing

	Consider instead defining $q_\tau$ to be the closed subset of $\mathbb{R}$ given by
	\begin{equation*}
		q_\tau = [\inf\{y \; ; \; P(Y_1 - Y_0 \leq y) \geq \tau\}, \inf\{y \; ; \; P(Y_1 - Y_0 \leq y) > \tau\}]
	\end{equation*}
	Note that this $q_\tau$ is the singleton $\inf\{y \; ; \; P(Y_1 - Y_0 \leq y) \geq \tau\}$, unless $P(Y_1 - Y_0 \leq \cdot)$ is flat when equal to $\tau$, in which case it equals the $\tau$-level set $\{y \; ; \; P(Y_1 - Y_0 \leq y) = \tau\}$. (Compare \cite{ehm2016quantiles}, who define the $\tau$-th quantile equivalently as $q_\tau = [\sup\{y \; ; \; P(Y_1 - Y_0 \leq y) < \tau\}, \sup\{y \; ; \; P(Y_1 - Y_0 \leq y) \leq \tau\}]$.) Let $Q_{I, \tau}$ denote the identified set of $q_\tau$ as defined in this remark. Lemma \ref{Lemma: identification, pseudo quantile, alternative definition} in appendix \ref{Appendix: identification} shows that under assumptions \ref{Assumption: setting} and \ref{Assumption: cost function} \ref{Assumption: cost function, CDF}, $q \in Q_{I, \tau}$ if and only if $\theta^L(q) \leq \tau \leq \theta^H(q)$.
\end{remark}

\subsection{Multiple treatment arms with exogenous treatment}

The identification results and estimators proposed above are easily extended to a setting with multiple treatment arms and exogenous treatment. Let the mutually exclusive treatment arms indexed by $d \in \{0, 1, \ldots, J\}$, with $d = 0$ indicating control. Let $Y_d$ be the potential outcome with treatment $d$, $D_d$ equal one if the unit has treatment $d$ and zero otherwise. The observed outcome is
\begin{align*}
	Y = \sum_{d=0}^J D_d Y_d
\end{align*}
Let $D = (D_0, D_1, \ldots, D_J)$ and assume
\begin{align*}
	(Y_0, Y_1, \ldots, Y_J) \perp D \mid X
\end{align*}
Note that the marginal distributions of $Y_d \mid X = x$, denoted $P_{d \mid x}$, are identified with the relation
\begin{align*}
	E_{P_{d \mid x}}[f(Y_d)] = E[f(Y_d) \mid X = x] = \frac{E[f(Y) D_d \mid X = x]}{P(D_d = 1 \mid X = x)}
\end{align*}

Let $\gamma_d = g(\theta_d, \eta_d)$ where $\theta_d = E[c(Y_d, Y_0)]$. Consider estimating the sharp identified set for $(\gamma_1, \ldots \gamma_J)$. For example, an RCT with two treatment arms may have similar average treatment effects. The treatment arms may be further distinguished by comparing $P(Y_1 - Y_0 > 0)$ with $P(Y_2 - Y_0 > 0)$, or $\text{Cov}(Y_1 - Y_0, Y_0)$ with $\text{Cov}(Y_2 - Y_0, Y_0)$.



Let $\theta_{d,x} = E[c(Y_1,Y_0) \mid X = x]$. The sharp identified set for $(\theta_{1,x}, \ldots, \theta_{J,x})$ is given by
\begin{align*}
	[\theta_{1,x}^L, \theta_{1,x}^H] \times \ldots \times [\theta_{J,x}^L, \theta_{J,x}^H]
\end{align*}
where $\theta_{d,x}^L = \theta^L(P_{d \mid x}, P_{0 \mid x})$ and $\theta_{d,x}^H = \theta^H(P_{d \mid x}, P_{0 \mid x})$ as in section \ref{Section: identification}.\footnote{This follows from existing results and the \textit{gluing lemma}, found in \cite{villani2009optimal} (pp. 11-12).}
The sharp identified set for $\theta_d$ is $[\theta_d^L, \theta_d^H]$ where $\theta_d^L = \sum_x s_x \theta_{d,x}^L$ and $\theta_d^H = \sum_x s_x \theta_{d,x}^H$, and the sharp identified set for $(\gamma_1, \ldots \gamma_J)$ is
\begin{align*}
	[\gamma_1^L, \gamma_1^H] \times \ldots \times [\gamma_J^L, \gamma_J^H]
\end{align*}
Sample analogues $(\hat{\gamma}_1^L, \hat{\gamma}_1^H, \ldots, \hat{\gamma}_J^L, \hat{\gamma}_J^H)$ can be formed just as in section \ref{Section: estimators}. Under natural adjustments to assumptions \ref{Assumption: cost function}, \ref{Assumption: parameter, function of moments}, and \ref{Assumption: full differentiability}, the same arguments work to show
\begin{align*}
	\sqrt{n}((\hat{\gamma}_1^L, \hat{\gamma}_1^H, \ldots, \hat{\gamma}_J^L, \hat{\gamma}_J^H) - (\gamma_1^L, \gamma_1^H, \ldots, \gamma_J^L, \gamma_J^H))
\end{align*}
is asymptotically Gaussian and the bootstrap consistently estimates its asymptotic distribution.



	\section{Conclusion}
\label{Section: conclusion}


This paper studies a large class of causal parameters that depend on a moment of the joint distribution of potential outcomes. The sharp identified set of such parameters is characterized with optimal transport. Estimators based on this identification are $\sqrt{n}$-consistent and converge in distribution under mild assumptions, and inference procedures based on the bootstrap are straightforward and computationally convenient.

	\newpage

	\bibliography{OTJointPO_bibliography}