# A second order finite volume IMEX Runge-Kutta scheme for two dimensional PDEs in finance

José Germán López-Salas, María Suárez Taboada, Manuel Jesús Castro Díaz,

Ana Ferreiro Ferreiro, José Antonio García Rodríguez

► To cite this version:

José Germán López-Salas, María Suárez Taboada, Manuel Jesús Castro Díaz, Ana Ferreiro Ferreiro, José Antonio García Rodríguez. A second order finite volume IMEX Runge-Kutta scheme for two dimensional PDEs in finance. 2024. hal-04720735

HAL Id: hal-04720735

https://hal.science/hal-04720735v1

Preprint submitted on 4 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A second order finite volume IMEX Runge-Kutta scheme for two dimensional PDEs in finance*

J. G. López-Salasª, M. Suárez-Taboadaª, M. J. Castroᵇ,*, A. M. Ferreiro-Ferreiroª, J. A. García-Rodríguezª

## Abstract

In this article we present a novel and general methodology for building second order finite volume implicit-explicit (IMEX) numerical schemes for solving two dimensional financial parabolic PDEs with mixed derivatives. In particular, applications to basket and Heston models are presented. The obtained numerical schemes have excellent properties and are able to overcome the well-documented difficulties related with numerical approximations in the financial literature. The methods achieve true second order convergence with non-regular initial conditions. Besides, the IMEX time integrator allows to overcome the tiny time-step induced by the diffusive term in the explicit schemes, also providing very accurate and non-oscillatory approximations of the Greeks. Finally, in order to assess all the aforementioned good properties of the developed numerical schemes, we compute extremely accurate semi-analytic solutions using multi-dimensional Fourier cosine expansions. A novel technique to truncate the

<footnote>ªDepartment of Mathematics, Faculty of Informatics and CITIC, Campus Elviña s/n, 15071-A Coruña (Spain)</footnote>

<footnote>ᵇDepartment of Análisis Matemático, Facultad de Ciencias, University of Málaga, Campus de Teatinos s/n, Málaga, 29080-Andalucía (Spain)</footnote>

<footnote>*M.J Castro research has been partially supported by the Spanish Government and FEDER through the coordinated Research project RTI2018-096064-B-C1, the Junta de Andalucía research project P18-RT-3163, the Junta de Andalucía-FEDER-University of Málaga research project UMA18-FEDERJA-16 and the University of Málaga. The other authors research has been partially supported by the Spanish MINECO under research project number PDI2019-108584RB-I00 and by the grant ED431G 2019/01 of CITIC, funded by Consellería de Educación, Universidade e Formación Profesional of Xunta de Galicia and FEDER.</footnote>

<footnote>*Corresponding author.</footnote>

<footnote>Email addresses: jose.lsalas@udc.es (J. G. López-Salas), maria.suarez3@udc.es (M. Suárez-Taboada), mjcastro@uma.es (M. J. Castro), ana.fferreiro@udc.es (A. M. Ferreiro-Ferreiro), jose.garcia.rodriguez@udc.es (J. A. García-Rodríguez)</footnote>
Fourier series for basket options is presented and it is efficiently implemented using multi-GPUs.

*Keywords:* finite volume method, option pricing, advection-diffusion, Runge-Kutta, IMEX, Fourier COS method.

2020 MSC: 65M08,91G20,91G80,65L04,65L12

# 1. Introduction

The goal of this article is to develop a novel methodology for pricing financial options. A general framework for building second order numerical schemes for financial advection-diffusion-reaction PDEs with mixed derivatives is designed. The technique is based in the use of finite volume implicit-explicit (IMEX) Runge-Kutta numerical methods. This approach is applied for developing a second order numerical method for solving basket and Heston models in two spatial dimensions. To our knowledge, this is the first time that finite volume IMEX Runge-Kutta numerical schemes have been successfully applied to the numerical solution of parabolic financial PDEs. The interplay between finite volume methods and IMEX time integrators is an extremely powerful combination in mathematical finance. On the one hand, finite volume methods are well suited to treat convective terms. In fact, it allows to build schemes with true second order convergence, even in the presence of non-regular initial or boundary conditions, which is a well-known difficulty in the financial literature. On the other hand, IMEX time integrator makes possible to overcome the tiny time steps induced by the spatial semi-discretization of the diffusive terms. This is also of paramount importance, otherwise the problem is intractable from the computational point of view. On top of that, the presented numerical schemes provide excellent finite difference approximations of the Greeks. Actually, these approximations are free of the spurious oscillations which typically appear when no strongly stable time integrators, such as Crank-Nicolson, are considered. Last but not least, in order to assess the robustness, accuracy and order of convergence of the proposed finite volume IMEX numerical schemes,
the multidimensional Fourier cosine method is implemented. A novel way to truncate the series expansion for basket options is presented. Additionally, an efficient parallel implementation in a multi-GPU environment allows to recover extraordinarily accurate option prices.

Option pricing models are of great importance in the financial industry. Several types of stochastic models can be found in the literature for modelling the price of an option. These models are formulated as stochastic differential equations (SDEs) or systems of SDEs. They are based on Brownian motions (or Wiener processes) that are used to simulate the evolution of the underlying assets and sometimes also their volatilities.

Option pricing methods can be mainly classified into two categories, deterministic and stochastic algorithms. The majority of the first ones deal with the solution of PDEs, while the second ones employ the Monte Carlo method based on the formulation of the model in terms of expectations. As we have said before, in this article we will focus on numerical methods to solve PDEs. In next paragraphs we point out the main drawbacks of Monte Carlo with respect to PDE based numerical methods.

Firstly, we note that Monte Carlo rate of convergence is very slow, of order $O(1/\sqrt{S})$ for all dimensions, $S$ being the number of simulations (see [1]). Secondly, if the risk neutral distribution of the underlying at the expiry of the option is known explicitly, there is no need to discretize the SDE to perform Monte Carlo simulation. In that situation, Monte Carlo will be faster than PDE methods, specially in high dimensions, since it does not suffer the curse of dimensionality. However, in several option pricing models this is not the case, so that the SDEs have to be discretized. In this procedure we must choose a suitable time step $\Delta t$ for the discretization of the SDE, so that the discretization error is of the same order of magnitude as the standard error of the simulation. Therefore, when increasing the number of simulations $S$, one must increase the number of time steps so that the time discretization error remains of the same order as the standard error of the simulation. Monte Carlo theory suggests that $\Delta t$ should
be proportional to $\frac{1}{\sqrt{s}}$. Thus, in a simulation with $S$ paths we take $\Delta t = \kappa \frac{1}{\sqrt{s}}$ for some constant of proportionality $\kappa$. Summing up, in order to increase Monte Carlo accuracy one should need to increase the number of simulations which in turn requires a decrease in the discretization time step, and leads to a very slow convergence.

Additionally, in the context of options pricing, standard Monte Carlo methods have difficulties in providing accurate sensitivities (Greeks) of the option prices with respect to the underlying stocks or parameters, which are important for trading purposes. These Greeks are directly given by the discretization of the partial derivatives involved in the PDE or can be easily obtained by numerical differentiation formulas. Also, early exercise features like in American or Bermudan options do not fit into the forward simulation paradigm of pure Monte Carlo and the combination with regression techniques turns out very costly (see [2] for details), whereas numerical methods for complementarity PDE problems are much more straightforward.

In this work we will focus in pricing options under financial models with two spatial variables. More precisely, the models we will be dealing with can be written as a 2D advection-diffusion-reaction PDE of the general form

$$ \frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} + b \frac{\partial u}{\partial y} + c \frac{\partial^2 u}{\partial x^2} + d \frac{\partial^2 u}{\partial y \partial x} + e \frac{\partial^2 u}{\partial y^2} + f = 0, \quad (1) $$

where $a, b, c, d, e$ are functions of $t, x$ and $y$, and $f = f(t, x, y, u)$.

Standard numerical methods for solving parabolic PDEs in finance, of the type (1), are finite difference (FD) based, see [3, 4, 5, 6], or finite elements (FE) based, see [7, 3, 8, 9]. Another possibility is to use Exponentially Fitted (EF) finite differences schemes (see [4]) or Alternate Directions (ADI) finite differences methods (see [10, 11]).

Nevertheless, standard FD or FE numerical schemes for financial PDEs have to face also several mathematical challenges. The first one is due to the lack of smoothness of the terminal conditions imposed by most of the traded option contracts. A financial option is a contract that gives the holder the right to trade in the future at a previously agreed exercise price. However, the option holder
is under no obligation and will only exercise if it is worth it. This optionality
in the exercise is represented mathematically in the payoff function of the con-
tract, commonly by means of the max function, i.e max(·, 0). As a consequence,
terminal conditions of PDE models, and their derivatives, are commonly non-
smooth. This non-smooth data originates noise in the numerical solution. As
a result, the expected rate of convergence is seriously degraded. Besides, the
numerical computation of derivatives of the numerical solution originates poor
estimates of the Greeks with oscillations. Several strategies were designed in
the literature in order to overcome these issues. For instance, smoothing meth-
ods for the terminal conditions were proposed. Heston and Zhou [12] replace
mesh nodal values with an average value computed over a close space. This
technique is natural in finite volume methods. In [6], Tavella and Randall pro-
posed to work with meshes such that the strike price is halfway between grid
points. Another method consists of projecting the terminal conditions onto a
set of basis functions, see [13]. Nevertheless, in general a smoothing technique
must be combined with thorough time discretizations in order to recover sec-
ond order convergence for option pricing problems with discontinuous payoff
functions (or its derivatives). Along this line, Rannacher [14, 13] proposed a
numerical scheme in which a finite number of fully implicit time steps are per-
formed before the Crank-Nicolson time discretization. This approach, known as
Rannacher timestepping, along with a suitable smoothing technique of the ter-
minal condition, preserves the quadratic convergence rate of the Crank-Nicolson
scheme (see [15, 16] for more details).

Another mathematical challenge in the numerical solution of this kind of
advection-diffusion PDEs arises when the advection term becomes larger than
the diffusion one, or when we have degenerated diffusion. Under this situa-
tion, instabilities show up, because the problem becomes more hyperbolic. In
order to face these problems, several techniques have been introduced. One
way to overcome these instability phenomena is to avoid centered schemes and
to consider upwind discretizations of advection terms, thus taking information
upstream and not downstream. We refer the reader to [3] for further details in
the mathematical finance field. Although upwind schemes are stable, they are
less accurate than the centered ones. Another approach to face instabilities in
convection-dominated diffusion problems is the semi-Lagrangian time stepping
method, that integrate the equations by searching backward in time the position
of the flow. Interpolation is needed in this tracing procedure along characteristic
curves. Several works in the literature combine the method of characteristics
with FE or FD procedures in option pricing problems. In [17], Forsyth et al.
price American Asian options under jump diffusion. In [8] the authors price
Asian options under a semi-Lagrangian finite element method. In [18], Ferreiro
et al. present a semi-Lagrangian finite difference method for pricing business
companies. The main disadvantages of semi-Lagrangian methods is the diffi-
culty to build high order numerical schemes. In fact, as stated in [17], due to
the non-smoothness of the payoff, “even if second order time stepping methods
are used, observed convergence as the mesh and time step is refined occurs at a
sub-second order rate”. Besides, the computational cost of this strategy can be
high due to the backward search in time. In fact, at each time step an initial
value problem for each point of the spatial mesh has to be solved, and then the
position of the point has to be traced backward in time. Finally, interpolation
with high enough order has to be performed. Another approach to properly solve
convection-dominated diffusion problems is to consider finite volume discretiza-
tion methods. Finite volume methods naturally handle non smooth payoffs and
the cases where the PDE becomes diffusion degenerated. Seminal works on ap-
plying finite volume method to option pricing problems are due to Forsyth and
Zvan, see [19, 20], where vertex type finite volume methods were applied for
problems in two spatial dimensions. More recently, in [21], finite volume meth-
ods were proposed for convection dominated pricing problems. More precisely,
the authors applied the Kurganov-Tadmor scheme proposed in [22, 23] to option
pricing problems written in conservative form. Later, in [24] the authors propose
a second order improvement of [21] with appropriate time methods and slope
limiters. Finally in [25] the authors apply a third order Kurganov-Levy scheme
presented in [26] based in the CWENO reconstructions presented in [27]. All
these works only deal with Black-Scholes PDEs in dimension one or Asian PDE problems. Therefore, they do not solve general two-dimensional pricing problems with mixed derivatives. In these works, the ordinary differential equations (ODEs) obtained after semi-discretizing in space are stiff. In fact, the authors use explicit-schemes in time. The main disadvantage of these approaches is that time steps have to be very small to maintain stability. Typically, Von Neumann stability analysis demands $\Delta t \le \frac{(\Delta x)^2}{2\eta}$, being $\eta$ the diffusion velocity.

Pareschi and Russo proposed in [28] implicit-explicit (IMEX) Runge-Kutta method for general hyperbolic systems of conservation laws with stiff relaxation terms. The idea is to apply an implicit discretization to the source terms and an explicit one to the nonstiff term. In [29, 30], the IMEX time marching scheme was used in mathematical finance only combined with finite differences and finite elements, respectively.

