# Variance of product of k correlated random variables

What is the variance of the product of $k$ correlated random variables?

More information on this topic than you probably require can be found in Goodman (1962): "The Variance of the Product of K Random Variables", which derives formulae for both independent random variables and potentially correlated random variables, along with some approximations. In an earlier paper (Goodman, 1960), the formula for the product of exactly two random variables was derived, which is somewhat simpler (though still pretty gnarly), so that might be a better place to start if you want to understand the derivation.

For completeness, though, it goes like this.

### Two variables

Assume the following:

• $$x$$ and $$y$$ are two random variables
• $$X$$ and $$Y$$ are their (non-zero) expectations
• $$V(x)$$ and $$V(y)$$ are their variances
• $$\delta_x = (x-X)/X$$ (and likewise for $$\delta_y$$)
• $$D_{i,j} = E \left[ (\delta_x)^i (\delta_y)^j\right]$$
• $$\Delta_x = x-X$$ (and likewise for $$\Delta_y$$)
• $$E_{i,j} = E\left[(\Delta_x)^i (\Delta_y)^j\right]$$
• $$G(x)$$ is the squared coefficient of variation: $$V(x)/X^2$$ (likewise for $$G(Y)$$)

Then: $$V(xy) = (XY)^2[G(y) + G(x) + 2D_{1,1} + 2D_{1,2} + 2D_{2,1} + D_{2,2} - D_{1,1}^2]$$ or equivalently:

$$V(xy) = X^2V(y) + Y^2V(x) + 2XYE_{1,1} + 2XE_{1,2} + 2YE_{2,1} + E_{2,2} - E_{1,1}^2$$

### More than two variables

The 1960 paper suggests that this an exercise for the reader (which appears to have motivated the 1962 paper!).

The notation is similar, with a few extensions:

• $$(x_1, x_2, \ldots x_n)$$ be the random variables instead of $$x$$ and $$y$$
• $$M = E\left( \prod_{i=1}^k x_i \right)$$
• $$A = \left(M / \prod_{i=1}^k X_i\right) - 1$$
• $$s_i$$ = 0, 1, or 2 for $$i = 1, 2, \ldots k$$
• $$u$$ = number of 1's in $$(s_1, s_2, \ldots s_k)$$
• $$m$$ = number of 2's in $$(s_1, s_2, \ldots s_k)$$
• $$D(u,m) = 2^u - 2$$ for $$m=0$$ and $$2^u$$ for $$m>1$$,
• $$C(s_1, s_2, \ldots, s_k) = D(u,m) \cdot E \left( \prod_{i=1}^k \delta_{x_i}^{s_i} \right)$$
• $$\sum_{s_1 \cdots s_k}$$ indicates summation of the $$3^k - k -1$$ sets of $$(s_1, s_2, \ldots s_k)$$ where $$2m + u > 1$$

Then, at long last:

$$V\left(\prod_{i=1}^k x_i\right) = \prod X_i^2 \left( \sum_{s_1 \cdots s_k} C(s_1, s_2 \ldots s_k) - A^2\right)$$

See the papers for details and slightly more tractable approximations!

• please note, that the above answer from Matt Krause contains a mistake as well as the paper itself. In the definition of the function C(s1,...,sk) it should be a product instead of a sum. – Nicolas Gisler Nov 19 '15 at 13:59
• Could you elaborate a little bit more..? "Because I - an anonymous person from the Internet - say so" is not really an answer... – Tim Nov 19 '15 at 14:02
• If you try to get the variance var(x*y) for independent random variables, via the formula for arbitrary k you can see that only a product and not a sum gives you the correct answer. In addition, if you look at the paper you can see it as well, on page 59 of the paper (at least in my version) he used a product instead of a sum. – Nicolas Gisler Nov 19 '15 at 14:08
• For the case of two random variables, an easier-to-read formula for the variance of the product of two correlated random variables can be found in this answer by @macro. This answer also points out the essential problem in $$V(xy) = X^2V(y) + Y^2V(x) + 2XYE_{1,1} + 2XE_{1,2} + 2YE_{2,1} + E_{2,2} - E_{1,1}^2,$$ viz., the thicket of notation conceals the essential fact that there are terms in it whose value cannot be determined unless we know cov$(x^2,y^2)$, or enough about the joint density of the two random variables to determine this quantity. – Dilip Sarwate Nov 19 '15 at 15:46
• An edit suggestion, that should really have been a comment, suggested that the original paper contained a typo where a sum and product were mixed up and this answer should be amended. See stats.stackexchange.com/review/suggested-edits/83662 – Silverfish Nov 19 '15 at 16:23

