Mixing
Compute $mathbb{E}[text{max}(x,0) text{max}(y,0)]$ where $(x,y)$ is jointly gaussian with given covariance...
$begingroup$
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_ysqrt{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{2rho mu_xmu_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.
probability-theory normal-distribution
edited Jan 31 at 12:37
Did
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
$endgroup$
add a comment |
$begingroup$
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_ysqrt{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{2rho mu_xmu_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.
probability-theory normal-distribution
edited Jan 31 at 12:37
Did
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
$endgroup$
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
1
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13
add a comment |
$begingroup$
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_ysqrt{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{2rho mu_xmu_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.
probability-theory normal-distribution
edited Jan 31 at 12:37
Did
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
$endgroup$
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_ysqrt{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{2rho mu_xmu_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.
probability-theory normal-distribution
probability-theory normal-distribution
edited Jan 31 at 12:37
Did
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
edited Jan 31 at 12:37
Did
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
edited Jan 31 at 12:37
Did
249k23227466
edited Jan 31 at 12:37
Did
249k23227466
edited Jan 31 at 12:37
Did
249k23227466
249k23227466
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
asked Jul 3 '17 at 5:16
Adel BibiAdel Bibi
208113
208113
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
1
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13
add a comment |
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
1
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
1
1
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13
add a comment |
2 Answers
2
active
oldest
votes
$begingroup$
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]=
intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty
left(mu_1+s_1 xiright)left(mu_2+s_2 etaright) cdot rho(xi,eta) dxi deta
end{eqnarray}
where
begin{equation}
rho(xi,eta) := frac{1}{2pi sqrt{1-rho^2}} expleft[
-frac{1}{2} frac{1}{1-rho^2}
left( xi^2+eta^2 - 2 rho xi etaright)
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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty 1 cdot rho(xi,eta) dxi deta=\
&&Tleft(-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)+1right)\
&&I_{0,1}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty eta cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi eta cdot rho(xi,eta) dxi deta=\
&&-rho cdot Tleft(-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)+1right)}{2 sqrt{2 pi } s_1}+frac{1}{4} rho
left(text{erf}left(frac{mu_2}{sqrt{2} s_2}right)+1right)-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-1right) 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 Tleft(0,-frac{rho}{sqrt{1-rho^2}},0right) + frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)\
&=& s_1 s_2 left(
frac{rho}{2pi} mbox{arctan}(frac{rho}{sqrt{1-rho^2}}) +
frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)
end{eqnarray}
as it should be.
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
add a comment |
$begingroup$
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}{2pi sqrt{1-rho^2}} expleft[ -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):=intlimits_0^infty x rho(x,y) dx &=& frac{1}{2pi sqrt{1-rho^2}} intlimits_0^infty expleft[-frac{1}{2} frac{1}{(1-rho^2)} ((x-rho y)^2 + y^2(1-rho^2))right]\
&=&frac{1}{2pi sqrt{1-rho^2}}e^{-frac{1}{2} y^2} intlimits_{-rho y}^infty (x+ rho y) e^{-frac{1}{2} frac{1}{1-rho^2} x^2} dx \
&=& frac{1}{2pi} 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 $yin(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}
intlimits_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 ^2right)^{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])
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
add a comment |
Your Answer
StackExchange.ifUsing("editor", function () {
return StackExchange.using("mathjaxEditing", function () {
StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
});
});
}, "mathjax-editing");
StackExchange.ready(function() {
var channelOptions = {
tags: "".split(" "),
id: "69"
};
initTagRenderer("".split(" "), "".split(" "), channelOptions);
StackExchange.using("externalEditor", function() {
// Have to fire editor after snippets, if snippets enabled
if (StackExchange.settings.snippets.snippetsEnabled) {
StackExchange.using("snippets", function() {
createEditor();
});
}
else {
createEditor();
}
});
function createEditor() {
StackExchange.prepareEditor({
heartbeatType: 'answer',
autoActivateHeartbeat: false,
convertImagesToLinks: true,
noModals: true,
showLowRepImageUploadWarning: true,
reputationToPostImages: 10,
bindNavPrevention: true,
postfix: "",
imageUploader: {
brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
allowUrls: true
},
noCode: true, onDemand: true,
discardSelector: ".discard-answer"
,immediatelyShowMarkdownHelp:true
});
}
});
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2344676%2fcompute-mathbbe-textmaxx-0-textmaxy-0-where-x-y-is-jointly-g%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
2 Answers
2
active
oldest
votes
2 Answers
2
active
oldest
votes
active
oldest
votes
active
oldest
votes
$begingroup$
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]=
intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty
left(mu_1+s_1 xiright)left(mu_2+s_2 etaright) cdot rho(xi,eta) dxi deta
end{eqnarray}
where
begin{equation}
rho(xi,eta) := frac{1}{2pi sqrt{1-rho^2}} expleft[
-frac{1}{2} frac{1}{1-rho^2}
left( xi^2+eta^2 - 2 rho xi etaright)
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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty 1 cdot rho(xi,eta) dxi deta=\
&&Tleft(-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)+1right)\
&&I_{0,1}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty eta cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi eta cdot rho(xi,eta) dxi deta=\
&&-rho cdot Tleft(-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)+1right)}{2 sqrt{2 pi } s_1}+frac{1}{4} rho
left(text{erf}left(frac{mu_2}{sqrt{2} s_2}right)+1right)-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-1right) 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 Tleft(0,-frac{rho}{sqrt{1-rho^2}},0right) + frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)\
&=& s_1 s_2 left(
frac{rho}{2pi} mbox{arctan}(frac{rho}{sqrt{1-rho^2}}) +
frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)
end{eqnarray}
as it should be.
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
add a comment |
$begingroup$
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]=
intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty
left(mu_1+s_1 xiright)left(mu_2+s_2 etaright) cdot rho(xi,eta) dxi deta
end{eqnarray}
where
begin{equation}
rho(xi,eta) := frac{1}{2pi sqrt{1-rho^2}} expleft[
-frac{1}{2} frac{1}{1-rho^2}
left( xi^2+eta^2 - 2 rho xi etaright)
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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty 1 cdot rho(xi,eta) dxi deta=\
&&Tleft(-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)+1right)\
&&I_{0,1}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty eta cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi eta cdot rho(xi,eta) dxi deta=\
&&-rho cdot Tleft(-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)+1right)}{2 sqrt{2 pi } s_1}+frac{1}{4} rho
left(text{erf}left(frac{mu_2}{sqrt{2} s_2}right)+1right)-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-1right) 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 Tleft(0,-frac{rho}{sqrt{1-rho^2}},0right) + frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)\
&=& s_1 s_2 left(
frac{rho}{2pi} mbox{arctan}(frac{rho}{sqrt{1-rho^2}}) +
frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)
end{eqnarray}
as it should be.
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
add a comment |
$begingroup$
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]=
intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty
left(mu_1+s_1 xiright)left(mu_2+s_2 etaright) cdot rho(xi,eta) dxi deta
end{eqnarray}
where
begin{equation}
rho(xi,eta) := frac{1}{2pi sqrt{1-rho^2}} expleft[
-frac{1}{2} frac{1}{1-rho^2}
left( xi^2+eta^2 - 2 rho xi etaright)
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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty 1 cdot rho(xi,eta) dxi deta=\
&&Tleft(-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)+1right)\
&&I_{0,1}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty eta cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi eta cdot rho(xi,eta) dxi deta=\
&&-rho cdot Tleft(-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)+1right)}{2 sqrt{2 pi } s_1}+frac{1}{4} rho
left(text{erf}left(frac{mu_2}{sqrt{2} s_2}right)+1right)-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-1right) 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 Tleft(0,-frac{rho}{sqrt{1-rho^2}},0right) + frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)\
&=& s_1 s_2 left(
frac{rho}{2pi} mbox{arctan}(frac{rho}{sqrt{1-rho^2}}) +
frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)
end{eqnarray}
as it should be.
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
$endgroup$
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]=
intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty
left(mu_1+s_1 xiright)left(mu_2+s_2 etaright) cdot rho(xi,eta) dxi deta
end{eqnarray}
where
begin{equation}
rho(xi,eta) := frac{1}{2pi sqrt{1-rho^2}} expleft[
-frac{1}{2} frac{1}{1-rho^2}
left( xi^2+eta^2 - 2 rho xi etaright)
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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty 1 cdot rho(xi,eta) dxi deta=\
&&Tleft(-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)+1right)\
&&I_{0,1}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty eta cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi cdot rho(xi,eta) dxi deta=\
&&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)+1right)+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}:=intlimits_{-frac{mu_1}{s_1}}^infty
intlimits_{-frac{mu_2}{s_2}}^infty xi eta cdot rho(xi,eta) dxi deta=\
&&-rho cdot Tleft(-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)+1right)}{2 sqrt{2 pi } s_1}+frac{1}{4} rho
left(text{erf}left(frac{mu_2}{sqrt{2} s_2}right)+1right)-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-1right) 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 Tleft(0,-frac{rho}{sqrt{1-rho^2}},0right) + frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)\
&=& s_1 s_2 left(
frac{rho}{2pi} mbox{arctan}(frac{rho}{sqrt{1-rho^2}}) +
frac{rho}{4} + frac{sqrt{1-rho^2}}{2pi}right)
end{eqnarray}
as it should be.
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
edited Jan 31 at 18:15
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
answered Jan 31 at 17:52
PrzemoPrzemo
4,65811032
4,65811032
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
add a comment |
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
"a claim has been made that the general case is somehow beyond reach" On this page? Where?
$endgroup$
– Did
Jan 31 at 20:18
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
@ Did In your answer in here math.stackexchange.com/questions/381161/… .
$endgroup$
– Przemo
Feb 1 at 9:51
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
$begingroup$
No. Please learn to read what is written and not what you imagine.
$endgroup$
– Did
Feb 5 at 20:29
add a comment |
$begingroup$
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}{2pi sqrt{1-rho^2}} expleft[ -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):=intlimits_0^infty x rho(x,y) dx &=& frac{1}{2pi sqrt{1-rho^2}} intlimits_0^infty expleft[-frac{1}{2} frac{1}{(1-rho^2)} ((x-rho y)^2 + y^2(1-rho^2))right]\
&=&frac{1}{2pi sqrt{1-rho^2}}e^{-frac{1}{2} y^2} intlimits_{-rho y}^infty (x+ rho y) e^{-frac{1}{2} frac{1}{1-rho^2} x^2} dx \
&=& frac{1}{2pi} 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 $yin(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}
intlimits_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 ^2right)^{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])
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
add a comment |
$begingroup$
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}{2pi sqrt{1-rho^2}} expleft[ -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):=intlimits_0^infty x rho(x,y) dx &=& frac{1}{2pi sqrt{1-rho^2}} intlimits_0^infty expleft[-frac{1}{2} frac{1}{(1-rho^2)} ((x-rho y)^2 + y^2(1-rho^2))right]\
&=&frac{1}{2pi sqrt{1-rho^2}}e^{-frac{1}{2} y^2} intlimits_{-rho y}^infty (x+ rho y) e^{-frac{1}{2} frac{1}{1-rho^2} x^2} dx \
&=& frac{1}{2pi} 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 $yin(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}
intlimits_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 ^2right)^{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])
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
$endgroup$
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
add a comment |
$begingroup$
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}{2pi sqrt{1-rho^2}} expleft[ -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):=intlimits_0^infty x rho(x,y) dx &=& frac{1}{2pi sqrt{1-rho^2}} intlimits_0^infty expleft[-frac{1}{2} frac{1}{(1-rho^2)} ((x-rho y)^2 + y^2(1-rho^2))right]\
&=&frac{1}{2pi sqrt{1-rho^2}}e^{-frac{1}{2} y^2} intlimits_{-rho y}^infty (x+ rho y) e^{-frac{1}{2} frac{1}{1-rho^2} x^2} dx \
&=& frac{1}{2pi} 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 $yin(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}
intlimits_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 ^2right)^{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])
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
$endgroup$
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}{2pi sqrt{1-rho^2}} expleft[ -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):=intlimits_0^infty x rho(x,y) dx &=& frac{1}{2pi sqrt{1-rho^2}} intlimits_0^infty expleft[-frac{1}{2} frac{1}{(1-rho^2)} ((x-rho y)^2 + y^2(1-rho^2))right]\
&=&frac{1}{2pi sqrt{1-rho^2}}e^{-frac{1}{2} y^2} intlimits_{-rho y}^infty (x+ rho y) e^{-frac{1}{2} frac{1}{1-rho^2} x^2} dx \
&=& frac{1}{2pi} 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 $yin(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}
intlimits_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 ^2right)^{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])
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
answered Jan 31 at 12:05
PrzemoPrzemo
4,65811032
4,65811032
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
add a comment |
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
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.
$endgroup$
– Did
Jan 31 at 12:32
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
@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.
$endgroup$
– Przemo
Jan 31 at 12:40
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
Quote from main: "Price's theorem can nicely handle the zero mean case with an elegant solution."
$endgroup$
– Did
Jan 31 at 13:26
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
$begingroup$
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.
$endgroup$
– Przemo
Jan 31 at 17:57
add a comment |
Thanks for contributing an answer to Mathematics Stack Exchange!
- Please be sure to answer the question. Provide details and share your research!
But avoid …
- Asking for help, clarification, or responding to other answers.
- Making statements based on opinion; back them up with references or personal experience.
Use MathJax to format equations. MathJax reference.
To learn more, see our tips on writing great answers.
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
StackExchange.ready(
function () {
StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmath.stackexchange.com%2fquestions%2f2344676%2fcompute-mathbbe-textmaxx-0-textmaxy-0-where-x-y-is-jointly-g%23new-answer', 'question_page');
}
);
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Sign up or log in
StackExchange.ready(function () {
StackExchange.helpers.onClickDraftSave('#login-link');
});
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
Required, but never shown
$begingroup$
The integral seems doable as an iterated integral and using u substitution $u = x^2$ and $u = y^2$ for each iterated integral.
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 5:50
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:09
$begingroup$
I said doable, not easy
$endgroup$
– Jonathan Davidson
Jul 3 '17 at 6:10
1
$begingroup$
@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.
$endgroup$
– Adel Bibi
Jul 3 '17 at 6:13