Operator splitting schemes for the two-asset Merton jump-diffusion model

TL;DR

This paper develops and compares seven operator splitting schemes for efficiently solving the complex two-asset Merton jump-diffusion PIDE, enabling accurate valuation of rainbow options with nonlocal and mixed derivative features.

Contribution

It introduces and evaluates seven novel operator splitting schemes tailored for the two-asset Merton jump-diffusion PIDE, enhancing numerical stability and efficiency.

Findings

Seven schemes exhibit stable convergence in numerical experiments.

Explicit treatment of the integral part improves computational efficiency.

Performance varies depending on option type and scheme choice.

Abstract

This paper deals with the numerical solution of the two-dimensional time-dependent Merton partial integro-differential equation (PIDE) for the values of rainbow options under the two-asset Merton jump-diffusion model. Key features of this well-known equation are a two-dimensional nonlocal integral part and a mixed spatial derivative term. For its efficient and stable numerical solution, we study seven recent and novel operator splitting schemes of the implicit-explicit (IMEX) and the alternating direction implicit (ADI) kind. Here the integral part is always conveniently treated in an explicit fashion. The convergence behaviour and the relative performance of the seven schemes are investigated in ample numerical experiments for both European put-on-the-min and put-on-the-average options.

Tables

Table 1. Table 1: Parameter sets for the two-asset Merton jump-diffusion model and European option.

	${\sigma }_1$	${\sigma }_2$	$\rho$	$\lambda$	${\gamma }_1$	${\gamma }_2$	$\hat{\rho }$	${\delta }_1$	${\delta }_2$	$r$	$K$	$T$
Set 1	0.12	0.15	0.30	0.60	-0.10	0.10	-0.20	0.17	0.13	0.05	100	1
Set 2	0.30	0.30	0.50	2	-0.50	0.30	-0.60	0.40	0.10	0.05	40	0.5
Set 3	0.20	0.30	0.70	8	-0.05	-0.20	0.50	0.45	0.06	0.05	40	1

Equations

$$S_{\tau }^{(1)}$$

$$S_{\tau }^{(2)}$$

$${\kappa }_i={\mathbb{E}}_{\mathbb{Q}}[e^{Y^{(i)}}-1](i=1,2).$$

$${\Sigma }_Y=\left(\begin{array}{cc}{\delta }_1^2& \hat{\rho }{\delta }_1{\delta }_2\\ \hat{\rho }{\delta }_1{\delta }_2& {\delta }_2^2\end{array}\right),$$

$$M_{\mathbb{Q},Y}(u)={\mathbb{E}}_{\mathbb{Q}}\left[e^{u^{\top }Y}\right]=\exp \left({\gamma }_1{u}_1+{\gamma }_2{u}_2+\frac{1}{2}\left({\delta }_1^2{u}_1^2+2\hat{\rho }{\delta }_1{\delta }_2{u}_1{u}_2+{\delta }_2^2{u}_2^2\right)\right),u=\left(\begin{array}{c}{u}_1\\ u_2\end{array}\right),$$

$${\kappa }_i=e^{{\gamma }_i+\frac{1}{2}{\delta }_i^2}-1(i=1,2).$$

$$\frac{\partial v}{\partial t}=$$

$$-(r+\lambda )v+\lambda {\int }_0^{\infty }{\int }_0^{\infty }v(s_1{y}_1,s_2{y}_2,t)f(y_1,y_2)\mathrm{d}{y}_1\mathrm{d}{y}_2$$

$$f(y_1,y_2)=\frac{1}{2\pi {\delta }_1{\delta }_2\sqrt{1-{\hat{\rho }}^2}{y}_1{y}_2}\exp \left(-\frac{{\left(\frac{\ln (y_1)-{\gamma }_1}{{\delta }_1}\right)}^2+{\left(\frac{\ln (y_2)-{\gamma }_2}{{\delta }_2}\right)}^2-2\hat{\rho }\left(\frac{\ln (y_1)-{\gamma }_1}{{\delta }_1}\right)\left(\frac{\ln (y_2)-{\gamma }_2}{{\delta }_2}\right)}{2(1-{\hat{\rho }}^2)}\right)$$

$$v(s_1,s_2,0)=\varphi (s_1,s_2),$$

$$\frac{{\partial }^{2}v}{\partial s_1^2}=0\text{leavevmode}\text{leavevmode}(\text{if}\text{leavevmode}{s}_1=S_{\max })\text{mbox}and\frac{{\partial }^{2}v}{\partial s_2^2}=0\text{leavevmode}\text{leavevmode}(\text{if}\text{leavevmode}{s}_2=S_{\max }).$$

$$\mathcal{D}v=\frac{1}{2}{\sigma }_1^2{s}_1^2\frac{{\partial }^{2}v}{\partial s_1^2}+\rho {\sigma }_1{\sigma }_2{s}_1{s}_2\frac{{\partial }^{2}v}{\partial s_1\partial s_2}+\frac{1}{2}{\sigma }_2^2{s}_2^2\frac{{\partial }^{2}v}{\partial s_2^2}+(r-\lambda {\kappa }_1)s_1\frac{\partial v}{\partial s_1}+(r-\lambda {\kappa }_2)s_2\frac{\partial v}{\partial s_2}-(r+\lambda )v.$$

