# Log normal simulation.

I want to calculate numerically the expectation of a lognormal random variable $Y=e^X$, where $X$ is normally distributed with mean $m$ and variance $V$.

The expectation is known as $e^{m+\frac{1}{2}V}$. When it comes to simulation , we can generate $N$ random numbers $\{Y_{k}\}_{k=1}^{N}$ centered normally distributed, and calculate : $\frac{1}{N}\sum_{k=1}^{N}{e^{m+\sqrt{V}Y_k}}$.

When $m$ and $V$ are relatively small, we can replicate the expected value. When $m$ and $V$ are very high, we get ridiculously high number.

What is the best method to reduce that kind of numerical errors?

Thanks.

• You'd need to use a different sampling scheme, because in that case the variance becomes high. The idea is called importance sampling; the idea is that in this case you are computing $\int_{-\infty}^\infty e^x f(x) dx$, and you want to rewrite this as $\int_{-\infty}^\infty e^x \frac{f(x)}{g(x)} g(x) dx$, in such a way that if $Z$ has pdf $g$ then $\frac{e^Z f(Z)}{g(Z)}$ has a small variance. – Ian May 26 '16 at 17:13
• What is a good choice of g? – Canardini May 26 '16 at 17:17
• The theoretically optimal choice of $g$ would be proportional to $e^x f(x)$. In this case $\frac{e^Z f(Z)}{g(Z)}$ would be a constant r.v. However, this cannot work because the normalization constant for this pdf is exactly the quantity you want to compute. But it hints at what you want to do: you want the behavior to be like $e^x f(x)$. So for example you could find the maximum of $e^x f(x)$ and take $g$ to be a Gaussian which is peaked at that maximum instead of at $m$. – Ian May 26 '16 at 17:21
• How about resampling the lognormal numbers to match first and second moments ? – Canardini May 26 '16 at 17:26
• I'm not following; are you suggesting an approach for computing the mean, or asking a related question about simultaneously computing two moments? – Ian May 26 '16 at 17:27

## 1 Answer

Of course, it depends on what high means when you say $\ldots$

"When $m$ and $V$ are very high, we get ridiculously high number."

Nonetheless, I could not replicate your concerns. I did two simulations ($N = 10\,000$ each) using MATLAB, one with $m$ varying from $-70$ to $70$, and another with $v$ varying from $0.5$ to $30$. The two figures show the expected values, computed from the simulated random variables, versus $\mathrm{e}^{m + \frac{1}{2}v}$. The axes are log-scaled to better see the deviations from the diagonal.

Left figure: If $v$ is fixed and $m$ varies from $-70$ to $70$, expected values are on the diagonal.

Right figure: If $m$ is fixed and $v$ is varied from $0.5$ to $30$, expected values scatter around the diagonal, and it does so the more, the larger $v$ is.

• In my simulation, I have m=-320.31 and V=639.02. I am expecting 0.449328, and I am getting 4.68e-103 for 10k paths. My simulations give me a mean of -320.74 and variance 633.02, which is pretty good. Switch to the exponential kills the process – Canardini May 26 '16 at 18:58
• Wow! If one bears in mind that $\exp(\sqrt{v})$ is the geometric standard deviation, I wonder what quantity in the universe could fluctuate by a factor of $\exp(\sqrt{640}) \approx 10^{11}$. – Björn Friedrich May 26 '16 at 19:27
• The largest reported value that I have seen in the literature was $v = 12.26$ (i.e. $\exp(\sqrt{v}) = 33.15$) in a measurement of species abundance in birds. – Björn Friedrich May 26 '16 at 19:35
• Let assume a process $r_t=x_t+\phi_t$ where $dx_t=-0.01x_t dt+0.2dW_t$. We have $e^{-\int_{0}^{T}{\phi ds}}=e^{-0.02*T}*e^{-0.5*V(T)}$ , with $V(T)=\frac{0.2*0.2}{0.01*0.01}(T-2\frac{1-e^{T0.01}}{0.01}+\frac{1-e^{T*0.02}}{0.02})$. I am trying to get $E(e^{-\int_{0}^{40}{r(s) ds}})$ – Canardini May 26 '16 at 19:46