International Journal of Computer Mathematics ISSN: 0020-7160 (Print) 1029-0265 (Online) Journal homepage: http://www.tandfonline.com/loi/gcom20 BENCHOP – The BENCHmarking project in option pricing Lina von Sydow, Lars Josef Höök, Elisabeth Larsson, Erik Lindström, Slobodan Milovanović, Jonas Persson, Victor Shcherbakov, Yuri Shpolyanskiy, Samuel Sirén, Jari Toivanen, Johan Waldén, Magnus Wiktorsson, Jeremy Levesley, Juxi Li, Cornelis W. Oosterlee, Maria J. Ruijter, Alexander Toropov & Yangzhang Zhao To cite this article: Lina von Sydow, Lars Josef Höök, Elisabeth Larsson, Erik Lindström, Slobodan Milovanović, Jonas Persson, Victor Shcherbakov, Yuri Shpolyanskiy, Samuel Sirén, Jari Toivanen, Johan Waldén, Magnus Wiktorsson, Jeremy Levesley, Juxi Li, Cornelis W. Oosterlee, Maria J. Ruijter, Alexander Toropov & Yangzhang Zhao (2015) BENCHOP – The BENCHmarking project in option pricing, International Journal of Computer Mathematics, 92:12, 2361-2379, DOI: 10.1080/00207160.2015.1072172 To link to this article: http://dx.doi.org/10.1080/00207160.2015.1072172 Accepted author version posted online: 24 Submit your article to this journal Aug 2015. Published online: 21 Sep 2015. Article views: 473 View related articles View Crossmark data Citing articles: 2 View citing articles Full Terms & Conditions of access and use can be found at http://www.tandfonline.com/action/journalInformation?journalCode=gcom20 Download by: [Bibliotheek TU Delft] Date: 06 January 2016, At: 06:27 International Journal of Computer Mathematics, 2015 Vol. 92, No. 12, 2361–2379, http://dx.doi.org/10.1080/00207160.2015.1072172 † BENCHOP – The BENCHmarking project in option pricing a∗ a a b Lina von Sydow , Lars Josef Höök , Elisabeth Larsson , Erik Lindström , a c a d Slobodan Milovanovic´ , Jonas Persson , Victor Shcherbakov , Yuri Shpolyanskiy , c e,f g b h Samuel Sirén , Jari Toivanen , Johan Waldén , Magnus Wiktorsson , Jeremy Levesley , h i,j i d h Juxi Li , Cornelis W. Oosterlee , Maria J. Ruijter , Alexander Toropov and Yangzhang Zhao a b Department of Information Technology, Uppsala University, Uppsala, Sweden; Centre for Mathematical c Sciences, Lund University, Lund, Sweden; SUNGARD FRONT ARENA, Stockholm, Sweden d www.sungard.com/frontarena; ORC GROUP, Stockholm, Sweden, www.orc-group.com and ITMO e University, St Petersburg, Russia; Department of Mathematical Information Technology, University of f Jyväskylä, Jyväskylä, Finland; The Institute for Computational and Mathematical Engineering (ICME), g Stanford University, Stanford, CA, USA; Haas School of Business, University of California Berkeley, h i Berkeley, CA, USA; Department of Mathematics, University of Leicester, Leicester, UK; Centrum j Wiskunde & Informatica, Amsterdam, The Netherlands; Delft Institute of Applied Mathematics, Delft University of Technology, Delft, The Netherlands (Received 30 January 2015; revised version received 18 June 2015; accepted 23 June 2015 ) The aim of the BENCHOP project is to provide the finance community with a common suite of benchmark problems for option pricing. We provide a detailed description of the six benchmark problems together with methods to compute reference solutions. We have implemented fifteen different numerical methods for these problems, and compare their relative performance. All implementations are available on line and can be used for future development and comparisons. Keywords: option pricing; numerical methods; benchmark problem; Monte Carlo method; Fourier- method; finite difference method; radial basis function 2010 AMS Subject Classifications: 65-02; 91G60; 91G20 1. Introduction The research on numerical methods for option pricing problems has been extensive over the last decades and there is now a plethora of methods targeting various types of options. However, there is a lack of cross comparisons between methods and a similar lack of common benchmarks to evaluate new approaches. The aim of BENCHOP is to provide the finance community with a set of common benchmark problems that can be used both for comparisons between methods and for evaluation of new methods. Furthermore, in order to facilitate comparisons, MATLAB implementations of a wide range of existing methods for each benchmark problem will be made available through the BENCHOP web site www.it.uu.se/research/project/compfin/benchop. *Corresponding author. Email: 2362 L. von Sydow et al. We also aim for BENCHOP to serve as a takeoff for future development of methods in option pricing. We expect future papers in the field to use the BENCHOP codes and problems to evaluate performance. In this way, we can contribute to a more uniform and comparable evaluation of the relative strengths and weaknesses of proposed methods. The benchmark problems have been chosen in such a way as to be relevant both for practi- tioners and researchers. They should also be possible to implement with a reasonable effort. We have selected problems with respect to a number of features that may be numerically challenging. These are early exercise properties, barriers, discrete dividends, local volatility, stochastic volatil- ity, jump diffusion, and two underlying assets. We have also included evaluation of hedging parameters in one of the problems, as this adds additional difficulties. In this paper, we present the benchmark problems with sufficient detail so that other people can solve them in the future. We also provide analytical solutions where such are available or methods for computing accurate reference solutions otherwise. Each problem is solved using MATLAB implementations of a number of already existing numerical methods, and timing results are provided as well as error plots. For details of the methods, we refer to the origi- nal papers and additional notes at the BENCHOP web site. The codes are not fully optimized, and the numerical results should not be interpreted as competition scores. We rather see it as a synoptical exposition of the qualities of the different methods. In Section 2, we state and motivate the benchmark problems while the numerical methods are briefly presented in Section 3. Section 4 is dedicated to the presentation of the numerical results and finally in Section 5, we discuss the results. In Appendix 1, we present how the reference values are computed for the different problems and in Appendix 2, we discuss how the local volatility surface is computed for one of the problems. 2. Benchmark problems In this section, we state each of the six benchmark problems. In the mathematical formulations, we let S represent the actual (stochastic) asset price realization, whereas s is the asset price variable in the PDE formulation of the problem, t is the time (with t = 0 representing today), r is the risk free interest rate, σ is the volatility, W is a Wiener process, u is the option price as a function of s, K is the strike price, and T is the time of maturity. The payoff function φ(s) is the value of the option at time T. In Problem 4, V is the stochastic variance variable, and v is the variance value in the PDE-formulation. In practice, the asset value today, S0, is a known quantity, while the strike price K can take on different values. In the benchmark problem descriptions below we have chosen to fix all parameters, and then solve for different values of S0 to simulate pricing of options that are ‘in the money’, ‘at the money’ and ‘out of the money’. The initial values are not given in the problem descriptions, but for each table and figure in the numerical results section, the values that were used are listed. 2.1 Problem 1: The Black–Scholes–Merton model for one underlying asset The celebrated Black–Scholes–Merton [4,39] option pricing model, developed in the early 1970’s, is arguably the most successful quantitative model ever introduced in social sciences, even initiating the new field of Financial Engineering, which occupies thousands of researchers in financial institutions and universities across the world. A key property of the model is that by building on so-called no-arbitrage arguments, it allows the price of plain vanilla call and put options to be calculated using variables that are either Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 International Journal of Computer Mathematics 2363 directly observable or can be easily estimated. The model is still widely used as a benchmark, although more advanced models have been developed over the years to take into account real- world features of asset price dynamics, such as jumps and stochastic volatility (see below). The Black–Scholes–Merton model has the advantage that closed form solutions exist for prices, as well as for hedging parameters, for some types of options. It has therefore been exten- sively used to test numerical methods that are then applied to more advanced problems. The computation of the hedging parameters (Greeks) is included in this benchmark problem as they are of significant practical interest and can be expensive and/or difficult to compute for some numerical methods. Mathematical formulation SDE-setting: dS = rS dt + σS dW. (1) 2 ∂u 1 ∂ u ∂u 2 2 PDE-setting: + σ s + rs − ru = 0. (2) 2 ∂t 2 ∂s ∂s Deliverables The pricing problem should be solved for three types of options; (a) a European call option, (b) an American put option, and (c) a barrier option. For the European option also the most common 2 2 hedging parameters = ∂u/∂s, Ŵ = ∂ u/∂s and V = ∂u/∂σ should be computed. Parameter and problem specifications For this problem we have two sets of model parameters, representing less and more numerically challenging situations, respectively. Standard parameters: σ = 0.15, r = 0.03, T = 1.0, and K = 100. (3) Challenging parameters: σ = 0.01, r = 0.10, T = 0.25, and K = 100. (4) The three types of options are characterized by their exercise properties and payoff functions. (a) European call : φ(s) = max(s − K, 0). (b) American put : φ(s) = max(K − s, 0), u(s, t) ≥ φ(s), 0 ≤ t ≤ T. { max(s − K, 0), 0 ≤ s < B (c) Barrier call up-and-out : φ(s) = , B = 1.25K. 0, s ≥ B 2.2 Problem 2: The Black–Scholes–Merton model with discrete dividends A shortcoming of the classical Black-Scholes formula is that it is only valid if the underlying stock does not pay dividends, invalidating the approach for many stocks in practice. In some special cases, for example, when dividend yields are constant and paid continuously over time, closed form solutions can be derived for dividend paying stocks too, see [39]. Usually, however, numerical methods are needed to calculate the option’s value. In practice, dividends are paid at discrete points in time, and the size of the dividend payments depends on the performance of the firm. For example, a firm whose performance has been poor may be capital constrained and therefore choose not to make a dividend payment, as may a company that needs its capital for a new investment opportunity. Fairly advanced stochastic modeling may therefore be needed in practice to capture the dividend dynamics of a company. In numerical tests, it is common to abstract away from these issues and simply assume that the firm makes discrete proportional dividend payments (i.e. has a constant dividend yield). Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 2364 L. von Sydow et al. Mathematical formulation SDE-setting: dS = rS dt + σS dW − δ(t − τ)D S dt. (5) 2 ∂u 1 ∂ u ∂u 2 2 PDE-setting: + σ s + rs − ru = 0, (6) 2 ∂t 2 ∂s ∂s In the SDE case, the (single) dividend at time τ enters explicitly, whereas in the PDE case, it is implicitly taken into account by enforcing − + u(s, τ ) = u(s(1 − D), τ ). (7) Deliverables Prices should be computed for a) a European call option and b) an American call option. Parameter and problem specifications The dividend is defined by τ = 0.4 and D = 0.03. We use the standard parameters (3) except for the expiration time, set to T = 0.5, together with standard payoff function φ(s) = max(s − K, 0) and European and American exercise properties respectively, see Problem 1. 2.3 Problem 3: The Black–Scholes–Merton model with local volatility As mentioned earlier, the Black–Scholes–Merton model with a constant volatility does not repro- duce market prices very well in practice. One discrepancy is the so-called volatility smile (which after the October 1987 crash is known to have turned into a smirk). If the implied volatility— the volatility in the Black–Scholes–Merton model that is consistent with the observed option price—is calculated for several options with the same exercise date but different strike prices, all options should under the classical assumptions of Black–Scholes–Merton have the same implied volatility. Instead, when plotted against the different strike prices, the curve is usually that of a U-shaped smile (or an L-shaped smirk). As discussed in [10], an approach to address this discrepancy between model and data is to assume local volatility, that is, to allow the volatility of the underlying asset to depend instanta- neously on the stock price s, and time t, generating a whole volatility surface. It is shown in [10] how to reverse engineer such a volatility surface from observed option prices. Given a volatility surface, the general Black–Scholes–Merton no-arbitrage approach can be used to derive the option price, although closed form solutions will typically no longer exist. We provide two volatility surfaces with different properties in order to see how the numerical methods handle such variable volatility coefficients. Mathematical formulation The equations for the option price are identical to (1) and (2) except that here σ = σ(s, t). Deliverables The price for a European call option should be computed in each case. Parameter and problem specifications The first local volatility surface is given by an explicit function 2 (s/100 − 1.2) σI(s, t) = 0.15 + 0.15(0.5 + 2t) . (8) 2 (s/100) + 1.44 The second local volatility surface σII (s, t) is based on market data and does not have an explicit form. The local surface is computed from a parametrization of the implied volatility data. The exact steps in the computation and the specific parametrization are given in Appendix 2. In both cases, we use K = 100 and r = 0.03, but the expiration times are different, with T = 1 for σI , and T = 0.5 for σII . The payoff function for a European call option is as before φ(s) = max(s − K, 0). Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 International Journal of Computer Mathematics 2365 2.4 Problem 4: The Heston model for one underlying asset The local volatility model allows for perfect matching of prices of European-style options but, just like the Black–Scholes–Merton model, also has its weaknesses. It does not perform very well for path dependent options and, moreover, there is clear evidence that in practice the volatility of asset prices is in itself random, beyond what can be simply be described as a function of time and underlying strike price [9,32,45]. The Heston model [20] assumes that in addition to the risk-factor that drives the value of the underlying asset, there is a another risk-factor that determines the underlying’s instantaneous variance, V. The PDE formulation of the model is therefore two-dimensional. Note that in contrast to the previous models, the market in Heston’s model is incomplete, and therefore additional assumptions about the market price of volatility risk are needed to determine the option price. The specific assumptions in [20] leads to the model below. Mathematical formulation SDE-setting: √ dS = rS dt + VS dW1, √ (9) dV = κ(θ − V) dt + σ V dW2, where W1 and W2 have correlation ρ. PDE-setting: 2 2 2 ∂u 1 ∂ u ∂ u 1 ∂ u ∂u ∂u 2 2 + vs + ρσvs + σ v + rs + κ(θ − v) − ru = 0. (10) 2 2 ∂t 2 ∂s ∂s∂v 2 ∂v ∂s ∂v Deliverables The price for a European call option should be computed. Parameter and problem specifications The model parameters are here given by r = 0.03, κ = 2, θ = 0.0225, σ = 0.25, ρ = −0.5, K = 100, and T = 1. The payoff function for the European call option is φ(s, v) = max(s − K, 0). With these parameters, the Feller condition is satisfied. 2.5 Problem 5: The Merton jump diffusion model for one underlying asset The Merton model [40] addresses another difference between real world asset price dynamics and the (local volatility) Black–Scholes–Merton model. What was identified early on, is that stock prices occasionally experience dramatic movements over very short time periods, that is, they sometimes ‘jump’. Such jumps make return distributions heavier-tailed than for pure diffusion processes, also in line with empirical observations and, as in the Heston model, causes the market to be incomplete, necessitating additional assumptions to price the option. The assumption used in [40] is that the underlying stock price follows a jump-diffusion process, where there is no risk- premium associated with jump risk. Under these conditions, the option price can be computed from a Partial-Integro Differential Equation (PIDE). Mathematical formulation SDE-setting: dS = (r − λξ)S dt + σS dW + S dQ, (11) where Q is a compound Poisson process with intensity λ > 0 and jump ratios that are √ 2 2 −(log y−γ ) /2δ log-normally distributed as p(y) = 1/ 2πyδe [50]. Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 2366 L. von Sydow et al. PIDE-setting: ∫ 2 ∞ ∂u 1 ∂ u ∂u 2 2 + σ s + (r − λξ)s − (r + λ)u + λ u(sy, τ )p(y) dy = 0. (12) 2 ∂t 2 ∂s ∂s 0 Deliverables The price for a European call option should be computed. Parameter and problem specifications 2 γ +δ /2 The parameters to use are r = 0.03, λ = 0.4, γ = −0.5, δ = 0.4, ξ = e − 1, σ = 0.15, K = 100, and T = 1. The payoff function is φ(s) = max(s − K, 0). 2.6 Problem 6: The Black–Scholes–Merton model for two underlying assets As an example of an option with more than one underlying, we use a spread option, which for K = 0 is called a Margrabe option [38]. This classic rainbow option has a payoff function that depends on two underlying assets, so the option price dynamic therefore depends on two risk- factors (as long as the two stocks’ returns are not perfectly correlated). In contrast to the Heston and Merton models, the model is still within the class of complete market models, and the option price is therefore completely determined without further assumptions. The reason that the market is complete in this case, in contrast to the other multi risk-factor models we have introduced, is that two underlying risky assets may be used in the formation of a hedging portfolio, whereas only one such asset is available with the Heston and Merton Models. Mathematical formulation SDE-setting: dS1 = rS1 dt + σ1S1 dW1, (13) dS2 = rS2 dt + σ2S2 dW2, where W1 and W2 have correlation ρ. PDE-setting: 2 2 2 ∂u 1 ∂ u ∂ u 1 ∂ u ∂u ∂u 2 2 2 2 + σ 1 s1 2 + ρσ1σ2s1s2 + σ2 s2 2 + rs1 + rs2 − ru = 0. (14) ∂t 2 ∂s 1 ∂s1∂s2 2 ∂s2 ∂s1 ∂s2 Deliverables The price for a European spread call option should be computed. Parameter and problem specifications The model parameters to use are r = 0.03, σ1 = σ2 = 0.15, ρ = 0.5, K = 0, and T = 1. The payoff function for the European call spread option is φ(s1, s2) = max(s1 − s2 − K, 0). 3. Numerical methods In Table 1, we display all methods that we have used. We also provide references to the original papers describing the methods. More information about the particular implementations used here can be found at www.it.uu.se/research/project/compfin/benchop. 4. Numerical results For each benchmark problem, we have decided on three (or five) evaluation points si (or (si, sj)). Each method must be tuned such that it delivers a solution u(si) with a relative error less than Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 International Journal of Computer Mathematics 2367 Table 1. List of methods used with abbreviations, marker symbol used in figures, and references. Abbr. Symbol Method References Monte Carlo methods MC Monte Carlo with Euler-Maruyama in time [16] MC-S Monte Carlo with analytical solution/Euler- [2,16,17,36,42,50] Maruyama/quadratic scheme in time and stratified sampling QMC-S Quasi Monte Carlo with analytical solution/Euler- [2,16,17,21,36,42,44,50] Maruyama/quadratic scheme in time, stratified sampling, and precomputed quasi random numbers Fourier methods FFT Fourier method with FFTs [6,31,33] FGL Fourier method with Gauss-Laguerre quadrature [1 (Page 890),7,18,31], [34 (Section 2.1–2.2), 35] COS Fourier method based on Fourier cosine series and [11,12,52,53] the characteristic function. Finite difference methods FD Finite differences on uniform grids with Rannacher [23,51,59,64, Chapter 78] smoothed CN in time FD-NU Finite differences on quadratically refined grids with [22,49,51,55,56] Rannacher smoothed CN / IMEX-CNAB in time FD-AD Adaptive finite differences with discontinuous [19,22,26,46,47,62] Galerkin / BDF-2 in time Radial basis function methods RBF Global radial basis functions with non-uniform [19,48] nodes and BDF-2 in time RBF-FD Radial basis functions generated finite differences [13,14,19,41,58,63,65] with BDF-2 in time RBF-PUM Radial basis functions partition of unity method with [19,54,57] BDF-2 in time RBF-LSML Least-squares multi-level radial basis functions with [19,27,28] BDF-2 in time RBF-AD Adaptive RBFs with CN in time [8,24,43] RBF-MLT Multi-level radial basis functions treating time as a [25,30,60] spatial dimension −4 10 in these points. We have not put any restrictions on the error in the rest of the domain. Due to this freedom, some codes have been tuned to narrowly target these points, while others (sometimes automatically) are tuned to achieve an evenly distributed error. Then the codes are run (on the same computer system) and the execution times are recorded. Each code is run four times, and the execution time reported in Tables 2–5 is the average of the last three runs. This is because the first time a MATLAB script is executed in a session it takes a bit longer. In the tables, we also show the approximate number of correct digits p in the result. This quantity is computed as p = [− log 10 er], where er is the maximum relative error and [·] indicates rounding. The maximum relative error is computed as ∣ ∣ ref ∣ u(s i) − u (si) ∣ ∣ ∣ er = max . ∣ ref ∣ i u (s i) For the Monte Carlo methods where the error is not deterministic, the errors are averaged over −4 the different runs. In some cases, a method was not able to reach a relative error of 10 within −3 reasonable time (1 hour), but a lower target 10 was attainable. These results are marked with a * in Tables 2–5. The execution time we report in the tables for Monte Carlo methods is the Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 2368 L. von Sydow et al. −4 Table 2. Problem 1. Computational time to compute a solution u that has a relative error < 10 at t = 0 and s = 90, 100, 110 for the standard parameters and at s = 97, 98, 99 for the challenging parameters. Standard parameters Challenging parameters Method (a) European (b) American (c) Up-and-out (a) European (b) American (c) Up-and-out ∗ MC 5.7e + 02 (3) – – × – – MC-S 1.5e + 01 (4) × × 1.6e + 01 (4) 9.8e − 02 (16) × QMC-S 5.7e − 01 (5) × – 6.3e − 01 (6) 1.2e + 02 (16) – FFT 1.3e − 03 (6) – – 1.3e − 03 (7) – – FGL 3.5e − 03 (14) 7.6e − 01 (5) – 1.8e − 02 (13) 2.2e + 00 (16) – COS 1.8e − 04 (5) 2.7e − 02 (4) 1.8e − 02 (4) 2.6e − 04 (4) 2.3e − 01 (5) 2.5e − 04 (4) FD 1.8e − 02 (4) 7.6e − 02 (4) 2.5e − 02 (4) 5.1e + 01 (4) 1.1e − 02 (11) 5.1e + 01 (4) FD-NU 9.2e − 03 (4) 5.8e − 02 (4) 1.6e − 02 (4) 5.0e − 02 (5) 6.0e − 03 (11) 5.2e − 02 (5) FD-AD 9.7e − 03 (4) 4.3e − 02 (4) 9.6e − 03 (4) 1.0e + 01 (4) 9.1e − 03 (4) 2.0e + 00 (4) RBF 6.2e − 02 (4) 4.6e + 00 (4) 1.4e − 01 (4) 7.7e + 01 (4) – 6.6e + 01 (4) RBF-FD 2.9e − 01 (4) 1.3e + 00 (4) 2.8e − 01 (4) 3.3e + 01 (5) 9.6e − 01 (4) 5.1e + 00 (4) RBF-PUM 2.8e − 02 (4) 3.6e + 00 (4) 5.4e − 02 (4) 3.4e + 00 (5) 4.5e + 00 (4) 1.9e + 00 (4) RBF-LSML 4.2e − 02 (4) – 3.0e − 02 (4) 7.5e + 00 (5) – – RBF-AD 7.9e − 01 (4) 1.7e + 01 (5) 2.4e + 01 (5) 3.3e + 00 (4) 1.2e + 00 (7) 1.6e + 01 (4) ∗ ∗ RBF-MLT 1.6e + 01 (4) – 2.4e + 02 (4) 3.3e + 02 (4) – 1.9e + 03 (3) Note: The numbers within parentheses indicate the approximate number of correct digits in the result. A ‘ − ’ indicates not implemented, while ‘ × ’ means implemented, but not accurate. 2 ∂u ∂ u ∂u Table 3. Problem 1. Computational time to compute hedging parameters = ∂s , Ŵ = ∂s2 and V = ∂σ that have a −4 relative error < 10 at t = 0 and s = 90, 100, 110 for the standard parameters and at s = 97, 98, 99 for the challenging parameters. Standard parameters Challenging parameters 2 2 ∂u ∂ u ∂u ∂u ∂ u ∂u Method = Ŵ = V = = Ŵ = V = 2 2 ∂s ∂s ∂σ ∂s ∂s ∂σ ∗ ∗ MC-S 3.0e + 00 (5) 5.7e + 01 (4) 1.7e + 01 (4) 3.3e + 00 (5) × 1.8e + 01 (3) ∗ ∗ QMC-S 5.7e − 01 (5) 1.0e + 00 (4) 6.3e − 01 (5) 6.3e − 01 (6) × 6.7e − 01 (4) FFT 1.3e − 03 (6) 1.2e − 03 (5) 2.1e − 03 (5) 1.2e − 03 (7) 1.2e − 03 (5) 2.5e − 03 (6) FGL 3.2e − 03 (14) 2.9e − 03 (14) 2.9e − 03 (14) 1.8e − 02 (14) 1.8e − 02 (11) 2.2e − 02 (11) COS 1.8e − 04 (4) 2.4e − 04 (5) 2.6e − 04 (5) 2.8e − 04 (4) 4.4e − 04 (5) 4.6e − 04 (5) FD 1.9e − 02 (5) 1.4e − 02 (5) 2.8e − 02 (4) 4.9e + 01 (5) 2.5e + 02 (4) 5.0e + 02 (5) FD-NU 7.8e − 03 (4) 8.9e − 03 (4) 1.8e − 02 (4) 6.5e − 02 (4) 5.0e − 01 (4) 1.6e + 00 (4) FD-AD 1.0e − 02 (4) 9.9e − 03 (4) 2.4e − 02 (4) 1.0e + 01 (4) 4.9e + 01 (4) 9.0e + 01 (4) RBF 6.6e − 02 (4) 6.9e − 02 (4) 7.9e − 02 (4) 9.3e + 01 (4) 3.1e + 02 (5) 7.0e + 02 (4) RBF-FD 3.0e − 01 (4) 2.5e − 01 (4) 5.6e − 01 (4) 3.4e + 01 (4) 3.7e + 01 (4) 1.0e + 02 (4) RBF-PUM 2.7e − 02 (4) 3.2e − 02 (4) 1.0e − 01 (4) 3.1e + 00 (5) 6.2e + 00 (4) 1.1e + 01 (4) RBF-LSML 1.3e − 01 (4) 2.7e − 01 (4) 7.1e − 02 (5) – – – RBF-AD 3.4e + 00 (4) 3.4e + 01 (4) 5.0e + 01 (4) 3.6e + 00 (5) 1.8e + 01 (4) 2.1e + 01 (4) ∗ ∗ RBF-MLT 1.8e + 01 (4) 1.8e + 01 (4) 3.2e + 01 (4) 4.2e + 02 (4) 3.4e + 02 (3) × Note: The numbers within parentheses indicate the approximate number of correct digits in the result. A ‘ − ’ indicates not implemented, while ‘ × ’ means implemented, but not accurate. time to compute the result for one evaluation point, whereas for other methods it is the time to compute the result for all evaluation points. In order to see how the errors behave away from the evaluation points, solutions are also plotted for a range of values in Figures 1–8. The figures show the absolute errors evaluated at the integer values between s = 60 and s = 160. No figure is shown for the second local volatility case in Problem 3, because there the local volatility result is only valid for a particular S0, not Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016 International Journal of Computer Mathematics 2369 −4 Table 4. Problems 2 and 3. Computational time to compute a solution u that has a relative error < 10 at t = 0 and s = 90, 100, 110. Discrete dividends Local volatility Method European call American call Smooth Implied ∗ MC-S 6.0e + 01 (3) – × × QMC-S 1.6e + 00 (4) – – – FFT 1.5e − 03 (7) – – – FGL 3.6e − 03 (14) 8.1e − 02 (6) – – COS 8.8e − 04 (5) 2.4e − 03 (4) 1.7e − 02 (4) – FD 2.0e − 02 (4) 1.5e − 02 (4) 2.2e − 02 (4) 1.2e + 00 (4) FD-NU 1.6e − 02 (4) 1.5e − 02 (4) 3.7e − 02 (4) 8.9e − 01 (4) FD-AD 2.1e − 02 (4) 2.3e − 02 (5) 3.6e − 02 (4) 8.8e − 01 (4) RBF 2.3e − 01 (4) 1.1e − 01 (4) 5.5e − 02 (4) 2.4e + 00 (4) RBF-FD 4.2e − 01 (4) 2.7e + 00 (4) 2.2e + 01 (4) 1.1e + 02 (4) RBF-PUM 3.3e − 02 (4) 3.0e − 02 (4) 1.4e − 01 (4) 9.0e − 01 (4) RBF-LSML 5.0e − 01 (4) 1.2e + 00 (4) 1.7e − 01 (4) 4.0e + 00 (4) Note: The numbers within parentheses indicate the approximate number of correct digits in the result. A ‘ − ’ indicates not implemented, while ‘ × ’ means implemented, but not accurate. −4 Table 5. Problems 4, 5, and 6. Computational time to compute a solution u that has a relative error < 10 at t = 0 −4 and s = 90, 100, 110 for the Heston and Merton models, and to compute a solution u that has a relative error < 10 at t = 0 and (s1, s2) = (100, 90), (100, 100), (100, 110), (90, 100), (110, 100) for the spread option. Method Heston Merton Spread ∗ MC-S × 1.6e + 01 (4) 2.9e + 01 (4) QMC-S – 3.4e − 01 (4) 1.8e + 00 (5) FFT 3.3e − 03 (5) 2.1e − 03 (5) – FGL 3.8e − 03 (14) 3.1e − 03 (13) 3.1e − 03 (14) COS 3.4e − 04 (4) 2.2e − 04 (4) 1.5e − 03 (4) FD – – – FD-NU 4.3e + 00 (4) 1.4e − 01 (4) 7.4e + 01 (4) FD-AD – – 4.7e + 01 (4) RBF 1.9e + 01 (4) – 7.7e + 01 (4) RBF-FD – – 2.2e + 03 (4) RBF-PUM 4.3e + 00 (4) – 1.3e + 01 (5) Note: The numbers within parentheses indicate the approximate number of correct digits in the result. A ‘ − ’ indicates not implemented, while ‘ × ’ means implemented, but not accurate. over a range. The vertical axis range in the figures is adjusted to the values that are plotted, but −20 the lower limit is not allowed to be lower than 10 . Errors falling below that value are not visible in the figures. The experiments have been performed on the Tintin cluster at Uppsala Multidisciplinary Cen- ter for Advanced Computational Science (UPPMAX), Uppsala University. The cluster consists of 160 dual AMD Opteron 6220 (Bulldozer) nodes. All codes are implemented (serially) in MATLAB. The names of the respective codes for each problem are indicated by the boldfaced heading over each (group of) plot(s). This generic name is then combined with an acronym for the particular method as for example BSeuCallUI_RBF.m. 5. Discussion Monte Carlo methods. MC methods are easy to implement in any number of dimensions, but the √ slow convergence rate, O(1/ N) for standard MC, makes it computationally expensive to reach Downloaded by [Bibliotheek TU Delft] at 06:27 06 January 2016