Nonparametric Point Identification of Treatment Effect Distributions via Rank Stickiness

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.

50,494 characters · 17 sections · 23 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractTreatment effect distributions are not identified without restrictions on the joint distribution of potential outcomes. Existing approaches either impose rank preservation---a strong assumption---or derive partial identification bounds that are often wide. We show that a single scalar parameter, rank stickiness, suffices for nonparametric point identification while permitting rank violations. The identified joint distribution---the coupling that maximizes average rank correlation subject to a relative entropy constraint, which we call the Bregman-Sinkhorn copula---is uniquely determined by the marginals and rank stickiness. Its conditional distribution is an exponential tilt of the marginal with a Bregman divergence as the exponent, yielding closed-form conditional moments and rank violation probabilities; the copula nests the comonotonic and Gaussian copulas as special cases. The empirical Bregman-Sinkhorn copula converges at the parametric $\sqrt{n}$-rate with a Gaussian process limit, despite the infinite-dimensional parameter space. We apply the framework to estimate the full treatment effect distribution, derive a variance estimator for the average treatment effect tighter than the Fr\'{e}chet--Hoeffding and Neyman bounds, and extend to observational studies under unconfoundedness.

\smallKeywords--- treatment effect distribution, identification, heterogeneity, imputation, nonparametric copula, sensitivity analysis.

Introduction

When the goal of a policy evaluation is to understand who benefits and who is harmed, the Average Treatment Effect (ATE) or Quantile Treatment Effects (QTE) are insufficient---what is needed is the entire Treatment Effect Distribution (TED). Consider a job training program: a positive ATE on earnings may obscure the fact that gains are concentrated among high-skilled workers while low-skilled workers experience adverse effects. The TED thus reveals the complete pattern of individual heterogeneity in program impacts.

However, the TED is not identified without additional assumptions: only one potential outcome is observed per unit, so the joint distribution $(Y_0, Y_1) \sim \pi$ is never observed. The existing literature faces a sharp trade-off between the strength of assumptions and informativeness of conclusions. At one extreme, manski1990nonparametric,manski1997monotone derive partial identification bounds under minimal assumptions, but these worst-case bounds are often wide blundell2007changes. Sharper bounds exploit the coupling structure: HeckmanSmithClements1997 and fan2010sharp employ Fr\'{e}chet--Hoeffding copula bounds; kim2014identifying tighten these under shape restrictions; firpo2019partial derive bounds on functionals of the joint distribution. At the other extreme, rank preservation---the requirement that each individual's rank in the outcome distribution is unchanged by treatment, $$ (Y_0 - Y_0') \cdot (Y_1 - Y_1') > 0 \quad \text{a.s.\ for independent draws } (Y_0, Y_1), (Y_0', Y_1') \sim \pi \;, $$ ---yields point identification: the TED reduces to the quantile-by-quantile difference of the two marginals firpo2007efficient; conditional variants have been studied by chernozhukov2005iv and athey2006identification. Yet HeckmanSmithClements1997 argue this assumption is implausible in many empirical settings; indeed, under existing approaches, even mild departures from rank preservation undermine nonparametric point identification of the TED.

This paper shows that point identification does not require rank preservation. We introduce a single structural parameter---rank stickiness $\rho$---that governs the degree of departure from comonotonicity and suffices for nonparametric point identification of the TED given the marginals. Rank stickiness captures the empirical regularity that treatment shifts the distribution of outcomes while individual ranks exhibit persistence, a phenomenon well-documented in earnings mobility both within individual lifetimes kopczuk2010earnings and across generations chetty2014united. We treat $\rho$ as a sensitivity parameter: rather than assuming a single value, we characterize how the identified TED varies with $\rho$.

Formally, we model the joint distribution of $(Y_0, Y_1) \sim \pi$ as the coupling that maximizes average rank correlation subject to a relative entropy constraint parameterized by $\rho \in (0,1]$: $$ \max_{\pi \in \Pi(\mu, \nu)~ \text{s.t.}~ \mathsf{KL}(\pi \mid \mu \otimes \nu) \leq R(\rho)} ~\E_{(Y_0, Y_1) \sim \pi} \E_{(Y_0', Y_1') \sim \pi} (Y_0 - Y_0') \cdot(Y_1 - Y_1')\;, $$ where $\Pi(\mu, \nu)$ is the set of couplings with marginals $\mu, \nu$, and $R(\rho)$ is a monotone function of $\rho$ that calibrates the degree of rank stickiness. We call the solution the Bregman-Sinkhorn copula. The conditional distribution of $Y_1 \mid Y_0$ takes the form of an exponential tilting of the treatment marginal, with exponent given by a Bregman divergence---yielding closed-form conditional moments and rank violation probabilities. The copula nests the comonotonic copula as $\rho \to 1$ and generalizes the Gaussian copula to the nonparametric setting.

The Bregman-Sinkhorn copula is nonparametrically identified by $(\mu, \nu, \rho)$ alone under mild regularity conditions on the marginals. Although rank violations occur with non-vanishing probability at every $\rho \in (0,1)$, the conditional rank of $Y_1$ given $Y_0$ remains sticky around the identity as $\rho \to 1$, with an explicit Bregman-divergence characterization of the conditional distribution. The relaxation of rank preservation incurs no loss in statistical efficiency: the empirical Bregman-Sinkhorn copula converges at the parametric $\sqrt{n}$-rate with a Gaussian process limit for any $\rho \in (0,1]$, despite the infinite-dimensional parameter space.

Applied to causal inference, efficient estimation of $\pi_\rho$ yields the entire treatment effect distribution $\cL(Y_1 - Y_0)$ and any functional of the joint distribution $(Y_0, Y_1)$. We derive a closed-form variance estimator for the ATE that is tighter than the Fr\'{e}chet--Hoeffding and Neyman bounds, and verify this in simulations that also examine sensitivity of the estimated TED to $\rho$. Empirically, we find the comonotonic copula ($\rho = 1$) simultaneously underestimates treatment effect heterogeneity and overestimates the variance of the ATE estimator; the Bregman-Sinkhorn copula corrects both. The framework extends to observational studies: under unconfoundedness, a covariate-specific stickiness parameter identifies the conditional joint distribution stratum by stratum, and the per-stratum transport map is computed via a gradient flow whose driving term can be estimated via binary classification on the observed data.

