(Sharp) Bounds on $E(XYZ)$ given all the bivariate marginals Suppose $X,Y,Z$ are all real-valued random variables. Suppose I know the joint marginal distributions of $(X,Y)$, $(Y,Z)$ and $(X,Z)$. I want to find bounds on $E(XYZ)$.
In the case of bounding $E(XY)$ given marginals for $X$ and $Y$, the result is given by the Hardy-Littlewood inequalities, using symmetric decreasing/increasing rearrangements (see, for example, the version given here).
Translating the Hardy-Littlewood idea to this more general setting seems difficult.
I wonder if there is something possible if we assume that $Y$ and $Z$, say, have finite support. In this case, the Hardy-Littlewood result can give bounds on
$$E(XY|Z=z)~,$$
by applying the increasing decreasing rearrangements to the conditional distributions of $X$ and $Y$. But these bounds are not sharp, since the rearrangements need not respect the joint distribution of $(X,Y)$. Some kind of "constrained rearrangement" could work here - but I am not sure what that would look like...
EDIT: We can assume that the bivariate marginals are consistent with some (unknown) joint distribution. As pointed out in the comments, this has some non-trivial implications. A source that characterizes this is: https://perso.math.u-pem.fr/kloeckner.benoit/papiers/JointLaws.pdf
Another source that discusses some implications is Section 3.4.3 in "Multivariate Models and Dependence Concepts" by Harry Joe. Plus, some bounds for the joint distribution are given at the start of Section 3.4.1, in Theorem 3.11:

where $F_{ij}$ are the bivariate marginals. These bounds are not generally CDFs, however, and so, are not sharp. Nevertheless, maybe these bounds can still yield sharp bounds on $E(XYZ)$.
 A: Update: The bounds below are not sharp after all, and I now think that the only sharp bounds can be combinatorial. For example, in the notation below, let $p_{ac}=p_{ad}=p_{bc}=p_{bd}=1/4$, and:

*

*let $X_a$ have probability $1/4$ for each of $\pm1, \pm2$

*let $X_b$ have probability $1/4$ for each of $\pm3, \pm4$

*let $X_c$ have probability $1/4$ for each of $\pm1, \pm3$

*let $X_d$ have probability $1/4$ for each of $\pm2, \pm4$
Then in the bounds below, all of the $W$'s are zero, and all of the terms with integrals are positive. So the argument below gives a positive bound for $E[XYZ]$. But in fact $X_{ac}=\pm1$, $X_{ad}=\pm2$, $X_{bc}=\pm3$, $X_{bd}=\pm4$, so $E[XYZ]=0$.
The same argument also works if $X_a,X_b,X_c,X_d$ are continuous approximations to the above discrete distributions. In any case the bounds below are not sharp.
For simplicity, I'll analyze the case with $X$ arbitrary, $Y$ binary, $Z$ binary. (I'll also be happy to analyze one example of this sort, if someone can provide a good one.)
The inputs to the problem are four values, four probabilities and four quantile functions:

*

*The two possible values $a\ge b$ for $Y$, and the two possible values $c\ge d$ for $Z$.

*The probability $p_{ac}=P[Y=a,Z=c]$ and similarly $p_{ad}$, $p_{bc}$, $p_{bd}$.

*The quantile functions $Q_a,Q_b,Q_c,Q_d$ for $X|Y=a$, $X|Y=b$, $X|Z=c$, $X|Z=d$.

The outputs of the problem will be four further quantile functions:

*

*The quantile functions $Q_{ac},Q_{ad},Q_{bc},Q_{bd}$ for $X$ under joint hypotheses on $Y$ and $Z$.

