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*}$$
 A: Another approach: $Z=X-\theta Y$ and $Y$ are independent because they are jointly normal with covariance $0$.
$EX^{2}Y^{2}= E[(Z+\theta Y)^{2}Y^{2}]=E[Z^{2}Y^{2}+2\theta ZY^{3}+\theta^{2}Y^{4}]$ which is easy to compute because of independence.
A: Here is one approach:
Recall that if $(X,Y)$ is bivariate normal with $X\sim N(\mu_X,\sigma_X^2),Y\sim N(\mu_Y,\sigma_Y^2), \text{Corr}(X,Y)=\rho$, then
$$Y|(X=x)\sim N\left(\mu_Y+\rho{\sigma_Y \over \sigma_X}(x-\mu_X), (1-\rho^2)\sigma^2_Y\right).$$
Now note that $$E[(XY)^2]=\text{Var}(XY)+(E[XY])^2.$$
By Eve's law,
$$\begin{align}\text{Var}(XY)&=E[X^2\text{Var}(Y|X)]+\text{Var}(XE[Y|X])\\
&=E[X^2(1-\rho^2)\sigma^2_Y]+\text{Var}\left(X\left(\mu_Y+\rho{\sigma_Y \over \sigma_X}(X-\mu_X)\right)\right)\end{align},$$
which under your parameters would simplify to $(1-\theta^2)E[X^2]+\theta^2\text{Var}(X^2)=1+\theta^2$  (since $X^2\sim \chi^2_1$).
Further, $$(E[XY])^2=(\text{Cov}(X,Y)+E[X]E[Y])^2,$$
which under your parameters would simplify to $\theta^2.$
The final expression works out to $1+2\theta^2.$
A: $$\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