At a broader level, access to the joint distribution of potential outcomes goes well beyond what conditional average treatment effects can deliver athey2017state,athey2019machine: the conditional mean identifies which subgroups benefit on average FarrellLiangMisra2021,athey2021policy, but the TED reveals the probability that a given individual is harmed by treatment---information essential for targeting and policy design.

We conclude this section by collecting notation used throughout. Let $\cY$ be a compact subset of $\R^d$, $\sP(\cY)$ the set of Borel probability measures on $\cY$, and $\Pi(\mu, \nu)$ the set of couplings with marginals $\mu, \nu \in \sP(\cY)$. We denote by $\cL(Y)$ the law of a random variable $Y$ and by $T \sharp \mu$ the pushforward of $\mu$ under a measurable map $T$, defined by $(T \sharp \mu)(B) = \mu(T^{-1}(B))$ for every Borel set $B$. For scalar-valued potential outcomes $(Y_0, Y_1) \sim \pi \in \Pi(\mu, \nu)$, we write $F(y_0) = \mu((-\infty, y_0])$ and $G(y_1) = \nu((-\infty, y_1])$ for the marginal distribution functions, and $F^{-1}$, $G^{-1}$ for the corresponding quantile functions. The stickiness parameter is $\rho \in (0, 1]$, and we set $\epsilon(\rho) := (1-\rho)/\rho$ unless otherwise stated.

Rank Violations and Optimal Transport

For a given individual, let $Y_0$ and $Y_1$ denote the potential outcomes under control and treatment, respectively; the individual treatment effect is $Y_1 - Y_0$. Since $Y_0$ and $Y_1$ are never jointly observed for the same unit, the law of $Y_1 - Y_0$ is not identified without additional restrictions on the coupling $(Y_0, Y_1) \sim \pi \in \cP(\cY \times \cY)$.

A canonical identifying restriction for scalar outcomes is strict rank preservation: $$ (Y_0, Y_1) \sim \pi = (F^{-1}, G^{-1})\sharp U, \qquad U \sim \mathrm{Unif}(0,1) \;. $$ This is a strong restriction: it implies $Y_1 = G^{-1} \circ F(Y_0)$ almost surely. Equivalently, if a unit is at quantile $u$ of the control-outcome distribution, that unit is also at quantile $u$ of the treatment-outcome distribution. HeckmanSmithClements1997 argue that this assumption is often implausible. We therefore consider a weaker notion, which we call average rank preservation.

Optimal Transport and Average Rank Preservation

We define average rank preservation in the general multivariate setting ($\cY = \R^d$).

definition[Average rank preservation] A coupling $\pi \in \Pi(\mu, \nu)$ satisfies average rank preservation if \begin{align*} \sR(\pi) := \frac{1}{2} \E_{(Y_0, Y_1) \sim \pi} \E_{(Y'_0, Y'_1) \sim \pi} \langle Y_0 - Y'_0 , Y_1 - Y'_1 \rangle > 0 \;, \end{align*} where $(Y_0, Y_1)$ and $(Y'_0, Y'_1)$ are independent draws from $\pi$. This quantity is the expected inner product between the control-outcome gap and the treatment-outcome gap for two independent units.

The following notion ties the average rank preservation measure to the Kullback-Leibler divergence of $\pi$ from the independence coupling.

definition[$\rho$-rank preservation] For $\rho \in (0, 1]$ and $\epsilon(\rho) := (1-\rho)/\rho$, a coupling $\pi \in \Pi(\mu, \nu)$ satisfies $\rho$-rank preservation if \begin{align*} \sR(\pi) > \epsilon(\rho) \, \mathsf{KL}\bigl( \pi \mid \mu \otimes \nu \bigr) \;. \end{align*} $\rho$-rank preservation generalizes average rank preservation: the case $\rho = 1$ reduces to the condition $\sR(\pi) > 0$.

We now formalize the connection between rank-preservation and optimal transport (OT).

proposition[Rank correlation and entropic OT] Suppose $\mu, \nu \in \sP(\cY)$ with bounded second moments. \begin{enumerate} • Suppose $\mu$ is absolutely continuous with respect to the Lebesgue measure. Then the unique maximizer of $\sR(\pi)$ over $\Pi(\mu, \nu)$ is the Monge--Kantorovich coupling $\pi_1 = (\mathop{id}, \nabla\phi_1)\sharp\mu$, where $\nabla\phi_1$ is the Brenier map from $\mu$ to $\nu$. • For each $\rho \in (0,1)$ and $\epsilon(\rho) := (1-\rho)/\rho$, maximizing $\sR(\pi) - \epsilon(\rho) \,\mathsf{KL}(\pi\mid \mu\otimes\nu)$ over $\Pi(\mu,\nu)$ is equivalent to solving the entropic OT problem $$ \min_{\pi \in \Pi(\mu, \nu)} \int_{\cY \times \cY} \frac{1}{2}\| y_0 - y_1 \|^2 \dd \pi(y_0, y_1) + \epsilon(\rho) \, \mathsf{KL}(\pi \mid \mu \otimes \nu) \;. $$ \end{enumerate}

Rank Violations

Consider the scalar case with $\cY = \R$. The coupling maximizing $\rho$-rank preservation exhibits rank violations for every $\rho \in (0, 1)$.

For the case $\rho = 1$, maximizing average rank preservation $\sR(\pi)$ over $\Pi(\mu, \nu)$ uniquely identifies the Monge--Kantorovich coupling $\pi_1 = (\mathop{id}, G^{-1} \circ F)\sharp\mu$. This coupling exhibits perfect rank preservation: if a unit is at quantile $u$ of the control-outcome distribution with $Y_0 = F^{-1}(u)$, that unit is also at quantile $u$ of the treatment-outcome distribution with $Y_1 = G^{-1}(u)$.

For the case $\rho \in (0, 1)$, we will prove in Theorem (ref) that under mild regularity conditions on $\mu, \nu$, maximizing $\rho$-rank preservation over $\Pi(\mu, \nu)$ uniquely identifies a nonparametric coupling $\pi_\rho$

align[align omitted — 141 chars of source]

