We have two random variables $X$ and $Y$ which are log normally distributed, with suitable parameters, what is the expected value for $\max(X,Y)$?

Given,

$$ X=e^{\mu+\sigma Z_{1}};\quad Y=e^{\nu+\tau Z_{2}};\quad Z_{1}\sim N(0,1);Z_{2}\sim N(0,1); $$

$Z_{1}, Z_{2}$ can be assumed independent if it simplifies matters.

We need to find an expression for $$E[\text{max}(X,Y)]$$

Please note I have reached the step below, but am unsure how to proceed further.

Steps Tried

\begin{eqnarray*} E\left[\max\left(X,Y\right)\right]=\int_{0}^{\infty}xf_{Y}\left(x\right)F_{X}\left(x\right)dx+\int_{0}^{\infty}yf_{X}\left(y\right)F_{Y}\left(y\right)dy \end{eqnarray*} \begin{eqnarray*} \int_{0}^{\infty}xf_{Y}\left(x\right)F_{X}\left(x\right)dx{\displaystyle =\int_{0}^{\infty}\frac{1}{\tau}\phi\left(\frac{\ln x-\nu}{\tau}\right)}\Phi\left(\frac{\ln x-\mu}{\sigma}\right)dx \end{eqnarray*} \begin{eqnarray*} {\displaystyle =\int_{0}^{\infty}\frac{1}{\tau}\phi\left(\frac{\ln x-\nu}{\tau}\right)}\Phi\left(\frac{\ln x-\mu}{\sigma}\right)dx\quad\text{, Substitution }u=\left(\frac{\ln x-\nu}{\tau}\right)\Rightarrow du=\frac{1}{x\tau}dx \end{eqnarray*} \begin{eqnarray*} {\displaystyle =\int_{-\infty}^{\infty}e^{u\tau+\nu}\phi\left(u\right)}\Phi\left(\frac{u\tau+\nu-\mu}{\sigma}\right)du \end{eqnarray*} \begin{eqnarray*} {\displaystyle =e^{\nu}\int_{-\infty}^{\infty}e^{u\tau}\phi\left(u\right)}\Phi\left(\frac{u\tau+\nu-\mu}{\sigma}\right)du \end{eqnarray*}

Related Question

Please note, an earlier question considers the case where there is only one source of randomness. This earlier question was mis-phrased due to my limited knowledge; but still provides an interesting and instructive solution.

Expected Value of Maximum of Two Lognormal Random Variables with One Source of Randomness

Please let me know of any other suggestions …

WLOG, we assume that $\sigma>0$ and $\tau >0$. Then \begin{align} E(\max(X, Y)) &= \iint\limits_{\mu+\sigma z_1 \ge \nu + \tau z_2}e^{\mu+\sigma z_1}\frac{1}{2\pi}e^{-\frac{1}{2}(z_1^2+z_2^2)}dz_1 dz_2 \\ &\quad + \iint\limits_{\mu+\sigma z_1 \le \nu + \tau z_2}e^{\nu+\tau z_2}\frac{1}{2\pi}e^{-\frac{1}{2}(z_1^2+z_2^2)}dz_1 dz_2. \end{align} Making the substitutions \begin{align*} x = \frac{\sigma z_1 - \tau z_2}{\sqrt{\sigma^2+\tau^2}} \mbox{ and } y = \frac{\tau z_1 + \sigma z_2}{\sqrt{\sigma^2+\tau^2}}, \end{align*} we obtain that \begin{align*} & \ \iint\limits_{\mu+\sigma z_1 \ge \nu + \tau z_2}e^{\mu+\sigma z_1}\frac{1}{2\pi}e^{-\frac{1}{2}(z_1^2+z_2^2)}dz_1 dz_2\\ =&\ \int_{-\infty}^{\infty} dy \int_{\frac{\nu-\mu}{\sqrt{\sigma^2+\tau^2}}}^{\infty}e^{\mu+\sigma \frac{\sigma x + \tau y}{\sqrt{\sigma^2+\tau^2}}}\frac{1}{2\pi}e^{-\frac{1}{2}(x^2+y^2)}dxdy\\ =&\ \int_{-\infty}^{\infty} dy \int_{\frac{\nu-\mu}{\sqrt{\sigma^2+\tau^2}}}^{\infty}e^{\mu+\frac{1}{2}\sigma^2}\frac{1}{2\pi}e^{-\frac{1}{2}\left[\big(x-\frac{\sigma^2}{\sqrt{\sigma^2+\tau^2}}\big)^2+\big(y-\frac{\sigma\tau}{\sqrt{\sigma^2+\tau^2}}\big)^2\right]}dxdy\\ =&\ \int_{\frac{\nu-\mu}{\sqrt{\sigma^2+\tau^2}}}^{\infty}e^{\mu+\frac{1}{2}\sigma^2}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}\big(x-\frac{\sigma^2}{\sqrt{\sigma^2+\tau^2}}\big)^2}dx\\ =&\ \int_{-\infty}^{-\frac{\nu-\mu}{\sqrt{\sigma^2+\tau^2}}}e^{\mu+\frac{1}{2}\sigma^2}\frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}\big(x+\frac{\sigma^2}{\sqrt{\sigma^2+\tau^2}}\big)^2}dx\\ =&\ e^{\mu+\frac{1}{2}\sigma^2} N\left(\frac{-\nu+\mu+\sigma^2}{\sqrt{\sigma^2+\tau^2}} \right), \end{align*} where $N$ is the cumulative distribution function of a standard normal random variable. Similarly, \begin{align*} \iint\limits_{\nu + \tau z_2 \ge \mu+\sigma z_1}e^{\nu+\tau z_2}\frac{1}{2\pi}e^{-\frac{1}{2}(z_1^2+z_2^2)}dz_1 dz_2 &= e^{\nu+\frac{1}{2}\tau^2} N\left(\frac{-\mu+\nu+\tau^2}{\sqrt{\sigma^2+\tau^2}} \right). \end{align*} That is, \begin{align*} E(\max(X, Y)) &= e^{\nu+\frac{1}{2}\tau^2} N\left(\frac{-\mu+\nu+\tau^2}{\sqrt{\sigma^2+\tau^2}} \right) + e^{\mu+\frac{1}{2}\sigma^2} N\left(\frac{-\nu+\mu+\sigma^2}{\sqrt{\sigma^2+\tau^2}} \right). \end{align*}