$$(s_{1,i},s_{2,j})\in [0,S_{\max }]\times [0,S_{\max }](0\le i\le m_1,0\le j\le m_2),$$

$$u^{\prime }(s_i)\approx {\omega }_{i,-1}u(s_{i-1})+{\omega }_{i,0}u(s_i)+{\omega }_{i,1}u(s_{i+1}),$$

$${\omega }_{i,-1}=\frac{-h_{i+1}}{h_i(h_i+h_{i+1})},{\omega }_{i,0}=\frac{h_{i+1}-h_i}{h_i{h}_{i+1}},{\omega }_{i,1}=\frac{h_i}{h_{i+1}(h_i+h_{i+1})}$$

$$u^{\prime \prime }(s_i)\approx {\omega }_{i,-1}u(s_{i-1})+{\omega }_{i,0}u(s_i)+{\omega }_{i,1}u(s_{i+1}),$$

$${\omega }_{i,-1}=\frac{2}{h_i(h_i+h_{i+1})},{\omega }_{i,0}=\frac{-2}{h_i{h}_{i+1}},{\omega }_{i,1}=\frac{2}{h_{i+1}(h_i+h_{i+1})}.$$

$$A^{(D)}V(t)$$

$$A^{(D)}=A^{(M)}+A_1+A_2,$$

$$A^{(M)}$$

$$A_1$$

$$A_2$$

$$X_k=\mathrm{diag}(s_{k,0},s_{k,1},\dots ,s_{k,m_k})$$

$$\mathcal{J}(s_1,s_2):=$$

$$=$$

$$\overline{\mathcal{J}}(x_1,x_2):=\lambda {\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\overline{v}({\xi }_1,{\xi }_2,t)\overline{f}({\xi }_1-x_1,{\xi }_2-x_2)\mathrm{d}{\xi }_1\mathrm{d}{\xi }_2=\lambda (\overline{v}\star \star \overline{f})(x_1,x_2),$$

$${\overline{J}}_{k,l}=\lambda \sum _{j=-{\check{M}}_2+1}^{{\check{M}}_2}\sum _{i=-{\check{M}}_1+1}^{{\check{M}}_1}{\overline{V}}_{i,j}{\overline{F}}_{i-k,j-l},$$

$${\overline{f}}_{i-k,j-l}$$

$${\overline{A}}^{(J)}=\lambda \overline{F}$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}\overline{F}=\sum _{d=-M_2+1}^{M_2-1}E(M_2,d)\otimes {\overline{F}}_d^1$$

$${\overline{F}}_d^1=\sum _{c=-M_1+1}^{M_1-1}E(M_1,c){\overline{F}}_{c,d}$$

$$\overline{J}=$$

Full text

1 Introduction

The inherent risk diversification of multi-name or rainbow options, for which the payoff depends on multiple underlying asset prices, has made them increasingly popular over the recent years. Margrabe [31] proposed such a rainbow option, which allows to exchange one asset for another, and derived a valuation formula under the Black–Scholes framework. The term rainbow option was later introduced by Rubinstein [35], who referred to the number of underlying assets as the number of colors of the option. The development of fast, accurate and stable methods for obtaining the fair values of options forms a central topic in the field of computational finance. The fair value of an option is equal to its expected discounted payoff under a risk-neutral measure. Acquiring a semi-closed analytical valuation formula, as for example Black & Scholes [5] and Margrabe [31] did, would be most favorable, but this is often not possible. Accordingly, there is a strong demand for reliable numerical methods that can efficiently approximate fair option values.

For the numerical approximation of the fair value of an option, three main approaches are considered in the literature up to now. Firstly, Monte Carlo methods estimate the expected discounted payoff value by computing sample means. The second way is by numerical integration, which is employed in for example the Carr–Madan method [9] and the COS method by Fang & Oosterlee [13]. The third main approach is to numerically solve a time-dependent partial differential equation (PDE) that holds for the option value. This approach was introduced in computational finance by Schwartz [38], who considered a finite difference discretization of the celebrated Black–Scholes PDE. The development, analysis and application of numerical methods for time-dependent PDEs has since then been widely studied in the literature for valuing many types of options under a broad range of asset pricing models, see for example [21, 39, 41].

It is well-known in financial markets that sudden, large changes in asset prices can occur. In view of this common phenomenon, Merton [32] added a jump term to the original stochastic process for the asset price assumed by Black & Scholes [5], yielding a so-called jump-diffusion model. Under such an advanced asset price process, the option value satisfies a one-dimensional time-dependent partial integro-differential equation (PIDE). The integral part, which stems from the jump term, runs over the whole asset price domain, and hence, is nonlocal.

In the present paper we consider the extension of the prominent Merton jump-diffusion model to two assets where the jumps in the asset prices occur simultaneously. This yields a two-dimensional time-dependent PIDE for the values of two-asset options. The efficient and stable numerical solution of this equation poses two principal challenges: (a) an effective treatment of the two-dimensional PDE part, which involves a mixed derivative term, and (b) an effective treatment of the two-dimensional integral part, which is nonlocal.

