|
| 1 | +> ## Important Information |
| 2 | +> This file has a lot of Latex and GitHub currently cannot render it on Markdown files. You can read all the math clearly as a [webpage](https://numericalalgorithmsgroup.github.io/NAGJavaExamples/nearest_correlation_matrices) or access this as a regular github [repository](https://github.com/numericalalgorithmsgroup/NAGJavaExamples/tree/main/nearest_correlation_matrices). |
| 3 | +> |
| 4 | +> The source of this example can be found [here](https://github.com/numericalalgorithmsgroup/NAGJavaExamples/blob/main/nearest_correlation_matrices/NcmNag.java). |
| 5 | +> |
| 6 | +> See the top directory of this repository for instructions to set up the [NAG Library for Java](https://github.com/numericalalgorithmsgroup/NAGJavaExamples). |
| 7 | +
|
| 8 | +# Fast Implied Volatilities using the NAG Library |
| 9 | + |
| 10 | +The Black-Scholes formula for the price of a European call option is |
| 11 | + |
| 12 | +$$P = 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),$$ |
| 13 | + |
| 14 | +where $$T$$ is the time to maturity, $$S_0$$ is the spot price of the underlying asset, $$K$$ is the strike price, $$r$$ is |
| 15 | +the interest rate and $$\sigma$$ is the volatility. A similar formula applies for European put options. |
| 16 | + |
| 17 | +An important problem in finance is to compute the implied volatility, $$σ$$, given values for $$T$$, $$K$$, $$S_0$$, |
| 18 | +$$r$$ and $$P$$. An explicit formula for $$\sigma$$ is not available. Furthermore, $$\sigma$$ cannot be directly measured from |
| 19 | +financial data. Instead, it must be computed using a numerical approximation. Typically, multiple values |
| 20 | +of the input data are provided, so the Black-Scholes formula must be solved many times. |
| 21 | + |
| 22 | +As shown in the figure below, the volatility surface (a three-dimensional plot of how the volatility varies |
| 23 | +according to the price and time to maturity) can be highly curved. This makes accurately computing |
| 24 | +the implied volatility a difficult problem. |
| 25 | + |
| 26 | +<div style="text-align: center;"> |
| 27 | + <img src="impvolsurf.png" width=500 /> |
| 28 | +</div> |
| 29 | + |
| 30 | +Before introducing our new NAG Library routine, let’s demonstrate how one might naively |
| 31 | +compute implied volatilities using a general purpose root finder. |
| 32 | + |
| 33 | +Let's generate some input data using a random number generator from the NAG Library: |
| 34 | + |
| 35 | +```java |
| 36 | + G05KG g05kg = new G05KG(); |
| 37 | + G05SQ g05sq = new G05SQ(); |
| 38 | + int ifail = 0; |
| 39 | + int lstate = 17; |
| 40 | + int[] state = new int[lstate]; |
| 41 | + int n = 10000; // This is the number of volatilities we will be computing |
| 42 | + |
| 43 | + double[] p = new double[n]; |
| 44 | + double[] k = new double[n]; |
| 45 | + double[] s0 = new double[n]; |
| 46 | + double[] t = new double[n]; |
| 47 | + double[] r = new double[n]; |
| 48 | + |
| 49 | + g05kg.eval(1, 0, state, lstate, ifail); |
| 50 | + g05sq.eval(n, 3.9, 5.8, state, p, ifail); |
| 51 | + g05sq.eval(n, 271.5, 272.0, state, k, ifail); |
| 52 | + g05sq.eval(n, 259.0, 271.0, state, s0, ifail); |
| 53 | + g05sq.eval(n, 0.016, 0.017, state, t, ifail); |
| 54 | + g05sq.eval(n, 0.017, 0.018, state, r, ifail); |
| 55 | +``` |
| 56 | + |
| 57 | +We have chosen the limits of the various uniform distributions above to ensure the input data takes |
| 58 | +sensible values. |
| 59 | + |
| 60 | +There are various standard root finding techniques that we could use to compute implied volatilities, |
| 61 | +a common example being bisection. The NAG Library routine ```C05AW```, is a general |
| 62 | +purpose root finder based on the secant method. It uses a *callback*, with data passed in via a communication object: |
| 63 | + |
| 64 | +```java |
| 65 | + public static class BlackScholes extends C05AW.Abstract_C05AW_F { |
| 66 | + public double eval() { |
| 67 | + double[] price = new double[1]; |
| 68 | + int ifail = 1; |
| 69 | + |
| 70 | + S30AA s30aa = new S30AA(); |
| 71 | + s30aa.eval("C", 1, 1, new double[]{this.getRUSER()[1]}, this.getRUSER()[2], new double[]{this.getRUSER()[3]}, this.getX(), this.getRUSER()[4], 0.0, price, 1, ifail); |
| 72 | + |
| 73 | + ifail = s30aa.getIFAIL(); |
| 74 | + |
| 75 | + if (ifail != 0) |
| 76 | + price[0] = 0.0; |
| 77 | + |
| 78 | + return price[0] - this.getRUSER()[0]; |
| 79 | + } |
| 80 | + } |
| 81 | +``` |
| 82 | + |
| 83 | + ```C05AW``` operates on scalars, so we need to call the routine |
| 84 | +once for every volatility we want to compute. We will time the computation and count how many errors are caught: |
| 85 | + |
| 86 | +```java |
| 87 | + int[] iuser = new int[5]; |
| 88 | + double[] ruser = new double[5]; |
| 89 | + int errorcount = 0; |
| 90 | + double sigma; |
| 91 | + |
| 92 | + C05AW c05aw = new C05AW(); |
| 93 | + BlackScholes blackScholes = new BlackScholes(); |
| 94 | + |
| 95 | + long tic = System.currentTimeMillis(); |
| 96 | + for (i = 0; i < n; i++) { |
| 97 | + //System.out.println("Info: i = " + i); |
| 98 | + ruser[0] = p[i]; |
| 99 | + ruser[1] = k[i]; |
| 100 | + ruser[2] = s0[i]; |
| 101 | + ruser[3] = t[i]; |
| 102 | + ruser[4] = r[i]; |
| 103 | + |
| 104 | + ifail = 1; |
| 105 | + c05aw.eval(0.15, 1.0e-14, 0.0, blackScholes, 500, iuser, ruser, ifail); |
| 106 | + |
| 107 | + sigma = c05aw.getX(); |
| 108 | + ifail = c05aw.getIFAIL(); |
| 109 | + |
| 110 | + if ((sigma < 0.0) || (ifail != 0)) { |
| 111 | + errorcount++; |
| 112 | + } |
| 113 | + } |
| 114 | + long toc = System.currentTimeMillis(); |
| 115 | + |
| 116 | + System.out.println("Using a general purpose root finder:"); |
| 117 | + System.out.printf(" Time taken: %.5f seconds\n", (toc-tic)/1000.0); |
| 118 | + System.out.printf(" There were %d failures\n", errorcount); |
| 119 | +``` |
| 120 | + |
| 121 | +Can a bespoke implied volatility routine do better? Our new routine at Mark 27.1 is called ```S30AC```. We call it as follows: |
| 122 | + |
| 123 | +```s30ac.eval("C", n, p, k, s0, t, r, sigma_arr, 2, ivalid, ifail);``` |
| 124 | + |
| 125 | +The return argument ```ivalid``` is an array recording any data points for which the volatility could not be computed. The argument ```mode``` allows us to select which algorithm to use – more on that in a moment, but |
| 126 | +for now we choose ```mode=2```. This selects the algorithm of Jäckel (2015), a very accurate method based |
| 127 | +on third order Householder iterations. |
| 128 | + |
| 129 | +Here is the call surrounded by some timing code: |
| 130 | + |
| 131 | +```java |
| 132 | + S30AC s30ac = new S30AC(); |
| 133 | + double[] sigma_arr = new double[n]; |
| 134 | + int[] ivalid = new int[n]; |
| 135 | + |
| 136 | + tic = System.currentTimeMillis(); |
| 137 | + ifail = 0; |
| 138 | + s30ac.eval("C", n, p, k, s0, t, r, sigma_arr, 2, ivalid, ifail); |
| 139 | + toc = System.currentTimeMillis(); |
| 140 | + |
| 141 | + System.out.println("S30AC with mode = 2 (Jäckel algorithm):"); |
| 142 | + System.out.printf(" Time taken: %.5f seconds\n", (toc-tic)/1000.0); |
| 143 | + System.out.printf(" There were %d failures\n", nonZeroLength(ivalid)); |
| 144 | +``` |
| 145 | + |
| 146 | +The new routine is several orders of magnitude faster than the root finder, with no failures reported. We could try |
| 147 | +tweaking the convergence parameters and iteration limits in ```C05AW```, and we could certainly |
| 148 | +improve the way data is passed through the callback, but we are unlikely to match the |
| 149 | +performance of ```S30AC```. |
| 150 | + |
| 151 | +Recently NAG embarked upon a collaboration with mathematicians at Queen Mary University of |
| 152 | +London, who have been developing an alternative algorithm for computing implied volatilities. The new |
| 153 | +algorithm (based on Glau et. al. (2018)) uses Chebyshev interpolation to remove branching and give |
| 154 | +increased SIMD performance. We access it by setting ```mode=1``` in the call to ```S30AC```: |
| 155 | + |
| 156 | +```java |
| 157 | + tic = System.currentTimeMillis(); |
| 158 | + ifail = 0; |
| 159 | + s30ac.eval("C", n, p, k, s0, t, r, sigma_arr, 1, ivalid, ifail); |
| 160 | + toc = System.currentTimeMillis(); |
| 161 | + |
| 162 | + System.out.println("S30AC with mode = 1 (Glau algorithm):"); |
| 163 | + System.out.printf(" Time taken: %.5f seconds\n", (toc-tic)/1000.0); |
| 164 | + System.out.printf(" There were %d failures\n", nonZeroLength(ivalid)); |
| 165 | +``` |
| 166 | + |
| 167 | +Depending on your system, you should find that, for similar accuracy, there is a modest speedup over the Jäckel algorithm. Our numerical experiments have shown that for very small arrays (containing fewer than 100 elements) the Jäckel algorithm actually |
| 168 | +outperforms that of Glau et.al., but for larger arrays the converse is true. As vector units continue |
| 169 | +to improve in the future, we expect the performance of the highly vectorizable Glau et.al. approach to |
| 170 | +improve similarly. |
| 171 | + |
| 172 | +So far, we have been computing implied volatilities with a relative accuracy as close as possible to |
| 173 | +double precision. However, in some applications implied volatilities are only required with a few decimal |
| 174 | +places of precision. One advantage of the Glau et.al. algorithm is that it can be run in a lower accuracy |
| 175 | +mode, aiming only for seven decimal places of accuracy. This is accessed by setting ```mode=0``` in the call |
| 176 | +to ```S30AC```. It roughly doubles the speed of the routine: |
| 177 | + |
| 178 | +```java |
| 179 | + tic = System.currentTimeMillis(); |
| 180 | + ifail = 0; |
| 181 | + s30ac.eval("C", n, p, k, s0, t, r, sigma_arr, 0, ivalid, ifail); |
| 182 | + toc = System.currentTimeMillis(); |
| 183 | + |
| 184 | + System.out.println("S30AC with mode = 0 (lower accuracy Glau algorithm):"); |
| 185 | + System.out.printf(" Time taken: %.5f seconds\n", (toc-tic)/1000.0); |
| 186 | + System.out.printf(" There were %d failures\n", nonZeroLength(ivalid)); |
| 187 | +``` |
| 188 | + |
| 189 | +The charts below summarize the results, using timings collected on an Intel Skylake machine. We can see that the Glau et.al. algorithm outperforms the Jäckel algorithm for large arrays but not for small arrays. Note that the general purpose root finder |
| 190 | +is omitted here as it is so much slower ```S30AC```. |
| 191 | + |
| 192 | +<div style="text-align: center;"> |
| 193 | + <img src="graphs.PNG" width=800 /> |
| 194 | +</div> |
| 195 | + |
| 196 | +In summary, NAG’s new state-of-the art algorithm can be run in three different modes, according to |
| 197 | +the length of the input arrays and the required accuracy. For more information, and to access the NAG |
| 198 | +Library, go to: https://www.nag.co.uk/content/nag-library. |
| 199 | + |
| 200 | +### References |
| 201 | + |
| 202 | +P. Jäckel (2015). Let’s be rational. *Wilmott* 2015, 40-53. |
| 203 | + |
| 204 | +K. Glau, P. Herold, D. B. Madan, C. Pötz (2019). The Chebyshev method for the implied volatility. |
| 205 | +*Journal of Computational Finance*, 23(3). |
0 commit comments