Help with calculating a complicated expectation

I'm trying to use Mathematica to help me simplify an expression involving expectations of products of powers for correlated, zero-mean, unit-variance normal variates. Specifically, for $(Z_1,Z_2,Z_3,Z_4)\sim N(0,\Sigma)$ where $\Sigma=(\rho_{ij})$ is the correlation matrix, I'd like to express $\mathbb{E}[(Z_1+Z_2+Z_3+Z_4)^4]$ in terms of the $\rho_{ij}$'s.

The results for the expectations of products of the $Z$'s up to 4th degree are known, for instance $\mathbb{E}(Z_j^4)=3$, $\mathbb{E}(Z_j^2Z_k^2)=1+2\rho_{jk}^2$ and $\mathbb{E}(Z_j^2Z_kZ_l)=2\rho_{jk}\rho_{jl}+\rho_{kl}$, etc. My goal is to substitute these expressions into the expansion of $\mathbb{E}[(Z_1+Z_2+Z_3+Z_4)^4]$ and simplify.

Apparently this is a very tedious task by hand. Is it even possible to do this in Mathematica? I'm OK with having to write some code for this, since my ultimate aim is to generalize this to $\mathbb{E}[(Z_1+\cdots+Z_d)^4]$, where $d>4$.

Any help is appreciated.

• Does Expectation[(Z1 + Z2 + Z3 + Z4)^4, {Z1, Z2, Z3, Z4} \[Distributed] MultinormalDistribution[{0, 0, 0, 0}, Table[p[i, j], {i, 1, 4}, {j, 1, 4}]]] not do what you want? – Patrick Stevens Sep 9 '15 at 22:39
• Thanks! I didn't realize it does that automatically without having to substitute expressions... – David L Sep 9 '15 at 22:45
• Would you like me to post it as an answer? (Or you could answer yourself, if you like.) – Patrick Stevens Sep 9 '15 at 22:45

Expectation[