We follow the well-known and versatile method of lines (MOL) approach [19], whereby the PIDE is first discretized in the spatial variables and subsequently in the temporal variable. Spatial discretization of the integral part gives rise to a very large, dense matrix that is computationally expensive to work with. However, by a change of variables, the double integral can be turned into a two-dimensional cross-correlation. Spatial discretization of this new integral version leads to a matrix having a block-Toeplitz structure for which matrix-vector products can be computed efficiently using the fast Fourier transform (FFT). For one-dimensional PIDEs, the application of FFT in approximating the integral part has been advocated in [1, 2, 12]. In particular, Almendral & Oosterlee [1] observed that the matrix corresponding to the jump term is Toeplitz and can be embedded in a circulant matrix, so that matrix-vector products can be computed exactly by means of FFT. In [10, 37] the FFT technique has next been applied in the case of two-dimensional PIDEs with a two-dimensional integral part. Salmi, Toivanen & von Sydow [37] considered a useful embedding of a block-Toeplitz matrix into a block-circulant one. Here we shall employ a similar idea, which was introduced in a different context in Barrowes, Teixeira & Kong [3].

After spatial discretization of the two-dimensional Merton PIDE one arrives at a very large, stiff system of ordinary differential equations (ODEs). The application of common implicit temporal discretization methods, such as the Crank–Nicolson scheme, requires the solution of very large, dense linear systems in every time step. The main aim of the present paper is to study a variety of recent and novel operator splitting schemes that effectively overcome this.

For one-dimensional option valuation PIDEs, Cont & Voltchkova [11] proposed an implicit-explicit (IMEX) Euler time-stepping scheme, where the PDE part is treated implicitly and the integral part explicitly. This IMEX Euler scheme is only first-order consistent. Various authors subsequently considered higher-order IMEX schemes for PIDEs in finance, see for example [8, 14, 30]. Salmi & Toivanen [36] recommended the IMEX CNAB scheme, where the PDE part is handled in a Crank–Nicolson manner and the integral part is treated in a two-step Adams–Bashforth fashion. This second-order two-step IMEX scheme has next been applied to two-dimensional option valuation PIDEs in for example [23, 37].

An alternative approach was proposed by Tavella & Randall [39] and d’Halluin, Forsyth & Vetzal [12]. These authors considered the Crank–Nicolson scheme together with a fixed-point iteration on the integral part for the solution of the linear system in each time step. This approach has been applied to two-dimensional option valuation PIDEs in Clift & Forsyth [10], including the two-dimensional Merton PIDE. It is found that in general two to three iterations are sufficient to reach a satisfactory error level. When the number of fixed-point iterations is frozen, a particular one-step IMEX scheme is obtained, see Section 4 below.

For multidimensional option valuation PDEs (without integral part), operator splitting schemes are widely considered in the literature. Notably, alternating direction implicit (ADI) methods, which combine implicit unidirectional correction stages with explicit stages involving the mixed derivative terms, are highly effective, see for example [16, 21, 23, 25, 40]. Hundsdorfer & in ’t Hout [18] recently studied a novel class of multistep stabilizing correction methods for PDEs with mixed derivative terms and showed that such methods can be competitive with the well-established one-step modified Craig–Sneyd (MCS) scheme [25] in the application to the familiar two-dimensional Heston PDE [17].

In the past years, ADI methods have been extended and applied to two-dimensional PIDEs in finance. In [28] Strang splitting was employed between the PDE and integral parts and an ADI scheme was used for the PDE subproblems, see also [27]. The direct extension of ADI methods to PIDEs has recently been considered in [24, 29]. Here the implicit unidirectional stages are combined with explicit stages involving both the mixed derivative term and the integral part. In ’t Hout & Toivanen [24] defined three adaptations of the MCS scheme to two-dimensional PIDEs that are all second-order consistent. They found that the adaptation where the integral part is treated in a two-step Adams–Bashforth fashion is preferred in the application to the well-known two-dimensional Bates PIDE [4], which merges the Heston PDE and the one-dimensional Merton PIDE.

An outline of the rest of this paper is as follows. Section 2 describes the asset price dynamics under the two-asset Merton jump-diffusion model and introduces the two-dimensional PIDE that holds for the value of a rainbow option under this model. Section 3 deals with the first step in the MOL approach, the spatial discretization of the two-dimensional Merton PIDE. Here particular attention is given to a suitable discretization of the double integral part, which enables application of the FFT algorithm from [3]. Section 4 concerns the second step in the MOL approach, the temporal discretization of the semidiscretized PIDE. We formulate seven recent and novel operator splitting schemes of both the IMEX and the ADI kind, which all treat the integral part in an explicit fashion. In Section 5 extensive numerical experiments are presented. Here two types of rainbow options are considered: a European put-on-the-min option and a European put-on-the-average option. For all seven splitting schemes, the behaviour of the temporal discretization error is investigated in detail. Moreover, employing the recent semi-closed analytical formula derived in Boen [6], also the total discretization error is studied in the case of the put-on-the-min option. The final Section 6 gives conclusions.

2 The two-asset Merton jump-diffusion model

Let $(\Omega,\mathcal{F},\{\mathcal{F}_{\tau}\},\mathbb{Q})$ be a probability measure space, with $\mathbb{Q}$ a risk-neutral measure. Under a two-asset jump-diffusion model, the asset prices under $\mathbb{Q}$ are given by

[TABLE]

where

• •

$r$ is the risk-free interest rate,

• •

$\sigma_{i}$ ($i=1,2$) is the instantaneous volatility of asset $i$, conditional on the event that no jumps occur,

• •

$W_{\tau}=(W^{(1)}_{\tau},W^{(2)}_{\tau})^{\top}$ consists of two correlated standard Brownian motions, with correlation coefficient $\rho$,

