Inference in partially identified moment models via regularized optimal transport

Inference in partially identified moment models via regularized optimal transport

\onehalfspacing

Introduction

Many quantities of interest in economics are not point-identified under realistic modeling choices. Two broad reasons can account for this: incomplete data, where relevant variables are only partially observed, and incomplete models, where economic theory does not pin down a unique data-generating mechanism. In such settings, the available data and modeling assumptions generally imply a set of parameter values rather than a unique point. The partial identification literature, such as the pioneering work by manski1990nonparametric, studies what can be learned about economically relevant parameters when the available information is insufficient for point identification; see also the handbook by manski2003partial and the survey by tamer2010partial.

Partial identification also often arises when the joint distribution of the data is unknown, even though the marginal distributions are observed. We consider a generalized method of moments (GMM) model where the parameter $\theta_0 \in \Theta \subset \mathbb{R}^k$ satisfies moment conditions

but the joint distribution $\pi_0$ of $(X,Y)$ is only known to lie in the set of couplings $\Pi(\mu,\nu)$ with observable marginals $\mu$ and $\nu$. To characterize the sharp identified set, fix a candidate $\theta$ and consider the set \[ \nu_{\Pi}(\theta) = \left\{\operatorname{\mathbb{E}}_\pi[\phi(X,Y,\theta)] : \pi \in \Pi(\mu,\nu)\right\} \] of all moment values generated by couplings consistent with the observed marginals. Then $\theta$ is compatible with the data if and only if $0 \in \nu_{\Pi}(\theta)$, or equivalently,

The last equality follows by norm duality and Sion's minimax theorem: see (ref) for details. The inner minimization is an optimal transport (OT) problem with linear costs $u'\phi(x,y,\theta)$, producing the coupling that minimizes the moment violation in direction $u$. The outer maximization identifies the worst-case direction. Thus, $\theta$ belongs to the sharp identified set exactly when this max-min value equals zero.

As a toy example, consider a randomized controlled trial (RCT) with potential outcomes $Y(0) \sim \mu$ and $Y(1) \sim \nu$. The share of units that benefit from treatment \( \theta = \mathbb P_\pi\left(Y(1) \ge Y(0)\right) \) is not point-identified since only the marginals $\mu,\nu$ of $\pi$ are observed. The moment function \( \phi(Y(0),Y(1),\theta)=\mathbf{1}\{Y(1) \ge Y(0)\} - \theta \) is one-dimensional, and $\theta$ belongs to the sharp identified set $\Theta_{I,0}$ if and only if there exists a coupling $\pi \in \Pi(\mu,\nu)$ such that \( \operatorname{\mathbb{E}}_{\pi}\left[\phi(Y(0),Y(1),\theta)\right] = 0. \) Equivalently, $\Theta_{I,0} = [\underline{\theta}, \overline{\theta}]$, where $\underline{\theta} = \min_{\pi \in \Pi(\mu,\nu)} \mathbb P_\pi\left(Y(1) \ge Y(0)\right)$ and $\overline{\theta} = \max_{\pi \in \Pi(\mu,\nu)} \mathbb P_\pi\left(Y(1) \ge Y(0)\right)$. This logic also extends to vector-valued or implicitly defined parameters.

In this paper, we develop a complete methodology for identification, estimation, and inference in the OT-based partially identified GMM model (ref). First, we characterize the sharp identified set for the parameter of interest using OT, as discussed above. The characterization provides both a geometric interpretation via support functions and informs an estimation procedure for the identified set.

Second, the implied estimation procedure involves solving an empirical version of the classical OT problem, which is known to be sensitive to sampling noise and having a slow convergence rate and a nonstandard limit distribution. To overcome this challenge, we employ entropic regularization, which penalizes the negative entropy of the joint distribution. Entropic OT has been widely used in statistics and machine learning since the seminal works by cuturi2013sinkhorn and galichon2022cupid, see the literature review below. This regularization yields a strictly convex problem that restores the usual $\sqrt{n}$-convergence and asymptotic normality, albeit at a cost of introducing a regularization bias. It is also computationally attractive, admitting fast implementation via the Sinkhorn algorithm for the inner (entropic OT) problem and the projected gradient ascent for the outer problem in the max-min representation (ref).

Third, we develop a procedure to test hypotheses and construct confidence regions for the identified set. To this end, we establish a uniform central limit theorem (CLT) for the entropic OT value uniformly in direction $u\in \mathbb{B}$ and parameter $\theta\in\Theta$. We build on mena2019statistical and goldfeld2024statistical and extend their framework to accommodate arbitrary smooth cost functions. To the best of our knowledge, we are the first to establish a CLT of such generality for the entropic OT value.

We then apply the functional delta method to the $\max$ functional to obtain the asymptotic distribution of our test statistic, relying on a result by franguridi2025set. This distribution depends on the unknown argmax set and is not available in closed form. Moreover, the standard bootstrap may be invalid when the hypothesized parameter value is on the boundary of the identified set, which occurs when the argmax set is not a singleton. We resolve this issue by employing the bootstrap for directionally differentiable functionals of fang2019inference. The resulting bootstrap-based test controls size locally uniformly and can be inverted to obtain a confidence region.

