Compute $\mathbb{E}[\text{max}(x,0) \text{max}(y,0)]$ where $(x,y)$ is jointly gaussian with given covariance and nonzero mean

Question

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.

Answer 1

Since a claim has been made that the general case is somehow beyond reach I decided to post it below. I am not going to present all the details of the computation because it took me half of a day to carry out it in Mathematica and I have not saved the intermediate steps. Yet let me assure you that I only used integration by parts and definitions of special functions. In some cases certain integrals were not reducible to known special functions so I just left those terms aside. Later it appeared that even those hard terms are all expressible through a generalized Owen's T function Generalized Owen's T function which is defined as a joint Gaussian probability in the following way: \begin{equation} T(h,a,b):= {\bf P}\left(X>h \quad \wedge \quad a X+b > Y > 0 \left. \right| X = N(0,1) , Y=N(0,1) \right) \end{equation}

Back to our question. Clearly the expectation in question equals: \begin{eqnarray} &&{\mathbf E}\left[\mbox{max}(X,0)\mbox{max}(Y,0)\right]= \int\limits_{-\frac{\mu_1}{s_1}}^\infty \int\limits_{-\frac{\mu_2}{s_2}}^\infty \left(\mu_1+s_1 \xi\right)\left(\mu_2+s_2 \eta\right) \cdot \rho(\xi,\eta) d\xi d\eta \end{eqnarray} where \begin{equation} \rho(\xi,\eta) := \frac{1}{2\pi \sqrt{1-\rho^2}} \exp\left[ -\frac{1}{2} \frac{1}{1-\rho^2} \left( \xi^2+\eta^2 - 2 \rho \xi \eta\right) \right] \end{equation} is the joint pdf of a bivariate Gaussian distribution with means zero and variances one. Now the onle thing we need to do is to compute four two dimensional integrals being expectations of $1$,$\xi$,$\eta$ and of $\xi \cdot \eta$. Here I only state the result and then as usual verify it numerically. We have: \begin{eqnarray} {\mathbf E}\left[\mbox{max}(X,0)\mbox{max}(Y,0)\right]= \mu_1 \mu_2 I_{0,0} + \mu_1 s_2 I_{0,1} + \mu_2 s_1 I_{1,0} + s_1 s_2 I_{1,1} \end{eqnarray} where \begin{eqnarray} &&I_{0,0}:=\int\limits_{-\frac{\mu_1}{s_1}}^\infty \int\limits_{-\frac{\mu_2}{s_2}}^\infty 1 \cdot \rho(\xi,\eta) d\xi d\eta=\\ &&T\left(-\frac{\mu_2}{s_2},\frac{\rho }{\sqrt{1-\rho ^2}},\frac{\mu_1}{\sqrt{1-\rho ^2} s_1}\right)+\frac{1}{4} \left(\text{erf}\left(\frac{\mu_2}{\sqrt{2} s_2}\right)+1\right)\\ &&I_{0,1}:=\int\limits_{-\frac{\mu_1}{s_1}}^\infty \int\limits_{-\frac{\mu_2}{s_2}}^\infty \eta \cdot \rho(\xi,\eta) d\xi d\eta=\\ &&\frac{e^{-\frac{\mu_2^2}{2 s_2^2}} \left(\text{erf}\left(\frac{\mu_1 s_2-\mu_2 \rho s_1}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)+1\right)+\rho e^{-\frac{\mu_1^2}{2 s_1^2}} \text{erfc}\left(\frac{\mu_1 \rho s_2-\mu_2 s_1}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)}{2 \sqrt{2 \pi }}\\ &&I_{1,0}:=\int\limits_{-\frac{\mu_1}{s_1}}^\infty \int\limits_{-\frac{\mu_2}{s_2}}^\infty \xi \cdot \rho(\xi,\eta) d\xi d\eta=\\ &&\frac{e^{-\frac{\mu_1^2}{2 s_1^2}} \left(\text{erf}\left(\frac{\mu_2 s_1-\mu_1 \rho s_2}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)+1\right)+\rho e^{-\frac{\mu_2^2}{2 s_2^2}} \text{erfc}\left(\frac{\mu_2 \rho s_1-\mu_1 s_2}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)}{2 \sqrt{2 \pi }}\\ &&I_{1,1}:=\int\limits_{-\frac{\mu_1}{s_1}}^\infty \int\limits_{-\frac{\mu_2}{s_2}}^\infty \xi \eta \cdot \rho(\xi,\eta) d\xi d\eta=\\ &&-\rho \cdot T\left(-\frac{\mu_2}{s_2},-\frac{\rho }{\sqrt{1-\rho ^2}},-\frac{\mu_1}{\sqrt{1-\rho ^2} s_1}\right)-\frac{\mu_1 \rho e^{-\frac{\mu_1^2}{2 s_1^2}} \left(\text{erf}\left(\frac{\mu_2 s_1-\mu_1 \rho s_2}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)+1\right)}{2 \sqrt{2 \pi } s_1}+\frac{1}{4} \rho \left(\text{erf}\left(\frac{\mu_2}{\sqrt{2} s_2}\right)+1\right)-\frac{\mu_2 \rho e^{-\frac{\mu_2^2}{2 s_2^2}} \text{erfc}\left(\frac{\mu_2 \rho s_1-\mu_1 s_2}{\sqrt{2-2 \rho ^2} s_1 s_2}\right)}{2 \sqrt{2 \pi } s_2}+\frac{\sqrt{1-\rho ^2} \exp \left(\frac{\mu_1^2 s_2^2-2 \mu_1 \mu_2 \rho s_1 s_2+\mu_2^2 s_1^2}{2 \left(\rho ^2-1\right) s_1^2 s_2^2}\right)}{2 \pi } \end{eqnarray}

