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!

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

