BENCHOP – SLV: the BENCHmarking project in Option Pricing – Stochastic and Local Volatility problems

ABSTRACT

In the recent project BENCHOP – the BENCHmarking project in Option Pricing we found that Stochastic and Local Volatility problems were particularly challenging. Here we continue the effort by introducing a set of benchmark problems for this type of problems. Eight different methods targeted for the Stochastic Differential Equation (SDE) formulation and the Partial Differential Equation (PDE) formulation of the problem, as well as Fourier methods making use of the characteristic function, were implemented to solve these problems. Comparisons are made with respect to time to reach a certain error level in the computed solution for the different methods. The implemented Fourier method was superior to all others for the two problems where it was implemented. Generally, methods targeting the PDE formulation of the problem outperformed the methods for the SDE formulation. Among the methods for the PDE formulation the ADI method stood out as the best performing one.

Keywords:

• Option pricing
• stochastic and local volatility
• numerical methods
• benchmark problem
• stochastic differential equation
• partial differential equation
• characteristic function

2010 AMS SUBJECT CLASSIFICATIONS:

• 65-02
• 91G60
• 91G20

1. Introduction

The aim of the original BENCHOP project described in [34] was to provide researchers and practitioners in finance with a set of common benchmark problems for comparisons between methods and for evaluation of new methods. A wide range of existing numerical methods for each benchmark problem was implemented and the relative performance of the methods was compared. The problems were selected with respect to features that may be numerically challenging such as early exercise properties, barriers, discrete dividends, local volatility, stochastic volatility, jump diffusion, and two underlying assets.

In [34] we found that stochastic and local volatility problems were particularly challenging. Hence, we here continue the effort to provide benchmark problems and performance comparisons between numerical methods at the research frontier for this type of problems. Stochastic and local volatility models provide an alternative approach to jump-diffusion models, and have proven effective in matching the rich asset price structure observed in derivatives markets. As is the case for jump-diffusion models, the market is typically incomplete in stochastic volatility models (in absence of assets that are directly traded on volatility risk), so additional assumptions are needed to uniquely determine asset prices. In [19], for example, it is assumed that the volatility risk-premium is proportional to the volatility level.

We present the benchmark problems with sufficient detail so that others can solve them in the future. We also provide analytical solutions where such are available. Each problem is solved using MATLAB implementations of a number of numerical methods, and error plots as a function of computational time are provided. For details of the methods, we refer to the original papers, listed in Table 1. The codes are not fully optimized, and the numerical results should not be interpreted as competition scores. We rather see it as a general indication of their performance in terms of obtained accuracy versus computational time that could be expected in the computing environment under consideration. In order to facilitate future comparisons, MATLAB p-codes of the here presented methods for each benchmark problem will be made available through the BENCHOP web site http://www.it.uu.se/research/scientific_computing/project/compfin/benchop.

Table 1. List of methods used with abbreviations, marker symbol used in figures, references, and whether the MATLAB-implementation makes use of parallellism.

Abbreviations	Symbol	Method	References	Parallellism
Methods for SDE formulation
MCA		Monte Carlo simulation with control and antithetic variables	[4,13,30]	Yes
MLMC		Multi-level Monte Carlo	[12]	No
SGBM		Stochastic grid bundling method	[10,22]	Yes
mSABR		The multi-step Monte Carlo simulation of the SABR model	[24]	Yes
Fourier method
FGL		Fourier method with Gauss–Laguerre quadrature	[1,5,6,14,23,25]	No
Methods for PDE formulation
ADI		Modified Craig–Sneyd alternating direction implicit finite difference method ($\theta =1/3$)	[17,20,21,36]	No
RBFFD		Radial basis function generated finite differences	[3,11,28,29]	Yes
RBFPUM		Radial basis function partition of unity methods	[28,31–33]	No

The paper is outlined as follows. In Section 2 we state 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. Finally, in Appendix 2 we specify the contributions from the authors of the paper.

2. Benchmark problems

2.1. Problem 1 – SABR stochastic-local volatility model

