The McCormick martingale optimal transport

Under a risk-neutral measure $\pi$ with discrete-time processes, the underlying assets $X$ and $Y$ are typically represented as:

$$X_{t+1}-X_{t}=X_{t}W_{t+1},\qquad Y_{t+1}-Y_{t}=Y_{t}Z_{t+1}.$$

(3.1)

Here, $\{W_{t}\}^{N-1}_{t=1}$ and $\{Z_{t}\}^{N-1}_{t=1}$ denote random returns for $X$ and $Y$, respectively, which are the only sources of randomness in the system (3.1). A standard assumption posits that $\{W_{t}\}^{N-1}_{t=1}$ are independent and identically distributed, as commonly imposed in binomial tree models. $\{Z_{t}\}^{N-1}_{t=1}$ also satisfy the same assumption. However, it should be noted that for each time $1\leq t\leq N-1$, $W_{t}$ and $Z_{t}$ are allowed to be correlated. For example, a parsimonious assumption is that the correlation between $W_{t}$ and $Z_{t}$ is a constant. Under these conditions, $X_{t+1}$ is independent of $Y_{t}$, conditional on $X_{t}$. This equivalently implies that $\pi(dx_{t+1}|x_{t},y_{t})=\pi(dx_{t+1}|x_{t})$, since knowledge of $Y_{t}$ does not aid in predicting $W_{t+1}$ once $X_{t}$ is known. It is crucial to emphasize that $Y_{t}$ and $X_{t+1}$ are not independent if $X_{t}$ is unknown.

Aligned with the nonanticipativity condition that originates from stochastic difference (or differential) equations, Lassalle, (2013) imposes that a transport plan $\pi$ should satisfy

$$\pi(dx_{t+1}|x_{1:t},y_{1:t})=\pi(dx_{t+1}|x_{1:t}),\quad t=1,...,N-1,\quad\text{$\pi$-a.s.},$$

(3.2)

which is equivalent to

$$\pi(dy_{t}|x_{1:N})=\pi(dy_{t}|x_{1:t}),\quad t=1,...,N-1,\quad\text{$\pi$-a.s.},$$

(3.3)

see Eckstein and Pammer, (2024, Remark 2.2). The property (3.2) (or (3.3)) is known as the causality condition. A transport plan that satisfies (3.2) (or (3.3)) is termed causal. If we interchange the positions of $X$ and $Y$ the following condition holds:

$$\pi(dx_{t}|y_{1:N})=\pi(dx_{t}|y_{1:t}),\quad t=1,...,N-1,\quad\text{$\pi$-a.s.},$$

(3.4)

then the transport plan is anticausal. Couplings that are both causal and anticausal are referred to as bicausal. For later use, denote ${\mathcal{M}}_{bc}(\bar{\mu},\bar{\nu})\subset{\mathcal{M}}(\bar{\mu},\bar{\nu})$ as the set of all bicausal and martingale transport plans, with $\mu_{t}$ and $\nu_{t}$ as marginals for each $t$.

Motivated by the nonanticipativity condition, we propose to further restrict the MOT problem (2.2) by considering bicausal plans:

$$\sup_{\pi\in{\mathcal{M}}_{bc}(\bar{\mu},\bar{\nu})}\int c(x_{1:N},y_{1:N})\pi(dx_{1:N},dy_{1:N})\,\text{ and }\inf_{\pi\in{\mathcal{M}}_{bc}(\bar{\mu},\bar{\nu})}\int c(x_{1:N},y_{1:N})\pi(dx_{1:N},dy_{1:N}).$$

(3.5)

Before delving into the main results, we discuss some properties for the bicausal MOT (3.5). The first consequence of the bicausal condition is the equivalence of the martingale property under individual filtrations and the joint filtration. The proof of Proposition 3.1 directly applies the causality condition (3.2) and is therefore omitted.

