Tuesday, April 08, 2014

5 Minutes of Xtend

There is a relatively new JVM based language, Xtend. Their homepage says "JAVA 10, TODAY!", so I thought I would give it a try, I was especially interested in operator overloading support, and the fact that it compiles to Java code, not Java byte code.

Unfortunately, after 5 minutes with it, and pasting some non Java code in an xtend file, Eclipse hangs forever, even on restart. After creating another workspace, just to trash the new workspace a similar way. This is quite incredible for a nearly 2 years old project, on eclipse.org.

Saturday, April 05, 2014

Building a more accurate basis point volatility formula

P. Jaeckel has defied the limits of accuracy with his latest Black-Scholes volatility solver, managing to also improve performance compared to his earlier solver "By Implication". Out of a silly exercise, I decided to try my hand for a more accurate Normal (or basis point) volatility solver.

In reality, the problem is much simpler in the Bachelier/Normal model. A very basic analysis of Bachelier formula shows that the problem can be reduced to a single variable, as Choi et al explain in their paper. So the problem is not really one of solving, but one of approximating (the inverse of) a function.

The first step to build that function is to actually have a highly accurate slow solver as reference. This is quite easy, I just started with Choi formula and used Halley's method to refine. In reality, Halley's method is already a bit overkill on this problem: it works impressively well, 1 iteration is enough to have an insane level of accuracy, only noticeable when one works in high precision arithmetic (for example 50 digits). For double precision, Newton's method would actually be enough - I initially thought that my Halley's implementation did not work as it produced the exact same output as Newton in double precision. Li proposes the use of the SOR method, which for this exercise, behaves very much like Newton's method.

I then followed the logic from Choi et al, but working directly with in-the-money call options instead of straddles. Straddles sound neat at first (hides that we work in-the-money), but it's actually useless for the algorithm. Choi et al. ignore half of the straddle range when they use their eta transform in the paper. One other change is the mapping itself, I found a better mapping for the call options (but not that far of Choi initial idea). Finally, because I am lazy, I did not go to the pain of finding a good rational fraction approximation along with the square root problem they describe, I just tried a Chebyshev polynomial.

Unfortunately, a single Chebyshev polynomial does not work well: even with a very large (1000) degree it's not very precise, so much that I thought that my transform was garbage. I had noticed by mistake, that on another part (negative) of the interval, the Chebyshev polynomial worked actually very well to approximate something related to the volatility of another option. Suddendly came to me the idea of, like Johnson does in his Faddeeva package, using N Chebyshev polynomials on N small intervals. This is like the big heavy hammer for which everything looks like nails, but it's actually very fast to evaluate as the degree of each polynomial can then be low (7), plus a table lookup (could be coded as switch statements if one really cares about such details). The slowest part is actually the call to the log function.

The final bit is the use of a Taylor approximation for my -u/log(1-u) transform as it is not all that accurate in double precision when u is near 0. And that produces the following graph

It is interesting to note that "solving" the b.p. vol is 10x faster than solving the Black vol.

I wrote a small paper around all this where you'll find the details as well as Matlab code.

Wednesday, March 19, 2014

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.