How to prove Isserli's theorem $ E(x_1x_2x_3x_4) = E(x_1x_2)E(x_3x_4) + E(x_1x_3)E(x_2x_4) + E(x_1x_4)E(x_2x_3) $?

Question

I am having trouble understanding and proving Isserli's theorem for n=4:

$$ E(x_1x_2x_3x_4) = E(x_1x_2)E(x_3x_4) + E(x_1x_3)E(x_2x_4) + E(x_1x_4)E(x_2x_3) $$

My attempt goes as follows:

Assuming

$$ P(x_i) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x_i^2}{2}}$$

Then

$$\small{ E(x_1x_2x_3x_4) = \int_{x_1=-\infty}^\infty \int_{x_2=-\infty}^\infty \int_{x_3=-\infty}^\infty \int_{x_4=-\infty}^\infty (x_1x_2x_3x_4) P(x_1x_2x_3x_4) dx_4dx_3dx_2dx_1}$$

$$\small{= \int_{x_1=-\infty}^\infty \int_{x_2=-\infty}^\infty \int_{x_3=-\infty}^\infty \int_{x_4=-\infty}^\infty x_1x_2x_3x_4 \frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{x_3^2}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{x_4^2}{2}}dx_4dx_3dx_2dx_1}$$

$$=\small{ \int_{x_1=-\infty}^\infty x_1\frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}}\int_{x_2=-\infty}^\infty x_2\frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}}\int_{x_3=-\infty}^\infty x_3\frac{1}{\sqrt{2\pi}}e^{-\frac{x_3^2}{2}}\left(\int_{x_4=-\infty}^\infty x_4 \frac{1}{\sqrt{2\pi}}e^{-\frac{x_4^2}{2}}dx_4\right)dx_3dx_2dx_1}$$

$$\small{= E(x_4) \int_{x_1=-\infty}^\infty x_1\frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}}\int_{x_2=-\infty}^\infty x_2\frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}}\left(\int_{x_3=-\infty}^\infty x_3\frac{1}{\sqrt{2\pi}}e^{-\frac{x_3^2}{2}}dx_3\right)dx_2dx_1}$$

$$\small{= E(x_3)E(x_4) \int_{x_1=-\infty}^\infty x_1\frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}}\left(\int_{x_2=-\infty}^\infty x_2\frac{1}{\sqrt{2\pi}}e^{-\frac{x_2^2}{2}}dx_2\right)dx_1}$$

$$\small{= E(x_2)E(x_3)E(x_4) \int_{x_1=-\infty}^\infty x_1\frac{1}{\sqrt{2\pi}}e^{-\frac{x_1^2}{2}}dx_1}$$

$$\small{= E(x_1)E(x_2)E(x_3)E(x_4)}$$

Since each of $x_1, x_2, x_3, x_4$ is assumed to be Gaussian independent random variables, I can't see where does the above calculation go wrong in the context of Isserli's theorem?

Answer 1

Isserli's Theorem equates various multivariate cumulants of a zero-mean multivariate Normal distribution. However, most of this theorem is a universal result, applicable to all multivariate distributions (with finite relevant moments), which takes on a particularly simple form for zero-mean distributions and an even simpler form for zero-mean Normal distributions. These simplifications arise from the very definitions of "zero-mean" (section 2 below) and "Multivariate Normal," (section 3 below) because those definitions essentially specify that particular cumulants are zero.

---

1. Generalities

There is a fully general relationship between moments and cumulants for any $n$-variate distribution of variables $(X_1,X_2,X_3,\ldots,X_n)$ (whose relevant moments are defined and finite).

Given such a distribution and for any subset of the indexes $B\subset\{1,2,3,\ldots,n\},$ $B = \{t_{1_1},t_{i_2},\ldots, t_{i_p}\},$ let $\kappa(B)$ be the coefficient of $i^pt_{i_1}t_{i_2}\cdots t_{i_p}$ (where $i^2 = -1$) in the power series expansion of the expectation

$$K(t_1,t_2,\ldots, t_n) = \log \left(E\left[\exp\left(it_1X_1+it_2X_2+\cdots +it_nX_n\right)\right]\right).$$

The argument of the log is the characteristic function of the distribution and the left hand side is called the cumulant generating function.

Then

$$E[X_1X_2X_3\cdots X_n] = \sum_{\pi}\prod_{B\in\pi}\kappa(B)\tag{*}$$

where $\pi$ ranges over all partitions of $\{1,2,\ldots, n\}.$ There are, for example, 15 partitions of $\{1,2,3,4\}.$ Four of them appear in the formulas below. The other 11 partitions all include singleton subsets, such as $\{\{1\},\ \{2,3,4\}\},$ $\{\{1\},\ \{2\},\ \{3,4\}\},$ and $\{\{1\},\ \{2\},\ \{3\},\ \{4\}\}.$

This is a purely formal combinatorial result having essentially nothing to do with distributions or multivariate Normal variables. It's a matter of manipulating formal power series for the exponential and the logarithm.

In particular, $E[X_i] = \kappa(\{i\})$ are the means and, when $n=2,$ $(*)$ becomes

$$E[X_iX_j] = \kappa(\{i,j\}) + \kappa(\{i\})\kappa(\{j\}) = \kappa(\{i,j\}) + E[X_i]E[X_j]$$

Upon solving for the cumulant this reveals the $\kappa(\{i,j\})$ to be the covariances.

---

2. Zero-mean distributions