Finally, our estimator applies broadly to settings where one observes marginal distributions but lacks the knowledge of the joint distribution. In (ref), we discuss four examples: fixed effects panel logit with attrition and refreshment samples, nonlinear treatment effects, nonparametric instrumental variables (IV) without a large support condition, and the Euler equation with repeated cross-sections. In the first example, we exploit the panel structure by fixing the joint distribution of retainers and solving OT only for attriters, which tightens the bounds for the common slope coefficient and average marginal effects (AME).

\paragraph{Related literature.} First, our paper contributes to the extensive literature on partially identified models. Among others, imbens2004confidence and stoye2009more study confidence sets for partially identified parameters with uniform coverage and optimality properties, beresteanu2008asymptotic and bontemps2012set develop asymptotic theory and geometric characterizations for partially identified models using random sets and support functions, and romano2010inference provide general procedures for inference on identified sets defined via minimization of a criterion function. For an overview of the broader literature on partial identification and inference, we refer to the review by canay2017practical. The closest antecedent is beresteanu2011sharp, which characterizes sharp identification regions in models with convex moment predictions using the theory of random sets and support-function representations. Our characterization of the identified set is in this spirit but does not rely on representations via random sets.

Second, there has been a growing literature that applied OT methods to economic problems. The paper closest to ours is fan2025partial, which studies partial identification of a finite-dimensional parameter defined by a moment equality model with incomplete data via a classical OT representation of the identified set. In contrast, our paper employs entropic OT and also develops the methodology for estimation and inference, and hence goes beyond identification. More broadly, OT has been deployed in a variety of economic applications, such as discrete choice models chiong2016duality, covariate matching for causal effects gunsilius2021matching, nonlinear difference-in-differences for multivariate counterfactuals torous2024optimal, policy learning in matching markets hazard2025whom, and combining stated and revealed preferences meango2025combining. Other works include voronin2025generalized, which introduces a generalized version of OT for estimation in a large class of partially identified models, and schennach2025optimally, which uses OT for estimation and inference in the overidentified GMM with measurement errors.

Third, we employ entropic regularization for the classical OT, which yields a strictly convex program that can be solved efficiently by the Sinkhorn algorithm. Entropic OT is now widely used in high-dimensional statistics and machine learning. In seminal works, cuturi2013sinkhorn introduces Sinkhorn distances as a fast approximation to Wasserstein distances based on entropic regularization, and galichon2022cupid show how entropic regularization of social surplus yields efficient algorithms for estimating matching models. We also contribute to the statistical theory of the entropic OT and rely on two important prior works. mena2019statistical establish asymptotic normality of the entropic OT cost for quadratic cost functions, and goldfeld2024statistical extend this analysis using empirical process theory to derive semiparametric efficiency bounds and bootstrap inference procedures. Our CLT extends these results by allowing for arbitrary smooth cost functions, which is essential when the cost itself encodes economically meaningful restrictions.

Finally, our analysis for the panel logit with attrition and refreshment is closely related to recent contributions on panel surveys with nonignorable attrition and refreshment in franguridi2024estimation,franguridi2024robust. They develop computationally feasible and robust procedures for estimation and inference in this setting under the assumption of additively separable attrition that restores point identification. In contrast, we do not impose any assumptions on attrition, and hence our model is partially identified even with refreshment samples. We then show how to characterize and estimate bounds on the common slope parameter and the AME (see, e.g., davezies2021identification), as well as conduct inference in this setting.

The remainder of the paper is organized as follows. (ref) introduces the partially identified OT-based GMM model, and characterizes the sharp identified set for the parameter of interest. (ref) develops procedures for estimation and inference using entropic regularization. (ref) discusses several economic models that fit our general framework. (ref) conducts a Monte Carlo simulation for the fixed effects panel logit with attrition and refreshment. (ref) concludes. (ref) contains proofs of all the theoretical results, and (ref) provides additional details for the fixed effects panel logit with attrition and refreshment.

Setup and partial identification

Model and sharp identified set

Let $X$ and $Y$ be random vectors in $\mathbb{R}^d$ with the joint distribution $\pi_0$ identified only up to the class $\Pi(\mu,\nu)$ of joint distributions with marginals $\mu$ and $\nu$. The true value $\theta_0$ of a parameter of interest $\theta\in\Theta \subset \mathbb{R}^k$ uniquely satisfies the set of moment conditions

where $\phi$ is a $p$-dimensional moment function. Since $\pi_0$ is only partially identified, so is the parameter $\theta$ in general. We allow for arbitrary $\dim(\phi)=p$ and $\dim(\theta)=k$ as long as the identified set is nonempty and compact. For $p>k$, the additional moments introduce extra directions to detect violations and may tighten the identified set.

For each fixed parameter value $\theta$, consider the set of moment predictions \[ \nu_{\Pi}(\theta) = \left\{\operatorname{\mathbb{E}}_{\pi} \phi(X,Y,\theta): \pi\in\Pi \right\}, \] i.e., the set of all moment vectors consistent with the point-identified set of distributions $\Pi = \Pi(\mu,\nu)$. Because expectation is linear and $\Pi$ is convex, $\nu_{\Pi}(\theta)$ is convex. The sharp identified set is then \[ \Theta_{I,0} = \left\{\theta\in\Theta: 0\in\nu_{\Pi}(\theta)\right\} = \left\{\theta\in\Theta: D_0(\theta)=0\right\}, \] where $D_0(\theta)=d(0,\nu_{\Pi}(\theta))$ denotes the Euclidean distance from the origin to the set $\nu_{\Pi}(\theta)$. Identification is equivalent to the origin lying in the set of moment predictions beresteanu2011sharp. Note that, although $\nu_{\Pi}(\theta)$ is convex for each $\theta$, the identified set $\Theta_{I,0}$ need not itself be convex.

