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!

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.