Alternative solution.

Note that \begin{align*} \max(X, Y) &= Y + \max(X-Y, \ 0)\\ &= Y + Y\max\Big(\frac{X}{Y}-1, \ 0\Big). \end{align*} Let $P$ be the given probability measure, and $Q$ be a probability measure that is defined by the Randon-Nykodim derivative \begin{align*} \frac{dQ}{dP} = e^{-\frac{1}{2}\tau^2 + \tau Z_2}. \end{align*} Then \begin{align*} \widehat{Z}_1 = Z_1, \mbox{ and } \widehat{Z}_2 = Z_2 - \tau \end{align*} are two standard independent normal random variables under the probability measure $Q$. Moreover, \begin{align*} Z \equiv \frac{X}{Y} &= e^{\mu-\nu - \tau^2 + \sigma \widehat{Z}_1 - \tau \widehat{Z}_2}\\ &= e^{\mu-\nu - \tau^2 + \sqrt{\sigma^2+\tau^2} \xi}, \end{align*} where \begin{align*} \xi = \frac{\sigma \widehat{Z}_1 - \tau \widehat{Z}_2}{\sqrt{\sigma^2+\tau^2}} \end{align*} is a standard normal random variable. Let \begin{align*} \lambda = \frac{-\mu+ \nu + \tau^2}{\sqrt{\sigma^2+\tau^2}}, \end{align*} and $E_Q$ be the expectation operator under the measure $Q$. Then, \begin{align*} E\big(\max(X, Y) \big) &= E(Y)+E\big(Y\max(Z-1, \ 0) \big)\\ &=E(Y)+e^{\nu+\frac{1}{2}\tau^2} E_Q\big(\max(Z-1, \ 0) \big). \end{align*} Furthermore, \begin{align*} E_Q\big(\max(Z-1, \ 0) \big) &= E_Q\left(Z\,\mathbb{I}_{Z \ge 1}\right) - E_Q\left(\mathbb{I}_{Z \ge 1}\right)\\ &=\int_{\lambda}^{\infty} \frac{1}{\sqrt{2\pi}}e^{\mu-\nu -\tau^2 + \sqrt{\sigma^2+\tau^2} z -\frac{1}{2}z^2} dz - N(-\lambda)\\ &=e^{\mu-\nu + \frac{1}{2} (\sigma^2-\tau^2)} \int_{\lambda}^{\infty} \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}\big(z-\sqrt{\sigma^2+\tau^2}\big)^2} dz - N(-\lambda)\\ &=e^{\mu-\nu + \frac{1}{2} (\sigma^2-\tau^2)}N\Big(\sqrt{\sigma^2+\tau^2} -\lambda\Big) - N(-\lambda). \end{align*} That is, \begin{align*} E\big(\max(X, Y) \big) &= e^{\nu + \frac{1}{2}\tau^2} \bigg[1 + e^{\mu-\nu + \frac{1}{2} (\sigma^2-\tau^2)}N\Big(\sqrt{\sigma^2+\tau^2} -\lambda\Big) - N(-\lambda) \bigg]\\ &=e^{\nu + \frac{1}{2}\tau^2} N(\lambda) + e^{\mu + \frac{1}{2} \sigma^2}N\Big(\sqrt{\sigma^2+\tau^2} -\lambda\Big)\\ &= e^{\nu + \frac{1}{2}\tau^2}N\left(\frac{-\mu+ \nu + \tau^2}{\sqrt{\sigma^2+\tau^2}}\right) + e^{\mu + \frac{1}{2} \sigma^2}N\left(\frac{-\nu+\mu+\sigma^2}{\sqrt{\sigma^2+\tau^2}} \right). \end{align*}

Comments. You can make this question more general by assuming that \begin{align*} X=e^{\mu+\sigma Z_1} \quad \mbox{and} \quad Y = e^{\nu+\tau\big(\rho Z_1 + \sqrt{1-\rho^2} Z_2\big)}, \end{align*} where $Z_1$ and $Z_2$ are two independent standard normal random variables, and $|\rho|<1$. In this case, the measure $Q$ can be defined by \begin{align*} \frac{dQ}{dP} = e^{-\frac{1}{2}\tau^2 + \tau \big(\rho Z_1 + \sqrt{1-\rho^2} Z_2\big)}. \end{align*} Moreover, \begin{align*} \widehat{Z}_1 = Z_1 - \rho\tau, \mbox{ and } \widehat{Z}_2 = Z_2 - \sqrt{1-\rho^2}\tau \end{align*} are two standard independent normal random variables under the probability measure $Q$. See also this question at https://quant.stackexchange.com/questions/21361/how-to-use-a-change-of-numeraire-to-price-this-option/21375#21375.