Integral of bivariate normal distribution function with respect to itself Define $F: \mathbb{R}^2 \rightarrow \mathbb{R}$ by 
\begin{align*}
F(x,y)=\int_{-\infty}^{x} \int_{-\infty}^y \frac{1}{2\pi \sqrt{1-\rho^2}}  
exp\left(\frac{-u^2-v^2+2\rho uv}{2(1-\rho^2)}\right) dudv,
\end{align*}
where $\rho$ is a constant with $-1\leq \rho \leq 1$.
Could you help me to prove that $\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} F(x,y) dF(x,y)=F(0,0)$.
I'm trying to solve this by substitution $dF(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) dxdy$, separating $x$ and $y$ variable and doing integration by part but couldn't work.
 A: Note that the integral equals
$$
   \int_{-\infty}^\infty \int_{-\infty}^\infty F(x,y)\mathrm{d}F(x,y) = \int_{-\infty}^\infty \int_{-\infty}^\infty \Pr\left(X_1 \leq x, Y_1 \leq y\right)\mathrm{d}F(x,y) = \Pr\left(X_1 \leq X_2, Y_1\leq Y_2\right)
$$
where $(X_1,Y_1)$ and $(X_2,Y_2)$ are iid random vectors from binormal distribution with zero means, unit variances and the same correlation coefficient $\rho$. Using iid standard normal variables $\{Z_k\}_{k=1}^6$, $(X_1,Y_1) \stackrel{d}{=} \left(Z_1, \rho Z_1 + \sqrt{1-\rho^2} Z_2 \right)$ and $(X_2,Y_2) \stackrel{d}{=} (Z_3, Z_3 \rho + \sqrt{1-\rho^2} Z_4)$. 
Hence
$$
 \Pr\left(X_1 <X_2, Y_1<Y_2\right) = \Pr\left( Z_1 < Z_3, Z_1 \rho + \sqrt{1-\rho^2} Z_2 < Z_3 \rho + \sqrt{1-\rho^2} Z_4 \right) = \Pr\left((Z_1-Z_3) < 0, (Z_1-Z_3)\rho + \sqrt{1-\rho^2} (Z_2-Z_4) < 0\right)
$$
Since $Z_1 - Z_2 \stackrel{d}{=} \sqrt{2} Z_5$ and $Z_3 - Z_4 \stackrel{d}{=} \sqrt{2} Z_6$ we have
$$ \begin{eqnarray}
   \Pr\left(X_1 \leq X_2, Y_1\leq  Y_2\right) &=& \Pr\left(Z_5 \leq  0, Z_5 \rho + Z_6 \sqrt{1-\rho^2 } \leq  0 \right) \\ &=& \Pr\left(X_3 \leq  0, Y_3 \leq  0\right) = F(0,0)
\end{eqnarray}$$
where we used that $(Z_5, Z_5 \rho + Z_6 \sqrt{1-\rho^2} )$ is equal in distribution to $(X_1, Y_1)$.
