# 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)$?

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?

• You assume the RVs are independent/uncorrelated. The joint distribution should have been in its general form consisting of the general covariance matrix etc. Commented May 18, 2023 at 14:30
• @gunes thank you. My statistics is a little rusty... could you help link to the full general formula for the multivariate joint distribution, please? Commented May 18, 2023 at 14:33
• It is crucial that $(x_1,x_2,x_3,x_4)$ follow a zero-mean multivariate Normal distribution!
– whuber
Commented May 18, 2023 at 15:05
• @whuber thank you, i am still trying to work this out. I gather that what I have attempted above was for the special/trivial case of mutually independent normally-distributed random variables? And I assume that Isserli's theorem applies to this independent variables case as well? But according to my workings, $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)$ so wouldn't this yield $E(x_1x_2)E(x_3x_4) + E(x_1x_3)E(x_2x_4) + E(x_1x_4)E(x_2x_3) = 3E(x_1x_2x_3x_4)$ contrary to the claim of the theorem? Commented May 18, 2023 at 15:15
• When the variables are independent and have zero means, all the expectations factor and every expectation is (by assumption) zero, so you have demonstrated that sums and products of zeros are zero. Without loss of generality you can standardize all the variables (both sides of the equation are multiplied by the same constant), whence the bivariate expectations are the correlation coefficients. The content of this theorem is that a particular fourth multivariate moment is a particularly simple function of those correlations.
– whuber
Commented May 18, 2023 at 15:17

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.

• thank you very much! I understand now the steps following $E(x_1x_2...x_n)=\sum_\pi\prod_B\kappa(B)$ to $\sum_\pi\prod_B\kappa(B)$ being nonzero only for those terms with no singletons and no partitions containing more than 2 members, thus leaving only $E(x_1x_2...x_n)=E(x_1x_2)E(x_3x_4) + E(x_1x_3)E(x_2x_4) + E(x_1x_4)E(x_2x_3)$. I don't fully understand how the combinatorial interpretation of the coefficients of the cumulant generating function comes about in $\sum_\pi\prod_B\kappa(B)=E(x_1x_2...x_n)$. Could you elaborate in a few lines how this combinatorial equivalence arises, please? Commented May 18, 2023 at 23:10
• That's a somewhat complicated subject. It merits an entire chapter in Kendall & Stuart's classic treatise on statistics, for instance. But it really comes down to the rules for working with power series along with the (very simple) power series for exp and log.
– whuber
Commented May 19, 2023 at 14:11
• @whuber, are you by chance aware of a multivariate (bivariate is enough) distribution with tractable fourth order cumulant (other than the "trivial" multivariate normal)? Commented Oct 17, 2023 at 9:54
• @ChristophHanck Could you elaborate on what you mean by "tractable"?
– whuber
Commented Oct 17, 2023 at 11:49
• Thanks, fair enough - is "closed form and ideally not too complicated a function of, say, the underlying parameters of the distribution" better? Essentially, I want to replace the multivariate normal in this answer with some other (mean zero if possible) distribution stats.stackexchange.com/questions/628868/… but haven't found anything yet. Commented Oct 17, 2023 at 12:11

$$\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.

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.