Numerical analysis of American option pricing in a two-asset jump-diffusion model

American options play a crucial role in financial markets and are widely traded for both hedging and speculative purposes. Unlike European options, which can only be exercised at expiration, American options allow exercise at any time before expiration. While this flexibility enhances their appeal across equities, commodities, and bonds, it also introduces substantial mathematical and computational challenges due to the absence of closed-form solutions in most cases.

From a mathematical perspective, American option pricing is formulated as a variational inequality due to its optimal stopping feature. The problem becomes particularly challenging when incorporating jumps, as the governing equations involve nonlocal integro-differential terms [70, 60]. These complexities necessitate advanced numerical methods to achieve accurate valuations.

A common approach to tackling the variational inequality in American option pricing is to reformulate it as a partial (integro-)differential complementarity problem [16, 71, 17, 33, 13]. This reformulation captures the early exercise feature through time-dependent complementarity conditions, effectively handling the boundary between the exercise and continuation regions. Finite difference techniques are widely used to solve these problems, resulting in nonlinear discretized equations that require iterative solvers such as the projected successive over-relaxation method (PSOR) [23] and penalty methods [71]. For models incorporating jumps, fixed-point iterations can be used to address the integral terms, as shown in [17] for two-asset jump-diffusion models. Additionally, efficient operator splitting schemes—including implicit-explicit and alternating direction implicit types—have recently been proposed for American options under the two-asset Merton jump-diffusion model [13], while the method of lines has been applied to the same problem under the two-asset Kou model [37].

An alternative approach to pricing American options is to approximate them using Bermudan options, a discrete-time counterpart, and analyze the convergence of the solution as the early exercise intervals shrink to zero. In this context, Fourier-based COS methods [62] have been applied to 2-D Bermudan options. However, for this application, the COS method requires particularly careful boundary tracking despite its high-order convergence for piecewise smooth problems. Notably, as observed in [32, 28], the COS method may produce negative option prices for short maturities, akin to cases where early exercise opportunities become more frequent. This phenomenon is closely tied to its potential loss of monotonicity, a limitation highlighted in previous discussions. Consequently, this approach may struggle to ensure convergence from Bermudan to American options.

This paper bridges the gap in existing numerical methods by introducing an efficient, straightforward-to-implement monotone integration scheme for the variational inequalities governing American options under the two-asset Merton jump-diffusion model. Our approach simultaneously handles cross-derivative terms and nonlocal integro-differential terms, simplifying the design of monotone schemes and ensuring convergence to the viscosity solution. In doing so, we address key computational challenges in current numerical techniques.

The main contributions of our paper are outlined below.

• (i)

  We present the localized variational inequality for pricing American options under the two-asset Merton jump-diffusion model, formulated on an infinite domain in log-price variables with a finite interior and artificial boundary conditions. Using a probabilistic technique, we demonstrate that the difference between the solutions of the localized and full-domain variational inequalities decreases exponentially with respect to the log-price domain size. In addition, we establish that the localized variational inequality satisfies a comparison result.

• (ii)

  We develop a monotone scheme for the variational inequality that explicitly enforces the inequality constraint. We solve a 2-D Partial Integro-Differential Equation (PIDE) at each timestep to approximate the continuation value, followed by an intervention action applied at the end of the timestep. Using the closed-form Fourier transform of the Green’s function, we derive an infinite series representation where each term is non-negative. This enables the direct approximation of the PIDE’s solutions via 2-D convolution integrals, using a monotone numerical integration method.

• (iii)

  We implement the monotone integration scheme efficiently by exploiting the Toeplitz matrix structure and using Fast Fourier Transforms (FFTs) combined with circulant convolution. The implementation process includes expanding the inner summation’s convolution kernel into a circulant matrix, followed by transforming the double summation kernel into a circulant block structure. This allows the circulant matrix-vector product to be efficiently computed as a circulant convolution using 2-D FFTs.

• (iv)

  We prove that the proposed monotone scheme is both $\ell_{\infty}$-stable and consistent in the viscosity sense, ensuring pointwise convergence to the viscosity solution of the variational inequality as the discretization parameter approaches zero.

• (v)

  Extensive numerical results demonstrate strong agreement with benchmark solutions from published test cases, including those obtained via operator splitting methods, establishing our method as a reference for validating numerical techniques.

While this work focuses on American option pricing under a correlated two-asset Merton jump-diffusion model, the core methodology—particularly the infinite series representation of the Green’s function, where each term is non-negative—can be extended to other stochastic control problems. One such application is asset allocation with a stock index and a bond index. For discrete rebalancing, the 2-D extension is straightforward, as time advancement can be handled using a similar infinite series representation. For continuous rebalancing, a natural approach is to start from the discrete setting and leverage insights from [69] to analyze the limit as the rebalancing interval approaches zero.

