Élien Aur

Jacopo Alfonsi

Benjamin Corbetta
email: j.corbetta@zeliade.com

Jourdain

SAMPLING OF ONE-DIMENSIONAL PROBABILITY MEASURES IN THE CONVEX ORDER AND COMPUTATION OF ROBUST OPTION PRICE BOUNDS

Keywords: convex order, martingale optimal transport, robust option price bounds, sampling techniques, linear programming AMS Subject Classification (2010): 91G60, 90C08, 60G42

, the set of probability measures on the real line with a finite first order moment is a complete lattice for the increasing and the decreasing convex orders. For µ and ν in this set, this enables us to define a probability measure µ ∨ ν (resp. µ ∧ ν) greater than µ (resp. smaller than ν) in the convex order. We give efficient algorithms permitting to compute µ ∨ ν and µ ∧ ν (and therefore µI ∨ νJ and µI ∧ νJ ) when µ and ν have finite supports. Last, we illustrate by numerical experiments the resulting sampling methods that preserve the convex order and their application to approximate martingale optimal transport problems and in particular to calculate robust option price bounds.

Introduction and motivations

We start with a practical motivation. Let us consider a trader who sells and hedges exotic options on an asset (S t ) t≥0 . To calculate option prices, the trader typically picks her favourite model, then calibrates this model to vanilla options prices and uses this calibrated model ( S t ) t≥0 to calculate exotic option prices. Then, a natural question is to know what is the range of prices that can be attained by this procedure. To make things simple, we assume zero interest rates and consider a trader who wants to deal an exotic option that pays c(S T 1 , S T 2 ) at time T 2 (with T 1 < T 2 and c : R × R → R) and observes on the market European call option prices at times T 1 and T 2 for all strikes. This amounts to know exactly µ and ν, the respective probability laws of S T 1 and S T 2 under all risk-neutral measures, and we suppose that S T 1 and S T 2 are respectively distributed according to µ and ν (this is achieved by the local volatility model and extensions). Then, all the models that are compatible with the absence of arbitrage are described by the set of martingale couplings

Π M (µ, ν) = π ∈ Π(µ, ν) : ∀x ∈ R, R |y|π Y |X (x, dy) < ∞ and R yπ Y |X (x, dy) = x ,
where π Y |X denotes a Markov kernel such that π(dx, dy) = µ(dx)π Y |X (x, dy). Notice that the joint law of ( S T 1 , S T 2 ) belongs to Π M (µ, ν). The range of prices of the exotic option is given by the set R×R c(x, y)π(dx, dy), π ∈ Π M (µ, ν) and we are interested in calculating the upper and lower bounds of this set, which are usually called the robust option price bounds. This problem is precisely the martingale optimal transport (MOT) problem studied by [START_REF] Beiglböck | Model-independent bounds for option prices-a mass transport approach[END_REF]. They have shown that the upper (resp. lower) bound is also the cheapest (resp. most expensive) price of a superhedging (resp. subhedging) strategy from the dual formulation of the problem. In practice, the trader is interested in calculating these bounds and comparing them to her model price, in order to evaluate the model risk. Unless for particular payoff functions c, she has to use a Monte-Carlo method and approximate the exotic price in her model by 1

I I i=1 c( S i T 1 , S i T 2 )
, where ( S i ) 1≤i≤I are independent copies of S. Let us note that by construction, the empirical measures µ I = 1

I I i=1 δ S i T 1 and ν I = 1 I I i=1 δ S i T 2
approximate µ and ν. The goal of the paper is to give a way to calculate from µ I and ν I the robust price option bounds. Our approach is to consider the corresponding discrete martingale optimal transport problem and to rely on the powerful linear programming library that have been developed. To do so, we need to modify slightly the empirical measures µ I and/or ν I in order to recover the convex order, as we explain now.

We now present the mathematical framework and consider X and Y , two random variables on R, with respective probability distributions µ and ν. The random variable X is said to be smaller than Y in the convex order if E[φ(X)] ≤ E[φ(Y )] for every convex function φ : R → R, provided that both expectations exist. In this case, we use the notation X ≤ cx Y or µ ≤ cx ν, since the convex order only involves the probability distributions. Theorem 8 in [START_REF] Strassen | The existence of probability measures with given marginals[END_REF] ensures that, when R |y|ν(dy) < ∞, Π M (µ, ν) = ∅ ⇐⇒ µ ≤ cx ν. In this work, we consider the approximation of the probability measures µ and ν by probability measures with finite supports µ I = I i=1 p i δ x i and ν J = J j=1 q j δ y j , with I, J ∈ N * , x i , y j ∈ R p i , q j > 0 for any i, j and I i=1 p i = J j=1 q j = 1. Note that we simply have p i = 1/I and q j = 1/J for empirical measures. We consider the approximation of the MOT problem: to minimize (or maximize)

I i=1 J j=1 r ij c(x i , y j ) (1.1)
under the constraints

r ij ≥ 0, I i=1 r ij = q j , J j=1 r ij = p i and J j=1
r ij y j = p i x i .

Thus, we can use a linear programming library to solve this problem numerically. The key issue to run these algorithms is the existence of such matrices (r ij ) 1≤i≤I,1≤j≤J , that amounts to the existence of a martingale coupling between µ I and ν J . Otherwise, the set of matrices satisfying the constraints is void. By Strassen's theorem, this is equivalent to have µ I ≤ cx ν J , which motivates the interest of preserving the convex order when approximating µ and ν.

Up to now, we have dealt with only two measures µ ≤ cx ν. More generally, if µ 1 ≤ cx . . . ≤ cx µ are probability measures, one may be interested in constructing approximations that preserve the convex order. In the same manner, such approximations would enable to tackle numerically multi-marginal martingale optimal transport problems, see again [START_REF] Beiglböck | Model-independent bounds for option prices-a mass transport approach[END_REF].

Up to our knowledge, few studies consider the problem of preserving the convex order while approximating two probability measures. We mention the thesis of [START_REF] Baker | Martingales with specified marginals[END_REF] who proposes the following construction for µ ≤ cx ν that are integrable. Let F µ (x) = µ((-∞, x]) and F ν (x) = ν((-∞, x]) be the cumulative distribution functions and for p ∈ (0, 1), F -1 µ (p) = inf{x ∈ R : F µ (x) ≥ p} and F -1 ν (p) = inf{x ∈ R : F ν (x) ≥ p} their left-continuous and non-decreasing generalized inverses also called quantile functions. For I ≥ 1, setting μI = 1

I I i=1 δ I i I i-1 I F -1 µ (u)du
and νI = 1

I I i=1 δ I i I i-1 I F -1 ν (u)du , (1.2)
one has by Theorem 2.4.11 in [START_REF] Baker | Martingales with specified marginals[END_REF] that μI ≤ cx νI (this result is consequence of the characterization of the convex order in terms of the quantile functions that we recall below in Theorem 2.1 (v)). This construction is easy to implement and can be obviously generalized to the multi-marginal case, but it requires an explicit calculation of the integrals of the quantile functions. Notice that nevertheless, these integrals may be deduced by calibration if one observes for i ∈ {1, . . . , I -1} the price C X i of the European call option with strike K X i = F -1 µ (i/I) written on the asset X as well as the current price s 0 = E[X] of this asset. Since C X i = E[(X -F -1 µ (i/I)) + ] =

1 i I
F -1 µ (p)dp -1 -i I F -1 µ (i/I), we get

I i I i-1 I F -1 µ (p)dp = I(C X i-1 -C X i )+(I -i+1)K X i-1 -(I -i)K X i (with the convention C X I = 0)
for i = 2, . . . , I and I 1 I 0 F -1 µ (p)dp = I(s 0 -C X 1 ) -(I -1)K X 1 . The dual quantization introduced by [START_REF] Pagès | Intrinsic stationarity for vector quantization: foundation of dual quantization[END_REF] gives another way to preserve the convex order in dimension one (see the remark after Proposition 10 in [START_REF] Pagès | Intrinsic stationarity for vector quantization: foundation of dual quantization[END_REF]). Suppose that µ and ν have a bounded support included in [x 1 , x I ] and that x 1 < • • • < x I . The dual quantization of µ is defined by

μI = µ({x 1 })δ x 1 + I-1 i=1 (x i ,x i+1 ] x -x i x i+1 -x i µ(dx)δ x i+1 + (x i ,x i+1 ] x i+1 -x x i+1 -x i µ(dx)δ x i .
This is the law

X = X1 X=x 1 + I-1 i=1 1 X∈(x i ,x i+1 ] x i+1 1 U ≤ X-x i x i+1 -x i + x i 1 U > X-x i x i+1 -x i
, where

U is sampled independently from X according to the uniform distribution on [0, 1]. Let φ : R → R be a convex function. We define φ :

[x 1 , x I ] → R by φ(x) = x-x i x i+1 -x i φ(x i+1 ) + x i+1 -x x i+1 -x i φ(x i ) for x ∈ [x i , x i+1 ], i = 1, . . . , I -1.
This is a convex function, and we have by )]. Note that both methods have serious drawbacks to be used in practice. Baker's construction requires to calculate explicitly the integrals of the quantile functions, which is often not possible. Dual quantization assumes probability distributions with bounded support, which is not satisfied by most of the asset models. On the other side, there are many recent works that deals with robust option price bounds for specific options. Without being exhaustive, [START_REF] Hobson | Robust bounds for forward start options[END_REF], [START_REF] Hobson | Robust price bounds for the forward starting straddle[END_REF] characterize the bounds for forward start options. Henry-Labordère and Touzi (2016) describe explicitly the bounds for an abstract class of payoff functions c satisfying the so-called Spence-Mirrlees condition, for which the left and right curtain couplings introduced by [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF] are optimal. These works crucially exploits particular properties of the payoff functions. In contrast, the numerical method that we propose in this paper to approximate robust option price bounds works for general payoff functions.