proposition[Rank violations] Let $\cY \in \R$ be a bounded interval, and let $\mu$ and $\nu$ be absolutely continuous with respect to the Lebesgue measure with strictly positive densities. For $\rho \in (0, 1)$, consider the coupling $(Y_0, Y_1) \sim \pi_\rho \in \Pi(\mu, \nu)$ uniquely defined in (ref). Then for any $u \in (0, 1)$, $$ \PP\Bigl( Y_1 \neq G^{-1}(u) \mid Y_0 = F^{-1}(u) \Bigr) > 0 \;. $$

Bregman-Sinkhorn Copula

Before formally deriving the nonparametric identification result under regularity conditions, we introduce the key properties of the Bregman-Sinkhorn copula, taking uniqueness of the solution to (ref) as given; the proof of uniqueness is deferred to Theorem (ref).

Sinkhorn Copula

Let $\epsilon(\rho) := (1-\rho)/\rho \in (0, \infty)$. The Sinkhorn dual potentials $\phi_\rho, \psi_\rho : \cY \to \R$ satisfy the system of equations

align[align omitted — 355 chars of source]

Under regularity conditions, the pair $(\phi_\rho, \psi_\rho)$ is uniquely identified by the marginals $\mu, \nu$ and the stickiness parameter $\rho$ (to be shown in Theorem (ref)), up to an additive constant $(\phi_\rho,\psi_\rho) \mapsto (\phi_\rho-c, \psi_\rho+c)$.

With these convex dual functions, we define the Sinkhorn copula $\pi_{\rho} \in \Pi(\mu, \nu)$ as

align[align omitted — 209 chars of source]

The name Sinkhorn copula comes from the fact that in the scalar-valued case,

align*[align* omitted — 233 chars of source]

where the density of the copula satisfies

align[align omitted — 220 chars of source]

Bregman Divergence and Conditional Exponential Family

The name Bregman-Sinkhorn copula reflects the following structure: under $\pi_\rho$, the conditional distribution of $Y_1$ given $Y_0$ is an exponential tilt of the marginal $\nu = \cL(Y_1)$, with a Bregman divergence as the exponent.