Let's now specialize to a zero-mean distribution. This greatly simplifies $(*),$ because any term for a partition containing a singleton will be a multiple of the expectation $\kappa(\{i\})=E[X_i]=0$ and thereby drop out; and the covariances likewise reduce to the expected products, $\operatorname{Cov}(X_i,X_j) = E[X_iX_j] = \kappa(\{i,j\}).$

With $n=4$ the equation $(*)$ reduces to

$$\begin{aligned} &E[X_1X_2X_3X_4]\\ &= \kappa(\{1,2,3,4\}) + \kappa(\{1,2\})\kappa(\{3,4\}) + \kappa(\{1,3\})\kappa(\{2,4\}) + \kappa(\{1,4\})\kappa(\{2,3\})\\ &= \kappa(\{1,2,3,4\}) + E[X_1X_2]E[X_3X_4] + E[X_1X_3]E[X_2X_4] + E[X_1X_4]E[X_2X_3]. \end{aligned}$$

---

3. Multivariate Normal distributions (with zero means)

Finally, the multivariate Normal distribution is defined by having all higher-order cumulants (with three or more $t_{i}$ involved) equal to zero. That is, its cumulant generating function is merely quadratic. Consequently $\kappa(\{1,2,3,4\})=0$ and the $n=4$ case (as well as all other cases) of Isserli's theorem is immediate, QED.

Answer 2

$\newcommand{\bt}{\boldsymbol t}$

Without attempting to the reach the generalities as in whuber's answer, below is a proof for $n = 4$ by directly differentiating the MGF $M(\bt)$ of $(X_1, X_2, X_3, X_4) \sim N_4(0, \Sigma)$. For this quartic Gaussian r.v., with $Q = \Sigma^{-1}$, we have \begin{align} M(\bt) = M(t_1, t_2, t_3, t_4) = \exp(-\bt'Q\bt/2), \end{align} and $E[X_1X_2X_3X_4]$ is given by \begin{align} E[X_1X_2X_3X_4] = \left.\frac{\partial^4 M(\bt)}{\partial t_4\partial t_3\partial t_2\partial t_1}\right|_{\bt = 0}. \tag{1} \end{align}

To evaluate the order 4 partial derivative in $(1)$, we need to apply the following equations repeatedly: \begin{align} & \frac{\partial \bt'Q\bt}{\partial t_i} = e_i'Q\bt, \; i = 1, 2, 3, 4. \tag{2} \\ & \frac{\partial e_i'Q\bt}{\partial t_j} = e_i'Qe_j, \; i, j \in \{1, 2, 3, 4\}. \tag{3} \\ \end{align} In $(2)$, $e_i$ is the $4$-long column vector with $i$-th entry $1$ and all the other entries $0$.

The rest of calculations are based on $(2), (3)$ and the chain rule: \begin{align} & M_1(\bt) := \frac{\partial M(\bt)}{\partial t_1} = -M(\bt)e_1'Q\bt, \tag{4} \\ & M_2(\bt) := \frac{\partial M^2(\bt)}{\partial t_2\partial t_1} = \frac{\partial M_1(\bt)}{\partial t_2} =M(\bt)e_2'Q\bt\cdot e_1'Q\bt - M(\bt)e_1'Qe_2 \tag{5} \end{align} Note that $(5)$ gives the expression of second joint moments $E[X_iX_j] = \left.\frac{\partial M^2(\bt)}{\partial t_j\partial t_i}\right|_{\bt = 0} = e_i'Qe_j$. Continue:
\begin{align} & M_3(\bt) := \frac{\partial^3 M(\bt)}{\partial t_3\partial t_2\partial t_1} = \frac{\partial M_2(\bt)}{\partial t_3} \\ =& -M(\bt)e_3'Q\bt\cdot e_2'Q\bt\cdot e_1'Q\bt + M(\bt)e_2'Qe_3\cdot e_1'Q\bt + M(\bt)e_2'Q\bt\cdot e_1'Qe_3 \\ &+ M(\bt)e_3'Q\bt \cdot e_1'Qe_2. \tag{6} \end{align} Finally, we can get $(1)$ from $(6)$: note that the partial derivative of the first term in the right-hand side of $(6)$ is $0$ when evaluated at $\bt = 0$ because there are more than one factor of the type $e_i'Q\bt$. It thus follows that \begin{align} & E[X_1X_2X_3X_4] = \left.\frac{\partial^4 M(\bt)}{\partial t_4\partial t_3\partial t_2\partial t_1}\right|_{\bt = 0} = \left.\frac{\partial M_3(\bt)}{\partial t_4}\right|_{\bt = 0} \\ =& e_2'Qe_3 \cdot e_1'Qe_4 + e_2'Qe_4 \cdot e_1'Qe_3 + e_3'Qe_4 \cdot e_1'Qe_2 \\ =& E[X_2X_3]E[X_1X_4] + E[X_2X_4]E[X_1X_3] + E[X_3X_4]E[X_1X_2]. \end{align} This completes the proof.

Answer 3

Let $x \sim N(0, \Sigma)$ and let $\Sigma = LL^T$ be the Cholesky decomposition of $\Sigma$, that is, $L$ is lower triangular. Then, we can write $x = L z$ where $z \sim N(0, I_4)$. Now write \begin{align*} \prod_{i=1}^4 x_i = \prod_{i=1}^4 (Lz)_i = \prod_{i=1}^4 \sum_{j=1}^i L_{ij} z_j. \end{align*} Next, switch the order of $\prod$ and $\sum$ (that is expend the product over summation). Then take expectation of both sides and use linearity of expectation and the fact that $z_j$ are zero mean and independent so many terms in the expansion will be zero.