using µ ≤ cx ν that E[φ( X)] = E[ φ(X)] ≤ E[ φ(Y )] = E[φ( Ŷ
The paper is structured as follows. In Section 2, we recall the characterization of the convex, increasing convex and decreasing convex orders in dimension one, and give a way to test if two probability measures with finite support are in the convex order. Then, in Section 3, we shift the approximating measures µ I and ν J so that they have the same mean and show that they are then in the convex order for I, J large enough, under suitable assumptions. However, this method is not reliable in practice since the dimension I × J of the linear programming problem (1.1) can then be too large for this problem to be solved quickly. It is thus crucial to use non asymptotic methods that gives approximating measures in the convex order for any I, J ∈ N * . In Section 4, we propose such a method based on the following result stated in [START_REF] Kertz | Stochastic and convex orders and lattices of probability measures, with a martingale interpretation[END_REF]: the set of probability measures with a finite first order moment is a lattice for the increasing and the decreasing convex orders. Given µ and ν in this set, we define µ ∨ ν (resp. µ ∧ ν) as the supremum (resp. infimum) of µ and ν for the decreasing convex order when R xµ(dx) ≤ R xν(dx) and for the increasing convex order otherwise. This way, µ ∨ ν (resp. µ ∧ ν) is greater then µ (resp. smaller than ν) in the convex order. In Section 5, we give efficient algorithms permitting to compute µ ∨ ν and µ ∧ ν (and therefore µ I ∨ ν J and µ I ∧ ν J ) when µ and ν are convex combinations of Dirac masses. We prove that, when µ ≤ cx ν, µ I ∨ ν J and µ I ∧ ν J respectively converge weakly to ν and µ as I, J → ∞, which shows that either µ I and µ I ∨ ν J or µ I ∧ ν J and ν J respectively approximate µ and ν and are in the convex order. Section 6 presents the numerical approximation of different MOT problems. We first consider academic examples where the MOT is known explicitly. Then, we address some practical examples with two or three marginals, when the marginal laws are the ones given by the Black and Scholes model.

This paper only deals with the one-dimensional problem. Of course, the same problem occurs in higher dimension when one consider d ≥ 2 assets instead of one. Thus, for two probability measures µ and ν on R d such that µ ≤ cx ν, and for discrete probability measures µ I and ν J that respectively approximate µ and ν, we want to slightly modify either µ I or ν J to recover the convex order between the approximating measures. Unfortunately, we can no longer rely on the lattice structure to construct such a modification since it has been shown by [START_REF] Müller | Stochastic order relations and lattices of probability measures[END_REF] that the set of probability measures with a given mean is not a lattice for the convex order when d ≥ 2. In a companion paper [START_REF] Alfonsi | Sampling of probability measures in the convex order by Wasserstein projection and approximation of martingale optimal transport problems[END_REF], we propose a suitable modification of µ I that is obtained by the mean of a projection for the Wasserstein distance and that can be efficiently computed.

Characterization of the convex, decreasing convex and increasing orders

Let P(R) denote the set of probability measures on the real line. We define

P -(R) = {µ ∈ P(R) : R x -µ(dx) < ∞}, P + (R) = {µ ∈ P(R) : R x + µ(dx) < ∞}, P 1 (R) = {µ ∈ P(R) : R |x|µ(dx) < ∞} = P -(R) ∩ P + (R).
For µ ∈ P(R), we set

F µ (x) = µ((-∞, x]), F µ (x-) = lim y→x-F µ (y) = µ((-∞, x)). For t ∈ R, we define ϕ µ (t) = t -∞ F µ (x)dx. Fubini's theorem ensures that for t ∈ R, ϕ µ (t) = R 1 {x≤t} R 1 {y≤x} µ(dy) = R (t -y) + µ(dy) (2.1)
where the right-hand side is finite iff µ ∈ P -(R). In a symmetric way, we define Fµ

(x) = µ([x, +∞)), φµ (t) = +∞ t
Fµ (x)dx and notice that

φµ (t) = +∞ t µ((x, +∞))dx = R (y -t) + µ(dy) (2.2)
where the right-hand side is finite iff µ ∈ P + (R). Let us note that µ((x, +∞)) = µ([x, +∞)) dx-a.e., and equation (2.2) is written with the open interval to have the same convention as [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF]. Denoting by μ the image of µ by x → -x, one has

∀x ∈ R, F μ(x) = μ((-∞, x]] = µ([-x, +∞)) = Fµ (-x) and thus ∀t ∈ R, φµ (t) = ϕ μ(-t).
(2.3) Last, we define the function

π µ (t) = R |t -x|µ(dx) = ϕ µ (t) + φµ (t),
that is usually called the potential of µ. It is finite when

µ ∈ P 1 (R). Since t - R xµ(dx) = ϕ µ (t) -φµ (t) = 2ϕ µ (t) -π µ (t) = π µ (t) -2 φµ (t), we have π µ (t) = 2ϕ µ (t) -t + R xµ(dx) = 2 φµ (t) + t - R xµ(dx).
One important particularity of the dimension 1 is the following result (see e.g. Theorems 3.A.1., 3.A.2. and 3.A.5. in Shaked and Shanthikumar (2007)).

Theorem 2.1. Let µ, ν ∈ P 1 (R). The following conditions are equivalent:

(i) µ ≤ cx ν, (ii) R xµ(dx) = R xν(dx) and ∀t ∈ R, π µ (t) ≤ π ν (t), (iii) R xµ(dx) = R xν(dx) and ∀t ∈ R, ϕ µ (t) ≤ ϕ ν (t), (iv) R xµ(dx) = R xν(dx) and ∀t ∈ R, φµ (t) ≤ φν (t), (v) 1 0 F -1 µ (p)dp = 1 0 F -1 ν (p)dp and ∀q ∈ (0, 1), 1 q F -1 µ (p)dp ≤ 1 q F -1 ν (p)dp.
Thus, to check if two probability measures are in the convex order, it is sufficient to focus on the the family of convex functions φ(x) = |t -x|, for t ∈ R. Thanks to this result, we obtain an interesting corollary that gives a necessary and sufficient condition for two probability measures with finite support to be in the convex order.

Corollary 2.2. Let µ = I i=1 p i δ x i and ν = J j=1 q j δ y j be two probability measures on R. Without loss of generality, we assume that

x 1 < • • • < x I , y 1 < • • • < y J and p 1 p I q 1 q J > 0.
Then we have µ ≤ cx ν if, and only if (i) y 1 ≤ x 1 and y J ≥ x I , (ii) for all j such that x 1 ≤ y j ≤ x I , π µ (y j ) ≤ π ν (y j ), (iii) I i=1 p i x i = J j=1 q j y j . By using (iii) and the link between π µ and ϕ µ , this corollary is still true if we replace (ii) by one of the following conditions:

(ii ) for all j such that x 1 ≤ y j ≤ x I , ϕ µ (y j ) ≤ ϕ ν (y j ), (ii ) for all j such that x 1 ≤ y j ≤ x I , φµ (y j ) ≤ φν (y j ).

Proof. The necessary condition is obvious. Let us check the sufficient condition. Since y 1 = min 1≤j≤J y j and y J = max 1≤j≤J y j , we obtain

∀t ≤ y 1 , π ν (t) = J j=1 q j y j -t and ∀t ≥ y M , π ν (t) = t - J j=1 q j y j . Similarly, we have π µ (t) = I i=1 p i x i -t for t ≤ x 1 and π µ (t) = t -I i=1 p i x i for t ≥ x I . By (i) and (iii), we get π µ (t) = π ν (t) for t ∈ (-∞, y 1 ] ∪ [y M , +∞) and π µ (t) = t - I i=1 p i x i = J j=1 q j (t -y j ) ≤ J j=1 q j |t -y j | = π ν (t)
for t ∈ (-∞, x 1 ] ∪ [x M , +∞) by Jensen's inequality. Moreover, since π µ is convex and π ν is affine on [y j ; y j+1 ] for every j = 1, . . . , M -1, we get from (ii) that π µ (t) ≤ π ν (t) for every t ∈ [y j ; y j+1 ] for every j = 1, . . . , M -1. Thus, we have π µ (t) ≤ π ν (t) for all t ∈ R.

In dimension one, the increasing and decreasing convex orders may be defined as follows.

Definition 2.3. For µ, ν ∈ P(R), we say that µ is smaller than ν in the increasing (resp. decreasing) convex order and denote µ ≤ icx ν (resp. µ ≤ dcx ν) if R φ(x)µ(dx) ≤ R φ(x)ν(dx) for each increasing (resp. decreasing) convex function φ : R → R such that the integrals make sense.

Lemmas 2.2 and 2.5 of [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF] give the following characterization of the increasing and decreasing convex orders.

Theorem 2.4. Let µ, ν ∈ P + (R) (resp. P -(R)). The following statements are equivalent :

(i) µ ≤ icx ν (resp. µ ≤ dcx ν), (ii) ∀t ∈ R, φµ (t) ≤ φν (t) (resp. ϕ µ (t) ≤ ϕ ν (t)), (iii) ∀q ∈ [0, 1], 1 q F -1 µ (p)dp ≤ 1 q F -1 ν (p)dp (resp. q 0 F -1 µ (p)dp ≥ q 0 F -1 ν (p)dp), (iv) μ ≤ dcx ν (resp. μ ≤ icx ν).
Note that the equivalence between (ii) and (iii) is a direct consequence that ϕ µ is the Legendre transform of the convex function q → q 0 F -1 µ (p)dp, and conversely, see e.g. Lemma A.22 of [START_REF] Föllmer | Stochastic finance, An introduction in discrete time[END_REF].

An asymptotic approach

Now, we turn to our practical problem. We assume that µ ≤ cx ν with µ, ν ∈ P 1 (R). We consider two i.i.d. samples X 1 , . . . , X I with distribution µ and Y 1 , . . . , Y J with distribution ν, and we set μI = 1

I I i=1 δ X i and νJ = 1 J J j=1 δ X j .
At least when µ and ν have densities, the empirical means XI = 1 I I i=1 X i and ȲJ = 1 J J j=1 Y j are almost surely distinct, and we cannot have μI ≤ cx νJ . However, we know that when I, J → +∞, both empirical means converge almost surely to R xµ(dx). A natural idea is then to slightly modify μI , or νJ , or both empirical measures so that they both have the same mean. The hope is then that for I and J large enough, the modified empirical measures would be naturally in the convex order. In this paragraph, we will take

µ I,J = 1 I I i=1 δ X i + ȲJ -XI , (3.1)
and keep νJ . Other choices are of course possible and some of them are presented in Subsection 4.2.

For a probability measure µ on the real line, let F -1 µ (0+) and F -1 µ (1-) respectively denote the left-hand and the right-hand limits of the function F -1 µ as p → 0 and p → 1. We have the following result.

Proposition 3.1. Let µ, ν ∈ P 1 (R) be such that µ ≤ cx ν and F -1 ν (0+) < F -1 µ (0+) ≤ F -1 µ (1-) < F -1 ν (1-).
The probability measure µ I,J defined by (3.1) (resp. νJ ) converges weakly to µ (resp. ν), almost surely when I, J → +∞. Suppose moreover that there exists ε > 0 such that for every t ∈

[F -1 µ (0+), F -1 µ (1-)] t -∞ (F ν (x) -F µ (x))dx ≥ ε. (3.2)
Then, almost surely, there exists M such that for all I, J ≥ M , µ I,J ≤ cx νJ .

(3.3)

Proof. Let f : R → R be a continuous bounded function. The strong law of large number gives the (almost sure) weak convergence of μI (resp. νJ ) towards µ (resp. ν), as well as the almost sure convergence of ȲJ -XI towards 0 as I, J → ∞. Now, we use that f is uniformly continuous on

[F -1 µ (0+) -1, F -1 µ (1-) + 1] to obtain that 1 I I i=1 f (X i + ȲJ -XI ) → I,J→+∞ f (x)µ(dx)
almost surely. Thus, µ I,J converges weakly to µ, almost surely. We now turn to the convex order. Let δ > 0 be such that

∀t ∈ [F -1 µ (0+) -δ, F -1 µ (1-) + δ], t -∞ (F ν (x) -F µ (x))dx ≥ ε/2. We take t ∈ (-∞, F -1 µ (0+) -δ) such that t -∞ F ν (x)dx < ε 6 , which implies that ∀t ∈ [F -1 µ (0+) -δ, F -1 µ (1-) + δ], t t (F ν (x) -F µ (x))dx ≥ ε 3 .
From the weak convergence of νJ (resp. µ I,J ) to ν (resp. µ), we obtain that dx-a.e., F νJ (x) → F ν (x) (resp. F µ I,J (x) → F µ (x)). Then, the dominated convergence theorem gives the pointwise convergence of t F µ (x)dx). Since these functions are nondecreasing and the limit is continuous with respect to t, we obtain the uniform convergence on each compact set, and deduce that, almost surely, there exists M such that

∀I, J ≥ M, ∀t ∈ [F -1 µ (0+) -δ, F -1 µ (1-) + δ], t t (F νJ (x) -F µ I,J (x))dx ≥ ε 4 .
Taking also M large enough so that for I, J ≥ M , | ȲJ -XI | < δ, we have F µ I,J (x) = 0 and

F νJ (x) -F µ I,J (x) ≥ 0 for x ≤ F -1 µ (0+) -δ and deduce that ∀t ∈ (-∞, F -1 µ (1-) + δ], t -∞ (F νJ (x) -F µ I,J (x))dx ≥ 0.
Last, for I, J ≥ M , F µ I,J (x) = 1 for x ≥ F -1 µ (1-)+δ and since lim t→+∞ ϕ νJ (t)-ϕ µ I,J (t) = 0, we get that for t

≥ F -1 µ (1-) + δ, t -∞ (F νJ (x) -F µ I,J (x))dx = +∞ t
(1 -F νJ (x))dx ≥ 0. We conclude by using Theorem 2.1 (iii).

We now give a sufficient condition for (3.2). Lemma 3.2. Let µ, ν be two probability measures on the real line with µ ≤ cx ν and

F -1 ν (0+) < F -1 µ (0+) ≤ F -1 µ (1-) < F -1 ν (1-). Suppose that ∃x 0 ∈ R, ∀x ∈ (-∞, x 0 ], F ν (x) -F µ (x) ≥ 0 and ∀x ∈ [x 0 , +∞), F ν (x) -F µ (x) ≤ 0.
Then, (3.2) holds.

Proof. We know that

F ν (x) -F µ (x) = 0 when x / ∈ [F -1 ν (0+), F -1 ν (1-)), F ν (x) -F µ (x) = F ν (x) > 0 when x ∈ (F -1 ν (0+), F -1 µ (0+)), F ν (x) -F µ (x) = F ν (x) -1 < 0 when x ∈ [F -1 µ (1-), F -1 ν (1-))
, and by the equality of the means we know

F -1 ν (1-) -∞
(F ν (x)-F µ (x))dx = 0. In particular, we necessarily have

x 0 ∈ [F -1 µ (0+), F -1 µ (1-)) and we get that t → t -∞ (F ν (x) -F µ (x))dx is nondecreasing on (∞, x 0 ], nonincreasing on [x 0 , ∞), increas- ing on [F -1 ν (0+), F -1 µ (0+)], decreasing on [F -1 µ (1-), F -1 ν (1-)]. Thus F -1 µ (0+) -∞ (F ν (x) - F µ (x))dx > 0, F -1 µ (1-) -∞ (F ν (x) -F µ (x))dx > 0 and (3.2) holds for ε = min F -1 µ (0+) -∞ (F ν (x) -F µ (x))dx > 0, F -1 µ (1-) -∞ (F ν (x) -F µ (x))dx .
In particular, we see that under the assumptions of Lemma 3.2, ϕ ν (t) > ϕ µ (t) for any t ∈ (F -1 ν (0+), F -1 ν (1-)), which means that (µ, ν) is irreducible, see Definition A.3 of [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF].

Proposition 3.1 gives a framework under which the empirical measures, up to a modification that equalizes their means, are asymptotically in the convex order. The assumptions of Proposition 3.1 are quite restrictive and could be refined. Nonetheless, it is not very convenient in practice to reach the convex order only asymptotically. If we have in mind to solve the discrete MOT problem (1.1), which is a linear programming problem in dimension I × J. To fix the ideas, if we need I, J 10 3 for the modified empirical measures to be in the convex order, the resolution of the discrete MOT problem (1.1) is already too greedy in time and memory. Thus, it would be much more convenient if we could guarantee that the approximating measure are in the convex order. This is why we prefer to focus on non asymptotic methods that ensures the convex order for any I, J.

4.

A non asymptotic approach based on the lattice structure 4.1. The lattice structure for the increasing and decreasing convex orders. It has been observed by Kertz and Rösler that for any a ∈ R,

η ∈ P 1 (R) such that R xη(dx) = a
is a lattice for the convex order (Proposition 1.6 of [START_REF] Kertz | Stochastic and convex orders and lattices of probability measures, with a martingale interpretation[END_REF]) and even a complete lattice (see top of p162 of [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF]). By Proposition 4.5 of [START_REF] Müller | Stochastic order relations and lattices of probability measures[END_REF], this property is no longer true in higher dimension. Thus, if µ, ν ∈ P 1 (R) share the same expectation, there is a unique probability measure µ [START_REF] Kertz | Stochastic and convex orders and lattices of probability measures, with a martingale interpretation[END_REF], φµ∧cxν is the greatest convex function below φµ ∧ φν and φµ∨cxν = φµ ∨ φν .

∧ cx ν ∈ P 1 (R) (resp. µ ∨ cx ν ∈ P 1 (R)) such that µ ∧ cx ν ≤ cx µ, µ ∧ cx ν ≤ cx ν (resp. µ ≤ cx µ ∨ cx ν, µ ≤ cx µ ∨ cx ν) and η ≤ cx µ ∧ cx ν for all η ∈ P 1 (R) such that η ≤ cx µ and η ≤ cx ν (resp. µ ∨ cx ν ≤ cx η for all η ∈ P 1 (R) such that µ ≤ cx η and ν ≤ cx η). By Proposition 1.6 of
To deal with the case when µ, ν

∈ P 1 (R) are such that R xµ(dx) = R xν(dx)
, we are going to use the complete lattice property of P + (R) (resp. P -(R)) endowed with the increasing (resp. decreasing) convex order stated in Theorem 3.4 of [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF] (resp. Theorem 3.7 (a) of [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF]). For µ, ν [START_REF] Kertz | Stochastic and convex orders and lattices of probability measures, with a martingale interpretation[END_REF], for µ, ν ∈ P + (R), φµ∧ icx ν is the greatest convex function below φµ ∧ φν and φµ∨ icx ν = φµ ∨ φν . By Theorem 3.7 (a) of [START_REF] Kertz | Complete lattices of probability measures with applications to martingale theory[END_REF], for µ, ν ∈ P -(R), µ ∧ dcx ν = μ ∧ icx ν and µ ∨ dcx ν = μ ∨ icx ν, so that by (2.3), ϕ µ∧ dcx ν is the greatest convex function below ϕ µ ∧ ϕ ν and ϕ µ∨ dcx ν = ϕ µ ∨ ϕ ν . Let us now compute the expectations of

∈ P + (R) (resp. P -(R)), let µ ∧ icx ν, µ ∨ icx ν ∈ P + (R) (resp. µ ∧ dcx ν, µ ∨ dcx ν ∈ P -(R)) satisfy the above properties of µ ∧ cx ν, µ ∨ cx ν but with ≤ cx replaced by ≤ icx (resp. ≤ dcx ). By Proposition 1.4 of
µ ∧ icx ν, µ ∨ icx ν, µ ∧ dcx ν and µ ∨ dcx ν when µ, ν ∈ P 1 (R). Lemma 4.1. Let µ, ν ∈ P 1 (R). Then µ ∧ icx ν, µ ∨ icx ν, µ ∧ dcx ν, µ ∨ dcx ν ∈ P 1 (R) and R xµ ∨ icx ν(dx) = R xµ(dx) ∨ R xν(dx), R xµ ∨ dcx ν(dx) = R xµ(dx) ∧ R xν(dx), R xµ ∧ icx ν(dx) = R xµ(dx) ∧ R xν(dx), R xµ ∧ dcx ν(dx) = R xµ(dx) ∨ R xν(dx). Proof. By (2.1), we have t -ϕ µ (t) = R xµ(dx) - R (x -t) + µ(dx)
, where the last term tends to 0 as t → +∞ by Lebesgue theorem. Since, in the same way, t -ϕ ν (t) tends to

R xν(dx), we deduce that t -ϕ µ ∨ ϕ ν (t) tends to R xµ(dx) ∧ R xν(dx). On the other hand, by (2.1) with µ replaced by µ ∨ dcx ν, t -ϕ µ ∨ ϕ ν (t) = t - R (t -x) + µ ∨ dcx ν(dx) = R (t ∧ x)µ ∨ dcx ν(dx) For t ≥ 0, the right-hand side is equal to - R x -µ ∨ dcx ν(dx) + R (x + ∧ t)µ ∨ dcx ν(dx), where the second term converges to R x + µ ∨ dcx ν(dx) as t → ∞ by monotone convergence. Therefore R |x|µ ∨ dcx ν(dx) < ∞ and R xµ ∨ dcx ν(dx) = R xµ(dx) ∧ R xν(dx).
The fact that µ ∨ icx ν ∈ P 1 (R) and the formula giving its expectation are deduced from the equality µ ∨ icx ν = μ ∨ dcx ν. For µ ∧ icx ν and µ ∧ dcx ν, this is an easy consequence of the inverse transform sampling method combined with the more convenient characterization of these probability measures based on the quantile functions stated in Lemma 4.2 just below.

Lemma 4.2.

For µ, ν ∈ P + (R), ∀q ∈ [0, 1], 1 q F -1 µ∧ icx ν (p)dp = 1 q F -1 µ (p)dp ∧ 1 q F -1 ν (p)dp,
and for µ, ν

∈ P -(R), ∀q ∈ [0, 1], q 0 F -1 µ∧ dcx ν (p)dp = q 0 F -1 µ (p)dp ∨ q 0 F -1 ν (p)dp.
Moreover, for all p ∈ (0, 1),

F -1 µ∧ icx ν (p), F -1 µ∧ dcx ν (p) ∈ {F -1 µ (p), F -1 ν (p)}. Proof.
The first statement is deduced from the second one and the equality µ ∧ dcx ν = μ ∧ icx ν. For f : R → (-∞, +∞], we define f * (y) = sup x∈R xy -f (x) the Fenchel-Legendre transform of f and recall that f * * := (f * ) * is the greatest convex function below f . Thus, we have ϕ µ∧ dcx ν = (ϕ µ ∧ ϕ ν ) * * . This gives ϕ * µ∧ dcx ν = (ϕ µ ∧ ϕ ν ) * . By using a standard property of the Fenchel-Legendre transform, we have (ϕ µ ∧ ϕ ν ) * = ϕ * µ ∨ ϕ * ν , and we deduce that [START_REF] Föllmer | Stochastic finance, An introduction in discrete time[END_REF].

q 0 F -1 µ∧ dcx ν (p)dp = q 0 F -1 µ (p)dp ∨ q 0 F -1 ν (p)dp for all q ∈ [0, 1] by Lemma A.23 of
Let us check now that

F -1 µ∧ icx ν (p) ∈ {F -1 µ (p), F -1 ν (p)} for p ∈ (0, 1). The function q → g(q) := 1 q F -1 µ (p)dp - 1 q F -1 ν (p)dp
is locally bounded with a locally bounded derivative on (0, 1) so that the distribution derivative of g + (q) is equal to 1 {g(q)>0} (F -1 ν (q) -F -1 µ (q)) and the one of

q → 1 q F -1 µ∧ icx ν (p)dp = 1 q F -1 µ (p)dp -g + (q)
is equal to -1 {g(q)≤0} F -1 µ (q) -1 {g(q)>0} F -1 ν (q). Therefore dq a.e. on (0, 1), F -1 µ∧ icx ν (q) = 1 {g(q)≤0} F -1 µ (q) + 1 {g(q)>0} F -1 ν (q). Let now p ∈ (0, 1). There is a sequence (q n ) n∈N of elements in (0, p) such that lim n→∞ q n = p and F -1

µ∧ icx ν (q n ) ∈ {F -1 µ (q n ), F -1 ν (q n )} for all n ∈ N. Either F -1 µ∧ icx ν (q n ) = F -1 µ (q n
) for infinitely many n and, by left-continuity of the quantile functions, F -1

µ∧ icx ν (p) = F -1 µ (p) or F -1 µ∧ icx ν (q n ) = F -1 ν (q n ) for infinitely many n and F -1 µ∧ icx ν (p) = F -1 ν (p). Similarly, we show that F -1 µ∧ dcx ν (p) ∈ {F -1 µ (p), F -1 ν (p)} by considering the derivative of q → q 0 F -1 µ (p)dp - q 0 F -1 ν (p)dp + + q 0 F -1 ν (p)dp. Definition 4.3. For µ, ν ∈ P 1 (R), we define µ ∧ ν, µ ∨ ν ∈ P 1 (R) by µ ∧ ν = 1 { R xµ(dx)≤ R xν(dx)} µ ∧ dcx ν + 1 { R xµ(dx)> R xν(dx)} µ ∧ icx ν, (4.1) µ ∨ ν = 1 { R xµ(dx)≤ R xν(dx)} µ ∨ dcx ν + 1 { R xµ(dx)> R xν(dx)} µ ∨ icx ν, (4.2) By Lemma 4.1, R xµ∧ν(dx) = R xν(dx) and R xµ∨ν(dx) = R xµ(dx)
. By Lemma 1.5 of [START_REF] Kertz | Stochastic and convex orders and lattices of probability measures, with a martingale interpretation[END_REF], for fixed a ∈ R, the convex, increasing convex and decreasing convex orders coincide on {η ∈ P

1 (R) such that R xη(dx) = a}. Therefore µ ∧ ν ≤ cx ν and µ ≤ cx µ ∨ ν, (4.3) and, when R xµ(dx) = R xν(dx), µ ∧ ν = µ ∧ cx ν and µ ∨ ν = µ ∨ cx ν. Remark 4.4. Let µ, ν ∈ P 1 (R). Then ν ≤ cx µ∨ν ⇐⇒ R xν(dx) = R xµ(dx). Let ν be the image of ν by t → t+ R xµ(dx)- R xν(dx). One has ϕ ν (t) = ϕ ν t + R xν(dx) - R xµ(dx) and φ ν (t) = φν t + R xν(dx) - R xµ(dx) for all t ∈ R. When R xµ(dx) ≤ R xν(dx), then ϕ ν ≤ ϕ ν and µ ∨ ν = µ ∨ dcx ν so that ϕ µ∨ν = ϕ µ ∨ ϕ ν ≤ ϕ µ ∨ ϕ ν = ϕ µ∨ ν . When R xµ(dx) > R xν(dx), then φν ≤ φ ν and µ ∨ ν = µ ∨ icx ν so that φµ∨ν = φµ ∨ φν ≤ φµ ∨ φ ν = φµ∨ ν . Therefore, in both cases, µ ≤ cx µ ∨ ν ≤ cx µ ∨ ν.
In the same way, µ∧ν

≤ cx µ if and only if R xµ(dx) = R xν(dx). Let µ denote the image of µ by t → t + R xν(dx) - R xµ(dx). One has F -1 µ (p) = F -1 µ (p) + R xν(dx) - R xµ(dx) for p ∈ (0, 1). When R xµ(dx) ≤ R xν(dx), then F -1 µ ≤ F -1 µ and µ ∧ ν = µ ∧ dcx ν, so that, by Lemma 4.2, q 0 F -1 µ∧ν (p)dp = q 0 F -1 µ (p)dp ∨ q 0 F -1 ν (p)dp ≤ q 0 F -1 µ (p)dp ∨ q 0 F -1 ν (p)dp = q 0 F -1 µ∧ν (p)dp for q ∈ [0, 1]. When R xµ(dx) > R xν(dx), then F -1 µ > F -1 µ and µ ∧ ν = µ ∧ icx ν, so that, again by Lemma 4.2, 1 q F -1 µ∧ν (p)dp = 1 q F -1 µ (p)dp ∧ 1 q F -1 ν (p)dp ≥ 1 q F -1 µ (p)dp ∧ 1 q F -1 ν (p)dp = 1 q F -1 µ∧ν (p)dp for q ∈ [0, 1]. With Theorem 2.4 (iii), we deduce that in both cases, µ ∧ ν ≤ cx µ ∧ ν ≤ cx ν.
Let us now check that the diameter of the set {µ, ν, µ∨ dcx ν, µ∨ icx ν, µ∧ dcx ν, µ∧ icx ν} (resp.

{µ, ν, µ ∨ dcx ν, µ ∧ dcx ν}, {µ, ν, µ ∨ icx ν, µ ∧ icx ν}) when µ, ν ∈ P 1 (R) (resp. P -(R), P + (R)))
in Wasserstein distance is equal to the Wasserstein distance between µ and ν.

