Pricing European and American options under Heston model using discontinuous Galerkin finite elements

Abstract

This paper deals with pricing of European and American options, when the underlying asset price follows Heston model, via the interior penalty discontinuous Galerkin finite element method (dGFEM). The advantages of dGFEM space discretization with Rannacher smoothing as time integrator with nonsmooth initial and boundary conditions are illustrated for European vanilla options, digital call and American put options. The convection dominated Heston model for vanishing volatility is efficiently solved utilizing the adaptive dGFEM. For fast solution of the linear complementary problem of the American options, a projected successive over relaxation (PSOR) method is developed with the norm preconditioned dGFEM. We show the efficiency and accuracy of dGFEM for option pricing by conducting comparison analysis with other methods and numerical experiments.

Introduction

In 1973, Black and Scholes [5] proposed a celebrated model under which the underlying price follows a geometric Brownian motion with constant volatility, and evaluate European options via an analytical formula. In spite of its simplicity and mathematical tractability, the Black–Scholes model is revealed to be unsuccessful, e.g., in predicting the volatility smirk due to the assumption of constant volatility. This important shortcoming can be overcame by regarding volatility as a source of randomness as in Heston model. Heston model treats the volatility as a square-root process, while also providing analytical prices for some options such as European vanilla options.

Heston model is a two-dimensional reaction–convection–diffusion (RCD) partial differential equation (PDE) with variable coefficients [20], [35]. The diffusion matrix contains the cross-diffusion term as a result of the correlation between the volatility and the underlying security. The commonly used method for option pricing under the Heston model is the finite difference method (FDM). It is relatively simple to implement, but the accuracy of the method is limited by time and space. For the effective numerical solution of these systems, FDMs with alternating direction implicit (ADI) time-stepping methods [29], [30], and high-order compact finite difference scheme in space [12] have been developed. When the boundary conditions, the shape of the domain, or data have limited regularity, finite elements methods (FEMs) perform better than FDMs [7], [52]. Radial basis [4] and spectral methods [41] are more accurate than the classical FDMs and FEMs, but they require the inversion of full system matrices. Fourier-based integration methods (COS) make use of the characteristic function, i.e., the Fourier transform of the probability density function of the underlying stock price of Heston model [20], which are fast and accurate [15]. Recently wavelet methods are developed using the B-splines basis [39] and Shannon wavelets inverse Fourier methods [40]. Due to the local nature of the wavelet basis, they are more robust and efficient than the COS methods. For American options, the value function satisfies a parabolic partial differential variational inequality system due to the early exercise constraint. FDMs are combined with operator splitting [8], [19], [27], [28], [38], [45] to the linear complementarity problem (LCP). FEM is also applied for solving American options under Heston model [16], [33], [55].

This paper deals with the numerical computation of European and American option pricing problems under the Heston model with discontinuous Galerkin finite elements method (dGFEM). The basis functions in dGFEM are discontinuous along the inter-element boundaries in contrast to the classical FEMs. The dGFEM has a number of desirable properties like the weakly enforcement of the boundary conditions, $hp$ (space and order) adaptivity and parallelization. In contrast to the stabilized continuous Galerkin finite elements methods for PDE with convection term like the Heston model, discontinuous Galerkin (dG) methods produce stable discretization without the need for extra stabilization strategies and damp the unphysical oscillations for convection dominated problems. The dG combines the best properties of the finite volume and continuous finite elements methods. Finite volume methods can only use lower degree polynomials, and continuous finite elements methods require higher regularity due to the continuity requirements. The disadvantage of the dG methods is higher condition numbers of the matrices of the associated algebraic systems [51] than for finite elements or finite differences. Several preconditioning techniques are developed as a remedy. One of them is applied here for American options.

We apply the symmetric interior penalty Galerkin (SIPG) method [2], [44] with upwinding for the convective part [3] to Heston model for some European options, such as vanilla and digital call options, and for American options. The SIPG method is the most used and popular method among the dG methods. On the other hand, other types of dGFEMs are applied to option pricing; pricing of European and American options with Constant Elasticity of Variance (CEV) model [37] and American options for the Black–Scholes equation [48], nonsymmetric interior penalty Galerkin (NIPG) method for pricing of European options under Black–Scholes [23], Heston models [24], and Asian options [25], [26]. An interesting feature of the Heston model is the occurrence of sharp layers or discontinuities for vanishing volatility. In the case of European call options, for instance, extremely high foreign interest rates make the problem convection dominated [30] when the underlying volatility is very small. In these cases, the naive approach is to refine the spatial mesh uniformly, which increases the degrees of freedom and refines the mesh unnecessarily in regions where the solutions are smooth. In some FDMs and FEMs [10], [29], [30], predefined nonuniform grids are used to resolve the discontinuities or the sharp layers accurately. In practice, the location of the interior or boundary layers for convection dominated problems are usually not known a priori. Adaptive methods can detect the layers using a posteriori error estimators and refine the mesh locally. Due to the local nature of the basis functions of the dGFEMs, the sharp layers and the singularities of the solution can be detected easily using the adaptive techniques [46]. We applied adaptive dGFEM based on a posteriori error estimator for an accurate and efficient solution of convection dominated Heston model for European options.

Option pricing models have nonsmooth initial data, with the discontinuous first derivatives of the payoff functions. The most popular time discretization method in option pricing is the Crank–Nicolson (CN) method, which leads undesired oscillations for nonsmooth initial data. The instability of the CN method is remedied by applying in the first four steps the implicit backward Euler (BE) method and then continuing with CN method as time integrator, known as Rannacher smoothing [43]. We apply the Rannacher smoothing for the Heston model with combination of SIPG and CN. The effect of the Rannacher smoothing is clearly visible for digital call options. For American options, the LCP is solved by projected successive over relaxation (PSOR) [27], [28] or by projected Gauss–Seidel (PSG) method. At each time step, a large full LCP has to be solved [8]. Because the condition number of the dGFEM discretized matrices increases rapidly for finer meshes, the convergence of the PSOR method slows down. Using the norm preconditioner in [17] designed for dGFEM discretization of RCD equations, we show that the convergence of PSOR method can be accelerated. Numerical results show that the number of iterations at each time step is essentially mesh independent. To the best of our knowledge, the PSOR with a matrix preconditioner is used for the first time in the evaluation of American options. We further present the detailed comparison of the SIPG for European and American options with radial basis functions with the partition of unity methods (RBF-PUM) [36], and with the NIPG [24].

The paper is organized as follows: In the next section, we introduce the Heston model for option pricing and give strong and variational forms of the underlying PDE. The space discretization via SIPG method and time discretization by CN method with Rannacher smoothing are described in Section 3. The American option pricing as a LCP is formulated in Section 4 with the preconditioned PSOR method. In Section 5, we give numerical orders of convergence for a test problem with a known solution, and we report on numerical results for European call as well as digital call and American put options, convection dominated problem for European call options using adaptive dGFEM. At the end of Section 5 we compare the dGFEM with the RBF-PUM methods for European and American options. The paper ends with some conclusions.