Let $\mathbb{B}$ be the unit ball in $\mathbb{R}^p$. By norm duality and Sion's minimax theorem, the distance $D_0(\theta)$ admits the following max-min representation

The inner minimization over couplings $\pi$ is an OT problem whose cost is the directional moment $u'\phi$. The outer maximization selects the direction $u$ in which the model slackness is largest. Since $D_0(\theta) = \max_{u\in \mathbb{B}}c_{\theta}(u)$ and $\theta \in \Theta_{I,0}$ if and only if $D_0(\theta) = 0$, we have that $\theta$ belongs to the sharp identified set exactly when $\max_{u\in \mathbb{B}}c_{\theta}(u) = 0$, or equivalently, when $c_{\theta}(u) \le 0$ for all $u \in \mathbb{B}$. This max-min representation informs our estimation and inference procedure based on estimating $c_{\theta}(u)$ for $u\in \mathbb{B}$ and checking whether the maximum is close to zero: see (ref).

Note that our setup is not a standard moment inequality model because the inner minimization over $\pi$ depends on the direction $u$ and delivers a direction-dependent evaluation of the set of predicted moments rather than a fixed collection of inequalities. It is also different from intersection bounds that aggregate separate scalar constraints.

To make the informal derivation above rigorous, we impose the following mild assumptions to ensure that the sharp identified set is nonempty and compact. Both of these properties are crucial for the Hausdorff consistency of the associated estimator, see (ref).

Toy example

To illustrate the characterization above, let us revisit the toy example in the introduction. Consider a RCT with potential outcomes $Y(0)\sim N(0,1)$ and $Y(1)\sim N(2,1)$. The parameter of interest is the share of units that benefit from treatment, \[ \theta=\mathbb P_\pi\left(Y(1)>Y(0)\right)=\operatorname{\mathbb{E}}_\pi\left[\mathbf{1}\{Y(1)>Y(0)\}\right]. \] The corresponding scalar moment condition is \( \operatorname{\mathbb{E}}_\pi\left[\phi\left(Y(0),Y(1),\theta\right)\right]=0, \) where $\phi\left(Y(0),Y(1),\theta\right)=\mathbf{1}\{Y(1)>Y(0)\}-\theta$. The identified set is then \[ \Theta_{I,0} :=\left\{\theta\in\mathbb{R}:\ \max_{u\in[-1,1]} c_\theta(u)=0\right\},\quad\text{where } c_\theta(u)=\min_{\pi\in\Pi} \operatorname{\mathbb{E}}_{\pi} \left[u\phi(Y(0),Y(1),\theta)\right]. \] This characterization has a simple interpretation as a two-sided game. Since $\phi$ is scalar, $u$ simply flips the sign of the moment: $u=1$ tests whether $\operatorname{\mathbb{E}}_\pi[\phi\left(Y(0),Y(1),\theta\right)]$ is positive under some $\pi\in\Pi$, while $u=-1$ tests the opposite. The adversary chooses the least favorable coupling under each sign. Thus, \[ \max_{u\in\{-1,0,1\}} c_\theta(u) = \max\left\{0,\,\,\, \min_{\pi\in\Pi} \mathbb P_{\pi}(Y(1)>Y(0)) - \theta,\,\,\, \theta - \min_{\pi\in\Pi} \mathbb P_{\pi}(Y(1)>Y(0))\right\} \] is the worst of these two one-sided checks, and equals $0$ at $u=0$. If neither side can produce a positive value, then $\theta$ is in the identified set $\Theta_{I,0}$.

For Gaussian marginals with common variance $\sigma^2$ and means $\mu_1\ge \mu_0$, the classical sharp bounds of makarov1982estimates yield \[ 1-2\Phi\left(\frac{-(\mu_1-\mu_0)}{2\sigma}\right) \le \mathbb{P}_\pi\left(Y(1)\ge Y(0)\right) \le\ 1, \] as also discussed in firpo2019partial. With $\sigma=1$, $\mu_0=0$, and $\mu_1=2$, the identified set is \( 0.69 \le \mathbb{P}_\pi\left(Y(1)\ge Y(0)\right) \le 1, \) which matches our characterization.

Figure (ref) plots $u\mapsto c_\theta(u)$ for $u\in[-1,1]$. The curve is piecewise linear, anchored at $u=1$ by the lowest feasible beneficiary share minus $\theta$ ($\min_{\pi\in\Pi} \mathbb P_{\pi}(Y(1)>Y(0))-\theta=0.69-\theta$), and at $u=-1$ by $\theta$ minus the highest feasible beneficiary share ($\theta-\max_{\pi\in\Pi} \mathbb P_{\pi}(Y(1)>Y(0))=\theta-1$). The identification check reduces to verifying whether this curve stays weakly below zero. For $\theta<0.69$, the right endpoint is above zero, whereas for $\theta\in[0.69,1]$, the whole curve remains nonpositive, and hence $\Theta_{I,0}=[0.69,1]$.