The extension to the 2-D Kou model presents significant additional challenges due to its piecewise-exponential structure, leading to complex multi-region double integrals. We leave this extension for future work and refer the reader to Subsection 3.2, where we discuss these difficulties and outline a potential neural network-based approach.

The remainder of the paper is organized as follows. In Section 2, we provide an overview of the two-asset Merton jump-diffusion model and present the corresponding variational inequality. We then define a localized version of this problem, incorporating boundary conditions for the sub-domains. Section 3 introduces the associated Green’s function and its infinite series representation. In Section 4, we describe a simple, yet effective, monotone integration scheme based on a composite 2-D quadrature rule. Section 5 establishes the mathematical convergence of the proposed scheme to the viscosity solution of the localized variational inequality. Numerical results are discussed in Section 6, and finally, Section 7 concludes the paper and outlines directions for future research.

For the underlying process $(X_{t},Y_{t}),~{}t\in[0,T]$, let $(a,b)$ be the state of the system. We denoted by $v^{\prime\prime}(a,b,t)$ the time-$t$ no-arbitrage price of a two-asset American option contract with maturity $T$ and payoff $\hat{v}(a,b)$. It is established that $v^{\prime\prime}(\cdot)$ is given by the optimal stopping problem [41, 45, 54, 59, 42, 70]

Here, $\gamma$ represents a stopping time; $\mathbb{E}^{x,y}_{t}$ denotes the conditional expectation under the risk-neutral measure $\mathfrak{Q}$, conditioned on $(X_{t},Y_{t})=(a,b)$. We focus on the put option case, where the payoff function $\hat{v}^{\prime}(\cdot)$ is bounded and continuous.

The methods of variational inequalities, originally developed in [11], are widely used for pricing American options, as evidenced by [70, 42, 60, 59], among many other publications. The value function $v^{\prime\prime}(\cdot)$, defined in (2.3), is known to be non-smooth, which prompts the use of the notion of viscosity solutions. This approach provides a powerful means for characterizing the value functions in stochastic control problems [20, 22, 21, 31, 59].

It is well-established that the value function $v^{\prime\prime}(\cdot)$, defined in (2.3), is the unique viscosity solution of a variational inequality as noted in [58, 60, 59]. While the original references describe the variational inequality using the spatial variables $(a,b)$, our approach employs a logarithmic transformation for theoretical analysis and numerical method development. Specifically, with $\tau=T-t$, and given positive values for $a$ and $b$, we apply the transformation $x=\ln(a)\in(-\infty,\infty)$ and $y=\ln(b)\in(-\infty,\infty)$. With ${\bf{x}}=(x,y,\tau)$, we define $v^{\prime}({\bf{x}})\equiv v^{\prime}(x,y,\tau)=v^{\prime\prime}(e^{x},e^{y},T-t)$ and $\hat{v}(x,y)=\hat{v}^{\prime}(e^{x},e^{y})$. Consequently, $v^{\prime}(\cdot)$ is the unique viscosity solution of the variational inequality given by

Here, the differential and jump operators $\mathcal{L}(\cdot)$ and $\mathcal{J}(\cdot)$ are defined as follows

where $f(s_{x},s_{y})$ is the joint probability density function of $\left(\ln(\xi_{\scalebox{0.6}{$x$}}),\ln(\xi_{\scalebox{0.6}{$y$}})\right)$.

Under the log transformation, the formulation (2.4) is posed on an infinite spatial domain $\mathbb{R}^{2}$. For problem statement and convergence analysis of numerical schemes, we define a localized pricing problem with a finitem, open, spatial interior sub-domain, denoted by $\mathbb{D}_{\scalebox{0.7}{\text{in}}}\subset\mathbb{R}^{2}$. More specifically, with $x_{\scalebox{0.55}{$\min$}}<0<x_{\scalebox{0.55}{$\max$}}$ and $y_{\scalebox{0.55}{$\min$}}<0<y_{\scalebox{0.55}{$\max$}}$, where $x_{\scalebox{0.55}{$\min$}}$, $x_{\scalebox{0.55}{$\max$}}$, $|y_{\scalebox{0.55}{$\min$}}|$, and $y_{\scalebox{0.55}{$\max$}}$ are sufficiently large, $\mathbb{D}_{\scalebox{0.7}{\text{in}}}$ and its complement $\mathbb{D}_{\scalebox{0.7}{\text{out}}}^{\scalebox{0.6}{$\infty$}}$ are respectively defined as follows