proposition[Bregman divergence representation] Consider $\cY = \R^d$. Let $(Y_0, Y_1) \sim \pi_\rho$ be the Sinkhorn copula defined in (ref), then for any $y_0 \in \cY$, the conditional distribution of $Y_1 \mid Y_0 = y_0$ is given by \begin{align} \PP( \dd y_1 \mid Y_0 = y_0 )= \frac{\exp\Bigl( \frac{ - B_{\psi_\rho}(y_1 \mid \nabla \psi^\ast_\rho(y_0))}{\epsilon(\rho)} \Bigr) \dd \nu(y_1) }{ \int_{\cY} \exp\Bigl( \frac{ - B_{\psi_\rho}(y'_1 \mid \nabla \psi^\ast_\rho(y_0))}{\epsilon(\rho)} \Bigr) \dd \nu(y_1') } \;, \end{align} where $B_{\psi_\rho}(y_1 \mid y'_1) := \psi_\rho(y_1) - \psi_\rho(y'_1) - \langle y_1 - y'_1, \nabla\psi_\rho(y'_1) \rangle$ is the Bregman divergence associated with $\psi_\rho$, and $\psi_\rho^\ast (y_0) := \sup_{y_1 \in \cY} \{ \langle y_0, y_1 \rangle - \psi_\rho(y_1) \}$ is the convex conjugate of $\psi_\rho$.

proof[Proof of Proposition (ref)] By definition of $\pi_\rho$ in (ref) $$ \PP( \dd y_1 \mid Y_0 = y_0 ) := \exp\Bigl( \frac{\langle y_0, y_1 \rangle - \phi_{\rho}(y_0) - \psi_{\rho}(y_1)}{\epsilon(\rho)} \Bigr) \dd \nu(y_1) \;. $$ By the Sinkhorn equation $\phi_\rho(y_0) = \epsilon(\rho) \log \int_{\cY} \exp\Bigl( \frac{\langle y_0, y_1 \rangle - \psi_{\rho}(y_1)}{\epsilon(\rho)} \Bigr) \dd \nu(y_1)$, we have \begin{align*} \PP( \dd y_1 \mid Y_0 = y_0 ) = \frac{\exp\Bigl( \frac{\langle y_0, y_1 \rangle - \psi_{\rho}(y_1)}{\epsilon(\rho)} \Bigr) \dd \nu(y_1) }{ \int_{\cY} \exp\Bigl( \frac{\langle y_0, y'_1 \rangle - \psi_{\rho}(y'_1)}{\epsilon(\rho)} \Bigr) \dd \nu(y_1') } = \frac{\exp\Bigl( \frac{ - B_{\psi_\rho}(y_1 \mid \nabla \psi^\ast_\rho(y_0))}{\epsilon(\rho)} \Bigr) \dd \nu(y_1) }{ \int_{\cY} \exp\Bigl( \frac{ - B_{\psi_\rho}(y'_1 \mid \nabla \psi^\ast_\rho(y_0))}{\epsilon(\rho)} \Bigr) \dd \nu(y_1') } \;. \end{align*} Here the last step follows from the fact that $\nabla \psi_\rho^\ast(y_0) = \argmax_{y_1} \{ \langle y_0, y_1 \rangle - \psi_\rho(y_1) \}$ and $\psi^\ast_\rho(y_0) = \langle y_0, \nabla \psi_\rho^\ast(y_0) \rangle - \psi_\rho(\nabla \psi_\rho^\ast(y_0))$, and thus $$\langle y_0, y_1 \rangle - \psi_\rho(y_1) - \psi^\ast_{\rho}(y_0) = - B_{\psi_\rho}(y_1 \mid \nabla \psi^\ast_\rho(y_0)) \;.$$

The conditional distribution of $Y_1 \mid Y_0 = y_0$ under the Sinkhorn copula is a member of the natural exponential family. Equivalently, letting $g(y_1) := \dd \nu(y_1) / \dd y_1$ denote the density of $\nu$, the conditional density of $Y_1 \mid Y_0 = y_0$ is

align[align omitted — 195 chars of source]

This is a member of the natural exponential family with natural parameter $y_0 / \epsilon(\rho)$ and log-partition function $\phi_\rho(y_0) / \epsilon(\rho)$. Consequently, the conditional mean and variance of $Y_1 \mid Y_0 = y_0$ are given by the first and second derivatives of $\phi_\rho$.

proposition[Conditional moments] Let $(Y_0, Y_1) \sim \pi_\rho$, the Bregman-Sinkhorn copula, and let $\phi_\rho$ be a solution to the Sinkhorn equation (ref), then, for any $y_0 \in \cY$, \begin{align*} \E[Y_1 \mid Y_0 = y_0] &= \nabla \phi_\rho(y_0) \;, \\ \mathop{Cov}[Y_1 \mid Y_0 = y_0 ] &= \epsilon(\rho) \; \nabla^2 \phi_\rho(y_0) \;. \end{align*}

Connection to Gaussian and Fréchet-Hoeffding Copulas

We restrict to the scalar case $\cY = \R$ and establish connections to the Fr\'{e}chet--Hoeffding and Gaussian copulas.

Consider the conditional Bregman-Sinkhorn copula: for any $\rho \in (0, 1]$ and $\epsilon \in [0, 1)$, the conditional density of $G(Y_1) = u_1 \mid F(Y_0) = u_0$ is given by

align[align omitted — 250 chars of source]

Here $\psi_\rho$ is the solution to the Sinkhorn system, and $\psi_\rho^\ast$ is its convex conjugate.

The Bregman-Sinkhorn copula nests two classical copulas as special cases:

(i) The Bregman-Sinkhorn copula generalizes the Gaussian copula to the nonparametric setting. The Bregman divergence reduces to a Mahalanobis distance when $\psi_\rho(x) = \rho x^2/2$, a quadratic function. With this parametric choice of $(\psi_{\rho}, \psi^\ast_\rho)$, $G = F = \Phi$ the standard Gaussian CDF, and the choice $\epsilon = (1-\rho^2)/\rho$, (ref) reduces to the conditional distribution of the Gaussian copula with correlation $\rho \in [0, 1)$

align*[align* omitted — 188 chars of source]

with joint copula density

align[align omitted — 241 chars of source]

(ii) The Bregman-Sinkhorn copula reduces to the comonotonic copula in the limit when $\rho \to 1$ and $\epsilon \to 0$. In this case, $\psi^\ast_1 = G^{-1} \circ F$ that achieves the Fréchet-Hoeffding copula, namely

align*[align* omitted — 74 chars of source]

Rank Stickiness

We now study the rank stickiness property of the nonparametric Bregman-Sinkhorn copula. The conditional rank distribution $G(Y_1) \mid F(Y_0) = u_0$ admits an analytic formula via the Bregman divergence representation in Proposition (ref). We characterize the rank distortion of the Bregman-Sinkhorn copula compared to the comonotonic copula, in the limiting case $\rho \to 1$ and $\epsilon \to 0$.

proposition[Rank stickiness] Consider the Bregman-Sinkhorn copula $(G(Y_1), F(Y_0)) \sim c_{\mathrm{B}}^{\rho, \epsilon}$ defined in (ref) with $\epsilon \in [0, 1)$ and $\rho \in (0, 1]$. Assume that potential $\psi_\rho$ is strictly convex and $\mu, \nu$ have strictly positive and smooth densities supported on a bounded interval. Then the following rank stickiness properties hold for any fixed $u_0 \in [0, 1]$ and $\rho \in (0, 1]$ and a sequence $\epsilon_n \to 0$: \begin{align*} \E[ G(Y_1) \mid F(Y_0) = u_0] &= G \circ \dot\psi^\ast_{\rho} \circ F^{-1} (u_0) + o_n(1) \\ \mathop{Var}[ G(Y_1) \mid F(Y_0) = u_0 ] &= o_n(1) \end{align*}

remarkWe call this rank stickiness for the following reason: when $\rho=1$ and $\epsilon_n$ small, the conditional rank of $G(Y_1) \mid F(Y_0) = u_0$ is concentrated around $G \circ \dot\psi^\ast_{1} \circ F^{-1} (u_0) = G \circ (G^{-1} \circ F) \circ F^{-1} (u_0) = u_0$, with bias and variance both of the order $o_n(1)$. Together with Proposition (ref), this shows that the Bregman-Sinkhorn copula exhibits rank stickiness in the sense that the conditional rank of $Y_1$ given $Y_0$ is near the true rank when $\rho \to 1$, while still allowing for rank violations with non-vanishing probability when $\rho \in (0, 1)$.

Identification and Inference

Identification

We now prove that the Bregman-Sinkhorn copula is nonparametrically identifiable and unique under the following assumptions, without imposing any parametric restrictions on the marginals $\mu, \nu$ or the copula $\pi_\rho$. We establish the result via an asymmetric contraction argument suited for the Sinkhorn system (ref).

assumption\begin{enumerate} • $\cY \subset \R^d$ is compact. • $\mu, \nu$ are absolutely continuous w.r.t. the Lebesgue measure and have strictly positive densities. \end{enumerate}

theorem[Nonparametric identifiability] Under Assumption (ref), the $\rho$-Bregman-Sinkhorn copula $\pi_\rho$ that solves (ref) is uniquely identifiable given the marginals $\mu, \nu$, for any $\rho\in (0, 1]$.

proof[Proof of Theorem (ref)] The case $\rho = 1$ follows from Brenier's theorem brenier1991polar. We hence consider $\rho \in (0,1)$ and denote $\epsilon := \epsilon(\rho) > 0$. By standard Sinkhorn theory for strictly positive continuous costs on compact spaces, there exists a bounded continuous pair $(\phi_\rho, \psi_\rho)$ solving \begin{align*} \phi_\rho(y_0) &= \epsilon \log \int_{\cY} \exp\Bigl( \frac{\langle y_0, y_1 \rangle - \psi_{\rho}(y_1)}{\epsilon} \Bigr) \dd \nu(y_1) \;, \\ \psi_\rho(y_1) &= \epsilon \log \int_{\cY} \exp\Bigl( \frac{\langle y_1, y_0 \rangle - \phi_\rho(y_0)}{\epsilon} \Bigr) \dd \mu(y_0) \;. \end{align*} It therefore suffices to prove uniqueness of such a pair up to additive constants: namely, the pair $(\phi_\rho, \psi_\rho)$ is the unique solution to the above system up to the transformation $(\phi_\rho, \psi_\rho) \mapsto (\phi_\rho - c, \psi_\rho + c)$ for any $c \in \R$. Let $(\phi^0, \psi^0)$ and $(\phi^1, \psi^1)$ be two bounded continuous solutions. Define \begin{align*} h(y_1) := \psi^1(y_1) - \psi^0(y_1) \;, \qquad M_\psi = \|h\|_\infty := \sup_{y_1 \in \cY} |h(y_1)| \;. \end{align*} We claim $h$ must be constant; suppose for contradiction that $h$ is non-constant. For $t \in [0,1]$, define \begin{align*} \psi^t := (1-t)\psi^0 + t\psi^1 \;, \qquad \phi^t(y_0) := \epsilon \log \int_{\cY} \exp\Bigl( \frac{\langle y_0, y_1 \rangle - \psi^t(y_1)}{\epsilon} \Bigr) \dd \nu(y_1) \;. \end{align*} Since $\psi^0,\psi^1$ are bounded and continuous and $\cY$ is compact, differentiation is justified, and for each $y_0 \in \cY$, \begin{align*} \phi^1(y_0) - \phi^0(y_0) = \int_0^1 \frac{\partial \phi^t(y_0)}{\partial t} \dd t = -\int_0^1 \int_{\cY} h(y_1) w_t(y_0,\dd y_1) \dd t \;, \end{align*} where \begin{align*} w_t(y_0,\dd y_1) := \frac{\exp\bigl( (\langle y_0,y_1\rangle - \psi^t(y_1))/\epsilon \bigr)}{\int_{\cY} \exp\bigl( (\langle y_0,y'_1\rangle - \psi^t(y'_1))/\epsilon \bigr) \dd \nu(y'_1)} \dd \nu(y_1) \end{align*} is a probability measure on $\cY$. Because $h$ is continuous and non-constant on the compact set $\cY$, there exist $\delta_\psi \in (0,1)$ and a nonempty open set \begin{align*} A_\psi := \{y_1 \in \cY: |h(y_1)| < (1-\delta_\psi) M_\psi\} \end{align*} with $\nu(A_\psi) > 0$. Next, let \begin{align*} R := \sup_{y_0,y_1 \in \cY} |\langle y_0,y_1\rangle| < \infty \;, \qquad B_\psi := \max\{\|\psi^0\|_\infty, \|\psi^1\|_\infty\} \;. \end{align*} Then for every $t \in [0,1]$ and $y_0,y_1 \in \cY$, \begin{align*} \exp\Bigl( \frac{-R-B_\psi}{\epsilon} \Bigr) \leq \exp\Bigl( \frac{\langle y_0,y_1\rangle - \psi^t(y_1)}{\epsilon} \Bigr) \leq \exp\Bigl( \frac{R+B_\psi}{\epsilon} \Bigr) \;. \end{align*} Hence, for every measurable $A \subset \cY$, \begin{align*} w_t(y_0,A) \geq \exp\Bigl( -\frac{2(R+B_\psi)}{\epsilon} \Bigr) \nu(A) =: c_\psi \nu(A) \;, \end{align*} uniformly in $y_0$ and $t$. Therefore, \begin{align*} \bigl|\phi^1(y_0)-\phi^0(y_0)\bigr| &\leq \int_0^1 \int_{\cY} |h(y_1)| w_t(y_0,\dd y_1) \dd t \\ &\leq \int_0^1 (1-\delta_\psi)M_\psi w_t(y_0,A_\psi) + M_\psi w_t(y_0,\cY \setminus A_\psi) \dd t \\ &\leq \bigl( 1 - \delta_\psi c_\psi \nu(A_\psi) \bigr) M_\psi \;. \end{align*} Taking the supremum over $y_0$ gives \begin{align} \|\phi^1 - \phi^0\|_\infty \leq \alpha_\psi \|\psi^1 - \psi^0\|_\infty \;, \qquad \alpha_\psi := 1 - \delta_\psi c_\psi \nu(A_\psi) < 1 \;. \end{align} The same interpolation applied to the second fixed-point equation gives, for each $y_1 \in \cY$, \begin{align*} \psi^1(y_1) - \psi^0(y_1) = -\int_0^1 \int_{\cY} \bigl(\phi^1(y_0) - \phi^0(y_0)\bigr)\, v_t(y_1, \dd y_0)\, \dd t \;, \end{align*} where $v_t(y_1, \dd y_0)$ is the analogous probability measure with $(\mu,\phi)$ in place of $(\nu,\psi)$. Since $v_t(y_1,\cdot)$ is a probability measure, \begin{align*} \|\psi^1 - \psi^0\|_\infty \leq \|\phi^1 - \phi^0\|_\infty \;. \end{align*} Combining with (ref) yields \begin{align*} \|\psi^1 - \psi^0\|_\infty \leq \alpha_\psi \|\psi^1 - \psi^0\|_\infty \;. \end{align*} Since $\alpha_\psi < 1$ and $M_\psi > 0$ ($\psi^0 \neq \psi^1$), this is a contradiction. Therefore $h$ must be constant, i.e., $\psi^1 - \psi^0 \equiv c$ for some $c \in \R$. Substituting $\psi^1 = \psi^0 + c$ into the first fixed-point equation gives, for every $y_0 \in \cY$, \begin{align*} \phi^1(y_0) = \epsilon \log \int_{\cY} \exp\Bigl( \frac{\langle y_0,y_1\rangle - \psi^0(y_1) - c}{\epsilon} \Bigr) \dd \nu(y_1) = \phi^0(y_0) - c \;. \end{align*} Hence any two solutions differ by additive constants of opposite sign. Finally, the coupling defined by (ref) is invariant under the transformation $(\phi,\psi) \mapsto (\phi-c, \psi+c)$. Therefore the resulting coupling $\pi_\rho$ is uniquely determined by $(\mu,\nu)$. This proves identifiability for $\rho \in (0,1)$, and the case $\rho=1$ was treated above.