Estimation and inference

Suppose now that we have access to random samples $X_1,\dots, X_n$ and $Y_1,\dots, Y_m$ from distributions $\mu$ and $\nu$, respectively. For simplicity, we let $n=m$ and assume that the two samples are independent, although these assumptions can be relaxed at the expense of heavier notation. The goal of this section is to describe our estimator of the sharp identified set $\Theta_I$ and a testing procedure for the hypothesis $H_0:\theta=\theta_0$.

Estimation

The characterization of the sharp identified set in (ref) directly informs an estimation procedure. Denote by $\hat\mu$, $\hat\nu$ the empirical distributions based on samples $(X_i)$ and $(Y_j)$, and let $\hat\Pi = \Pi(\hat\mu,\hat\nu)$ be the set of joint distributions with marginals $\hat\mu$ and $\hat\nu$. The sample analog of $c_{\theta,0}(u)$ is then

This is the value of the empirical OT problem with the cost function $(x,y) \mapsto u'\phi(x,y,\theta)$. This is an infinite-dimensional linear program that is known to be computationally challenging, sensitive to sampling noise, and having nonstandard convergence rates when the (effective) dimensions of $X$ and $Y$ are greater than $4$, see, e.g., cuturi2013sinkhorn,hundrieser2024unifying. We therefore suggest using entropic regularization -- a classical technique for improving analytical and computational properties of OT that was introduced in cuturi2013sinkhorn and a working paper version of galichon2022cupid. This constitutes adding a term to the cost function that penalizes deviations from the independence distribution $\hat\mu\otimes\hat\nu$, viz.,

where $\operatorname{KL}$ is the Kullback-Leibler divergence,

The advantage of such regularization is that the program becomes strictly convex, much less sensitive to sampling noise, and regains standard asymptotics as we show in the next subsection. Importantly, the regularized program can be solved using explicit iterations of the Sinkhorn algorithm, see (ref). Although regularization introduces (small) bias, explicit debiasing procedures are available in the literature, see, e.g., pooladian2022debiaser and references therein.

The identified set of interest then becomes

where $\varepsilon>0$ is a penalty parameter and

is the population distance statistic. Throughout the rest of the paper, we drop the subscript $\varepsilon$ for brevity.

We define our estimator of $\Theta_I$ as

where $\eta_n>0$ is a tuning parameter and the (sample) distance statistic is

We impose the following assumptions.

(ref) imposes weak separation of the identified set by the criterion function. (ref) below implies that (ref) holds with $r_n = C n^{-1/2}$ for some constant $C$. Finally, (ref) requires picking $\eta_n$ that converges to zero sufficiently slowly.

The following theorem establishes the convergence of our estimator to the sharp identified set in the Hausdorff distance $d_H$.

Inference

Now our goal is to develop a test of the hypothesis $H_0:\theta=\theta_0$. To this end, we use the rescaled distance $\sqrt{n}\cdot \hat D(\theta_0)$ as the test statistic. We characterize its asymptotic behavior by first establishing a novel CLT for the regularized OT value $\hat c_{\theta_0}(u)$, uniformly in $u\in \mathbb{B}$ and $\theta_0\in\Theta$, and then applying the functional delta method to obtain the limiting distribution of our test statistic. Since this distribution does not have a simple form, we employ a result in franguridi2025set to establish the validity of the bootstrap for directionally differentiable functionals of fang2019inference.

Our first result is the uniform CLT for the regularized OT value, which we derive under the following assumptions.

(ref) helps establish uniform bounds on the optimal potentials and their derivatives. At the expense of more complicated proofs, this assumption can be relaxed to bounds on the tails of $\mu,\nu$ related to the cost function, similar to goldfeld2024statistical. (ref) can be relaxed to accommodate dependence within samples, e.g., when the samples are stationary $\beta$-mixing processes, as well as dependence across samples, see Remark 11 in goldfeld2024statistical.

Applying the functional delta method to this uniform convergence statement with the functional $\chi(c) :=\max_{u\in \mathbb{B}} c(u)$ leads to the following result.

The asymptotic distribution of the distance statistic is neither available in closed form, nor is easy to simulate from, because it depends on the unknown features of the data-generating process such as the argmax set $U_c(\theta_0)$. Moreover, standard bootstrap often fails to control size uniformly in partially identified models andrews2009validity,andrews2009invalidity. This failure occurs when the parameter is on the boundary of the parameter space andrews2000inconsistency. In our model, this corresponds to the case where the hypothesized value $\theta_0$ is on the boundary of the identified set or, equivalently, where $U_c(\theta_0)$ is not a singleton, such as $\theta\in\{0.69,1.00\}$ in (ref).

To overcome this challenge, we make use of the bootstrap for directionally differentiable functionals of fang2019inference. Our testing procedure is described in (ref) and depends on an additional tuning parameter $\iota_n$. We impose the following assumption.

