Two dimensional Cantelli's inequality, $\Pr[X > 0, Y > 0]$ Cantelli's inequality says $\Pr[X > 0] \ge \frac{E[X^+]^2}{E[X^2]}$.
If $E[X]=0$ and $E[X^2]=1$, the fourth moment bound, $E|X|\ge\frac{E[X^2]^{3/2}}{E[X^4]^{1/2}}$, can give
$$\Pr[X > 0] \ge E[X^+]^2 = E\left[\frac{X+|X|}{2}\right]^2 = \frac14 E[|X|]^2 \ge \frac1{4E[X^4]}.$$
(With more care this can be improved to $\frac{2 \sqrt{3}-3}{E[X^4]}\approx\frac1{2.1547 E[X^4]}$.)
If we instead have two random variables $X,Y$, I'm wondering if there is a way to lower bound $\Pr[X > 0 \text{ and } Y > 0]$ using the moments of $(X, Y)$?
Berge's 2-dimensional inequality and the multivariate Chebyshev inequality both concern the norm $(X^2+Y^2)^{3/2}$, which treats positive and negative values the same.
If $(X,Y)$ was a sum of smaller increments, I could use Berry Esseen to bound $\Pr[X > 0 \text{ and } Y > 0]$ in terms of $E[(X^2+Y^2)^{3/2}]$.
Using the Cauchy Schwarz method from Cantelli's proof, we can easily show $\Pr[X>0, Y>0] \ge \frac{E[X^+Y^+]^2}{E[X^2Y^2]}$ (where $X^+ = X[X > 0]$), however unlike in the case of Cantelli, we can't just say $E[X^+Y^+] \ge E[XY]$, since $E[XY] = E[X^+Y^+ - X^+Y^- - X^-Y^+ + X^-Y^-]$ and $E[X^-Y^-]$ could perhaps be big.
Edit: I expect the marginals $X$ and $Y$ to be close to symmetric, so $\Pr[X>0]$ and $\Pr[Y>0]$ are not larger than 1/2 individually. This prevents me from using a union bound.
 A: I found the book "Probability Inequalities in Multivariate Distributions" by Y. L. Tong.
In chapter 7 multiple results are given, one of which is the following:
Assume $E[X]=E[Y]=0$ and $E[X^2]=E[Y^2]=1$.
Let $\rho=E[XY]\in(-1,1)$, then
$$\Pr[X>-1,Y>-1]
\ge 1-2\big(\frac{1+\sqrt{1-\rho}+\rho}{3+\rho}\big)^2
\ge\frac19+\frac{4\rho}{27}.
$$
That's the kind of result I was looking for, except I would have preferred $>0$ rather than $>-1$.
However, it's not too hard to come up with examples showing that won't work as long as we're not using higher moments of $X$ and $Y$.
I also realized that one can expand the proof by Cantelli as follows:
Let $X^+ = X[X>0]$ be the positive part of a random variable.
Now, using Cauchy Schwarz:
$$
E[X^+Y^+]
=E[XY[X>0,Y>0]]
\le E[(XY)^2]^{1/2}\Pr[X>0,Y>0]^{1/2}
\implies \Pr[X>0,Y>0]\ge\frac{E[X^+Y^+]^2}{E[(XY)^2]}.
$$
Unfortunately It's not clear to me how to bound $E[X^+Y^+]$ using higher moments of $X$, $Y$ and $XY$.
A: There are several lower bounds one could use. For instance, we have for any two events $A$ and $B$, 
\begin{equation}
P(A \cap B) \ge P(A) + P(B) - 1
\end{equation}
In our case, $A = \{X > 0\}$ and $B = \{Y > 0\}$.
You can then apply the Cantelli inequality for each term. If $A$ and $B$ are independent, you do get a more useful bound, obviously.
A: Edit: The answer below is probably broken in multiple ways, so please take it lightly.
However, I still conjecture that when $X$ and $Y$ matches the moments of gaussians with $E[XY]=\rho$ (and so the fourth moments are $E[X^2Y^2]=1+2\rho^2$,  $E[X^3Y]=3\rho$ and $E[X^4]=3$), we have $\Pr[X>0\wedge Y>0]\sim2\rho^4$.
In particular $\Pr[X>0\wedge Y>0]>0$ whenever $\rho>0$.
I would love to see a correct proof of this.