The Stochastic Alpha Beta Rho (SABR) model [18] is an established Stochastic Differential Equation (SDE) system which is often used for interest rates and FX modelling in practice. A key feature of the model is that it matches the observed dynamic behaviour of the volatility smile, namely that when the price of the underlying decreases, the volatility smile shifts to lower prices, and vice versa. The SABR model is based on a parametric local volatility component in terms of a model parameter, β. The formal definition of the SABR model reads $\begin{array}{rl}\mathrm{d}{F}_t& ={\sigma }_t{F}_t^{\beta }\mathrm{d}{W}_t^F,\\ \mathrm{d}{\sigma }_t& =\alpha {\sigma }_t\mathrm{d}{W}_t^{\sigma },\end{array}$ where $F_t=S_t\exp \left(r\right(T-t\left)\right)$ denotes the forward value of the underlying asset $S_t$, with r the interest rate, $S_0$ the spot price and T the contract's final time. The quantity ${\sigma }_t$ denotes the stochastic volatility, $W_t^F$ and $W_t^{\sigma }$ are two correlated Brownian motions with constant correlation coefficient ρ (i.e. $W_t^F{W}_t^{\sigma }=\rho t$). The free model parameters are $\alpha >0$ (the volatility of the volatility), $0\le \beta \le 1$ (the elasticity) and ρ (the correlation coefficient). The corresponding Partial Differential Equation (PDE) for the valuation of options is given by $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\frac{1}{2}{\sigma }^2{s}^{2\beta }\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{s}^2}+\rho \alpha s^{\beta }{\sigma }^2\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }\sigma \mathrm{\partial }s}+\frac{1}{2}{\alpha }^2{\sigma }^2\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{\sigma }^2}-ru=0,$ for s>0, $\sigma >0$ and $0\le t<T$.

Deliverables: The problem should be solved for a European call option with payoff $max(S-K,0)$ at t=T, with three strikes $\begin{array}{rl}K& =S_0\exp (0.1\times \sqrt{T}\times \delta ),\\ \delta & =-1.0,0.0,\mathrm{1.0.}\end{array}$

Parameter and problem specifications:

• Set I ([16]): T=2, r=0, $S_0=0.5$, ${\sigma }_0=0.5$, $\alpha =0.4$, $\beta =0.5$ and $\rho =0$.

• Set II ([7]): T=10, r=0, $S_0=0.07$, ${\sigma }_0=0.4$, $\alpha =0.8$, $\beta =0.5$ and $\rho =-0.6$.

2.2. Problem 2 – quadratic local stochastic volatility model

The Quadratic Local Stochastic Volatility (QLSV) model, can be viewed as a generalization of Heston's stochastic volatility model [19]. In the QLSV model, the square root term is multiplied by a quadratic function in the underlying asset. When the function is a constant, Heston's original model is obtained. The additional degrees of freedom provided by the quadratic function allows for improved volatility surface calibration.

We use the formulation in, e.g. [26] and define the following SDE: $\begin{array}{rl}\mathrm{d}{S}_t& =r{S}_t\mathrm{d}t+\sqrt{V_t}f\left(S_t\right)\mathrm{d}{W}_t^1,\\ \mathrm{d}{V}_t& =\kappa (\eta -V_t)\mathrm{d}t+\sigma \sqrt{V_t}\mathrm{d}{W}_t^2,\end{array}$ with $f\left(s\right)=\frac{1}{2}\alpha s^2+\beta s+\gamma$. The PDE for the valuation of options is then given by $\frac{\mathrm{\partial }u}{\mathrm{\partial }t}+\frac{1}{2}f(s{)}^{2}v\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{s}^2}+\rho \sigma f(s)v\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }s\mathrm{\partial }v}+\frac{1}{2}{\sigma }^{2}v\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{v}^2}+rs\frac{\mathrm{\partial }u}{\mathrm{\partial }s}+\kappa (\eta -v)\frac{\mathrm{\partial }u}{\mathrm{\partial }v}-ru=0,$ for s>0, v>0 and $0\le t<T$.

We consider a case when the Feller condition, $2\kappa \eta \ge {\sigma }^2$ is not satisfied (as is often the case in practice). The Feller condition [9] ensures that the volatility process, $V_t$, remains strictly positive. If the condition is violated, the process may reach the boundary $V_t=0$. In this case, additional specification of the behaviour at the boundary is needed to ensure (weak) positivity, and the problem is also numerically more challenging (see [8,27]).

Deliverables: The problem should be solved for