Since the jump operator $\mathcal{J}(\cdot)$ is non-local, computing the integral (2.5) for $(x,y)\in\mathbb{D}_{\scalebox{0.7}{\text{in}}}$ typically requires knowledge of $v(\cdot)$ within the infinite outer boundary sub-domain $\mathbb{D}_{\scalebox{0.7}{\text{out}}}^{\scalebox{0.6}{$\infty$}}$. Therefor, appropriate boundary conditions must be established for $\mathbb{D}_{\scalebox{0.7}{\text{out}}}^{\scalebox{0.6}{$\infty$}}$. In the following, we define the definition domain and its sub-domains, discuss boundary conditions, and investigate the impact of artificial boundary conditions on $v(\cdot)$.

A proof of Lemma 2.1 is provided in Appendix A. This lemma establishes that, in log-price coordinates, the error between the localized and full-domain variational inequality solutions decays exponentially as the spatial interior domain size increases, i.e. the truncation error is of order $\mathcal{O}(e^{-A})$ as $A\to\infty$.¹¹1Equivalently, if $A^{\prime}$ represents the spatial domain size in the original price scale for both assets, the error is $\mathcal{O}(1/A^{\prime})$. The result is a local pointwise estimate, indicating that the localization error is more pronounced near the boundary. This rapid decay implies that smaller computational domains can be used, thereby significantly reducing computational costs.

The boundedness conditions $\|\hat{v}(\cdot)\|_{\infty}<\infty$ and $\|\hat{p}(\cdot)\|_{\infty}<\infty$ are satisfied for standard put options and Dirichlet boundary conditions. However, for payoffs that grow unbounded as $x,y\to+\infty$ (e.g. standard calls), a more refined analysis is required. Large-deviation techniques have been employed in pure-diffusion models (see [10]) to estimate the probability of large excursions beyond the truncated domain. These techniques can also be adapted to certain (finite-activity) jump processes (see, for example, [44]). A rigorous extension of these ideas to jump-diffusion models with unbounded payoffs, while theoretically possible, is considerably more involved and lies beyond the scope of this work. We therefore concentrate on bounded-payoff options—such as puts—where $\|\hat{v}(\cdot)\|_{\infty}$ remains finite.

For the remainder of the analysis, we choose the Dirichlet condition based on discounted payoff as follows

While more sophisticated boundary conditions might involve the asymptotic properties of the variational inequality (2.4a) as $x,y\to-\infty$ or $x,y\to\infty$, our observations indicate that these sophisticated boundary conditions do not significantly impact the accuracy of the numerical solution within $\Omega_{\scalebox{0.7}{\text{in}}}$. This will be illustrated through numerical experiments in Subsection 6.2.4.

We now write (2.1) in a compact form, which includes the terminal and boundary conditions in a single equation. We let $Dv({\mathbf{x}})$ and $D^{2}v({\mathbf{x}})$ represent the first-order and second-order partial derivatives of $v\left({\mathbf{x}}\right)$. The variational inequality (2.9) can be expressed compactly as

where

For a locally bounded function $\psi:\mathbb{D}\rightarrow\mathbb{R}$, where $\mathbb{D}$ is a closed subset of $\mathbb{R}^{n}$, we recall its upper semicontinuous (u.s.c. in short) and the lower semicontinuous (l.s.c. in short) envelopes given by

In the context of numerical solutions to degenerate parabolic equations in finance, the convergence to viscosity solutions is ensured when the scheme is stable, consistent, and monotone, provided that a comparison result holds [43, 8, 7, 9, 5, 6, 4]. Specifically, stability, consistency and monotonicity facilitate the identification of u.s.c. subsolutions and l.s.c. supersolutions through the respective use of $\limsup$ and $\liminf$ of the numerical solutions as the discretization parameter approaches zero.

Suppose $\underline{v}(\cdot)$ and $\overline{v}(\cdot)$ respectively denote such subsolution and supersolution within a region, referred to as $\mathcal{S}$, where $\mathcal{S}=\mathbb{S}\times[0,T]$ for an open set $\mathbb{S}\subseteq\mathbb{R}^{2}$. By construction using $\limsup$ for $\underline{v}(\cdot)$ and $\liminf$ for $\overline{v}(\cdot)$, and the nature of $\limsup\geq\liminf$, we have $\underline{v}({\mathbf{x}})\geq\overline{v}({\mathbf{x}})$ for all ${\mathbf{x}}\in\mathcal{S}$. If a comparison result holds in $\mathcal{S}$, it means that $\underline{v}({\mathbf{x}})\leq\overline{v}({\mathbf{x}})$ for all ${\mathbf{x}}\in\mathcal{S}$. Therefore, $v({\mathbf{x}})=\underline{v}({\mathbf{x}})=\overline{v}({\mathbf{x}})$ is the unique, continuous viscosity solution within the region $\mathcal{S}$.