(ref) posits that the maxima of $c_{\theta_0}$ are well-separated. This assumption is equivalent to the (super)gradient of $c_{\theta_0}$ being bounded away from zero on the complement of the argmax set $U_c(\theta_0)$. It also suffices for this assumption to hold in a small neighborhood around $U_c(\theta_0)$ rather than on the entire ball $\mathbb{B}$. (ref) requires $\iota_n$ to converge to zero slower than $n^{-1/2}$. This guarantees that the enlarged argmax $\hat U_n$ converges in the Hausdorff distance to the true argmax $U_c(\theta_0)$. If $U_c(\theta_0)$ is known to be a singleton, we can take $\hat U_n$ to be the sample argmax of $\hat c_{\theta_0}$.

It is straightforward to establish that under (ref), our testing procedure controls size locally uniformly in the sense of Corollary 3.2 of fang2019inference. We refer the reader to Theorem 4 of franguridi2025set for details.

Numerical implementation

Our procedure requires computing the distance statistic (ref), which amounts to solving two nested optimization problems.

The inner problem (entropic OT) is solvable by a very fast and simple numerical procedure called the Sinkhorn algorithm described in Algorithm (ref).

The outer problem is the maximization of a concave function $\hat c_\theta(u)$ (the output of (ref)) over the unit ball $u\in\mathbb{B}$. We solve this problem using projected gradient ascent described in (ref). Of course, many other off-the-shelf algorithms are available for constrained concave optimization; for a comprehensive review, see bubeck2015convex.

Illustrative examples

This section presents four empirical examples that illustrate how our methods can address partial identification problems arising from missing or incomplete data. The examples span different areas of econometrics: panel data, causal inference, and macro-finance. In each case, we show how the inability to observe certain joint distributions leads to partial identification of parameters of interest and how our methodology can be applied to characterize the identified sets in various contexts.

Fixed effects panel logit with attrition and refreshment

This is our leading example and forms the basis for the Monte Carlo simulations in (ref). Panel logit models are widely used in empirical studies to analyze binary outcomes while controlling for individual heterogeneity through fixed effects. However, attrition is a common issue in panel data, where units may drop out of the sample in subsequent periods for reasons potentially correlated with outcomes of interest. When a refreshment sample is available, the model can be partially identified via OT. With an adjustment that exploits the panel structure, our methodology can be used to derive tight bounds for the common slope parameter and the AME.

Common slope parameter

The common slope parameter is point identified under complete data or when attrition is independent of outcomes conditional on observables, but becomes partially identified under unrestricted attrition. For notation simplicity, consider a static panel logit model with two periods $T=2$,

where $\theta$ captures the effect of covariates on the outcome, $\alpha_i$ represents individual fixed effects, and $\varepsilon_{it}$ follows a standard logistic distribution.

The standard approach to eliminate the incidental parameter $\alpha_i$ is to condition on the sufficient statistic $S_i = Y_{i1} + Y_{i2}$. For individuals with $S_i = 1$ (switchers), the conditional log-likelihood is

with the corresponding conditional score function $s(Y_{i1},Y_{i2},X_{i1},X_{i2};\theta)$ and the moment condition

The explicit form of the score is given in (ref). To translate this moment condition into our OT framework, we embed the event $\{Y_{i1}+Y_{i2}=1\}$ into the cost function and define

A naive approach to partially identify $\theta$ would use only the marginal distributions from period 1 (original sample) and period 2 (refreshment sample). However, we can achieve tighter bounds by exploiting the panel structure. Since we observe both periods for retainers, we fix their joint distribution $f_{1,2\mid\text{ret}}$ and apply OT only to couple the attriter distributions $f_{1\mid\text{att}}$ and $f_{2\mid\text{att}}$ across periods. While $f_{1\mid\text{att}}$ is directly observed, $f_{2\mid\text{att}}$ is unobserved because attriters are missing in period 2, but it can be recovered from the observed refreshment and retainer samples via the law of total probability. See (ref) for details.

Let $p$ denote the retention rate. The attriter contribution to the moment bounds is

Since $\operatorname{\mathbb{E}}_{f_{1,2\mid\text{ret}}}[\phi(Y_1,Y_2,X_1,X_2;\theta)]$ can be computed directly from the observed data for retainers, the overall bounds are

The identified set for $\theta$ is therefore $\Theta_{I,0} = \{\theta :\, \underline{\nu}(\theta) \leq 0 \leq \overline{\nu}(\theta)\}.$ (ref) in (ref) presents our estimation procedure.

Our estimator naturally extends to dynamic panel logit models where lagged dependent variables appear as regressors. It is possible to incorporate the moment conditions developed by honore2024moment as the cost function, with a similar but more complicated partition structure separating retainers and attriters to achieve efficiency gains.

AME

We now consider estimation of the AME, which captures the average change in outcome probability induced by a marginal change in a covariate and provides a directly interpretable measure for empirical work. Specifically, the AME of covariate $j$ at period $\tau$ is defined as \[ \delta_{\tau,j} = \theta_j\operatorname{\mathbb{E}}[\Lambda(X_{\tau}'\theta + \alpha)(1-\Lambda(X_{\tau}'\theta + \alpha))]. \] The AME is partially identified even without attrition due to the incidental parameters problem combined with the nonlinear structure of the logit model. Under unrestricted attrition, this identification issue becomes more severe as the joint distribution of outcomes across periods is no longer observable.

