Q:

Variance of XY for 2 correlated normal variables

Question 5 in Allan Glut's "An Intermediate Course in Probability" asks "Suppose that X and Y have a joint normal distribution with $E(X)=E(Y)=0$ and $Var(X)=\sigma_x^2$, $Var(Y)=\sigma_y^2$, and a correlation coefficient $\rho$. Compute $E(XY)$ and $Var(XY)$. One may use the fact that $X$ and a suitable linear combination of $X$ and $Y$ are independent." 
I found $E(XY)=\rho\sigma_x\sigma_y$ easily enough but I'm stumped on $Var(XY)$. I can set $U=-\rho X+\sigma_x^2Y$ and then $U$ and $X$ are independent and $Var(U)=\sigma_x^4\sigma_y^2-\rho^2\sigma_x^2$. Then I can do:
$$Var(XY)=E(X^2Y^2)-(E(XY))^2=E\left(X^2\left(\frac{\rho X-U}{\sigma_x^2}\right)^2\right)-\rho^2\sigma_x^2\sigma_y^2$$
$$=\frac{1}{\sigma_x^4}\left(\rho^2E(X^4)-2E(X^3U)+E(X^2U^2)\right)-\rho^2\sigma_x^2\sigma_y^2$$
And since $U$ and $X$ are independent, this gives:
$$=\frac{1}{\sigma_x^4}\left(\rho^2E(X^4)-2E(X^3)E(U)+E(X^2)E(U^2)\right)-\rho^2\sigma_x^2\sigma_y^2$$
And using moment-generating functions, I found $E(X^3)=0$ and $E(X^4)=3\sigma_x^4$, which gives $\sigma_x^2\sigma_y^2(1-\rho^2)+2\rho^2$. The answer in the book is just $\sigma_x^2\sigma_y^2(1-\rho^2)$. 
Is there an easier way to do this, i.e., without moment-generating functions?  And any idea what I'm doing wrong? 

A:

If $A$ and $B$ are independent random variables, then $\operatorname{Var}(A - B) = \operatorname{Var}(A) + \operatorname{Var}(-B) = \operatorname{Var}(A) + \operatorname{Var}(B)$. This means $\operatorname{Var}(U) = \rho^2 \sigma_x^2 + \sigma_x^4\sigma_y^2$.
Also, $\operatorname{cov}(X, U) = -\rho \sigma_x^2 + \rho \sigma_x^3 \sigma_y$, so $U$ and $X$ are not independent. 
Instead, you need to chose $U = X - \frac{\sigma_x}{\sigma_y \rho} Y$. Then $X$ and $U$ are independent and you should get the desired result. 
You don't need the moment generating function of the normal distribution to calculate its moments.  The odd moments vanish by symmetry and the even moments can be calculated via partial integration (Note that $x^{2n} \exp(-x^2/2) = x^{2n - 1} \cdot x\exp(-x^2/2)$).