Correlated joint normal distribution: calculating a probability

Given

$$ f_{XY}(x,y) = \frac{1}{2\pi \sqrt{1-\rho^2}} \exp \left( -\frac{x^2 +y^2 - 2\rho xy}{2(1-\rho^2)} \right) $$

$Y = Z\sqrt{1-\rho^2} + \rho X$

And

$$ f_{XZ}(x,z) = \frac{1}{2\pi } \exp \left( -\frac{x^2 +z^2}{2} \right) $$

Show that $P(X>0,Y>0)= \frac{1}{4}+\frac{1}{2\pi}(\arcsin \rho) $

I'm supposed to use the fact that X and Z are independent standard normal random variables, but I don't quite understand how. Any help would be greatly appreciated.

Solution 1:

Here's a solution that only uses linear algebra and geometry.

If $\pmatrix{X\\ Y}$ is bivariate normal with mean $\pmatrix{0\\0}$ and covariance matrix $\Sigma=\pmatrix{1&\rho\\\rho&1}$, then $\pmatrix{U\\V}=\Sigma^{-1/2} \pmatrix{X\\Y}$ is bivariate normal with mean $\pmatrix{0\\0}$ and covariance matrix $\pmatrix{1&0\\ 0&1}.$ That is, $U$ and $V$ are independent, standard normal random variables.

The illustration below shows that the probability that $\pmatrix{X\\Y}$ lies in the upper quadrant (in blue), is the same as the probability that $\pmatrix{U\\V}$ lies in the wedge (in orange). Since the distribution of $\pmatrix{U\\V}$ is rotationally invariant, simple geometry gives $\mathbb{P}(X>0,Y>0)={\theta\over 2\pi}.$

With $v=\Sigma^{-1/2}\pmatrix{0\\1}$ and $w=\Sigma^{-1/2}\pmatrix{1\\0},$ we have $\cos(\theta)=\langle v,w\rangle /\|v\| \|w\|.$ But \begin{eqnarray*} \langle v,w\rangle &=& (0\ 1)\,\Sigma^{-1}\pmatrix{1\\0}=-\rho/(1-\rho^2)\\[5pt] \|v\|^2&=&(0\ 1)\,\Sigma^{-1}\pmatrix{0\\1}=1/(1-\rho^2)\\[5pt] \|w\|^2&=&(1\ 0)\,\Sigma^{-1}\pmatrix{1\\0}=1/(1-\rho^2), \end{eqnarray*} so that $\cos(\theta)=-\rho.$ Putting it all together gives $$\mathbb{P}(X>0,Y>0)={\arccos(-\rho)\over 2\pi}.$$

Solution 2:

First, for the connections among $X, Y,$ and $Z$, please see this page on how to generate two correlated normal random variables. I will leave the verification of means, variances, and correlation to you. Also, I assume $p$ in your (edited) post should be $\rho$, and you should change that.

Here is a brief simulation in R, for the case $\rho = 0.8,$ that might be helpful.

 m = 10^6;  rho = .8;  dlt = sqrt(1 - rho^2)
 x = rnorm(m);  z = rnorm(m)  # indep std norm
 y = dlt*z + rho*x
 mean(x > 0 & y > 0)
 ## 0.397643
 1/4 + 1/(2*pi)*asin(rho)
 ## 0.3975836  # agrees to 3 places
 acos(-rho)/(2*pi)  # appended later in view of @Byron's Answer
 ## 0.3975836

Here is a plot of 30,000 (out of a million) $(X, Y)$ pairs simulated.