Fast and Accurate Implied Volatility Solver

The calibration of a stochastic volatility model or a volatility surface parameterization (like SVI) involves minimizing the model options volatilities against market options volatilities. Often, the model computes an option price, not an implied volatility. It is therefore useful to have a fast way to invert that option price to get back the implied volatility that corresponds to it. Furthermore during the calibration procedure, the model option price can vary widely: it is convenient to have a robust implied volatility solver.

Another more basic use of implied volatility solvers, is for the computation of Black-Scholes greeks for a given market option price.

A few years ago, P. Jaeckel produced a much better algorithm than a simple Newton or Brent solver, in his paper "By Implication". There is also a much simpler algorithm from Li, based on SOR and a good initial guess for SOR, which I found to be actually quite robust and fast. Now P. Jaeckel has a newer algorithm, faster and more accurate, close to double accuracy.

I have tested those on 1 million options of random volatility (between 0.001 and 3), random strikes (N deviations with a cap at 1M) for a few expiries.

In terms of performance, the results are independent of expiries, but in terms of accuracy, the new Jaeckel algorithm is particularly more accurate for the long-term options (5y and 15y).

Algorithm   Expiry  Vol RMSE  Price RMSE  Time
Jaeckel2014 5y      3.8E-16   2.1E-16     1.8s
Jaeckel2006 5y      1.3E-10   1.1E-10     3.0s
Jaeckel2014 15y     2.0E-12   2.4E-16     1.8s
Jaeckel2006 15y     1.5E-7    8.0E-11     2.7s

The maximum error from Jaeckel2006 is around 1e-8, while the one from Jaeckel2014 is 2e-15 (for very large unrealistic strike)

As a comparison, the simpler Li SOR-TS algorithm follows the given price accuracy independently of the expiry; I have tested with 1E-12. The error in implied volatility will be slightly higher: different close vols can give the same price due to the maximum achievable accuracy of the Black-Scholes formula with double numbers, even with a good cumulative normal distribution implementation. Its performance is however dependent on the number of deviations considered: closer to ATM means faster for Li algorithm.

Algorithm   Expiry  Deviation  Vol RMSE  Price RMSE  Time
Jaeckel2014 1y      5          4.2E-16   2.0E-17     1.8s
LiSORTS     1y      5          8.5E-9    1.9E-13     2.1s
Jaeckel2014 1y      3          3.1E-16   5.9E-17     1.7s
LiSORTS     1y      3          4.3E-12   1.6E-13     1.3s

Actually, I have cheated in my Li SOR-TS implementation: I have reused the idea from P. Jaeckel to compute the Black-Scholes price with erfc_x (unscaled erfc) instead of erfc. This simple change divides the number of exp calls by 2. Without this trick, for 5 deviations, SOR-TS took 3.6s (almost twice the time).

I would not be surprised if this was the main performance improvement between the two Jaeckels.