When the joint distributions and, consequently, conditional kernels $\mu(dx_{t+1}|x_{1:t})$ are provided, the causality condition (3.2) (or (3.3)) imposes a linear constraint on $\pi$; see Backhoff-Veraguas et al., (2017, Proposition 2.4 (3)). However, this is no longer applicable when the conditional kernels are not predetermined. The corresponding equivalent condition is outlined in Proposition 3.2. The proof follows a rationale similar to that of Backhoff-Veraguas et al., (2017, Proposition 2.4 (3)), and we include it in the Appendix for the sake of completeness.

To observe that the causality constraint is no longer linear in $\pi$, one can think about the discrete case. The conditional mass function $\pi(x_{t+1:N}|x_{1:t})=\pi(x_{1:t},x_{t+1:N})/\pi(x_{1:t})$, where the denominator $\pi(x_{1:t})>0$ is not uniquely determined. In the subsequent section, we demonstrate that this constraint is bilinear in $\pi$ and non-convex. This characteristic poses a primary challenge in the bicausal MOT and makes the program difficult to solve. Consequently, we propose adopting the McCormick relaxation (McCormick,, 1976) for the bilinear constraint. The resulting relaxed problem is convex, offering a computationally more tractable approach while maintaining the capacity to enhance bounds in various scenarios.

To elucidate the main concept, we examine the case with discrete marginals, interpreting $\pi(x_{1:N},y_{1:N})$ as probability mass functions. The causality constraint (3.2) (or (3.3)) can be expressed equivalently as

$$\pi(x_{1:N},y_{t})\pi(x_{1:t})=\pi(x_{1:t},y_{t})\pi(x_{1:N}),\quad t=1,\ldots,N-1.$$

(4.1)

We stress that (4.1) holds pointwise, while (3.2) (or (3.3)) holds $\pi$-a.s. Indeed, if $\pi(x_{1:N})>0$, then $\pi(x_{1:t})>0$. Both conditional probabilities in (3.3) are well-defined. Then (4.1) holds. If $\pi(x_{1:N})=0$, then $\pi(dy_{t}|x_{1:N})$ is not defined. (3.3) is not required to hold. However, we note that $\pi(x_{1:N},y_{t})=0$ since $\pi(x_{1:N})=0$. Therefore, (4.1) is valid as $\pi(x_{1:N},y_{t})\pi(x_{1:t})=\pi(x_{1:t},y_{t})\pi(x_{1:N})=0$.

(4.1) essentially constitutes a bilinear equality constraint on $\pi$. Consequently, it is computationally challenging to find exact upper/lower bounds in the bicausal MOT problem (3.5). In the optimization literature, bilinear programming falls under quadratically constrained quadratic programming (QCQP). Several algorithms, including branch and cut with reformulation-linearization techniques (Audet et al.,, 2000), and semi-definite programming relaxations (Anstreicher,, 2009; Sherali and Fraticelli,, 2002), have been proposed. In this study, we opt for the classic McCormick relaxation in McCormick, (1976) for several reasons. It is widely employed as a built-in routine in optimization software and is also straightforward to implement manually. Moreover, the dual formulation remains tractable, allowing for the demonstration of the impact of causality constraints on hedging strategies.

McCormick relaxations construct the convex and concave envelopes as follows. The bilinear term on the left-hand side of (4.1) is substituted with $w_{1}(x_{1:N},y_{t})=\pi(x_{1:N},y_{t})\pi(x_{1:t})$. For probability $\pi(x_{1:N},y_{t})$, we suppose the upper bound $U_{N,t}(x_{1:N},y_{t})$ and lower bound $L_{N,t}(x_{1:N},y_{t})$ are known. The corresponding upper and lower bounds for $\pi(x_{1:t})$ are denoted as $U_{t,0}(x_{1:t})$ and $L_{t,0}(x_{1:t})$, respectively. Here, the subscripts in $L_{t,s}$ and $U_{t,s}$, $0\leq t,s\leq N$, indicate that different functions are used when the variables $(x,y)$ have different indices. However, for the sake of brevity, we denote $L_{t,s}=:L$ and $U_{t,s}=:U$, since the indices in $(x,y)$ imply the functions used. Then we simply write

