Solution (source code)

= Solution

Write $s=\sqrt{1-\rho^2}>0$. The inverse change of variables is $y=\rho x+sz$, whose <Jacobian determinant> has absolute value $s$. Completing the square gives
$$
\frac{x^2-2\rho xy+y^2}{1-\rho^2}=x^2+z^2.
$$
The <change of variables formula> therefore transforms the joint <probability density function> to
$$
g(x,z)=s f(x,\rho x+sz)
=\frac1{2\pi}e^{-(x^2+z^2)/2}
=\left(\frac{e^{-x^2/2}}{\sqrt{2\pi}}\right)
 \left(\frac{e^{-z^2/2}}{\sqrt{2\pi}}\right).
$$
It factors into two normalized <standard normal distribution> densities. Thus \b[$X$ and $Z$ are <independent> $N(0,1)$ variables]; factorization proves the requested <independence>, rather than simply zero <covariance>.

The <standard Gaussian random vector> $(X,Z)$ has a rotation-invariant density. In <polar coordinates>, $X=r\cos\theta$, $Z=r\sin\theta$, its density with the area factor is $(2\pi)^{-1}re^{-r^2/2}$, so the angular coordinate is uniform on an interval of length $2\pi$.

Put $\alpha=\arcsin\rho\in(-\pi/2,\pi/2)$, so $s=\cos\alpha$ and
$$
Y=\rho X+sZ=r\sin(\theta+\alpha).
$$
In the right half-plane, choose $-\pi/2<\theta<\pi/2$. The additional condition $Y>0$ then requires $-\alpha<\theta<\pi/2$. This wedge has angle $\pi/2+\alpha$; its boundary has <probability> zero. Consequently
$$
\boxed{\mathbb P(X>0,Y>0)=\frac14+\frac{\arcsin\rho}{2\pi}.}
$$
At $\rho=0$ this is $1/4$, as required by <independence>. The limits as $\rho\to1$ and $\rho\to-1$ are $1/2$ and zero, respectively. The angular calculation is also the geometric <basis> of the <Gaussian sign-correlation identity>.