Another approach:
Let $Z=X+Y$ and $T=X-Y$,
then
$$\begin{align}
\Pr[X>0 \wedge Y>0]
&=\Pr[Z>0 \wedge  |T|<Z]
\\&=\int[z>0]\Pr[|T|<z]dP(z)
\\&\ge\int[z > s]\left(1-\frac{s^2}{z^2}\right)dP(z),
\end{align}$$
By Chebyshev's inequality, where $s^2=E[T^2]=E[X^2]+E[Y^2]-2E[XY]=2-2\rho$.
Now let define
$$\begin{align}
Q_{l,r}(x)=\frac{(l-x)^2 (x-s) \left(
r (2 l r-4 r^2+3 l s-5 r s) +( -l (r+2 s)+r (3 r+4 s))x 
\right)}{r^3 (l-r)^3}
\end{align},$$
where $l,r\in\mathbb R$.
Then we have
$$[z > s]\left(1-\frac{s^2}{z^2}\right) \ge Q_{l,r}(x).$$
For all $x$ and $l<s<r$.
See the picture for how $Q$ embraces the function.

(In this case $l=0$, $s=1$ and $r=2$.)
We thus have
$$\int[z > s]\left(1-\frac{s^2}{z^2}\right)dP(z)
\ge E[Q_{l,r}(Z)].$$
Optimizing over $l$ and $r$ we get some ugly expression in $s$, $\rho$, $\kappa=E[(X+Y)^3]$ and $\gamma=E[(X+Y)^4]$.
The resulting formula is a bit wildy.
But if I didn't make a mistake, the result equals
$$\begin{align}
\frac{\left(8 \rho^2+\sqrt{2} \sqrt{1-\rho} \gamma +8 \rho\right)^4}{4 (\rho-1)^2 \gamma^3}+O\left(\left(\frac{1}{\gamma}\right)^4\right)
\end{align}$$
for large $\gamma$.
In the "normal" case, $\gamma= 12 (1 + 2 \rho + \rho^2)$ and $\kappa= 0$
we get
$$\frac{16 \rho^4 (2 \rho+1) \sqrt{\frac{\rho+1}{\rho (7 \rho+4)+1}}}{\left(\rho \sqrt{1-\rho}+\sqrt{1-\rho}+\sqrt{(\rho+1) (\rho (7 \rho+4)+1)}\right)^3}=2 \rho^4-8 \rho^5+O\left(\rho^6\right).$$
The function looks like this:

I wonder if this $\Pr[X>0,Y>0]\sim \rho^4$ is the best one can hope for.
Again in the case where $X$ and $Y$ are normal distributed, with the same moments as above, we have $\Pr[X>0,Y>0]=\frac{2 \tan ^{-1}\left(\frac{\rho}{\sqrt{1-\rho^2}}\right)+\pi }{4 \pi }=\frac{1}{4}+\frac{\rho}{2\pi}$, which is of course much better.
I tried solving the problem numerically.
That is,

Maximize
$$\left\langle \left[\begin{smallmatrix}
1 & E[X] & E[X^2] & E[X^3] & E[X^4] \\
E[Y] & E[XY] & E[X^2Y] & E[X^3Y] & 0 \\
E[Y^2] & E[XY^2] & E[X^2Y^2] & 0 & 0 \\
E[Y^3] & E[XY^3] & 0 & 0 & 0 \\
E[Y^4] & 0 & 0 & 0 & 0
\end{smallmatrix}\right], M \right\rangle$$ over matrices $M\in\mathbb   R^{4\times 4}$ such
that for all $x,y\in\mathbb R$: $[1, x, x^2, x^3, x^4]^T M [1, y, y^2, y^3, y^4] \le [x > 0 \wedge y > 0]$.

To handle the $\forall$ I sampled 10,000 pairs $x,y\sim N(0,5)$ and solved the resulting linear program.
The result was the following curve:

Here the moments $E[X^i Y^j]$ are again taken from the corresponding multi-normal with covariance $[\begin{smallmatrix}1 & \rho\\\rho & 1\end{smallmatrix}]$.
The x-axis is $\rho$ and the y-axis is the probability lower bound.
This plots seems to indicate that the analytical bound isn't that far off from what's theoretically possibly to say with access only to the given moments.
To get a bound $>0$ even when $\rho=0$ or negative, it seems we need to use sixth moments. It seems (again from numerical experiments) that in this case we can get $\Pr[X>0,Y>0] \ge c + \Omega(\rho)$ for some small $c\approx 0.02$:

