Variance of product of k correlated random variables

Question

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

Answer 1

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!

Answer 2

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*}

Answer 3

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.