davezies2021identification show how to construct outer bounds on the AME without attrition. They show that $\delta_{\tau,j}$ belongs to the interval $\tilde{\delta} \pm \bar{b}$, where $\tilde{\delta} = \operatorname{\mathbb{E}}[p(X_{1:T},S,\theta_0)]$ and $\bar{b} = \operatorname{\mathbb{E}}[a(X_{1:T},S,\theta_0)]$ for functions $p$ and $a$ constructed from a degree-$(T+1)$ Chebyshev polynomial. The explicit formulas for $p$, $a$, and the associated coefficients are collected in (ref), including closed-form expressions for the case $T=2$.

Under unrestricted attrition, we extend the partition approach for $\theta$ to bound the AME. First, we compute the identified set $\widehat\Theta_I$ for the common slope parameters $\theta$ as in (ref) and construct a finite grid $\left\{\theta^{(g)}\right\}\subset \widehat\Theta_I$. Second, for each grid point $\theta^{(g)}$ we plug it into the Chebyshev approximation and, using our partition of retainers and attriters, compute the corresponding AME bounds $\left[\underline{\delta}(\theta^{(g)}), \overline{\delta}(\theta^{(g)})\right]$ via OT. Finally, we profile over $\theta$ by taking the union to obtain the identified set for the AME: \( \bigcup_g \left[\underline{\delta}(\theta^{(g)}), \overline{\delta}(\theta^{(g)})\right]. \) See (ref) for details.

Although this grid-based profiling yields conservative AME bounds due to the Chebyshev polynomial approximation and the two-stage approach that first estimates bounds on $\theta$ and then bounds on $\delta$, it delivers a significant improvement in computational efficiency over alternative methods such as Hankel moment matrix positivity bounds or a single-stage approach that estimates both $\theta$ and $\delta$ simultaneously.

Nonlinear treatment effects

In causal inference, researchers face the fundamental problem of missing data because they observe either the potential outcome under treatment $Y(1)$ or under control $Y(0)$ for each individual, but never both simultaneously. This example demonstrates how OT methods provide identified sets when the object of interest is a functional of the joint distribution of potential outcomes. See also for example fan2025partial.

Consider a binary treatment $D\in\{0,1\}$, outcome $Y$, and pre-treatment covariates $X$. Let $Y(d)$ denote the potential outcome under treatment status $d$, with the observed outcome being $Y = DY(1) + (1-D)Y(0)$. We impose the standard assumptions $(Y(0),Y(1)) \perp D \mid X$ (strong ignorability) and $0 < e(X) < 1$ almost surely (monotonicity), where $e(X) = P(D=1\mid X)$ is the propensity score.

The parameter of interest is \[ \theta = \operatorname{\mathbb{E}}\left[h(Y(1),Y(0))\right], \] where $h: \mathbb{R}^2 \to \mathbb{R}$ is a known function of both potential outcomes. This encompasses parameters such as the average treatment effect (ATE) $h(y_1,y_0)=y_1-y_0$, distributional treatment effects $h(y_1,y_0)=\mathbf{1}\{y_1-y_0\leq t\}$, the variance of treatment effects $h(y_1,y_0)=(y_1-y_0)^2$, and the share of units that benefit from treatment $h(y_1,y_0)=\mathbf{1}\{y_1 > y_0\}$. While the ATE is point identified, parameters involving the joint distribution of $(Y(1), Y(0))$ are generally only partially identified since we never observe both potential outcomes for the same individual.

To implement our OT approach, we employ propensity score stratification. By the propensity score theorem of rosenbaum1983central, we have $(Y(0),Y(1)) \perp D \mid e(X)$. Rather than conditioning on the covariate vector $X$ directly, we discretize the propensity score into $K$ strata, $\left\{I_k\right\}_{k=1}^K$, with equal probability $\pi > 0$ for each stratum. The propensity score stratification reduces the curse of dimensionality when $X$ is high-dimensional, ensures sufficient sample sizes, and facilitates the implementation of OT estimation within each stratum.

Within each stratum $k$, the conditional distributions $F_{Y(d) \mid e(X) \in I_k}(y)$ for $d \in \{0,1\}$ are point identified. However, the joint distribution $F_{Y(1), Y(0) \mid e(X) \in I_k}(y_1, y_0)$ remains unknown, leading to partial identification of stratum-specific parameters $\theta_k = \operatorname{\mathbb{E}}[h(Y(1),Y(0)) \mid e(X) \in I_k]$. The identified set for $\theta_k$ is given by the OT bounds \( \Theta_k = \left[\underline\theta_k, \overline\theta_k\right], \) where the cost function is $\phi(y_1,y_0) = h(y_1,y_0)$,

and $F_{d \mid k} = F_{Y(d) \mid e(X) \in I_k}$ denotes the marginal distribution of $Y(d)$ in stratum $k$. The parameter $\theta$ is then bounded by $\theta \in \pi\left[ \sum_{k=1}^K \underline\theta_k,\, \sum_{k=1}^K \overline\theta_k\right]$.