• a European call option with payoff $max(S-K,0)$ and K=100,

• a Double-no-touch option paying 1 if $L<S_t<U$ (for all t) and 0 else with L=50, U=150,

and three spot values: $(S_0,V_0)=(S_0,0.114)$ for $S_0=75,100,125$.

Parameter and problem specifications: We consider

• a Heston model with $\alpha =0$, $\beta =1$, $\gamma =0$,

• a QLSV model with $\alpha =0.02$, $\beta =0$, $\gamma =0$.

For both models we use the parameter set (see [26]): $T=1,r=0,\kappa =2.58,\eta =0.043,\sigma =1,\rho =-\mathrm{0.36.}$

2.3. Problem 3 – Heston–Hull–White model

The Heston–Hull–White (HHW) model is a hybrid asset price model combining the Heston stochastic volatility model, and Hull–White stochastic interest rate model, see e.g. [15,17]. With such an approach the skew pattern for equity can be matched, while still allowing for stochastic interest rates. We define the HHW SDE: $\begin{array}{rl}\mathrm{d}{S}_t& =R_t{S}_t\mathrm{d}t+\sqrt{V_t}{S}_t\mathrm{d}{W}_t^1,\\ \mathrm{d}{V}_t& =\kappa (\eta -V_t)\mathrm{d}t+{\sigma }_1\sqrt{V_t}\mathrm{d}{W}_t^2,\\ \mathrm{d}{R}_t& =a\left(b\right(t)-R_t)\mathrm{d}t+{\sigma }_2\mathrm{d}{W}_t^3,\end{array}$ and the corresponding HHW PDE: $\begin{array}{rl}\frac{\mathrm{\partial }u}{\mathrm{\partial }t}& +\frac{1}{2}{s}^{2}v\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{s}^2}+\frac{1}{2}{\sigma }_1^{2}v\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{v}^2}+\frac{1}{2}{\sigma }_2^2\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }{r}^2}\\ & +{\rho }_{12}{\sigma }_{1}sv\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }s\mathrm{\partial }v}+{\rho }_{13}{\sigma }_{2}s\sqrt{v}\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }s\mathrm{\partial }r}+{\rho }_{23}{\sigma }_1{\sigma }_2\sqrt{v}\frac{{\mathrm{\partial }}^{2}u}{\mathrm{\partial }v\mathrm{\partial }r}\\ & +rs\frac{\mathrm{\partial }u}{\mathrm{\partial }s}+\kappa (\eta -v)\frac{\mathrm{\partial }u}{\mathrm{\partial }v}+a(b\left(t\right)-r)\frac{\mathrm{\partial }u}{\mathrm{\partial }r}-ru=0,\end{array}$ for s>0, v>0, and $0\le t<T$. Again, we will consider a case when the Feller condition is violated, see Section 2.2.

Deliverables: The problem should be solved for a European call option with payoff $max(S-K,0)$, K=100, and three spot values: $(S_0,V_0,R_0)=(S_0,0.04,0.10)$ for $S_0=75,100,125$.

Parameter and problem specifications: We use the parameter set (cf. [2,17]): $\begin{array}{rl}T& =10,\kappa =0.5,\eta =0.04,{\sigma }_1=1,{\sigma }_2=0.09,{\rho }_{12}=-0.9,{\rho }_{13}=0,\\ {\rho }_{23}& =0,a=0.08\mathrm{and}b\left(t\right)\equiv \mathrm{0.10.}\end{array}$

3. Numerical methods

In Table 1 we display all the methods that we have used. The methods are developed from the SDE or PDE formulation of the problem respectively, or make use of the characteristic function for the Fourier method. The methods are organized according to this in the table. To have a uniform presentation in Section 4 we use the same abbreviation and marker symbol for a given method in all figures in this section. Finally, Table 1 provides references to the original papers describing the methods and also indicates whether the implementation is employing parallel constructions provided in Matlab Parallel Computing Toolbox or not. From the references in the table it should be clear which flavour of a method is actually used, when there are several similar approaches available.

It should be noted that even if a method is not parallelized here, it may still be parallelizable with good scalability.

4. Numerical results