rho =.; mu1 =.; mu2 =.; s1 =.; s2 =.;
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)];
T[h_, a_, b_] := 
  NIntegrate[(E^(-(b^2/2) - xi b h - 1/2 (1 + xi^2) h^2)) /(
     2 (1 + xi^2) \[Pi]) - 
     b  /(2 Sqrt[2] Sqrt[ \[Pi]]) (
      xi  Erfc[(h + xi (b + xi h))/(Sqrt[2] Sqrt[1 + xi^2])])/ ((1 + 
        xi^2)^(3/2)) E^(-(b^2/(2 + 2 xi^2))), {xi, 0, a}, 
    WorkingPrecision -> 20] + Erfc[h/Sqrt[2]] Erf[b/Sqrt[2]] 1/4;

{mu1, mu2, s1, s2} = RandomReal[{0, 1}, 4, WorkingPrecision -> 50];
rho = RandomReal[{0, Sqrt[s1 s2]}, WorkingPrecision -> 50];

l1 = NIntegrate[{1, y, x, x y} myrho[x, y], {x, -mu1/s1, 
    Infinity}, {y, -mu2/s2, Infinity}];


l2 = {1/4 (1 + Erf[mu2/(Sqrt[2] s2)]) + 
    T[-(mu2/s2), rho/Sqrt[1 - rho^2], mu1/(Sqrt[1 - rho^2] s1)],
   1/(2 Sqrt[
     2 \[Pi]]) (E^(-(mu2^2/(
        2 s2^2))) (1 + 
         Erf[(-mu2 rho s1 + mu1 s2)/(Sqrt[2 - 2 rho^2] s1 s2)]) + 
      E^(-(mu1^2/(2 s1^2)))
        rho Erfc[ (-mu2 s1 + mu1 rho s2)/(Sqrt[2 - 2 rho^2] s1 s2)]),
   1/(2 Sqrt[
     2 \[Pi]]) (E^(-(mu1^2/(
        2 s1^2))) (1 + 
         Erf[(-mu1 rho s2 + mu2 s1)/(Sqrt[2 - 2 rho^2] s1 s2)]) + 
      E^(-(mu2^2/(2 s2^2)))
        rho Erfc[ (-mu1 s2 + mu2 rho s1)/(Sqrt[2 - 2 rho^2] s1 s2)]),
    Sqrt[1 - rho^2]/(2 \[Pi]) E^((
     mu2^2 s1^2 - 2 mu1 mu2 rho s1 s2 + mu1^2 s2^2)/(
     2 (-1 + rho^2) s1^2 s2^2)) + 
    1/4 rho (1 + Erf[mu2/(Sqrt[2] s2)]) - ( mu1 rho)/(
     2 Sqrt[2 \[Pi]]
       s1) (1 + 
       Erf[(mu2 s1 - mu1 rho s2)/(Sqrt[2 - 2 rho^2] s1 s2)]) E^(-(
      mu1^2/(2 s1^2))) - ( mu2 rho )/(2 Sqrt[2 \[Pi]] s2)
      Erfc[(mu2 rho s1 - mu1 s2)/(Sqrt[2 - 2 rho^2] s1 s2)] E^(-(
      mu2^2/(2 s2^2))) - 
    rho T[-(mu2/s2), -(rho/Sqrt[1 - rho^2]), -(mu1/(
       Sqrt[1 - rho^2] s1))]};
MatrixForm[{l1, l2}]

It would be nice to do some additional sanity checks like finding some limiting cases. I will attempt to do it asap.

Update: Let us check the case $\mu_1=\mu_2=0$. In here we have: \begin{eqnarray} {\mathbf E}\left[\mbox{max}(X,0)\mbox{max}(Y,0)\right] &=& s_1 s_2 I_{1,1}\\ &=&s_1 s_2 \left( -\rho \cdot T\left(0,-\frac{\rho}{\sqrt{1-\rho^2}},0\right) + \frac{\rho}{4} + \frac{\sqrt{1-\rho^2}}{2\pi}\right)\\ &=& s_1 s_2 \left( \frac{\rho}{2\pi} \mbox{arctan}(\frac{\rho}{\sqrt{1-\rho^2}}) + \frac{\rho}{4} + \frac{\sqrt{1-\rho^2}}{2\pi}\right) \end{eqnarray} as it should be.

Answer 2

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])