The bounds derived above are sharp for many functions $h$ of interest. As shown in fan2017partial, when $h$ is either a supermodular functional or an indicator of a nondecreasing transformation, the sharp bounds are given by the Fr\'echet-Hoeffding bounds, which are achieved by comonotone and counter-comonotone couplings (corresponding to perfect positive and negative dependence, respectively).

Nonparametric IV without large support

IV methods are commonly used in causal inference when strong ignorability fails. The classical control function approach to IV requires a large support condition for point identification. However, this assumption often fails in practice when the treatment is discrete or has limited variation, such as binary treatments in angrist1996identification and discrete treatment intensity in angrist1995two. This example shows how our methodology can deliver meaningful bounds even when the large support assumption fails.

Consider the standard nonparametric IV model \[ Y = g(X,V), \quad X = h(Z,W), \quad (W,V)\perp Z, \] where $Y$ is the outcome, $X$ is the endogenous treatment, $Z$ is the instrument, and $(V,W)$ are the unobserved errors. We want to estimate a generic object of interest \(\theta=\operatorname{\mathbb{E}}[\Lambda(g(X,V))]\) beyond the local average treatment effects (LATE), e.g., $\Lambda(y) = y$ for the ATE, or $\Lambda(y) = \mathbf{1}\{y \leq t\}$ for the distributional treatment effect. The classical control variable formula requires not only treatment monotonicity ($w \mapsto h(z,w)$ strictly increasing) but also large support ($\operatorname{\text{supp}}(R\mid V=v) = \operatorname{\text{supp}}(R)$ for all $v$, where $R = F_{X\mid Z}(X)$), see, e.g., Section 4.1 in gunsilius2025primer. When this large support condition fails, $\theta$ becomes partially identified.

To formalize this setting, we impose two assumptions. First, similar to the classical setup, we assume treatment monotonicity: $w \mapsto h(z,w)$ is strictly increasing for each $z$, so under a monotone transformation, we can define the control variable $W=F_{X\mid Z}(X)$. Second, we also assume outcome monotonicity: $v \mapsto g(x,v)$ is strictly increasing, so define $V = F_{Y\mid X}(Y) \sim \text{Uniform}[0,1]$. Note that $W$ and $V$ are, in general, not independent. If $W$ has limited variation, such as discrete or censored treatment, the large support condition may fail, i.e., $\operatorname{\text{supp}}(W\mid V=v) \subsetneq \operatorname{\text{supp}}(W)$ for some $v$.

In many empirical applications, the instrument $Z$ is supported on a finite number of values. For example, angrist1995two use quarter-of-birth dummies as instruments for schooling, and card1995geographic uses a binary instrument based on proximity to four-year colleges. For each $z$, the marginals $F_{W}$ and $F_V$ can be recovered from the data, and the moment function is \[ \phi(w,v;z) = \Lambda\left(g\left(F_{X\mid Z=z}^{-1}(w), v\right)\right). \] Then the sharp identified set for $\theta$ is the interval $\left[\underline\theta, \overline\theta\right]$ with \[ \underline\theta = \inf_{F\in\Pi(F_{W}, F_V)} \sum_z\operatorname{\mathbb{E}}_F [\phi(W,V;z)]\Pr(Z=z),\quad \bar\theta = \sup_{F\in\Pi(F_{W}, F_V)} \sum_z\operatorname{\mathbb{E}}_F [\phi(W,V;z)]\Pr(Z=z). \]

Euler equation estimation with repeated cross-sections

A prominent example in macro-finance is estimating the discount factor $\beta$ and risk aversion $\gamma$ from the constant relative risk aversion (CRRA) Euler equation \[ \operatorname{\mathbb{E}}\left[\beta\left(C_{i,t+1}/C_{it}\right)^{-\gamma}R_{t+1} - 1\mid \mathcal I_{it}\right] = 0, \] where $C_{it}$ is individual $i$'s consumption at time $t$, $R_{t+1}$ is the common asset return, and $\mathcal I_{it}$ is the information set.

In practice, this single nonlinear conditional moment is converted into an overidentified unconditional GMM by introducing a $k$-dimensional vector of instruments with $k\ge2$, \( Z_{it} = \left(Z_{t}^A, Z_{it}^I\right), \) where $Z_t^A$ are lagged macro variables (such as GDP growth rates and interest rates), and $Z_{it}^I$ are lagged individual variables (such as demographics, as well as prior income and consumption). The moment condition then becomes \[ \operatorname{\mathbb{E}}\left[Z_{it}\left(\beta\left(C_{i,t+1}/C_{it}\right)^{-\gamma}R_{t+1} - 1\right)\right] = 0. \] Under a standard rank condition, this equation can be used to estimate $(\beta,\gamma)$.

To estimate the Euler equation, one would ideally use data that preserve cross-sectional heterogeneity while providing a sufficient sample size, which in practice motivates the use of repeated cross-sections or short rotating panels. Much of the early empirical literature, however, relied on aggregate time-series consumption and return data, as in hansen_singleton1982. Because the Euler equation is nonlinear in consumption, Jensen's inequality implies that the nonlinear moment evaluated at aggregate consumption would differ from the cross-sectional average of individual-level terms. This mismatch can induce systematic bias in GMM estimates. Subsequent work emphasized granular data to retain heterogeneity. For example, dynan_skinner_zeldes2004 exploit the panel structure of the Panel Study of Income Dynamics (PSID), but the PSID has a relatively small cross-sectional sample size. In contrast, many large-scale household surveys, such as the Consumer Expenditure Survey, offer rich cross-sectional coverage via repeated cross-sections or short rotating panels, making them well-suited for our OT-based approach. See also liu2023full on estimating full structural models with repeated cross-sections.

