Q:

Bias of moment estimator of lognormal distribution

I am doing some numerical experiment that consists in sampling a lognormal distribution $X\sim\mathcal{LN}(\mu, \sigma)$, and trying to estimate the moments $\mathbb{E}[X^n]$ by two methods:

Looking at the sample mean of the $X^n$
Estimating $\mu$ and $\sigma^2$ by using the sample means for $\log(X), \log^2(X)$, and then using the fact that for a lognormal distribution, we have $\mathbb{E}[X^n]=\exp(n \mu + (n \sigma)^2/2)$.

The question is:
I find experimentally that the second method performs much better then the first one, when I keep the number of samples fixed, and increase $\mu, \sigma^2$ by some factor T. Is there some simple explanation for this fact? 
I'm attaching a figure in which the x-axis is T, while the y axis are the values of $\mathbb{E}[X^2]$ comparing the true values of $\mathbb{E}[X^2] = \exp(2 \mu + 2 \sigma^2)$ (orange line), to the estimated values. method 1 - blue dots, method 2 - green dots. y-axis is in log scale

EDIT:
Below is a minimal Mathematica code to produce the results for one T, with the output:
   ClearAll[n,numIterations,sigma,mu,totalTime,data,rmomentFromMuSigma,rmomentSample,rmomentSample]
(* Define variables *)
n=2; numIterations = 10^4; sigma = 0.5; mu=0.1; totalTime = 200;
(* Create log normal data*)
data=RandomVariate[LogNormalDistribution[mu*totalTime,sigma*Sqrt[totalTime]],numIterations];

(* the moment by theory:*)
rmomentTheory = Exp[(n*mu+(n*sigma)^2/2)*totalTime];

(*Calculate directly: *)
rmomentSample = Mean[data^n];

(*Calculate through estimated mu and sigma *)
muNumerical = Mean[Log[data]]; (*numerical \[Mu] (gaussian mean) *)
sigmaSqrNumerical = Mean[Log[data]^2]-(muNumerical)^2; (* numerical gaussian variance *)
rmomentFromMuSigma = Exp[ muNumerical*n + (n ^2sigmaSqrNumerical)/2];

(*output*)
Log@{rmomentTheory, rmomentSample,rmomentFromMuSigma}

Output:
(*Log of {analytic, sample mean of r^2, using mu and sigma} *)
{140., 91.8953, 137.519}

above, the second result is the sample mean of $r^2$, which is below the two other results

A:

There is something puzzling in those results since

the first method provides an unbiased estimator of $\mathbb{E}[X^2]$, namely$$\frac{1}{N}\sum_{i=1}^N X_i^2$$has $\mathbb{E}[X^2]$ as its mean. Hence the blue dots should be around the expected value (orange curve);
the second method provides a biased estimator of $\mathbb{E}[X^2]$, namely$$\mathbb{E}[\exp(n \hat\mu + n^2 \hat{\sigma}^2/2)]>\exp(n \mu + (n \sigma)^2/2)$$when $\hat\mu$ and $\hat\sigma²$ are unbiased estimators of $\mu$ and $\sigma²$ respectively, and it is thus strange that the green dots are aligned with the orange curve.

but they are due to the problem and not to the numerical computations: I repeated the experiment in R and got the following picture with the same colour code and the same sequence of $\mu_T$'s and $\sigma_T$'s, which represents each estimator divided by the true expectation:

Here is the corresponding R code:
moy1=moy2=rep(0,200)
mus=0.14*(1:200)
sigs=sqrt(0.13*(1:200))
tru=exp(2*mus+2*sigs^2)
for (t in 1:200){
x=rnorm(1e5)
moy1[t]=mean(exp(2*sigs[t]*x+2*mus[t]))
moy2[t]=exp(2*mean(sigs[t]*x+mus[t])+2*var(sigs[t]*x+mus[t]))}

plot(moy1/tru,col="blue",ylab="relative mean",xlab="T",cex=.4,pch=19)
abline(h=1,col="orange")
lines((moy2/tru),col="green",cex=.4,pch=19)

Hence there is indeed a collapse of the second empirical moment as $\mu$ and $\sigma$ increase that I would attribute to the enormous increase in the variance of the said second empirical moment as $\mu$ and $\sigma$ increase.

My explanation of this curious phenomenon is that, while $\mathbb{E}[X^2]$ obviously is the mean of
  $X^2$, it is not a central value: actually the median of $X^2$ is equal to $e^{2\mu}$. When representing the random variable $X^2$ as $\exp\{2\mu+2\sigma\epsilon\}$ where $\epsilon\sim\mathcal{N}(0,1)$, it is clear that,
  when $\sigma$ is large enough, the random variable $\sigma\epsilon$ is almost never of the
  magnitude of $\sigma^2$. In other words if
  $X$ is $\mathcal{LN}(\mu,\sigma)$ $$\begin{align*}\mathbb{P}(X^2>\mathbb{E}[X^2])&=\mathbb{P}(\log\{X^2\}>2\mu+2\sigma^2)\\&=\mathbb{P}(\mu+\sigma\epsilon>\mu+\sigma^2)\\&=\mathbb{P}(\epsilon>\sigma)\\
&=1-\Phi(\sigma)\end{align*}$$
  which can be arbitrarily small.

A:

I thought I'd throw up some figs showing that both user29918 and Xi'an's plots are consistent. Fig 1 plots what user29918 did, and Fig 2 (based on same data), does what Xi'an did for his plot. Same result, different presentation.
What's happening is that as T increases, the variances becomes huge and the estimator $\frac{1}{n} \sum_i x_i^2$ becomes like trying to estimate the population mean of the Powerball Lotto by buying Lotto tickets! A large percentage of the time, you will underestimate the payoff (because no sample observation hits the jackpot) and a tiny percentage of the time, you will massively overestimate the payoff (because there's a jackpot winner in the sample). The sample mean is an unbiased estimate but it's not expected to be precise, even with thousands and thousands of draws! In fact, as it becomes harder and harder to win the lotto, your sample mean will be below the population mean the vast majority of the time.
Further Comments:

An unbiased estimator does not mean the estimator is expected to be close! The blue dots need not be near the expectation. Eg. a single observation chosen at random gives an unbiased estimate of the population mean, but that estimator would not be expected to be close.
The issue is coming up as the variance is becoming absolutely astronomical. As the variance goes batshit, the estimate for the first method is being driven be just a few observations. You also start having a tiny, tiny probability of an INSANELY, INSANELY, INSANELY big number... 
This is an intuitive explanation. Xi'an has a more formal derivation. His result $P(X^2 > E[X^2]) = 1 - \Phi(\sigma)$ implies that as $\sigma$ gets large, it becomes incredibly unlikely to ever draw an observation above the mean, even with thousands of observations. My language of "winning the lotto" refers to an event where $X^2 > E[X^2]$.