We assume the inputs are consistent. We abbreviate $p_a=p_{ac}+p_{ad}$ and $W_a=p_aE[X|Y=a]=p_a\int_0^1 Q_a(q)dq$, etc., where the latter is a weighted expectation.
We want to maximize $E[XYZ]$, and we can represent that as
$$E[XYZ]=ac F_{ac}+ad F_{ad}+bcF_{bc}+bdF_{bd}$$
where the four variables $F_{ac}, F_{ad}, F_{bc}, F_{bd}$ represent
\begin{align}
F_{ac}=p_{ac}\,&E[X|Y=a,Z=c]\\
F_{ad}=p_{ad}\,&E[X|Y=a,Z=d]\\
F_{bc}=p_{bc}\,&E[X|Y=b,Z=c]\\
F_{bd}=p_{bd}\,&E[X|Y=b,Z=d]\\
\end{align}
These variables jointly satisfy four equalities
\begin{align}
F_{ac} + F_{ad} =W_a\\
F_{bc} + F_{bd} =W_b\\
F_{ac} + F_{bc} =W_c\\
F_{ad} + F_{bd} =W_d\\
\end{align}
(though only three equalities are independent) and each variable satisfies four inequalities, such as:
\begin{align}
p_a\int_{0}^{p_{ac}/p_a} Q_a(q)\,dq \ \le F_{ac}\  \le
p_a\int_{0}^{p_{ac}/p_a} Q_a(1-q)\,dq \\
p_c\int_{0}^{p_{ac}/p_c} Q_c(q)\,dq \ \le F_{ac}\  \le
p_c\int_{0}^{p_{ac}/p_c} Q_c(1-q)\,dq \\
\end{align}
Using the equalities, we can rewrite $E[XYZ]$
$%=a c F_{ac} + a d F_{ad} + b c F_{bc} + b d F_{bd}$
as:
\begin{align}
E[XYZ] %&=a c F_{ac} + ad(W_a-F_{ac})+bc(W_c-F_{ac}) + \frac{bd}{2}(W_b+W_d-W_a-W_c+2F_{ac})\\
&=adW_a+bcW_c+\frac{bd}{2}(W_b+W_d-W_a-W_c) +(a-b)(c-d)F_{ac} \\
%E[XYZ] &= ac(W_a-F_{ad})+ad F_{ad}+\frac{bc}{2}(W_b+W_c-W_a-W_d+2F_{ad}) + b d(W_d-F_{ad})\\
 &= acW_a+bdW_d+\frac{bc}{2}(W_b+W_c-W_a-W_d)+(a-b)(d-c)F_{ad}\\
%E[XYZ] &= a c(W_c-F_{bc}) + \frac{ad}{2}(W_a+W_d-W_b-W-c+2 F_{bc}) + b c F_{bc} + b d(W_b-F_{bc})$\\
&=acW_c + bdW_b + \frac{ad}{2}(W_a+W_d-W_b-W_c) +(b-a)(c-d)F_{bc}\\
%E[XYZ]&=\frac{a c}{2}(W_a+W_c-W_b-W_d+2F_{bd}) + a d(W_d-F_{bd}) + b c(W_b-F_{bd}) + b d F_{bd}\\
&=ad W_d + bc W_b +\frac{a c}{2}(W_a+W_c-W_b-W_d) + (b-a)(d-c)F_{bd}\\
\end{align}
So the maximum possibility for $E[XYZ]$ is the minimum of is bounded above by:
\begin{align}
\Big\{adW_a+bcW_c+\frac{bd}{2}(W_b+W_d-W_a-W_c) + (a-b)(c-d)&
p_a\int_{0}^{p_{ac}/p_a} Q_a(1-q)\,dq,\\
adW_a+bcW_c+\frac{bd}{2}(W_b+W_d-W_a-W_c) + (a-b)(c-d)&
p_c\int_{0}^{p_{ac}/p_c} Q_c(1-q)\,dq,\\
acW_a+bdW_d+\frac{bc}{2}(W_b+W_c-W_a-W_d)+(a-b)(d-c)&
p_a\int_{0}^{p_{ad}/p_a} Q_a(q)\,dq,\\
acW_a+bdW_d+\frac{bc}{2}(W_b+W_c-W_a-W_d)+(a-b)(d-c)&
p_d\int_{0}^{p_{ad}/p_a} Q_d(q)\,dq,\\
acW_c + bdW_b + \frac{ad}{2}(W_a+W_d-W_b-W_c) +(b-a)(c-d)&
p_b\int_{0}^{p_{bc}/p_b} Q_b(q)\,dq,\\
acW_c + bdW_b + \frac{ad}{2}(W_a+W_d-W_b-W_c) +(b-a)(c-d)&
p_c\int_{0}^{p_{bc}/p_c} Q_c(q)\,dq,\\
ad W_d + bc W_b +\frac{a c}{2}(W_a+W_c-W_b-W_d) + (b-a)(d-c)&
p_b\int_{0}^{p_{bd}/p_b} Q_b(1-q)\,dq,\\
ad W_d + bc W_b +\frac{a c}{2}(W_a+W_c-W_b-W_d) + (b-a)(d-c)&
p_d\int_{0}^{p_{bd}/p_d} Q_d(1-q)\,dq\Big\}\\
\end{align}
