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.
 A: 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}.$$ 
A: 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.