Limit Theorem

Suppose $(Y_{0,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \mu$ and $(Y_{1,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \nu$ are independent samples. Let $\mu_n, \nu_n$ be the corresponding empirical measures, $\mu_n := \frac{1}{n} \sum_{i=1}^n \delta_{Y_{0,i}}$ and $\nu_n := \frac{1}{n} \sum_{i=1}^n \delta_{Y_{1,i}}$. Let $F_n, G_n$ be the empirical CDFs of $\mu_n, \nu_n$, and let $(\phi_n, \psi_n)$ be the empirical Sinkhorn potentials constructed from $\mu_n, \nu_n$ as in (ref).

We establish central limit theorems (CLT) for the empirical potentials $\phi_n \oplus \psi_n (y_0, y_1) := \phi_n(y_0) + \psi_n(y_1)$ constructed from the marginal empirical measures $\mu_n, \nu_n$. These complement the identifiability result of (ref): the Bregman-Sinkhorn copula is not only nonparametrically identified but also estimable at the parametric $\sqrt{n}$-rate, with a Gaussian process limit distribution.

The following two theorems use different proof strategies. For $\rho = 1$, the empirical rank-preserving map $G_n^{-1} \circ F_n$ is Hadamard differentiable as a composition of quantile functional, and the result follows from the functional delta method van1996weak. For $\rho \in (0,1)$, the Sinkhorn dual potentials are defined implicitly by the fixed-point system (ref); the main technical work is to establish (i) the Fr\'echet derivative of the Sinkhorn operator with respect to $h = \phi \oplus \psi$, and (ii) that the linearized operator (the derivative) is a bounded linear isomorphism on a gauge-centered Banach subspace, using a variant of the argument in Theorem (ref). The CLT then follows by applying the implicit function theorem in Banach spaces. We defer the proofs to the appendix.

We first establish the limit theorem for the rank-preserving case $\rho=1$, where $G_n^{-1} \circ F_n$ and $G^{-1} \circ F$ denote the empirical and population couplings, respectively.

theorem[CLT for empirical approximation, $\rho = 1$] Let $\cY \subset \R$ be a bounded interval. Denote $\ell^\infty(\cY)$ as the space of bounded measurable functions. Suppose $(Y_{0,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \mu$ and $(Y_{1,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \nu$ are independent samples. Assume $F$ and $G$ are continuously differentiable on $\cY$ with densities $f,g$, and \begin{align*} 0 < \inf_{y \in \cY} f(y) \leq \sup_{y \in \cY} f(y) < \infty \;, \qquad 0 < \inf_{y \in \cY} g(y) \leq \sup_{y \in \cY} g(y) < \infty \;. \end{align*} Then, in $\ell^\infty(\cY)$, we have the following limit: \begin{align*} \sqrt{n} \, \bigl( G_n^{-1} \circ F_n(y) - G^{-1} \circ F(y) \bigr) \rightsquigarrow \frac{\sqrt{2}\,\bB(F(y))}{g(G^{-1}(F(y)))} \;, \end{align*} for a standard Brownian bridge $\bB$.

Now we establish the limit theorem for the $\rho \in (0, 1)$ case. We prove for the general multivariate case $\cY \subset \R^d$ with $d \geq 1$ under Assumption (ref). Limit theorems for entropic OT have been established in the recent literature gonzalez2022weak, goldfeld2024limit, here we provide a self-contained analysis.

theorem[CLT for empirical approximation, $\rho \in (0,1)$] Let Assumption (ref) hold and consider $\rho \in (0,1)$. Denote $\cC(\cY)$ as the space of continuous functions. Suppose $(Y_{0,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \mu$ and $(Y_{1,i})_{i=1}^n \stackrel{i.i.d.}{\sim} \nu$ are independent samples. Let $(\phi_\rho,\psi_\rho)$ be Sinkhorn dual potentials for $(\mu,\nu)$, and $(\phi_n,\psi_n)$ the corresponding empirical dual potentials for $(\mu_n,\nu_n)$, as in (ref). Then, in $C(\cY) \oplus C(\cY) \subset C(\cY \times \cY)$, \begin{align*} \sqrt{n} \, \bigl( \phi_n \oplus \psi_n - \phi_\rho \oplus \psi_\rho \bigr) \rightsquigarrow \bZ_\rho := \bZ_{0, \rho} \oplus \bZ_{1, \rho} \;, \end{align*} where $(\bZ_{0,\rho}, \bZ_{1,\rho})$ is a jointly mean-zero Gaussian process in $C(\cY) \times C(\cY)$.

Application to Causal Inference

This section applies the Bregman-Sinkhorn copula to three problems in causal inference: (i) Section (ref) on estimation of the treatment effect distribution (TED), (ii) Section (ref) on estimation of the variance of the ATE estimator under rank violations, and (iii) Section (ref) on extending the framework to observational studies via covariate adjustment.

We adopt the potential outcomes framework. Let $(Y_0, Y_1) \in \cY \times \cY$ denote the potential outcomes under control and treatment, respectively, and let $T \in \{0, 1\}$ be a binary treatment indicator. The observed outcome is $Y_i = T_i Y_{1,i} + (1-T_i) Y_{0,i}$, and we observe $n$ i.i.d.\ copies $\{(Y_i, T_i)\}_{i=1}^n$. The object of interest is the treatment effect distribution (TED) $\cL(Y_1 - Y_0)$ and functionals thereof, including the variance of the Horvitz-Thompson estimator of the ATE.

Treatment Effect Distribution

We first work in the experimental setting where $T \sim \mathrm{Bern}(1/2)$ is independent of $(Y_0, Y_1)$. The joint distribution of potential outcomes is assumed to follow the Bregman-Sinkhorn copula: $(Y_0, Y_1) \sim \pi_\rho \in \Pi(\mu, \nu)$, where $\mu$ and $\nu$ are the marginal distributions of $Y_0$ and $Y_1$, respectively.

The Bregman-Sinkhorn copula $\widehat{\pi}_\rho$ is estimated by solving the Sinkhorn system with empirical marginals

align*[align* omitted — 186 chars of source]

where $\sC := \{i : T_i = 0\}$ and $\sT := \{i : T_i = 1\}$ are the control and treatment groups, with sizes $n_{\sC}$ and $n_{\sT}$, respectively. Since $T \perp\!\!\!\perp (Y_0, Y_1)$, we have $\cL(Y \mid T=0) = \cL(Y_0)$ and $\cL(Y \mid T=1) = \cL(Y_1)$, so the empirical marginals $\mu_n, \nu_n$ are consistent for $\mu, \nu$.

The Sinkhorn system (ref) is solved iteratively until convergence, yielding empirical dual potentials $(\phi_n, \psi_n)$. The empirical copula $\widehat{\pi}_\rho$ is then obtained by substituting $(\phi_n, \psi_n)$ and $(\mu_n, \nu_n)$ into (ref).

The TED is then estimated by pushing $\widehat{\pi}_\rho$ forward through the map $(y_0, y_1) \mapsto y_1 - y_0$: $$ \widehat{\cL_\rho(Y_1 - Y_0)} := (y_1 - y_0) \sharp \widehat{\pi}_\rho \in \sP(\cY) \;. $$

figure[figure omitted — 1,315 chars of source]

(ref) illustrates the imputed TED $\widehat{\cL_\rho(Y_1 - Y_0)}$ under two data-generating processes (DGPs): a Bregman-Sinkhorn copula DGP with $(F(Y_0), G(Y_1)) \sim c_{\mathsf{S}}^{\rho_1=0.99}$ as in (ref), and a Gaussian copula DGP with $(F(Y_0), G(Y_1)) \sim c_{\mathsf{G}}^{\rho_2=0.90}$ as in (ref). Both DGPs exhibit non-trivial rank violations, as evidenced by (ref). The top row of (ref) displays the sampling distribution of the ATE estimator $\widehat{\tau} = n_{\sT}^{-1} \sum_{i \in \sT} Y_i - n_{\sC}^{-1} \sum_{i \in \sC} Y_i$ across 1000 runs of the treatment assignment $T$, holding the potential outcomes fixed; the population ATE is a small positive constant in both cases.

For each DGP, $n = 2000$ i.i.d.\ draws of $(Y_0, Y_1, T)$ are generated, and the imputed TED $\widehat{\cL_\rho(Y_1 - Y_0)}$ is computed for stickiness levels $\rho \in \{0.90, 0.95, 0.99, 0.999, 1.00\}$. Rows 2--3 of (ref) compare the imputed TED $\widehat{\cL_\rho(Y_1 - Y_0)}$ against the true, unobserved $\cL_\rho(Y_1 - Y_0)$ for each DGP across these values of $\rho$.

Both DGPs are calibrated so that the population ATE $\tau := n^{-1}\sum_{i=1}^n (Y_{1,i} - Y_{0,i})$ is small and positive, yet the TED is bimodal: approximately $60\%$ of units experience a negative individual treatment effect and $40\%$ experience a positive one. The imputed TED $\widehat{\cL_\rho(Y_1 - Y_0)}$ recovers this bimodal structure across all stickiness levels $\rho$, and a well-chosen $\rho$ (e.g., $\rho = 0.99$) closely approximates the true population TED. By contrast, the ATE estimator fails to capture this heterogeneity and, under the couplings considered here, exhibits substantial variance, rendering ATE-based inference uninformative.

The comonotonic coupling ($\rho = 1$) underestimates treatment effect heterogeneity by producing an overly concentrated TED. The Bregman-Sinkhorn copula with $\rho < 1$ permits rank violations, interpolating between the comonotonic coupling ($\rho = 1$, which over-shrinks the TED) and the independence coupling ($\rho \to 0$, which over-disperses it), and recovering the full heterogeneity at intermediate values of $\rho$.

Variance Functional of ATE Estimator

The Horvitz-Thompson estimator of the ATE and its exact variance are

align*[align* omitted — 205 chars of source]

Since $\mathop{Var}(\widehat{\tau})$ depends on the unobserved joint distribution of $(Y_0, Y_1)$, it must be bounded or imputed under structural assumptions. Under the structural assumption $(Y_0, Y_1) \sim \pi_\rho$, the Bregman-Sinkhorn copula, the conditional moment formulas of Proposition (ref) yield a tractable expression for imputation.

proposition[Variance functional] Under the Bregman-Sinkhorn copula $(Y_0, Y_1) \sim \pi_\rho$, \begin{align*} \E\bigl[(Y_{0,i} + Y_{1,i})^2\bigr] = \E\bigl[ \bigl(Y_{0,i} + \dot\phi_\rho(Y_{0,i})\bigr)^2 + \epsilon(\rho)\,\ddot\phi_\rho(Y_{0,i}) \bigr] = \E\bigl[ \bigl(Y_{1,i} + \dot\psi_\rho(Y_{1,i})\bigr)^2 + \epsilon(\rho)\,\ddot\psi_\rho(Y_{1,i}) \bigr] \;. \end{align*} The resulting plug-in estimator of $\mathop{Var}(\widehat{\tau})$ is \begin{align*} \widehat{\mathop{Var}}_{\mathsf{SB}}(\widehat{\tau}) := \frac{1}{n^2} \sum_{i=1}^n \Bigl\{ \bigl[Y_i + \dot\phi_\rho(Y_i)\bigr]^2 + \epsilon(\rho)\,\ddot\phi_\rho(Y_i) \Bigr\} T_i + \Bigl\{ \bigl[Y_i + \dot\psi_\rho(Y_i)\bigr]^2 + \epsilon(\rho)\,\ddot\psi_\rho(Y_i) \Bigr\} (1-T_i) \;. \end{align*}

proofApply Proposition (ref) and the law of total expectation.

In practice, $\dot\phi_\rho$ and $\ddot\phi_\rho$ are computed from the empirical dual potentials $\phi_n, \psi_n$ obtained by running the Sinkhorn algorithm on $\mu_n, \nu_n$. By Proposition (ref), the empirical counterparts are (with $\dot\psi_\rho$ and $\ddot\psi_\rho$ defined analogously):

align*[align* omitted — 529 chars of source]

We benchmark $\widehat{\mathop{Var}}_{\mathsf{SB}}$ against two competitors. The first is the Fréchet-Hoeffding bound, corresponding to the special case $\rho = 1$ (rank preservation):

align*[align* omitted — 235 chars of source]

The second is the classical upper bound of neyman1923application, which relies only on the marginal variances:

align*[align* omitted — 293 chars of source]

figure[figure omitted — 1,157 chars of source]

(ref) reports a numerical comparison of all three bounds under two DGPs: (i) a Bregman-Sinkhorn copula DGP with $(F(Y_0), G(Y_1)) \sim c_{\mathsf{S}}^{\rho_0=0.95}$ and pronounced rank violations, and (ii) a Gaussian copula DGP with $(F(Y_0), G(Y_1)) \sim c_{\mathsf{G}}^{\rho_0=0.95}$ with moderate rank violations. Each bound is evaluated across 500 independent draws of $T$, with the potential outcomes held fixed, and compared against the true $\mathop{Var}(\widehat{\tau})$ computed from the unobserved $(Y_0, Y_1)$. The Bregman-Sinkhorn bound $\widehat{\mathop{Var}}_{\mathsf{SB}}$ is substantially tighter than both $\widehat{\mathop{Var}}_{\mathsf{FH}}$ and $\widehat{\mathop{Var}}_{\mathsf{N}}$, which overestimate the true variance across all replications.

The comonotonic coupling ($\rho = 1$) errs in both directions: it underestimates treatment effect heterogeneity by producing an overly concentrated TED, and overestimates the variance of the ATE estimator. The Bregman-Sinkhorn copula corrects both simultaneously: it recovers the full heterogeneity of the TED and yields a tighter variance estimator for the ATE.

Covariate Adjustment

In observational studies, treatment assignment $T$ may depend on covariates $X$, so $(Y_0, Y_1) \not\!\perp\!\!\!\perp T$. The Bregman-Sinkhorn framework extends naturally to this setting by conditioning on $X$.

assumption[Unconfoundedness] $(Y_0, Y_1) \perp\!\!\!\perp T \mid X = x$ for every $x \in \cX$.

assumption[Conditional Bregman-Sinkhorn copula] For every $x \in \cX$, the conditional distribution of $(Y_0, Y_1)$ given $X = x$ is a Bregman-Sinkhorn copula with stickiness parameter $\rho_x \in (0, 1]$: \begin{align*} (Y_0, Y_1) \mid X = x \;\sim\; \pi_{\rho_x}(\cdot \mid X = x) \in \Pi(\mu_x, \nu_x), \end{align*} where $\mu_x$ and $\nu_x$ denote the conditional distributions of $Y_0$ and $Y_1$ given $X = x$, respectively.

Assumption (ref) generalizes the conditional rank preservation assumption of chernozhukov2005iv and athey2006identification, which corresponds to $\rho_x = 1$ for every $x$.

Under Assumption (ref), the conditional marginals $\mu_x = \cL(Y_0 \mid X = x) = \cL(Y \mid X = x, T = 0)$ and $\nu_x = \cL(Y_1 \mid X = x) = \cL(Y \mid X = x, T = 1)$ are identified from the observed data. Under Assumption (ref), the conditional joint distribution $\pi_{\rho_x}(\cdot \mid X = x)$ is then identified as $$ \pi_{\rho_x}(\cdot \mid X = x) := \argmax_{\pi \in \Pi(\mu_x, \nu_x)} \bigl\{ \sR(\pi) - \epsilon(\rho_x)\,\mathsf{KL}(\pi \mid \mu_x \otimes \nu_x) \bigr\} \;. $$ The conditional TED $\cL(Y_1 - Y_0 \mid X = x)$ is recovered by pushing $\pi_{\rho_x}(\cdot \mid X = x)$ forward through $(y_0, y_1) \mapsto y_1 - y_0$, and the marginal TED $\cL(Y_1 - Y_0)$ is obtained by integrating over the distribution of $X$.

When $\rho_x = 1$ for every $x$, the conditional Bregman-Sinkhorn copula reduces to the conditional rank-preserving coupling, yet the marginal distribution of $(Y_0, Y_1)$ need not be rank-preserving---a feature unavailable to any single unconditional Bregman-Sinkhorn copula. To see this, let $X \in \{0, 1\}$ with equal probability, $(Y_0, Y_1) \mid X = 0 \sim \frac{1}{2} \delta_{(0.5,\, 0.25)} + \frac{1}{2} \delta_{(0.75,\, 1)}$, and $(Y_0, Y_1) \mid X = 1 \sim \frac{1}{2} \delta_{(0.25,\, 0.5)} + \frac{1}{2} \delta_{(1,\, 0.75)}$: each conditional is rank-preserving, yet the marginal is not. The conditional Bregman-Sinkhorn mixture therefore generates a strictly richer class of marginal dependence structures, accommodating heterogeneous rank dependence across covariate strata.

We now present an algorithm that computes the conditional rank-preserving Bregman-Sinkhorn copula $\pi_{\rho_x \equiv 1}(\cdot \mid X = x)$ jointly over all $x \in \cX$, for $\cY = \mathbb{R}$. The algorithm generalizes the parabolic Monge--Amp\`ere PDE of DebLiang2025, which computes the optimal transport map between $\cL(Y_0 \mid X = x)$ and $\cL(Y_1 \mid X = x)$. The parabolic Monge--Amp\`ere PDE defines a gradient flow that converges to $\dot{\phi}_\infty(y; x)$ as $t \to \infty$: with initial condition $\dot{\phi}_0(y; x) = y$,

align*[align* omitted — 97 chars of source]

where $h(\cdot\,; x)$ is the score of the conditional log-density ratio,

align*[align* omitted — 117 chars of source]

Setting the right-hand side to zero yields static Monge--Amp\`ere equations Liang2025DS2,Liang2025DS1 that characterize the optimal transport map $\dot{\phi}_\infty(\cdot\,; x)$ for each $x$: $$ \dot{\phi}_\infty(\cdot\,; x) \sharp \cL(Y_0 \mid X = x) = \cL(Y_1 \mid X = x) \;. $$

Under Assumption (ref), $h(\cdot\,; x)$ is identified from the observed data $(Y, X, T)$:

align*[align* omitted — 540 chars of source]

The joint log-likelihood ratio $\log \frac{\PP( \dot\phi_t(Y; X) = y, X = x \mid T = 0)}{\PP(Y = y, X = x \mid T = 1)}$ is identified from the observed data and can be estimated by training a binary classifier to predict $T$ from the pooled sample $\{ (\dot\phi_t(Y_i; X_i), X_i) : T_i = 0 \}$ and $\{ (Y_i, X_i) : T_i = 1 \}$. Formal consistency guarantees for the resulting plug-in estimator of the conditional TED are left for future work.

Acknowledgements

TL gratefully acknowledges support from the NSF Career Award (DMS-2042473) and the Wallman Society of Fellows at the University of Chicago.

\ifbiblatex \else \fi