Expectation of the product of iid random variables If we have iid random variables $X_1,X_2,...,X_N$ with $\mathbb{E}X_i=\mu$, is it true that $\mathbb{E}\prod X_i=\mu^N$?
I had no doubt that this is true, until I tried it out with Python, using random.normalvariate() to generate the set of samples, and surprisingly found that the product of all these data points are generally a lot smaller than $\mu^N$.
For example, I used that function to generate 2 million (there's absolutely no need to have a dataset of this size, but I went for it anyways) data points that are, supposedly, distributed as $N(1, 0.2)$. I was hoping for their product to scatter around 1 as I repeat the trial, but instead I got numbers $\sim\pm10^{-18500}$ constantly.
For what it's worth, I've tried sample sizes ranging from 1 to 1000, and they all fell below their respective $\mu^N$, significantly -- the difference was visible on a log-scaled plot.
I suspected random.normalvariate() generated something whose PDF is not $N(1,0.2)$. But I plotted that 2 million data points and got a perfect bell-shape curve.
I suspected that the data are correlated among themselves. But I tried to compute $\mathbb{E}\prod X_i$ with $\text{corr}(X_i,X_j)=\rho$ and found that my calculation could not explain it. I'm not hundred percent confident in my calculation though.
And I tried to understand it intuitively, and had the following thought. Say we have a lot of $N(1,0.2)$ data points. It is conceivable that they lie symmetrically around 1. So, we can group the data points into pairs that are roughly like $\{(1-a_n),(1+a_n)\}$. This should be feasible when the sample size is large enough. But each pair has their product being less than 1. Therefore, the total product is a lot smaller than 1.
So, this seems to me a paradox. I can't dissuade myself from $\mathbb{E}\prod X_i=\mu^N$, but neither can I find the loophole in the thought above (or my empirical tests). I feel that I must have made a blatant mistake somewhere, but I can't locate that error. Please help me out if you know the answer!
 A: There are several things to note here.

*

*Multiplying lots of numbers in computer leads to rounding errors, especially if these numbers vary in magnitude, i.e. small numbers are multiplied by large numbers.

*You need to sample the $\prod X_i$, not $X_i$.

*In simulations you rely on law of large numbers. However the bigger the variance of the random variable you sample, the slower the convergence.

To see the point 3 let us calculate the variance of $\prod X_i$:
$$
Var(\prod X_i) = E(\prod X_i)^2 - (E\prod X_i)^2 = (EX_i^2)^N - (EX_i)^N. 
$$
For normal variable $N(\mu, \sigma^2)$ this expression turns into:
$$
(\sigma^2+\mu^2)^N-\mu^N \approx N\sigma^2
$$
if $\mu=1$. So the variance scales with $N$.
So here is an example in R which showcases this:
> set.seed(666)
> prod_sample <- sapply(1:1000, function(x)prod(rnorm(100, mean = 1, sd = sqrt(0.2))))
> quantile(prod_sample)
           0%           25%           50%           75%          100% 
-3.001294e+00 -5.287966e-07  5.889490e-11  2.814698e-06  7.951951e+01     
> mean(prod_sample)
[1] 0.1500762

This is for $N=100$ and the product is sampled 1000 times.
If we chose $N=3$, then the result is more aligned with the expectations:
> prod_sample <- sapply(1:1000, function(x)prod(rnorm(3, mean = 1, sd = sqrt(0.2))))
> quantile(prod_sample)
        0%        25%        50%        75%       100% 
-0.9823311  0.4047946  0.7918063  1.3744908  5.5923739 
> 
> mean(prod_sample)
[1] 1.000404

A: First, let's establish the correct identity.
When $X_1, \ldots, X_N$ are independent variables with finite expectations $\mu_i=X_i,$ then by laws of conditional expectation,
$$E\left[\prod_{i=1}^N X_i\right] =E\left[X_N E\left[\prod_{i=1}^{N-1} X_i \mid X_{N}\right]\right] = E\left[X_N \prod_{i=1}^{N-1} \mu_i\right] = \mu_N\prod_{i=1}^{N-1} \mu_i= \prod_{i=1}^N \mu_i$$
gives a proof by mathematical induction (beginning with the base case $N=1$ where $$E\left[\prod_{i=1}^N X_i\right] = E\left[X_1\right] = \mu_1 = \prod_{i=1}^N \mu_i$$ is trivially true).
Now, let's find an explanation for the simulation results.
The pairing argument in the question is an interesting one, because it shows that when multiplications by $(1-a)$ and $(1+a)$ occur in equal numbers (approximately $N/2$ each), the net product is $(1-a^2)^{N/2}\approx \exp(-Na^2/2).$  This suggests that when $N$ is sufficiently large, it's nearly certain that the product will be tiny--certainly less than the common mean of $1.$  The reason this is not a paradox is that there will be a vanishingly small--but still positive--probability of yielding a whopping big number on the order of $(1+a)^N \approx \exp(aN).$  This rare chance of a huge product balances out all the tiny products, keeping the mean at $1.$
It is not easy to analyze the product of many Normal variables.  Instead, we may gain insight from a simpler case.  Let $Y_1, Y_2, \ldots,$ be a sequence of independent Rademacher variables: that is, each of these has a $1/2$ chance of being either $1$ or $-1.$  Pick some number $0 \lt a \lt 1$ and define $X_i = 1 + aY_i,$ so that each $X_i$ has equal chances of being $1\pm a.$ Clearly $E[X_i] = 1 = \mu_i$ for all $i.$
Consider the product of the first $N$ of these $X_i.$  Suppose, in a simulation, that $k$ of these values equal $1-a$ and (therefore) the remaining $N-k$ of them equal $1+a.$  The product then is $(1-a)^k(1+a)^{N-k}.$ How small must $k$ be for this product to exceed $1$?
Given $N$ and $a,$ we must solve the inequality
$$(1-a)^k(1+a)^{N-k} \ge 1$$
for $k.$  By taking logarithms, this is equivalent to
$$k \le N \frac{\log(1+a)}{\log(1+a) - \log(1-a)}.$$
Because each $X_i$ has equal and independent chances of being $1\pm a,$ the distribution of $k$ is Binomial$(N, 1/2),$ which even for moderate sizes of $N$ ($N \ge 10$ is fine) is nicely approximated by a Normal$(N/2, \sqrt{N}/2)$ distribution.  Thus, the chance that the product is $1$ or greater will be close to the value of the standard Normal distribution at $Z$ (the tail area under the Bell Curve left of $Z$) where
$$Z = \frac{N \frac{\log(1+a)}{\log(1+a) - \log(1-a)} - \frac{N}{2}}{\sqrt{N}/2} = \text{constant}\times \sqrt{N}.$$
You can see where this is going!  As $N$ grows large, $Z$ is pushed further out to the left, making it less and less likely to observe any product greater than $1$ in a simulation.
In the question, where the standard deviation is $0.2,$ the value $a=0.2$ will closely reproduce the simulation behavior.  In this case the constant is
$$\text{constant} = \frac{2\log(1+a)}{\log(1+a)-\log(1-a)} - 1 = -0.100\ldots$$
Taking $N=2\times 10^6,$ for instance, as in the question, compute $Z \approx -142.$  The chance that $k$ is small enough to produce a value this negative is less than $10^{-10000}.$  You can't even represent that in double precision floats.  It would take far more than the age of the universe to create a simulation that had the remotest chance of producing such an imbalance between the $1+a$ and $1-a$ values that the product exceeds $1.$
In short, for all practical purposes, when $N$ is sufficiently large ($N \gg 5000$ will do when $a=1/5$), you will never observe a value above $1$ in this simulation, even though the mean of the product is $1.$
A: As an addendum to @whuber's explanation, consider that
\begin{align}
\prod_{i=1}^N \vert X_i\vert &= \exp\{ \sum\nolimits_{i=1}^N \overbrace{\ln \vert X_i\vert}^{Y_i} \}\\
&= \exp\{ \sum\nolimits_{i=1}^N Y_i \}\\
&= \exp\{ \sum\nolimits_{i=1}^N (Y_i - \tilde \mu)+N\tilde\mu\}\qquad\qquad\mathbb \ \ E[Y_i]=\tilde\mu\\
&= \exp\left\{\frac{\tilde\sigma\sqrt{N}}{\sqrt{N}\tilde\sigma} 
\sum\nolimits_{i=1}^N (Y_i - \tilde \mu)+N\tilde\mu\right\}\qquad\text{var}(Y_i)=\tilde\sigma^2\\
&\approx \exp\{\tilde\sigma\sqrt{N}\zeta+N\tilde\mu\}\qquad\quad\qquad\qquad\zeta\sim\mathcal N(0,1)\\
\end{align}
which shows that the expression tends to behave more and more erratically as $N$ increases. (For $|\mu|$ large enough, $\tilde \mu\approx\log|\mu|$.) Note that
$$\xi=\exp\{\tilde\sigma\sqrt{N}\zeta+N\tilde\mu\}$$
is a log-normal $\mathcal L\mathcal N(N\tilde\mu,N\tilde\sigma^2)$
random variable with
$$\mathbb E[\xi]=\exp\{N\tilde\mu+N\tilde\sigma^2/2)\qquad
\text{var}(\xi)=[\exp\{N\tilde\sigma^2/2\}-1]\exp\{2N\tilde\mu+N\tilde\sigma^2\}$$
As a last remark, not quite related with the spirit of the question, while $$\prod_{i=1}^NX_i$$ is indeed an unbiased estimator of $\mu^N$, a more efficient estimator based on the $N$ rv's $X_i$ would be
$$\mathbb E\left[\prod_{i=1}^NX_i\Big\vert\bar X_N\right]$$
thanks to the Rao-Blackwell theorem. I however do not see an easy way out for computing this conditional expectation (which should be a polynomial of degree $N$ in $\bar X_N$). But since $\bar X_N$ is minimal sufficient and complete there exists a single (UMVU) estimator.
