Problem
I am trying to simulate $\log E[e^X]$ on a computer, where $X$ is a random variable which I can simulate only numerically. Due to the fat-tailed properties of the exponential-transform, convergence is really bad. However, $X$ itself is sort-of-ish well-behaved, being visually close to normally distributed (with a bit of negative skewness).
If $X$ would be exactly normal, we could instead compute $\log E[e^X] = E[X] + Var[X]/2$, where $E[X]$ and $Var[X]$ have good numerical convergence properties in our case. Hence, I was wondering, if I could approximately evaluate $\log E[e^X]$ using a similar approach.
Unfortunately, just using the log-normal expected value turns out not to be a good enough approximation, due to the negative skewness of $X$. My question is, is there a good way to correct for the skewness in a similar approximation approach?
Previous approach
I tried using a Taylor-expansion based approach (https://en.wikipedia.org/wiki/Taylor_expansions_for_the_moments_of_functions_of_random_variables), approximating
\begin{align} \log E[e^X] & \approx \log \left\{ e^{\mu_X} + e^{\mu_X} \frac{Var[X]} {2} + e^{\mu_X} \frac{E[(X-\mu_X)^3]} {3!} + \dots \right\} \\ & = \mu_X + \log \left\{ 1 + \frac{Var[X]} {2} + \frac{E[(X-\mu_X)^3]} {3!} + \dots \right\}, \end{align} but it turns out that the approximation is really bad. E.g., for $X\sim N(0,\sigma^2)$, we have $\log E[e^X] = \sigma^2/2$, whereas the quadratic approximation would give $\log (1+\sigma^2/2)$, which is only good for $\sigma$ small. Adding higher order terms didn't help much either in a previous numerical exploration.