A Bachelier Normal Implied Volatility Solver

Sun, Jun 7, 2026

Around 2014, I proposed a few Bachelier implied volatility “solvers”. The first one was inspired by Steven G. Johnson’s Faddeeva package for the complex error function: it used a piecewise Chebyshev polynomial representation to have a near machine accurate representation of the Bachelier implied volatility. Why did I put “solver” in quotes? Because the problem can be reduced to a 1D function representation (as originally shown by Choi, Kim and Kwak in 2009). Then, I was encouraged by Gary Kennedy to find a simpler rational function approximation, I used 3 or 4 rational functions to cover the full range of the Bachelier implied volatility function.

Recently, I decided to revisit a little bit the idea, using simpler transformations for some zones to make the solver even faster (and slightly more accurate). Compared to Peter Jäckel’s Implied Normal Volatility solver, it is more than twice as fast. This is because there is a correction step involving the normal density function in P.J.’s solver. Interestingly, testing out initially Quantlib’s implementation gave me the impression it wasn’t so accurate. There may be something wrong there as a direct port of P.J. algo to Rust or Java was very accurate.

The paper on the new LFK2026 Normal implied volatlity solver is available at https://arxiv.org/abs/2605.18343 and minimal Java code here. I also give minimal Rust code in the below update (the performance increase is even more significant with Rust).

Update June 16 I have tried two new things:

1. Use of Fused-Multply-Add instructions (FMA) for the Horner evaluation. This speeds up the calculation even more: the new approximations are nearly 3x faster than Jäckel. Another advantage is reduced error overall. This is not so visible in terms of maximum error, the effect is more obvious if we look at the RMSE.

   

2. Use of corrected FMA Horner scheme. This allows to gain a few extra (typically two) machine epsilons in accuracy. The cost is decrease in performance (nearly 2x computational cost).

   

Overall it is quite challenging to increase accuracy further. While it is an interesting problem to look at, it is not necessarily so practical: is this gain of one to three machine epsilons really worth the cost?

Rust code for the LFK solvers and reproduction of the results presented in the paper (instructions in the README.md).