I have two normally distributed variables $X_1$ and $X_2$ with mean zero and covariance matrix $\Sigma$. I am interested in trying to calculate the value of $E[X_1^2 X_2^2]$ in terms of the entries of $\Sigma$.
I used the law of total probability to get $E[X_1^2 X_2^2] = E[X_1^2 E[X_2^2 | X_1]]$ but I'm not sure what the inner expectation reduces to. Is there another method here?
Thanks.
Edit: The variables are also multivariate normally distributed.
The expectation clearly is proportional to the product of the squared scale factors $\sigma_{11}\sigma_{22}$. The constant of proportionality is obtained by standardizing the variables, which reduces $\Sigma$ to the correlation matrix with correlation $\rho = \sigma_{12}/\sqrt{\sigma_{11}\sigma_{22}}$.
Assuming bivariate normality, then according to the analysis at https://stats.stackexchange.com/a/71303 we may change variables to
$$X_1 = X,\ X_2 = \rho X + \left(\sqrt{1-\rho^2}\right) Y$$
where $(X,Y)$ has a standard (uncorrelated) bivariate Normal distribution, and we need only compute
$$\mathbb{E}\left(X^2 (\rho X + \left(\sqrt{1-\rho^2}\right) Y)^2\right) = \mathbb{E}(\rho^2 X^4 + (1-\rho^2)X^2 Y^2 + c X^3 Y)$$
where the precise value of the constant $c$ does not matter. ($Y$ is the residual upon regressing $X_2$ against $X_1$.) Using the univariate expectations for the standard normal distribution
$$\mathbb{E}(X^4)=3,\ \mathbb{E}(X^2) = \mathbb{E}(Y^2)=1,\ \mathbb{E}Y=0$$
and noting that $X$ and $Y$ are independent yields
$$\mathbb{E}(\rho^2 X^4 + (1-\rho^2)X^2 Y^2 + c X^3 Y) = 3\rho^2 + (1-\rho^2) + 0 = 1 + 2\rho^2.$$
Multiplying this by $\sigma_{11}\sigma_{22}$ gives
$$\mathbb{E}(X_1^2 X_2^2) = \sigma_{11}\sigma_{22} + 2\sigma_{12}^2.$$
The same method applies to finding the expectation of any polynomial in $(X_1,X_2)$, because it becomes a polynomial in $(X, \rho X + \left(\sqrt{1-\rho^2}\right)Y)$ and that, when expanded, is a polynomial in the independent normally distributed variables $X$ and $Y$. From
$$\mathbb{E}(X^{2k}) = \mathbb{E}(Y^{2k}) = \frac{(2k)!}{k!2^k} = \pi^{-1/2} 2^k\Gamma\left(k+\frac{1}{2}\right)$$
for integral $k\ge 0$ (with all odd moments equal to zero by symmetry) we may derive
$$\mathbb{E}(X_1^{2p}X_2^{2q}) = (2q)!2^{-p-q}\sum_{i=0}^q \rho^{2i}(1-\rho^2)^{q-i}\frac{(2p+2i)!}{(2i)! (p+i)! (q-i)!}$$
(with all other expectations of monomials equal to zero). This is proportional to a hypergeometric function (almost by definition: the manipulations involved are not deep or instructive),
$$\frac{1}{\pi} 2^{p+q} \left(1-\rho ^2\right)^q \Gamma \left(p+\frac{1}{2}\right) \Gamma \left(q+\frac{1}{2}\right) \, _2F_1\left(p+\frac{1}{2},-q;\frac{1}{2};\frac{\rho ^2}{\rho ^2-1}\right).$$
The hypergeometric function times $\left(1-\rho ^2\right)^q$ is seen as a multiplicative correction for nonzero $\rho$.
