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:

1. Looking at the sample mean of the $X^n$
2. 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

• An unbiased estimator does not imply that the blue dots should be near the expected value (orange curve). An estimator can be unbiased if it has a high probability of being too low and small (perhaps vanishingly small) probability of being way too high. That is what's occurring as T increases and the variance goes batshit huge (see my answer). – Matthew Gunn May 26 '16 at 9:28
• For how to obtain unbiased estimators, please see stats.stackexchange.com/questions/105717. UMVUEs of the mean and variance are given in the answers and comments thereto. – whuber Nov 1 '16 at 19:35

There is something puzzling in those results since

1. 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);
2. 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.

• I'm also puzzled. I'm adding a minimal code with the results (Mathematica) – user29918 May 26 '16 at 8:11
• Ok. Thanks! Putting some numbers, I see now that my meager sample size was really not up for the task! – user29918 May 26 '16 at 14:06
• @user29918: Sorry, I do not see the sample size as the problem but rather the fact that the log-normal becomes very skewed when $\sigma$ grows to infinity as the mean to become useless. – Xi'an May 26 '16 at 14:58
• @Xi'an Good stuff! $P(X^2 > \mathbb{E}[X^2]) = 1 - \Phi(\sigma)$. That captures in precisely in equations what I was (rather imprecisely) trying to express in words, that as $\sigma$ increases, it becomes increasingly likely (and for large $\sigma$, near certain) that an observation is below the mean. Indeed the probability is so high that it's highly likely the entire sample is below the mean! – Matthew Gunn May 26 '16 at 18:50
• This type of asymptotic is not very helpful in that the number of simulations needed to approximate correctly the moments grows exponentially fast with $\sigma$. – Xi'an May 28 '16 at 15:58

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.

3. 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]$.