How to find a closed-form expression for $E[\max(X,Y)]$? Let $X$ and $Y$ be two random variables, with
$X\sim N(a,b)$ and $Y\sim N(c,d)$. Furthermore, $X$ and $Y$ are correlated with a correlation coefficient equal to $p$.
How do I find a closed-form expression for $E[\max(X,Y)]$?
Whichever way I look at it, I am not able to get rid of the randomness in the solution.
My solution.
$E[\max(X,Y)] = XP(X>Y) + YP(Y>X)$
But if I do this I will get an answer in terms of $X$ and $Y$ which is still a random quantity.
 A: Take the joint density
$$\exp\left(-z/2(1-p^2)/(2\pi bd \sqrt{(1-p^2})\right)$$ where
$z=(x-a)^2/b^2 -2p(x-a)(y-b)/(bd) + (y-c)^2/d^2$.
Multiply it by $x$ over the region where $x>y$ and by $y$ over the region where $y>x$ and compute the double integral over those regions. This will not have a closed form but can be integrated numerically.
A: First, technically speaking, given only $X \sim N(\mu_1, \sigma_1^2), Y \sim N(\mu_2, \sigma_2^2)$ and $\operatorname{Corr}(X, Y) = \rho$ does not imply the joint distribution of $(X, Y)$ is bivariate normal $N_2(\mu_1, \mu_2, \sigma_1^2, \sigma_2^2, \rho)$.  Check this answer for many counterexamples.  Therefore, the accepted answer actually imposed additional conditions that were not originally posted in the question.
Secondly, if we do have $(X, Y) \sim N_2(\mu_1, \mu_2, \sigma_1^2, \sigma_2^2, \rho)$, then $E[\max(X, Y)]$ indeed admits a closed-form expression, which is deduced below.
Write $E[\max(X, Y)] = \frac{1}{2}E[(X + Y) + |X - Y|]$.  Under the joint normality condition, we have
\begin{align}
& X + Y \sim N(\mu_1 + \mu_2, \sigma_1^2 + \sigma_2^2 + 2\rho\sigma_1\sigma_2), \\
& X - Y \sim N(\mu_1 - \mu_2, \sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2). 
\end{align}
It thus suffices to compute $E[Z]$ and $E[|Z|]$ for a univariate r.v. $Z \sim N(\mu, \sigma^2)$. Clearly, $E[Z] = \mu$, while $E[|Z|] = \sigma E[|Z'|]$ with $Z' \sim N(\mu\sigma^{-1}, 1)$, and
\begin{align}
& E[|Z'|] = \int_{-\infty}^\infty |x|\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}(x - \mu\sigma^{-1})^2}dx \\
=& \frac{1}{\sqrt{2\pi}}\int_0^\infty xe^{-\frac{1}{2}(x - \mu\sigma^{-1})^2}dx - 
\frac{1}{\sqrt{2\pi}}\int_{-\infty}^0 xe^{-\frac{1}{2}(x - \mu\sigma^{-1})^2}dx \\
=& \left[\frac{1}{\sqrt{2\pi}}e^{-\mu^2\sigma^{-2}/2} + \mu\sigma^{-1}\Phi(\mu\sigma^{-1})\right] - 
\left[-\frac{1}{\sqrt{2\pi}}e^{-\mu^2\sigma^{-2}/2} + \mu\sigma^{-1}\Phi(-\mu\sigma^{-1})\right] \\
=& \frac{2}{\sqrt{2\pi}}e^{-\mu^2\sigma^{-2}/2} + \mu\sigma^{-1}(2\Phi(\mu\sigma^{-1}) - 1). 
\end{align}
Denote $\sqrt{\sigma_1^2 + \sigma_2^2 - 2\rho\sigma_1\sigma_2}$ by $\sigma_0$, then
\begin{align*}
E[|X - Y|] = \sigma_0
\left(\frac{2}{\sqrt{2\pi}}e^{-(\mu_1 - \mu_2)^2\sigma_0^{-2}/2} + (\mu_1 - \mu_2)\sigma_0^{-1}(2\Phi((\mu_1 - \mu_2)\sigma_0^{-1}) - 1)\right)
\end{align*}
and
\begin{align}
 & E[\max(X, Y)] \\
=& \frac{1}{2}E[X + Y] + \frac{1}{2}E[|X - Y|] \\
=& \frac{\mu_1 + \mu_2}{2} + \frac{\sigma_0}{\sqrt{2\pi}}e^{-(\mu_1 - \mu_2)^2\sigma_0^{-2}/2} + \frac{1}{2}(\mu_1 - \mu_2)(2\Phi((\mu_1 - \mu_2)\sigma_0^{-1}) - 1) \\
=& \mu_2 + \frac{\sigma_0}{\sqrt{2\pi}}e^{-\frac{(\mu_1 - \mu_2)^2}{2\sigma_0^2}} + (\mu_1 - \mu_2)\Phi\left(\frac{\mu_1 - \mu_2}{\sigma_0}\right).
\end{align}
