Multivariate normal distribution moments I would like to evaluate the following higher order moments of a multivariate normal distribution in the case of mean $0$ and in the case of mean $\mu$:
\begin{equation}
E[X_i^{2 n}] \qquad E[(X_i^2 X_{i+1}^2)^n]
\end{equation}
In the $0$ mean case I understand from the Wick Theorem that we should have $E[X_i^{2 n}]= \frac{(2 n -1)!}{2^{n-1}(n-1)!}E[X_i^{2}]^n$ but I cannot obtain the combinatorial factors of the other. In the non-central case I am quite lost.
 A: Let $Z$ be a standard normal r.v. and set $\sigma_i\equiv\sigma_{ii}$. If $\mu_i\ne 0$,
\begin{align}
\mathsf{E}X_i^{2n}&=\mathsf{E}(\sigma_iZ+\mu_i)^{2n}=\sum_{k=0}^{2n}\binom{2n}{k} \sigma_i^{k}\mu_i^{2n-k}\mathsf{E}Z^k \\
&=\sum_{k=0}^{n}\binom{2n}{2k} \sigma_i^{2k}\mu_i^{2(n-k)}(2k-1)!!
\end{align}
because $\mathsf{E}Z^{2k}=(2k-1)!!$. When $\mu_i=0$,
$$
\mathsf{E}X_i^{2n}=\mathsf{E}(\sigma_iZ)^{2n}=\sigma_i^{2n}\mathsf{E}Z^{2n}=\sigma_i^{2n}(2n-1)!!.
$$

For the expectation of cross-products let $Z_1$ and $Z_2$ be independent standard normal r.v.s. Then $(X_i,X_j)\overset{d}{=}(v_iZ_1,v_{ij}Z_1+v_j Z_2)+(\mu_i,\mu_j)$, where
$$
\begin{bmatrix}
v_i & 0 \\
v_{ij} & v_j
\end{bmatrix}=\frac{1}{\sigma_i}
\begin{bmatrix}
\sigma_i^2 & 0 \\
\sigma_{ij} & \sqrt{\sigma_i^2\sigma_j^2-\sigma_{ij}^2}
\end{bmatrix}
$$
is the Cholesky decomposition of $\operatorname{Var}([X_i, X_j]^{\top})$. Using the multinomial theorem (when $\mu_i\ne 0$, $\mu_j\ne 0$, and $\sigma_{ij}\ne 0$),
$$
\mathsf{E}[X_iX_j]^{2n}=\sum_{k_1+\cdots+k_5=2n}\binom{2n}{k_1,\ldots,k_5}\prod_{l=1}^5 \alpha_l^{k_l}\times \mathsf{E}Z_1^{k_1+2k_2+k_3}\mathsf{E}Z_2^{k_1+k_4},
$$
where
$$
\begin{align}
\alpha_1&=v_iv_j, \quad \alpha_2=v_iv_{ij}, \\
\alpha_3&=v_i\mu_j+v_{ij}\mu_i, \\
\alpha_4&=v_j\mu_i, \quad \alpha_5=\mu_i\mu_j.
\end{align}
$$
When $\mu_i=\mu_j=0$ and $\sigma_{ij}\ne 0$, 
$$
\mathsf{E}[X_iX_j]^{2n}=\sum_{k=0}^n \alpha_1^{2k}\alpha_2^{2(n-k)}(2(2n-k))!!\,(2k-1)!!.
$$
A: (Sorry I don't have enough reputation to comment, but) 
What's the significance of the index $i$ here? Are you looking at a stochastic process? Otherwise, can we simply consider the bivariate case? If so, let's denote $i$ by $1$ and $j$ by $2$. 
If $X_1$ and $X_2$ are correlated, we can represent $X_2$ as a linear combination of $X_1$ and some independent $Y$. In any case, moments of products of normal variables (and their powers) can be found here, for example. 
