Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
92,316 characters · 21 sections · 34 citation commands
Estimating Functionals of the Joint Distribution of Potential Outcomes with Optimal Transport
\pagenumbering{gobble}
\pagenumbering{arabic}
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” 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 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). As argued in 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 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. 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. 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 galichon2017survey for a recent survey. Several recent working papers utilize optimal transport for issues related to casual inference, including inverse propensity weighting dunipace2021optimal, matching on covariates gunsilius2021matching, and obtaining counterfactual distributions torous2021optimal. In concurrent and highly complementary work, 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 ji2023model.
The remainder of this paper is organized as follows. Section (ref) formalizes the setting and introduces the class of parameters under study. Optimal transport is introduced in section (ref), and used in identification in section (ref). Section (ref) proposes the estimators and contains the asymptotic results. Section (ref) 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) discusses straightforward extensions, and section (ref) concludes.
Consider a potential outcomes framework with binary treatment, a binary instrument, and finitely supported covariates (imbens1994late, 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
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
It is assumed that the instrument itself does not affect the outcome.\footnote{\linespread{1}\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 instrument exclusion, also known as the 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 compliers.
Assumption (ref) formalizes the setting.
Assumption (ref) is essentially equivalent to assumption 2.1 in abadie2003semiparametric, with the addition that covariates are finitely supported. Instrument independence is sometimes referred to as 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 (ref) 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.
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. abadie2003semiparametric shows that assumption (ref) suffices to identify the marginal distributions of $Y_1$ and $Y_0$ for the subpopulation of compliers.
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}$,
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 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). 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
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\}$.
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
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) places restrictions on the cost function to ensure optimal transport can be used characterize and estimate the sharp identified set for $\theta$.
Assumption (ref) covers every example listed below. Continuous cost functions $c$ are given a unified analysis, but for reasons discussed in section (ref) 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) (ref) requires the cdfs $F_{d \mid x}$ be continuous. As discussed in section (ref), 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) hold 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) and (ref), the sharp identified set for $\theta$ is an interval $[\theta^L, \theta^H]$. Assumption (ref) contains conditions on $g$ and $\eta$.
Note that when $\theta$ itself is of interest, assumption (ref) is satisfied with $g(\theta, \eta) = \theta$. Assumption (ref) (ref) ensures the identified set for $\gamma$ is the interval $[\gamma^L, \gamma^H]$, and assumption (ref) (ref) is used to apply the delta method. It is straightforward to show assumption (ref) (ref) 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) (ref) to hold.
The following examples are intended both to fix ideas and illustrate the broad scope of the parameter class described above.
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 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)]$:
This minimization problem in (ref) is referred to as the 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)$:
The dual problem chooses a pair of functions in $\Phi_c$ to maximize the sum of the corresponding expectations:
When the cost function is lower semicontinuous and bounded from below, the primal problem is attained and strong duality holds:
The dual problem will be used to construct and analyze estimators. Indeed, the identification of $P_{d \mid x}$ in lemma (ref) 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
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 $c$-transforms of a function $\varphi$:
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 $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
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
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
Equation (ref) 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
This problem will be referred to as the 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 (ref) when considering identification, and as the functional given in (ref) 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.
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:
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), 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) (ref) holds ($c(y_1,y_0)$ is Lipschitz continuous and $\mathcal{Y}$ is compact), define
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) (ref) holds (the cost function is $c(y_1,y_0) = \mathbbm{1}\{y_1 - y_0 \leq \delta\}$) define
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) hold 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) and (ref), 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
Under assumptions (ref), (ref), and (ref), the identified set for $\gamma$ is $[\gamma^L, \gamma^H]$, with endpoints
The following theorem summarizes the discussion above. Let $\theta^L(\cdot, \cdot)$ and $\theta^H(\cdot, \cdot)$ be given by (ref) or (ref) depending on the cost function, and set
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 (ref) and (ref) 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$,
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.
Sample analogues of the expressions identifying $P_{1 \mid x}$, $P_{0 \mid x}$, and $s_x$ in lemma (ref) 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:
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:
In this notation, $s_x = P(X = x \mid D_1 > D_0)$ and its empirical analogue $\hat{s}_x$ are
The maps $P_{d \mid x}$ and their empirical analogues are
Under assumption (ref), $\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:
where $f_i = f(Y_i)$ and the weights $\omega_{d,x,i}$ can be computed directly from data:
Sample analogue estimators of $\gamma^L$ and $\gamma^H$ are based on equations (ref), (ref), (ref), (ref), and (ref). 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) is written
Here $\mathcal{F}_c$, $\mathcal{F}_c^c$, and the functions $\theta^L(\cdot)$, $\theta^H(\cdot)$ are defined according to the cost function:
The sample analogue estimators are given by
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 (ref): 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:
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 (ref) 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$.
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 (ref), (ref), (ref), (ref), and (ref) such that
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). Under assumption (ref), (ref), and (ref), $\mathcal{F}$ is a Donsker set and $T(\cdot)$ is continuous at $P$, which implies the esimators are consistent:
The map $T(\cdot)$ is not only continuous under assumptions (ref), (ref), and (ref), 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) 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), (ref), and (ref), but involves a more complex procedure described in section (ref). 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
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) states this high-level uniqueness condition, while the following lemma (ref) 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.
When treatment is exogenous, condition (ref) of lemma (ref) 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) can hold even when the conditions of lemma (ref) do not. For example, when interest is in the cumulative distribution function and assumption (ref) (ref) is satisfied, the dual problem is essentially optimizing over the difference of CDFs (see remark (ref)). 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 staudt2022uniqueness.
The following theorem gives the main weak convergence result.
To make use of the weak convergence result of theorem (ref) 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 vaart1997weak are computationally convenient for reasons discussed below. These procedures define a new map $\mathbb{P}_n^* \in \ell^\infty(\mathcal{F})$ pointwise with
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 efron1979bootstrap and the “Bayesian” bootstrap of rubin1981bayesian. Either bootstrap can be used to estimate the asymptotic distribution. The Bayesian bootstrap may be preferable in small samples for reasons discussed below.
The map $\mathbb{P}_n^*$ in (ref) 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
and the bootstrap analogue of $\hat{s}_x$ is
The maps $\hat{P}_{d \mid x}$ have bootstrap analogues
where $f_i = f(Y_i)$ and $\omega_{d,x,i}^*$ are bootstrap versions of the weights in (ref):
Finally, $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*})$ can be computed with
Under assumption (ref), estimating the distribution of $T_P'(\mathbb{G})$ is straightforward.
It is worth emphasizing the computationally convenience of the bootstrap $\mathbb{P}_n^*$ given in (ref) when treatment is exogenous. The weights given in display (ref) simplify to
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 (ref):
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 (ref) attempting to divide by zero. This problem cannot arise when using the Bayesian bootstrap suggested in (ref); 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$.
The solutions to optimal transport may not be unique as assumption (ref) requires. As emphasized in the statement of theorem (ref), assumption (ref) is not needed to obtain the asymptotic distribution of the estimators. However, without assumption (ref) the procedure suggested by lemma (ref) may not consistently estimate that limiting distribution. When in doubt, researchers can make use of an alternative procedure based on the results of 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$:
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:
and the maps
and
The alternative procedure uses the conditional law of
given the data, where $\hat{D}_4$ and $\hat{D}_3$ are matrices given by
Theorems (ref) and (ref) make it straightforward to conduct inference. For example, a simple confidence set for the identified set $[\gamma^L, \gamma^H]$ is given by
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) through (ref) hold, let $(\hat{\gamma}^{L*}, \hat{\gamma}^{H*}) = T(\mathbb{P}_n^*)$. When assumptions (ref) through (ref) hold but assumption (ref) 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
through simulation:
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,
Confidence sets for the parameter could be constructed following imbens2004confidence.
In this section I demonstrate the estimators in revisiting the famous National Supported Work Demonstration program (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 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):
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$:
as described in example (ref), 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.
This section briefly describes simple extensions.
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
where for a known set $A \subseteq \mathcal{X}$,
The identified set for $\gamma_A$ is straightforward to characterize and estimate. First note that
where $s_A = \sum_{x \in A} s_x$. The proof of theorem (ref) 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
and the sharp identified set for $\gamma_A$ is $[\gamma_A^L, \gamma_A^H]$ where
Let $\hat{s}_x$, $\hat{\theta}_x^L$, and $\hat{\theta}_x^H$ be as defined in section (ref). Let $\hat{s}_A = \sum_{x \in A} \hat{s}_x$ and
Under assumptions (ref), (ref), and (ref), $\sqrt{n}((\hat{\gamma}_A^L, \hat{\gamma}_A^H) - (\gamma_A^L, \gamma_A^H)$ will converge weakly. With assumption (ref) the straightforward bootstrap will consistently estimate its asymptotic distribution.
Example (ref) considers the parameter $q_\tau$ solving
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$:
Let $Q_{I, \tau}$ denote the sharp identified set for $q_\tau$.
Lemma (ref) 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$.
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
Let $D = (D_0, D_1, \ldots, D_J)$ and assume
Note that the marginal distributions of $Y_d \mid X = x$, denoted $P_{d \mid x}$, are identified with the relation
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
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).\footnote{This follows from existing results and the gluing lemma, found in 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
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). Under natural adjustments to assumptions (ref), (ref), and (ref), the same arguments work to show
is asymptotically Gaussian and the bootstrap consistently estimates its asymptotic distribution.
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.