• •

$\lambda$ is the jump intensity of the Poisson arrival process $N=\{N_{\tau},\tau\geq 0\}$,

• •

$Y^{(i)}_{k}(k=1,2,3,\ldots)$ are independent jump sizes that have an identical distribution to the random jump size $Y^{(i)}$ and, for any given $k$, the jumps $Y^{(i)}_{k}\ (i=1,2)$ occur simultaneously, driven by the same Poisson arrival process $N$, and are correlated with correlation coefficient $\widehat{\rho}$,

• •

$\kappa_{i}$ is the expected relative jump size

[TABLE]

The $W_{\tau}$, $N_{\tau}$ and $Y=(Y^{(1)},Y^{(2)})^{\top}$ are assumed to be independent of each other whenever $\tau\geq 0$.

In the two-dimensional version of the Merton model [32], the random vector $Y$ is bivariate normally distributed with mean $(\gamma_{1},\gamma_{2})^{\top}$ and covariance matrix

[TABLE]

where $\delta_{i}^{2}$ is the variance of $Y^{(i)}$ $(i=1,2)$. The moment generating function of $Y$ is then given by

[TABLE]

and the expected relative jump size is

[TABLE]

Under the two-asset Merton jump-diffusion model, the value $v=v(s_{1},s_{2},t)$ of a European-style222A European-style option can only be exercised at the maturity date $T$. option with maturity date $T>0$ and $s_{i}$ ($i=1,2$) representing the price of asset $i$ at time $\tau=T-t$, satisfies the following PIDE [10]:

[TABLE]

whenever $s_{1}>0$, $s_{2}>0$, $0<t\leq T$. Here

[TABLE]

is the probability density function of a bivariate lognormal distribution. The initial condition for (2) is given by

[TABLE]

where $\phi$ denotes the payoff function of the option, that is, the option value at expiry. Further, as the boundary condition, it holds that the PIDE (2) is itself fulfilled on the two sides $s_{1}=0$ and $s_{2}=0$, respectively.

3 Spatial discretization

For the numerical solution of the initial-boundary value problem for (2) we consider the well-known and versatile method of lines (MOL) approach, which consists of two consecutive steps [19]: first the PIDE (2) is discretized in space and subsequently in time. This section discusses the spatial discretization. In the next section we turn to the temporal discretization.

3.1 Convection-diffusion-reaction part

To render the numerical solution of (2) feasible, the spatial domain is truncated to a bounded set $[0,S_{\rm max}]\times[0,S_{\rm max}]$ with fixed value $S_{\rm max}$ taken sufficiently large. On the two far sides $s_{1}=S_{\rm max}$ and $s_{2}=S_{\rm max}$ a linear boundary condition is imposed, which is common in finance,

[TABLE]

In this subsection we describe the finite difference discretization of the convection-diffusion-reaction part of (2), given by

[TABLE]

The option value function $v$ will be approximated at a nonuniform, Cartesian set of spatial grid points,

[TABLE]

with $s_{1,0}=s_{2,0}=0$ and $s_{1,m_{1}}=s_{2,m_{2}}=S_{\rm max}$. The nonuniform grid in each spatial direction is chosen such that relatively many grid points lie in a region of financial and/or numerical interest. It is generated by applying a smooth transformation to an artificial uniform grid. In Section 5 more details shall be given of the specific type of spatial grids that is used.

Let $0=s_{0}<s_{1}<\cdots<s_{m}=S_{\rm max}$ be any given smooth, nonuniform, unidirectional spatial grid with mesh widths $h_{i}=s_{i}-s_{i-1}$ ($1\leq i\leq m$) and let $u:[0,S_{\rm max}]\to\mathbb{R}$ be any given smooth function. The following second-order central finite difference formulas are employed for the approximation of the first and second derivatives of $u$:

[TABLE]

with

[TABLE]

and

[TABLE]

with

[TABLE]

for $1\leq i\leq m-1$. If $i=0$, no finite difference formulas are required in view of the degeneracy of (2) at the zero boundaries. If $i=m$, the first derivative is approximated by the first-order backward finite difference formula and the second derivative vanishes by the linear boundary condition (3.1). Concerning the mixed derivative term in (2), this is discretized by successively applying the relevant finite difference formulas for the first derivative in the two spatial directions.

Let the vector $V(t)=(V_{0,0}(t),V_{1,0}(t),\ldots,V_{m_{1}-1,m_{2}}(t),V_{m_{1},m_{2}}(t))^{\top}$ where entry $V_{i,j}(t)$ denotes the semidiscrete approximation to $v(s_{1,i},s_{2,j},t)$ whenever $0\leq i\leq m_{1}$, $0\leq j\leq m_{2}$. The semidiscrete version of the convection-diffusion-reaction part $\mathcal{D}v$ can then be written as

[TABLE]

with matrix

[TABLE]

where

[TABLE]

Here $I_{k}$, $X_{k}$, $D_{k}^{(l)}$ are given $(m_{k}+1)\times(m_{k}+1)$ matrices for $k,l\in\{1,2\}$, where $I_{k}$ is the identity matrix, $X_{k}$ is the diagonal matrix

[TABLE]

