# 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.

• The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral. – Jonathan Davidson Jul 3 '17 at 5:50
• @JonathanDavidson Are you referring to the second integral in the posted question? If yes, keep in mind the PDF has non-zero mean. $f(x,y)$ has the form $\text{exp}( (\frac{x - \mu_x}{\sigma_x})^2 + .. )$. A proper substitution there would be $u = \frac{x - \mu_x}{\sigma_x}$. However, that will change the integral limits resulting in many difficulties down the road. – Adel Bibi Jul 3 '17 at 6:09
• I said doable, not easy – Jonathan Davidson Jul 3 '17 at 6:10
• @JonathanDavidson I see. Unfortunately, I did go through that route. Once the integral limits change, the second integral becomes nasty and not doable. It will involve integrating over the error function multiplied by a polynomial. As far as mathematica is concerned, it doesn't seem to have a closed solution in terms of known functions. – Adel Bibi Jul 3 '17 at 6:13

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: $$$$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)$$$$

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 $$$$\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]$$$$ 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.

• "a claim has been made that the general case is somehow beyond reach" On this page? Where? – Did Jan 31 at 20:18
• @ Did In your answer in here math.stackexchange.com/questions/381161/… . – Przemo Feb 1 at 9:51
• No. Please learn to read what is written and not what you imagine. – Did Feb 5 at 20:29

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: $$$$\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]$$$$ 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])


• Please compare your "Assume for simplicity that the variables have variance one and mean zero" to the mention that "The mean is not zero" (the question). The general case is not trivially reducible to the centered case. – Did Jan 31 at 12:32
• @Did Yes it is not. Yet there is no reason why we should not look at special cases . It always gives insight to analyze such cases rather than just saying "It is not feasible". As a matter of fact even in the general case there are certain terms which are doable and other that are not. It makes sense to separate those rather than saying ..well nothing can be done anyway. – Przemo Jan 31 at 12:40
• Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution." – Did Jan 31 at 13:26
• Alright, the way I phrased my answer was suggestive of what you stated namely that the general case could be reduced to this special one which, as you rightfully pointed out, is not true. Yet I do not take back that all calculations in here are pretty straightforward and boil down to integration by parts and differentiation with respect to a parameter (in case of the generalized Owen's T function). I had the gut feeling that that function is the only thing we need to to complete calculations of such questions. So far I was right in that. – Przemo Jan 31 at 17:57