In this section, we present the results using the different methods presented in Section 3 for the problems defined in Section 2. For all problems, we display errors in the computed solution as a function of wall clock time. The implementations were made in MATLAB and encrypted into p-code. The experiments were carried out at the Rackham Cluster at Uppsala Multidisciplinary Center for Advanced Computational Science at Uppsala University http://www.uppmax.uu.se/resources/systems/the-rackham-cluster/. The Rackham Cluster consists of 334 nodes where each node has two 10-core Intel Xeon V4 CPUs and at least 128 GB RAM memory. The experiments were performed on a dedicated node, i.e. using up to 20 MATLAB-workers. For the methods based on the SDE formulation of the problems, i.e. a Monte Carlo based solution method, the code was run three times and the mean of the results was reported. Finally, note that the axes in the figures are different for the different problems.

4.1. Results for SABR model

In Figures 1 and 2 we display the error $\mathrm{\Delta }{u}_{max}=\underset{K}{max}|u(S_0,{\sigma }_0,0)-u_{\mathrm{ref}}(S_0,{\sigma }_0,0)|,$ for the set of parameters (including $S_0$ and ${\sigma }_0$) defined in Section 2.1. Here u denotes the computed solution and $u_{\mathrm{ref}}$ the reference solution, as a function of wall clock time for the SABR stochastic-local volatility model. Figure 1 shows the results for Parameter Set I, and Figure 2 the results for Parameter Set II. For both parameter sets, the ADI method is the most favourable one to use – for a given computational time the error is always smallest for this method. Among the Monte Carlo based methods, mSABR is the least favourable one to use for Parameter Set I while it is quite favourable to use for Parameter Set II if the required error in the solution is not too small (slightly less than ${10}^{-3}$). For Parameter Set I MCA and MLMC behave similarly while MCA performs better than MLMC for Parameter Set II. For the RBF methods, RBFPUM performs better than RBFFD for relatively large errors ($\sim {10}^{-3}$). However, especially for Parameter Set I, RBFFD can provide smaller errors than RBFPUM does. For the SABR model, FGL was not implemented.

Figure 1. Results for the European call option under the SABR model, Parameter Set I. The reference values for $K_i=S_0\exp (0.1\times \sqrt{T}\times {\delta }_i)$, ${\delta }_i=-1.0,0.0,1.0$ are given by 0.221383196830866, 0.193836689413803, and 0.166240814653231.

Figure 2. Results for the European call option under the SABR model, Parameter Set II. The reference values for $K_i=S_0\exp (0.1\times \sqrt{T}\times {\delta }_i)$, ${\delta }_i=-1.0,0.0,1.0$ are given by 0.052450313614407, 0.046585753491306, and 0.039291470612989.

4.2. Results for QLSV model

The error (1) $\mathrm{\Delta }{u}_{max}=\underset{S_0}{max}|u(S_0,V_0,0)-u_{\mathrm{ref}}(S_0,V_0,0)|,$(1) as a function of wall clock time for the European call option under the Heston model is presented in Figure 3, and the Double-no-touch option under the same model in Figure 4. Here we have used the set of parameters (including $S_0$ and $V_0$) defined in Section 2.2. For the European call option, FGL outperforms all other methods. Already in less than 0.1 s, the error is below ${10}^{-7}$ for this method. The second best method is ADI which is also the method that is performing best for the Double-no-touch option. When it comes to the two RBF methods, RBFFD is the method of choice for this problem. For the European call option, the error for RBFPUM stays slightly below ${10}^{-1}$ and for the Double-no-touch option it stays around ${10}^{-2}$. RBFFD on the other provides an error that is slightly larger than ${10}^{-4}$ in less than 100 s. The SGBM method is the best performing Monte Carlo type method for the Heston model. Both for the European option and the Double-no-touch option the time to reach a certain accuracy for SGBM is well below the time for MCA. Finally, for the European option, the performance of RBFFD and SGBM is similar, while RBFFD gives better results than SGBM for the Double-no-touch option.

Figure 3. Results for the European call option under the Heston model. The reference values for $S_0=75$, 100, and 125 are given by 0.908502728459621, 9.046650119220969, and 28.514786399298796.

Figure 4. Results for the Double-no-touch option under the Heston model. The reference values for $S_0=75$, 100, and 125 are given by 0.834539127387590, 0.899829293208866, and 0.668399975738358.