and $D_{k}^{(l)}$ is the matrix representing numerical differentiation of order $l$ in the $k$-th spatial direction by the relevant finite difference formula above. The matrix $A^{(M)}$ corresponds to the mixed derivative term in (2) and $A_{k}$ corresponds to all derivative terms in the $k$-th spatial direction ($k=1,2$). Further, the reaction term has been distributed equally across $A_{1}$ and $A_{2}$.

3.2 Integral part

In order to discretize the integral part in the PIDE (2), we apply a well-known transformation to the log-price variable. Let $x_{i}=\ln(s_{i})$ be the log-price of the $i$-th asset and $\eta_{i}=\ln(y_{i})$ ($i=1,2$), then

[TABLE]

where $\overline{v}(\xi_{1},\xi_{2},t)=v(e^{\xi_{1}},e^{\xi_{2}},t)$ and $\overline{f}(\eta_{1},\eta_{2})=$ $f(e^{\eta_{1}},e^{\eta_{2}})e^{\eta_{1}}e^{\eta_{2}}$ is the probability density function of a bivariate normal distribution. Since $\overline{f}$ is a real-valued function, it is clear that the double integral in the log-price variable can be viewed as a two-dimensional cross-correlation (see e.g. [7]):

[TABLE]

where $\star\star$ denotes the two-dimensional cross-correlation operator and $\overline{\mathcal{J}}(x_{1},x_{2})=\mathcal{J}(e^{x_{1}},e^{x_{2}})$. For the semidiscretization, we restrict the log-price domain to $[-X_{\rm max},X_{\rm max}]\times[-X_{\rm max},X_{\rm max}]$ with $X_{\rm max}=\ln(S_{\rm max})$ and consider a $M_{1}\times M_{2}$ uniform grid $(x_{1,k},x_{2,l})=(k\Delta x_{1},l\Delta x_{2})$ for ${k=-\widecheck{M}_{1}+1,\ldots,\widecheck{M}_{1},\ l=-\widecheck{M}_{2}+1,\ldots,\widecheck{M}_{2}}$, where $\widecheck{M}_{i}=M_{i}/2$, $\widecheck{M}_{i}\Delta x_{i}=X_{\rm max}$ and $M_{i}$ is chosen to be a power of 2 such that the mesh width $\Delta x_{i}$ is smaller than the smallest mesh width in the nonuniform $\ln(s_{i})$-grid ($i=1,2$). The double integral $\overline{\mathcal{J}}(x_{1,k},x_{2,l})$ is approximated by

[TABLE]

with approximation $\overline{V}_{i,j}\approx\overline{v}(x_{1,i},x_{2,j},t)$ defined below and $\overline{F}_{i-k,j-l}=\overline{f}_{i-k,j-l}\Delta x_{1}\Delta x_{2}$ where

[TABLE]

Let $E(M,d)$ denote the $M\times M$ matrix with ones on the $d$-th off-diagonal and zeros elsewhere. The matrix

[TABLE]

with

[TABLE]

and

[TABLE]

is then an asymmetrical two-level block-Toeplitz matrix of size $(M_{1}M_{2})\times(M_{1}M_{2})$. Its structure is illustrated in Appendix A. Let the vectors

[TABLE]

then

[TABLE]

This matrix-vector product can be computed very efficiently using the FFT algorithm for matrix-vector multiplications with asymmetric multilevel block-Toeplitz matrices, introduced in Barrowes, Teixeira & Kong [3]. This fast algorithm computes the matrix-vector product using two FFTs and one inverse FFT (IFFT), reducing the computational cost to $\mathcal{O}(M_{1}M_{2}\log(M_{1}M_{2}))$. The algorithm from [3] essentially embeds the block-Toeplitz matrix in a circulant matrix before applying the FFT, conveniently avoiding any wrap-around effect, which would have been present when applying the FFT directly to the block-Toeplitz matrix. Further, we mention that the algorithm from [3] is more efficient than that employed in [37], as the former uses only one-dimensional FFTs, whereas the latter requires two-dimensional FFTs.

Note that the spatial grid on which we discretize the PIDE (2) is nonuniform, and in general does not lead to a uniform log-price grid. We therefore (bi)linearly interpolate the option value approximations between the two grids directly before and after application of the algorithm of [3]. Linearly interpolating $V$ to $\overline{V}$ (using Lagrange basis functions on the $s$-grid, as in e.g. [12]), yields

[TABLE]

whenever $s_{1,\Pi(i)}\leq e^{x_{1,i}}\leq s_{1,\Pi(i)+1}$ and $s_{2,\Omega(j)}\leq e^{x_{2,j}}\leq s_{1,\Omega(j)+1}$, where $\phi_{\Pi(i)},\psi_{\Omega(j)}\in[0,1]$ are interpolation weights . This can be written as

[TABLE]

where the matrix $\mathcal{X}$ contains the interpolation weights such that (3.2) holds. Similarly, linearly interpolating $\overline{J}$ back to the original $s$-grid, leads to approximations $J_{k,l}$ of $\mathcal{J}(s_{1,k},s_{2,l})$ given by

[TABLE]

whenever $e^{x_{1,\Lambda(k)}}\leq s_{1,k}\leq e^{x_{1,\Lambda(k)+1}}$ and $e^{x_{2,\Upsilon(l)}}\leq s_{2,l}\leq e^{x_{2,\Upsilon(l)+1}}$, where $\overline{\phi}_{\Lambda(k)},\overline{\psi}_{\Upsilon(l)}\in[0,1]$ are interpolation weights. This can be written as