Just to add to the awesome answer of Matt Krause (in fact easily derivable from there). If x, y are independent then, \begin{equation*} \begin{split} E_{1,1} &= E[(x-E[x])(y-E[y])] = Cov(x,y) = 0\\ E_{1,2} &= E[(x-E[x])(y-E[y])^2] \\ &= E[x-E(x)]E[(y-E[y])^2] \\ &= (E[x]-E[x])E[(y-E[y])^2]=0\\ E_{2,1} &= 0\\ E_{2,2} &= E[(x-E[x])^2(y-E[y])^2]\\ &= E[(x-E[x])^2]E[(y-E[y])^2\\ &= V[x]V[y]\\ V[xy] &= E[x]^2 V[y] + E[y]^2 V[x] + V[x]V[y] \end{split} \end{equation*}

• The result for the case of $n$ independent random variables has been discussed here. – Dilip Sarwate Nov 19 '15 at 15:24

In addition to the general formula given by Matt it may be worth noting that there is a somewhat more explicit formula for zero mean Gaussian random variables. It follows from Isserlis' theorem, see also Higher moments for the centered multivariate normal distribution.

Suppose that $(x_1, \ldots, x_k)$ follows a multivariate normal distribution with mean 0 and covariance matrix $\Sigma$. If the number of variables $k$ is odd, $E\left(\prod_i x_i\right) = 0$ and
$$V\left(\prod_i x_i\right) = E\left( \prod_i x_i^2\right) = \sum \prod \tilde{\Sigma}_{i,j}$$ where $\Sigma \prod$ means sum over all partitions of $\{1, \ldots, 2k\}$ into $k$ disjoint pairs $\{i, j\}$ with each term being a product of the corresponding $k$ $\tilde{\Sigma}_{i,j}$'s, and where $$\tilde{\Sigma} = \left( \begin{array}{cc} \Sigma & \Sigma \\ \Sigma & \Sigma \end{array} \right)$$ is the covariance matrix for $(x_1, \ldots, x_k, x_1, \ldots, x_k)$. If $k$ is even, $$V\left(\prod_i x_i\right) = \sum \prod \tilde{\Sigma}_{i,j} - \left(\sum \prod \Sigma_{i,j}\right)^2.$$ In the case $k = 2$ we get $$V(x_1x_2) = \Sigma_{1,1} \Sigma_{2,2} + 2 (\Sigma_{1,2})^2 - \Sigma_{1,2}^2 = \Sigma_{1,1} \Sigma_{2,2} + (\Sigma_{1,2})^2.$$ If $k = 3$ we get $$V(x_1x_2x_3) = \sum \Sigma_{i,j}\Sigma_{k,l}\Sigma_{r,t},$$ where there are 15 terms in the sum.

It is, in fact, possible to implement the general formula. The most difficult part appears to be the computation of the required partitions. In R, this can be done with the function setparts from the package partitions. Using this package it was no problem to generate the 2,027,025 partitions for $k = 8$, the 34,459,425 partitions for $k = 9$ could also be generated, but not the 654,729,075 partitions for $k = 10$ (on my 16 GB laptop).

A couple of other things are worth noting. First, for Gaussian variables with non-zero mean it should be possible to derive an expression as well from Isserlis' theorem. Second, it is unclear (to me) if the above formula is robust against deviations from normality, that is, if it can be used as an approximation even if the variables are not multivariate normally distributed. Third, though the formulas above are correct, it is questionable how much the variance tells about the distribution of the products. Even for $k = 2$ the distribution of the product is quite leptokurtic, and for larger $k$ it quickly becomes extremely leptokurtic.

• Neat approach! For what it's worth, the formula in my answer also has a combinatorial blow-up: the summation over C involves summing $O(3^k)$ terms. – Matt Krause Nov 22 '15 at 2:47