All in all, our goal is to first semi-discretize in space 2D Black-Scholes equations (having mixed derivatives) with second-order finite volume methods, thus properly treating convection terms and non-smooth payoffs. Later, we propose to integrate in time the resulting system of stiff ODEs by means of the second-order IMEX Runge-Kutta time marching scheme. In this way, we will apply an implicit discretization to the diffusion (stiff) part and an explicit one to the convection and source terms (non stiff). As a result, the time step restriction will depend only on the stability of the degenerate diffusion PDE, i.e., $\Delta t \le \frac{\Delta x}{\alpha}$, being $\alpha$ the convection velocity. The here proposed numerical methods will provide also good approximations of the Greeks (without oscillations).

The organization of this paper is as follows. In Section 2 we overview a couple of two-dimensional Black-Scholes PDE problems with huge interest in the financial industry. Section 3 presents the proposed finite volume IMEX Runge-Kutta scheme for general two-dimensional convection diffusion PDEs written in conservative form. Section 4 is devoted to an alternative Fourier numerical method designed to solve the PDE problems presented in Section 2. In order to obtain a high-accurate reference solution to validate the here proposed finite volume methods, a novel and very efficient multi-GPU implementation of the
Fourier method is carried out. In Section 5 numerical experiments are performed to price several options under the Black-Scholes and Heston models. A thorough inspection of the accuracy and order of convergence of the developed numerical methods is detailed. Finally, conclusions are drawn in Section 6.

Notations of the paper.

(i) $\ln(x)$ stands for the natural logarithm of $x \in \mathbb{R}_+$.

(ii) Re$(z)$ denotes the real part of the complex $z \in \mathbb{C}$; we set $\mathbf{i} := \sqrt{-1}$.

(iii) $N$ is the standard normal cumulative distribution function.

(iv) $\mathbf{x}^\mathrm{T}$ stands for the transpose of the vector $\mathbf{x} = (x_1, \dots, x_n) \in \mathbb{R}^n$.

## 2. Option pricing PDE models

The simplest option is a contract giving its owner the right (not the obligation) to buy (call option) or sell (put option) a particular asset at a fixed price at a specified time in the future. Therefore, the option contract has to specify three terms. The first one is the underlying asset; the price of this asset will be denoted by $s_t$. The second one is the fixed price, the so-called exercise price or strike price; it will be denoted by $K$. Finally, the given date is termed the maturity or the expiration date; it is often denoted by $T$. If the right to exercise the option is only permitted at maturity, the option is called European. This kind of options are basic derivatives, and are often called plain vanilla options. Option pricing is one of the many problems of financial mathematics. The price of the option at maturity is straightforward. If $s_T$ is the price of the underlying asset at maturity, the payoff of a call option is $\max(s_T - K, 0)$, whilst $\max(K - s_T, 0)$ is the benefit of the option's holder at the option's expiry. Nevertheless, determining the price that the buyer of the option has to pay to the seller at the initial time $t$, the so-called premium, is a challenging mathematical problem because $s_T$ is not known at time $t$. We will use $u$ to mean the value of the option. It is clear that $u$ is a function of the value of the underlying asset $s$ at time $t$. The goal of option pricing PDE models is to answer the following question: what is the market price $u(s, t)$ of the option at time $t$, $0 \le t < T$?
The theory of option pricing started in the seventies with the celebrated works of Black and Scholes [31] and Merton [32]. The Black-Scholes model assumes that the price of the underlying asset is the solution of the following stochastic differential equation

$$ \frac{ds_t}{s_t} = rdt + \sigma dW_t, $$

where $W$ is a standard Brownian motion, $r \in \mathbb{R}$ is the interest rate and $\sigma \in \mathbb{R}_+$ is the volatility of the underlying asset. Therefore, the vanilla call option price is given by the following expectation under the risk-neutral measure

$$ u(s_t, t) = e^{-r(T-t)} \mathbb{E}[\max(s_T - K, 0)]. $$

Besides $u$ solves the backward PDE

$$ \frac{\partial u}{\partial t} + \frac{1}{2}\sigma^2 s^2 \frac{\partial^2 u}{\partial s^2} + rS \frac{\partial u}{\partial s} - ru = 0, \quad t \in [0, T). $$

In this work we will work on deterministic methods for pricing options, thus proposing appropriate numerical methods for solving parabolic advection-diffusion PDEs. However, we will deal with more complicated contracts and models than the one described above, so that the resulting PDE becomes two-dimensional in space. More precisely, we are interested in pricing options on many underlyings, the so-called basket options, and also in pricing vanilla options under stochastic volatility models.

On the one hand, an arithmetic basket call option has a payoff of

$$ \max(c_1s_{1T} + c_2s_{2T} - K, 0), \quad c_1, c_2 \in \mathbb{R}. $$

The Black-Scholes model results in closed form solutions for pricing vanilla options, due to the fact that the underlying asset price at fixed time follows a lognormal distribution. Nevertheless, considering the Black-Scholes model for several underlying assets does not lead to closed form formulas for the price of a basket option. The obstacle is mainly the fact that the distribution of a weighted average of lognormals is not known. Thus, numerical methods are required to approximate the price of arithmetic basket options.
On the other hand, the Black-Scholes model assumes that $\sigma$, the volatility of the underlying, is constant. However, this model cannot recover the market volatility smile and skew, which are behaviours of the market implied volatility surface (see [33] and references therein for details). In order to overcome this limitation of constant volatility models, stochastic volatility models were proposed. They change our model for the evolution of the underlying asset by adding a new variable for the volatility which satisfy some SDE. The most popular and renowned of all stochastic volatility models is the Heston model proposed in [34]. Pricing vanilla options under this model also requires numerical methods if we do not want to depend on quasi-closed form solutions.

In the next two subsections the problems we want to solve in this article are stated explicitly.

## 2.1. Options on a basket of two assets

We consider a basket of two assets, whose prices are $s_1$ and $s_2$. Under the Black-Scholes model, the prices of the underlying assets follow the following system of stochastic differential equations:

$$ 
\begin{aligned}
\frac{ds_{1t}}{s_{1t}} &= (r - q_1)dt + \sigma_1 dW_{1t}, \quad s_{10} \text{ known,} \\
\frac{ds_{2t}}{s_{2t}} &= (r - q_2)dt + \sigma_2 dW_{2t}, \quad s_{20} \text{ known,}
\end{aligned}
\qquad (2) $$

where $W$ is a two dimensional correlated Brownian motion. We set

$$ W_{1t} = \bar{W}_{1t}, \quad W_{2t} = \rho \bar{W}_{1t} + \sqrt{1 - \rho^2} \bar{W}_{2t}, \qquad (3) $$

where $(\bar{W}_1, \bar{W}_2)$ is a two dimensional standard Brownian motion and $\rho \in (-1, 1)$ is the constant correlation parameter. Besides, $\sigma_1, \sigma_2 \in \mathbb{R}_+$ are the market volatilities of the assets $s_1$ and $s_2$, respectively. Finally, $r$ is the interest rate, and $q_1$ and $q_2$ are the dividend yields of the assets. The price of an option with maturity $T$ and payoff function $u_T(s_1, s_2)$ is

$$ u(s_{1t}, s_{2t}, t) = e^{-r(T-t)} \mathbb{E}^Q [u_T(s_{1T}, s_{2T})], $$
where $Q$ is the selected martingale measure. By applying the two-dimensional
Itô formula, we can find the PDE for the price of the option $u(s_1, s_2, t)$:

$$
\begin{align}
& \frac{\partial u}{\partial t} + \frac{1}{2}\sigma_1^2 s_1^2 \frac{\partial^2 u}{\partial s_1^2} + \frac{1}{2}\sigma_2^2 s_2^2 \frac{\partial^2 u}{\partial s_2^2} + \rho\sigma_1\sigma_2 s_1 s_2 \frac{\partial^2 u}{\partial s_1 \partial s_2} \nonumber \\
& \qquad + (r - q_1)s_1 \frac{\partial u}{\partial s_1} + (r - q_2)s_2 \frac{\partial u}{\partial s_2} - ru = 0, \quad t \in [0, T), \nonumber \\
& u(s_1, s_2, T) = u_T(s_1, s_2), \quad s_1, s_2 > 0. \tag{4}
\end{align}
$$

For convenience, a time reversal $\tau = T-t$ is performed. Doing so, the backward
PDE (4) is then transformed to the following forward PDE

$$
\begin{align}
& \frac{\partial u}{\partial \tau} - \frac{1}{2} \sigma_1^2 s_1^2 \frac{\partial^2 u}{\partial s_1^2} - \frac{1}{2} \sigma_2^2 s_2^2 \frac{\partial^2 u}{\partial s_2^2} - \rho \sigma_1 \sigma_2 s_1 s_2 \frac{\partial^2 u}{\partial s_1 \partial s_2} \nonumber \\
& \qquad - (r-q_1) s_1 \frac{\partial u}{\partial s_1} - (r-q_2) s_2 \frac{\partial u}{\partial s_2} + ru = 0, \quad \tau \in (0, T], \nonumber \\
& u(s_1, s_2, 0) = u_0(s_1, s_2), \quad s_1, s_2 > 0. \tag{5}
\end{align}
$$

In the sequel, the $\tau$ notation is dropped for simplicity and the forward time is
again written as $t$. In this work we consider the following payoff function for
the arithmetic basket call option

$$
u_0(s_1, s_2) = \max \left( \frac{1}{2}(s_1 + s_2) - K, 0 \right), \quad (6)
$$

where $K$ is the fixed strike price. As previously said, it is not possible to compute
the value of an arithmetic basket call option analytically. We will compare the
results of the proposed finite volume method with reference values calculated
with another method, namely the COS method.

In order to compute a numerical approximation of $u$, we need to truncate
the domain in the variables $s_1$ and $s_2$, so $u$ will be computed for $s_1, s_2 \in (0, \bar{S})$,
with $\bar{S}$ large enough. Besides, boundary conditions have to be imposed at the
artificial boundaries $s_1 = s_2 = \bar{S}$. The PDE (5) is degenerate on the boundaries
$s_1, s_2 = 0$ and is reduced into a Black-Scholes 1D PDE. For the arithmetic
basket call option with payoff function (6) the following Dirichlet and Neumann
boundary conditions can be used:

$$u(0, s_2, t) = C\left(\frac{1}{2}s_2, t, \sigma_1, q_1\right), \quad \frac{\partial u}{\partial s_1}(\bar{S}, s_2, t) = \frac{1}{2},$$

$$u(s_1, 0, t) = C\left(\frac{1}{2}s_1, t, \sigma_2, q_2\right), \quad \frac{\partial u}{\partial s_2}(s_1, \bar{S}, t) = \frac{1}{2},$$

where C is the Black-Scholes 1D price of the vanilla call European option, given by

$$C(S, t, \sigma, q) = Se^{-qt}N(d_1) - Ke^{-rt}N(d_2),$$

$$d_1 = \frac{\ln(\frac{S}{K}) + (r - q + \frac{\sigma^2}{2})t}{\sigma\sqrt{t}}, \quad d_2 = d_1 - \sigma\sqrt{t}.$$

## 2.2. Heston model

In the Heston model [34], there is an underlying asset s whose volatility is a stochastic process driven by a second Brownian motion:

$$ds_t = (r - q)s_t dt + \sqrt{v_t}s_t dW_{1t}, \quad s_0 \text{ known},$$

$$dv_t = \kappa(\theta - v_t)dt + \sigma\sqrt{v_t}dW_{2t}, \quad v_0 \text{ known}, \qquad (7)$$

where r, q, κ, θ, σ, s₀, v₀ are constant parameters in ℝ⁺ and W is a two dimensional correlated Brownian motion, settled as in (3). Moreover, r is the fixed interest rate, q is the constant dividend yield, κ is the mean-reversion speed for the variance, θ is the mean reversion level for the variance, σ is the volatility of the variance (the so-called volatility of the volatility) and v₀ is the initial level of the variance. The process v represents the variance of s and its stochastic differential equation is a version of the square root process described by Cox, Ingersoll and Ross [35].

The price of a derivative with payoff function $u_T(s)$ is the solution of the following backward PDE:

$$\frac{\partial u}{\partial t} + \frac{1}{2}s^2 v \frac{\partial^2 u}{\partial s^2} + \rho \sigma s v \frac{\partial^2 u}{\partial v \partial s} + \frac{1}{2}\sigma^2 v \frac{\partial^2 u}{\partial v^2} \\ + (r-q)s \frac{\partial u}{\partial s} + \kappa(\theta-v)\frac{\partial u}{\partial v} - ru = 0, \quad t \in [0, T), \\ u(s, v, T) = u_T(s), \quad s > 0. \qquad (8)$$
Once more, if we perform a time reversal $\tau = T - t$, the backward PDE (8) is transformed to the following forward PDE

