# construct volatility smile based on historic observations

So I calculated historic volatility/skewness/kurtosis for a commodity. I now would like to construct a volatility smile that reflects this historically realized distribution. I tried using some cornish-fisher like expansions but it has only a limited valid range and it seems in general a bit more complicated than necessary. Any ideas here how to do this?

IMHO the simplest way would be to: (1) fit a probability distribution to the $T$-period returns you've historically observed. This can be done by moment-matching the sample variance/skewness/kurtosis statistics you've already computed, or using kernel density estimation (2) compute European option prices by numerically integrating the $T$-period returns pdf (transforming it to prices if it's easier for you) (3) inverting the obtained prices through the BS formula to extract your BS implied volatility smile.
Be careful that the information you're using is provided under the real-world measure $\mathbb{P}$ and not in the risk-neutral world $\mathbb{Q}$. So when computing the option prices you should use a stochastic discount factor if you want to be rigorous. See discussions here for instance.