In Figures 5 and 6 the results for the European and Double-no-touch options under the QLSV model are presented. The Fourier method was not implemented for this model and the ADI method performed best for both types of options this time. To obtain low accuracy (${10}^{-1}$ and ${10}^{-2}$ for European and Double-no-touch options respectively) the RBFPUM is the second best method that reaches these accuracies in less than 1 s. To obtain higher accuracy (${10}^{-3}$ and ${10}^{-4}$ respectively) RBFFD is preferable between the two RBF methods. For the Double-no-touch option, MCA performs similarly as RBFFD while RBFFD outperforms MCA for the European call option.

Figure 5. Results for the European call option under the QLSV model. The reference values for $S_0=75$, 100, and 125 are given by 0.527472759419533, 8.902347915743665, and 29.159828965633729.

Figure 6. Results for the Double-no-touch option under the QLSV model. The reference values for $S_0=75$, 100, and 125 are given by 0.933800903110254, 0.914799140676374, and 0.592983062889906.

4.3. Results for HHW model

In Figure 7 we present the results from the Heston–Hull–White model which is the only three-dimensional model among our benchmark problems. The error (2) $\mathrm{\Delta }{u}_{max}=\underset{S_0}{max}|u(S_0,V_0,R_0,0)-u_{\mathrm{ref}}(S_0,V_0,R_0,0)|,$(2) as a function of wall clock time is presented for the set of parameters (including $S_0$, $V_0$, and $R_0$) defined in Section 2.3. Again, FGL outperforms all other methods and reaches $\mathrm{\Delta }{u}_{max}<{10}^{-7}$ in less than 0.1 s. ADI shows the same robust convergence behaviour as in the previous cases but this time MCA seems to perform equally well when considering accuracy versus computational time. It should be noted that the ADI code is a serial code while MCA is parallelized. Further, for MCA confidence intervals are not given. Finally, RBFPUM gives a fairly large error in this experiment.

Figure 7. Results for the European call option under the HHW model. The reference values for $S_0=75$, 100, and 125 are given by 35.437896876285350, 54.728065308229503, and 75.397596834993621.

5. Discussion

In this paper, we have defined a set of benchmarking problems in option pricing for stochastic and local volatility problems. Eight different numerical methods have been implemented and compared with respect to error versus computational time.

The Fourier method FGL was implemented only for the European call option with the Heston model and the Heston–Hull–White model. In both cases, this method was superior to all others and reached extremely high accuracy with an error less than ${10}^{-7}$ in less than 0.1 s. However, Fourier methods rely on the availability of the characteristic function of the underlying stochastic process or efficient approximations of the same.

For the problems when FGL was not implemented, the most efficient method was the ADI method. This was also together with RBFPUM and MCA the only method that was implemented for all problems. For the RBF methods, RBFPUM often reached a reasonable accuracy in a short time but had troubles with convergence with increasing computational time. RBFFD however, was in many cases the method that reached the smallest error after ADI (and FGL where applicable). The main challenge for the RBF methods was to apply appropriate boundary conditions in the volatility and interest rate dimensions. Some commonly used conditions introduced bias, while the option of leaving the boundary open in some cases lead to instability or inaccuracy.

For the methods based on the SDE formulation of the problems, it is difficult to draw any far-reaching conclusions. MLMC performed reasonably well for the SABR problems while mSABR gave reasonably small errors for Parameter Set II in a fairly short time. The mSABR method requires a smart application of a stochastic collocation-based method (see [16]) to make the algorithm affordable in terms of computational time. Once the so-called integrated variance and volatility processes are simulated, the SABR forward asset process has been simulated by an extended Log-Euler discretization scheme, called Log-Euler+. MCA also gave quite accurate results in a reasonable time, especially for the SABR problems. For the Heston model (the special case of QLSV), SGBM was the best performing Monte Carlo type method. For the 3D problem Heston–Hull–White, MCA and ADI showed comparable results in this setting. This indicates as expected that Monte Carlo methods probably will be more competitive in higher dimensions.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

The computations were performed on resources provided by SNIC through Uppsala Multidisciplinary Center for Advanced Computational Science (UPPMAX) under Project SNIC 2017/1-236.

Appendices

Appendix 1. The methods used for computing the reference values used in the comparisons