It is established that the full-domain variational inequality defined in (2.4) satisfies a comparison result in [64, 59, 3, 40]. Similarly, the localized variational inequality (2.11) also satisfies a comparison result, as detailed in the lemma below. We recall $\Omega^{\scalebox{0.5}{$\infty$}}_{\scalebox{0.7}{\text{out}}}$ defined in (2.1).

The proof of the comparison result follows a similar approach to that in [51], and is therefore omitted here for brevity.

We let $G(\eta_{x},\eta_{y},\Delta\tau)$ be the Fourier transform of $g(x,y,\Delta\tau)$ with respect to the spatial variables, i.e.

A closed-form expression for $G(\eta_{x},\eta_{y},\cdot)$ is given as follows [62]

Here, $\Gamma(\eta_{x},\eta_{y})=\iint_{\mathbb{R}^{2}}f(s_{x},s_{y})~{}e^{i(\eta_{x}s_{x}+\eta_{y}s_{y})}~{}ds_{x}ds_{y}$, where $f(s_{x},s_{y})$ is the joint probability density function of random variables $\xi_{x}$ and $\xi_{y}$ given in (2.2).

For convenience, we define ${\boldsymbol{z}}=[x,y]$, ${\boldsymbol{\eta}}=[\eta_{x},\eta_{y}]$ and ${\boldsymbol{s}}=[s_{x},s_{y}]$ are the column vectors, ${\boldsymbol{z}}\cdot{\boldsymbol{\eta}}$ is the dot product of vectors ${\boldsymbol{z}}$ and ${\boldsymbol{\eta}}$, ${\boldsymbol{z}}^{\top}$ is the transpose of a vector ${\boldsymbol{z}}$, and $\tilde{\boldsymbol{C}}$ is the covariance matrix of $x$ and $y$. The covariance matrix $\tilde{\boldsymbol{C}}$ and its inverse $\tilde{\boldsymbol{C}}^{-1}$ are respectively given as follows

For subsequent use, we express the function $G(\eta_{x},\eta_{y},\cdot)$ given in (3.3) in a compact matrix-vector form as follows

where $\Gamma({\boldsymbol{\eta}})=\int_{\mathbb{R}^{2}}f({\boldsymbol{s}})~{}e^{i{\boldsymbol{s}}\cdot{\boldsymbol{\eta}}}~{}d{\boldsymbol{s}}$, and $\tilde{{\boldsymbol{\beta}}}=\left[(r-\lambda\kappa_{x}-\frac{\sigma_{\scalebox{0.6}{$x$}}^{2}}{2}),~{}(r-\lambda\kappa_{y}-\frac{\sigma_{\scalebox{0.6}{$y$}}^{2}}{2})\right]$ is the column vector. For brevity, we use $\int_{{\boldsymbol{\eta}}\in\mathbb{R}^{2}}(\cdot)~{}d{\boldsymbol{\eta}}$ to represent the 2-D integral $\iint_{\mathbb{R}^{2}}(\cdot)~{}d\eta_{x}d\eta_{y}$.

A proof of Lemma 3.1 is provided in Appendix B. We emphasize that the infinite series representation in Lemma 3.1 does not rely on the specific form of the joint probability density function $f(\cdot)$, and thus it applies broadly to general two-asset jump-diffusion model. In the specific case of the two-asset Merton jump-diffusion model, where the joint probability density function $f(\cdot)$ is given by (2.2), the terms of the series can be explicitly evaluated, as detailed in the corollary below.

A proof of Corollary 3.1 is given in Appendix C.

In the 2-D Merton model of Corollary 3.1, $[\log(\xi_{x}),\log(\xi_{y})]$ are assumed to follow a bivariate normal distribution. Consequently, the sum of $k$ i.i.d. random vectors, each following this bivariate normal distribution, is also bivariate normal. Convolution with the 2-D diffusion kernel $\exp\left(\theta-\frac{(\cdots)^{\top}\boldsymbol{C}^{-1}(\cdots)}{2}\right)$, which is itself Gaussian, then reduces to a straightforward Gaussian–Gaussian double integral. Summing over all possible numbers of jumps according to the Poisson distribution gives an infinite series of tractable Gaussian double integrals, thus leading to the closed‐form solution presented in Corollary 3.1.