Lemma 4.5. For η ∈ {µ, ν, µ∨ dcx ν, µ∨ icx ν, µ∧ dcx ν, µ∧ icx ν}, ∀p ∈ (0, 1), F -1 µ (p)∧F -1 ν (p) ≤ F -1 η (p) ≤ F -1 µ (p) ∨ F -1 ν (p). Let ≥ 1.
Since, by Proposition 2.17 of [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF], the -Wasserstein distance between two probability measures η 1 , η 2 ∈ P(R) is given by

W (η 1 , η 2 ) = 1 0 |F -1 η 1 (p) -F -1 η 2 (p)| dp 1/
, one easily deduces that for all η 1 , η 2 ∈ {µ, ν, µ

∨ dcx ν, µ ∨ icx ν, µ ∧ dcx ν, µ ∧ icx ν}, W (η 1 , η 2 ) ≤ 1 0 (F -1 µ (p) ∨ F -1 ν (p) -F -1 µ (p) ∧ F -1 ν (p)) dp 1/ = W (µ, ν).
Proof. For η ∈ {µ ∧ icx ν, µ ∧ dcx ν} this is a consequence of the last statement in Lemma 4.2 and the left-continuity of the quantile functions. Let

η ∈ {µ ∨ icx ν, µ ∨ dcx ν}. Since for all t ∈ R, t -∞ F µ∧ dcx ν (x)dx = t -∞ F µ (x)dx ∨ t -∞ F ν (x)dx and +∞ t Fµ∧ icx ν (x)dx = +∞ t Fµ (x)dx ∨ +∞ t
Fν (x)dx a reasoning analogous to the proof of this last statement ensures that ∀x ∈ R,