$$ \frac{\partial u}{\partial \tau} - \frac{1}{2} s^2 v \frac{\partial^2 u}{\partial s^2} - \rho \sigma s v \frac{\partial^2 u}{\partial v \partial s} - \frac{1}{2} \sigma^2 v \frac{\partial^2 u}{\partial v^2} \\ - (r-q) s \frac{\partial u}{\partial s} - \kappa(\theta - v) \frac{\partial u}{\partial v} + ru = 0, \quad \tau \in (0, T], \\ u(s, v, 0) = u_0(s), \quad s > 0. \tag{9} $$

As before, the forward time $\tau$ will be written as $t$. The payoff function of an European call option is

$$ u_0(s) = \max(s - K, 0), \tag{10} $$

where $K$ is the strike. Since the natural logarithmic characteristic function of $s$, $E[e^{i\xi \ln s_T}]$, is known explicitly, the price of the call option can be also computed by means of Fourier methods.

The finite volume numerical approximation of $u$ will be computed in a truncated domain, $s \in (0, \bar{S})$ and $v \in (0, \bar{V})$. In this work, the following boundary conditions were imposed

$$ u(0, v, t) = 0, \quad \frac{\partial^2 u}{\partial s^2}(\bar{S}, v, t) = 0, \quad \frac{\partial^2 u}{\partial v^2}(s, \bar{V}, t) = 0. $$

At the boundary $v = 0$ no condition is specified. If the so-called Feller condition is satisfied, $2\kappa\theta > \sigma^2$, it follows that this is and outflow boundary.

### 3. Finite volume IMEX Runge-Kutta numerical method

Equation (1) can be written in conservative form in the following compact way:

$$ \frac{\partial}{\partial t} u(x, y, t) + \frac{\partial f_1}{\partial x}(u) + \frac{\partial f_2}{\partial y}(u) = \frac{\partial g_1}{\partial x}(u_x, u_y) + \frac{\partial g_2}{\partial y}(u_x, u_y) + h(u), \tag{11} $$

or

$$ \frac{\partial}{\partial t} u(x, y, t) + \operatorname{div} \mathbf{F}(u) = \operatorname{div} \mathbf{G}(\nabla u) + h(u), \tag{12} $$

with $\mathbf{F}(u) = (f_1(u), f_2(u))$ and $\mathbf{G}(\nabla u) = (g_1(\nabla u), g_2(\nabla u))$. More precisely, $\mathbf{F} : \mathbb{R}^3 \to \mathbb{R}^2$, $\mathbf{G} : \mathbb{R}^4 \to \mathbb{R}^2$ and $h : \mathbb{R}^3 \to \mathbb{R}$ are functions of $x, y, u, \nabla u$,
although in order to keep the notation simple some of these dependencies were omitted.

The numerical solution of equation (12) using a finite volume semi-implicit scheme is difficult because of the presence of the diffusive part. In this article we have considered the implicit-explicit (IMEX) Runge-Kutta time discretization numerical scheme proposed in [28]. This time-marching scheme play a major rule in the treatment of stiff systems of differential equations in the form

$$ \frac{\partial U}{\partial t} + E(U) = I(U), \qquad (13) $$

where $U = U(t) \in \mathbb{R}^N$ and $E, I: \mathbb{R}^N \to \mathbb{R}^N$, being $E$ the non-stiff term and $I$ the stiff part. Both parts will be handled simultaneously with the same IMEX solver.

The rest of the section is organized as follows. In Section 3.1 we describe the space discretization obtained by finite volume schemes. Section 3.2 is devoted to the IMEX Runge-Kutta scheme applied to the stiff system of differential equations (13) obtained by the finite volume space discretization.

### 3.1. Space discretization

Let $\Delta x, \Delta y$ be the mesh length in the $X, Y$ directions respectively. We define the grid points $x_i = i\Delta x, y_j = j\Delta y, i, j \in \mathbb{Z}$. Let $x_{i+1/2} = x_i + \frac{1}{2}\Delta x$ and $y_{j+1/2} = y_j + \frac{1}{2}\Delta y$. We consider rectangular finite volumes $V_{ij} = [x_{i-1/2}, x_{i+1/2}] \times [y_{j-1/2}, y_{j+1/2}]$, where $(x_i, y_j)$ is the center of the finite volume $V_{ij}$. The area of $V_{ij}$ will be denoted as $|V_{ij}|$, i.e. $|V_{ij}| = \Delta x \Delta y$. Besides, $\Gamma_{ij}$ represents the boundary of $V_{ij}$. Let $\bar{u}_{ij} = \frac{1}{|V_{ij}|} \int_{V_{ij}} u(x, y, t) dxdy$ be the volume average. The here described finite volume grid is sketched in Figure 1.

Integrating equation (12) in space on $V_{ij}$ and dividing by $|V_{ij}|$ we obtain the semi-discrete equation

$$ \frac{d\bar{u}_{ij}}{dt} = - \frac{1}{|V_{ij}|} \int_{V_{ij}} \mathrm{div} \mathbf{F}(u) dxdy + \frac{1}{|V_{ij}|} \int_{V_{ij}} \mathrm{div} \mathbf{G}(u_x, u_y) dxdy \\ + \frac{1}{|V_{ij}|} \int_{V_{ij}} h(u) dxdy. \qquad (14) $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Finite volume stencil.</figcaption></figure>

Applying the divergence theorem we can rewrite the first two volume integrals of (14) as line integrals

$$ \frac{d\bar{u}_{ij}}{dt} = - \frac{1}{|V_{ij}|} \oint_{\Gamma_{ij}} \mathbf{F}(u) \cdot \mathbf{n} \, d\gamma + \frac{1}{|V_{ij}|} \oint_{\Gamma_{ij}} \mathbf{G}(u_x, u_y) \cdot \mathbf{n} \, d\gamma \\ + \frac{1}{|V_{ij}|} \int_{V_{ij}} h(u) \, dxdy. \qquad (15) $$

The unknowns of the problem are the volume averages {$\bar{u}_{ij}(t)$}. Therefore, in order to convert (15) into a numerical scheme, we have to approximate on the right hand side of this equation with functions of {$\bar{u}_{ij}(t)$}. In the following two content blocks, the numerical treatment of the source and advection terms (Section 3.1.1) and the diffusion ones (Section 3.1.2) is discussed. The source and convective terms will be treated explicitly by IMEX Runge-Kutta, while the diffusion part will be managed implicitly in time.
### 3.1.1. Explicit part: advection-reaction space discretization

For the advective part one readily gets

$$ \oint_{\Gamma_{ij}} \mathbf{F}(u) \cdot \mathbf{n} \, d\gamma = \int_{\Gamma_{i+1/2,j}} f_1(u) \, d\gamma_j + \int_{\Gamma_{i-1/2,j}} f_1(u) \, d\gamma_j + \\ \qquad \int_{\Gamma_{i,j+1/2}} f_2(u) \, d\gamma_i + \int_{\Gamma_{i,j-1/2}} f_2(u) \, d\gamma_i, \quad (16) $$

where for second order schemes each line integral can be approximated by means of midpoint quadrature rule as follows

$$ \int_{\Gamma_{i \pm 1/2, j}} f_1(u) \, d\gamma_j \approx \pm \Delta y f_1(u(x_{i \pm 1/2}, y_j)), \\ \int_{\Gamma_{i,j \pm 1/2}} f_2(u) \, d\gamma_i \approx \pm \Delta x f_2(u(x_i, y_{j \pm 1/2})). $$

At this point the unknown function $u(x, y, t)$ is reconstructed by a piecewise polynomial using the volume averages $\{\bar{u}_{ij}(t)\}$. More precisely, starting from $\{\bar{u}_{ij}\}$, we compute a piecewise polynomial reconstruction

$$ \mathcal{R}(x, y) = \sum_{i,j} P_{ij}(x, y) \mathbf{1}_{ij}(x, y), $$

where $P_{ij}$ is a suitable polynomial satisfying some accuracy and non oscillatory property, and $\mathbf{1}_{ij}$ is the indicator function of the volume $V_{ij}$. Second order numerical schemes can be obtained by means of piecewise linear polynomials, although higher order schemes can be obtained by polynomials of higher order. Here we consider the natural extension of MUSCL reconstruction to 2D Cartesian grids (see [36]).

The flux functions at the midpoints of the boundaries of the volumes can be computed by using a suitable numerical flux function, consistent with the analytical flux. Hereafter we detail the approximation of $f_1(u(x_{i+1/2}, y_j))$, being the other fluxes approximated in the same way:

$$ f_1(u(x_{i+1/2}, y_j)) \approx \mathcal{F}_1(u_{i+1/2,j}^{-}, u_{i+1/2,j}^{+}). $$

The values $u_{i+1/2,j}^{\pm}$ are obtained from the reconstruction as

$$ u_{i+1/2,j}^{\pm} = \lim_{x \to x_{i+1/2}^{\pm}} \mathcal{R}(x, j\Delta y), $$
with $x$ in a normal line to the boundary $\Gamma_{i+1/2,j}$. For example the left reconstructed value at the edge $\Gamma_{i+1/2,j}$ is

$$u_{i+1/2,j}^{-} = \bar{u}_{i,j} + \frac{\Delta x}{2} u'_{i,j},$$

where the slope $u'_{i,j}$ is a first order approximation of the space derivative in the $X$ direction of $u(x, y, t)$ at point $(x_i, y_j)$ at every time $t$. This slope must satisfy the TVD property and thus we must use slope limiters. In our case we use the minmod limiter, where the slope is given by

$$u'_{ij} = \frac{1}{\Delta x} \minmod(\bar{u}_{i,j} - \bar{u}_{i-1,j}, \bar{u}_{i+1,j} - \bar{u}_{i,j}),$$

$$\minmod(a, b) = \begin{cases} \min(a, b) & \text{if } a, b > 0, \\ \max(a, b) & \text{if } a, b < 0, \\ 0 & \text{otherwise.} \end{cases}$$

All the calculations performed in this article have been carried out using the CIR numerical flux

$$\mathcal{F}_1(u_{i+1/2,j}^{-}, u_{i+1/2,j}^{+}) = \frac{1}{2}(f_1(u_{i+1/2,j}^{-}) + f_1(u_{i+1/2,j}^{+})) - \frac{1}{2}\alpha(u_{i+1/2,j}^{+} - u_{i+1/2,j}^{-}),$$

$$\text{where } \alpha = \left| f'_{1} \left( \frac{u_{i+1/2,j}^{-} + u_{i+1/2,j}^{+}}{2} \right) \right|.$$

Therefore, the line integral (16) of the convective term, is finally approximated as

$$\int_{\Gamma_{ij}} \mathbf{F}(u) \cdot \mathbf{n} \, d\gamma \approx \Delta y (\mathcal{F}_1(u_{i+1/2,j}^{-}, u_{i+1/2,j}^{+}) - \mathcal{F}_1(u_{i-1/2,j}^{-}, u_{i-1/2,j}^{+})) + \\ \Delta x (\mathcal{F}_2(u_{i,j+1/2}^{-}, u_{i,j+1/2}^{+}) - \mathcal{F}_2(u_{i,j-1/2}^{-}, u_{i,j-1/2}^{+})), \quad (17)$$

being $\mathcal{F}_1$ and $\mathcal{F}_2$ the numerical fluxes of the physical flux functions $f_1$ and $f_2$, respectively.

Finally, the volume integral of the source term can be discretized using the midpoint quadrature rule

$$\int_{V_{ij}} h(u) dxdy \approx |V_{ij}| h(\bar{u}_{ij}). \qquad (18)$$
### 3.1.2. Implicit part: diffusion space discretization

The diffusion part is:

$$ \int_{\Gamma_{ij}} \mathbf{G}(u_x, u_y) \cdot \mathbf{n} \, d\gamma = \int_{\Gamma_{i+1/2, j}} g_1(u_x, u_y) \mathrm{d}\gamma_j + \int_{\Gamma_{i-1/2, j}} g_1(u_x, u_y) \mathrm{d}\gamma_j + \\ \int_{\Gamma_{i,j+1/2}} g_2(u_x, u_y) \mathrm{d}\gamma_i + \int_{\Gamma_{i,j-1/2}} g_2(u_x, u_y) \mathrm{d}\gamma_i. \quad (19) $$

Before approximating each one of these line integrals, a suitable approximation of the partial derivatives $u_x$ and $u_y$ has to be built for the volume $V_{ij}$. With this aim, we build the second order Lagrange interpolating polynomial of $u$ centered in the volume $V_{ij}$. Let $L_{ij}^u$ be that polynomial. Considering the nine nodes of the dual mesh $\{x_{i+k}, y_{j+l}\}$, $k, l = -1, 0, 1$, and the averaged values $\{\bar{u}_{i+k,j+l}\}$ of the solution at each volume $V_{i+kj+l}$, this polynomial is given by:

$$ L_{ij}^{\bar{u}}(x, y) = \sum_{k,l=-1}^{1} \bar{u}_{i+k,j+l} \ell_{i+k}(x) \ell_{j+l}(y), $$

where $\ell_{i+k}$ and $\ell_{j+l}$ are the one dimensional Lagrange polynomial basis, i.e:

