Covariance of two geometric Brownian motions Assume we have two geometric Brownian motions
$$
dX_t = \mu X_t dt + \sigma X_t dW^1_t, \qquad \qquad dY_t = \mu Y_t dt + \sigma Y_t dW^2_t
$$
where the Wiener processes are correlated with $E[dW^1_t dW_t^2]  =\rho dt$. I want to show that 
$$
\text{Cov}(X_t, Y_t) = X_0Y_0e^{2\mu t}(e^{\rho\sigma^2t}-1).
$$
(This holds accoring to wikipedia )
My try:
If one takes the Ito solutions of the differentials, we have
$$
X_t = X_0e^{(\mu-\sigma^2/2)t+\sigma W^1_t}\\
Y_t = Y_0e^{(\mu-\sigma^2/2)t+\sigma W^2_t}
$$
Now caluclating $\text{Cov} = E[X_t Y_t]-E[X_t]E[Y_t]$, I start with $E[X_t Y_t]$:
$$
E[X_t Y_t] = X_0 Y_0 e^{(2\mu - \sigma^2)t} E[e^{\sigma(W_t^1 + W_t^2)}]
$$
Here is where I am not sure about how to calculate the last expectation. I can think of first deriving the joint distribution of the Wiener processes and then caluclate the integral but it seems like a huge detour. I would hope there is some more direct way, involving properties of the Brownian motion?
 A: Basically, you can simply apply Itô's formula:
$$\begin{align*} f(W_t^1, W_t^2)-f(0,0) &= \int_0^t f_x(W_s^1,W_s^2) \, dW_s^1 + \int_0^t f_y(W_s^1,W_s^2) \, dW_s^2 \\ &\quad +\frac{1}{2} \bigg( \int_0^t f_{xx} (W_s^1,W_s^2) \, \underbrace{dW_s^1 dW_s^1}_{ds} + 2 \int_0^t f_{xy}(W_s^1,W_s^2) \, \underbrace{dW_s^1 dW_s^2}_{\varrho \, ds} \\ &\quad + \int_0^t f_{yy}(W_s^1,W_s^2) \, \underbrace{dW_s^2 dW_s^2}_{ds} \bigg). \end{align*}$$
If we choose $f(x,y) := \exp(\sigma(x+y))$ and take the expectation on both sides, it follows that
$$\mathbb{E}e^{\sigma (W_t^1+W_t^2)}-1 = \sigma^2(\varrho+1) \int_0^t \mathbb{E}e^{\sigma (W_s^1+W_s^2)} \, ds.$$
Setting $\varphi(t) := \mathbb{E}e^{\sigma (W_t^1+W_t^2)}$ we find that $\varphi$ is a solution to the ordinary differential equation
$$\varphi'(t) = \sigma^2 (\varrho+1) \varphi(t) \qquad \varphi(0)=1.$$
This ODE can be solved explicitly,
$$\varphi(t) = \exp(\sigma^2(\varrho+1)t).$$
Combining the results, we conclude
$$\mathbb{E}(X_t\cdot Y_t) = X_0 Y_0 e^{2\mu t} e^{\sigma^2 \varrho t}.$$