[TABLE]

where the matrix $\overline{\mathcal{X}}$ contains the interpolation weights such that (3.2) holds. Consequently, we have

[TABLE]

The entire procedure is summarized in Algorithm 1.

Note that in the interpolation step from the uniform log-price grid to the nonuniform price grid, we obtain values $J_{i,j}$ for $i=1,\ldots,m_{1},\ j=1,\ldots,m_{2}$. On the boundaries $s_{1}=0$ and $s_{2}=0$, the two-dimensional Merton PIDE (2) reduces to the one-dimensional Merton PIDE and the jump part can be computed using a similar FFT technique (see f.i. [12]).

Denoting the approximation of $\mathcal{J}$ on the full grid (3.2) by $A^{(J)}V(t)$, we arrive at the following semidiscrete system of ODEs:

[TABLE]

for $0<t\leq T$, where $A=A^{(D)}+A^{(J)}=A^{(M)}+A_{1}+A_{2}+A^{(J)}$.

The maximum norm of the matrix $A^{(J)}$ can be bounded independently of the spatial grid by a moderate constant whenever $\lambda$ is moderate. Accordingly, the dense matrix $A^{(J)}$ constitutes a nonstiff part of the semidiscrete system, whereas the sparse matrix $A^{(D)}$ constitutes a stiff part.

4 Temporal discretization

For the temporal discretization of the semidiscrete system (3.5) we shall investigate seven operator splitting schemes. Each of these schemes conveniently treats the nonstiff integral part, given by the matrix $A^{(J)}$, in an explicit fashion.

Let integer $N\geq 1$ be given and step size $\Delta t=T/N$. Each of the following schemes defines an approximation $V^{n}$ to $V(t_{n})$ at the temporal grid point $t_{n}=n\Delta t$ successively for $n=1,2,\ldots,N$. The initial vector $V(0)=V^{0}$ is determined by the payoff function and will be specified in Section 5.

Crank–Nicolson Forward Euler (CNFE) scheme: In this basic method the convection-diffusion-reaction part, including the mixed derivative, is handled implicitly using the Crank–Nicolson scheme and the integral part is treated explicitly in a forward Euler manner:

[TABLE]

Method (4.1) is a one-step IMEX scheme. Due to the application of forward Euler, its order333This refers to the classical order of consistency, that is, for fixed nonstiff ODEs. is just equal to one. 2. 2.

Crank–Nicolson scheme with fixed-point iteration (CNFI): In Tavella & Randall [39] and d’Halluin, Forsyth & Vetzal [12] the Crank–Nicolson scheme has been proposed with a fixed-point iteration on the integral part:

[TABLE]

for $k=1,2,\ldots,l$ and $V^{n}=Y_{l}$ using the starting value $Y_{0}=V^{n-1}$. Clearly, with just one iteration, method (4.1) is obtained. Method (4.2) can be applied with a dynamic convergence criterion. Numerical experiments in [10, 12] reveal that this generally leads to $l=2$ or $l=3$ iterations. For comparable computational work to the other schemes, formulated below, we choose here a fixed number of $l=2$ iterations. Then method (4.2) can be regarded as a one-step IMEX scheme, where the integral part is treated by the explicit trapezoidal rule, also called the modified Euler method. It can be verified, by Taylor expansion, that its order is equal to two. 3. 3.

Implicit-Explicit Trapezoidal Rule (IETR): A different combination of the implicit trapezoidal rule (Crank–Nicolson) for the convection-diffusion-reaction part and its explicit variant for the integral part is considered by in ’t Hout [21]:

[TABLE]

Here the first two internal stages $Y_{0}$, $\widehat{Y}_{0}$ are explicit, whereas the third internal stage $Y_{1}$ is implicit. Method (4.3) is also a one-step IMEX scheme with order equal to two. This method is already well-known in the literature on the numerical solution of PDEs, without integral part, see e.g. Hundsdorfer & Verwer [19]. 4. 4.

Crank–Nicolson Adams–Bashforth (CNAB) scheme: In the CNAB method, the convection-diffusion-reaction part is again treated by the Crank–Nicolson scheme, but the integral part is now handled in a two-step Adams–Bashforth fashion:

[TABLE]

Method (4.4) forms a two-step IMEX scheme and is of order equal to two. This method has first been studied for the application to PIDEs in Salmi & Toivanen [36] and Salmi, Toivanen & von Sydow [37]. Previously, it has been investigated for the numerical solution of PDEs in e.g. [15, 19, 20].

It is well-known that the nonsmoothness of the initial (payoff) function can have an adverse effect on the convergence behaviour of the Crank–Nicolson scheme, which can be resolved by first applying, at initial time $t=0$, two time steps with step size $\Delta t/2$ using the implicit Euler method to define the approximation $V^{1}$ to $V(t_{1})$, see e.g. [21, 33, 34]. In the present application to PIDEs, the implicit Euler method can be computationally demanding. We therefore apply, for all four schemes above, the IMEX Euler scheme as the starting method:

[TABLE]

In the following, three recent operator splitting schemes for PIDEs are formulated that employ a subsequent, useful splitting of the two-dimensional convection-diffusion-reaction part, given by the matrix $A^{(D)}$.

One-step adaptation of the modified Craig–Sneyd (MCS) scheme: This method is a direct adaptation of a well-established one-step ADI scheme:

[TABLE]

where $\theta>0$ is a given parameter. The MCS scheme was introduced by in ’t Hout & Welfert [25] for the numerical solution of PDEs with mixed derivative terms. The above, direct adaptation to PIDEs has recently been studied in in ’t Hout & Toivanen [24]. Method (4.5) is of order two for any value $\theta$. Here we make the common choice $\theta=\frac{1}{3}$, which is motivated by stability and accuracy results in the literature for two-dimensional problems, see e.g. [22, 24, 25, 26]. It is readily seen that in (4.5) the integral part and mixed derivative term are both treated by the explicit trapezoidal rule. Notice that the explicit stages $\widehat{Y}_{0}$, $\widetilde{Y}_{0}$ can be combined, so that the integral part is evaluated twice per time step. The implicit stages $Y_{j}$, $\widetilde{Y}_{j}$ (for $j=1,2$) are often called stabilizing corrections and are unidirectional. The pertinent linear systems for these stages are tridiagonal, so that they can be solved very efficiently using an a priori $LU$ factorization. 2. 6.

Two-step adaptation of the MCS (MCS2) scheme: In [24] an alternative adaptation of the MCS scheme to PIDEs has been proposed, where the integral part is dealt with in a two-step Adams–Bashforth fashion:

[TABLE]

Method (4.6) is also of order two for any value $\theta$. We take again $\theta=\frac{1}{3}$ and, for starting (4.6), define $V^{1}$ by the one-step method (4.5). 3. 7.

Stabilizing correction two-step Adams-type (SC2A) scheme: In Hundsdorfer & in ’t Hout [18] a novel class of stabilizing correction multistep methods has recently been investigated for the numerical solution of PDEs. We select here a prominent member of this class, the two-step Adams-type scheme called SC2A, and consider its direct adaptation to PIDEs:

[TABLE]

with coefficients $(\widehat{b}_{1},\widehat{b}_{2})=\left(\frac{3}{2},-\frac{1}{2}\right)$ and $(\widecheck{b}_{1},\widecheck{b}_{2})=\left(\frac{3}{2}-\theta,-\frac{1}{2}+\theta\right)$. The integral part and mixed derivative term are now both handled by the two-step Adams–Bashforth scheme. Method (4.7) is also of order two for any value $\theta$. Following [18] we take $\theta=\frac{3}{4}$, which is motivated by stability and accuracy results. For starting (4.7), the one-step method (4.5) is used with $\theta=\frac{1}{3}$ to define $V^{1}$.

In each of the seven time-stepping schemes described above, any matrix-vector product with the matrix $A^{(J)}$ means applying Algorithm 1, hence requiring two FFTs and one IFFT. In the schemes (4.1), (4.4), (4.6), (4.7) only one such matrix-vector product arises per time step, whereas two such products appear in the schemes (4.2), (4.3), (4.5).

5 Numerical study

In this section we examine through extensive numerical experiments the convergence behaviour of the seven operator splitting schemes formulated in Section 4 in the numerical solution of the semidiscrete two-dimensional Merton PIDE. For the numerical study we consider the three parameter sets for the two-asset Merton jump-diffusion model and European-style option given in Table 1. Set 1 is taken from Clift & Forsyth [10]. Set 2 has the same diffusion parameters as in Zvan, Forsyth & Vetzal [42], and is complemented with jump parameters where the intensity $\lambda$ is taken larger than in Set 1. Set 3 is fully new and yields more intensive jumps. Here $\lambda T$ is rather large, which poses a particular numerical challenge. Notice that for all three sets the correlation coefficients $\rho$ and $\widehat{\rho}$ are both nonzero.

Two types of rainbow options are considered: a European put-on-the-min option and a European put-on-the-average option. Their payoff functions are, respectively,

[TABLE]

and

[TABLE]

Clearly, the two payoff functions $\phi$ are continuous, but not everywhere continuously differentiable. For each option a specific nonuniform spatial grid $\{(s_{1,i},s_{2,j}):0\leq i\leq m_{1},0\leq j\leq m_{2}\}$ is chosen, where the nonsmoothness of $\phi$ is taken into account. The number of grid points is taken to be the same in both directions, $m_{1}=m_{2}=m$. For the put-on-the-min option, relatively many points $s_{1,i}$ and $s_{2,j}$ are placed around the locations $s_{1}=K$ and $s_{2}=K$, respectively. For the put-on-the-average option, the payoff is nonsmooth on the line segment given by $s_{1}+s_{2}=2K$. Here the grid in each direction is taken to be uniform on the interval $[0,2K]$ and nonuniform outside, with relatively many points inside this interval. For both options, the nonuniform grids are generated by applying a smooth transformation to an artificial uniform grid, see e.g. [21]. For Sets 1, 2, 3 we heuristically choose $S_{\rm max}=5K,30K,50K$, respectively, in the case of the put-on-the-min and $S_{\rm max}=5K,15K,25K$, respectively, in the case of the put-on-the-average option. As an illustration, Figure 1 displays the two spatial grids if $m=50$, $K=100$, $S_{\rm max}=5K$. In Figure 2 the graphs of the two payoff functions are shown.

The initial vector $V(0)=V^{0}$ is defined by pointwise evaluation of the payoff function $\phi$ at the spatial grid points, except at those points that lie closest to the set of nonsmoothness of $\phi$, where cell averaging is applied. For $k=1,2$ let

