Compute $\mathbb{E}[\text{max}(x,0) \text{max}(y,0)]$ where $(x,y)$ is jointly gaussian with given covariance and nonzero mean 
Is there a closed form expression for the expectation of $g(x,y) = \text{max}(x,0) \text{max}(y,0)$ where $(x,y)$ is jointly gaussian with the bivariate gaussian distribution, when the mean is not zero and the covariance is non singular?

I tried to apply Price's theorem that can be written in the following form:
$$\frac{\partial^2 \mathbb{E}[g(x,y)]}{\partial \rho^2} = \mathbb{E} \left[ \frac{\partial^4 g(x,y)}{\partial x^2 \partial y^2} \right] = \mathbb{E}[\delta(x) \delta(y)] = \frac{1}{2 \pi \sigma_x \sigma_y\sqrt{1 - \rho^2}} \text{exp}\left(-\frac{1}{2(1 - \rho^2)} \left(\frac{\mu_x^2}{\sigma_x^2} +  \frac{\mu_y^2}{\sigma_y^2} - \frac{2\rho \mu_x\mu_y}{\sigma_x \sigma_y} \right)\right)$$
Unfortunately, integrating the last equality with respect to the correlation coefficient doesn't seem to have a closed form in terms of known functions. I have posted the integral question earlier before here.
Moreover, it doesn't seem to be possible to find the expectation directly as:
$$\int_0^{\infty} \int_0^{\infty} xy f(x,y) dx dy$$
EDIT: Price's theorem can nicely handle the zero mean case with an elegant solution.
 A: The calculation in question is pretty straightforward and can be done using elementary methods. Assume for simplicity that the variables have variance one and mean zero. Then the joint pdf reads:
\begin{eqnarray}
\rho(x,y) = \frac{1}{2\pi \sqrt{1-\rho^2}} \exp\left[ -\frac{1}{2} \frac{1}{(1-\rho^2)} (x^2+y^2-2 \rho x y )\right]
\end{eqnarray}
We integrate over $x$ first.
\begin{eqnarray}
I(y):=\int\limits_0^\infty x \rho(x,y) dx &=& \frac{1}{2\pi \sqrt{1-\rho^2}} \int\limits_0^\infty \exp\left[-\frac{1}{2} \frac{1}{(1-\rho^2)} ((x-\rho y)^2 + y^2(1-\rho^2))\right]\\
&=&\frac{1}{2\pi \sqrt{1-\rho^2}}e^{-\frac{1}{2} y^2} \int\limits_{-\rho y}^\infty (x+ \rho y) e^{-\frac{1}{2} \frac{1}{1-\rho^2} x^2} dx \\
&=& \frac{1}{2\pi} \sqrt{1-\rho^2} e^{-\frac{1}{2} y^2 \frac{1}{1-\rho^2}} + \frac{\rho  e^{-\frac{y^2}{2}} y}{2 \sqrt{2 \pi }} + \frac{\rho  e^{-\frac{y^2}{2}} y \text{erf}\left(\frac{\rho  y}{\sqrt{2-2 \rho ^2}}\right)}{2 \sqrt{2 \pi }}
\end{eqnarray}
Now we multiply the result by $y$ and integrate the whole thing over $y\in(0,\infty)$. Even at the first glance it is clear that the integrals from the first two terms on the right hand side are doable it is only the last integral that might cause difficulties. However even that integral is doable and it reads:
\begin{equation}
\int\limits_0^\infty y^2 e^{-\frac{1}{2}y^2} Erf[a y] dy = \frac{1}{\sqrt{\pi}} \left[ \frac{2 a}{1+2 a^2} + \sqrt{2} \arctan(\sqrt{2} a)\right]
\end{equation}
The result was derived by differentiating with respect to the parameter $a$.
The final result is as follows:
\begin{eqnarray}
\left< max(x,0) max(y,0) \right> = \frac{\left(1-\rho ^2\right)^{3/2}}{2 \pi } + \frac{\rho }{4} + \frac{\rho  \left(\sqrt{1-\rho ^2} \rho +\tan ^{-1}\left(\sqrt{\frac{\rho ^2}{1-\rho ^2}}\right)\right)}{2 \pi }
\end{eqnarray}
rho =.;
myrho[x_, y_] := 
  1/(2 Pi Sqrt[1 - rho^2]) Exp[-1/
      2 1/(1 - rho^2) (x^2 + y^2 - 2 rho x y)];
rho = RandomReal[{0, 1}, WorkingPrecision -> 50];
NIntegrate[x y myrho[x, y], {x, 0, Infinity}, {y, 0, Infinity}, 
 WorkingPrecision -> 20]
(1 - rho^2)^(3/2)/(2 \[Pi]) + rho/4 + (
 rho (rho Sqrt[1 - rho^2] + ArcTan[Sqrt[rho^2/(1 - rho^2)]]))/(
 2 \[Pi])