F η (x) ∈ {F µ (x), F ν (x)}. For p ∈ (0, 1), since {x ∈ R : F µ (x) ≥ p} ∩ {x ∈ R : F ν (x) ≥ p} = {x ∈ R : F µ (x) ∧ F ν (x) ≥ p} ⊂ {x ∈ R : F η (x) ≥ p}, {x ∈ R : F µ (x) ≥ p} ∪ {x ∈ R : F ν (x) ≥ p} = {x ∈ R : F µ (x) ∨ F ν (x) ≥ p} ⊃ {x ∈ R : F η (x) ≥ p}, we get F -1 µ (p) ∧ F -1 ν (p) ≤ F -1 η (p) ≤ F -1 µ (p) ∨ F -1 ν (p). 4.2.
Approximations in convex order. Let µ, ν ∈ P 1 (R) be two probability measures such that µ ≤ cx ν. We want to construct a couple of measures with finite supports and in convex order which approximate (µ, ν). We first generate ( x I i ) 1≤i≤I and ( y J j ) 1≤j≤J such that ( 1

I I i=1 δ x I i , 1 J J j=1 δ y J j
) approximates (µ, ν) in one of the following ways :

1. ( x I 1 , . . . , x I I ) = (X 1 , . . . , X I ) and ( y J J , . . . , y

J J ) = (Y 1 , . . . , Y J ) where ((X i , Y i )) i≥1 are random vectors i.i.d. according to µ ⊗ ν, 2. x I i = F -1 µ 2i-1 2I
for i ∈ {1, . . . , I} and y

J j = F -1 ν 2j-1 2J
for j ∈ {1, . . . , J}.

We may want either to simply keep these samples with, in general, distinct empirical means by setting a. (x I 1 , . . . , x I I ) = ( x I 1 , . . . , x I I ) and (y J 1 , . . . , y J J ) = ( y J 1 , . . . , y J J ), or to modify the approximate measures to ensure that they have the same expectation. We then propose one of the following modifications b. set (x I 1 , . . . ,

x I I ) = ( x I 1 + R xµ(dx) -1 I I i=1 x I i , . . . , x I I + R xµ(dx) -1 I I i=1 x I i ) and (y J 1 , . . . , y J J ) = ( y J 1 + R yν(dy) -1 J J j=1 y J j , . . . , y J J + R yν(dy) -1 J J j=1 y J j ). c. set (y J 1 , . . . , y J J ) = ( y J 1 + 1 I I i=1 x I i -1 J J k=1 y J k , . . . , y J J + 1 I I i=1 x I i -1 J J k=1 y J k ) and (x I 1 , . . . , x I I ) = ( x I 1 , . . . , x I I ), d. compute the estimators of the variances σ 2 x = 1 I-1 I i=1 ( x I i ) 2 -1 I(I-1) ( I i=1 x I i ) 2 and σ 2 y = 1 J-1 J j=1 ( y J j ) 2 -1 J(J-1) ( J j=1 y J j ) 2 and add J σ 2 x I σ 2 y +J σ 2 x 1 J J j=1 y J j -1 I I i=1 x I i to each x I i to obtain x I i and I σ 2 y I σ 2 y +J σ 2 x 1 I I i=1 x I i -1 J J j=1 y J j to each y J j to obtain y J j .
The first modification necessitates the knowledge of R xµ(dx) = R yν(dy) which is generally the case in financial applications, since this is the current price of the underlying asset. The second modification seems suitable for the deterministic initialization of points with I = J, since formally, because of the inequality µ ≤ cx ν, the tails of µ should be lighter than the tails of ν enabling a quicker convergence of 1

I I i=1 F -1 µ 2i-1 2I than of 1 I I i=1 F -1 ν 2i-1 2I
to the common expectation of µ and ν as I → ∞. When R y 2 ν(dy) < ∞, the convex combination of 1

I I i=1 X i and 1 J J j=1 Y j with minimal vari- ance is Var(Y 1 ) IVar(Y 1 ) + JVar(X 1 ) I i=1 X i + Var(X 1 ) IVar(Y 1 ) + JVar(X 1 ) J j=1 Y j
and is approximated by the common mean

σ 2 y I σ 2 y +J σ 2 x I i=1 x I i + σ 2 x I σ 2 y +J σ 2
x J j=1 y J j of the two samples (x I 1 , . . . , x I I ) and (y J 1 , . . . , y J J ) for the random initialization of points. We then set

µ I = 1 I I i=1 δ x I i and ν J = 1 J J j=1 δ y J j ,
and consider either (µ I , µ I ∨ ν J ) or (µ I ∧ ν J , ν J ) to approximate (µ, ν) in the convex order. The next proposition shows without giving any speed of convergence that the measures µ and ν are then well approximated asymptotically. Notice that, according to Remark 4.4, the measure µ I ∨ ν I obtained with the second modification dominates in the convex order the one obtained without modification.

Proposition 4.6. Let µ and ν be two probability measures on R such that µ ≤ cx ν and R |x|ν(dx) < ∞. As I, J → ∞, µ I and µ I ∨ ν J (resp. µ I ∧ ν J and ν J ) converges a.s. weakly to µ and ν.

From this proposition, it seems then natural to approximate the continuous MOT inf R×R c(x, y)π(dx, dy), π ∈ Π M (µ, ν) by the minimum of the discrete MOT (1.1) for the measures µ I ∧ ν J and ν J when I, J are large. We will do this with I = J in Section 6 on numerical examples. However, up to our knowledge, no convergence result is known in the literature for the MOT cost between probability measures μ and ν such that μ ≤ cx ν, in the limit μ → µ and ν → ν for the weak sense. We intend to address this important question in future research.

To prove Proposition 4.6, we need the following lemma to deal with the deterministic choice of points.

Lemma 4.7. Let µ be such that R |x|µ(dx) < ∞. As I → ∞, 1 I I i=1 δ F -1 µ ( 2i-1 2I ) converges weakly to µ, 1 I I i=1 F -1 µ 2i-1 2I
converges to R xµ(dx) and for all t ∈ R, 1

I I i=1 (t - F -1 µ 2i-1 2I ) + converges to R (t -x) + µ(dx).
Proof. The weak convergence follows from the continuity of (0, 1) p → F -1 µ (p) outside of a set with zero Lebesgue measure and the weak convergence of 1

I I i=1 δ 2i-1 2I
to the Lebesgue measure on (0, 1). We get by the weak convergence lim inf

I→∞ 1 I I i=1 (F -1 µ ((2i -1)/2I)) -≥ 1 0 (F -1 µ (p)) -dp.
From the inequality

1 I I i=1 (F -1 µ ((2i -1)/2I)) -≤ 1 2I 0 (F -1 µ (p)) -+ 1-1 2I 0 (F -1 µ (p)) -dp,
and since

1 0 |F -1 µ (p))|dp < ∞, we deduce that lim I→∞ 1 I I i=1 (F -1 µ ((2i -1)/2I)) -= 1 0 (F -1 µ (p)) -dp. (4.4)
Dealing in a symmetric way with the positive part ensures the convergence of the expectations. The inequality (t -x) + ≤ t + + x -for t, x ∈ R combined with (4.4) yields some uniform integrability property ensuring the last convergence.

Proof of Proposition 4.6. Let us deal with the second modification, the proof being similar for the first or third modifications or with no modification. We also deal with the random choice of points, the previous lemma replacing the strong law of large numbers for the deterministic choice. The a.s. weak convergence of the empirical measure µ I to µ as I → ∞ is a consequence of the strong law of large numbers. Still by the strong law of large numbers,

1 I I i=1 X i -1 J J j=1 Y j converges a.s. to R xµ(dx) - R xν(dx) = 0 as I, J → ∞. For t ∈ R, since x → (t -x) + is Lipschitz continuous with constant 1, R (t -x) + ν J (dx) - R (t -x) + ν(dx) ≤ 1 I I i=1 X i - 1 J J j=1 Y j + 1 J J j=1 (t -Y j ) + - R (t -x) + ν(dx) .
Dealing with the second term in the right-hand side by the strong law of large numbers, one deduces that ϕ ν J (t) converges a.s. to ϕ ν (t) as I, J → ∞. The monotonicity in t and the continuity of the limit in the same variable ensures that a.s. the convergence holds for all t ∈ R. In the same way a.s. for all t ∈ R, ϕ µ I (t) converges to ϕ µ (t). Since ϕ µ I ∨ν J (t) = max (ϕ µ I (t), ϕ ν J (t)), we deduce that as I, J → ∞, a.s., for all t ∈ R,

ϕ µ I ∨ν J (t) converges to max (ϕ µ (t), ϕ ν (t)) = ϕ ν (t) since µ ≤ cx ν. For t ∈ R and h > 0, 1 h (ϕ µ I ∨ν J (t) -ϕ µ I ∨ν J (t -h)) = 1 h t t-h F µ I ∨ν J (x)dx ≤ F µ I ∨ν J (t) ≤ 1 h t+h t F µ I ∨ν J (x)dx = 1 h (ϕ µ I ∨ν J (t + h) -ϕ µ I ∨ν J (t))
where the left-hand and the right-hand sides of the inequality respectively converge to

1 h R (t -x) + -(t -h -x) + ν(dx) = 1 h t t-h F ν (x)dx and 1 h R (t + h -x) + -(t -x) + ν(dx) = 1 h t+h t F ν (x)dx.
Letting h → 0, we conclude that a.s., F µ I ∨ν J (t) converges to F ν (t) for all t ∈ R such that ν({t}) = 0 which is equivalent to the weak convergence of µ I ∨ ν J to ν.

We now prove the a.s. weak convergence of µ I ∧ ν J to µ. We have by construction ϕ µ I ∧ν J (t) ≤ min(ϕ µ I (t), ϕ ν J (t)). This gives a.s, that for all t ∈ R, lim sup I,J→+∞ ϕ µ I ∧ν J (t) ≤ ϕ µ (t) since µ ≤ cx ν. Now, we observe that the convergence of ϕ µ I (resp. ϕ ν J ) to ϕ µ (resp. ϕ ν ) is uniform on R. In fact, ϕ µ I and ϕ µ are nondecreasing functions such that lim t→-∞ ϕ µ I (t) = lim t→-∞ ϕ µ (t) = 0 and ϕ µ is continuous. From the a.s. simple convergence, we deduce that a.s., ϕ µ I converges uniformly to ϕ µ on (-∞, T ], for any T > 0. The same reasoning gives that a.s., φµ I converges uniformly to φµ on [T, +∞). Since, by (2.1) and (2.2),

t - R xµ(dx) = R (t -x) + µ(dx) - R (x -t) + µ(dx) = ϕ µ (t) -φµ (t), t -1 I I i=1 X i = ϕ µ I (t) -φµ I (t) and 1 I I i=1 X i → I→+∞ R xµ(dx)
, we get that a.s., ϕ µ I converges uniformly to ϕ µ on R. The same argument works for the uniform convergence of ϕ ν J to ϕ ν . Thus, for ε > 0 there exists N such that for I, J ≥ N ,

∀t ∈ R, max(|ϕ µ I -ϕ µ (t)|, |ϕ ν J -ϕ ν (t)|) ≤ ε. We deduce ∀t ∈ R, ϕ µ I (t) = ϕ µ I (t) ∧ ϕ ν J (t) + ϕ µ I (t) ∨ ϕ ν J (t) -ϕ ν J (t) ≤ ϕ µ I (t) ∧ ϕ ν J (t) + ϕ µ (t) ∨ ϕ ν (t) -ϕ ν (t) + 2ε = ϕ µ I (t) ∧ ϕ ν J (t) + 2ε, since µ ≤ cx ν. Since ϕ µ I -2ε is a convex function below ϕ µ I ∧ ϕ ν J , we get that ∀I, J ≥ N, t ∈ R, ϕ µ I (t) -2ε ≤ ϕ µ I ∧ν J (t).
Since ε > 0 is arbitrary, this gives a.s, that for all t ∈ R, lim inf I,J→+∞ ϕ µ I ∧ν J (t) ≥ ϕ µ (t). From the a.s. convergence for all t ∈ R of ϕ µ I ∧ν J (t) to ϕ µ (t), we deduce as before the convergence for all t ∈ R such that µ({t}) = 0 of F µ I ∧ν J (t) to F µ (t).

Remark 4.8. The proof of Proposition 4.6 does not use the independence between (X i ) i≥1 and (Y i ) i≥1 when they are sampled. In fact, among the three proposed modifications, only the third one exploits the independence in order to get the best estimation of the common expectation R xµ(dx) = R yµ(dy) from the samples (X i ) 1≤i≤I and (Y j ) 1≤j≤J . An interesting thing to notice is that the construction of approximate probability measures that preserve the convex order can be easily extended to the multi-marginal case. More precisely, let ≥ 2, I 1 , . . . , I be positive integers and µ 1 , . . . , µ be probability measures on R such that µ 1 ≤ cx . . . ≤ cx µ and R |x|µ (dx) < ∞. We either set