When there is no analytical formula available we have used the ADI method with a very small tolerance to compute the reference values. The estimated error in these values is computed as the absolute difference between two subsequent solutions in a refinement sequence of ADI solutions. As with all numerical methods, there could be some bias in this.

A.1. Methods for problem 1 – SABR

Case I: To our knowledge, no analytical formula allowing for highly accurate evaluation of the European call option under the SABR model was previously available in the literature. For the parameters we use in Case I, with uncorrelated Brownian motions $W_S$ for the stock and $W_{\sigma }$ for the volatility, a new price formula has been derived. The full details of the derivations are found in [35]. If we denote the European call option price at time t for the SABR model with strike K and maturity T by ${\pi }_C^E(t,S_t,T,K)$, we have that ${\pi }_C^E(0,S_0,T,K)=S_0(1-F_S(K\left)\right)-K(1-F_K(K\left)\right),$ where the probabilities $F_S\left(K\right)$ and $F_K\left(K\right)$ are expressed through the following integrals: $\begin{array}{rl}{F}_K\left(K\right)& =\frac{1}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}\mathrm{\Re }\left(\frac{q(\alpha ,{\sigma }_0,T,x)}{{\overline{z}}^{\prime }+i{\omega }^{\prime }}\right)\mathrm{d}{\omega }^{\prime }\mathrm{d}x,\\ F_S\left(K\right)& =\frac{1}{\pi }{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\int }_0^{\mathrm{\infty }}\mathrm{\Re }\left(\frac{q(\alpha ,{\sigma }_0,T,x)}{({\overline{z}}^{\prime }+i{\omega }^{\prime })(1+{\overline{z}}^{\prime }+i{\omega }^{\prime }{)}^2}\right)\mathrm{d}{\omega }^{\prime }\mathrm{d}x,\end{array}$ where $q(\alpha ,{\sigma }_0,T,x)=\frac{{\mathrm{e}}^{-\left(1/2{\alpha }^{2}T\right)\left(\mathrm{arcCosh}({\mathrm{e}}^{-x}\lambda {\alpha }^2/{\sigma }_0^2+\cosh (x){)}^2-x^2+(x+T{\alpha }^2/2{)}^2\right)}}{\sqrt{2\pi {\alpha }^{2}T}},$ and the variable ${\overline{z}}^{\prime }$, which is associated with the integration path of an inverse Laplace transform, should be chosen such that $0<{\overline{z}}^{\prime }<\frac{1}{2}\left(x_1+\sqrt{x_1^2+x_2}\right),$ where $\begin{array}{rl}{x}_1& =\frac{1}{K}{\left(\frac{({\mathrm{e}}^x-1){\sigma }_0}{2\alpha }\right)}^2+\frac{S_0}{K}-1,\\ x_2& =\frac{4}{K}{\left(\frac{({\mathrm{e}}^x-1){\sigma }_0}{2\alpha }\right)}^2.\end{array}$

We calculate these integrals using Gauss–Hermite quadrature points for x and Gauss–Laguerre quadrature points for ${\omega }^{\prime }.$ Using 2000 points in the x and ${\omega }^{\prime }$ direction, i.e. a total of 4 million points we get the values.

Table A1. SABR European call prices using the parameters from case I. Error $\sim 1\mathrm{e}-14$.

K	$0.5\exp (-0.1\sqrt{2})$	0.5	$0.5\exp \left(0.1\sqrt{2}\right)$
Call price	0.221383196830866	0.193836689413803	0.166240814653231

Case II: The solution was obtained by an ADI method since the analytical method only works in the zero correlation case. The ADI method employed $890\times 445$ spatial grid points and 445 time steps.

Table A2. SABR European call prices using the parameters from case II. Error $\sim 1\mathrm{e}-6$.

K	$0.07\exp (-0.1\sqrt{10})$	0.07	$0.07\exp \left(0.1\sqrt{10}\right)$
Call price	0.052450313614407	0.046585753491306	0.039291470612989

A.2. Methods for problem 2 – QSLV

Heston European call: The solution was obtained using a Fourier–Gauss–Laguerre implementation with 1000 quadrature points.

Table A3. Heston European call prices. Error $\sim 1\mathrm{e}-16$.

$S_0$	75	100	125
Call price	0.908502728459621	9.046650119220969	28.514786399298796

