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!
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}}. $$