Expectation of the Truncated Lognormal Distribution

I derive the expected value of a random variable that is left-truncated and lognormally distributed.

Consider a random variable YY which is lognormally distributed with parameters μ\mu and σ\sigma:

Y=exp(X),XN(μ,σ2).(1) Y = \exp(X), \qquad X \sim \mathcal{N}(\mu, \sigma^2). \tag{1}

Now consider the left-truncated variable Y^\hat{Y}, defined as

Y^={Yif Yk>0,0else.(2) \hat{Y} = \begin{cases} Y & \text{if $Y \geq k \gt 0$,} \\ 0 & \text{else.} \end{cases} \tag{2}

The goal of this post is to derive the expected value of Y^\hat{Y}. Intuitively, we might expect this value to be greater than the expected value of YY, since left-truncation eliminates probability mass in the left tail of the distribution (Figure 11).

Figure 1. The probability density function of a random variable which is lognormally distributed with parameters μ=0\mu = 0 and σ=1\sigma = 1 and then truncated at k=0.3k = 0.3 (left) and k=1.5k = 1.5 (right).

In general, if ZZ is a random variable left-truncated at aa and if f(z)f(z) denotes its probability density function (PDF) and F(z)F(z) denotes its cumulative distribution function (CDF), then its expected value is

E[ZZ>a]=azg(z)dz1F(a),(3) \mathbb{E}[Z \mid Z \gt a] = \frac{\int_a^{\infty} z g(z) \text{d}z}{1 - F(a)}, \tag{3}

where

g(z)={f(z)z>k,0else.(4) g(z) = \begin{cases} f(z) & z \gt k, \\ 0 & \text{else}. \end{cases} \tag{4}

I think this makes intuitive sense if we re-arrange the terms. Consider this equation:

azg(z)dz=E[ZZ>a]P(Z>a).(5) \int_a^{\infty} z g(z) \text{d}z = \mathbb{E}[Z \mid Z \gt a] \mathbb{P}(Z \gt a). \tag{5}

Here, the left-hand side is not the desired expectation, as it would be strictly lower than the expected value of ZZ. This would not make sense, since the lower bound on the admissible values of ZZ is actually increasing as aa increases. So we need to adjust this integral by the probability that ZZ is greater than aa.

To compute the expectation in our case, we want to simplify the integral II,

I=logkexp{x}12πσ2exp{12[xμσ]2}dx.(6) I = \int_{\log k}^{\infty} \exp\left\{x\right\} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left\{ -\frac{1}{2} \left[ \frac{x - \mu}{\sigma} \right]^2 \right\} \text{d}x. \tag{6}

The lower bound is logk\log k, not kk, because we have used a change of variables, x=logyx = \log y, in order to express the expectation in terms of the density function of the normal distribution.

We can combine the two exponent terms and then simplify the expression to be again quadratic in xx. We can then separate the terms that are quadratic in xx from the terms that do not depend on xx:

exp{12[xμσ]2+x}=exp{12σ2[x2μ22xμ2σ2x]}=exp{12σ2[x2μ22xμ2σ2x+(2σ2μ+σ4)(2σ2μ+σ4)]}=exp{12σ2[(xμσ2)2(2σ2μ+σ4)]}=exp{12σ2[xμσ2]2}exp{μ+12σ2}.(7) \begin{aligned} &\exp\left\{ -\frac{1}{2} \left[ \frac{x - \mu}{\sigma} \right]^2 + x\right\} \\ &= \exp\left\{ -\frac{1}{2\sigma^2} \left[ x^2 - \mu^2 - 2x\mu - 2 \sigma^2 x \right]\right\} \\ &= \exp\left\{ -\frac{1}{2\sigma^2} \left[ x^2 - \mu^2 - 2x\mu - 2 \sigma^2 x + (2 \sigma^2 \mu + \sigma^4) - (2 \sigma^2 \mu + \sigma^4) \right]\right\} \\ &= \exp\left\{ -\frac{1}{2\sigma^2} \left[ (x - \mu - \sigma^2)^2 - (2 \sigma^2 \mu + \sigma^4) \right]\right\} \\ &= \exp\left\{ -\frac{1}{2 \sigma^2} \left[ x - \mu - \sigma^2 \right]^2 \right\} \exp\left\{ \mu + \frac{1}{2} \sigma^2 \right\}. \end{aligned} \tag{7}

The right term does not depend on xx and can thus be pulled out the integral, giving us

I=exp{μ+12σ2}logk12πσ2exp{12σ2[xμσ2]2}dx.(8) I = \exp\left\{ \mu + \frac{1}{2} \sigma^2 \right\} \int_{\log k}^{\infty} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left\{ -\frac{1}{2 \sigma^2} \left[ x - \mu - \sigma^2 \right]^2 \right\} \text{d}x. \tag{8}

Now note that the term inside the new integral IrI_r is the density function for a normally distributed random variable,

GN(μ+σ2,σ2).(9) G \sim \mathcal{N}(\mu + \sigma^2, \sigma^2). \tag{9}

And the new integral is simply the survival function:

Ir=logk12πσ2exp{12σ2[xμσ2]2}dx=P(Glogk).(10) I_r = \int_{\log k}^{\infty} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp\left\{ -\frac{1}{2 \sigma^2} \left[ x - \mu - \sigma^2 \right]^2 \right\} \text{d}x = \mathbb{P}(G \geq \log k). \tag{10}

And we can rewrite the survival function in terms of the standard normal distribution’s CDF, which I’ll denote with Φ\Phi:

P(Glogk)=1P(Glogk)=1P(Gμσσlogkμσσ)=1Φ(logkμσσ)=Φ(σlogkμσ).(11) \begin{aligned} \mathbb{P}(G \geq \log k) &= 1 - \mathbb{P}(G \leq \log k) \\ &= 1 - \mathbb{P}\left(\frac{G - \mu}{\sigma} - \sigma \leq \frac{\log k - \mu}{\sigma} - \sigma\right) \\ &= 1 - \Phi\left(\frac{\log k - \mu}{\sigma} - \sigma\right) \\ &= \Phi\left(\sigma - \frac{\log k - \mu}{\sigma}\right). \end{aligned} \tag{11}

The last step holds because of the symmetry of the normal distribution. This gives us:

I=exp{μ+12σ2}Φ(σlogkμσ).(12) I = \exp\left\{ \mu + \frac{1}{2} \sigma^2 \right\} \Phi\left(\sigma - \frac{\log k - \mu}{\sigma}\right). \tag{12}

Finally, we need the normalizing factor in Equation 33. Here, that is one minus the CDF of YY, which is lognormally distributed. This CDF can be expressed in terms of the CDF of the normal distribution:

F(y)=Φ(logyμσ).(13) F(y) = \Phi\left(\frac{\log y - \mu}{\sigma} \right). \tag{13}

Putting this together with Equation 1212, we have:

E[Y^]=E[YY>k]=exp{μ+12σ2}Φ(σlogkμσ)1Φ(logkμσ).(14) \mathbb{E}[\hat{Y}] = \mathbb{E}[Y \mid Y \gt k] = \frac{\exp\left\{ \mu + \frac{1}{2} \sigma^2 \right\} \Phi\left(\sigma - \frac{\log k - \mu}{\sigma}\right)}{1 - \Phi\left(\frac{\log k - \mu}{\sigma}\right)}. \tag{14}

And we’re done! We can simplify this a bit if we see that the exponent is the mean of YY—see this appendix for a derivation—and if we let uu denote the zz-score (logkμ)/σ(\log k - \mu) / \sigma. This gives us

E[YY>k]=E[Y]Φ(σu)1Φ(u).(15) \mathbb{E}[Y \mid Y \gt k] = \frac{ \mathbb{E}[Y] \Phi\left(\sigma - u\right) }{ 1 - \Phi\left(u\right) }. \tag{15}

This result is not obvious to me, and I haven’t seen the equation before—hence why I’m deriving it. So to sanity check it, I used Monte Carlo sampling. In Figure 22, I generated a normalized histogram from one million samples from the truncated lognormal distribution and compared this to its PDF. I then compared the empirical mean to the mean using Equation 1515. This suggests that we’ve derived the mean of the truncated lognormal distribution correctly.

Figure 2. Monte Carlo estimation of the lognormal distribution and its mean, compared with the mean derived in Equation 1515. The parameters used were μ=0\mu = 0, σ=1\sigma = 1, and k=1.5k = 1.5.

Here is the Python code used to generate this figure.

from scipy.stats import norm, lognorm

N = norm(0, 1).cdf

def trunc_lognorm_pdf(xx, mu, sigma, k):
    ln_dist = lognorm(scale=np.exp(mu), s=sigma)
    numer = ln_dist.pdf(xx)
    numer[xx <= k] = 0
    denom = 1 - ln_dist.cdf(k)
    return numer / denom

def trunc_lognorm_mean(mu, sigma, k):
    z = (np.log(k) - mu) / sigma
    return np.exp(mu + 0.5 * sigma**2) * N(sigma - z) / (1 - N(z))

def trunc_lognorm_rvs(mu, sigma, k, size=1):
    out = np.empty(0)
    while out.size != size:
        xx = lognorm(scale=np.exp(mu), s=sigma).rvs(size=size)
        out = np.append(out, xx[xx > k])[:size]
    return out

See this appendix for details on SciPy’s parameterization of its lognormal implementation (lognorm).