x k i = F -1 µ k 2i-1 I k for i ∈ {1, . . . , I k } and k ∈ {1, . . . , } or x k i = X k i
where the random vectors ((X 1 i , . . . , X i )) i≥1 are i.i.d. according to µ 1 ⊗ . . . ⊗ µ . Next we choose (x k i ) 1≤i≤I k ,1≤k≤ as one of the following vectors : -(

x k i ) 1≤i≤I k ,1≤k≤ , -( x k i -1 I k I k i=1 x k i + R xµ 1 (dx)) 1≤i≤I k ,1≤k≤ , -( x k i -1 I k I k i=1 x k i + 1 I 1 I 1 i=1 x 1 i ) 1≤i≤I k ,1≤k≤ , -x k i + j=1 I j / σ 2 j l=1 I l / σ 2 l ( 1
I j I j i=1 x j i -1 I k I k i=1 x k i ) 1≤i≤I k ,1≤k≤
, where σ 2 j = 1

I j -1 I j i=1 ( x j i ) 2 -1 I j (I j -1) ( I j i=1 x j i ) 2 and set µ k I k = 1 I k I k i=1 δ x k i for k ∈ {1, . . . , }.
For the first choice, the expectations of the measures µ k I k are in general distinct and, to obtain an increasing sequence of approximate probability measures for the convex order, we compute either

µ 1 I 1 , µ 1 I 1 ∨ µ 2 I 2 , (µ 1 I 1 ∨ µ 2 I 2 ) ∨ µ 3 I 3 , . . . , ((µ 1 I 1 ∨ µ 2 I 2 ) ∨ . . .) ∨ µ I or µ 1 I 1 ∧ (. . . ∧ (µ -1 I -1 ∧ µ I )), . . . , µ -2 I -2 ∧ (µ -1 I -1 ∧ µ I ), µ -1 I -1 ∧ µ I , µ I .
For the other choices, the expectations are equal and the order in which the suprema ∨ (resp. the infima ∧) are computed no longer matters. Indeed, for f, g, h : R → R, we have max(f, g, h) = max(max(f, g), h) and min(f, g, h) = min((min(f, g)), h)). One easily generalizes the proof of Proposition 4.6 to obtain the following result. Proposition 4.9. Let µ 1 , . . . , µ be probability measures on R such that µ

1 ≤ cx . . . ≤ cx µ and R |x|µ (dx) < ∞. Then, as I 1 , . . . , I → +∞, µ 1 I 1 , µ 1 I 1 , µ 1 I 1 ∨ µ 2 I 2 , (µ 1 I 1 ∨ µ 2 I 2 ) ∨ µ 3 I 3 , . . . , ((µ 1 I 1 ∨µ 2 I 2 )∨. . .)∨µ I (resp. µ 1 I 1 ∧(. . .∧(µ -1 I -1 ∧µ I )), . . . , µ -2 I -2 ∧(µ -1 I -1 ∧µ I ), µ -1 I -1 ∧ µ I , µ I ) converge a.
s. weakly to µ 1 , . . . , µ .

Algorithms and complexity

In this section, we present the algorithms to calculate µ I ∧ ν J and µ I ∨ ν J defined in (4.1) and (4.2) when µ I = I i=1 p i δ x i and ν J = J j=1 q j δ y j are probability measures with finite support. In this case, the probability measures µ I ∧ ν J and µ I ∨ ν J also have a finite support, and the number of points in their support is discussed in Subsection 5.3. When

x 1 < • • • < x I and y 1 < • • • < y J , the number of operations required to calculate µ I ∧ ν J or µ I ∨ ν J is in O(I + J).
Otherwise, the computation cost is in O(I log(I) + J log(J)) since one has first to sort the points. 5.1. Algorithm for the calculation of the supremum of two probability measures with finite support. Let µ I = I i=1 p i δ x i with (p 1 , . . . , p I ) ∈ (0, 1] I and ν J = J j=1 q j δ y j with (q 1 , . . . , q J ) ∈ (0, 1] J be two probability measures on the real line. Up to reducing the number of points and sorting them, we may suppose that x 1 < x 2 < . . . < x I (resp y 1 < y 2 < . . . < y J ). The equality µ I ∨ icx ν J = μI ∨ dcx νJ implies that it is enough to focus on the computation of µ I ∨ dcx ν J .

The functions ϕ µ I and ϕ ν J are piecewise affine with slopes taking respectively their values in {0, p 1 , p 1 +p 2 , . . . , I-1 i=1 p i , 1} and {0, q 1 , q 1 +q 2 , . . . , J-1 j=0 q j , 1}. As a consequence ϕ µ I ∨ ϕ ν J is piecewise affine with slopes taking values in the union of these two sets so that µ I ∨ν J is a probability measure supported on K ≤ I +J -1 points z 1 < z 2 < . . . < z K . Our aim is now to provide an efficient algorithm aimed at expliciting µ I ∨ dcx ν J = K k=1 k δ z k . Let w 1 < w 2 < . . . < w L with L ≥ I ∨ J be such that {w 1 , . . . , w L } = {x 1 , . . . , x I } ∪ {y 1 , . . . , y J }.

For convenience, we set w 0 = w 1 -1. Of course, ϕ µ I (w 0 ) = F µ I (w 0 ) = ϕ ν J (w 0 ) = F ν J (w 0 ) = 0. We now proceed inductively on ∈ {1, . . . , L} remarking that the functions ϕ µ I and ϕ ν J are affine with respective slopes F µ I (w -1 ) and F ν J (w -1 ) on [w -1 , w ] :

1. We first compute

ϕ µ I (w ) = ϕ µ I (w -1 ) + F µ I (w -1 )(w -w -1 ), F µ I (w ) = F µ I (w -1 ) + µ I ({w }) ϕ ν J (w ) = ϕ ν J (w -1 ) + F ν J (w -1 )(w -w -1 ), F ν J (w ) = F ν J (w -1 ) + ν J ({w }). 2. If (ϕ ν J (w ) -ϕ µ I (w ))(ϕ ν J (w -1 ) -ϕ µ I (w -1
)) < 0, then, computing the intersection point between ϕ µ I and ϕ ν J on [w -1 , w ] and the increase of the slope of ϕ µ I ∨ ϕ ν J at this point, we obtain that

µ I ∨ dcx ν J |ϕ ν J (w ) -ϕ µ I (w )|w -1 + |ϕ ν J (w -1 ) -ϕ µ I (w -1 )|w |ϕ ν J (w ) -ϕ µ I (w )| + |ϕ ν J (w -1 ) -ϕ µ I (w -1 )| = µ I ∨ dcx ν J ((w -1 , w )) = |ϕ µ I (w ) -ϕ µ I (w -1 ) -ϕ ν J (w ) + ϕ ν J (w -1 )| w -w -1 = |ϕ ν J (w ) -ϕ µ I (w )| + |ϕ ν J (w -1 ) -ϕ µ I (w -1 )| w -w -1 = |F ν (w -1 ) -F µ (w -1 )|.
Otherwise, we either have ∀t ∈ (w -1 , w )

ϕ µ I ∨ dcx ν J (t) = ϕ µ I (t) or ∀t ∈ (w -1 , w ) ϕ µ I ∨ dcx ν J (t) = ϕ ν J (t)
, and the slope of ϕ µ I ∨ ϕ ν J is constant on (w -1 , w ) so that µ I ∨ dcx ν J ((w -1 , w )) = 0.

3. The mass µ I ∨ dcx ν J ({w }) is obtained in the following way.

(i) If ϕ µ I (w ) = ϕ ν J (w ), then the slopes of ϕ µ I ∨ ϕ ν J on (w -1 , w ) and on (w , w +1 ) (convention w L+1 = w L + 1) are respectively equal to

F µ I (w -1 ) ∧ F ν J (w -1 ) and F µ I (w ) ∨ F ν J (w ) so that µ I ∨ dcx ν J ({w }) = F µ I (w ) ∨ F ν J (w ) -F µ I (w -1 ) ∧ F ν J (w -1 ). (ii) If ϕ µ I (w ) > ϕ ν J (w ) then ϕ µ I ∨ ϕ ν J is equal to ϕ µ I on a neighbourhood of w so that µ I ∨ dcx ν J ({w }) = µ I ({w }). (iii) If ϕ µ I (w ) < ϕ ν J (w ), then, in a symmetric way, µ I ∨ dcx ν J ({w }) = ν J ({w }).
Of course, the slope of ϕ µ I ∨ ϕ ν J is equal to 0 on (-∞, w 1 ) and to 1 on (w L , +∞) so that

µ I ∨ dcx ν J ((-∞, w 1 ) ∪ (w L , +∞)) = 0. Moreover, µ I ∨ dcx ν J ({w 1 }) = µ I ({w 1 }) ∨ ν J ({w 1 }).
On the other hand, by (2.1), ϕ µ I (w L ) = w L -I i=1 p i x i , ϕ ν J (w L ) = w L -J j=1 q j y j and when either I i=1 p i x i < J j=1 q j y j and µ I ({w L }) = 0 or I i=1 p i x i > J j=1 q j y j and ν J ({w L }) = 0, then µ I ∨ dcx ν J ({w L }) = 0. Otherwise, and in particular when I i=1 p i x i = J j=1 q j y j , µ I ∨ dcx ν J ({w L }) > 0. Algorithm 1 below implements these ideas to compute µ I ∨ dcx ν J . 5.2. Algorithm for the calculation of the infimum of two probability measures with finite support. We consider the same setting as in Subsection 5.1. The equality µ I ∧ icx ν J = μI ∧ dcx νJ implies that it is enough to focus on the computation of µ I ∧ dcx ν J . Recall that w 1 < w 2 < . . . < w L with L ≥ I ∨J are such that {w 1 , . . . , w L } = {x 1 , . . . , x I }∪ {y 1 , . . . , y J }. Since convex hulls will appear in the characterization, when d = 1, of the Wasserstein projections considered in Section 3, rather than using the characterization of µ I ∧ dcx ν J deduced from Lemma 4.2 which would lead to some algorithm similar to the previous one, we compute the convex hull ϕ µ I ∧ dcx ν J of ϕ µ I ∧ ϕ ν J . To do so, we observe that its epigraph is the convex hull of the union of the epigraphs of ϕ µ I and ϕ ν J . Thus, we can use Andrew's monotone chain convex hull algorithm. We describe here briefly how it works in our case. We exclude the trivial case L = 1 and assume that L ≥ 2. The principle is to calculate for ∈ {2, . . . L} the convex hull of {(y, w)

: w ∈ [w 1 , w ], y ≥ ϕ µ I ∧ ϕ ν J (w)}. It is fully described by indices i (1) < • • • < i (n ) such that this convex hull is equal to ∪ n -1 n=1 {(y, w) : w ∈ [w i (n) , w i (n+1) ], y ≥ w i (n+1) -w w i (n+1) -w i (n) ϕ µ I ∧ϕ ν J (w i (n) )+ w-w i (n) w i (n+1) -w i (n) ϕ µ I ∧ ϕ ν J (w i (n+1) )}. Obviously n 2 = 2, i 2 (1) = 1, i 2 (2) =
2 and from this initialization, we proceed by induction. To calculate the indices i +1 from the indices i , Andrew's monotone chain convex hull algorithm works as follows: {i +1 (1), . . . , i (n

+1 )} is equal to {1}∪{i (n) : 2 ≤ n ≤ n , ϕµ I ∧ϕν J (w +1 )-ϕµ I ∧ϕν J (w i (n) ) w +1 -w i (n) > ϕµ I ∧ϕν J (w i (n) )-ϕµ I ∧ϕν J (w i (n-1) ) w i (n) -w i (n-1)
} ∪ { + 1}. This means that we remove all the previous points that would give a nonincreasing slope. Once we have calculated the indices, the weights of the measure µ I ∧ dcx ν J are obtained as the difference of two successive slopes. This is presented in Algorithm 2.

Examples.

In this subsection, we discuss on different examples the number K of points in the support of µ I ∨ ν J (resp. µ I ∧ ν J ), in the case I i=1 p i x i = J j=1 q j y j . We show that the maximal value of K is I + J -1 (resp. I + J -2).

We start by studying examples for µ I ∨ ν J and we suppose for convenience that I ≤ J.

1. If I = 1, then necessarily µ I ≤ cx ν J so that µ I ∨ ν J = ν J and K = J, 2. If I ≥ 2, then the support of µ I ∨ν J contains at least the two distinct points x 1 ∧y 1 and x I ∨y J . Now either µ I ∨ν J ∈ {µ I , ν J } or neither µ I ≥ cx ν J nor ν J ≥ cx µ I which implies that K ≥ 3 : the sign of ϕ ν J -ϕ µ I changes on the interval (x 1 ∧ y 1 , x I ∨ y J ), which implies that the cumulative distribution function F equal to the slope of ϕ µ I ∨ ϕ ν J is not constant on this interval. In particular, when I ≥ 3 then K ≥ 3. Notice that K may be equal to 3 even if I is arbitrarily large. For instance, for the centered probability measures

µ I = 1 3 δ -1 + 2 3(I -1) I i=2 δ i-2 I-2
and ν I = 2 3(I -1)

I-1 i=1 δ 1-I+i I-2 + 1 3 δ 1 with I ≥ 3, µ I ∨ ν I = 1 3 (δ -1 + δ 0 + δ 1 ). Indeed ϕ µ I (t) = 1+t 3 for t ∈ [-1, 0] whereas ϕ ν I (0) = 2 3(I-1) I-1 i=1 i-1 I-2 = 1 3 so that ϕ ν I (t) = 1+2t 3 for t ∈ [0, 1].
With the convexity of ϕ µ I and ϕ ν I and the equality of these functions at the extreme points -1 and 1 (deduced at point 1 from the equality of the expectations), one concludes that ϕ µ I ∨ ϕ ν I (t) = (1+t) + +t + +(t-1) + 3 .

