I have a scenario, where I know that $\boldsymbol{z} \sim \mathcal{N}( \boldsymbol{\mu}, \boldsymbol{\Sigma} ) $ and I need to compute $ \mathbb{E}( z_i^2 z_j ) $, where $z_i$ is the i-th component of vector $\boldsymbol{z}$.
Is there any way of computing this, assuming that $\boldsymbol{z}$ is a random vector of dimension n>2, say e.g. n=5?
If there is not a general way of computing this, is there anything for the special case of $\boldsymbol{\mu}=\boldsymbol{0}$?
The closest thing I found so far is Isserlis' theorem, but it still seems to handle a slightly different case.
We can find $a,b$ such that $z_i$ and $az_i+bz_j$ are independent: just make the covariance $0$. [Since they are jointly normal this is enough to say they are independent]. Now $Ez_i^{2}(az_i+bz_j)=Ez_i^{2}E(az_i+bz_j)$. The right side can be computed and the left side is $Ez_i^{2}(az_i+bz_j)=aEz_i^{3}+bEz_i^{2} z_j$. From this you can compute $Ez_i^{2}z_j$.