However, extending this approach to the 2-D Kou model is much more challenging—even in the simplest case where $\log(\xi_{x})$ and $\log(\xi_{y})$ are independent. In this setting, the joint PDF of $\log(\xi_{x})$ and $\log(\xi_{y})$ is simply the product of two 1-D Kou PDFs, resulting in a four‐region, piecewise‐exponential density on $\mathbb{R}^{2}$. Summing $k$ i.i.d. random vectors, each following this independent 2-D Kou jump distribution, leads to multi‐region double integrals, since each dimension can exhibit positive or negative jumps. As $k$ grows, the number of piecewise‐exponential components increases combinatorially. Moreover, these sums produce factors that often combine polynomial or gamma‐type terms with exponentials in two variables.

When convolved with the 2-D Gaussian diffusion kernel, the resulting integrals become exponential–Gaussian double integrals, which generally do not reduce to standard special functions (unlike the 1-D case, which sometimes admits the $\mathrm{Hh}(\cdot)$ family from [1], as shown in [48]). Although one can still construct a series expansion over the Poisson‐distributed jumps, the individual terms would require intricate generalized special functions or numerous piecewise integrals, making the final expression far more cumbersome and less recognizable as a “closed‐form” solution.

A possible way forward is inspired by [28], where a single-hidden-layer neural network (NN) with Gaussian activation functions is used to approximate an unknown transition density or Green’s function via a closed-form expression of its Fourier transform. However, this approach can lead to potential loss of nonnegativity and requires retraining the neural network for different $\Delta\tau$.

For the general 2-D Kou model, an alternative approach is to instead approximate the joint PDF of $\log(\xi_{x})$ and $\log(\xi_{y})$ using a NN akin to that in [28], where Gaussian activation functions are employed. Once trained, the joint PDF is represented as a finite sum of 2-D Gaussians, effectively forming a Gaussian mixture model. Notably, the nonnegativity of this approximation can also be enforced, aligning with the monotonicity requirements of numerical schemes.²²2While we focus on the 2-D Kou model, this approach can, in principle, be applied to other 2-D jump-diffusion dynamics where the joint jump density can be effectively approximated by a Gaussian mixture.

This formulation transforms the convolution with the 2-D diffusion kernel into a tractable Gaussian–Gaussian double integral, allowing the same efficient techniques developed in this paper for the 2-D Merton model to be applied. Notably, because the NN does not depend on $\Delta\tau$, it can be trained only once and reused for all $\Delta\tau$, ensuring computational efficiency by eliminating the need for retraining at different time step sizes. We plan to explore this approach in future work in the context of fully 2-D asset allocation.

For subsequent analysis, we study the truncation error in the infinite series representation of the Green’s function $g(\cdot)$ as given in (3.6). Notably, this truncation error bound is derived independently of the specific form of the joint probability density function $f(\cdot)$, ensuring its applicability to a broad range of two-asset jump-diffusion models.

Specifically, for a fixed ${\boldsymbol{z}}=[x,y]\in\mathbb{R}^{2}$, we denote by $g({\boldsymbol{z}},\Delta\tau,K)$ an approximation of the Green’s function $g({\boldsymbol{z}},\Delta\tau)$ using the first $(K+1)$ terms from the series (3.6). As $K$ approaches $\infty$, the approximation $g({\boldsymbol{z}},\Delta\tau,K)$ becomes exact with no loss of information. However, with a finite $K$, the approximation incurs an error due to the truncation of the series. This truncation error can be bounded as follows:

Here, in (i), we apply the following fact: if $\omega$ denotes a complex number, then the modulus of the complex exponential is equivalent to the exponential of the real part of $\omega$, i.e $\left|e^{\omega}\right|=\exp(\Re(\omega))$ and $\left|\left(\Gamma\left({\boldsymbol{\eta}}\right)\right)^{K+1}\right|\leq\left(\int_{\mathbb{R}^{2}}f({\boldsymbol{s}})~{}\left|e^{i{\boldsymbol{s}}\cdot{\boldsymbol{\eta}}}\right|~{}d{\boldsymbol{s}}\right)^{K+1}\leq 1$, (ii) is due to the Chernoff-Hoeffding bound for the tails of a Poisson distribution $\mathrm{Poi}(\lambda\Delta\tau)$, which reads as $\mathbb{P}\left(\mathrm{Poi}(\lambda\Delta\tau)\geq k\right)\leq\frac{e^{-\lambda\Delta\tau}(e\lambda\Delta\tau)^{k}}{k^{k}}$, for $k>\lambda\Delta\tau$ [56].