3. As explained in Subsection 5.1, we necessarily have K ≤ I + J -1. For α ∈ ( 2 I , 2 I-1 ),

µ I = 1 I I i=1 δ x i with x i = (1 + α)i I + 1 - α 2 and ν I+1 = 1 I + 1 I i=0 δ i I ,
the cardinality of the support of µ I ∨ ν I+1 is 2I. Indeed, it is easy to check that the condition α < 2 I-1 is necessary and sufficient to ensure i-1 I < x i < i I for all i ∈ {1, . . . , I}. Under this condition, it is enough to check that for all i ∈ {1, . . . , I}, ϕ ν I+1 (x i ) > ϕ µ I (x i ) and for all i ∈ {1, . . . I -1},

ϕ µ I ( i I ) > ϕ ν I+1 ( i I ) to ensure that F indeed takes the 2I + 1 possible different values. Notice that on [ i-1 I , i I ], ϕ ν I+1 (t) = i I+1 t -i-1 2I and that on [x i , x i+1 ], ϕ µ I (t) = i I t -(1+α)(i+1) 2(I+1) + α 2 . The inequality ϕ ν I+1 (x i ) > ϕ µ I (x i ) multiplied by 2I(I+1) 2 i therefore writes (2-α(I -1))(I +1-i) > 0.
On the other hand, the inequality

ϕ µ I ( i I ) > ϕ ν I+1 ( i I ) multiplied by 2I 2 (I+1) i rewrites (αI -2)(I -i) > 0.
4. When the probability measures µ I and ν J are uniform :

p i = 1
I for all i ∈ {1, . . . , I} and q j = 1 J for all j ∈ {1, . . . , J}. Then K ≤ I + J -gcd(I, J) (in particular K ≤ I when J = I) since the cardinality of { i I : i ∈ {0, . . . , I}} ∪ { j J : j ∈ {0, . . . , J}} is I + J + 1 -gcd(I, J). Let I = I/ gcd(I, J) and J = J/ gcd(I, J). These integers are prime together and we have i I = j J ⇐⇒ iJ = jI ⇐⇒ i = κI and j = κJ , with κ ∈ {0, . . . , gcd(I, J)} when 0 ≤ i ≤ I and 0 ≤ j ≤ J. Therefore, the cardinal of { i I : i ∈ {0, . . . , I}} ∩ { j J : j ∈ {0, . . . , J}} is 1 + gcd(I, J), which gives the claim. We now discuss the number K of points in the support of µ I ∧ ν J and suppose again I ≤ J. We assume I ≥ 2, otherwise µ I ∧ ν J has clearly one element (K = 1) and is the Dirac mass at x 1 .

1 . If I = J, x i = y i and p i = q i for 1 ≤ i ≤ I, we have µ I ∧ ν J = µ I and thus K = I.

2 . If I = 2, x 1 ≤ min j=1,...,J y j and x 2 ≥ max j=1,...,J y j , then for any j ∈ {1, . . . , J}, we have

y j = x 2 -y j x 2 -x 1 x 1 + y j -x 1 x 2 -x 1 x 2 . Thus, if p 1 = J j=1 q j
x 2 -y j

x 2 -x 1 and p 2 = J j=1 q j y j -x 1

x 2 -x 1 , we have ν J ≤ cx µ I and thus K = J.

. We have

K ≤ I + J -2. It is clear that K ≤ I + J since the support of µ I ∧ ν J is included in {x 1 , . . . , x I } ∪ {y 1 , . . . , y J }. If x 1 < y 1 , then x 1 is not in the support of µ I ∧ ν J since ϕ µ I ∧ν J (y 1 ) = 0.
Thus, at most one element of {x 1 , y 1 } is in the support of µ I ∧ ν J . If x I > y J , then x I is not in the support of µ I ∧ ν J , otherwise we would have ϕ µ I ∧ν J (x I -) < 1 and thus

ϕ µ I ∧ν J (y J ) > ϕ µ I ∧ν J (x I ) + y J -x I = x I - I i=1 p i x i + y J -x I = y J - J j=1
q j y j = ϕ ν J (y J ).

Thus, at most one element of {x I , y J } is in the support of µ I ∧ ν J . We eventually get K ≤ I + J -2. This bound is attained by the previous example: for α ∈ ( 2 I , 2 I-1 ),

µ I = 1 I I i=1 δ x i with x i = (1+α)i I+1 -α 2 and ν I+1 = 1 I+1 I i=0 δ i I
, the cardinality of the support of µ I ∧ ν I+1 is 2I -1. Let us recall that we have already shown that for all i ∈ {1, . . . I}, x i ∈ ( i-1 I , i I ) and ϕ ν I+1 (x i ) > ϕ µ I (x i ) and, for all i ∈ {1, . . . I -1},

ϕ µ I ( i I ) > ϕ ν I+1 ( i I ). Thus, to prove that the support of µ I ∧ ν J is {x i , 1 ≤ i ≤ I} ∪ { i I , 1 ≤ i ≤ I -1}
, it is necessary and sufficient to prove that

ϕ µ I (x i ) -ϕ ν I+1 i-1 I x i -i-1 I < ϕ ν I+1 i I -ϕ µ I (x i ) i I -x i , i ∈ {1, . . . , I}, ϕ ν I+1 i I -ϕ µ I (x i ) i I -x i < ϕ µ I (x i+1 ) -ϕ ν I+1 i I x i+1 -i I , i ∈ {1, . . . , I -1}.
We observe that

ϕ µ I (x i ) = 1+α 2I(I+1) i(i -1) and ϕ ν I+1 i I = 1 2I(I+1) i(i + 1
). After simple but tedious calculations, we can check that the first condition is equivalent to 0 < (2 -α (I -1)) (I + 1 -i), which is true since α < 2 I-1 . Similarly, the second condition is equivalent to (I -i) (2 -αI) < 0, which is true since α > 2 I . 4 . If I ≥ 2, the support of µ I ∧ ν J has at least two elements, i.e. K ≥ 2. Since m := I i=1 p i x i = J j=1 q j y j and the elements x i (resp. y j ) are distinct, we have ϕ µ I (m) > 0 (resp. ϕ ν J (m) > 0). If K = 1, we would have then necessary µ I ∧ ν J = δ m , which would be in contradiction with the previous statement. Notice that K may be equal to 2 even if I is arbitrarily large. For I ≥ 2, we take when I -1 = 2k is even

µ I = 1 I -1 k i=1 1 2 δ -1-i k + δ -1+ i k + 1 2 δ 1 , ν I = 1 2 δ -1 + 1 I -1 k i=1 1 2 δ 1-i k + δ 1+ i k ,
and when I -1 = 2k + 1 is odd,

µ I = 1 I k i=0 1 2 δ -1-i k + δ -1+ i k + 1 2 δ 1 , ν I = 1 2 δ -1 + 1 I k i=0 1 2 δ 1-i k + δ 1+ i k .
By construction, the support of µ I (resp. ν I ) has exactly I points. Moreover, one easily checks that µ I ∧ ν I = 1 2 (δ -1 + δ 1 ).

5.4. Irreducible components. Specializing Theorem A.4 of [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF] and its proof to the case of probability measures, one has the following decomposition.

Proposition 5.1. Let µ ≤ cx ν be two distinct probability measures such that µ, ν ∈ P 1 (R).

There exists a countable family of disjoint intervals

(t n , t n ), 1 ≤ n ≤ N , with N ∈ N * ∪{∞} such that {t ∈ R, ϕ µ (t) < ϕ ν (t)} = ∪ N n=1 (t n , t n ). We have µ((t n , t n )) > 0 for any 1 ≤ n ≤ N . Let E = ∩ 1≤n≤N {(-∞, t n ] ∪ [t n , +∞)}, µ 0 (dx) = 1 E (x)µ(dx), µ n (dx) = 1 (t n ,tn) (x)µ(dx)
so that µ = N n=0 µ n . Then, there exists a unique decomposition ν = N n=0 ν n such that

ν n (R) = µ n (R) and µ n µ n (R) ≤ cx ν n µ n (R) , for 1 ≤ n ≤ N and, when µ 0 (R) > 0, µ 0 µ 0 (R) ≤ cx ν 0 µ 0 (R) .
It is given by ν 0 = µ 0 and, for 1 ≤ n ≤ N , by

ν n (dx) = 1 (t n ,tn) (x)ν(dx) + (F ν (t n ) -F µ (t n ))δ t n (dx) + (F µ (t n -) -F ν (t n -))δ tn (dx).
Last, if P is a martingale coupling between µ and ν, P (dx, dy) = µ 0 (dx)δ x (dy) + N n=1 µ((t n , t n ))P n (dx, dy), where P n is a martingale coupling between µ n µ n (R) and ν n µ n (R) .

The intervals ((t n , t n )) 1≤n≤N are called the irreducible components of the couple (µ, ν). According to Corollary A.7 of [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF], solving a martingale optimal transport problem between µ and ν is equivalent to solving the martingale optimal transport problem with the same payoff function between µ n µ n (R) and ν n µ n (R) for each 1 ≤ n ≤ N . This motivates our interest in studying the irreducible components of (µ I ∧ ν J , ν J ) and (µ I , µ I ∨ ν J ): to solve numerically the corresponding discrete MOT problem, this replace a linear programming problem by N linear programming problems of smaller dimension, where N is the corresponding number of irreducible components.

We now consider two probability measures with finite supports µ I = I i=1 p i δ x i and ν J = J j=1 q j δ y j with

x 1 < • • • < x I , y 1 < • • • < y J , p i > 0 for i ∈ {1, .
. . , I} and q j > 0 for j ∈ {1, . . . , J} and I i=1 p i x i ≤ J j=1 q j y j . Since, by Lemma 4.1,

R xµ I ∨ dcx ν J (dx) = R xµ I (dx) and R xµ I ∧ dcx ν J (dx) = R xν J (dx), µ I ≤ cx µ I ∨ dcx ν J and µ I ∧ dcx ν J ≤ cx ν J .
We are interested in describing precisely and computing the irreducible components of (µ I , µ I ∨ dcx ν J ) (resp. (µ I ∧ dcx ν J , ν J )). Of course, there are finitely many components. The computation of the irreducible components of (µ

I , µ I ∨ icx ν J ) (resp. (µ I ∧ icx ν J , ν J )) when I i=1 p i x i > J j=1 q j y j is easily deduced from the equality µ I ∨ icx ν J = μI ∨ dcx νJ (resp. µ I ∧ icx ν J = μI ∧ dcx νJ ).
Irreducible components of (µ I , µ I ∨ dcx ν J ) when I i=1 p i x i ≤ J j=1 q j y j . The irreducible components are defined as the largest intervals on which ϕ µ I < ϕ µ I ∨ dcx ν J . We have

µ I ∨ dcx ν J = K k=1 k δ z k with z 1 < • • • < z K , k > 0,
and Algorithm 1 calculates the points z k and the weights k . Among these points, we focus on the points such that

ϕ µ I (z k ) = ϕ µ I ∨ dcx ν J (z k ). Thus, we define v 1 < v 2 < . . . < v M such that {v 1 , . . . , v M } = {z k : 1 ≤ k ≤ K and ϕ µ I (z k ) = ϕ µ I ∨ dcx ν J (z k )}.
(5.1)

For t ≥ z K , ϕ µ I (t) ≤ ϕ µ I ∨ dcx ν J (t) = t -I i=1 p i x i . Since the function ϕ µ I is convex and equal to t -I i=1 p i x i for t ≥ x I , one deduces that ϕ µ I (t) = t -I i=1 p i x i for t ≥ z K . Since ϕ µ I (t) = ϕ µ I ∨ dcx ν J (t) = 0 for t ≤ z 1 and we have M ≥ 1 and all the irreducible components lie in (z 1 , z K ). If K = 1, µ I = µ I ∨ dcx ν J = δ x 1
and there is no irreducible component. Otherwise, we have M ≥ 2, which we assume now.

Remark 5.2. When M ≥ 3, looking at the irreducible components reduces the dimension of the discrete MOT problem. We notice that if I ≥ 2 and I i=1 p i x i = J j=1 q j y j , then

M = 2 =⇒ µ I ≤ cx ν J . In fact, M = 2 implies that ϕ µ I (z k ) < ϕ ν J (z k ) for 2 ≤ k ≤ K -1. Since 0 = ϕ µ I (z 1 ) ≥ ϕ ν J (z 1 ) and z K -I i=1 p i x i = ϕ µ I (z K ) ≥ ϕ ν J (z K ) ≥ z K -I i=1 p i x i + by Jensen's inequality, we get ϕ µ I (z 1 ) = ϕ ν J (z 1 ) and ϕ µ I (z K ) = ϕ ν J (z K ) and thus µ I ≤ cx ν J . We have ϕ µ I ≤ ϕ µ I ∨ ϕ ν J = ϕ µ I ∨ dcx ν J and, for any m ∈ {1, . . . , M }, ϕ µ I (v m ) = ϕ µ I ∨ dcx ν J (v m ).
We can therefore compare the derivatives and get