$$ \ell_{i+k}(x) = \prod_{\substack{p=i-1 \\ p \neq i+k}}^{i+1} \frac{x-x_p}{x_{i+k}-x_p}, \quad \ell_{j+l}(y) = \prod_{\substack{q=j-1 \\ q \neq j+l}}^{j+1} \frac{y-y_q}{y_{j+l}-y_q}. $$

Therefore, we use the approximations $u_x \approx \partial_x L_{ij}^{\bar{u}}$ and $u_y \approx \partial_y L_{ij}^{\bar{u}}$. The computation of these approximations of the partial derivatives of the solution is straightforward, since it just involves the computation of derivatives of one dimensional Lagrange polynomial basis.

As a summary, the line integrals in equation (19) are approximated as

$$ \int_{\Gamma_{i \pm 1/2, j}} g_1(u_x, u_y) \mathrm{d}\gamma_j \approx \int_{\Gamma_{i \pm 1/2, j}} g_1(\partial_x L_{ij}^{\bar{u}}, \partial_y L_{ij}^{\bar{u}}) \mathrm{d}\gamma_j \approx \\ \pm \Delta y g_1(\partial_x L_{ij}^{\bar{u}}(x_{i \pm 1/2}, y_j), \partial_y L_{ij}^{\bar{u}}(x_{i \pm 1/2}, y_j)), $$

$$ \int_{\Gamma_{i,j \pm 1/2}} g_2(u_x, u_y) \mathrm{d}\gamma_i \approx \int_{\Gamma_{i,j \pm 1/2}} g_2(\partial_x L_{ij}^{\bar{u}}, \partial_y L_{ij}^{\bar{u}}) \mathrm{d}\gamma_i \approx \\ \pm \Delta x g_2(\partial_x L_{ij}^{\bar{u}}(x_i, y_{j \pm 1/2}), \partial_y L_{ij}^{\bar{u}}(x_i, y_{j \pm 1/2})). $$
### 3.2. IMEX Runge-Kutta time discretization

An IMEX Runge-Kutta scheme consists of applying an implicit discretization to the diffusion terms (stiff terms) and an explicit one to the convective and source terms (non stiff terms). When applied to system (13) at second order it takes the form

$$U^{(k)} = U^n - \Delta t \sum_{l=1}^{k-1} \tilde{a}_{kl} E(U^{(l)}) + \Delta t \sum_{l=1}^{2} a_{kl} I(U^{(l)}), \quad (20)$$

$$U^{n+1} = U^n - \Delta t \sum_{k=1}^{2} \tilde{\omega}_k E(U^{(k)}) + \Delta t \sum_{k=1}^{2} \omega_k I(U^{(k)}), \quad (21)$$

where $U^n = (\bar{u}_{ij}^n)$, $U^{n+1} = (\bar{u}_{ij}^{n+1})$ are the vectors of unknowns volume averages at times $t^n$ and $t^{n+1}$, and $U^{(k)}$ and $U^{(l)}$ are the vector of unknowns at the stages $k, l$ of the IMEX Runge-Kutta scheme. The matrices $\bar{A} = (\tilde{a}_{ij})$, $\tilde{a}_{kl} = 0$ for $l \ge k$ and $A = (a_{kl})$ are $2 \times 2$ matrices such that the scheme is explicit in $E$ and implicit in $I$. The coefficient vectors $\tilde{w} = (\tilde{w}_1, \tilde{w}_2)$ and $w = (w_1, w_2)$ complete the IMEX Runge-Kutta scheme.

One should consider diagonally implicit Runge-Kutta (DIRK) schemes [37] for the implicit terms ($a_{kl} = 0$, for $l > k$) in order to solve efficiently the algebraic equations corresponding to the implicit part of the discretization at each time step. Besides, the use of a DIRK scheme for the implicit part $I$ guarantees that $F$ is always evaluated explicitly.

IMEX Runge-Kutta schemes can be represented by a double tableau in the usual Butcher notation,

$$\begin{array}{c|c|c} \tilde{c} & \tilde{A} & c \\ \hline \tilde{\omega} & & \omega \end{array}$$

where the coefficients vectors $\tilde{c} = (\tilde{c}_1, \tilde{c}_2)^T$ and $c = (c_1, c_2)^T$ used for the treatment of non autonomous systems, are given by the usual relation

$$\tilde{c}_k = \sum_{l=1}^{k-1} \tilde{a}_{kl}, \quad c_k = \sum_{l=1}^{k} a_{kl}. \quad (22)$$
The PDEs (5) and (9) we will be dealing with in this article, after the space discretization, give rise to a system of non-autonomous ODEs. Therefore, the constants $\tilde{c}_k$ and $c_k$ were omitted in the time marching scheme (20)-(21).

In this work we will consider the second order IMEX-SSP2(2,2,2) (L-stable scheme, see [28]), whose *tableaus* for the explicit (left) and implicit (right) parts are

$$ \begin{array}{c|cc|c|cc} 0 & 0 & 0 & \gamma & \gamma & 0 \\ 1 & 1 & 0 & 1-\gamma & 1-2\gamma & \gamma \\ \hline & 1/2 & 1/2 & & 1/2 & 1/2 \end{array} \qquad \gamma = 1 - \frac{1}{\sqrt{2}} $$

Typically, an explicit time integrator needs extremely small time steps due to quadratic CFL-like condition:

$$ 2\eta_1 \frac{\Delta t}{(\Delta x)^2} + 2\eta_2 \frac{\Delta t}{(\Delta y)^2} + \frac{1}{2}\eta_3 \frac{\Delta t}{\Delta x \Delta y} \le 1, \qquad (23) $$

$$ \alpha_1 \frac{\Delta t}{\Delta x} + \alpha_2 \frac{\Delta t}{\Delta y} \le 1, \qquad (24) $$

where

$$ \begin{aligned} \eta_1 &= \left| \frac{\partial g_1}{\partial u_x} \right|, \eta_2 = \left| \frac{\partial g_2}{\partial u_y} \right|, \eta_3 = \left| \frac{\partial g_1}{\partial u_y} + \frac{\partial g_2}{\partial u_x} \right|, \\ \alpha_1 &= \left| \frac{\partial f_1}{\partial u} \right|, \alpha_2 = \left| \frac{\partial f_2}{\partial u} \right|, \end{aligned} $$

for all volumes $V_{ij}$ and $\forall (x, y) \in \Gamma_{ij}$. However, IMEX only needs to satisfy the advection stability condition (24).

## 4. Reference solutions. Fourier Cosine Method

To the best of our knowledge, the most accurate alternative method to assess on the accuracy of our developed finite volume scheme is the COS method proposed by Fang and Oosterlee in [38]. This pricing method is based on Fourier-cosine series expansions. The keystone of Fourier methods is the fact that while probability densities are usually not known for many relevant asset processes, their Fourier transforms, known as the characteristic functions, are available in many relevant cases. Let us describe the method in the two dimensional case.
In this section $t \in [0, T)$ denotes the backward time, so we are interested in the price of the options before the maturity $T$ of the option. Let $\mathbf{x} = (x_1, x_2) \in \mathbb{R}^2$. Consider a two-dimensional pricing model whose risky assets $\mathbf{s} = (s_1, s_2)$ under a selected martingale measure $Q$ take the form

$$ \mathbf{s}_T = (s_{1T}, s_{2T}) = (e^{x_{1T}+x_{10}}, e^{x_{2T}+x_{20}}) = (s_{10}e^{x_{1T}}, s_{20}e^{x_{2T}}), $$

where for $i = 1, 2$, $x_{i0} = \ln(s_{i0})$, $s_{i0} > 0$, and $\mathbf{x}_T = (x_{1T}, x_{2T})$ is a random variable whose characteristic function $\varphi_{\mathbf{x}_T}$ is known. Let $\mathbf{x}_0 = (x_{10}, x_{20})$. The value at time $t$ of an European option with maturity $T$ and payoff function $u_T$ is given by

$$ u(\mathbf{s}_0, t) = e^{-r(T-t)} \mathbb{E}^Q [u_T(\ln(\mathbf{s}_T))] = e^{-r(T-t)} \int_{\mathbb{R}^2} u_T(\mathbf{x} + \mathbf{x}_0) \Phi_{\mathbf{x}_T}(\mathbf{x}) d\mathbf{x}, $$

where $r$ is the interest rate and $\Phi_{\mathbf{x}_T}$ is the probability density function of $\mathbf{x}_T = (x_{1T}, x_{2T})$. Since the tales of $\Phi_{\mathbf{x}_T}$ rapidly decay to zero, we can truncate the infinite integration range $\mathbb{R}^2$ to some domain $I = [a_1, b_1] \times [a_2, b_2] \subset \mathbb{R}^2$ without losing significant accuracy (see [38, 39] for details on how to reduce the computational domain), thus

$$ u(\mathbf{s}_0, t) \approx e^{-r(T-t)} \int_I u_T(\mathbf{x} + \mathbf{x}_0) \Phi_{\mathbf{x}_T}(\mathbf{x}) d\mathbf{x}. $$

At this point we use the Fourier-cosine approximation of $\Phi_{\mathbf{x}_T}$. Let $\mathbf{k} = (k_1, k_2) \in \mathbb{N}^2$, and for $i = 1, 2$ let

$$ \phi_{k_i}^{(i)}(x_i) = \begin{cases} 1, & \text{if } k_i = 0, \\ \sqrt{2} \cos \left( k_i \pi \frac{x_i - a_i}{b_i - a_i} \right), & \text{if } k_i > 0. \end{cases} $$

The cosine basis functions $\phi_k(\mathbf{x}) = \phi_{k_1}^{(1)}(x_1)\phi_{k_2}^{(2)}(x_2)$, are orthonormal in $I$ with respect to the probability density weight $\nu$ of the continuous uniform distribution in $I$, i.e $\nu(\mathbf{x}) = \nu_1(x_1)\nu_2(x_2)$, where for $i = 1, 2$

$$ \nu_i(x_i) = \begin{cases} \frac{1}{b_i - a_i}, & \text{if } a_i \le x_i \le b_i, \\ 0, & \text{if } x_i < a_i \text{ or } x_i > b_i. \end{cases} $$
If $\Phi_{x_T} \in L^2_\nu(I)$, it can be decomposed in its Fourier-cosine series

$$ \Phi_{x_T}(\mathbf{x}) = \sum_{\mathbf{k} \in \mathbb{N}^2} \alpha_{\mathbf{k}} \phi_{\mathbf{k}}(\mathbf{x}), $$

where the Fourier coefficients $\alpha_k$ are given by

$$ \alpha_k = \int_I \Phi_{x_T}(\mathbf{x}) \phi_k(\mathbf{x}) \nu(\mathbf{x}) d\mathbf{x} = \left( \prod_{i=1}^{2} \frac{1}{b_i - a_i} \right) \int_I \Phi_{x_T}(\mathbf{x}) \phi_k(\mathbf{x}) d\mathbf{x}. $$

Let $A_k = \int_I \Phi_{x_T}(\mathbf{x}) \phi_k(\mathbf{x}) d\mathbf{x}$. In practice, the Fourier coefficients are computed only on a finite set of multi-indices $\Gamma \subset \mathbb{N}^d$, therefore $\Phi_{x_T}(\mathbf{x}) \approx \sum_{\mathbf{k} \in \Gamma} \alpha_k \phi_k(\mathbf{x})$. Coming back to the option price, we have that

$$ \begin{aligned} u(s_0, t) &\approx e^{-r(T-t)} \sum_{\mathbf{k} \in \Gamma} \alpha_{\mathbf{k}} \int_I u_T(\mathbf{x} + \mathbf{x}_0) \phi_{\mathbf{k}}(\mathbf{x}) d\mathbf{x} \\ &= e^{-r(T-t)} \sum_{\mathbf{k} \in \Gamma} A_{\mathbf{k}} \int_I u_T(\mathbf{x} + \mathbf{x}_0) \phi_{\mathbf{k}}(\mathbf{x}) \nu(\mathbf{x}) d\mathbf{x}. \end{aligned} $$

Note that

$$ \beta_k(\mathbf{x}_0) = \int_I u_T(\mathbf{x} + \mathbf{x}_0) \phi_k(\mathbf{x}) \nu(\mathbf{x}) d\mathbf{x} = \left( \prod_{i=1}^{2} \frac{1}{b_i - a_i} \right) \int_I u_T(\mathbf{x} + \mathbf{x}_0) \phi_k(\mathbf{x}) d\mathbf{x}, $$

are the Fourier coefficients of the Fourier-cosine expansion of the payoff function $u_T$. Let $B_k(\mathbf{x}_0) = \int_I u_T(\mathbf{x} + \mathbf{x}_0) \phi_k(\mathbf{x}) d\mathbf{x}$. All in all we have

$$ u(s_0, t) \approx e^{-r(T-t)} \left( \prod_{i=1}^{2} \frac{1}{b_i - a_i} \right) \sum_{\mathbf{k} \in \Gamma} A_{\mathbf{k}} B_{\mathbf{k}}(\mathbf{x}_0). $$

The computation of the coefficients $B_k(\mathbf{x}_0)$ depends on the payoff function $u_T$ and only for a few of these financial functions the integral can be made analytically. For the here discussed arithmetic basket options with payoff function (6), where