$\displaystyle L(x_{1:N},y_{t})$

$\displaystyle\leq\pi(x_{1:N},y_{t})\leq U(x_{1:N},y_{t})\quad\text{ and }\quad L(x_{1:t})\leq\pi(x_{1:t})\leq U(x_{1:t}).$

Since $(\pi(x_{1:N},y_{t})-L(x_{1:N},y_{t}))(\pi(x_{1:t})-L(x_{1:t}))\geq 0$, expanding the expression gives

$$w_{1}(x_{1:N},y_{t})\geq L(x_{1:N},y_{t})\pi(x_{1:t})+\pi(x_{1:N},y_{t})L(x_{1:t})-L(x_{1:N},y_{t})L(x_{1:t}).$$

(4.2)

Similarly, the same procedure applied to products involving $U(x_{1:N},y_{t})-\pi(x_{1:N},y_{t})$ and $U(x_{1:t})-\pi(x_{1:t})$ yields three additional inequalities. We also replace the right-hand side of (4.1) with $w_{2}(x_{1:N},y_{t})=\pi(x_{1:t},y_{t})\pi(x_{1:N})$, and four corresponding inequalities are derived similarly. The causality constraint becomes $w_{1}(x_{1:N},y_{t})=w_{2}(x_{1:N},y_{t})$.

Together with the relaxation on anticausality and the martingale condition, we introduce McCormick martingale couplings as follows:

All of the inequalities in (4.4) are derived in the same way as (4.2). Note that (4.4) holds for all $x_{1:N},y_{t}$, instead of merely $\pi$-a.s. In condition (3), it is possible to simplify $w_{1}(x_{1:N},y_{t})=w_{2}(x_{1:N},y_{t})=:w(x_{1:N},y_{t})$. The term $w(\cdot)$ (or $w_{1}(\cdot)$ and $w_{2}(\cdot)$) also serves as a variable in the optimization process. However, in (4.4), we can eliminate the dependence on $w(\cdot)$ by imposing that the maximum of the lower bounds does not exceed the minimum of the upper bounds for $w(\cdot)$. Therefore, we only write $\pi\in{\mathcal{M}}(\bar{\mu},\bar{\nu};L,U)$ to signify that $\pi$ is a McCormick martingale coupling.

A simple choice of bounds is $L=0$ and $U(x_{1:N},y_{t})=\min\{\mu_{1}(x_{1}),\ldots,\mu_{N}(x_{N}),\nu_{t}(y_{t})\}$, which are derived from the marginal conditions and do not require additional modeling assumptions. Under this choice, ${\mathcal{M}}(\bar{\mu},\bar{\nu};L,U)$ is also not empty. Indeed, we can first construct (or rather arbitrarily pick) the marginal laws, $\mu\in\Pi(\overline{\mu})$ and $\nu\in\Pi(\overline{\nu})$, such that $X_{1:N}$ and $Y_{1:N}$ are martingales under $\mu$ and $\nu$, respectively. Then the independent coupling $\mu\otimes\nu$ serves as a McCormick martingale coupling. In the general finite and discrete cases, many LP solvers can check whether ${\mathcal{M}}(\bar{\mu},\bar{\nu};L,U)$ is non-empty.

We focus on addressing the minimization problem associated with the McCormick MOT in (4.5), as the maximization problem can be treated similarly:

$$\inf_{\pi\in{\mathcal{M}}(\bar{\mu},\bar{\nu};L,U)}\int c(x_{1:N},y_{1:N})\pi(dx_{1:N},dy_{1:N}).$$

(4.5)

In finite and discrete scenarios, the infimum in the primal problem is attained when the set ${\mathcal{M}}(\bar{\mu},\bar{\nu};L,U)$ is not empty. The dual problem is derived through classic finite-dimensional LP and is omitted here.