F µ I ∨ dcx ν J (v m -) ≤ F µ I (v m -), F µ I ∨ dcx ν J (v m ) ≥ F µ I (v m ). Let m ≤ M -1 and k m be such that v m = z km . If F µ I ∨ dcx ν J (v m ) = F µ I (v m ) then, since ϕ µ I ∨ dcx ν J ≥ ϕ µ I and the derivative F µ I ∨ dcx ν J of ϕ µ I ∨ ϕ ν J is constant and equal to F µ I ∨ dcx ν J (v m ) on (z km , z km+1 ) whereas the one F µ I of ϕ µ I is nondecreasing, one has ∀t ∈ [z km , z km+1 ], ϕ µ I ∨ dcx ν J (t) = ϕ µ I (t) and, in particular, k m+1 = k m + 1 and µ I ((v m , v m+1 )) = 0. In a symmetric way, if F µ I (v m+1 -) = F µ I ∨ dcx ν J (v m+1 -), then k m = k m+1 -1, µ I ((v m , v m+1 )) = 0 and ϕ µ I (t) = ϕ µ I ∨ dcx ν J (t) for t ∈ [v m , v m+1 ]. Thus, µ I ((v m , v m+1 )) > 0 implies that F µ I ∨ dcx ν J (v m ) > F µ I (v m ) and F µ I (v m+1 -) > F µ I ∨ dcx ν J (v m+1 -), which gives that ∀t ∈ (v m , v m+1 ), ϕ µ I ∨ dcx ν J (t) > ϕ µ I (t).
Proposition 5.3. We assume

I i=1 p i x i ≤ J j=1 q j y j so that µ I ≤ cx µ I ∨ dcx ν J . Let v 1 < v 2 < . . . < v M be defined by (5.1). We have {t ∈ R, ϕ µ I (t) < ϕ µ I ∨ dcx ν J (t)} = 1≤m≤M -1, µ I ((vm,v m+1 ))>0 (v m , v m+1 ),
i.e. the intervals (v m , v m+1 ) such that µ I ((v m , v m+1 )) > 0 are the irreducible components of (µ I , µ I ∨ dcx ν J ). Moreover µ 0 I = M m=1 µ I ({v m })δ vm and any martingale coupling P between µ I and µ I ∨ dcx ν J is such that for all m ∈ {1, . . . , M },

P ({(v m , v m )}) = µ I ({v m }) and for all m ∈ {1, . . . , M -1}, P ((v m , v m+1 ) × [v m , v m+1 ]) = µ I ((v m , v m+1 )).
Any martingale optimal transport problem between the marginals µ I and µ I ∨ dcx ν J may be decomposed into the martingale optimal transport problems with the same payoff function but between the marginals

1 µ I ((vm,v m+1 )) I i=1 p i 1 {vm<x i <v m+1 } δ x i and 1 µ I ((v m , v m+1 )) (F µ I ∨ dcx ν J (v m ) -F µ I (v m ))δ vm + K k=1 k 1 {vm<z k <v m+1 } δ z k + (F µ I (v m+1 -) -F µ I ∨ dcx ν J (v m+1 -))δ v m+1
for m such that µ I ((v m , v m+1 )) > 0. Doing so, we replace one linear programming problem by N linear programming problems of smaller size, where N = M -1 m=1 1 {µ I ((vm,v m+1 ))>0} is the number of irreducible components.

Setting

η m-= F µ I (v m -) -F µ I ∨ dcx ν J (v m -), η m = µ I ({v m }) and η m+ = F µ I ∨ dcx ν J (v m ) - F µ I (v m )
, Algorithm 1 can be modified in the following way to compute (v m , η m-, η m , η m+ ) 1≤m≤M . First initialize m ← 1 before the while loop. In this while loop,

1. if ( ϕ ν -ϕ µ )(ϕ ν -ϕ µ ) < 0, add v m ← z k , η m-← 1 {ϕν >ϕµ} k , η m ← 0, η m+ ← 1 {ϕν <ϕµ} k , m ← m + 1 just before k ← k + 1, 2. if ϕ ν -ϕ µ = 0, add v m ← w, η m-← (F µ -F ν ) + , η m ← F µ -F µ , η m+ ← ( F ν -F µ ) + , m ← m + 1 just before k ← k + 1. 3. if ( ϕ µ -ϕ ν )p i 1 { z=x i } > 0, add v m ← w, η m-← 0, η m ← p i , η m+ ← 0, m ← m + 1 just before k ← k + 1.
Irreducible components of (µ I ∧ dcx ν J , ν J ) when I i=1 p i x i ≤ J j=1 q j y j . We start with a simple observation.

Lemma 5.4. For j ∈ {1, . . . , J -1}, we have

∃t 0 ∈ (y j , y j+1 ), ϕ ν J (t 0 ) = ϕ µ I ∧ dcx ν J (t 0 ) ⇐⇒ ∀t ∈ [y j , y j+1 ], ϕ ν J (t) = ϕ µ I ∧ dcx ν J (t).
Proof. The function ϕ ν J is affine on [y j , y j+1 ]. Since the function ϕ µ I ∧ dcx ν J is below ϕ ν J and convex, the equality of both functions at t 0 implies that they coincide on [y j , y j+1 ]. Now, we assume that I i=1 p i x i ≤ J j=1 q j y j so that µ I ∧ dcx ν J ≤ cx ν J and we consider

v 1 < • • • < v M such that {v 1 , . . . , v M } = {y j : 1 ≤ j ≤ J and ϕ ν J (y j ) = ϕ µ I ∧ dcx ν J (y j )}, (5.2)
and we denote by j m the index such that v m = y jm . We have ϕ ν J (y 1 ) = ϕ µ I ∧ dcx ν J (y 1 ) = 0 and thus M ≥ 1. We have ϕ ν J (y J ) = y J -J j=1 q j y j and ϕ µ I ∧ dcx ν J (y J ) ≥ y J -J j=1 q j y j since µ I ∧ dcx ν J and µ I have the same means. This gives ϕ ν J (y J ) = ϕ µ I ∧ dcx ν J (y J ) and v M = y J . We consider M ≥ 2 and 1 ≤ m ≤ M -1. By Lemma 5.4, we either have

j m+1 = j m + 1 or j m+1 > j m + 1 and ϕ ν J (t) > ϕ µ I ∧ dcx ν J (t) for t ∈ (v m , v m+1 ). Since j m+1 = j m + 1 ⇐⇒ ν J (v m , v m+1 ) = 0,
we get the following result.

Proposition 5.5. We assume

I i=1 p i x i ≤ J j=1 q j y j so that µ I ∧ dcx ν J ≤ cx ν J . Let v 1 < v 2 < . . . < v M be defined by (5.2). We have {t ∈ R, ϕ µ I ∧ dcx ν J (t) < ϕ µ I (t)} = 1≤m≤M -1, ν J ((vm,v m+1 ))>0 (v m , v m+1 ). i.e. the intervals (v m , v m+1 ) such that ν J ((v m , v m+1 )) > 0 are the irreducible components of (µ I ∧ dcx ν J , ν J ). Moreover ν 0 J = M m=1 µ I ∧ dcx ν J ({v m })δ vm and any martingale cou- pling P between µ I ∧ dcx ν J and ν J is such that for all m ∈ {1, . . . , M }, P ({(v m , v m )}) = µ I ∧ dcx ν J ({v m }) and for all m ∈ {1, . . . , M -1}, P ((v m , v m+1 ) × [v m , v m+1 ]) = µ I ∧ dcx ν J ((v m , v m+1 )).
Any MOT problem between the marginals µ I ∧ dcx ν J = K k=1 k δ z k and ν J = J j=1 q j δ y j may be decomposed into the MOT problems with the same payoff function but between the marginals

1 K k=1 k 1 {vm<z k <v m+1 } K k=1 k 1 {vm<z k <v m+1 } δ z k and 1 K k=1 k 1 {vm<z k <v m+1 } (F ν J (v m ) -F µ I ∧ dcx ν J (v m ))δ vm + J j=1 q j 1 {vm<y j <v m+1 } δ y j + (F µ I ∧ dcx ν J (v m+1 -) -F ν J (v m+1 -))δ v m+1
for m such that ν J ((v m , v m+1 )) > 0. Doing so, we replace one linear programming problem by N linear programming problems of smaller size, where N = M -1 m=1 1 {ν J ((vm,v m+1 ))>0} is the number of irreducible components.

To compute the irreducible components, it is convenient to add two coordinates η k and F k ν to each element in the list L in Algorithm 2. Instead of adding ( w, ϕ µ ∧ ϕ ν , 0, 0) to the list, one should add ( w, ϕ µ ∧ ϕ ν , 0, 0, 1 { w=y j ,ϕν ≤ϕµ} q j , F ν ). At the end of the algorithm, one computes {v

1 , . . . , v m } by m = 1, for k = 1 to s(L) if η k > 0, then v m = z k , µ I ∧ dcx ν J ({v m }) = k , F µ I ∧ dcx ν J (v m -)-F ν J (v m -) = F k -k -F k ν +η k , F ν J (v m )-F µ I ∧ dcx ν J (v m ) = F k ν -F k , m = m + 1.

Numerical experiments

In this section, we show how the different constructions of probability measures presented in the paper can be used to approximate the solution of a martingale optimal transport problem. In each case, we have used as approximating measures the empirical measures obtained of i.i.d. samples. This corresponds to the most generic practical case for the financial application: martingale asset prices are simulated to compute option prices by the Monte-Carlo method, and we want to reuse these simulations to compute bounds on these option prices. To solve numerically (1.1),we have used the linear programming solver GLPK † . With this solver, we are limited to use samples with I = J ≈ 100: the dimension of the matrix r ij is then about 10000. Other solvers may be more efficient, but since the dimension of the linear problem is the square of the number of points, the size of the samples which can be treated by exact solvers is limited. Approximated solvers such as the one proposed by [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] with entropic regularization seem to be a promising way to overcome this issue. The discussion on the maximal size of the samples that is possible to handle is of course interesting, but outside the scope of the paper. Here, we rather want to compare the different constructions and illustrate their relevance.

Last, we only consider in our numerical experiments approximations constructed from i.i.d. samples of the measures µ and ν. When the quantile functions F -1 µ and F -1 ν are known explicitly, it is possible to consider approximations with more evenly distributed points such as (1.2) or ηI = 1

I I-1 i=0 δ F -1
η ((i+1/2)/I) , η ∈ {µ, ν}. The former approximation, which requires in addition to know explicitly the antiderivative of the quantile functions, preserves the convex order. The latter is likely to do the same, but one has to check this property by using for instance Corollary 2.2 before using the linear programming solver. This check has the same computational complexity as calculating the probability measure μI ∧ νI equal to μI when μI ≤ cx νI , which ensures in any case the convex ordering. Besides, for the financial application, the quantile functions coming from the marginal distributions of a local volatility model or a stochastic volatility model are not known explictly in general. Moreover, when considering large i.i.d. samples with I = J it is possible to use Baker's construction (1.2) on the re-ordered modifications of empirical measures. This approach is made precise in Example6.3 and leads to very evenly spaced points in Figure 7. We have noticed in our numerical experiments almost no difference for the same final sample size between this procedure and the direct use of (1.2) when the (integrated) quantile functions can be calculated explicitly. This is why we only present here approximations obtained from i.i.d. samples. 6.1. An example with an explicit MOT. We start with an example where the martingale optimal transport is known explicitly. Let > 2,

µ(dx) = 1 2 1 [-1,1] (x)dx and ν(dx) = 1 4 1 [-2,2] (x)dx.
We consider the following martingale optimal transport problem:

min π∈Π M (µ,ν) R×R
|y -x| π(dx, dy).

For any π ∈ Π M (µ, ν), we have [START_REF] Hobson | Robust bounds for forward start options[END_REF] in Example 6.1. Note that the optimal coupling is neither the left-curtain coupling π lc (dx, dy)

R×R |y -x| 2 π(dx, dy) = R y 2 ν(dy) - R x 2 µ(dx) = 1. For > 2,
= 1 2 1 [-1,1] (x)dx 1 4 δ -x 2 -3 2 (dy) + 3 4 δ 3 2 x+ 1 2 (dy) nor the right-curtain coupling π rc (dx, dy) = 1 2 1 [-1,1] (x)dx 3 4 δ 3 2 x-1 2 (dy) + 1 4 δ 3 2 -x 2
(dy) introduced by Beiglböck and Juillet (2016), Theorem 1.5, see also [START_REF] Henry-Labordère | An explicit martingale version of the onedimensional Brenier theorem[END_REF] for an explicit calculation of those couplings. For the implementation, we consider independent random samples X 1 , . . . , X I distributed according to µ and Y 1 , . . . , Y I distributed according to ν, with I = 100. We set 2. This modification requires to know explicitly the common mean of µ and ν. Note that this is usually the case for the financial application: discounted assets prices are martingales, and their mean is given by the present values. We consider = 2.3. We have run the linear programming solver to solve (1.1) with (µ I ∧ ν I , ν I ) and (µ I , µ I ∨ ν I ) in Figure 1 and with ( µ I ∧ ν I , ν I ) and ( µ I , µ I ∨ ν I ) in Figure 2. In those figures, we have drawn the points with positive probability under the optimal transport, i.e. with the notation of (1.1), the points (x i , y j ) such that r ij > 0, where r is the optimal solution of (1.1). We expect to recover the line segments y = x+1 and y = x-1 for x ∈ [-1, 1]. From those figures, we observe that the martingale optimal transport map is much better approximated when we equalize the means. This can be heuristically explained as follows. Since µ I and ν I have the same mean, we have µ I ∧ ν I ≤ cx ν I and µ I ∧ ν I ≤ cx µ I . The latter inequality ensures from Corollary 2.2 that the lowest (resp. highest) point weighted by µ I ∧ ν I is greater (resp. smaller) than the lowest (resp. highest) point weighted by µ I . Thus, µ I ∧ ν I mainly weights points in [-1, 1]. Instead, when we do not equalize the means, we only have µ I ∧ ν I ≤ cx ν I , but there is absolutely no reason to have µ I ∧ ν I ≤ cx µ I . Therefore, µ I ∧ ν I may weight points that are significantly outside [-1, 1]: we observe in Figure 1 (left) many points close to 2. A similar phenomenon happens for the supremum. When we equalize the means, we have both µ I ≤ cx µ I ∨ ν I and ν I ≤ cx µ I ∨ ν I , and the latter condition gives that µ I ∨ ν I mainly weights points across [-2, 2]. In contrast, we do not have ν I ≤ cx µ I ∨ ν I and the points weighted by µ I ∨ ν I may not span [-2, 2]. In Figure 1 (right), the highest point weighted by µ I ∨ ν I is 1.215876, which is far from 2.