$$ u_T(x_1, x_2) = \max \left( \frac{1}{2}(s_{10}, s_{20}) \cdot (e^{x_1}, e^{x_2}), 0 \right), $$

the payoff Fourier coefficients have to be computed numerically. Let us now discuss how to compute the coefficients $A_k$ efficiently. For the sake of brevity we restrict the notation to the two dimensional case. In order to alleviate the
notation, let us denote $\hat{k}_i = \frac{\pi}{b_i - a_i} k_i \in \mathbb{R}$ for $i = 1, 2$. In this case,

$$ \phi_k(\mathbf{x}) = \begin{cases} 1, & \text{if } \mathbf{k} = (0,0), \\ \sqrt{2} \cos(\hat{k}_1(x_1 - a_1)), & \text{if } \mathbf{k} = (k_1, 0), k_1 > 0, \\ \sqrt{2} \cos(\hat{k}_2(x_2 - a_2)), & \text{if } \mathbf{k} = (0, k_2), k_2 > 0, \\ 2 \cos(\hat{k}_1(x_1 - a_1)) \cos(\hat{k}_2(x_2 - a_2)), & \text{if } \mathbf{k} = (k_1, k_2), k_1, k_2 > 0. \end{cases} $$

Since

$$ \cos(\hat{k}_i(x_i - a_i)) = \operatorname{Re}(e^{-i a_i \hat{k}_i} e^{i x_i \hat{k}_i}), \text{ and} $$

$$ 2 \cos(\hat{k}_1(x_1 - a_1)) \cos(\hat{k}_2(x_2 - a_2)) = $$

$$ \cos(\hat{k}_1(x_1 - a_1) + \hat{k}_2(x_2 - a_2)) + \cos(\hat{k}_1(x_1 - a_1) - \hat{k}_2(x_2 - a_2)) = $$

$$ \operatorname{Re}(e^{-i(a_1,a_2)\cdot(\hat{k}_1,\hat{k}_2)} e^{i(x_1,x_2)\cdot(\hat{k}_1,\hat{k}_2)}) + \operatorname{Re}(e^{i(a_1,a_2)\cdot(-\hat{k}_1,\hat{k}_2)} e^{i(x_1,x_2)\cdot(\hat{k}_1,-\hat{k}_2)}), $$

and

$$ \int_I \Phi_{x_T}(\mathbf{x}) e^{i(x_1,x_2)\cdot(\hat{k}_1,\hat{k}_2)} d\mathbf{x} \approx \int_{\mathbb{R}^2} \Phi_{x_T}(\mathbf{x}) e^{i(x_1,x_2)\cdot(\hat{k}_1,\hat{k}_1)} d\mathbf{x} = \varphi_{x_T}(\hat{k}_1, \hat{k}_1), $$

it is clear that

$$ A_{(k_1, k_2)} \approx \begin{cases} \operatorname{Re}(\varphi_{x_T}(0,0)), & \text{if } \mathbf{k} = (0,0), \\ \sqrt{2} \operatorname{Re}(e^{-i a_1 \hat{k}_1} \varphi_{x_T}(\hat{k}_1,0)), & \text{if } \mathbf{k} = (k_1,0), k_1 > 0, \\ \sqrt{2} \operatorname{Re}(e^{-i a_2 \hat{k}_2} \varphi_{x_T}(0,\hat{k}_2)), & \text{if } \mathbf{k} = (0,k_2), k_2 > 0, \\ \operatorname{Re}(e^{-i(a_1,a_2)\cdot(\hat{k}_1,\hat{k}_2)} \varphi_{x_T}(\hat{k}_1,\hat{k}_2)) + \\ \operatorname{Re}(e^{i(a_1,a_2)\cdot(-\hat{k}_1,\hat{k}_2)} \varphi_{x_T}(\hat{k}_1,-\hat{k}_2)), & \text{if } \mathbf{k} = (k_1,k_2), k_1, k_2 > 0. \end{cases} \tag{25} $$

This spectral method suffers the curse of dimensionality: not only in the computation of the multidimensional integrals defining $B_k(\mathbf{x}_0)$, but also in the series expansion itself, especially if $\Phi_{x_T}$ is a heavy tailed density. A smart choice of the multi-indices coefficients $\Gamma$ will help to mitigate the curse of dimensionality of this algorithm.

In the Black-Scholes model (2) we have

$$ \varphi_{x_T}(\mathbf{u}) = \exp \left( iT\mu\mathbf{u}^\top - \frac{1}{2}T\mathbf{u}\Sigma\mathbf{u}^\top \right), \tag{26} $$
where $\boldsymbol{\mu} = (\mu_1, \mu_2)$ with $\mu_i = (r - q_i - \frac{1}{2}\sigma_i^2)$ and $\Sigma$ is the covariance matrix with elements $\Sigma_{i,j} = \sigma_i\sigma_j\rho_{ij}$. Let $c_1^{(i)} = T(r - \frac{1}{2}\sigma_i^2)$ and $c_2^{(i)} = T\sigma_i^2$ be the cumulants of the Brownian motion. Here, the way to reduce the computational domain is to consider (see [38])

$$a_i = c_1^{(i)} - L\sqrt{c_2^{(i)}}, \quad b_i = c_1^{(i)} + L\sqrt{c_2^{(i)}}, \quad L \in \mathbb{R}_{+}. \quad (27)$$

At this point let us discuss the computation of the coefficients $A_k$ in the general case $(k_1, k_2 > 0)$ for the Black-Scholes model. By plugging (26) into (25) considering (27), since $\frac{T\mu_i - a_i}{b_i - a_i} = \frac{1}{2}$, one readily obtains

$$A_{(k_1,k_2)} \approx \\ \mathrm{Re} \left( e^{-i(a_1,a_2) \cdot (\hat{k}_1, \hat{k}_2)} \varphi_{x_T} (\hat{k}_1, \hat{k}_2) \right) + \mathrm{Re} \left( e^{i(a_1,a_2) \cdot (-\hat{k}_1, \hat{k}_2)} \varphi_{x_T} (\hat{k}_1, -\hat{k}_2) \right) \\ = C_1 \mathrm{Re} \left( e^{i(T\mu_1 - a_1, T\mu_2 - a_2) \cdot (\hat{k}_1, \hat{k}_2)} \right) + C_2 \mathrm{Re} \left( e^{i(T\mu_1 - a_1, T\mu_2 - a_2) \cdot (\hat{k}_1, -\hat{k}_2)} \right) \\ = C_1 \cos \left( (k_1 + k_2) \frac{\pi}{2} \right) + C_2 \cos \left( (k_1 - k_2) \frac{\pi}{2} \right),$$

where $C_1, C_2 \in \mathbb{R}$ are given by

$$C_1 = \exp \left( -\frac{1}{2} T \left( \hat{k}_1, \hat{k}_2 \right) \Sigma \left( \hat{k}_1, \hat{k}_2 \right)^{\mathrm{T}} \right), \\ C_2 = \exp \left( -\frac{1}{2} T \left( \hat{k}_1, -\hat{k}_2 \right) \Sigma \left( \hat{k}_1, -\hat{k}_2 \right)^{\mathrm{T}} \right).$$

As a consequence, for the Black-Scholes model $A_{(k_1,k_2)} \neq 0$ if and only if $k_1 + k_2$ is even. Therefore, the here considered set of multi-indices is

$$\Gamma^{BS} = \{ (k_1, k_2) \in \mathbb{N}^2 : k_1 \le K_1, k_2 \le K_2, k_1 + k_2 \text{ even} \}.$$

This choice will allow us to overcome the curse of dimensionality in the series expansion. Thus, we will be able to compute, at the same cost, more accurate option prices. Regarding the numerical computation of the payoff coefficients $B_k(\mathbf{x}_0)$, we have used composed Gauss-Legendre quadrature rule with two thousand points on each spatial direction. In order to speedup the computations, we have taken advantage of high performance computing techniques. The loop for the multi-indices $k \in \Gamma$ was parallelized in CPU using four OpenMP threads.
Each one of those threads, is responsible for computing $A_k$ and $B_k(\mathbf{x}_0)$. The calculation of $B_k(\mathbf{x}_0)$ is very costly. In order to accelerate this computation, we have parallelized the quadrature rule using Nvidia CUDA code. This code was ran on four Nvidia Tesla V100.

<caption>**Algorithm 1** Sketch of parallel COS method</caption>

**Require:** $s_0 > 0$
1: $u(\mathbf{s}_0, t) = 0$
2: **parfor** $\mathbf{k} \in \Gamma^{BS}$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;This is a parallel OpenMP for loop
4: &nbsp;&nbsp;&nbsp;&nbsp;Compute $A_k$
5: &nbsp;&nbsp;&nbsp;&nbsp;Compute $B_k(\mathbf{x}_0)$
6: &nbsp;&nbsp;&nbsp;&nbsp;GPU kernel in NVIDIA V100
7: &nbsp;&nbsp;&nbsp;&nbsp;$u(\mathbf{s}_0, t) = A_k B_k(\mathbf{x}_0)$
8: &nbsp;&nbsp;&nbsp;&nbsp;Reduction variable
9: **end parfor**
10: $u(\mathbf{s}_0, t) = e^{-r(T-t)} \left( \prod_{i=1}^{2} \frac{1}{b_i - a_i} \right) u(\mathbf{s}_0, t)$

Finally, in order to price options under the Heston model using the COS method, several approaches are available. For example, one could use the two dimensional COS here presented, using the bivariate characteristic function presented in [39]. Nevertheless, in the Fourier framework, Heston problem can be reduced to a one-dimensional problem. Here we have considered the following characteristic function

$$ \varphi_{x_T}(\xi) = \exp\left(\mathbf{i}\xi(r-q)T + \frac{\nu_0}{\sigma^2}\left(\frac{1-e^{-D(\xi)T}}{1-G(\xi)e^{-D(\xi)T}}\right)(\kappa - \mathbf{i}\rho\sigma\xi - D(\xi))\right) \times \\ \exp\left(\frac{\kappa\theta}{\sigma^2}\left(T(\kappa - \mathbf{i}\rho\sigma\xi - D(\xi)) - 2\ln\left(\frac{1-G(\xi)e^{-D(\xi)T}}{1-G(\xi)}\right)\right)\right), $$

where $\nu_0$ is the initial variance and $D$, $G$ the following functions

$$ D(\xi) = \sqrt{(\kappa - \mathbf{i}\rho\sigma\xi)^2 + (\xi + \mathbf{i})\xi\sigma^2}, \quad G(\xi) = \frac{\kappa - \mathbf{i}\rho\sigma\xi - D(\xi)}{\kappa - \mathbf{i}\rho\sigma\xi + D(\xi)}. $$

Besides, the payoff coefficients for the vanilla call option with payoff $\max(s_T - K, 0)$ can be computed analytically (see [38] for details).
## 5. Numerical experiments

In this section, several numerical experiments are developed to assess the accuracy and performance of the new numerical methods proposed in previous sections for the discussed two dimensional problems.

Experiments are carried out on several option pricing problems. In subsection 5.1 basket call options over two underlyings following Black-Scholes model are priced. Vanilla call options under Heston stochastic volatility model are priced in subsection 5.2. Each subsection is organized as follows.

We start by writing the PDE model in conservative form. Then, we solve the PDE problems using the discussed finite volume IMEX Runge-Kutta numerical method. Both convection-dominated and diffusion-dominated scenarios are considered. As a result, surface prices are presented. In addition, plots for the numerical derivatives of the solution are presented. More precisely, delta and gamma Greeks are shown. All pictures, tables and errors in this whole section will be computed considering the numerical solution at the last time step where $t = T$, i.e, we approximate numerically the option prices “today”.

Later, reference solutions for the studied pricing problems are accurately computed at $t = T$ by means of the COS Fourier method explained in Section 4. They allow us to compute $L_1$ and $L_\infty$ error norms, and the following relative and mean absolute errors given by:

$$L_1 \text{ error} = \Delta x \Delta y \sum_{i,j} |\bar{u}_{i,j} - u(x_i, y_j)|,$$

$$L_\infty \text{ error} = \max_{i,j} |\bar{u}_{i,j} - u(x_i, y_j)|,$$

$$L_\infty \text{ relative error} = \frac{L_\infty \text{ error}}{\max_{i,j} |u(x_i, y_j)|},$$

$$\text{Mean absolute error} = \text{mean}_{i,j} |\bar{u}_{i,j} - u(x_i, y_j)|.$$

Besides, absolute errors of the finite volume IMEX Runge-Kutta method against the COS Fourier method are plotted in order to quickly comprehend the quality of the approximation and the regions with higher errors.
On top of that, orders of convergence are validated in-depth. Both explicit and IMEX Runge-Kutta time integrators were implemented over the same explained space discretization. In the explicit case second order Heun's ODE solver was considered. Finite volume with both time marching schemes achieves the announced second order accuracy in all regimes. IMEX and explicit methods yield similar results in terms of accuracy and convergence order. Nevertheless, the IMEX solver can take advantage of much larger times steps, thus offering much better performance in terms of execution times. The explicit method dramatically slows down the computation due to the presence of diffusive terms.

