This paper addresses a significant gap in rigorous numerical treatments for pricing American options under correlated two-asset jump-diffusion models using the viscosity solution approach, with a particular focus on the Merton model. The pricing of these options is governed by complex two-dimensional (2-D) variational inequalities that incorporate cross-derivative terms and nonlocal integro-differential terms due to the presence of jumps. Existing numerical methods, primarily based on finite differences, often struggle with preserving monotonicity in the approximation of cross-derivatives-a key requirement for ensuring convergence to the viscosity solution. In addition, these methods face challenges in accurately discretizing 2-D jump integrals. We introduce a novel approach to effectively tackle the aforementioned variational inequalities, seamlessly managing cross-derivative terms and nonlocal integro-differential terms through an efficient and straightforward-to-implement monotone integration scheme. Within each timestep, our approach explicitly tackles the variational inequality constraint, resulting in a 2-D Partial Integro-Differential Equation (PIDE) to solve. Its solution is then expressed as a 2-D convolution integral involving the Green's function of the PIDE. We derive an infinite series representation of this Green's function, where each term is strictly positive and computable. This series facilitates the numerical approximation of the PIDE solution through a monotone integration method, such as the composite quadrature rule. The proposed method is demonstrated to be both $\ell_{\infty} $-stable and consistent in the viscosity sense, ensuring its convergence to the viscosity solution of the variational inequality. Extensive numerical results validate the effectiveness and robustness of our approach, highlighting its practical applicability and theoretical soundness. 🔗 View origin paper >>