<s>[INST] Show the value of $\mathbb{E}(X^2Y^2)$ Suppose that $(X,Y)$ is distributed ranodm 2-vectors having the normal distribution with $\mathbb{E}X=\mathbb{E}Y=0,\mathbb{Var}(X)=\mathbb{Var}(Y)=1,$ and $\text{Cov}(X,Y)=\theta\in(-1,1).$ Show the value of  $\mathbb{E}(X^2Y^2).$

The joint pdf of $(X,Y)$ is $$p(x,y)=\frac{1}{2\pi\sqrt{1-\theta^2}}\exp\left\{-\frac{1}{2(1-\theta^2)} \left [  x^2-2\theta xy+y^2\right ]\right\}.$$
Then $$\mathbb{E}(X^2Y^2)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}x^2y^2\cdot p(x,y)\mathrm{d}x\mathrm{d}y.$$
Let $$\begin{cases}
 u=\frac{x-y}{\sqrt{2}} \\
 v=\frac{x+y}{\sqrt{2}} 
\end{cases}\Rightarrow \mathbb{E}(X^2Y^2)=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot p\left (\frac{u+v}{\sqrt{2}},\frac{u-v}{\sqrt{2}}\right )\mathrm{d}u\mathrm{d}v.$$
But I don't know how to deal with this integral
$$\begin{align*}
 &\Delta:=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot \exp\left \{-\frac{1}{2(1-\theta^2)}\left [  (1-\theta)u^2+(1+\theta)v^2\right ]\right \}\mathrm{d}u\mathrm{d}v\\
 &\quad=4\int_{0}^{\infty}\int_{0}^{\infty}\frac{(u^2-v^2)^2}{2}\cdot \exp\left \{-\frac{1}{2(1-\theta^2)}\left [  (1-\theta)u^2+(1+\theta)v^2\right ]\right \}\mathrm{d}u\mathrm{d}v,\quad \theta\in(-1,1).
\end{align*}$$
 [/INST] $$\begin{align*}
\Bbb E\left[X^2Y^2\right] &= \frac1{2\pi \sqrt{1-\theta^2}} \int\limits_0^\infty \int\limits_0^\infty \left(u^2-v^2\right)^2 e^{\Large-\frac{(1-\theta) u^2+(1+\theta)v^2}{2(1-\theta^2)}} \, du \, dv \\
\tag{1} &= \frac4\pi \int\limits_0^\infty \int\limits_0^\infty \left((1+\theta)s^2-(1-\theta)t^2\right)^2 e^{-(s^2+t^2)} \, ds \, dt \\[1ex]
\tag{2} &= \frac4\pi \int\limits_0^{\pi/2} \int\limits_0^\infty \bigg((1+\theta)^2 r^4\cos^4(\phi) - 2(1-\theta^2)r^4\cos^2(\phi)\sin^2(\phi) \\
&\qquad \qquad \qquad \qquad + (1-\theta)^2 r^4 \sin^4(\phi)\bigg) re^{-r^2} \, dr \, d\phi \\[1ex]
\tag{3} &= \frac4\pi \int\limits_0^{\pi/2} \bigg((1+\theta)^2 \cos^4(\phi) - 2(1-\theta^2) \cos^2(\phi) \sin^2(\phi) \\
&\qquad \qquad \qquad \qquad + (1-\theta)^2 \sin^4(\phi)\bigg) \, d\phi \\[1ex]
\tag{4} &= \frac4\pi \left(\frac{3\pi}{16} (1+\theta)^2 - \frac\pi8(1-\theta^2) + \frac{3\pi}{16} (1-\theta)^2\right) \\[1ex]
&= 1+2\theta^2
\end{align*}$$


*

*$(1)$ : substitute $(u,v)=\left(\sqrt{2(1+\theta)}\,s,\sqrt{2(1-\theta)}\,t\right)$

*$(2)$ : substitute $(s,t)=(r\cos(\phi),r\sin(\phi))$

*$(3)$ : integrate by parts;
$$\int\limits_0^\infty r^5 e^{-r^2} \, dr = 1$$

*$(4)$ : reduce powers via $\cos^2(\alpha)=\dfrac{1+\cos(2\alpha)}2$ and $\sin^2(\alpha)=\dfrac{1-\cos(2\alpha)}2$ and integrate

</s>