In the continuous scenario similar to Korman and McCann, (2015); Bogachev, (2022), we make the assumption that the domains ${\mathcal{X}}_{t}={\mathcal{Y}}_{t}=\mathbb{R}$ and the marginals $\mu_{t}$ and $\nu_{t}$ are absolutely continuous with respect to the Lebesgue measure on $\mathbb{R}$. Similarly, we consider couplings $\pi$ that are absolutely continuous with respect to the Lebesgue measure $\lambda$ on $\mathbb{R}^{N}\times\mathbb{R}^{N}$. When the context is evident, we denote $\lambda(dx_{1:N},dy_{1:N})$ by $d\lambda$ and omit the arguments. Let $f(x_{1:N},y_{1:N})$ represent the density of a coupling $\pi$. Although probability measures with densities are absolutely continuous with respect to the Lebesgue measure, we restrict them to the Borel $\sigma$-algebra on the Euclidean space and interpret them as Borel probability measures. Additionally, nonnegative stock prices result in zero densities in negative domains.

Due to absolute continuity, the causality constraint (3.2) (or (3.3)) transforms into the following density equality:

$$f(x_{1:N},y_{t})f(x_{1:t})=f(x_{1:t},y_{t})f(x_{1:N}),\quad\text{$\lambda$-a.e.},\quad t=1,\ldots,N-1.$$

(4.6)

Crucially, (4.6) holds $\lambda$-a.e., instead of merely $\pi$-a.s. It can be verified in the same spirit of (4.1), by noting that both sides of (4.6) are zero when regular conditional kernels are not defined.

With a slight abuse of notation, let $w_{1}(x_{1:N},y_{t})=f(x_{1:N},y_{t})f(x_{1:t})$ still be the auxiliary variable. With the same convention on omitting the subscripts as in the discrete case, we suppose that the upper bound $u(\cdot)$ and lower bound $l(\cdot)$ for the density $f(\cdot)$ are known, such that

$\displaystyle l(x_{1:N},y_{t})$

$\displaystyle\leq f(x_{1:N},y_{t})\leq u(x_{1:N},y_{t}),\;\text{$\lambda$-a.e.}\quad\text{ and }\quad l(x_{1:t})\leq f(x_{1:t})\leq u(x_{1:t}),\;\text{$\lambda$-a.e.},$

which are referred to as capacity/density constraints (Korman and McCann,, 2015; Bogachev,, 2022). Similarly, the inequality (4.2) becomes

$$w_{1}(x_{1:N},y_{t})\geq l(x_{1:N},y_{t})f(x_{1:t})+f(x_{1:N},y_{t})l(x_{1:t})-l(x_{1:N},y_{t})l(x_{1:t}),\quad\text{$\lambda$-a.e.}$$

(4.7)

Other inequalities can be derived in a similar fashion. It is important to emphasize that $l(\cdot)$ and $u(\cdot)$ are not necessarily density functions.

Denote $g_{t}$ as the density of $\mu_{t}$, and $\bar{g}:=(g_{1},\ldots,g_{N})$ as the individual density vectors. Similarly, let $h_{t}$ be the density of $\nu_{t}$, and $\bar{h}:=(h_{1},\ldots,h_{N})$ be the individual density vectors. We define the continuous version of McCormick martingale couplings similarly. Note that (4.9) holds $\lambda$-a.e.

With continuous densities, we denote the risk-neutral price of the exotic option as

$$P(f):=\int c(x_{1:N},y_{1:N})f(x_{1:N},y_{1:N})d\lambda.$$

(4.10)

The primal problem, incorporating McCormick relaxation and capacity constraints, can be expressed as

$$P=\inf_{f\in{\mathcal{M}}(\bar{g},\bar{h};l,u)}P(f).$$

(4.11)

This formulation is referred to as the continuous version of McCormick MOT.

Theorem 4.5 demonstrates the primal attainment of the problem (4.11).