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?