Finally, each subsection is closed with a table with some reference option prices. Since the here considered option prices are not known is closed form, an interested reader may find this table useful, particularly those developers of alternative option price calculators.

Note that the designed algorithms were implemented using C++ (GNU C++ compiler 9.3.0) and were ran in a machine with 128 GB of RAM, one AMD Ryzen 9 5950X 16-Core Processor, and four GPUs NVIDIA Volta V100. All codes were compiled using double precision arithmetic. Besides, for linear algebra operations Eigen library [40] was considered. For all experiments in this article a CFL of 0.5 is taken into account in the stability conditions.

### 5.1. Options on a basket of two assets

In this experiment we solve the basket European option pricing problem considering two underlying. The PDE model (5) can be written in the conservative form (11) as follows

$$ \frac{\partial u}{\partial t} + \frac{\partial f_1}{\partial s_1}(u) + \frac{\partial f_2}{\partial s_2}(u) = \frac{\partial g_1}{\partial s_1}(u_{s_1}, u_{s_2}) + \frac{\partial g_2}{\partial s_2}(u_{s_1}, u_{s_2}) + h(u), \quad (28) $$
where the functions $f_1, f_2, g_1, g_2$ and $h$ are given by:

$$f_1(u) = (\sigma_1^2 - r + q_1)s_1 u(s_1, s_2, t) + \frac{\rho}{2}\sigma_1\sigma_2 s_1 u(s_1, s_2, t),$$

$$f_2(u) = (\sigma_2^2 - r + q_2)s_2 u(s_1, s_2, t) + \frac{\rho}{2}\sigma_1\sigma_2 s_2 u(s_1, s_2, t),$$

$$g_1(u_{s_1}, u_{s_2}) = \frac{1}{2}\sigma_1^2 s_1^2 u_{s_1}(s_1, s_2, t) + \frac{\rho}{2}\sigma_1\sigma_2 s_1 s_2 u_{s_2}(s_1, s_2, t),$$

$$g_2(u_{s_1}, u_{s_2}) = \frac{1}{2}\sigma_2^2 s_2^2 u_{s_2}(s_1, s_2, t) + \frac{\rho}{2}\sigma_1\sigma_2 s_1 s_2 u_{s_1}(s_1, s_2, t),$$

$$h(u) = (\sigma_2^2 + \rho\sigma_1\sigma_2 + \sigma_1^2 + q_1 + q_2 - 3r)u(s_1, s_2, t).$$

The initial condition $u_0$, being equation (6) for the arithmetic basket call option, was averaged on each volume. Indeed, the initial volume averages are set as

$$\bar{u}_{i,j}^0 = \frac{1}{|V_{ij}|} \int_{V_{ij}} u_0(s_1, s_2) ds_2 ds_1,$$

where the mid-point quadrature formula is used to approximate the integral.

In this section we perform two basket option pricing tests with the market parameters of Table 1. The first set of parameters, denoted as Test 1, stands for a convection-dominated regime. Although nowadays this setup with an interest rate $r = 0.5$ is financially unrealistic, it is useful as a stress-test. The second group of parameters, labelled as Test 2, represents a diffusion-dominated setting. In both scenarios, in order to avoid noise coming from the artificial boundary conditions, the computational space domain was considered as $(s_1, s_2) \in [0, 5K] \times [0, 5K]$.

<table><caption>Table 1: Market data for basket options under Black-Scholes model.</caption><thead><tr><th></th><th>σ<sub>1</sub></th><th>σ<sub>2</sub></th><th>r</th><th>q<sub>1</sub></th><th>q<sub>2</sub></th><th>ρ</th><th>T</th><th>K</th></tr></thead><tbody><tr><td>Test 1</td><td>0.1</td><td>0.1</td><td>0.5</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Test 2</td><td>0.5</td><td>0.5</td><td>0.1</td><td>0.00</td><td>0.00</td><td>0.5</td><td>0.25</td><td>30</td></tr></tbody></table>

Surfaces of basket option prices at $t = T$ for Test 1 and Test 2 are shown in Figure 2. Cuts of these surfaces with planes parallel to $s_2 = 0$ are shown in Figure 3, along with COS solutions. All these plots were drawn using the solution obtained with the finite volume IMEX Runge-Kutta solver. Moreover,
first numerical derivative (delta) and second numerical derivative (gamma) of
the numerical approximation of the option prices are shown in Figures 4 and
5, respectively. The finite volume IMEX Runge-Kutta numerical scheme offers
high resolution approximation of basket option values. Furthermore, it provides
accurate and non-oscillatory approximations of the Greeks.

<figure><img src="image_3.png" /><figcaption>Figure 2: Basket option prices for Test 1 (left) and Test 2 (right) at $t = T$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Cuts of price surfaces of Figure 2. Numerical solution is shown in continuous line and COS solution with squares.</figcaption></figure>

At this point we compute the reference option prices by using the described
COS method. For a mesh of size $N_1 \times N_2 = 1600 \times 1600$, i.e., a mesh with 1600
discretization points in each space direction, the errors of the approximation
<figure><img src="image_5.png" /><figcaption>Figure 4: Basket option deltas for Test 1 (left) and Test 2 (right) at $t = T$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Basket option gammas for Test 1 (left) and Test 2 (right) at $t = T$.</figcaption></figure>

obtained with the finite volume IMEX Runge-Kutta method are collected in Table 2. Figure 6 presents contour plots for the surface defined by the absolute errors in the option prices at $t = T$. The here developed finite volumes (explicit or IMEX) schemes offer high-resolution approximations even at regions of discontinuities and non-smoothness in the initial condition. Therefore, there is no need to apply smoothing methods like the Rannacher technique.

Tables 3 and 4 record $L_1$ errors and $L_1$ orders of convergence for Test 1 and Test 2, respectively. The errors and orders of convergence are shown for both the IMEX and explicit finite volume numerical schemes. Both schemes
<table><caption>Table 2: Finite volume IMEX Runge-Kutta numerical errors against reference basket option prices computed with COS method at $t = T$.</caption><thead><tr><th colspan="3">Test 1</th></tr><tr><th>L<sub>∞</sub> error</th><th>L<sub>∞</sub> relative error</th><th>Mean absolute error</th></tr></thead><tbody><tr><td>1.0245 × 10<sup>-4</sup></td><td>1.3976 × 10<sup>-6</sup></td><td>3.5079 × 10<sup>-6</sup></td></tr></tbody><thead><tr><th colspan="3">Test 2</th></tr><tr><th>L<sub>∞</sub> error</th><th>L<sub>∞</sub> relative error</th><th>Mean absolute error</th></tr></thead><tbody><tr><td>2.3406 × 10<sup>-5</sup></td><td>1.9439 × 10<sup>-7</sup></td><td>2.7099 × 10<sup>-6</sup></td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 6: Contour plots for the surface of absolute errors of basket options prices for Test 1 (left) and Test 2 (right) at $t = T$.</figcaption></figure>

achieve second-order accuracy in the $L_1$ norm. Additionally, these tables present the time steps and the execution times for the two time integrators. As said before, IMEX method is able to converge using much larger time steps than the explicit one. The stability condition of the explicit scheme requires extremely small time steps, thus making the method useless in practice for refined grids in space. Therefore, IMEX offers much better performance in terms of execution times, and allows us to solve the PDE problems with refined meshes in space. In fact, IMEX method is able to run for grids finer than 800 × 800 grid, while the explicit method is not in reasonable computational times. In Figure 7, the natural logarithms of $L_1$ errors and execution times of Table 4 are shown for both the IMEX and explicit numerical schemes. IMEX clearly outperforms the
explicit method in this scenario with diffusion dominance, which is the typical situation in finance. Although for convective dominated problems in coarse grids both time marching schemes perform similarly, as soon as the mesh is refined IMEX is the only practical choice.

<table><caption>Table 3: $L_1$ errors and $L_1$ orders of convergence of the IMEX and explicit finite volume methods for Test 1.</caption><thead><tr><th colspan="5">IMEX</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>9.9867 × 10<sup>1</sup></td><td>--</td><td>2.05 × 10<sup>-2</sup></td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>50 × 50</td><td>3.3457 × 10<sup>1</sup></td><td>1.57</td><td>1.03 × 10<sup>-2</sup></td><td>3.9 × 10<sup>-2</sup></td></tr><tr><td>100 × 100</td><td>9.1341 × 10<sup>0</sup></td><td>1.87</td><td>5.13 × 10<sup>-3</sup></td><td>2.3 × 10<sup>-1</sup></td></tr><tr><td>200 × 200</td><td>2.3529 × 10<sup>0</sup></td><td>1.95</td><td>2.56 × 10<sup>-3</sup></td><td>1.8 × 10<sup>0</sup></td></tr><tr><td>400 × 400</td><td>5.6234 × 10<sup>-1</sup></td><td>2.06</td><td>1.28 × 10<sup>-3</sup></td><td>1.6 × 10<sup>1</sup></td></tr><tr><td>800 × 800</td><td>1.1257 × 10<sup>-1</sup></td><td>2.32</td><td>6.41 × 10<sup>-4</sup></td><td>1.5 × 10<sup>2</sup></td></tr></tbody></table><table><thead><tr><th colspan="5">Explicit</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>9.9752 × 10<sup>1</sup></td><td>--</td><td>3.56 × 10<sup>-2</sup></td><td>9.4 × 10<sup>-4</sup></td></tr><tr><td>50 × 50</td><td>3.3453 × 10<sup>1</sup></td><td>1.57</td><td>8.88 × 10<sup>-3</sup></td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>100 × 100</td><td>9.1503 × 10<sup>0</sup></td><td>1.87</td><td>2.22 × 10<sup>-3</sup></td><td>8.3 × 10<sup>-2</sup></td></tr><tr><td>200 × 200</td><td>2.3614 × 10<sup>0</sup></td><td>1.95</td><td>5.56 × 10<sup>-4</sup></td><td>1.4 × 10<sup>0</sup></td></tr><tr><td>400 × 400</td><td>5.6466 × 10<sup>-1</sup></td><td>2.06</td><td>1.39 × 10<sup>-4</sup></td><td>2.1 × 10<sup>1</sup></td></tr><tr><td>800 × 800</td><td>1.1304 × 10<sup>-1</sup></td><td>2.32</td><td>3.47 × 10<sup>-5</sup></td><td>3.6 × 10<sup>2</sup></td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Figure 7: Efficiency curve of IMEX and explicit time marching schemes for Test 1.</figcaption></figure>

Finally, Tables 5 and 6 show some basket option prices at $t = T$, for Test 1 and Test 2, respectively. These prices were obtained with COS method and
<table><thead><tr><th colspan="5">IMEX</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>9.6620 × 10<sup>1</sup></td><td>--</td><td>4.71 × 10<sup>-2</sup></td><td>4.5 × 10<sup>-3</sup></td></tr><tr><td>50 × 50</td><td>2.5178 × 10<sup>1</sup></td><td>1.94</td><td>2.35 × 10<sup>-2</sup></td><td>3.3 × 10<sup>-2</sup></td></tr><tr><td>100 × 100</td><td>6.4828 × 10<sup>0</sup></td><td>1.95</td><td>1.18 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-1</sup></td></tr><tr><td>200 × 200</td><td>1.6209 × 10<sup>0</sup></td><td>2.00</td><td>5.89 × 10<sup>-3</sup></td><td>1.2 × 10<sup>0</sup></td></tr><tr><td>400 × 400</td><td>3.9419 × 10<sup>-1</sup></td><td>2.03</td><td>2.94 × 10<sup>-3</sup></td><td>9.8 × 10<sup>0</sup></td></tr><tr><td>800 × 800</td><td>7.9229 × 10<sup>-2</sup></td><td>2.31</td><td>1.47 × 10<sup>-3</sup></td><td>8.5 × 10<sup>1</sup></td></tr></tbody></table>

<table><caption>Table 4: $L_1$ errors and $L_1$ orders of convergence of the IMEX and explicit finite volume methods for Test 2.</caption><thead><tr><th colspan="5">Explicit</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>9.0224 × 10<sup>1</sup></td><td>--</td><td>1.42 × 10<sup>-3</sup></td><td>1.7 × 10<sup>-2</sup></td></tr><tr><td>50 × 50</td><td>2.3440 × 10<sup>1</sup></td><td>1.94</td><td>3.56 × 10<sup>-4</sup></td><td>1.2 × 10<sup>-1</sup></td></tr><tr><td>100 × 100</td><td>5.9498 × 10<sup>0</sup></td><td>1.97</td><td>8.89 × 10<sup>-5</sup></td><td>1.8 × 10<sup>0</sup></td></tr><tr><td>200 × 200</td><td>1.4834 × 10<sup>0</sup></td><td>2.00</td><td>2.22 × 10<sup>-5</sup></td><td>3.0 × 10<sup>1</sup></td></tr><tr><td>400 × 400</td><td>3.5473 × 10<sup>-1</sup></td><td>2.06</td><td>5.56 × 10<sup>-6</sup></td><td>4.9 × 10<sup>2</sup></td></tr><tr><td>800 × 800</td><td>7.1095 × 10<sup>-2</sup></td><td>2.31</td><td>1.39 × 10<sup>-6</sup></td><td>7.9 × 10<sup>3</sup></td></tr></tbody></table>

