Expected value of $\max\{X,Y\}$ when $(X,Y)$ is bivariate normal 
Random vector $(X,Y)\sim N(0,0,1,1,\rho)$, that is to say, the density function of $(X,Y)$ is given by
  $$f(x,y)=\frac{1}{2\pi\sqrt{1-\rho^2}}\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2-2\rho xy+y^2)\right\}$$
  Prove that 
  $$E\max\{X,Y\}=\sqrt{\frac{1-\rho}{\pi}}.$$

I have trouble evaluating the integral:
$$\iint_{x\ge y}x\exp\left\{-\frac{1}{2(1-\rho^2)}(x^2-2\rho xy+y^2)\right\}dxdy$$
This can be turned into 
$$\int_{-\infty}^{\infty}\exp\left\{-\frac{y^2}{2}\right\}dy\int_{y}^{\infty}x\exp\left\{-\frac{(x-\rho y)^2}{2(1-\rho^2)}\right\}dx$$
However, when changing variable in the second part, I cannot get rid of the $\rho y \exp\{\cdots\}$and thus got stuck here.
Any hint or solutions are welcomed, thanks!
 A: Using the change of variable $x=z+\rho y$ in the inner integral yields the integral you are after as the sum of two terms. 
The inner integral of the first term involves the function $z\mathrm e^{-cz^2/2}$ for some positive $c$, which is easily integrated since it has a primitive proportional to $\mathrm e^{-cz^2/2}$. 
The second term is proportional to
$$
\int_\mathbb R\mathrm e^{-y^2/2}\left(\int_{(1-\rho)y}^\infty y\mathrm e^{-cz^2/2}\mathrm dz\right)\mathrm dy=\int_\mathbb R\mathrm e^{-cz^2/2}\left(\int_{-\infty}^{z/(1-\rho)}y\mathrm e^{-y^2/2}\mathrm dy\right)\mathrm dz,
$$
that is,
$$
-\int_\mathbb R\mathrm e^{-cz^2/2}\mathrm e^{-bz^2/2}\mathrm dz,
$$
for some positive $b$. This is the integral of a multiple of a gaussian density, hence has a well known value.

An easier road is to note that $2\max(X,Y)=X+Y+|X-Y|$, that $E[X]=E[Y]=0$ and that $X-Y$ is normal centered with variance $\sigma^2=2(1-\rho)$ hence
$$
2E[\max(X,Y)]=E[|X-Y|]=\sigma E[|Z|],
$$
where $Z$ is standard gaussian. A standard computation yields $E[|Z|]=\sqrt{2/\pi}$ hence
$$
E[\max(X,Y)]=\frac{\sigma}2\cdot\sqrt{\frac2\pi}=\sqrt{\frac{\sigma^2}{2\pi}}=\sqrt{\frac{1-\rho}{\pi}}.
$$
A: If you are interested further in this topic, there is a published paper by Nadarajah and Kotz that derives the pdf of max$(X,Y)$ as an Azzalini skew-Normal, and derives the expectation etc:

Nadarajah, S. and Kotz, S. (2008), "Exact Distribution of the Max/Min of Two Gaussian Random Variables", IEEE TRANSACTIONS ON VERY LARGE SCALE INTEGRATION (VLSI) SYSTEMS, VOL. 16, NO. 2, FEBRUARY 2008

Some of these results go back years.
If you are interested in the trivariate case, i.e. the pdf of max$(X,Y,Z)$, when $(X,Y,Z)$ ~ trivariate Normal, there is working paper / presentation available for download on the internet by Balakrishnan ... just google it and you will (hopefully) find it.

A Skewed Look at Bivariate and Multivariate Order Statistics
  Prof. N. Balakrishnan
  Working paper / presentation (2007)

