The Black-Scholes formula for the price of a European call option is

\[ C = S_0 \Phi \left( \frac {\ln \left( \frac {S_0}{K} \right) + \left[ r + \frac {\sigma^2}{2} \right] T}{\sigma \sqrt {T}} \right) – Ke^{rT} \Phi \left( \frac {\ln \left( \frac {S_0}{K} \right) + \left[ r – \frac {\sigma^2}{2} \right] T}{\sigma \sqrt {T}} \right) \]

where \(T\) is the time to maturity of the contract, \(S_0\) is the spot price of the underlying asset, \(K\) is the strike price of exercising the option, \(r\) is the interest rate and \(\sigma\) is the volatility. An important problem in finance is to compute the *implied volatility*, \(\sigma\), given values for \(T\), \(K\), \(S_0\), \(r\) and \(C\). An explicit formula for \(\sigma\) is not available. Furthermore, \(\sigma\) cannot be directly measured from financial data. Instead, it must be computed using a numerical approximation.

As shown in the figure, the *volatility surface* (a three-dimensional plot of how the volatility varies according to the price and time to maturity) can be highly curved. This makes accurately computing the volatility a difficult problem.

Mark 27.1 of the NAG Library contains a new routine, s30acf, for computing the implied volatility of a European option contract for arrays of input data.

This routine gives the user a choice of two algorithms. The first is the method of Jäckel (2015), which uses a third order Householder method to achieve close to machine accuracy for all but the most extreme inputs. This method is fast for short vectors of input data.

The second algorithm is based on that of Glau et al. (2018), with additional performance enhancements developed in a collaboration between NAG and mathematicians at Queen Mary University of London. This method uses Chebyshev interpolation and is designed for long vectors of input data, where vector instructions can be exploited. For applications in which accuracy to machine precision is not required, the algorithm can also be instructed to aim for accuracy to roughly single precision (approximately seven decimal places), giving even further performance improvements.