Multigrid Iterative Algorithms based on Compact Finite Difference Schemes and Hermite interpolation for Solving Regime Switching American Options
11 Multigrid Iterative Algorithms based on Compact Finite Difference Schemes and Hermite interpolation for Solving Regime Switching American Options
Chinonso I. Nwankwo a, *, Weizhong Dai b a Department of
Mathematics , Statistics, and Computer Science, University of Illinois at Chicago, Chicago, IL 60607, USA b Department of Mathematics and Statistics, Louisiana Tech University, Ruston LA 71272, USA * Corresponding author, [email protected], [email protected]
Abstract
We present multigrid iterative algorithms for solving a system of coupled free boundary problems for pricing American put options with regime-switching. The algorithms are based on our recent developed compact finite difference scheme coupled with Hermite interpolation for solving the ๐ coupled partial differential equations consisting of the asset, delta, gamma, and speed options. In the algorithms, we first use the Gauss-Seidel as a smoother, and then implement V-cycle and modified multigrid strategies for solving our discretized equations. Hermite interpolation with Newton interpolatory divided difference (as the basis) is used in estimating the coupled asset, delta, gamma, and speed options in the set of equations. A numerical experiment is performed with the two-regimes example and compared with other existing methods to validate the optimal strategy. Results show that these algorithms provide fast and efficient tools for pricing American put options with regime-switching. Keywords: multigrid methods,
American put options, Hermite and high order interpolation, optimal exercise boundary, compact finite difference method, logarithmic transformation Introduction
American style option valuation based on the regime-switching model involves ๐ coupled free boundary problem where ๐ represent the number of regimes. The closed-form of this model is difficult to find, hence, several numerical methods have been proposed (Khaliq and Liu, 2009; Nielsen et al., 2002; Zhang et al., 2013, Chiarella et al., 2016; Meyer and van der Hoek, 1997, Han and Kim, 2016; Egorova et al., 2016; Shang and Bryne, 2019). Moreover, it is always a challenging task to approximate the Greeks associated with this model especially, beyond two regimes. In our previous work (Nwankwo et al., 2019), we generated a system of ๐ coupled partial differential equations through logarithmic transformation and by taking further derivatives of the regime-switching model. We called this system of ๐ equations, the ๐ asset-delta-gamma-speed equations. By employing the compact finite difference method for discretizing the system of ๐ coupled equation and the Hermite interpolation for estimating the coupled regimes, we were able to approximate the ๐ asset options, option Greeks, and optimal exercise boundaries with high order accuracy using Gauss-Seidel iterative methods. it is imperative to mention that for each iteration using a front-fixing approach for solving the regime-switching model, we approximate the coupled regime with Hermite Interpolation which increases the overall computational burden. When a small step size and time step are used, this effect is more pronounced. Hence, there is a need to find a strategy that reduces the number of iterations and CPU time required to achieve numerical convergence. The motivation of this article is to develop a multigrid iterative algorithm to speed up our computation for solving the system of ๐ equations. As we know, the major advantage of the multigrid iterative method is that it can remove the high-frequency error from numerical solutions in a few iterations (using a good smoother) and transfers the low-frequency error to a coarse grid where the latter is cheaper and faster to remove. This improves the overall computational effort, speed, and efficiency of numerical approximations for large scale problems (Hafner and Konke, 2006). Because of this, the multigrid method has gained a broader application after the work of Fedorenko (1961), Bachvalov (1966), Brandt (1977), and Saad (2003) which has been applied to elliptic and time-dependent partial differential equations (Wang and Ge, 2018; Rosseel et al., 2007; Horton and Vandewalle, 1994; Lee and Meyer, 1979; Yavneh, 2006; Khelifi et al., 2014). The efficiency of the multigrid method with a compact finite difference scheme has further been investigated by researchers (Wang and Ge, 2018; Pardhanani et al., (1997); Zhang et al., (2002); Moghaderi and Dehghan, 2014; Ghaffar et al., 2016). Furthermore, Shieh (1985) and Guo et al. (2012) applied a multigrid method to obtain a numerical solution of the coupled partial differential equation based on domain partitioning. In option pricing, Urschel (2013) implemented an adaptive space-time multigrid approach to barrier option using implicit Euler and Crank-Nicholson method. Jeong et al. (2013) applied an adaptive multigrid method on the Black Scholes equation with Crank-Nicholson discretization. Clarke and Parrott (1999) implemented a multigrid method in discrete linear complementary problem for pricing American options with stochastic volatility. In this study, we will implement the multigrid method to a compact finite difference scheme coupled with the Hermite interpolation for pricing American options with regime-switching. The rest of the paper is organized as follows. In section 2, we present our mathematical model. In section 3, we present the compact finite difference scheme and estimate the coupled asset, delta, gamma, and speed options in the set of ๐ equations using the Hermite interpolation with the Newton interpolatory divided difference formula as the basis. In section 4, we implement multigrid strategies and their algorithms for solving the derived numerical scheme and obtain the option values, optimal exercise boundary, and the Greeks in each regime. In section 5, we investigate and compare the numerical performance of our present algorithms using the two-regimes example and conclude the paper in section 6. Mathematical Model
Our mathematical model is based on the American put option written on the asset ๐ ๐ก with strike price ๐พ and expiration time ๐. Let ๐ ๐ (๐, ๐ก) denote the option price and ๐ = ๐ โ ๐ก . Then, ๐ ๐ (๐, ๐) satisfies the coupled free boundary value problem: โ ๐๐ ๐ (๐, ๐)๐๐ + 12 ๐ ๐ ๐ ๐ ๐ (๐, ๐)๐๐ + ๐ ๐ ๐ ๐๐ ๐ (๐, ๐)๐๐ โ ๐ ๐ ๐ ๐ (๐, ๐) + โ ๐ ๐๐๐โ ๐ [๐ ๐ (๐, ๐) โ ๐ ๐ (๐, ๐)] = 0,for ๐ > ๐ ๐(๐) (๐), (1) ๐ ๐ (๐, ๐) = ๐พ โ ๐, for ๐ < ๐ ๐(๐) (๐). (2) Here, the initial and boundary conditions are given as: ๐ ๐ (๐, 0) = ๐๐๐ฅ( ๐พ โ ๐, 0), ๐ ๐(๐) (0) = ๐พ; (3๐) ๐ ๐ (๐ ๐(๐) , ๐) = ๐พ โ ๐ ๐(๐) (๐), ๐ ๐ (0, ๐) = ๐พ, ๐ ๐ (โ, ๐) = 0, ๐๐๐ ๐ ๐ (๐ ๐(๐) , ๐) = โ1, (3๐) where ๐ ๐(๐) (๐) is the optimal exercise boundary for the ๐ ๐กโ regime. We further fix the free boundary challenge by employing a front-fixing logarithmic transformation (Egorova, 2016) on multi-variable domains as ๐ฅ ๐ = ln ๐๐ ๐(๐) (๐) = ln ๐ โ ln ๐ ๐(๐) (๐) , ๐ = 1,2,โโโ, ๐ผ, (4๐) The transformed ๐ option value functions ๐ ๐ (๐ฅ ๐ , ๐) are related to the original ๐ option value functions ๐ ๐ (๐, ๐ก) by the dimensionless transformation ๐ ๐ (๐ฅ ๐ , ๐) = ๐ ๐ (๐, ๐), ๐ = 1,2,โโโ, ๐ผ. (4๐) Using this transformation and eliminating the first-order derivative by taking further derivatives, we then obtain a set of ๐ coupled partial differential equations, which we call the asset, delta, gamma, and speed PDE equations and are given as follows: ๐๐ ๐ ๐๐ โ 12 ๐ ๐ ๐ ๐ ๐๐ฅ ๐2 โ (๐ โฒ๐(๐) ๐ ๐(๐) + ๐ ๐ โ ๐ ๐ + (๐ ๐ โ ๐ ๐๐ )๐ ๐ โ โ ๐ ๐๐๐โ ๐ ๐ ๐ = 0, (5๐) ๐๐ ๐ ๐๐ โ 12 ๐ ๐ ๐ ๐ ๐๐ฅ ๐2 โ (๐ โฒ๐(๐) ๐ ๐(๐) + ๐ ๐ โ ๐ ๐ ๐ ๐๐ฅ ๐2 + (๐ ๐ โ ๐ ๐๐ )๐ ๐ โ โ ๐ ๐๐๐โ ๐ ๐ ๐ = 0, (5๐) ๐๐ ๐ ๐๐ โ 12 ๐ ๐ ๐ ๐ ๐๐ฅ ๐2 โ (๐ โฒ๐(๐) ๐ ๐(๐) + ๐ ๐ โ ๐ ๐ ๐ ๐๐ฅ ๐2 + (๐ ๐ โ ๐ ๐๐ )๐ ๐ โ โ ๐ ๐๐๐โ ๐ ๐ ๐ = 0, (5๐) ๐๐ ๐ ๐๐ โ 12 ๐ ๐ ๐ ๐ ๐๐ฅ ๐2 โ (๐ โฒ๐(๐) ๐ ๐(๐) + ๐ ๐ โ ๐ ๐ ๐ ๐๐ฅ ๐2 + (๐ ๐ โ ๐ ๐๐ )๐ ๐ โ โ ๐ ๐๐๐โ ๐ ๐ ๐ = 0, (5๐) where ๐ = 1,2,โโโ, ๐ผ, ๐ฅ ๐ โ [0, โ) and the initial and boundary conditions for ๐ ๐ (๐ฅ ๐ , ๐) , ๐ ๐ (๐ฅ ๐ , ๐) , ๐ ๐ (๐ฅ ๐ , ๐) , and ๐ ๐ (๐ฅ ๐ , ๐) are defined as : ๐ ๐(๐) (0) = ๐พ, ๐ ๐ (๐ฅ ๐ , 0) = 0, ๐ ๐ (๐ฅ ๐ , 0) = 0, ๐ ๐ (๐ฅ ๐ , 0) = 0, ๐ ๐ (๐ฅ ๐ , 0) = 0; (6๐) ๐ ๐ (0, ๐) = ๐พ โ ๐ ๐(๐) (๐), ๐ ๐ (0, ๐) = โ๐ ๐(๐) (๐); (6๐) ๐ ๐ (0, ๐) = โ๐ ๐(๐) (๐), ๐ ๐ (0, ๐) = โ๐ ๐(๐) (๐); (6๐) ๐ ๐ (โ, ๐) = 0; ๐ ๐ (โ, ๐) = 0; ๐ ๐ (โ, ๐) = 0, ๐ ๐ (โ, ๐) = 0. (6๐) Numerical Discretization and Interpolation.
To solve the above system of ๐ PDEs that consist of the asset, delta, gamma and speed options in a uniform grid [0, โ) ร [0 ๐] for each regime and take into account, the relationship among the regimesโ intervals using cubic interpolation, the infinite domain is truncated to a finite large domain with an estimated boundary (๐ฅ ๐ ) ๐ far enough that the error generated due to the truncation is negligible (Egorova et al., 2016; Kangro and Nicolaides, 2000) . Letting ๐ and ๐ represent the node points in the ๐ ๐กโ and ๐ ๐กโ regimesโ intervals, respectively, and M and N represent the numbers of grid points and time steps, respectively, we have (๐ฅ ๐ ) ๐ = ๐โ, (๐ฅ ๐ ) ๐ = ๐โ, โ = (๐ฅ ๐ ) ๐ ๐ = (๐ฅ ๐ ) ๐ ๐ , ๐, ๐ โ [0, M]; (7๐) ๐ = ๐๐ , ๐ ๐ = ๐๐, ๐ โ [0, N]. (7๐) It should be pointed out that we choose the same length of interval for all regimes so that we have a uniform grid size.
In this section, we find the numerical solutions of the ๐ asset options, option Greeks, and optimal exercise boundaries which we denote as (๐ข ๐ ) ๐๐ , (๐ค ๐ ) ๐๐ , (๐ฆ ๐ ) ๐๐ , (๐ง ๐ ) ๐๐ , (ฮ ๐ ) ๐๐ , (ฮ ๐ ) ๐๐ , (ฮ ๐ ) ๐๐ , and ๐ ๐(๐)๐ . Compact Finite Difference Scheme
In our previous study , the set of ๐ equations in (5) were discretized using compact and Crank-Nicolson scheme in space and time, respectively. To develop a compact finite difference scheme in space for the ๐ asset option at (๐ฅ ๐ ) = 0 , we first employed a compact finite difference formula as described in the following lemma. Lemma.
Assume ๐(๐ฅ) โ ๐ถ [๐ฅ , ๐ฅ ] , then it holds
74 ๐ โฒโฒ (๐ฅ ) + 34 ๐ โฒโฒ (๐ฅ ) = 5โ [๐(๐ฅ ) โ ๐(๐ฅ )] โ 5โ ๐ โฒ (๐ฅ ) โ โ4 ๐ (3) (๐ฅ ) + โ6 ๐ (3) (๐ฅ ) + ๐(โ ). (8) Proof.
Along the lines of Nwankwo et al. (2019) and Hirsh (1975). Using (8) for the second-order derivative term in (1a), we obtained โ 12 ๐ [74 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐2 + 34 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐2 ]= 5๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) โ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )โ โ 1โ ๐๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐ ]โ ๐ โ2 [14 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐3 + 16 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐3 ] + ๐(โ ). (9) The first-order derivative in (9) at (๐ฅ ๐ ) was then evaluated based on (6b) as ๐๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐ โ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) = ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) โ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) = โ๐พ. (10) To evaluate the third-order derivative term in (9) at (๐ฅ ๐ ) , we let ๐ฅ ๐ โ 0 + in (1b) and discretized ๐๐ ๐ /๐๐ก. This gives ๐ โ8 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐3 = โ4 ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1 ) โ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐ )๐ โ โ(๐ฝ ๐ ) ๐+12 ๐ ((๐ฅ ๐ ) , ๐ ๐+12 ) + โ( ๐ ๐ โ ๐ ๐๐ )4 ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) โ โ4 โ ๐ ๐๐๐โ ๐ ๐ ๐ ((๐ฅ ๐ ) ๐ โ |๐=0 , ๐ ๐+1/2 ) + ๐(๐ ), (11๐) and similarly, ๐ โ12 ๐ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 )๐๐ฅ ๐3 = โ4 ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1 ) โ ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐ )๐ โ โ(๐ฝ ๐ ) ๐+1/2 ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) + โ( ๐ ๐ โ ๐ ๐๐ )4 ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) โ โ4 โ ๐ ๐๐๐โ ๐ ๐ ๐ ((๐ฅ ๐ ) ๐ โ |๐=1 , ๐ ๐+1/2 ) + ๐(๐ ), (11๐) where (๐ฝ ๐ ) ๐+1/2 โก 2(๐ ๐(๐)๐+1 โ ๐ ๐(๐)๐ )๐(๐ ๐(๐)๐+1 + ๐ ๐(๐)๐ ) + ๐ ๐ โ ๐ ๐ = ๐ ๐โ . (12) Here, (๐ฅ ๐ ) ๐ โ |๐=0 and (๐ฅ ๐ ) ๐ โ |๐=1 are the locations in the space for the ๐ ๐กโ equation corresponding to (๐ฅ ๐ ) and (๐ฅ ๐ ) in the ๐ ๐กโ equation, respectively. For the term ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) and ๐ ๐ ((๐ฅ ๐ ) , ๐ ๐+1/2 ) in (1b), we employed a fourth-order approximation (Adam, 1975; Adam, 1977, Liao and Khaliq, 2009; Dremkova and Ehrhardt, 2011) as ๐ ๐ ((๐ฅ ๐ ) ๐ , ๐ ๐+1/2 )= 3โ [๐ ๐ ((๐ฅ ๐ ) ๐+2 , ๐ ๐+1/2 ) โ ๐ ๐ ((๐ฅ ๐ ) ๐ , ๐ ๐+1/2 )] โ 4๐ ๐ ((๐ฅ ๐ ) ๐+1 , ๐ ๐+1/2 )โ ๐ ๐ ((๐ฅ ๐ ) ๐+2 , ๐ ๐+1/2 ) + ๐(โ ). (13) Substituting (11) and (13) into (9) and applying the Crank-Nicholson method in time , we then obtained the compact finite difference scheme at (๐ฅ ๐ ) = 0 as ๐ ๐1 (๐ข ๐ ) + ๐ ๐1 (๐ข ๐ ) = (๐ ) ๐๐+1/2 , (14) with the truncation error of ๐(๐ + โ ) . Here, ๐ โ |๐ = 0 ,1 indicate the values of ๐ โ given at (๐ฅ ๐ ) and (๐ฅ ๐ ) , respectively. At each interior grid point, (๐ฅ ๐ ) ๐ = 1,2, โฆ , ๐ โ 1 , using the compact finite difference scheme (Yan et al., 2019) as
112 ๐ โฒโฒ (๐ฅ ๐โ1 ) + 1012 ๐ โฒโฒ (๐ฅ ๐ ) + 112 ๐ โฒโฒ (๐ฅ ๐+1 ) = 1โ [๐(๐ฅ ๐โ1 ) โ 2๐(๐ฅ ๐ ) + ๐(๐ฅ ๐+1 )] + ๐(โ ), (15) for (1a)-(1d), we obtained ๐ ๐1 (๐ข ๐ ) ๐โ1๐+1 + ๐ ๐1 (๐ข ๐ ) ๐๐+1 + ๐ ๐1 (๐ข ๐ ) ๐+1๐+1 = (๐ ๐๐ข ) ๐๐+1/2 , (16๐) ๐ ๐1 (๐ค ๐ ) ๐โ1๐+1 + ๐ ๐1 (๐ค ๐ ) ๐๐+1 + ๐ ๐1 (๐ค ๐ ) ๐+1๐+1 = (๐ ๐๐ค ) ๐๐+1/2 , (16๐) ๐ ๐1 (๐ฆ ๐ ) ๐โ1๐+1 + ๐ ๐1 (๐ฆ ๐ ) ๐๐+1 + ๐ ๐1 (๐ฆ ๐ ) ๐+1๐+1 = (๐ ๐๐ฆ ) ๐๐+1/2 , (16๐) ๐ ๐1 (๐ง ๐ ) ๐โ1๐+1 + ๐ ๐1 (๐ง ๐ ) ๐๐+1 + ๐ ๐1 (๐ง ๐ ) ๐+1๐+1 = (๐ ๐๐ง ) ๐๐+1/2 , ๐ = 1,2, โฆ , ๐ โ 1. (16๐) To further ensure that our coefficient matrix entries preserve both symmetricity and positive definiteness, we multiply (14) by ๐ ๐1 / ๐ ๐1 and obtain ๐ฬ ๐1 (๐ข ๐ ) + ๐ ๐1 (๐ข ๐ ) = (๐ฬ ) ๐๐+1/2 , (16๐) where ๐ฬ ๐1 = ๐ ๐1 ๐ ๐1 ๐ ๐1 , ๐ฬ = ๐ ๐1 (๐ ) ๐๐+1/2 ๐ ๐1 . (16๐) Hence, we obtain a discretized system for the asset, delta, gamma, and speed option equations as follows: ๐ด ๐ ๐ = ๐ ๐๐ข , ๐ต ๐ ๐ = ๐ ๐๐ค , ๐ต ๐ ๐ = ๐ ๐๐ฆ , ๐ต ๐ ๐ = ๐ ๐๐ง , (17) where ๐ด ๐ = [ ๐ฬ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 โฑ โฑ โฑ ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ] ๐ต ๐ = [ ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 โฑ โฑ โฑ ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ๐ ๐1 ] ๐ ๐๐ข = [ (๐ฬ ) ๐๐+1/2 (๐ ) ๐๐+1/2 (๐ ) ๐๐+1/2 ...(๐ ๐โ1๐ข ) ๐๐+1/2 ] ๐ ๐๐ค = [ (๐ ) ๐๐+1/2 โ ๐ ๐1 (๐ค ๐ ) (๐ ) ๐๐+1/2 (๐ ) ๐๐+1/2 ...(๐ ๐โ1๐ค ) ๐๐+1/2 ] โฆ ๐ ๐๐ง = [ (๐ ) ๐๐+12 โ ๐ ๐1 (๐ง ๐ ) (๐ ) ๐๐+12 (๐ ) ๐๐+12 ...(๐ ๐โ1๐ง ) ๐๐+12 ] . (18) Here, (๐ ) ๐๐+1/2 = ๐ ๐2 (๐ข ๐ ) + ๐ ๐2 (๐ข ๐ ) + ๐พ๐๐ ๐ + โ6 [(๐ค ๐ ) โ (๐ค ๐ ) ] + โ๐(๐ ๐ โ ๐ ๐๐ )12 [(๐ค ๐ ) + (๐ค ๐ ) ]+ ๐(๐ฝ ๐ ) ๐+1/2 ๐ ) + (๐ค ๐ ) ] + 15[(๐ค ๐ ) + (๐ค ๐ ) ] + 9[(๐ค ๐ ) + (๐ค ๐ ) ]โ 6[(๐ค ๐ ) + (๐ค ๐ ) ]]โ โ๐(๐ฝ ๐ ) ๐+1/2
24 [14[(๐ฆ ๐ ) + (๐ฆ ๐ ) ] โ 5[(๐ฆ ๐ ) + (๐ฆ ๐ ) ] โ 2[(๐ฆ ๐ ) + (๐ฆ ๐ ) ]] + โ๐24 โ ๐ ๐๐๐โ ๐ [3((๐ค ๐ ) ๐ โ |๐=0๐+1 + (๐ค ๐ ) ๐ โ |๐=0๐ ) โ 2((๐ค ๐ ) ๐ โ |๐=1๐+1 + (๐ค ๐ ) ๐ โ |๐=1๐ )]+ ๐8 โ ๐ ๐๐ [7((๐ข ๐ ) ๐ โ |๐=0๐+1 + (๐ข ๐ ) ๐ โ |๐=0๐ ) + 3((๐ข ๐ ) ๐ โ |๐=1๐+1 + (๐ข ๐ ) ๐ โ |๐=1๐ )], ๐โ ๐ (19๐) (๐ ๐๐ข ) ๐๐+1/2 = ๐ ๐2 (๐ข ๐ ) ๐โ1๐ + ๐ ๐2 (๐ข ๐ ) ๐๐ + ๐ ๐2 (๐ข ๐ ) ๐+1๐ + ๐ ๐ (๐ฝ ๐ ) ๐+1/2
24 [[(๐ค ๐ ) ๐โ1๐+1 + (๐ค ๐ ) ๐โ1๐ ] + 10[(๐ค ๐ ) ๐๐+1 + (๐ค ๐ ) ๐๐ ] + [(๐ค ๐ ) ๐+1๐+1 + (๐ค ๐ ) ๐+1๐ ]]+ ๐24 โ ๐ ๐๐ [[(๐ข ๐ ) ๐ โ โ1๐+1 + (๐ข ๐ ) ๐ โ โ1๐ ] + 10[(๐ข ๐ ) ๐ โ ๐+1 + (๐ข ๐ ) ๐ โ ๐ ] ๐โ ๐ + [(๐ข ๐ ) ๐ โ +1๐+1 + (๐ข ๐ ) ๐ โ +1๐ ]] , (19๐) (๐ ๐๐ค ) ๐๐+1/2 = ๐ ๐2 (๐ค ๐ ) ๐โ1๐ + ๐ ๐2 (๐ค ๐ ) ๐๐ + ๐ ๐2 (๐ค ๐ ) ๐+1๐ + ๐ ๐ ๐ ) ๐+1/2 [[(๐ข ๐ ) ๐โ1๐+1 + (๐ข ๐ ) ๐โ1๐ ] โ 2[(๐ข ๐ ) ๐๐+1 + (๐ข ๐ ) ๐๐ ] + [(๐ข ๐ ) ๐+1๐+1 + (๐ข ๐ ) ๐+1๐ ]]+ โ ๐ ๐๐ ๐โ ๐ [[(๐ค ๐ ) ๐ โ โ1๐+1 + (๐ค ๐ ) ๐ โ โ1๐ ] + 10[(๐ค ๐ ) ๐ โ ๐+1 + (๐ค ๐ ) ๐ โ ๐ ] + [(๐ค ๐ ) ๐ โ +1๐+1 + (๐ค ๐ ) ๐ โ +1๐ ]] , (19๐) (๐ ๐๐ฆ ) ๐๐+1/2 = ๐ ๐2 (๐ฆ ๐ ) ๐โ1๐ + ๐ ๐2 (๐ฆ ๐ ) ๐๐ + ๐ ๐2 (๐ฆ ๐ ) ๐+1๐ + ๐ ๐ ๐ ) ๐+1/2 [[(๐ค ๐ ) ๐โ1๐+1 + (๐ค ๐ ) ๐โ1๐ ] โ 2[(๐ค ๐ ) ๐๐+1 + (๐ค ๐ ) ๐๐ ] + [(๐ค ๐ ) ๐+1๐+1 + (๐ค ๐ ) ๐+1๐ ]]+ ๐24 โ ๐ ๐๐๐โ ๐ [[(๐ฆ ๐ ) ๐ โ โ1๐+1 + (๐ฆ ๐ ) ๐ โ โ1๐ ] + 10[(๐ฆ ๐ ) ๐ โ ๐+1 + (๐ฆ ๐ ) ๐ โ ๐ ] + [(๐ฆ ๐ ) ๐ โ +1๐+1 + (๐ฆ ๐ ) ๐ โ +1๐ ]] , (19๐) (๐ ๐๐ง ) ๐๐+1/2 = ๐ ๐2 (๐ง ๐ ) ๐โ1๐ + ๐ ๐2 (๐ง ๐ ) ๐๐ + ๐ ๐2 (๐ง ๐ ) ๐+1๐ + ๐ ๐ ๐ ) ๐+1/2 [[(๐ฆ ๐ ) ๐โ1๐+1 + (๐ฆ ๐ ) ๐โ1๐ ] โ 2[(๐ฆ ๐ ) ๐๐+1 + (๐ฆ ๐ ) ๐๐ ] + [(๐ฆ ๐ ) ๐+1๐+1 + (๐ฆ ๐ ) ๐+1๐ ]]+ ๐24 โ ๐ ๐๐๐โ ๐ [[(๐ง ๐ ) ๐ โ โ1๐+1 + (๐ง ๐ ) ๐ โ โ1๐ ] + 10[(๐ง ๐ ) ๐ โ ๐+1 + (๐ง ๐ ) ๐ โ ๐ ] + [(๐ง ๐ ) ๐ โ +1๐+1 + (๐ง ๐ ) ๐ โ +1๐ ]], (19๐) ๐ ๐1 = 7 + โ4 + 54 ๐ ๐ + 5โ4 ๐ ๐ + (7 + โ)๐8 (๐ ๐ โ ๐ ๐๐ ), ๐ ๐1 = 34 โ 54 ๐ ๐ + 3๐8 (๐ ๐ โ ๐ ๐๐ ); (20๐) ๐ ๐1 = 1012 + ๐ ๐ ๐ โ ๐ ๐๐ ), ๐ ๐1 = 112 โ ๐ ๐ ๐ โ ๐ ๐๐ ); (20๐) ๐ ๐2 = 7 + โ4 โ 54 ๐ ๐ โ 5โ4 ๐ ๐ โ (7 + โ)๐8 (๐ ๐ โ ๐ ๐๐ ), ๐ ๐2 = 34 + 54 ๐ ๐ โ 3๐8 (๐ ๐ โ ๐ ๐๐ ); (20๐) ๐ ๐2 = 1012 โ ๐ ๐ ๐ โ ๐ ๐๐ ), ๐ ๐2 = 112 + ๐ ๐ ๐ โ ๐ ๐๐ ) ; (20๐) ๐ ๐ ๐ โ ๐ ๐๐ )4 + 5๐
2โ . (20๐) where ๐ โ represents the location for the ๐ ๐กโ regime corresponding to (๐ฅ ๐ ) ๐ , and the truncation error is ๐(๐ + โ ) . The optimal exercise boundary and the initial and boundary conditions for each regime are calculated as ๐ ๐(๐) ๐+1 = ๐พ โ (๐ข ๐ ) , (๐ค ๐ ) = โ๐ ๐(๐) ๐+1 , (๐ฆ ๐ ) = โ๐ ๐(๐) ๐+1 , (๐ง ๐ ) = โ๐ ๐(๐) ๐+1 ; (21 ) (๐ข ๐ ) ๐๐+1 = 0, (๐ค ๐ ) ๐๐+1 = 0, (๐ฆ ๐ ) ๐๐+1 = 0, (๐ง ๐ ) ๐๐+1 = 0; (22 ) (๐ข ๐ ) ๐0 = (๐ค ๐ ) ๐0 = (๐ฆ ๐ ) ๐0 = (๐ง ๐ ) ๐0 = 0, ๐ = 1,2,โโโ, ๐. (23) Let the approximate solutions of the theta, delta decay, and color options for each regime be given as ๐๐ ๐ ((๐ฅ ๐ ) ๐ , ๐ ๐ )๐๐ โ (ฮ ๐ ) ๐๐ , ๐๐ ๐ ((๐ฅ ๐ ) ๐ , ๐ ๐ )๐๐ โ (ฮ ๐ ) ๐๐ , ๐๐ ๐ ((๐ฅ ๐ ) ๐ , ๐ ๐ )๐๐ โ (ฮ ๐ ) ๐๐ ; (24) respectively. For ๐ = 1 , we approximate these three Greeks using first-order backward finite differences (ฮ ๐ ) ๐1 โ (๐ข ๐ ) ๐1 โ (๐ข ๐ ) ๐0 ๐ , (ฮ ๐ ) ๐1 โ (๐ค ๐ ) ๐1 โ (๐ค ๐ ) ๐0 ๐ , (ฮ ๐ ) ๐1 โ (๐ฆ ๐ ) ๐1 โ (๐ฆ) ๐0 ๐ . (25๐) Subsequently, we use the second-order backward finite difference approximations as (ฮ ๐ ) ๐๐+1 = 3(๐ข ๐ ) ๐๐+1 โ 4(๐ข ๐ ) ๐๐ + (๐ข ๐ ) ๐๐โ1
2๐ , (ฮ ๐ ) ๐๐+1 = 3(๐ค) ๐๐+1 โ 4(๐ค ๐ ) ๐๐ + (๐ค ๐ ) ๐๐โ1
2๐ ; (25๐) (ฮ ๐ ) ๐๐+1 = 3(๐ฆ) ๐๐+1 โ 4(๐ฆ ๐ ) ๐๐ + (๐ฆ ๐ ) ๐๐โ1
2๐ . (25๐)
The initial conditions of the theta, delta decay and color options for each regime are calculated as (ฮ ๐ ) ๐0 = 0, (ฮ ๐ ) ๐0 = 0, (ฮ ๐ ) ๐0 = 0, ๐ = 0,1,โโโ, ๐. (26) Hermite Interpolation
The relationship between the fixed interval (and the mesh) for the ๐ ๐กโ the regime and the fixed interval (and the mesh) for the ๐ ๐กโ regime after the logarithmic transformation is as follows: (๐ฅ ๐ ) ๐ = (๐ฅ ๐ ) ๐ โ ln ๐ ๐(๐) (๐ ๐ )๐ ๐(๐) (๐ ๐ ) . (27) If ๐ ๐(๐) (๐ ๐ ) = ๐ ๐(๐) (๐ ๐ ) , then the fixed interval for the ๐ ๐กโ regime overlaps completely with the fixed interval for the ๐ ๐กโ regime. Hence, (๐ฅ ๐ ) ๐ = (๐ฅ ๐ ) ๐ and we solve only one interval. However, if ๐ ๐(๐) (๐ ๐ ) โ ๐ ๐(๐) (๐ ๐ ), we need to solve on multi intervals and there will exist three cases. We refer the reader to the work of Egorova et al. (2016) for further reading. We only focus on the case when (๐ฅ ๐ ) ๐ < (๐ฅ ๐ ) ๐ โ = (๐ฅ ๐ ) ๐ โ ln ๐ ๐(๐) (๐ ๐ )๐ ๐(๐) (๐ ๐ ) < (๐ฅ ๐ ) ๐+1 , (28) which requires interpolating between nodes. Here (๐ฅ ๐ ) ๐ โ โ [(๐ฅ ๐ ) ๐ , (๐ฅ ๐ ) ๐+1 ] and (๐ข ๐ ) ๐ โ ๐ , (๐ค ๐ ) ๐ โ ๐ , (๐ฆ ๐ ) ๐ โ ๐ and (๐ง ๐ ) ๐ โ ๐ need to be evaluated using high order interpolation and based on (๐ข ๐ ) ๐๐ , (๐ค ๐ ) ๐๐ , (๐ฆ ๐ ) ๐๐ and (๐ง ๐ ) ๐๐ . We employ the Hermite interpolation with a Newton basis to estimate these coupled regimes in the set of ๐ equations. The choice of using the Newton interpolatory divided difference formula as the basis provides faster convergence to our numerical solution than the initial choice of Lagrange basis in our earlier work. For the readerโs interest, please see the work of Powell (1981), Egecioglu et al. (1989), and Burden et al. (2016) for further explanation and proof of this interpolation technique. Approximating the coupled regimes in the set of ๐ equations using the Hermite interpolation with Newton basis is described as follows. Let ๐ โ โ [๐, ๐ + 1] be the point in the interval of ๐ ๐กโ regime that we need to approximate the asset, delta, gamma, and speed function. Hence, (๐ข ๐ ) ๐ โ ๐ = (๐ข ๐ ) ๐๐ + (๐ค ๐ ) ๐๐ [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ]โ + [(๐ข ๐ ) ๐+1๐ โ (๐ข ๐ ) ๐๐ โ (๐ค ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] โ+ [(๐ค ๐ ) ๐+1๐ โ (๐ค ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐+1 ]โ , (29๐) (๐ค ๐ ) ๐ โ ๐ = (๐ค ๐ ) ๐๐ + (๐ฆ ๐ ) ๐๐ [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ]โ + [(๐ค ๐ ) ๐+1๐ โ (๐ค ๐ ) ๐๐ โ (๐ฆ ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] โ+ [(๐ฆ ๐ ) ๐+1๐ โ (๐ฆ ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐+1 ]โ โฒ (29๐) (๐ฆ ๐ ) ๐ โ ๐ = (๐ฆ ๐ ) ๐๐ + (๐ง ๐ ) ๐๐ [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ]โ + [(๐ฆ ๐ ) ๐+1๐ โ (๐ฆ ๐ ) ๐๐ โ (๐ง ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] โ+ [(๐ง ๐ ) ๐+1๐ โ (๐ง ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐+1 ]โ , (29๐) (๐ง ๐ ) ๐ โ ๐ = (๐ง ๐ ) ๐๐ + (๐งฬ ๐ ) ๐๐ [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ]โ + [(๐ง ๐ ) ๐+1๐ โ (๐ง ๐ ) ๐๐ โ (๐งฬ ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] โ+ [(๐งฬ ๐ ) ๐+1๐ โ (๐งฬ ๐ ) ๐๐ ][(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐ ] [(๐ฅ ๐ ) ๐ โ โ (๐ฅ ๐ ) ๐+1 ]โ . (29๐) Moreover, it is worth noting that the derivative of the speed option, (๐งฬ ๐ ) ๐๐ , is employed in the cubic Hermite interpolations. To evaluate (๐งฬ ๐ ) ๐๐ with fourth-order accuracy, we obtain it from the speed option PDE by further taking derivative. Multigrid Method
Generally, the process of approximating the ๐ discrete system of the asset, delta, gamma, and speed equations in (17) involves the following steps. First, we start our iteration on a fine uniform grid 1 ๐ด ๐,โ ๐ โ = ๐ ๐,โ๐ข , ๐ต ๐,โ ๐ โ = ๐ ๐,โ๐ค , ๐ต ๐,โ ๐ โ = ๐ ๐,โ๐ฆ , ๐ต ๐,โ ๐ โ = ๐ ๐,โ๐ง , (30) where h is the step size of the finest grid. Next, we relax (30) ๐ times using an iterative method that is a good smoother and compute their residual as follows: ๐ ๐,โ๐ข = ๐ ๐,โ๐ข โ ๐ด ๐,โ ๐ โ , ๐ ๐,โ๐ค = ๐ ๐,โ๐ค โ ๐ด ๐,โ ๐ โ , (31๐) ๐ ๐,โ๐ฆ = ๐ ๐,โ๐ฆ โ ๐ด ๐,โ ๐ โ , ๐ ๐,โ๐ง = ๐ ๐,โ๐ง โ ๐ด ๐,โ ๐ โ . (31๐) After computing the residual, we transfer the latter to a coarse grid using the restriction operator ๐ ๐,2โ๐ข = โ โ2โ ๐ ๐,โ๐ข , ๐ ๐,2โ๐ค = โ โ2โ ๐ ๐,โ๐ค , ๐ ๐,2โ๐ฆ = โ โ2โ ๐ ๐,โ๐ฆ , ๐ ๐,2โ๐ง = โ โ2โ ๐ ๐,โ๐ง . (32) The restriction operator can be calculated from the prolongation matrix as follows: ๐ซ = 12 [ 2 1 1 2 โฑ โฑ 1 1 2 1 1 2 ] , โ โ2โ = 12 (๐ซ ) ๐ . (33) However, in our case, to avoid matrix-vector multiplication and properly include the Neumann boundary effect on the error and the residual for the ๐ asset options, we employ the method described in the work of Briggs et. al (2000) and obtain ๐ ๐,2โ๐ข directly from ๐ ๐,โ๐ข as follows: (๐ ๐,2โ๐ข ) ๐๐+1 = (๐ ๐,โ๐ข ) + 2(๐ ๐,โ๐ข ) + (๐ ๐,โ๐ข ) โฎ (๐ ๐,2โ๐ง ) ๐๐+1 = (๐ ๐,โ๐ง ) + 2(๐ ๐,โ๐ง ) + (๐ ๐,โ๐ง ) (๐ ๐,2โ๐ข ) = 2(๐ ๐,โ๐ข ) + (๐ ๐,โ๐ข ) (๐ ๐,2โ๐ค ) = (๐ ๐,โ๐ค ) , (๐ ๐,2โ๐ฆ ) = (๐ ๐,โ๐ฆ ) , (๐ ๐,2โ๐ ) = (๐ ๐,โ๐ ) ; (35๐) (๐ ๐,2โ๐ข ) ๐/2๐+1 = (๐ ๐,โ๐ข ) ๐๐+1 , (๐ ๐,2โ๐ค ) ๐/2๐+1 = (๐ ๐,โ๐ค ) ๐๐+1 ; (36๐) (๐ ๐,2โ๐ฆ ) ๐/2๐+1 = (๐ ๐,โ๐ฆ ) ๐๐+1 , (๐ ๐,2โ๐ ) ๐/2๐+1 = (๐ ๐,โ๐ ) ๐๐+1 . (36๐)
Next, we obtain the exact solution of the defect equations 2 ๐ด ๐,2โ ๐ ๐,2โ๐ข = ๐ ๐,2โ๐ข , ๐ต ๐,2โ ๐ ๐,2โ๐ค = ๐ ๐,2โ๐ค , ๐ต ๐,2โ ๐ ๐,2โ๐ฆ = ๐ ๐,2โ๐ฆ , ๐ต ๐,2โ ๐ ๐,2โ๐ง = ๐ ๐,2โ๐ง , (37) where ๐ is the error term on the coarse grid. We then transfer the error term to the fine grid using prolongation (interpolation) operators ๐ ๐,โ๐ข = ๐ซ ๐ ๐,2โ๐ข , ๐ ๐,โ๐ค = ๐ซ ๐ ๐,2โ๐ค , ๐ ๐,โ๐ฆ = ๐ซ ๐ ๐,2โ๐ฆ , ๐ ๐,โ๐ง = ๐ซ ๐ ๐,2โ๐ง . (38) Finally, we correct our approximation of the asset, delta, gamma, and speed options ๐ โ = ๐ โ + ๐ ๐,โ๐ข , ๐ โ = ๐ โ + ๐ ๐,โ๐ค , ๐ โ = ๐ โ + ๐ ๐,โ๐ฆ , ๐ โ = ๐ โ + ๐ ๐,โ๐ง , (39) and relax (30) ๐ times using a smoothing method with our new approximation as the initial guess. Smoothing Methods
Good smoothers are components of effective multigrid methods. This is because few iterations with those relaxation methods remove the high-frequency error. In this work, we use Gauss-Seidel iteration as the smoother and rearrange the compact scheme discretization of our system of partial differential equations in a pointwise manner as follows: (๐ข ๐๐+1 ) = (๐ ๐๐ข ๐ ) โ ๐ ๐1 (๐ข ๐๐ ) ๐ ๐1 , (40๐) (๐ข ๐๐+1 ) ๐๐+1 = (๐ ๐๐ข ๐ ) ๐๐+1/2 โ ๐ ๐1 (๐ข ๐๐+1 ) ๐โ1๐+1 โ ๐ ๐1 (๐ข ๐๐ ) ๐+1๐+1 ๐ ๐1 , (40๐) (๐ค ๐๐+1 ) ๐๐+1 = (๐ ๐๐ค ๐ ) ๐๐+1/2 โ ๐ ๐1 (๐ค ๐๐+1 ) ๐โ1๐+1 โ ๐ ๐1 (๐ค ๐๐ ) ๐+1๐+1 ๐ ๐1 , (41) (๐ฆ ๐๐+1 ) ๐๐+1 = (๐ ๐๐ฆ ๐ ) ๐๐+1/2 โ ๐ ๐1 (๐ฆ ๐๐+1 ) ๐โ1๐+1 โ ๐ ๐1 (๐ฆ ๐๐ ) ๐+1๐+1 ๐ ๐1 , (42) (๐ง ๐๐+1 ) ๐๐+1 = (๐ ๐๐ง ๐ ) ๐๐+1/2 โ ๐ ๐1 (๐ง ๐๐+1 ) ๐โ1๐+1 โ ๐ ๐1 (๐ง ๐๐ ) ๐+1๐+1 ๐ ๐1 , ๐ = 1,2, โฏ , ๐ โ 1. (43) In our previous work, we presented a Gauss-Seidel algorithm for solving our system of equation, hence, we skip it in this section.
Multigrid Algorithms
V-cycle multigrid method : This is an extension of two grids method as described in (30)-(39) by continuous restriction of the residual and relaxation of the defect equation into a pre-defined coarsest grid in a recursive manner. The error is then solved in an exact form in the coarsest grid after which it is 3 interpolated and smoothed to the original fine grid where a better approximation is obtained. An algorithm for obtaining the numerical solutions of the optimal exercise boundary, asset option, and the option Greeks in each regime using the V-cycle multigrid method is described below.
Algorithm 1.
An algorithm based on the V-cycle multigrid Initialize ๐ ๐(๐)๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ฏ ๐ ) โ๐ , (๐ฑ ๐ ) โ๐ , and (๐ช ๐ ) โ๐ for ๐ = 1,2, โฆ , ๐ผ . ๐๐จ๐ซ ๐ง = ๐ ๐ญ๐จ ๐ Compute (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , and (๐ ๐ ) โ๐ for ๐ = 1,2, โฆ , ๐ผ and ๐ โ ๐ based on (29) ๐๐๐ก ๐ ๐(๐)๐+1(๐๐ก=0) = ๐ ๐(๐)๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) =(๐ ๐ ) โ๐ , and (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ . ๐ฐ๐ก๐ข๐ฅ๐ ๐ญ๐ซ๐ฎ๐ Compute (๐ ๐ ) โ๐+1(๐๐ก) , โฆ, and (๐ ๐ ) โ๐+1(๐๐ก) for ๐ = 1,2, โฆ , ๐ผ and ๐ โ ๐ based on (29) Relax ๐ด ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1/3) = ๐ ๐,โ๐ข , โฆ, and ๐ต ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1/3) = ๐ ๐,โ๐ง , and ๐ ๐(๐)๐+1(๐๐ก+1/3) ๐ times using either on (40)-(43) with the initial guess (๐ ๐ ) โ๐+1(๐๐ก) , โฆ, and (๐ ๐ ) โ๐+1(๐๐ก) , and ๐ ๐(๐)๐+1(๐๐ก) Compute ๐ ๐,โ๐ข , โฆ, and ๐ ๐,โ๐ง based on (31) If max |๐ ๐(๐)๐+1(It+1/3) โ ๐ ๐(๐)๐+1(It+1/3) | < ๐ and max |๐ ๐,โ๐ข | < ๐ set (๐ ๐ ) โ๐ = (๐ ๐ ) โ๐+1 , โฆ, and (๐ ๐ ) โ๐ = (๐ ๐ ) โ๐+1 , and ๐ ๐(๐)๐ = ๐ ๐(๐)๐+1 . Compute ( ๐ฏ ๐ ) โ๐+1 ,( ๐ฑ ๐ ) โ๐+1 , and ( ๐ช ๐ ) โ๐+1 . break else , go to 10 If โ is the coarsest grid, go to 17, else Compute ๐ ๐,2โ๐ข , โฆ , ๐ ๐,2โ๐ง based on (32) or (34). Let ๐ โ be the step size of the coarsest grid. For ๐ = 1,2, โฆ , ๐ โ 1, t hen compute recursively: ๏ท Relax ๐ด ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ข = ๐ ๐,2 ๐ โ๐ข , โฆ, and ๐ต ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ง = ๐ ๐,2 ๐ โ๐ง with an initial guess of ๐ ๐,2 ๐ โ๐ข =โฏ = ๐ ๐,2 ๐ โ๐ง = 0 ๏ท Compute ๐ ๐,2 ๐+1 โ๐ข , โฆ , ๐ ๐,2 ๐+1 โ๐ง Solve ๐ด ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ข = ๐ ๐,2 ๐ โ๐ข , โฆ, and ๐ต ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ง = ๐ ๐,2 ๐ โ๐ง based on Compute ๐ฬ ๐,2 ๐โ1 โ๐ข ,..., ๐ฬ ๐,2 ๐โ1 โ๐ง from ๐ ๐,2 ๐ โ๐ข , โฆ , ๐ ๐,2 ๐ โ๐ง based on (38). For ๐ = ๐ โ 1, ๐ โ 2, โฆ ,2, c ompute recursively: ๏ท Calculate ๐ ๐,2 ๐ โ๐ข = ๐ ๐,2 ๐ โ๐ข + ๐ฬ ๐,2 ๐ โ๐ข , โฆ, and ๐ ๐,2 ๐ โ๐ง = ๐ ๐,2 ๐ โ๐ง + ๐ฬ ๐,2 ๐ โ๐ง ๏ท R elax ๐ด ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ข = ๐ ๐,2 ๐ โ๐ข , โฆ, and ๐ต ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ง = ๐ ๐,2 ๐ โ๐ง ๏ท Compute ๐ ๐,2 ๐โ1 โ๐ข = โฏ = ๐ ๐,2 ๐โ1 โ๐ง . Compute (๐ ๐ ) โ๐+1(๐๐ก+2/3) = (๐ ๐ ) โ๐+1(๐๐ก+1/3) + ๐ ๐,โ๐ข , ..., and (๐ ๐ ) โ๐+1(๐๐ก+2/3) = (๐ ๐ ) โ๐+1(๐๐ก+1/3) +๐ ๐,โ๐ง based on (39) Compute (๐ ๐ ) โ๐+1 , (๐ ๐ ) โ๐+1 , (๐ ๐ ) โ๐+1 , and (๐ ๐ ) โ๐+1 for ๐ = 1,2, โฆ , ๐ผ, and ๐ โ ๐ based on (29). Relax ๐ด ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1) = ๐ ๐,โ๐ข , โฆ, and ๐ต ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1) = ๐ ๐,โ๐ง and ๐ ๐(๐)๐+1(๐๐ก+1) ๐ times based on (40)-(43) with the initial guess (๐ ๐ ) โ๐+1(๐๐ก+2/3) , โฆ , (๐ ๐ ) โ๐+1(๐๐ก+2/3) , and ๐ ๐(๐)๐+1(๐๐ก+2/3) endif endwhile 21. endfor Modified cycle multigrid method (M-cycle): Here, we implement a modified multigrid strategy presented in the work of Hafner and Konke (2006). They mentioned that it is necessary for balancing of computational effort in each mesh. In this method, rather than recursive restriction and smoothing across each coarse grid, it is done once on a coarse grid and interpolated back to the fine grid in an increasing fashion starting from the coarsest grid. We observe during the numerical experiment that interpolating with cubic or high order interpolation speeds up convergence with this method. Another distinctness of this approach is that rather than solving the defect equation in exact form, a smoothing process is carried out with many iteration steps. The number of inner iteration steps in each grid which is given as follow: ๐ ๐ = ๐(๐ โ ๐) , (44) where ๐ = 0,1, โฆ , ๐ and ๐ is a constant factor . ๐ = 0 and ๐ = ๐ represent the coarsest and the finest grid, respectively . An algorithm for obtaining the numerical solutions of the optimal exercise boundary, asset option, and the option Greeks in each regime using the M-cycle multigrid method is described below.
Algorithm 2.
An algorithm based on the Modified-cycle multigrid (M-cycle) Initialize ๐ ๐(๐)๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ฏ ๐ ) โ๐ , (๐ฑ ๐ ) โ๐ , and (๐ช ๐ ) โ๐ for ๐ = 1,2, โฆ , ๐ผ . ๐๐จ๐ซ ๐ง = ๐ ๐ญ๐จ ๐ Compute (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐ , and (๐ ๐ ) โ๐ for ๐ = 1,2, โฆ , ๐ผ and ๐ โ ๐ based on (29) ๐๐๐ก ๐ ๐(๐)๐+1(๐๐ก=0) = ๐ ๐(๐)๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ , (๐ ๐ ) โ๐+1(๐๐ก=0) =(๐ ๐ ) โ๐ , and (๐ ๐ ) โ๐+1(๐๐ก=0) = (๐ ๐ ) โ๐ . ๐ฐ๐ก๐ข๐ฅ๐ ๐ญ๐ซ๐ฎ๐ Compute (๐ ๐ ) โ๐+1(๐๐ก) , โฆ, and (๐ ๐ ) โ๐+1(๐๐ก) for ๐ = 1,2, โฆ , ๐ผ and ๐ โ ๐ based on (29) Relax ๐ด ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1/๐) = ๐ ๐,โ๐ข , โฆ, and ๐ต ๐,โ (๐ ๐ ) โ๐+1(๐๐ก+1/๐) = ๐ ๐,โ๐ง , and ๐ ๐(๐)๐+1(๐๐ก+1/๐) ๐ times using either on (40)-(43) with the initial guess (๐ ๐ ) โ๐+1(๐๐ก) , โฆ, and (๐ ๐ ) โ๐+1(๐๐ก) , and ๐ ๐(๐)๐+1(๐๐ก) Compute ๐ ๐,โ๐ข , โฆ, and ๐ ๐,โ๐ง based on (31). If max |๐ ๐(๐)๐+1(It+1/๐) โ ๐ ๐(๐)๐+1(It+1/๐) | < ๐ and max |๐ ๐,โ๐ข | < ๐ Set (๐ ๐ ) โ๐ = (๐ ๐ ) โ๐+1 , โฆ, and (๐ ๐ ) โ๐ = (๐ ๐ ) โ๐+1 , and ๐ ๐(๐)๐ = ๐ ๐(๐)๐+1 . Compute ( ๐ฏ ๐ ) โ๐+1 ,( ๐ฑ ๐ ) โ๐+1 , and ( ๐ช ๐ ) โ๐+1 . break else , go to 10 For ๐ = ๐ โ 1, ๐ โ 2, โฆ ,1, compute recursively:
Compute ๐ ๐,2 ๐ โ๐ข , โฆ , ๐ ๐,2 ๐ โ๐ง Relax ๐ด ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ข = ๐ ๐,2 ๐ โ๐ข , โฆ, and ๐ต ๐,2 ๐ โ ๐ ๐,2 ๐ โ๐ง = ๐ ๐,2 ๐ โ๐ง with an initial guess of ๐ ๐,2 ๐ โ๐ข = โฏ =๐ ๐,2 ๐ โ๐ง = 0 Compute ๐ ๐,โ๐ข ,..., ๐ ๐,โ๐ง from ๐ ๐,2 ๐ โ๐ข , โฆ , ๐ ๐,2 ๐ โ๐ง with cubic or high order interpolation. Compute (๐ ๐ ) โ๐+1[๐๐ก+(๐โ๐+1/๐)] = (๐ ๐ ) โ๐+1[๐๐ก+(๐โ๐/๐)] + ๐ ๐,โ๐ข , โฆ, a nd ( ๐ ๐ ) โ๐+1 [ ๐๐ก+(๐โ๐/๐)+1 ] = ( ๐ ๐ ) โ๐+1 [ ๐๐ก+(๐โ๐/๐) ] + ๐ ๐,โ๐ง based on (39) Compute (๐ ๐ ) โ๐+1 , (๐ ๐ ) โ๐+1 , (๐ ๐ ) โ๐+1 , and (๐ ๐ ) โ๐+1 for ๐ = 1,2, โฆ , ๐ผ and ๐ โ ๐ based on (29) Relax ๐ด ๐,โ (๐ ๐ ) โ๐+1[๐๐ก+(๐โ๐+1/๐)] = ๐ ๐,โ๐ข , โฆ, and ๐ต ๐,โ ( ๐ ๐ ) โ๐+1 [ ๐๐ก+ ( ๐โ๐+1/๐ )] = ๐ ๐,โ๐ง and ๐ ๐(๐)๐+1[๐๐ก+(๐โ๐+1/๐)] ๐ times based on (40)-(43) endif 20. endwhile 21. endfor Numerical Experiment
To check the performance of our proposed algorithms, we test their accuracies and present numerical examples. In the numerical example, we only consider American put options pricing problem with two regimes. The numerical experiment was carried out on the mesh with a uniform grid size. The numerical code was written with MATLAB 2019a on Intel Core i5-3317U CPU 1.70GHz 64-bit ASUS Laptop.
Two Regimes Example
In this example, we consider a two-regimes problem with the following data:
๐พ = 9, ๐ = 1, ๐ = [โ6 69 โ 9] , ๐ = [0.100.05] , ๐ = [0.800.30] , ๐ = 10 โ8 . (45) In our computation, we chose the interval ๐ โค 3. Our time step ๐ was chosen such that ๐ = โ . We label results based on multigrid methods as follows: ๏ท ๐(๐บ๐: ๐บ๐) โ
V-cycle multigrid with Gauss-Seidel smoother, ๏ท ๐(๐บ๐: โฏ : ๐บ๐) โ
M-cycle multigrid with Gauss-Seidel smoother. We compared the multigrid methods with MTree (Liu, 2010), IMS1, IMS2 (Khaliq and Liu, 2009), and MOL (Chiarella et al., 2016). which we label as FF-CS1 and listed the results in Tables 1 and 2. The plots of the option price, option Greeks, and optimal exercise boundary for each regime were displayed in Figs. 1 and 2. From Tables 1-2, based on the comparison between the multigrid methods and other existing methods, we can observe that the results from the multigrid methods are very close to MTree (Liu, 2010) and MOL (Chiarella et al., 2016). It is important to point out that Chiarella et al. (2016) mentioned that MTree data was used as the benchmark result in the work of Khaliq and Liu (2009). 6
Fig. 1.
Optimal exercise boundaries, asset options, and option Greeks with V-cycle multigrid method in Algorithm 1 ( ๐ = ๐, โ = 0.00625 ) . Fig. 2.
Optimal exercise boundaries, asset options, and option Greeks with M-cycle multigrid method in Algorithm 2 ( ๐ = ๐, โ = 0.00625 ) . Table 1.
Comparison of the American put option price in regime 1. S MTree IMS1 IMS2 MOL
๐(1,1)
๐(2,2,2,2) โ = 0.0125 0.00625 0.0125 0.00625
Table 2.
Comparison of the American put option price in regime 2. S MTree IMS1 IMS2 MOL
๐(1,1)
๐(2,2,2,2) โ = 0.0125 0.00625 0.0125 0.00625 ๐ norm per each iteration at the final time level and displayed the plot in Fig. 3. Next, we computed the convergent factor of the multigrid methods. For the numerical convergent factor at each time level, we estimated it as follows (Gander and Neumuller, 2014): ๐ = max โ๐ ๐+1 โ โ๐ ๐ โ , (46) where ๐ ๐ is the residual at ๐ iteration for a given time level. The result of the numerical convergent factor was displayed in Fig. 4. In Table 3, we confirmed the h-independency of the M-cycle multigrid using the measured convergent factor. We further verified the global maximum iteration and the average CPU time of our previous method with the V-cycle and M-cycle methods in Algorithms 1 and 2 and displayed the results in Table 4. 9 Fig. 3. ๐ norm of the normalized residue and log scale of the absolute error in Regime 1 ( โ =0.0125; ๐ = 2; ๐(2,2,2,2); ๐(1,1)) . Fig. 4.
The convergent factor of the multigrid methods ( โ = 0.0125; ๐ = 2; ๐(2,2,2,2); ๐(1,1)) . Table 3.
Confirming the โ independency of the M-cycle using the measured convergent factor when ๐ =๐ with different grid sizes and fixed time step ๐ = 0.00015625 . โ ๐(๐บ๐: ๐บ๐: ๐บ๐: ๐บ๐)
Table 4.
Comparison of the global maximum iterations and average CPU time(s) where โ = 0.0125 , inner smoothing maximum number of iterations = 2 and ๐ = 2 . Method
๐(1,1) ๐(10,10)
๐(1,1,1,1) ๐(2,2,2,2) ๐(1,1,1,1,1) ๐(2,2,2,2,2)
Previous Method CPU time(s)
Moreover, it can be seen from Table 4 that the number of iterations required to achieve numerical convergence with a tolerance of ๐ = 10 โ8 is smaller in the multigrid methods when compared with our previous method which is based on classical Gauss-Seidel iteration. From Fig. 4, we can further observe that the convergence factor of the V-cycle multigrid is closer to one than the M-cycle multigrid. Moreover, M-cycle multigrid requires less CPU time than the V-cycle multigrid. This implies that the M-cycle multigrid performs better than the V-cycle multigrid in terms of numerical convergence. One can further observe that the CPU time for the multigrid methods is twice smaller than our previous method; a clear indication of the good performance of the multigrid methods over our previous method. Moreover, the M-cycle method in Algorithm 2 has less CPU time when compared with the V-cycle method in Algorithm 1. Comparing the M-cycle method based on the number of grid sequence from Table 4, the M-cycle multigrid with five grid sequence has less CPU time than that of four grid sequence, a performance that could be very useful when modeling a large-scale problem or system with large state space that requires small step size to approximate the numerical solution of a function. Finally, we can observe that the measured convergent factor in Table 3 is almost independent of the varying step sizes for a fixed time step; an indication of the h-independent convergence of the multigrid. Conclusion
We have presented two multigrid strategies in Algorithms 1 and 2 for solving American put options with regime-switching and investigated their numerical performances with our previous method (that is based on the classical Gauss-Seidel iterative method) and other existing methods. The numerical investigation is carried out with a two-regimes example. From the numerical results, the capability of the multigrid methods can easily be observed. The global maximum iteration required to achieve numerical convergence with ๐ = 10 โ8 is much smaller with multigrid methods when compared with our previous method that is based on the Gauss-Seidel iteration. In terms of error per iteration, it can be observed that the error reduces faster in multigrid methods when compared with our previous method. Moreover, the error reduces much faster in the M-cycle method when compared with the V-cycle method. The multigrid methods require less iteration and CPU time when compared with our previous method. One can also observe that the M-cycle method requires the least CPU time and have a convergent factor less than the V-cycle method which is a clear indication of its superior performance over the V-cycle method. Furthermore, M-cycle multigrid with five grid sequence has less CPU time than that of four grid sequence. In general, our result shows that the multigrid methods could be very useful for American options valuation with switching regimes, especially when solving problems with a large finite state space. 1 References Adam, Y. (1975) A Hermitian finite difference method for the solution of parabolic equations. Computers and Mathematics with Applications, 1:393-406. 2.
Adam, Y. (1977) Highly accurate compact implicit methods and boundary conditions. Journal of Computational Physics, 24:10-22. 3.
Bachvalov. N. S. (1966) On the convergence of a relaxation method with natural constraints on the elliptic Operator. USSR Comput. Math. and Math. Phys. 6:101-135. 4.
Brandt, A. (1977) Multi-level adaptive solutions to boundary-value problems. Numer. Math. 31:333-390. 5.
Brigg, W. L., Henson, V. E., McCormick, S. F. (2000) A Multigrid tutorial. SIAM. 6.
Burden, R. L., Faires, D. J., and Burden, A. M. (2016) Numerical analysis. Cengage Learning, Boston, MA. 7.
Chiarella, C., Nikitopoulos Sklibosios, C., Schlogl, E., and Yang, H. (2016) Pricing American options under regime switching using method of lines. SSRN Electronic Journal. 8.
Ciesielski, P. E., Fulton, S. R., and Schubert, W. H. (1985) Multigrid solution of an elliptic boundary value problem from tropical cyclone theory. Monthly Weather Review, 116:797-807. 9.
Clarke, N., and Parrott. K. (1999) Multigrid for American option pricing with stochastic volatility. Applied Mathematical Finance
Dremkova, E., and Ehrhardt, M. (2011) A high-order compact method for non-linear Black-Scholes option pricing equation of American options. International Journal of Computational Mathematics, 88:2782-2797. 11.
Egecioglu, O., Gallopoulos, E., and Koc, C. K. (1989) Fast computation of divided differences and parallel Hermite interpolation. Journal of Complexity, 5:417-437. 12.
Egorova, V. N., Company, R., and Jรณdar, L. (2016) A new efficient numerical method for solving American option under regime switching model. Computers and Mathematics with Applications, 71:224โ237. 13.
Fedorenko, R. P. (1961) A relaxation method for solving elliptic difference equations. USSR Comput. Math. and Math. Phys. 1:1092-1096. 14.
Gander, M., and Neumuller, M. (2014) Analysis of a time multigrid algorithm for DG-discretization in time. arXiv Preprint: arXiv:1409.5254v2 [math.NA] Ghaffar, F., Badshah, N., Islam, I., and Khan, M. A. (2016) Multigrid method based on transformation-free high-order scheme for solving 2D Helmholtz equation on nonuniform grids. Advances in Difference Equations. 16.
Guo, Z., Mi, J., and Grant, P.S. (2012) An implicit parallel multigrid computing scheme to solve coupled thermal-solute phase-field equations for dendrite evolution. Journal of Computational Physics, 231 (1781-1796). 17.
Hafner, S., and Konke, C. (2006) Multigrid preconditioned conjugate gradient method in the mechanical analysis of heterogeneous solids. 17th International Conference on the Application of Computer Science and Mathematics in Architecture and Civil Engineering. 18.
Han, Y., and Kim, G. (2016) Efficient lattice method for valuing of options with barrier in a regime switching model. Discrete Dynamics in Nature and Society. 19.
Hirsh, R. S. (1975) Higher order accurate difference solutions of fluid mechanics by a compact differencing technique. Journal of Computational Physics, 19:90-109. 20.
Horton, G., and VandeWalle, S. (1995) A space-time multigrid method for parabolic partial differential equations. SIAM Journal of Scientific Computing, 16:848โ864. 21.
Jeong, D., Kim, J., and Wee, I. S. (2009) An accurate and efficient numerical method for Black-Scholes equations. Commun. Korean Math. Soc. 24:617โ628. 22.
Kangro, R., and Nicolaides, R. (2000) Far field boundary conditions for Black--Scholes equations. SIAM Journal on Numerical Analysis, 38:1357-1368. 23.
Khaliq, A. Q. M., and Liu, R. H. (2009) New numerical scheme for pricing American option with regime-switching. International Journal of Theoretical and Applied Finance, 12:319-340. 24.
Khelifi, S. C., Mechitoua, N., Hulsemann, F., and Magoules, F. (2014) A hybrid multigrid method for convectionโdiffusion problems. Journal of Computational and Applied Mathematics, 259:711โ719. 25.
Lee, H. N., and Meyer, R. E. (1979) On time dependent multi-grid numerical technique. Comp. & Maths. with Appl., 6:61-65. 26.
Liao, W., and Khaliq, A. Q. M. (2009) High-order compact scheme for solving nonlinear Black-Scholes equation with transaction cost. International Journal of Computer Mathematics, 86:1009-1023. 27.
Liu, R. H. (2010). Regime-switching recombining tree for option pricing. International Journal of Theoretical and Applied Finance, 13:479-499. 28.
Meyer, G. H., and van der Hoek, J. (1997) The evaluation of American options with the method of lines. 3
Moghaderi, H., and
Dehghan, M. (2014) A multigrid compact finite difference method for solving the one-dimensional nonlinear sine-Gordon equation. Math. Meth. Appl. Sci., 38:3901โ3922. 30.
Nielsen, B. F., Skavhaug O., and Tveito A. (2002) A penalty and front-fixing methods for the numerical solution of American option problems. Journal of Computational Finance, 5:69-97. 31.
Nwankwo, C., Dai, W., and Liu, R. (2019) Compact finite difference scheme with Hermite interpolation for pricing American put options based on regime switching model. arXiv Preprint: arXiv:1908.04900v6 [q-fin.CP]. 32.
Pardhanani, A. L., Spotz, W. F., and Carey G. F. (1997) A stable multigrid strategy for convection-diffusion using high order compact discretization. Electronic Transaction on Numerical Analysis, 6:211-223. 33.
Powell, M. J. D. (1981) Approximation theory and methods. Cambridge University Press, Cambridge. 34.
Rosseel, E., Boonen, T., and VandeWalle, S. (2007) Algebraic multigrid for stationary and timeโdependent partial differential equations with stochastic coefficients. Numer. Linear Algebra Appl., 15:141โ163. 35.
Saad, Y. (2003) Iterative methods for sparse linear systems. SIAM. 36.
Shang, Q., and Bryne, B. (2019) An efficient lattice search algorithm for the optimal exercise boundary in American options. SSRN Electronic Journal. 37.
Shieh, A. S. (1985) On the solution of coupled system of PDE by a multigrid method. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 4:527-530. 38.
Trottenberg, U., Oosterlee, C., and Schuller, A. (2001) Multigrid. Academic Press, New York. 39.
Urschel, J. C. (2013) A space-time multigrid method for the numerical valuation of barrier options. Communications in Mathematical Finance, 2:1-20. 40.
Wang, Y., and Ge, Y. (2018) High-order compact difference scheme and multigrid method for solving the 2D elliptic problems. Mathematical Problem in Engineering. 41.
Yavneh, I. (2006). Why multigrid methods are so efficient. Comput. Sci. Engrg. 8:12โ22. 42.
Yan, Y., Dai, W., Wu, L., and Zhai, S. (2019) Accurate gradient preserved method for solving heat conduction equations in double layers. Applied Mathematics and Computation, 354:58-85. 43.