with the finite volume IMEX Runge-Kutta solver executed over a space grid of size 1600 × 1600. In the COS method, in order to truncate the series expansion, we considered $K_1 = K_2 = 51$ and we set $L = 12$ to cut the tales of the unknown joint density of the underlying assets.
<table><caption>Table 5: Comparison between COS and finite volume IMEX Runge-Kutta numerical solutions for Test 1 at $t = T$. At each cell of the table the COS price is the first one, while the finite volume price is the second one.</caption><thead><tr><th>s<sub>1</sub></th><th>s<sub>2</sub></th><th>20.03125</th><th>40.03125</th><th>60.03125</th><th>80.03125</th></tr></thead><tbody><tr><td rowspan="2">20.03125</td><td>1.1495452217 × 10<sup>-11</sup></td><td>3.5571109320</td><td>13.5563429224</td><td>23.5563429224</td></tr><tr><td>6.1707381613 × 10<sup>-12</sup></td><td>3.5571143743</td><td>13.5563429224</td><td>23.5563429224</td></tr><tr><td rowspan="2">40.03125</td><td>3.5571109320</td><td>13.5563429224</td><td>23.5563429224</td><td>33.5563429224</td></tr><tr><td>3.5571143743</td><td>13.5563429224</td><td>23.5563429224</td><td>33.5563429224</td></tr><tr><td rowspan="2">60.03125</td><td>13.5563429224</td><td>23.5563429224</td><td>33.5563429224</td><td>43.5563429224</td></tr><tr><td>13.5563429224</td><td>23.5563429224</td><td>33.5563429224</td><td>43.5563429224</td></tr><tr><td rowspan="2">80.03125</td><td>23.5563429224</td><td>33.5563429224</td><td>43.5563429224</td><td>53.5563429224</td></tr><tr><td>23.5563429224</td><td>33.5563429224</td><td>43.5563429224</td><td>53.5563429224</td></tr></tbody></table>

<table><caption>Table 6: Comparison between COS and finite volume IMEX Runge-Kutta numerical solutions for Test 2 at $t = T$. At each cell of the table the COS price is the first one, while the finite volume price is the second one.</caption><thead><tr><th>s<sub>1</sub></th><th>s<sub>2</sub></th><th>20.296875</th><th>39.046875</th><th>57.796875</th><th>76.546875</th></tr></thead><tbody><tr><td rowspan="2">20.296875</td><td>9.9328285143 × 10<sup>-2</sup></td><td>2.7953041634</td><td>10.1409569225</td><td>19.2023244718</td></tr><tr><td>9.9312989509 × 10<sup>-2</sup></td><td>2.7953229793</td><td>10.1409654135</td><td>19.2023248970</td></tr><tr><td rowspan="2">39.046875</td><td>2.7953041634</td><td>10.0983319903</td><td>19.1914361881</td><td>28.5400105675</td></tr><tr><td>2.7953229793</td><td>10.0983391832</td><td>19.1914362212</td><td>28.5400103410</td></tr><tr><td rowspan="2">57.796875</td><td>10.1409569225</td><td>19.1914361881</td><td>28.5396217167</td><td>37.9127244740</td></tr><tr><td>10.1409654135</td><td>19.1914362212</td><td>28.5396214779</td><td>37.9127244317</td></tr><tr><td rowspan="2">76.546875</td><td>19.2023244718</td><td>28.5400105675</td><td>37.9127244740</td><td>47.2875872432</td></tr><tr><td>19.2023248970</td><td>28.5400103410</td><td>37.9127244317</td><td>47.2875872518</td></tr></tbody></table>
## 5.2. Heston model

In this experiment we price a vanilla call option under the Heston stochastic volatility model. The PDE model (9) can be written in the conservative form (11) as follows:

$$ \frac{\partial u}{\partial t} + \frac{\partial f_1}{\partial s}(u) + \frac{\partial f_2}{\partial v}(u) = \frac{\partial g_1}{\partial s}(u_s, u_v) + \frac{\partial g_2}{\partial v}(u_s, u_v) + h(u), \quad (29) $$

where the functions $f_1, f_2, g_1, g_2$ and $h$ are given by:

$$ f_1(u) = (v - r + q)su(s, v, t), $$

$$ f_2(u) = \left( \rho\sigma v - \kappa(\theta - v) + \frac{1}{2}\sigma^2 \right) u(s, v, t), $$

$$ g_1(u_s, u_v) = \frac{1}{2} s^2 v u_s(s, v, t) + \rho \sigma s v u_v(s, v, t), $$

$$ g_2(u_s, u_v) = \frac{1}{2} \sigma^2 v u_v(s, v, t), $$

$$ h(u) = (v - 2r + q + \kappa + \rho\sigma)u(s, v, t). $$

The initial condition $u_0$, being equation (10) for the European call option, was averaged on each volume. The initial volume averages are set as

$$ \bar{u}_{i,j}^0 = \frac{1}{|V_{ij}|} \int_{V_{ij}} u_0(s) ds dv. $$

In this part we perform two vanilla call option pricing tests with the market parameters of Table 7. The first set of parameters, labelled as Test 3, is taken from [10]. The second group of parameters, denoted as Test 4, is a variation of Test 3, swapping the interest rate $r$ and the volatility of the volatility $\sigma$. In both tests, in order to minimize the numerical errors coming from the artificial boundary conditions, the computational space domain was considered as $(s, v) \in [0, 800] \times [0, 4]$.

Figure 8 displays the finite volume IMEX Runge-Kutta approximations of the exact option price functions $u$ at $t = T$ corresponding to the two cases of Table 7. Cuts of these price surfaces with planes parallel to $v = 0$ are shown in Figure 9 in combination with COS solutions. Additionally, Figures 10 and 11 present the numerical Greeks of the numerical approximation, deltas and
<table><caption>Table 7: Market data for vanilla call options under Heston model.</caption><thead><tr><th></th><th>σ</th><th>r</th><th>q</th><th>κ</th><th>θ</th><th>ρ</th><th>T</th><th>K</th></tr></thead><tbody><tr><td>Test 3</td><td>0.3</td><td>0.025</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Test 4</td><td>0.025</td><td>0.3</td><td>0.0</td><td>1.5</td><td>0.04</td><td>-0.9</td><td>0.25</td><td>100</td></tr></tbody></table>

gammas, respectively. Once again, the here proposed numerical scheme provides accurate and non-oscillatory approximations of the Greeks, even at regions of discontinuities and non-smoothness in the initial condition.

<figure><img src="image_9.png" /><figcaption>Figure 8: Vanilla call option prices for Test 3 (left) and Test 4 (right) at $t = T$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Cuts of prices surfaces of Figure 8. Numerical solution is shown in continuous line and COS solution with squares.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Vanilla call option deltas for Test 3 (left) and Test 4 (right) at $t = T$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11: Vanilla call option gammas for Test 3 (left) and Test 4 (right) at $t = T$.</figcaption></figure>

At this stage, using the described COS method, we compute reference call option prices. For a mesh of size $3200 \times 3200$, the errors of the approximation obtained with the finite volume IMEX Runge-Kutta method are presented in Table 8. Figure 8 displays contour plots for the surface defined by the absolute errors in the call option prices at $t = T$. These plots show that the here developed finite volume schemes provide high-quality approximations even at the region of non-smoothness in the payoff.

In Tables 9 and 10, the $L_1$ errors and the $L_1$ orders of the explicit and IMEX methods for the Test 3 and Test 4 are displayed, respectively. Again, both
<table><caption>Table 8: Finite volume IMEX Runge-Kutta numerical errors against reference call option prices computed with COS method at $t = T$.</caption><thead><tr><th colspan="3">Test 3</th></tr><tr><th>L<sub>∞</sub> error</th><th>L<sub>∞</sub> relative error</th><th>Mean absolute error</th></tr></thead><tbody><tr><td>2.0377 × 10<sup>-4</sup></td><td>3.2852 × 10<sup>-7</sup></td><td>3.1711 × 10<sup>-5</sup></td></tr></tbody></table><table><thead><tr><th colspan="3">Test 4</th></tr><tr><th>L<sub>∞</sub> error</th><th>L<sub>∞</sub> relative error</th><th>Mean absolute error</th></tr></thead><tbody><tr><td>1.5196 × 10<sup>-4</sup></td><td>1.9969 × 10<sup>-5</sup></td><td>2.6785 × 10<sup>-7</sup></td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 12: Contour plots for the surface of absolute errors of call option prices for Test 3 (left) and Test 4 (right) at $t = T$.</figcaption></figure>

numerical schemes achieve second-order accuracy in the $L_1$ norm. Even if we impose time-dependent boundary conditions, we do not observe the troublesome order reduction effect described in other papers. As said previously, IMEX is able to converge using much larger time steps than the explicit method, thus it consumes much less computing time.

Finally, Tables 11 and 12 show some vanilla call option prices at $t = T$ under Heston stochastic volatility model, for Test 3 and Test 4, respectively. These prices were obtained with COS method and with the finite volume IMEX Runge-Kutta solver executed over a grid with size 3200 × 3200.
<table><thead><tr><th colspan="5">IMEX</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>3.1356 × 10<sup>2</sup></td><td>--</td><td>3.85 × 10<sup>-3</sup></td><td>1.6 × 10<sup>-2</sup></td></tr><tr><td>50 × 50</td><td>1.0516 × 10<sup>2</sup></td><td>1.58</td><td>1.92 × 10<sup>-3</sup></td><td>9.7 × 10<sup>-2</sup></td></tr><tr><td>100 × 100</td><td>2.5715 × 10<sup>1</sup></td><td>2.03</td><td>9.61 × 10<sup>-4</sup></td><td>7.4 × 10<sup>-1</sup></td></tr><tr><td>200 × 200</td><td>6.2028 × 10<sup>0</sup></td><td>2.05</td><td>4.80 × 10<sup>-4</sup></td><td>7.4 × 10<sup>0</sup></td></tr><tr><td>400 × 400</td><td>1.4586 × 10<sup>0</sup></td><td>2.09</td><td>2.40 × 10<sup>-4</sup></td><td>6.8 × 10<sup>1</sup></td></tr><tr><td>800 × 800</td><td>2.9016 × 10<sup>-1</sup></td><td>2.33</td><td>1.20 × 10<sup>-4</sup></td><td>6.4 × 10<sup>2</sup></td></tr></tbody></table>

<table><caption>Table 9: $L_1$ errors and $L_1$ orders of convergence of the IMEX and explicit finite volume methods for Test 3.</caption><thead><tr><th colspan="5">Explicit</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>3.1312 × 10<sup>2</sup></td><td>--</td><td>1.92 × 10<sup>-4</sup></td><td>1.9 × 10<sup>-1</sup></td></tr><tr><td>50 × 50</td><td>1.0489 × 10<sup>2</sup></td><td>1.58</td><td>4.81 × 10<sup>-5</sup></td><td>2.7 × 10<sup>0</sup></td></tr><tr><td>100 × 100</td><td>2.5649 × 10<sup>1</sup></td><td>2.03</td><td>1.20 × 10<sup>-5</sup></td><td>4.4 × 10<sup>1</sup></td></tr><tr><td>200 × 200</td><td>6.1864 × 10<sup>0</sup></td><td>2.05</td><td>3.00 × 10<sup>-6</sup></td><td>7.3 × 10<sup>2</sup></td></tr><tr><td>400 × 400</td><td>1.4543 × 10<sup>0</sup></td><td>2.09</td><td>7.52 × 10<sup>-7</sup></td><td>1.2 × 10<sup>4</sup></td></tr><tr><td>800 × 800</td><td>2.8924 × 10<sup>-1</sup></td><td>2.33</td><td>1.87 × 10<sup>-7</sup></td><td>2.0 × 10<sup>5</sup></td></tr></tbody></table>

<table><thead><tr><th colspan="5">IMEX</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>8.9656 × 10<sup>1</sup></td><td>--</td><td>3.87 × 10<sup>-3</sup></td><td>1.5 × 10<sup>-2</sup></td></tr><tr><td>50 × 50</td><td>2.4203 × 10<sup>1</sup></td><td>1.89</td><td>1.94 × 10<sup>-3</sup></td><td>9.0 × 10<sup>-2</sup></td></tr><tr><td>100 × 100</td><td>9.3022 × 10<sup>0</sup></td><td>1.38</td><td>9.69 × 10<sup>-4</sup></td><td>7.0 × 10<sup>-1</sup></td></tr><tr><td>200 × 200</td><td>1.2578 × 10<sup>0</sup></td><td>2.89</td><td>4.84 × 10<sup>-4</sup></td><td>6.9 × 10<sup>0</sup></td></tr><tr><td>400 × 400</td><td>2.9883 × 10<sup>-1</sup></td><td>2.07</td><td>2.42 × 10<sup>-4</sup></td><td>6.8 × 10<sup>1</sup></td></tr><tr><td>800 × 800</td><td>5.9846 × 10<sup>-2</sup></td><td>2.32</td><td>1.21 × 10<sup>-4</sup></td><td>6.3 × 10<sup>2</sup></td></tr></tbody></table>