µ I = 1 I I i=1 δ X i , ν I = 1 I I i=1 δ Y i , µ I = 1 I I i=1 δ X i -XI , ν I = 1 I I i=1 δ Y i -ȲI ,
To complement this discussion, we now compare the value of the discrete optimal cost to the continuous one, i.e. to R×R |y -x| π (dx, dy) = 1. On the run of Figures 1 and2, we have obtained respectively 0.9294, 0.9101, 1.0773 and 1.0773 for (µ I ∧ ν I , ν I ), (µ I , µ I ∨ ν I ), ( µ I ∧ ν I , ν I ) and ( µ I , µ I ∨ ν I ). On 100 independent samples, we have obtained respectively the mean costs 0.7506, 0.7319, 1.0020, 1.0020 with the standard deviations 0.2148, 0.2148, 0.1400 and 0.1400. These results confirm that the cost is much better approximated in this example when we equalize the mean. Let us note here that the costs (and the optimal probability measure) obtained with ( µ I ∧ ν I , ν I ) and ( µ I , µ I ∨ ν I ) are the same because we already have on our samples µ I ≤ cx ν I , which is due to the fact that ν largely dominates µ for the convex order. From now on, we will only present numerical results for which we equalize the means. 6.2. The at-the-money discrete Asian call option. We consider the payoff function c(x, y) = (x + y) + = 2((x + y)/2) + , which corresponds to the case of a discret Asian call option where the average is calculated on two dates. We consider probability distributions µ ≤ cx ν with zero mean, so that the call is at-the-money.

Before to consider shifted lognormal laws, we first study the case

µ(dx) = 1 2 1 [-1,1] dx, ν(dx) = 1 12 1 [-3,-1] (x) + 1 3 1 [-1,1] (x) + 1 12 1 [1,3] (x) dx,
because we know then an explicit optimal martingale transport. If X is a uniform random variable on [-1, 1] and ξ an independent random variable such that P(ξ = -2) = P(ξ = 2) = 1/6 and P(ξ = 0) = 2/3, X + ξ ∼ ν. Since E[ξ] = 0, this shows that µ ≤ cx ν. We are interested in the optimal coupling π ∈ Π M (µ, ν) that minimizes R×R (x + y) + π(dx, dy). By Jensen's inequality, we have satisfies this condition and is therefore optimal. We observe that π (dx, dy) has the rightmonotone property for x ∈ [-1, 0] (see Definition 1.4 of [START_REF] Beiglböck | On a problem of optimal transport under marginal martingale constraints[END_REF]) and the left-monotone property for x ∈ [0, 1]. Let us note that we have not shown that π is the unique optimal martingale coupling. We also mention that, by similar arguments, π is also an optimal martingale coupling that minimizes R×R (3x -y) + π(dx, dy). We have plotted in Figure 3 the points weighted by the discrete MOT problem with the measures µ

R×R (x + y) + µ(dx)π Y |X (x, dy) ≥ R x + R yπ Y |X (x, dy) + µ(dx) = R (2x) + µ(dx) = 1 2 .

The equality in

I = 1 I I i=1 δ X i -XI and ν I = 1 I I i=1 δ Y i -ȲI ,
where X 1 , . . . , X I and Y 1 , . . . , Y I are sampled independently according to µ and ν. Note that for I = 100, like the one plotted in Figure 3, most of the samples already satisfy µ I ≤ cx ν I . As expected, the plotted points are close to the segment lines y = -x and y = 3x for x ∈ [-1, 1]. The corresponding cost is 0.5085, which is not so far from the theoretical cost 1/2 obtained for the continuous MOT. On 100 independent runs, we have obtained a mean cost of 0.5106 with a standard deviation equal to 0.0326. Now, we carry on our numerical study and consider lognormal distributions (Black and Scholes model). Namely, we assume that X i

(d) = exp σ X G -1 2 σ 2 X -1 and Y i (d) = exp σ Y G -1 2 σ 2 Y -1, with G ∼ N (0, 1
), σ X = 0.24 and σ Y = 0.28. We have estimated for I = 100 the probability that µ I ≤ cx ν I on 10 5 runs: we have obtained P( µ I ≤ cx ν I ) ≈ 0.4601, with a 95% confidence interval [0.4570, 0.4632]. Therefore, it is crucial to use either µ I ∧ ν I or µ I ∨ ν I on this example to recover discrete probability measures in the convex order. In Figure 4, we have plotted the points weighted by the solution to the discrete MOT minimization of c(x, y) = (x + y) + for ( µ I ∧ ν I , ν I ) and ( µ I , µ I ∨ ν I ) in a case where µ I ≤ cx ν I . Though being different, the two graphs are very similar. The costs are the same up to 9 digits (compared to 12 digits when µ I ≤ cx ν I ) and are equal to 0.1809. Thus, in our experiments, we have not observed any important differences between the MOT problems on ( µ I ∧ ν I , ν I ) and ( µ I , µ I ∨ ν I ), and we will work later on with ( µ I ∧ ν I , ν I ). Now, let us comment qualitatively Figure 4. We have plotted, as in Figure 3 the segment lines y = -x and y = 3x for x ∈ [-1, 1]. As one may expect, there are still many points on the first segment line, which is the boundary of the set of points where the cost function vanishes. Instead, the other segment line was more related to the specific laws that we have considered in Figure 3. In Figure 4, there are points close to y = 3x at the origin, but then they move away from this line. Points at the top right or at bottom left are mostly gathered around the line y = x, which is due to the martingale constraint. This is particularly noticeable on our example for the points at the bottom left.

6.3. Lognormal distributions with c(x, y) = |x-y| . We consider the cases = 2.1 and = 1.9, and take again µ and ν as the lognormal distributions of exp σ X G -1 2 σ 2 X -1 and exp σ Y G -1 2 σ 2 Y -1, with G ∼ N (0, 1), σ X = 0.24 and σ Y = 0.28. In Figure 5, we have plotted the weighted points of the discrete MOT for the minimization problem. We still observe for = 2.1 two curves that do not cross like in Figure 2, but the curves obtained for the lognormal distributions are quite far from the lines y = x + 1 and y = x -1 obtained for the uniform distribution. In Figure 6 we have plotted the weighted points for the maximization program. Interestingly, the optimal minimizing coupling for = 2.1 (resp. = 1.9) is very close to the optimal maximizing coupling for = 1.9 (resp. = 2.1). The graphs that we have produced up to now depend on the samples X 1 , . . . , X I and Y 1 , . . . , Y I . From a run to another, we observe the same patterns but with different points. One may like to have a discrete MOT with smaller statistical error. It is possible to do this by using Baker's construction on probability measures with finite supports. Namely, .

By Theorem 2.4.11 of Baker (2012), we have μI ≤ cx νI since µ I ∧ ν I ≤ cx ν I . We have plotted in Figure 7 the same discrete MOT as in Figure 5, but using (μ I , νI ) with I = 100 and I = 10000. The plot is clearly less noisy. This method can also be used to reduce variance on the calculation of the optimal cost. On the example of Figure 2, with 100 samples, I = 100 and I = 10000, we have a mean cost equal to 0.9981 with a standard deviation equal to 0.0148. This is to compare with the mean cost of 1.0020 and the standard deviation of 0.1400 obtained with ( µ I ∧ ν I , ν I ) and I = 100. 6.4. An example with three marginal laws. We first write the discrete MOT problem in this case. We have to minimize (or maximize) For a solution to exist, the measures µ = I i=1 p i δ x i , ν = J j=1 q j δ y j and η = K k=1 s k δ z k have to satisfy µ ≤ cx ν ≤ cx η.

For i = 1, . . . , I, we consider independent samples X i

(d) = exp σ X G -1 2 σ 2 X -1, Y i (d) = exp σ Y G -1 2 σ 2 Y -1 and Z i (d)
= exp σ Y G -1 2 σ 2 Y -1, with G ∼ N (0, 1), σ X = 0.24, σ Y = 0.28, σ Z = 0.32. In the financial application, the variables X + 1, Y + 1 and Z + 1 represent the values of an asset at three different dates t 1 < t 2 < t 3 in a Black-Scholes model, and we are interested in calculating price bounds for the option that pays Z -X+Y This has two advantages. First, as discussed in the example of Figure 7, the result is less noisy. Second, the dimension of r in the problem (6.4) is fixed and equal to (I ) 3 , while we do not know a priori the number of points weighted by µ I ∧ ( ν I ∧ η I ) and ν I ∧ η I , see Subsection 5.3. Thus, we keep the control on the time needed by the linear programming solver.

We have plotted in Figure 8 the points weighted by the discrete MOT that minimizes the cost. As one may expect, many points are on the hyperplane z = x+y 2 which is the boundary of the set of points where the cost function vanishes. With this example, we have obtained a cost of 0.0303, which is a lower bound for the price. Running the maximization program on the same sample, we have obtained 0.0856 for the price upper bound. As expected, the Black-Scholes price 0.0681 is between these bounds: the position of the Black-Scholes price with respect to the bounds is an indication on how conservative the model is to price the

  t

  ,y) = |x-y|^{2.3}

Figure 1 .

 1 Figure 1. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) with c(x, y) = |y -x| 2.3 for (µ I ∧ ν I , ν I ) at left and (µ I , µ I ∨ ν I ) at right, with I = 100. The segments y = x ± 1, x ∈ [-1, 1] indicate the support of π given by (6.1).

Figure 2 .

 2 Figure 2. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) with c(x, y) = |y -x| 2.3 for ( µ I ∧ ν I , ν I ) at left and ( µ I , µ I ∨ ν I ) at right, with I = 100. The segments y = x ± 1, x ∈ [-1, 1] indicate the support of π given by (6.1).

  Jensen's inequality is equivalent to have either π Y |X (x, [-x, +∞)) = 1 or π Y |X (x, (-∞, -x]) = 1, µ(dx)-a.e.. The martingale kernel π Y |X (x, dy) = 1 2 (δ -x (dy) + δ 3x (dy)) (6.2)

Figure 3 .

 3 Figure3. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) with c(x, y) = (y + x) + for ( µ I ∧ ν I , ν I ), with I = 100. The segments y = 3x and y = -x, x ∈ [-1, 1] indicate the support of π given by (6.2).

Figure 4 .

 4 Figure 4. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) with c(x, y) = (y + x) + for ( µ I ∧ ν I , ν I ) at left and ( µ I , µ I ∨ ν I ) at right, with I = 100. The segments y = 3x and y = -x, x ∈ [-1, 1] indicate the support of π given by (6.2).

Figure 5 .

 5 Figure5. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) for ( µ I ∧ ν I , ν I ) with c(x, y) = |y -x| 2.1 at left and c(x, y) = |y -x| 1.9 at right, with I = 100.

Figure 6 .

 6 Figure6. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT maximization (1.1) for ( µ I ∧ ν I , ν I ) with c(x, y) = |y -x| 2.1 at left and c(x, y) = |y -x| 1.9 at right, with I = 100.

rFigure 7 .r

 7 Figure7. Points (x i , y j ) such that r ij > 0, where r is the optimal solution of the MOT minimization (1.1) with c(x, y) = |y -x| 2.1 at left and c(x, y) = |y -x| 1.9 at right for (μ I , νI ) defined by (6.3) with I = 100 and I = 10000.

2+Figure 8 .

 8 Figure 8. Points (x i , y j , z k ) such that r ijk > 0, where r is the optimal solution of the MOT minimization (6.4) with c(x, y, z) = z -x+y 2 + for

  Jensen's inequality gives = 1). The equality condition in Jensen's inequality gives that |y -x| is constant π(dx, dy)-a.s. and thus equal to 1. The martingale condition then gives π = π . Let us note that for < 2, the same arguments with R×R |y -x| π(dx, dy) ≤ R×R |y -x| 2 π(dx, dy) 2 = 1 gives that π is the unique martingale coupling that maximizes R×R |y -x| π(dx, dy). The optimality of π for = 1 was obtained by

	We observe that
	π (dx, dy) =	1 2	1 [-1,1] (x)	δ x+1 (dy) + δ x-1 (dy) 2	dx		(6.1)
	achieves this lower bound. Let π ∈ Π M (µ, ν) be an optimal coupling (i.e.	R×R |y -
	x| π(dx, dy)
	|y -x| π(dx, dy) ≥	|y -x| 2 π(dx, dy)	2	= 1.
	R×R		R×R
	† https://www.gnu.org/software/glpk/

This research benefited from the support of the "Chaire Risques Financiers", Fondation du Risque, the French National Research Agency under the program ANR-12-BS01-0019 (STAB) and was completed after the hiring of Jacopo Corbetta by Zeliade Systems. 1

Algorithm 1 Calculate z 1 < z 2 < . . . < z K and ( 1 , . . . , K ) ∈ (0, 1] K s.t.

Remove the last element of L. end if if 1 = 0 then Remove the first element of L. end if option.