In an extreme case, suppose we observe only repeated cross-sections, so that the marginal distributions \( f_{C_{it},Z_{it}^I}(c,z^I) \text{ and } f_{C_{i,t+1}}(\tilde c), \) are known, in contrast to their joint law. Let $\theta = (\beta,\gamma)'$. Then the parameter is only partially identified with the identified set \( \Theta = \left\{\theta :\, \underline\nu(\theta)\le0\le \overline\nu(\theta)\right\}, \) where \[ \underline\nu(\theta) = \inf_{f\in\Pi\left(f_{C_{it},Z_{it}^I}, f_{C_{i,t+1}}\right)} \operatorname{\mathbb{E}}_f\left[\phi(Z_{it},C_{it},C_{i,t+1},R_{t+1};\theta)\right], \] and $\overline\nu(\theta)$ is the supremum of the same expression. Here $\Pi(\cdot)$ is the set of all couplings consistent with the observed marginals, and the moment function is \[ \phi(z,c,\tilde c,r;\theta) = z'\left(\beta\left(\tilde c/c\right)^{-\gamma}r - 1\right). \]

With short rotating panels, where households are observed for only several periods, the OT problem becomes more complex: one can exploit the limited longitudinal links to tighten the identified set while using OT for the remaining unlinked portions of the data.

Monte Carlo simulation

In this section, we give a small illustration of the performance of our inference procedure.

The data-generating process is the fixed effects panel logit with attrition and refreshment as described in (ref) and defined in (ref). The true parameter is $\theta_0=(1.0,2.0)'$. The covariate vector $X_{it}=(X_{it,1},X_{it,2})'$ consists of two independent components $X_{it,1}$ and $X_{it,2}$ that have discrete uniform distributions on three-valued sets $\mathcal{X}_1 = \{0.42, 0.55, 0.60\}$ and $\mathcal{X}_2 = \{0.54, 0.65, 0.72\}$, respectively, and are independent across units $i$ and time $t$. The fixed effects $\alpha_i$ are drawn i.i.d.\ from the standard normal distribution. The idiosyncratic error terms $\varepsilon_{it}$ are drawn i.i.d.\ from the standard logistic distribution with zero mean and unit scale. The sizes of both the first-period and the refreshment samples are $n_\text{org}=n_\text{ref}=15000$. The attrition rate is $1-p=10\%$, and the units drop out of the sample completely at random.

We conduct $400$ simulations. For each simulation, we construct the 90% confidence region by inverting the test in (ref) on a grid of hypothesized values $\theta^* \in [-0.5,2.5]\times [0.5,3.5]$ centered at the true value $\theta_0=(1.0,2.0)'$. We use (ref) to solve the discrete entropic OT problem on a grid, where each of the two marginals is defined on $2\times 3\times 3$ values in the joint support $\{0,1\}\times \mathcal{X}_1\times\mathcal{X}_2$ of $(y_{it},x_{it1},x_{it2})$. We then use grid search on the unit circle to calculate the distance $\hat D(\theta^*)$. We set the entropic regularization parameter $\varepsilon=0.1$ and the tuning parameter $\iota_n=0.05 n^{-1/2}\log n$ for constructing the enlarged argmax set in the bootstrap procedure.

(ref) shows the simulation results. Each value $\theta^*$ in the grid is colored according to the value of $\hat D(\theta^*)$ (left panel) or the proportion of times it is covered by the 90% confidence region (right panel). The identified set $\Theta_{I,\varepsilon}$ is indicated by the black contour line. We see that the distance statistic tracks the identified set and the confidence region performs reasonably well.

Conclusion

This paper develops a methodology for estimation and inference in GMM where the distribution of the data is identified only up to its marginals. We characterize the identified set for the parameter of interest using tools from convex analysis and OT. The practical implementation of classical OT is hindered by both theoretical and computational limitations. To overcome these issues, we rely on the regularized (entropic) version of OT. The resulting OT-based characterization directly informs an estimator and a test statistic for conducting inference and constructing confidence regions.

We establish a central limit theorem for the entropic OT value under smooth cost functions and use it to show $\sqrt{n}$-consistency and asymptotic normality of our proposed statistic. We then obtain valid critical values via the bootstrap for directionally differentiable functionals developed in fang2019inference.

Our estimation and inference methodology is generic and computationally efficient. It is also highly relevant for applied work, since many important economic questions, ranging from the effects of policy interventions to the dynamics of household behavior, involve parameters that are characterized via OT-based partially identified GMM.

Our framework admits several promising theoretical extensions. First, it naturally extends to settings with more than two marginals, such as panel data with multiple waves or repeated cross-sections over several periods (multi-marginal OT). Second, it would be of interest to extend it to GMM models with conditional moment restrictions. Finally, when the moment conditions underidentify the parameter even under point identification of the data distribution, the interaction between these two sources of partial identification is theoretically challenging and deserving of further study. We leave these extensions for future work.