Sampling of probability measures in the convex order by Wasserstein projection

Motivated by the approximation of Martingale Optimal Transport problems, we study sampling methods preserving the convex order for two probability measures $\mu$ and $\nu$ on $\mathbb{R}^d$, with $\nu$ dominating $\mu$. When $(X_i)_{1\le i\le I}$ (resp. $(Y_j)_{1\le j\le J}$) are i.i.d. according $\mu$ (resp. $\nu$), the empirical measures $\mu_I$ and $\nu_J$ are not in the convex order. We investigate modifications of $\mu_I$ (resp. $\nu_J$) smaller than $\nu_J$ (resp. greater than $\mu_I$) in the convex order and weakly converging to $\mu$ (resp. $\nu$) as $I,J\to\infty$. In dimension 1, according to Kertz and R\"osler (1992), the set of probability measures with a finite first order moment is a lattice for the increasing and the decreasing convex orders. From this result, we can define $\mu\vee\nu$ (resp. $\mu\wedge\nu$) that is greater than $\mu$ (resp. smaller than $\nu$) in the convex order. We give efficient algorithms permitting to compute $\mu\vee\nu$ and $\mu\wedge\nu$ when $\mu$ and $\nu$ are convex combinations of Dirac masses. In general dimension, when $\mu$ and $\nu$ have finite moments of order $\rho\ge 1$, we define the projection $\mu\curlywedge_\rho \nu$ (resp. $\mu\curlyvee_\rho\nu$) of $\mu$ (resp. $\nu$) on the set of probability measures dominated by $\nu$ (resp. larger than $\mu$) in the convex order for the Wasserstein distance with index $\rho$. When $\rho=2$, $\mu_I\curlywedge_2 \nu_J$ can be computed efficiently by solving a quadratic optimization problem with linear constraints. It turns out that, in dimension 1, the projections do not depend on $\rho$ and their quantile functions are explicit, which leads to efficient algorithms for convex combinations of Dirac masses. Last, we illustrate by numerical experiments the resulting sampling methods that preserve the convex order and their application to approximate Martingale Optimal Transport problems.