Heston Double-no-touch: The solution was obtained by an ADI method since analytical methods only work in the European case. The ADI method employed $1514\times 758$ spatial grid points and 757 time steps.

Table A4. Heston Double-no-touch prices with lower barrier L=50 and upper barrier U=150. Error $\sim 1.1\mathrm{e}-7$.

$S_0$	75	100	125
Call price	0.834539127387590	0.899829293208866	0.668399975738358

QSLV European call: The solution was obtained by an ADI method since analytical methods only work in the zero correlation case. The ADI method employed $1486\times 744$ spatial gridpoints and 743 time steps.

Table A5. Case quadratic stochastic-Local-Volatility European call prices. Error $\sim 2.4\mathrm{e}-6$.

$S_0$	75	100	125
Call price	0.527472759419533	8.902347915743665	29.159828965633729

QSLV Double-no-touch: The solution was obtained by an ADI method using $2120\times 1061$ spatial grid points and 1060 time steps.

Table A6. quadratic stochastic-Local-Volatility Double-no-touch prices with lower barrier L=50 and upper barrier U=150. Error $\sim 7.3\mathrm{e}-8$.

$S_0$	75	100	125
Call price	0.933800903110254	0.914799140676374	0.592983062889906

A.3. Methods for problem 3 – HHW

The solution was computed using a Fourier–Gauss–Laguerre implementation with 1000 quadrature points.

Table A7. Heston European call prices. Error $\sim 1\mathrm{e}-16$.

$S_0$	75	100	125
Call price	35.437896876285350	54.728065308229503	75.397596834993621

Appendix 2. The contribution from the authors

Apart from what is listed in Table 9 all authors took part in reading and correcting of the final version of the paper.

Table A8. Name and contribution of authors.

Name	Contributions
Lina von Sydow	Project leader
	Wrote Sections 1, 3, 4, 5, and Appendix 2
	Took part in writing the rest of the paper
	One of two behind code for RBFFD (SABR and QLSV models)
	Took part in meetings A, B, and C
Slobodan Milovanović	Performed all numerical experiments
	Created all plots in Section 4
	Took part in writing the paper
	One of two who wrote unit tests for contributed codes
	Took part in the computation of reference values
	One of two behind code for RBFFD (SABR and QLSV models)
	Took part in meetings A and B
Elisabeth Larsson	Deputy project leader
	Took part in the writing of Appendix 1
	Took part in writing the rest of the paper
	Took part in the computation of reference values
	One of two behind code for RBFPUM (all models)
	Took part in meetings A, B, and C
Karel In 't Hout	Lead the project of defining benchmark problems
	Took part in the writing of Section 2
	Took part in the computation of reference values
	One of two behind code for ADI (all models)
	Took part in meetings A, B, and C
Magnus Wiktorsson	Took part in the writing of Appendix 1
	Took part in the computation of reference values
	Behind code for FGL (Heston European call option and HHW model)
	Behind code for MCA (all models)
	Took part in meeting A
Cornelis W. Oosterlee	Took part in defining benchmark problems
	Took part in the writing of Section 2
	One of two behind code for mSABR (SABR model)
	One of two behind code for SGBM (QLSV model)
	Took part in meetings A, B, and C
Victor Shcherbakov	One of two who wrote unit tests for contributed codes
	One of two behind code for RBFPUM (all models)
	Took part in meetings A and B
	Took part in the computation of reference values
Maarten Wyns	One of two behind code for ADI (SABR and QLSV models)
	Took part in meeting B
Alvaro Leitao	One of two behind code for mSABR (SABR model)
	Behind code for MLMC (SABR model)
	Took part in meeting C
Shashi Jain	One of two behind code for SGBM (QLSV model)
Tinne Haentjens	One of two behind code for ADI (HHW model)
Johan Waldén	Took part in the writing of Sections 1 and 2
	Took part in meeting A

During the time for the project, three meetings were held:

1. A kick-off meeting in Uppsala in August 2016, dedicated for the set-up of the project.

2. An intermediate planning meeting in Austin in connection with the SIAM Conference on Financial Mathematics & Engineering in November 2016.

3. A final planning meeting in Leiden in connection with the Lorentz center workshop in Applied Mathematics Techniques for Energy Markets in Transition in September 2017.

In Table 9 it is presented at what meetings the authors took part.