[TABLE]

with $s_{k,-1/2}=0$ and $s_{k,m+1/2}=S_{\rm max}$. Then we define [21]

[TABLE]

whenever the cell

[TABLE]

has a nonempty intersection with the set of points where $\phi$ is nonsmooth. For both payoffs under consideration, the above integral is readily calculated.

In the following we investigate, for all seven operator splitting methods formulated in Section 4 for the semidiscrete system (3.5), the temporal discretization error at $t=T$ on a natural region of financial interest (ROI),

[TABLE]

Clearly, the temporal discretization error is measured in the maximum norm, which is the most relevant norm in finance. Numerical experiments reveal that the evaluation of the two-dimensional integral part is computationally the dominant part in each time step, compare also [37]. In view of this, for a fair comparison of the different time-stepping methods, we consider the three schemes (4.2), (4.3), (4.5) with $N^{\prime}=N$ time steps and the four schemes (4.1), (4.4), (4.6), (4.7) with $N^{\prime}=2N$ time steps, as the latter four schemes require only one evaluation of the integral part per time step, whereas the former three schemes employ two such evaluations. A reference value for $V(T)$ has been computed by applying the MCS2 scheme with $3000$ time steps.

Figure 3 displays the temporal errors $\widehat{E}^{ROI}(m,N)$ for $m=150$ and a range of values $10\leq N\leq 1000$ for the put-on-the-min option (left column) and put-on-the-average option (right column) and the three parameter sets given in Table 1. For all seven splitting methods, the positive result holds that the temporal errors always lie below a moderate value and show a regular, monotonically decreasing behaviour as $N$ increases. The CNFE scheme (4.1) shows an order of convergence just equal to one, as expected. Each of the other six splitting schemes (4.2)–(4.7) attains a favourable order of convergence equal to two. Comparing the different second-order methods, the MCS2 scheme yields the smallest temporal error constant in all experiments. It is further interesting to notice that for Set 3, with more intensive jumps, each of the two-step schemes CNAB, MCS2, SC2A outperforms all of the one-step schemes under consideration. The range of the temporal error constants for the different second-order methods is also the largest for this set. We performed additional numerical experiments, choosing a coarser ($m=75$) as well as a finer ($m=250$) spatial grid. The obtained results were almost identical to those displayed in Figure 3, indicating that the observed temporal convergence orders are valid in a favourable, stiff sense. Recall that the CNFI scheme has been applied with a fixed number of $l=2$ iterations. We also carried out experiments for this scheme with $l=3$ iterations, taking now $N^{\prime}=\lceil 2N/3\rceil$ time steps for a fair comparison. This did not lead to a significant improvement in accuracy, however, and did not alter the above conclusions.

A semi-closed analytical formula has recently been derived in Boen [6] for the value of a European put-on-the-min option under the two-asset Merton model. This formula is stated in Appendix B. Hence, in this case, we can also study the total discretization error at $t=T$ on the region of interest,

[TABLE]

Figure 4 displays the total errors $E^{ROI}(m,N)$ for a range of values $10\leq m\leq 500$ and the three parameter sets given in Table 1, where we have chosen $N=\lceil m/3\rceil$ and $N^{\prime}$ as above. Clearly, with all splitting methods except CNFE, a second-order convergence behaviour is found, which is as desired. The total errors obtained with these six time-stepping methods are visually almost identical, indicating that the spatial discretization error dominates in these experiments.

6 Conclusion

In this paper we have studied seven different operator splitting schemes when applied to the two-dimensional Merton PIDE with a focus on IMEX and ADI methods. Each of the considered schemes conveniently treats the nonstiff, nonlocal integral part in an explicit fashion. Any matrix-vector product arising from the semidiscretization of the integral part is computed efficiently by means of FFT using Algorithm 1. Through ample numerical experiments, we have investigated the convergence behaviour of the seven schemes and examined their relative performance. Except for the first-order CNFE scheme, all of the considered splitting schemes showed a stiff order of temporal convergence equal to two, which is as desired. The two-step MCS2 scheme with $\theta=\frac{1}{3}$ stood out as yielding the smallest temporal error constant in all the conducted experiments. Together with the positive results recently derived in [24], we therefore recommend this splitting scheme for the efficient and stable temporal discretization of two-dimensional PIDEs in financial option valuation.

Acknowledgements

The authors acknowledge the support of the Research Fund (BOF) of the University of Antwerp (41/FA070300/3/FFB150337).

Figures

Appendix A The block-Toeplitz structure of $\overline{A}^{(J)}$

As an illustration, if $M_{1}=M_{2}=4$, then $\overline{F}$ in the matrix $\overline{A}^{(J)}=\lambda\overline{F}$ is of the following form:

This is a block-Toeplitz matrix with Toeplitz blocks $\overline{F}^{1}_{d}$, ${d=-3,\ldots,3}$:

[TABLE]

Appendix B Semi-closed analytical formula for put-on-the-min option

Let $M(x_{1},x_{2},\Sigma)$ denote the bivariate normal cumulative distribution function, evaluated at $(x_{1},x_{2})$, with mean $(0,0)^{\top}$ and covariance matrix $\Sigma$. Under the two-asset Merton jump-diffusion model, the value of a European put-on-the-min option at inception is given in Boen [6]:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

and

[TABLE]