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}