Median and mean of the sample mean of i.i.d. log-normal Let $y:=\frac1n\sum_{i=1}^n x_i$, where $\{x_i\}_{i=1}^n$ is a set of i.i.d. random variables, and every $x_i$ has a lognormal distribution $x_i \sim\text{Lognormal}(\mu,\sigma^2)$. Let $\text{Med}[y]$ be the median of $y$. Is the following inequality true $\forall (n,\mu,\sigma)$?
$$\text{Med}[y]<\mathbf E[y]$$ 

Motivation: I am computing the sample mean of the lognormal random variables via Monte Carlo. The sample mean seems tend to concentrate below the mean for large $\sigma$. I am wondering whether this is true for all cases. It is true for a sigle sample. However, there is no explicit formula for the distribution of the mean of finite number of --- not even two -- samples i.i.d. lognormal variables. I have no idea how to prove it.
 A: I would like to provide an asymptotic solution, as the sample size $n\rightarrow\infty$.
Using the Fenton-Wilkson empirical asymptotic result in [1], we know that the sample sum of i.i.d. lognormals $$X_i\sim \text{LN}(\mu,\sigma^2)$$ can be approximated by lognormal with different mean and variance $Y=\sum_{i=1}^nX_i$ is distributed as $$Y\sim \text{LN}(\mu_n,\sigma^2_n)$$ where $\sigma_{n}^{2}=\log\left(\frac{e^{\sigma^{2}}-1}{n}+1\right)$  and $\mu_{n}=\log\left(n\cdot e^{\mu}\right)+\frac{1}{2}\left(\sigma_{n}^{2}-\sigma^{2}\right)$.
We know from properties of lognormal distributions that $\text{Median}(Y)=e^{\mu_n}$ and $\text{Mean}(Y)=\exp\big(\mu_n+\frac{\sigma^2_n}{2}\big)$and therefore $\text{Median}(Y)\le \text{Mean}(Y)$ is asymptotically proved since $\sigma^2_n\rightarrow 0$. However this argument holds only when the Fenton-Wilkson approximation works within the order of accuracy of $O\Big(e^{\frac{\sigma^2_n}{2}}\Big)$, which is not always the case. If we know the exact distribution of $Y$, this argument can be modified to see if such an argument holds in general.
[1] Cobb, Barry R., R. Rumi, and Antonio Salmerón. "Approximating the distribution of a sum of log-normal random variables." Statistics and Computing 16.3 (2012): 293-308. http://leo.ugr.es/pgm2012/submissions/pgm2012_submission_6.pdf
A: Here is a positive result for the average of two iid lognormals.  Recall that the average of a lognormal distributed as $LN(\mu,\sigma)$ is $\exp(\mu+\sigma^2/2)$.

Let $y_2$ be the average of two iid lognormals each distributed as
  $LN(\mu,\sigma)$.  Then:
For small $\sigma$, $Med[y_2] \approx \exp(\mu+\sigma^2/4)$, which is less
  than $E[y_2]$.
For large $\sigma$, $Med[y_2] \approx \exp(\mu+.545\sigma)$, which is less
  than $E[y_2]$.
From those we can establish that $Med[y_2] < E[y_2]$ for all positive $\sigma$.

The approximations are valid to second order in $\sigma$, in the sense that
$$P\left[y_2 < \exp(\sigma^2/4)\right] = \frac{1}{2} + O(\sigma^3)$$
$$P\left[y_2 < \exp(0.545\sigma)\right] = \frac{1}{2} + O(\frac{1}{\sigma^3})$$
where .545 is an abbreviation for $\ln(F_1^{-1}(\frac{1}{\sqrt{2}}))$, and $F_{\sigma}$ is the CDF of $LN[0,\sigma]$, and we take $\mu=0$ for simplicity.
Proof:
We prove the first approximation from the general convolution formula
$\begin{align}
P[y_2 < m]
 &= \int\int_{\frac{x_1+x_2}2<m}dF(x_1)dF(x_2)\ \ \text{for arbitrary i.i.d. with CDF }F \\
 &= \int_{x_1} F(2m-x_1)\, dF(x_1)\ \ \text{by integrating over }x_2 \\
 &= \int_{0}^{1} F\big(2m - F^{-1}(p)\big)dp\ \ \text{where } p=F(x_1) \\
\end{align}$
which in this particular case gives
$$P[y_2 < m] = \int_{0}^{F_\sigma(2m)} F_\sigma\big(2m - F_\sigma^{-1}(p)\big)dp.$$
$F_\sigma(2m)$ is the largest possible $p$ such that $(x_1+x_2)/2<m$ for some positive $x_2$.

The graph shows the line $1-p$ and the integrand for $\sigma=1/2, \ m=\exp(\sigma^2/4)$.  The integrand approaches the line as $\sigma$ goes to 0.  The median is $m$ iff the areas under the two curves are equal.
We use Mathematica to differentiate the formula for $P[y_2 < m]$ under the integral sign and at the limit of integration, and this computes that the coefficients in the power series are $1/2, 0, 0$.
F[u_] = Simplify[CDF[LogNormalDistribution[0, sigma], u], Assumptions -> u > 0]
x1 = u /. Solve[p == F[u], u][[1]]
x2 = u /. Solve[q == F[u], u][[1]]
qsol = q /. Solve[(x1 + x2)/2 == Exp[sigma^2/4], q][[1]]
plim = F[2 Exp[sigma^2/4]]
Limit[plim, sigma -> 0]
Limit[D[plim, sigma], sigma -> 0]
Limit[D[plim, {sigma, 2}], sigma -> 0]
Integrate[Limit[qsol, sigma -> 0], {p, 0, 1}]
Integrate[Limit[D[qsol, sigma], sigma -> 0], {p, 0, 1}]
Integrate[Limit[D[qsol, {sigma, 2}], sigma -> 0], {p, 0, 1}]
Plot[{qsol /. sigma -> 1/2, 1 - p, 1}, {p, 0, 1}, PlotRange -> {0, 1}, 
 AxesLabel -> {"p=F[x1]", "q=F[x2]"}, PlotStyle -> {, , {Black, Thin}}]

We prove the second approximation from
$$\frac{1}{2} < P{\Large[}y_2 < e^{0.545\sigma}{\Large]} < F_\sigma(2 e^{0.545\sigma})^2\ \text{ and } 
\lim_{\sigma\rightarrow\infty}F_\sigma(2 e^{0.545\sigma})^2=\frac{1}{2}$$
and the fact (not proved here) that $P[y_2 < e^{0.545\sigma}]$ has eventually monotonic derivatives.
The inequality comes from 
\begin{align}
\frac{1}{2}
&= F_\sigma(e^{0.545\sigma})^2 \\
&= P[x_1 < e^{0.545\sigma}\ \& \ x_2 < e^{0.545\sigma}] \\
&< P[\ \ \ \ \ \ \ x_1 + x_2 < 2e^{0.545\sigma} \ \ \ \ \ \ \ \ \, ] \\
&< P[x_1 < 2e^{0.545\sigma}\ \& \ x_2 < 2e^{0.545\sigma}] \\
&= F_\sigma(2 e^{0.545\sigma})^2 \\
\end{align}
We use Mathematica to show that the limit of the right-hand side is 1/2.
x3 = u /. Solve[1/Sqrt[2] == F[u], u][[1]]
Limit[F[2 x3]^2, sigma -> Infinity] // FullSimplify
N[x3]