<table><caption>Table 10: $L_1$ errors and $L_1$ orders of convergence of the IMEX and explicit finite volume methods for Test 4.</caption><thead><tr><th colspan="5">Explicit</th></tr><tr><th>N<sub>1</sub> × N<sub>2</sub></th><th>L<sub>1</sub> error</th><th>Order</th><th>Δt</th><th>Time (s)</th></tr></thead><tbody><tr><td>25 × 25</td><td>9.0576 × 10<sup>1</sup></td><td>--</td><td>1.99 × 10<sup>-4</sup></td><td>1.8 × 10<sup>-1</sup></td></tr><tr><td>50 × 50</td><td>2.4440 × 10<sup>1</sup></td><td>1.89</td><td>4.99 × 10<sup>-5</sup></td><td>2.7 × 10<sup>0</sup></td></tr><tr><td>100 × 100</td><td>9.3648 × 10<sup>0</sup></td><td>1.38</td><td>1.25 × 10<sup>-5</sup></td><td>4.4 × 10<sup>2</sup></td></tr><tr><td>200 × 200</td><td>1.2671 × 10<sup>0</sup></td><td>2.89</td><td>3.12 × 10<sup>-6</sup></td><td>7.1 × 10<sup>2</sup></td></tr><tr><td>400 × 400</td><td>3.0089 × 10<sup>-1</sup></td><td>2.07</td><td>7.79 × 10<sup>-7</sup></td><td>1.2 × 10<sup>4</sup></td></tr><tr><td>800 × 800</td><td>6.0411 × 10<sup>-2</sup></td><td>2.32</td><td>1.95 × 10<sup>-7</sup></td><td>1.9 × 10<sup>5</sup></td></tr></tbody></table>
<table><caption>Table 11: Comparison between COS and finite volume IMEX Runge-Kutta numerical solutions for Test 3 at $t = T$. At each cell of the table the COS price is the first one, while the finite volume price is the second one.</caption><thead><tr><th>v</th><th>75.125</th><th>100.125</th><th>125.125</th><th>150.125</th></tr></thead><tbody><tr><td rowspan="2">0.200625</td><td>0.4317502658</td><td>8.5901684247</td><td>27.6694906563</td><td>51.1935422402</td></tr><tr><td>0.4316035999</td><td>8.5901562104</td><td>27.6695002425</td><td>51.1935390259</td></tr><tr><td rowspan="2">0.400625</td><td>1.8664224003</td><td>11.8552660514</td><td>30.0081668834</td><td>52.3493028388</td></tr><tr><td>1.8662948639</td><td>11.8552481800</td><td>30.0081846067</td><td>52.3493135590</td></tr><tr><td rowspan="2">0.600625</td><td>3.3658574222</td><td>14.3630707461</td><td>32.1382630288</td><td>53.7378969543</td></tr><tr><td>3.3657449349</td><td>14.3630493039</td><td>32.1382791463</td><td>53.7379138145</td></tr><tr><td rowspan="2">0.800625</td><td>4.7800189282</td><td>16.4716211703</td><td>34.0691648265</td><td>55.1769879733</td></tr><tr><td>4.7799153837</td><td>16.4715967918</td><td>34.0691774343</td><td>55.1770065036</td></tr></tbody></table>

<table><caption>Table 12: Comparison between COS and finite volume IMEX Runge-Kutta numerical solutions for Test 4 at $t = T$. At each cell of the table the COS price is the first one, while the finite volume price is the second one.</caption><thead><tr><th>v</th><th>75.125</th><th>100.125</th><th>125.125</th><th>150.125</th></tr></thead><tbody><tr><td rowspan="2">0.200625</td><td>1.3840492642</td><td>12.2239735395</td><td>33.1373172158</td><td>57.4478630579</td></tr><tr><td>1.3839721771</td><td>12.2239654379</td><td>33.1372986255</td><td>57.4478503229</td></tr><tr><td rowspan="2">0.400625</td><td>3.2984130167</td><td>15.2478827153</td><td>34.8602149676</td><td>58.0809835379</td></tr><tr><td>3.2983708067</td><td>15.2478733955</td><td>34.8602005816</td><td>58.0809679113</td></tr><tr><td rowspan="2">0.600625</td><td>5.0070610101</td><td>17.6182306257</td><td>36.6338455366</td><td>59.0749794899</td></tr><tr><td>5.0070098251</td><td>17.6182226102</td><td>36.6338456927</td><td>59.0749764400</td></tr><tr><td rowspan="2">0.800625</td><td>6.5349321045</td><td>19.6278242193</td><td>38.3205477349</td><td>60.2165220825</td></tr><tr><td>6.5348770744</td><td>19.6278146746</td><td>38.3205529668</td><td>60.2165267905</td></tr></tbody></table>
## 6. Conclusions

We have shown that finite volume IMEX Runge-Kutta numerical schemes are suitable for solving convection-diffusion PDE option pricing problems, even in the presence of mixed derivatives in the diffusive operator. This opens the door to the application of these schemes to numerous models in finance, even those giving rise to non-linear PDEs.

The obtained numerical schemes are highly efficient. On the one hand, large time steps can be used, avoiding the need to use small time steps enforced by the diffusion stability condition that appears when explicit schemes are considered. On the other hand, the schemes are second order accurate. This fact is of paramount importance in order to obtain accurate approximations of the Greeks without oscillations. Besides, we have shown that second order accuracy is preserved even when non smooth initial conditions (payoffs) are considered, which is the usual situation in finance. Thus, additional special techniques, like the Rannacher time stepping, do not need to be taken into account. Moreover, the here developed option price calculators can be extended to build numerical solvers with higher order convergence.

Last, but not least, in this work we provide an alternative way to compute very accurate approximations of the prices of arithmetic basket call options by means of a highly accurate and efficient multidimensional COS Fourier method. These Fourier semi-analytical solutions can be also very valuable benchmarks for a broader audience that work in the development of high order schemes in the general parabolic setting, even outside the financial world.

## References

[1] P. Glasserman, Monte Carlo Methods in Financial Engineering, Springer, New York, 2003.

[2] E. Gobet, J. López-Salas, C. Vázquez, Quasi-Regression Monte-Carlo Scheme for Semi-Linear PDEs and BSDEs with Large Scale Parallelization
on GPUs, Archives of Computational Methods in Engineering 27 (2020)
889–921.

[3] Y. Achdou, O. Pironneau, Frontiers in Applied Mathematics Computational Methods for Option Pricing, SIAM, Philadelphia, 2005.

[4] D. Duffy, Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach, Wiley, New Jersey, 2006.

[5] P. Wilmott, Paul Wilmott on Quantitative Finance, Wiley, West Sussex, 2006.

[6] D. Tavella, C. Randall, Pricing Financial Instruments: The Finite Difference Method, John Wiley & Sons, New York, 2000.

[7] Y. Achdou, O. Pironneau, Volatility smile by multilevel least square, International Journal of Theoretical and Applied Finance 05 (06) (2002) 619–643.

[8] A. Bermúdez, M. R. Nogueiras, C. Vázquez, Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristics/Finite Elements. Part I: Time Discretization, SIAM Journal on Numerical Analysis 44 (5) (2006) 1829–1853.

[9] A. Bermúdez, M. R. Nogueiras, C. Vázquez, Numerical Analysis of Convection-Diffusion-Reaction Problems with Higher Order Characteristics/Finite Elements. Part II: Fully Discretized Scheme and Quadrature Formulas, SIAM Journal on Numerical Analysis 44 (5) (2006) 1854–1876.

[10] K. J. I. T. Hout, S. Foulon, ADI finite difference schemes for option pricing in the Heston model with correlation, International Journal of Numerical Analysis and Modeling 7 (2) (2010) 303–320.

[11] C. Hendricks, C. Heuer, M. Ehrhardt, M. Günther, High-order adi finite difference schemes for parabolic equations in the combination technique with
application in finance, Journal of Computational and Applied Mathematics 316 (2017) 175-194.

[12] S. Heston, G. Zhou, On the Rate of Convergence of Discrete-Time Contingent Claims, Mathematical Finance 10 (2000) 53-75.

[13] R. Rannacher, Finite element solution of diffusion problems with irregular data, Numerische Mathematik 43 (1984) 309-327.

[14] M. Luskin, R. Rannacher, W. Wendland, On the smoothing property of the crank-nicolson scheme, Applicable Analysis 14 (2) (1982) 117-135.

[15] M. Giles, R. Carter, Convergence analysis of Crank-Nicolson and Rannacher time-marching, Journal of Computational Finance 9 (2005).

[16] D. M. Pooley, K. R. Vetzal, P. A. Forsyth, Convergence remedies for non-smooth payoffs in option pricing, Journal of Computational Finance 6 (4) (2002).

[17] Y. d'Halluin, P. A. Forsyth, G. Labahn, A semi-Lagrangian approach for American Asian options under jump diffusion, SIAM Journal on Scientific Computing 27 (1) (2006) 315-345.

[18] D. Castillo, A. Ferreiro, J. García-Rodríguez, C. Vázquez, Numerical methods to solve PDE models for pricing business companies in different regimes and implementation in GPUs, Applied Mathematics and Computation 219 (24) (2013) 11233-11257.

[19] R. Zvan, P. A. Forsyth, K. R. Vetzal, A Finite Volume Approach For Contingent Claims Valuation, IMA Journal of Numerical Analysis 21 (3) (2001) 703-731.

[20] R. Zvan, P. A. Forsyth, K. R. Vetzal, Robust numerical methods for PDE models of Asian options, Journal of computational Finance 1 (2) (1997) 39-78.
[21] G. I. Ramírez-Espinoza, M. Ehrhardt, Conservative and finite volume methods for the convection-dominated pricing problem, Advances in Applied Mathematics and Mechanics 5 (6) (2013) 759–790.

[22] H. Nessyahu, E. Tadmor, Non-oscillatory central differencing for hyperbolic conservation laws, Journal of Computational Physics 87 (2) (1990) 408–463.

[23] A. Kurganov, E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations, Journal of Computational Physics 160 (1) (2000) 241–282.

[24] O. Bhatoo, A. A. I. Peer, E. Tadmor, D. Y. Tangman, A. A. E. F. Saib, Efficient conservative second-order central-upwind schemes for option-pricing problems, Journal of computational Finance 22 (5) (2019) 39–78.

[25] O. Bhatoo, A. A. I. Peer, E. Tadmor, D. Y. Tangman, A. A. E. F. Saib, Conservative Third-Order Central-Upwind Schemes for Option Pricing Problems, Vietnam Journal of Mathematics 47 (4) (2019) 813–8332.

[26] A. Kurganov, D. Levy, A Third-Order Semidiscrete Central Scheme for Conservation Laws and Convection-Diffusion Equations, SIAM Journal on Scientific Computing 22 (4) (2000) 1461–1488.

[27] D. Levy, G. Puppo, G. Russo, Compact central WENO schemes for multidimensional conservation laws, SIAM Journal on Scientific Computing 22 (2) (2000) 656–672.

[28] L. Pareschi, G. Russo, Implicit–Explicit Runge–Kutta Schemes and Applications to Hyperbolic Systems with Relaxation, Journal of Scientific Computing 25 (2005) 129–155.

[29] S. Salmi, J. Toivanen, IMEX schemes for pricing options under jump-diffusion models, Applied Numerical Mathematics 84 (2014) 33–45.

[30] J. de Frutos, Implicit–explicit Runge–Kutta methods for financial derivatives pricing models, European Journal of Operational Research 171 (2006) 991–1004.
[31] F. Black, M. Scholes, The Pricing of Options and Corporate Liabilities, The Journal of Political Economy 81 (3) (1973) 637–654.

[32] R. C. Merton, Theory of rational option pricing, The Bell Journal of Economics and Management Science 4 (1) (1973) 141–183.

[33] J. Gatheral, The Volatility Surface, John Wiley & Sons, New Jersey, 2006.

[34] S. L. Heston, A Closed-Form Solution for Options with Stochastic Volatility with Applications to Bond and Currency Options, The Review of Financial Studies 6 (2) (1993) 327–343.

[35] J. C. Cox, J. E. Ingersoll, S. A. Ross, A theory of the term structure of interest rates, Econometrica 53 (1985) 385–407.

[36] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method, Journal of Computational Physics 32 (1) (1979) 101–136.

[37] E. Hairer, G. Wanner, Solving Ordinary Differential Equations, Vol.2. Stiff and Differential-algebraic Problems, Springer-Verlag, New York, 1987.

[38] F. Fang, C. W. Oosterlee, A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions, SIAM Journal on Scientific Computing 31 (2) (2008) 826–848.

[39] M. J. Ruijter, C. W. Oosterlee, Two-Dimensional Fourier Cosine Series Expansion Method for Pricing Financial Options, SIAM Journal on Scientific Computing 34 (5) (2012) B642–B671.

[40] G. Guennebaud, B. Jacob, et al., Eigen v3.4.0, http://eigen.tuxfamily.org (2010).