Compact Finite Difference Scheme with Hermite Interpolation for Pricing American Put Options Based on Regime Switching Model
11 Compact Finite Difference Scheme with Hermite Interpolation for Pricing American Put Options Based on Regime Switching Model
Chinonso I. Nwankwo a , Weizhong Dai a, *, Ruihua Liu b a Department of Mathematics and Statistics, Louisiana Tech University, Ruston LA 71272, USA b Department of Mathematics, University of Dayton, 300 College Park, Dayton, OH 45469, USA * Corresponding author, [email protected]
Abstract
We consider a system of coupled free boundary problems for pricing American put options with regime-switching. To solve this system, we first employ the logarithmic transformation to map the free boundary for each regime to multi-fixed intervals and then eliminate the first-order derivative in the transformed model by taking derivatives to obtain a system of partial differential equations which we call the asset-delta-gamma-speed equations. As such, the fourth-order compact finite difference scheme can be used for solving this system. The influence of other asset, delta, gamma, and speed options in the present regime is estimated based on Hermite interpolations. Finally, the numerical method is tested with several examples. Our results show that the scheme provides an accurate solution that is fast in computation as compared with other existing numerical methods.
Keywords:
American put options with regime switching, logarithmic transformation, optimal exercise boundary, compact finite difference method, Hermite interpolation Introduction
The well-known Black-Scholes model has been used over decades in options valuation. This model constructs a delta hedging portfolio with an assumption of the frictionless market, no-arbitrage, and constant risk-free interest and volatility (Ugur, 2006). To remove this ideal assumption and reproduce the actual market price, risk, behavior, and dynamics, researchers have proposed several improvements by including stochastic volatility (Chockalingam and Muthuraman, 2011; Düring and Fournié, 2012; Garnier and Sølna, 2017; Huang et al., 2011; Hull and White, 1987; Ikonen and Toivanen, 2007; Zhylyevskyy, 2009), jump-diffusion (Bingham, 2006; Chen et al., 2019; Cont and Tankov, 2004; Guoqing and Hanson, 2006; Kou, 2002), and regime-switching (Company et al., 2016a; Egorova et al., 2016; Huang et al., 2011; Khaliq and Liu, 2009; Mamon and Rodrigo, 2005) in the pricing models. The regime-switching model for American option valuation, first introduced by Hamilton (1989), has gained broader interest after the seminal work of Buffington and Elliot (2002). It defines a finite number of market states known as regimes. Each regime has its own set of market variables, and the market randomly switches among different regimes (Chiarella et al., 2016). The model for option valuation with regime-switching involves a system of partial differential equations with free boundaries for which the analytical solution is very difficult to obtain in general. Thus, some works in the literature have proposed numerical techniques for solving the option pricing equation with regime-switching. Among them, the commonly known numerical methods are the penalty method (Khaliq and Liu, 2009; Nielsen et al., 2002; Zhang et al., 2013), the method of line (MOL) (Chiarella et al., 2016; Meyer and van der Hoek, 1997), the lattice method (Han and Kim, 2016; Shang and Bryne, 2019), the fast Fourier transform (Boyarchenko and Levendorskii, 2008; Liu et al., 2006), and the front-fixing techniques (Egorova et al., 2016). The lattice-based method is more common among practitioners. However, tracking the optimal exercise boundary can be a challenge (Shang and Bryne, 2019). Fast Fourier transform method is efficient in solving the European options (Chiarella et al., 2016).
The penalty method removes the free boundary by introducing a penalty term (Khaliq and Liu, 2009). The MOL method calculates the asset and delta options and the optimal exercise boundary simultaneously during computation. Meyer and van der Hoek (1997) pointed out that there are still some complications with the MOL method due to the singularity of the solution and infinite interval. The front-fixing technique (Blackwell and Hogan, 1994; Company et al., 2016b; Company et al., 2016c; Landau, 1950; Mitchell and Vynnycky, 2009; Mitchell and Vynnycky) was first applied by Egorova et al. (2016) to the regime-switching model. To the best of our knowledge, the above methods provide up to second-order accurate solutions. The motivation of this research is to propose a higher accurate front-fixing numerical method for solving the regime-switching pricing model. To this end, we first use a logarithmic transformation to map the free boundary for each regime to multi-fixed intervals and then eliminate the first-order derivative in the transformed model by taking derivatives to obtain a system of partial differential equations which we call the asset-delta-gamma-speed equations. As such, the fourth-order compact finite difference can be used for solving this system. The influence of other asset, delta, gamma, speed options in the present regime is estimated based on Hermite interpolations. Finally, the numerical scheme is solved using either the Gauss-Seidel or Newton iterative method, which predicts the optimal exercise boundary, option value, and option Greeks in each regime. The rest of the paper is organized as follows. In section 2, we consider a regime-switching model and its transformations. We transform the model to obtain coupled partial differential equations for option values, delta, gamma, and speed options in each regime. In section 3, we develop an accurate numerical method and its algorithms for solving these equations and obtaining the option values, optimal exercise boundary, and the Greeks in each regime. In section 4, we test our algorithms using examples with two, four, eight, and sixteen regimes. We conclude the paper in section 5. Regime Switching Model and its Transformations
Regime Switching Model
Let us consider a continuous-time Markov chain whose states are labeled as 𝑚 = 1,2,∙∙∙, 𝐼.
Let
𝑄 =(𝑞 𝑚𝑙 ) 𝐼×𝐼 represent the generator matrix with the entry elements 𝑞 𝑚𝑙 satisfying the condition below (Norris, 1998): 𝑞 𝑚𝑚 = − ∑ 𝑞 𝑚𝑙 , 𝑙≠𝑚 𝑞 𝑚𝑙 ≥ 0, 𝑓𝑜𝑟 𝑙 ≠ 𝑚, 𝑙 = 1,2,∙∙∙, 𝐼. (1) Assuming a risk-neutral measure (Elliot et al., 2007), the underlying asset follows a stochastic process 𝑑𝑆 𝑡 = 𝑆 𝑡 (𝑟 𝛼 𝑡 𝑑𝑡 + 𝜎 𝛼 𝑡 𝑑𝐵 𝑡 ), 0 ≤ 𝑡 < ∞, (2) where 𝑟 𝛼 𝑡 and 𝜎 𝛼 𝑡 are the interest rate and volatility of the asset, respectively, and are dependent on the Markov chain state with 𝑟 𝛼 𝑡 |𝛼 𝑡 = 𝑟 𝑚 , 𝜎 𝛼 𝑡 |𝛼 𝑡 = 𝜎 𝑚 , 𝑚 = 1,2,∙∙∙, 𝐼. ( ) We consider an American put option written on the asset 𝑆 𝑡 with strike price 𝐾 and expiration time 𝑇. Let 𝑉 𝑚 (𝑆, 𝑡) denote the option price and 𝜏 = 𝑇 − 𝑡 . Then, 𝑉 𝑚 (𝑆, 𝜏) satisfies the following parabolic PDEs with free boundaries (Khaliq and Liu, 2009): − 𝜕𝑉 𝑚 (𝑆, 𝜏)𝜕𝜏 + 12 𝜎 𝑆 𝜕 𝑉 𝑚 (𝑆, 𝜏)𝜕𝑆 + 𝑟 𝑚 𝑆 𝜕𝑉 𝑚 (𝑆, 𝜏)𝜕𝑆 − 𝑟 𝑚 𝑉 𝑚 (𝑆, 𝜏) + ∑ 𝑞 𝑚𝑙𝑙≠𝑚 [𝑉 𝑙 (𝑆, 𝜏) − 𝑉 𝑚 (𝑆, 𝜏)] = 0,for 𝑆 > 𝑠 𝑓(𝑚) (𝜏), (4) 𝑉 𝑚 (𝑆, 𝜏) = 𝐾 − 𝑆, for 𝑆 < 𝑠 𝑓(𝑚) (𝜏). (5) Here, the initial and boundary conditions are given as: 𝑉 𝑚 (𝑆, 0) = 𝑚𝑎𝑥( 𝐾 − 𝑆, 0), 𝑠 𝑓(𝑚) (0) = 𝐾; (6) 𝑉 𝑚 (𝑠 𝑓(𝑚) , 𝜏) = 𝐾 − 𝑠 𝑓(𝑚) (𝜏), 𝑉 𝑚 (0, 𝜏) = 𝐾, 𝑉 𝑚 (∞, 𝜏) = 0, 𝜕𝜕𝑆 𝑉 𝑚 (𝑠 𝑓(𝑚) , 𝜏) = −1, (7) where 𝑠 𝑓(𝑚) (𝜏) is the optimal exercise boundary for the 𝑚 𝑡ℎ regime. Logarithmic Transformation
To fix the free boundary challenge, we employ a transformation (Egorova et al., 2016; Wu and Kwok, 1997) on multi-variable domains as 𝑥 𝑚 = ln 𝑆𝑠 𝑓(𝑚) (𝜏) = ln 𝑆 − ln 𝑠 𝑓(𝑚) (𝜏), 𝑚 = 1,2,∙∙∙, 𝐼, (8) where the variable 𝑥 𝑚 exists in a positive domain 𝑥 𝑚 ∈ (0, ∞) if 𝑆 > 𝑠 𝑓(𝑚) . The transformed 𝑚 option value functions 𝑈 𝑚 (𝑥 𝑚 , 𝜏) are related to the original 𝑚 option value functions 𝑉 𝑚 (𝑆, 𝜏) by the dimensionless transformation 𝑈 𝑚 (𝑥 𝑚 , 𝜏) = 𝑉 𝑚 (𝑆, 𝜏), 𝑚 = 1,2,∙∙∙, 𝐼. (9𝑎) Applying this transformation, we obtain the following relations: 𝜕𝑥 𝑚 𝜕𝑆 = 1𝑆 , 𝜕𝑥 𝑚 𝜕𝜏 = − 𝑠 ′𝑓(𝑚) (𝜏)𝑠 𝑓(𝑚) (𝜏) , 𝜕𝑉 𝑚 𝜕𝑆 = 1𝑆 𝜕𝑈 𝑚 𝜕𝑥 𝑚 ; (9𝑏) 𝜕 𝑉 𝑚 𝜕𝑆 = 1𝑆 (− 𝜕𝑈 𝑚 𝜕𝑥 𝑚 + 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 ), 𝜕𝑉 𝑚 𝜕𝜏 = 𝜕𝑈 𝑚 𝜕𝜏 − 𝑠 ′𝑓(𝑚) (𝜏)𝑠 𝑓(𝑚) (𝜏) 𝜕𝑈 𝑚 𝜕𝑥 𝑚 . (9𝑐) Because our interest is to also calculate speed, delta decay, and color options, we differentiate further to obtain higher derivatives of the 𝑚th option value function as 𝜕 𝑉 𝑚 𝜕𝑆 = 1𝑆 (2 𝜕𝑈 𝑚 𝜕𝑥 𝑚 − 3 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 + 𝜕 𝑈 𝑚 𝜕𝑥 𝑚3 ), 𝜕 𝑉 𝑚 𝜕𝑆𝜕𝜏 = 1𝑆 ( 𝜕 𝑈 𝑚 𝜕𝑥 𝑚 𝜕𝜏 − 𝑠 ′𝑓(𝑚) (𝜏)𝑠 𝑓(𝑚) (𝜏) 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 ) ; (9𝑑) 𝜕 𝑉 𝑚 𝜕𝑆 𝜕𝜏 = 1𝑆 ( 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 𝜕𝜏 − 𝜕 𝑈 𝑚 𝜕𝑥 𝑚 𝜕𝜏 + 𝑠 ′𝑓(𝑚) (𝜏)𝑠 𝑓(𝑚) (𝜏) 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 − 𝑠 ′𝑓(𝑚) (𝜏)𝑠 𝑓(𝑚) (𝜏) 𝜕 𝑈 𝑚 𝜕𝑥 𝑚3 ). (9𝑒) Let 𝑙 represent the coupled regime(s) in the 𝑚 free boundary PDEs. The former also has a variable 𝑥 𝑙 = ln 𝑆𝑠 𝑓(𝑙) (𝜏) = ln 𝑆 − ln 𝑠 𝑓(𝑙) (𝜏), 𝑙 ≠ 𝑚, 𝑙 = 1,2,∙∙∙, 𝐼. (10) Eliminating 𝑆 in the 𝑙 𝑡ℎ and 𝑚 𝑡ℎ equations, we obtain 𝑥 𝑙 = 𝑥 𝑚 − ln 𝑠 𝑓(𝑙) (𝜏)𝑠 𝑓(𝑚) (𝜏). (11) Substituting (8) into (4) and (5) (i.e.,
𝑆 = 𝑠 𝑓 ( 𝑚 ) ( 𝜏 ) 𝑒 𝑥 𝑚 ) , the model can be changed to 𝜕𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝜏 − 12 𝜎 𝜕 𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚 + 𝑟 𝑚 𝑈 𝑚 (𝑥 𝑚 , 𝜏)− ∑ 𝑞 𝑚𝑙𝑙≠𝑚 [𝑈 𝑙 (𝑥 𝑚 , 𝜏) − 𝑈 𝑚 (𝑥 𝑚 , 𝜏)] = 0, for 𝑥 𝑚 > 0; (12) 𝑈 𝑚 (𝑥 𝑚 , 𝜏) = 𝐾 − 𝑆 = 𝐾 − 𝑠 𝑓(𝑚) (𝜏)𝑒 𝑥 𝑚 , for 𝑥 𝑚 < 0; (13) where the initial condition (5) is changed to 𝑈 𝑚 (𝑥 𝑚 , 0) = max(𝐾 − 𝐾𝑒 𝑥 𝑚 , 0) = 0, 𝑥 𝑚 ≥ 0, 𝑠 𝑓(𝑚) (0) = 𝐾. (14) By letting 𝑥 𝑚 → 0 − , we obtain from (13) that 𝑈 𝑚 (0, 𝜏) = 𝐾 − 𝑠 𝑓(𝑚) (𝜏) . Thus, together with (7), we obtain the boundary condition for (12) as 𝑈 𝑚 (0, 𝜏) = 𝐾 − 𝑠 𝑓(𝑚) (𝜏), 𝑈 𝑚 (∞, 𝜏) = 0. (15) To apply the high order compact finite difference method, we further transform the system in (12)-(15) by eliminating the first-order derivative. To this end, we let 𝑊 𝑚 (𝑥 𝑚 , 𝜏) represent the derivative of the option value in each regime known as the delta option and given as 𝑊 𝑚 (𝑥 𝑚 , 𝜏) = 𝜕𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚 , for 𝑥 𝑚 > 0 and 𝑥 𝑚 < 0. (16) Differentiating (12) and (13) with respect to 𝑥 𝑚 , respectively , we generate a system of partial differential equations in terms of delta option as 𝜕𝑊 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝜏 − 12 𝜎 𝜕 𝑊 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 + 𝑟 𝑚 𝑊 𝑚 (𝑥 𝑚 , 𝜏)− ∑ 𝑞 𝑚𝑙𝑙≠𝑚 (𝑊 𝑙 (𝑥 𝑚 , 𝜏) − 𝑊 𝑚 (𝑥 𝑚 , 𝜏)) = 0, for 𝑥 𝑚 > 0; (17) 𝑊 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) 𝑒 𝑥 𝑚 , for 𝑥 𝑚 < 0; (18) where the initial condition for (17) is obtained based on (14) as : 𝑊 𝑚 (𝑥 𝑚 , 0) = 0, 𝑥 𝑚 ≥ 0. (19) By letting 𝑥 𝑚 → 0 − in (18) together with (7), we obtain the boundary condition for (17) as 𝑊 𝑚 (0, 𝜏) = −𝑠 𝑓(𝑚) , 𝑊 𝑚 (∞, 𝜏) = 0. (20) It is important to point out that for 𝑊 𝑚 (𝑥 𝑚 , 𝜏) at 𝑥 𝑚 = 0 when 𝜏 = 0 , its value obtained based on the initial condition and the boundary condition are different. This happens in many PDE problems. We are mostly concerned with what happens for 𝑊 𝑚 (𝑥 𝑚 , 𝜏) and other functions in 𝑥 𝑚 ≥ 0 when 𝜏 > 0 . Here, we set 𝑊 𝑚 (0, 0) = 0 . We have used the average of the limit from the left and right as the value of 𝑊 𝑚 (0, 0) . We have also used the smoothstep method (Bravo and Mcgraw, 2015) to approximate 𝑊 𝑚 (0, 0). However, when comparing them with our choice of 𝑊 𝑚 (0, 0) = 0, we found no significant difference. Furthermore, we let 𝑌 𝑚 (𝑥 𝑚 , 𝜏) represent the derivative of the delta option in the 𝑚 𝑡ℎ regime known as the gamma option and obtain 𝑌 𝑚 (𝑥 𝑚 , 𝜏) = 𝜕𝑊 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚 = 𝜕 𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 , for 𝑥 𝑚 > 0 and 𝑥 𝑚 < 0. (21) Differentiating (17) and (18) with respect to 𝑥 𝑚 , respectively, we generate a system of gamma option PDEs for each regime of the form 𝜕𝑌 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝜏 − 12 𝜎 𝜕 𝑌(𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑊 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 + 𝑟 𝑚 𝑌 𝑚 (𝑥 𝑚 , 𝜏)− ∑ 𝑞 𝑚𝑙𝑙≠𝑚 (𝑌 𝑙 (𝑥 𝑚 , 𝜏) − 𝑌 𝑚 (𝑥 𝑚 , 𝜏)) = 0, for 𝑥 𝑚 > 0; (22) 𝑌 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) 𝑒 𝑥 𝑚 , for 𝑥 𝑚 < 0; (23) where the initial condition for (22) is obtained based on (19) as 𝑌 𝑚 (𝑥 𝑚 , 0) = 0, 𝑥 𝑚 ≥ 0. (24) By letting 𝑥 𝑚 → 0 − in (23) together with (7), we obtain the boundary condition for (22) 𝑌 𝑚 (0, 𝜏) = −𝑠 𝑓(𝑚) , 𝑌 𝑚 (∞, 𝜏) = 0. (25) Finally, we let 𝑍 𝑚 (𝑥 𝑚 , 𝜏) represent the derivative of the gamma option known as the speed option and obtain 𝑍 𝑚 (𝑥 𝑚 , 𝜏) = 𝜕𝑌 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚 = 𝜕 𝑊 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 = 𝜕 𝑈 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 , for 𝑥 𝑚 > 0 and 𝑥 𝑚 < 0, (26) Differentiating (22), (23) with respect to 𝑥 𝑚 , we generate a system of speed option PDEs as 𝜕𝑍 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝜏 − 12 𝜎 𝜕 𝑍(𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑌 𝑚 (𝑥 𝑚 , 𝜏)𝜕𝑥 𝑚2 + 𝑟 𝑚 𝑍 𝑚 (𝑥 𝑚 , 𝜏)− ∑ 𝑞 𝑚𝑙𝑙≠𝑚 (𝑍 𝑙 (𝑥 𝑚 , 𝜏) − 𝑍 𝑚 (𝑥 𝑚 , 𝜏)) = 0, for 𝑥 𝑚 > 0; (27) 𝑍 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) 𝑒 𝑥 𝑚 , for 𝑥 𝑚 < 0; (28) where the initial and boundary conditions for (27) are given as : 𝑍 𝑚 (𝑥 𝑚 , 0) = 0, 𝑥 𝑚 ≥ 0, 𝑍 𝑚 (0, 𝜏) = −𝑠 𝑓(𝑚) , 𝑍 𝑚 (∞, 𝜏) = 0. (29) Thus, a set of asset-delta-gamma-speed option PDEs in each regime can be written as follows: 𝜕𝑈 𝑚 𝜕𝜏 − 12 𝜎 𝜕 𝑈 𝑚 𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑚 + (𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑈 𝑚 − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑈 𝑙 = 0, (30𝑎) 𝜕𝑊 𝑚 𝜕𝜏 − 12 𝜎 𝜕 𝑊 𝑚 𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑈 𝑚 𝜕𝑥 𝑚2 + (𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑊 𝑚 − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑊 𝑙 = 0, (30𝑏) 𝜕𝑌 𝑚 𝜕𝜏 − 12 𝜎 𝜕 𝑌 𝑚 𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑊 𝑚 𝜕𝑥 𝑚2 + (𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑌 𝑚 − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑌 𝑙 = 0, (30𝑐) 𝜕𝑍 𝑚 𝜕𝜏 − 12 𝜎 𝜕 𝑍 𝑚 𝜕𝑥 𝑚2 − (𝑠 ′𝑓(𝑚) 𝑠 𝑓(𝑚) + 𝑟 𝑚 − 𝜎 𝑌 𝑚 𝜕𝑥 𝑚2 + (𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑍 𝑚 − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑍 𝑙 = 0, (30𝑑) where 𝑚 = 1,2,∙∙∙, 𝐼, 𝑥 𝑚 > 0, and the initial and boundary conditions for 𝑈 𝑚 (𝑥 𝑚 , 𝜏) , 𝑊 𝑚 (𝑥 𝑚 , 𝜏) , 𝑌 𝑚 (𝑥 𝑚 , 𝜏) , and 𝑍 𝑚 (𝑥 𝑚 , 𝜏) are given as : 𝑈 𝑚 (𝑥 𝑚 , 0) = 0, 𝑊 𝑚 (𝑥 𝑚 , 0) = 0, 𝑌 𝑚 (𝑥 𝑚 , 0) = 0; (31𝑎) 𝑍 𝑚 (𝑥 𝑚 , 0) = 0, 𝑠 𝑓(𝑚) (0) = 𝐾, 𝑥 𝑚 ≥ 0; (31𝑏) 𝑈 𝑚 (0, 𝜏) = 𝐾 − 𝑠 𝑓(𝑚) (𝜏), 𝑊 𝑚 (0, 𝜏) = 𝑌 𝑚 (0, 𝜏) = 𝑍 𝑚 (0, 𝜏) = −𝑠 𝑓(𝑚) (𝜏); (31𝑐) 𝑈 𝑚 (∞, 𝜏) = 0; 𝑊 𝑚 (∞, 𝜏) = 0, 𝑌 𝑚 (∞, 𝜏) = 0, 𝑍 𝑚 (∞, 𝜏) = 0. (31𝑑) On the other hand, for 𝑥 𝑚 < 0 , 𝑈 𝑚 (𝑥 𝑚 , 𝜏) = 𝐾 − 𝑠 𝑓(𝑚) (𝜏)𝑒 𝑥 𝑚 , 𝑊 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) (𝜏)𝑒 𝑥 𝑚 ; (32𝑎) 𝑌 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) (𝜏)𝑒 𝑥 𝑚 , 𝑍 𝑚 (𝑥 𝑚 , 𝜏) = −𝑠 𝑓(𝑚) (𝜏)𝑒 𝑥 𝑚 . (32𝑏) Numerical Formulation
To solve the above asset-delta-gamma-speed option PDEs, we first design a uniform grid [0, ∞) × [0 𝑇] for each regime taking into consideration how the 𝑚 𝑡ℎ regime’s interval relates to the 𝑙 𝑡ℎ regime’s interval using the Hermite interpolation technique. The infinite boundary is replaced with the far estimate boundary (Egorova et al., 2016; Kangro and Nicolaides, 2000; Toivanen, 2010), which we denote as (𝑥 𝑚 ) 𝑀 . Representing 𝑖 as the node point in the 𝑚 𝑡ℎ regime’s interval, 𝑗 as the node point in the 𝑙 𝑡ℎ regime’s interval and 𝑛 as the time level. For given positive integers M and N representing the numbers of grid points and time steps, respectively, we have (𝑥 𝑚 ) 𝑖 = 𝑖ℎ, (𝑥 𝑙 ) 𝑗 = 𝑗ℎ, 𝜏 𝑛 = 𝑛𝑘, ℎ = (𝑥 𝑚 ) 𝑀 𝑀 , 𝑘 = 𝑇𝑁 , 𝑖, 𝑗 ∈ [0, M], 𝑘 ∈ [0, N]. (33)
We denote the numerical solutions of 𝑈 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 ), 𝑈 𝑙 ((𝑥 𝑙 ) 𝑗 , 𝜏 𝑛 ), 𝑊 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 ), 𝑊 𝑙 ((𝑥 𝑙 ) 𝑗 , 𝜏 𝑛 ), 𝑌 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 ) , 𝑌 𝑙 ((𝑥 𝑙 ) 𝑗 , 𝜏 𝑛 ) , 𝑍 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 ), 𝑍 𝑙 ((𝑥 𝑙 ) 𝑗 , 𝜏 𝑛 ), 𝑠 𝑓(𝑚) ( 𝜏 𝑛 ) , and 𝑠 𝑓(𝑙) ( 𝜏 𝑛 ) as (𝑢 𝑚 ) 𝑖𝑛 , (𝑢 𝑙 ) 𝑗𝑛 , (𝑤 𝑚 ) 𝑖𝑛 , (𝑤 𝑙 ) 𝑗𝑛 , (𝑦 𝑚 ) 𝑖𝑛 , (𝑦 𝑙 ) 𝑗𝑛 , (𝑧 𝑚 ) 𝑖𝑛 , (𝑧 𝑙 ) 𝑗𝑛 , 𝑠 𝑓(𝑚)𝑛 , and 𝑠 𝑓(𝑙)𝑛 , respectively. Compact Finite Difference Scheme
In the numerical discretization for the asset, delta, gamma and speed options in each regime, the higher-order compact finite difference method is used in space, while the second-order Crank-Nicolson method is used in time. To develop a compact finite difference scheme in space at (𝑥 𝑚 ) = 0 , we first derive a compact finite difference formula as described in the following lemma. Lemma . Assume 𝑓(𝑥) ∈ 𝐶 [𝑥 , 𝑥 ] , then it holds
74 𝑓 ′′ (𝑥 ) + 34 𝑓 ′′ (𝑥 ) = 5ℎ [𝑓(𝑥 ) − 𝑓(𝑥 )] − 5ℎ 𝑓 ′ (𝑥 ) − ℎ4 𝑓 (3) (𝑥 ) + ℎ6 𝑓 (3) (𝑥 ) + 𝑂(ℎ ). (34) Proof . Applying the Taylor expansion at 𝑥 , we obtain ℎ12 𝑓 ′′ (𝑥 ) − ℎ12 𝑓 ′′ (𝑥 ) = ℎ
12 𝑓 (3) (𝑥 ) + ℎ
24 𝑓 (4) (𝑥 ) + ℎ
72 𝑓 (5) (𝑥 ) + ⋯, (35𝑎)
53 [𝑓(𝑥 ) − 𝑓(𝑥 )ℎ ] − 53 𝑓 ′ (𝑥 ) − 5ℎ6 𝑓 ′′ (𝑥 ) = 5ℎ
18 𝑓 (3) (𝑥 ) + 5ℎ
72 𝑓 (4) (𝑥 ) + ℎ
72 𝑓 (5) (𝑥 ) + ⋯ . (35𝑏) Eliminating the fifth-order derivative by subtracting (35b) from (35a), we obtain − 53 [𝑓(𝑥 ) − 𝑓(𝑥 )ℎ ] + 53 𝑓 ′ (𝑥 ) + ℎ12 𝑓 ′′ (𝑥 ) + 9ℎ12 𝑓 ′′ (𝑥 ) = − 7ℎ
36 𝑓 (3) (𝑥 ) − ℎ
36 𝑓 (4) (𝑥 ) + 𝑂(ℎ ). (36) On the other hand, we have (Hirsh, 1975)
23 𝑓 ′ (𝑥 ) + 13 𝑓 ′ (𝑥 ) = [𝑓(𝑥 ) − 𝑓(𝑥 )ℎ ] − ℎ6 𝑓 ′′ (𝑥 ) + 𝑂(ℎ ), (37𝑎) for which we multiply it by ℎ /6 , differentiate twice, and then rearrange. This gives ℎ6 𝑓 ′′ (𝑥 ) − ℎ6 𝑓 ′′ (𝑥 ) = − 2ℎ
18 𝑓 (3) (𝑥 ) − ℎ
18 𝑓 (3) (𝑥 ) − ℎ
36 𝑓 (4) (𝑥 ) + 𝑂(ℎ ). (37𝑏) We then subtract (37b) from (36) to eliminate the fourth-order derivative and obtain − 53 [𝑓(𝑥 ) − 𝑓(𝑥 )ℎ ] + 53 𝑓 ′ (𝑥 ) + 7ℎ12 𝑓 ′′ (𝑥 ) + 3ℎ12 𝑓 ′′ (𝑥 ) = − 3ℎ
36 𝑓 (3) (𝑥 ) + ℎ
18 𝑓 (3) (𝑥 ) + 𝑂(ℎ ). (38) Dividing (38) by ℎ/3 and rearranging the terms give Eq. (34) and hence the proof has been completed. We now use (34) for the second-order derivative term in (30a) and obtain − 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 ] + 𝑂(ℎ ). (39) To evaluate the first-order derivative in (39) at (𝑥 𝑚 ) , we use (15), (20) and obtain 𝜕𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )𝜕𝑥 𝑚 − 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) = 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) = −𝐾. (40𝑎) To evaluate the third-order derivative term in (39) at (𝑥 𝑚 ) , we let 𝑥 𝑚 → 0 + in (30b) and discretize 𝜕𝑊 𝑚 /𝜕𝑡. This gives 𝜎 ℎ8 𝜕 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )𝜕𝑥 𝑚3 = ℎ4 [𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1 ) − 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛 )𝑘 − 𝜔 𝑛+1/2 𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) + ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑊 𝑙 ((𝑥 𝑚 ) 𝑗 ∗ |𝑖=0 , 𝜏 𝑛+1/2 )] + 𝑂(𝑘 ). (40𝑏) Equivalently, the third-order derivative at (𝑥 𝑚 ) is evaluated as follows: 𝜎 ℎ12 𝜕 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )𝜕𝑥 𝑚3 = ℎ6 [𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1 ) − 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛 )𝑘 − 𝜔 𝑛+1/2 𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) + ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − ∑ 𝑞 𝑚𝑙𝑙≠𝑚 𝑊 𝑙 ((𝑥 𝑚 ) 𝑗 ∗ |𝑖=1 , 𝜏 𝑛+1/2 )] + 𝑂(𝑘 ). (40𝑐) Here, (𝑥 𝑚 ) 𝑗 ∗ |𝑖=0 and (𝑥 𝑚 ) 𝑗 ∗ |𝑖=1 are the locations in the space for the 𝑙 𝑡ℎ equation corresponding to (𝑥 𝑚 ) and (𝑥 𝑚 ) in the 𝑚 𝑡ℎ equation, respectively, and (𝜔 𝑚 ) 𝑛+1/2 ≡ 2(𝑠 𝑓(𝑚)𝑛+1 − 𝑠 𝑓(𝑚)𝑛 )𝑘(𝑠 𝑓(𝑚)𝑛+1 + 𝑠 𝑓(𝑚)𝑛 ) + 𝑟 𝑚 − 𝜎
0 For the term 𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) in (30b), we employ a fourth-order approximation (Adam, 1975; Adam, 1976, Liao and Khaliq, 2009; Dremkova and Ehrhardt, 2011) as 𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )= 3ℎ [𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )] − 4𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )+ 𝑂(ℎ ). (42) Substituting (40) and (42) into (39), we obtain − 12 𝜎 [74 𝜕 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )𝜕𝑥 𝑚2 + 34 𝜕 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )𝜕𝑥 𝑚2 ]= 5𝜎 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )ℎ ] + 5𝜎 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝐾ℎ ]+ ℎ4 [[𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1 ) − 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛 )]𝑘 ] − ℎ6 [[𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1 ) − 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛 )]𝑘 ]− 34 (𝜔 𝑚 ) 𝑛+1/2 (𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ))+ ℎ24 (𝜔 𝑚 ) 𝑛+1/2 (14𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) + 3𝑌 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ))+ ℎ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )4 𝑈 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 ) − ℎ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )4 𝐾 − ℎ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )6 𝑊 𝑚 ((𝑥 𝑚 ) , 𝜏 𝑛+1/2 )− ℎ12 ∑ 𝑞 𝑚𝑙𝑙≠𝑚 (3𝑊 𝑙 ((𝑥 𝑚 ) 𝑗 ∗ |𝑖=0 , 𝜏 𝑛+1/2 ) − 2𝑊 𝑙 ((𝑥 𝑚 ) 𝑗 ∗ |𝑖=1 , 𝜏 𝑛+1/2 )) + 𝑂(ℎ ). (43) Thus, applying the Crank-Nicholson method in time for (43) , we obtain the compact finite difference scheme at (𝑥 𝑚 ) = 0 as 𝑚 ) − (𝑢 𝑚 ) 𝑘 ] + 34 [(𝑢 𝑚 ) − (𝑢 𝑚 ) 𝑘 ] − 5𝜎
4ℎ [(𝑢 𝑚 ) − (𝑢 𝑚 ) ℎ − (𝑢 𝑚 ) ]− 5𝜎
4ℎ [(𝑢 𝑚 ) − (𝑢 𝑚 ) ℎ − (𝑢 𝑚 ) ] − 𝐾 [5ℎ( 𝑟 𝑚 − 𝑞 𝑚𝑚 )4 + 5𝜎
2ℎ ]+ (𝑟 𝑚 − 𝑞 𝑚𝑚 )8 [7 + ℎ[(𝑢 𝑚 ) + (𝑢 𝑚 ) ] + 3[(𝑢 𝑚 ) + (𝑢 𝑚 ) ]] − ℎ6 [(𝑤 𝑚 ) − (𝑤 𝑚 ) 𝑘 ]− ℎ(𝑟 𝑚 − 𝑞 𝑚𝑚 )12 [(𝑤 𝑚 ) + (𝑤 𝑚 ) ]− (𝜔 𝑚 ) 𝑛+1/2 𝑚 ) + (𝑤 𝑚 ) ] + 3[(𝑤 𝑚 ) + (𝑤 𝑚 ) ] + 3[(𝑤 𝑚 ) + (𝑤 𝑚 ) ]]+ ℎ(𝜔 𝑚 ) 𝑛+1/2
24 [14[(𝑦 𝑚 ) + (𝑦 𝑚 ) ] + 3[(𝑦 𝑚 ) + (𝑦 𝑚 ) ]] − ℎ24 ∑ 𝑞 𝑚𝑙𝑙≠𝑚 [3((𝑤 𝑙 ) 𝑗 ∗ |𝑖=0𝑛+1 + (𝑤 𝑙 ) 𝑗 ∗ |𝑖=0𝑛 ) − 2((𝑤 𝑙 ) 𝑗 ∗ |𝑖=1𝑛+1 + (𝑤 𝑙 ) 𝑗 ∗ |𝑖=1𝑛 )]− 18 ∑ 𝑞 𝑚𝑙 [7((𝑢 𝑙 ) 𝑗 ∗ |𝑖=0𝑛+1 + (𝑢 𝑙 ) 𝑗 ∗ |𝑖=0𝑛 ) + 3((𝑢 𝑙 ) 𝑗 ∗ |𝑖=1𝑛+1 + (𝑢 𝑙 ) 𝑗 ∗ |𝑖=1𝑛 )] 𝑙≠𝑚 = 0, (44) 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 (Zhao et al., 2007; Liao and Khaliq, 2009; Cao et al., 2011; Gao and Sun, 2013; Yan et al., 2019) as
112 𝑓 ′′ (𝑥 𝑖−1 ) + 1012 𝑓 ′′ (𝑥 𝑖 ) + 112 𝑓 ′′ (𝑥 𝑖+1 ) = 1ℎ [𝑓(𝑥 𝑖−1 ) − 2𝑓(𝑥 𝑖 ) + 𝑓(𝑥 𝑖+1 )] + 𝑂(ℎ ), (45) for (30a)-(30d), we obtain
112 [(𝑢 𝑚 ) 𝑖−1𝑛+1 − (𝑢 𝑚 ) 𝑖−1𝑛 𝑘 ] + 1012 [(𝑢 𝑚 ) 𝑖𝑛+1 − (𝑢 𝑚 ) 𝑖𝑛 𝑘 ] + 112 [(𝑢 𝑚 ) 𝑖+1𝑛+1 − (𝑢 𝑚 ) 𝑖+1𝑛 𝑘 ]− 𝜎 [(𝑢 𝑚 ) 𝑖−1𝑛+1 − 2(𝑢 𝑚 ) 𝑖𝑛+1 + (𝑢 𝑚 ) 𝑖+1𝑛+1 ] − 𝜎 [(𝑢 𝑚 ) 𝑖−1𝑛 − 2(𝑢 𝑚 ) 𝑖𝑛 + (𝑢 𝑚 ) 𝑖+1𝑛 ]− (𝜔 𝑚 ) 𝑛+1/2
24 [[(𝑤 𝑚 ) 𝑖−1𝑛+1 + (𝑤 𝑚 ) 𝑖−1𝑛 ] + 10[(𝑤 𝑚 ) 𝑖𝑛+1 + (𝑤 𝑚 ) 𝑖𝑛 ] + [(𝑤 𝑚 ) 𝑖+1𝑛+1 + (𝑤 𝑚 ) 𝑖+1𝑛 ]]+ (𝑟 𝑚 − 𝑞 𝑚𝑚 )24 [[(𝑢 𝑚 ) 𝑖−1𝑛+1 + (𝑢 𝑚 ) 𝑖−1𝑛 ] + 10[(𝑢 𝑚 ) 𝑖𝑛+1 + (𝑢 𝑚 ) 𝑖𝑛 ] + [(𝑢 𝑚 ) 𝑖+1𝑛+1 + (𝑢 𝑚 ) 𝑖+1𝑛 ]]− ∑ 𝑞 𝑚𝑙
24 [[(𝑢 𝑙 ) 𝑗 ∗ −1𝑛+1 + (𝑢 𝑙 ) 𝑗 ∗ −1𝑛 ] + 10[(𝑢 𝑙 ) 𝑗 ∗ 𝑛+1 + (𝑢 𝑙 ) 𝑗 ∗ 𝑛 ] + [(𝑢 𝑙 ) 𝑗 ∗ +1𝑛+1 + (𝑢 𝑙 ) 𝑗 ∗ +1𝑛 ]] 𝑙≠𝑚 = 0, (46𝑎)
112 [(𝑤) 𝑖−1𝑛+1 − (𝑤 𝑚 ) 𝑖−1𝑛 𝑘 ] + 1012 [(𝑤 𝑚 ) 𝑖𝑛+1 − (𝑤 𝑚 ) 𝑖𝑛 𝑘 ] + 112 [(𝑤 𝑚 ) 𝑖+1𝑛+1 − (𝑤 𝑚 ) 𝑖+1𝑛 𝑘 ]− 𝜎 [(𝑤 𝑚 ) 𝑖−1𝑛+1 − 2(𝑤 𝑚 ) 𝑖𝑛+1 + (𝑤 𝑚 ) 𝑖+1𝑛+1 ] − 𝜎 [(𝑤 𝑚 ) 𝑖−1𝑛 − 2(𝑤 𝑚 ) 𝑖𝑛 + (𝑤 𝑚 ) 𝑖+1𝑛 ]− (𝜔 𝑚 ) 𝑛+1/2 [[(𝑢 𝑚 ) 𝑖−1𝑛+1 + (𝑢 𝑚 ) 𝑖−1𝑛 ] − 2[(𝑢 𝑚 ) 𝑖𝑛+1 + (𝑢 𝑚 ) 𝑖𝑛 ] + [(𝑢 𝑚 ) 𝑖+1𝑛+1 + (𝑢 𝑚 ) 𝑖+1𝑛 ]]+ (𝑟 𝑚 − 𝑞 𝑚𝑚 )24 [[(𝑤 𝑚 ) 𝑖−1𝑛+1 + (𝑤 𝑚 ) 𝑖−1𝑛 ] + 10[(𝑤 𝑚 ) 𝑖𝑛+1 + (𝑤 𝑚 ) 𝑖𝑛 ] + [(𝑤 𝑚 ) 𝑖+1𝑛+1 + (𝑤 𝑚 ) 𝑖+1𝑛 ]]− ∑ 𝑞 𝑚𝑙 𝑙≠𝑚 [[(𝑤 𝑙 ) 𝑗 ∗ −1𝑛+1 + (𝑤 𝑙 ) 𝑗 ∗ −1𝑛 ] + 10[(𝑤 𝑙 ) 𝑗 ∗ 𝑛+1 + (𝑤 𝑙 ) 𝑗 ∗ 𝑛 ] + [(𝑤 𝑙 ) 𝑗 ∗ +1𝑛+1 + (𝑤 𝑙 ) 𝑗 ∗ +1𝑛 ]]= 0, (46𝑏)
112 [(𝑦 𝑚 ) 𝑖−1𝑛+1 − (𝑦 𝑚 ) 𝑖−1𝑛 𝑘 ] + 1012 [(𝑦 𝑚 ) 𝑖𝑛+1 − (𝑦 𝑚 ) 𝑖𝑛 𝑘 ] + 112 [(𝑦 𝑚 ) 𝑖+1𝑛+1 − (𝑦 𝑚 ) 𝑖+1𝑛 𝑘 ] − 𝜎 [(𝑦 𝑚 ) 𝑖−1𝑛+1 − 2(𝑦 𝑚 ) 𝑖𝑛+1 + (𝑦 𝑚 ) 𝑖+1𝑛+1 ] − 𝜎 [(𝑦 𝑚 ) 𝑖−1𝑛 − 2(𝑦 𝑚 ) 𝑖𝑛 + (𝑦 𝑚 ) 𝑖+1𝑛 ]− (𝜔 𝑚 ) 𝑛+1/2 [[(𝑤 𝑚 ) 𝑖−1𝑛+1 + (𝑤 𝑚 ) 𝑖−1𝑛 ] − 2[(𝑤 𝑚 ) 𝑖𝑛+1 + (𝑤 𝑚 ) 𝑖𝑛 ] + [(𝑤 𝑚 ) 𝑖+1𝑛+1 + (𝑤 𝑚 ) 𝑖+1𝑛 ]]+ (𝑟 𝑚 − 𝑞 𝑚𝑚 )24 [[(𝑦 𝑚 ) 𝑖−1𝑛+1 + (𝑦 𝑚 ) 𝑖−1𝑛 ] + 10[(𝑦 𝑚 ) 𝑖𝑛+1 + (𝑦 𝑚 ) 𝑖𝑛 ] + [(𝑦 𝑚 ) 𝑖+1𝑛+1 + (𝑦 𝑚 ) 𝑖+1𝑛 ]]− ∑ 𝑞 𝑚𝑙 𝑙≠𝑚 [[(𝑦 𝑙 ) 𝑗 ∗ −1𝑛+1 + (𝑦 𝑙 ) 𝑗 ∗ −1𝑛 ] + 10[(𝑦 𝑙 ) 𝑗 ∗ 𝑛+1 + (𝑦 𝑙 ) 𝑗 ∗ 𝑛 ] + [(𝑦 𝑙 ) 𝑗 ∗ +1𝑛+1 + (𝑦 𝑙 ) 𝑗 ∗ +1𝑛 ]]= 0, (46𝑐)
112 [(𝑧 𝑚 ) 𝑖−1𝑛+1 − (𝑧 𝑚 ) 𝑖−1𝑛 𝑘 ] + 1012 [(𝑧 𝑚 ) 𝑖𝑛+1 − (𝑧 𝑚 ) 𝑖𝑛 𝑘 ] + 112 [(𝑧 𝑚 ) 𝑖+1𝑛+1 − (𝑧 𝑚 ) 𝑖+1𝑛 𝑘 ]− 𝜎 [(𝑧 𝑚 ) 𝑖−1𝑛+1 − 2(𝑧 𝑚 ) 𝑖𝑛+1 + (𝑧 𝑚 ) 𝑖+1𝑛+1 ] − 𝜎 [(𝑧 𝑚 ) 𝑖−1𝑛 − 2(𝑧 𝑚 ) 𝑖𝑛 + (𝑧 𝑚 ) 𝑖+1𝑛 ]− 𝜔 𝑛+1/2 [[(𝑦 𝑚 ) 𝑖−1𝑛+1 + (𝑦 𝑚 ) 𝑖−1𝑛 ] − 2[(𝑦 𝑚 ) 𝑖𝑛+1 + (𝑦 𝑚 ) 𝑖𝑛 ] + [(𝑦 𝑚 ) 𝑖+1𝑛+1 + (𝑦 𝑚 ) 𝑖+1𝑛 ]]+ (𝑟 𝑚 − 𝑞 𝑚𝑚 )24 [[(𝑧 𝑚 ) 𝑖−1𝑛+1 + (𝑧 𝑚 ) 𝑖−1𝑛 ] + 10[(𝑧 𝑚 ) 𝑖𝑛+1 + (𝑧 𝑚 ) 𝑖𝑛 ] + [(𝑧 𝑚 ) 𝑖+1𝑛+1 + (𝑧 𝑚 ) 𝑖+1𝑛 ]]− ∑ 𝑞 𝑚𝑙 𝑙≠𝑚 [[(𝑧 𝑙 ) 𝑗 ∗ −1𝑛+1 + (𝑧 𝑙 ) 𝑗 ∗ −1𝑛 ] + 10[(𝑧 𝑙 ) 𝑗 ∗ 𝑛+1 + (𝑧 𝑙 ) 𝑗 ∗ 𝑛 ] + [(𝑧 𝑙 ) 𝑗 ∗ +1𝑛+1 + (𝑧 𝑙 ) 𝑗 ∗ +1𝑛 ]]= 0, (46𝑑) 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 ; (47𝑎 ) (𝑢 𝑚 ) 𝑀𝑛+1 = 0, (𝑤 𝑚 ) 𝑀𝑛+1 = 0, (𝑦 𝑚 ) 𝑀𝑛+1 = 0, (𝑧 𝑚 ) 𝑀𝑛+1 = 0; (47𝑏 ) (𝑢 𝑚 ) 𝑖0 = (𝑤 𝑚 ) 𝑖0 = (𝑦 𝑚 ) 𝑖0 = (𝑧 𝑚 ) 𝑖0 = 0, 𝑖 = 1,2,∙∙∙, 𝑀. (47𝑐) Let the approximate solutions of the theta, delta decay, and color options for each regime be given as 𝜕𝑈 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 )𝜕𝜏 ≈ (Θ 𝑚 ) 𝑖𝑛 , 𝜕𝑊 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 )𝜕𝜏 ≈ (Κ 𝑚 ) 𝑖𝑛 , 𝜕𝑌 𝑚 ((𝑥 𝑚 ) 𝑖 , 𝜏 𝑛 )𝜕𝜏 ≈ (Γ 𝑚 ) 𝑖𝑛 ; (48) respectively. For 𝑛 = 1 , we approximate these three Greeks using first-order backward finite differences (Θ 𝑚 ) 𝑖1 ≈ (𝑢 𝑚 ) 𝑖1 − (𝑢 𝑚 ) 𝑖0 𝑘 , (Κ 𝑚 ) 𝑖1 ≈ (𝑤 𝑚 ) 𝑖1 − (𝑤 𝑚 ) 𝑖0 𝑘 , (Γ 𝑚 ) 𝑖1 ≈ (𝑦 𝑚 ) 𝑖1 − (𝑦) 𝑖0 𝑘 . (49𝑎) Subsequently, we use the second-order backward finite difference approximations as (Θ 𝑚 ) 𝑖𝑛+1 ≈ 3(𝑢 𝑚 ) 𝑖𝑛+1 − 4(𝑢 𝑚 ) 𝑖𝑛 + (𝑢 𝑚 ) 𝑖𝑛−1
2𝑘 , (Κ 𝑚 ) 𝑖𝑛+1 ≈ 3(𝑤) 𝑖𝑛+1 − 4(𝑤 𝑚 ) 𝑖𝑛 + (𝑤 𝑚 ) 𝑖𝑛−1
2𝑘 ; (49𝑏) (Γ 𝑚 ) 𝑖𝑛+1 ≈ 3(𝑦) 𝑖𝑛+1 − 4(𝑦 𝑚 ) 𝑖𝑛 + (𝑦 𝑚 ) 𝑖𝑛−1
2𝑘 . (49𝑐)
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,∙∙∙, 𝑀. (50) Hermite Interpolation
To evaluate (𝑢 𝑙 ) 𝑗 ∗ 𝑛 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 , and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 in (46), we need to consider the relationship between the fixed interval (and the mesh) for the 𝑙 𝑡ℎ regime and the fixed interval (and the mesh) for the 𝑚 𝑡ℎ regime after the logarithmic transformation. If 𝑠 𝑓(𝑙) (𝜏 𝑛 ) = 𝑠 𝑓(𝑚) (𝜏 𝑛 ) , then the fixed interval for the 𝑙 𝑡ℎ regime overlaps completely with the fixed interval for the 𝑚 𝑡ℎ regime. Hence, (𝑥 𝑚 ) 𝑖 = (𝑥 𝑙 ) 𝑗 . If 𝑠 𝑓(𝑙) (𝜏 𝑛 ) ≠ 𝑠 𝑓(𝑚) (𝜏 𝑛 ) , then there are three possible cases as shown in Fig. 1. Fig. 1 . Relationship between the 𝑙 𝑡ℎ and 𝑚 𝑡ℎ intervals and the location of the (𝑥 𝑚 ) 𝑖 in the 𝑙 𝑡ℎ interval. 4 Fig 1a shows that there exists a possibility for (𝑥 𝑚 ) 𝑖 corresponding to (𝑥 𝑙 ) 𝑗 ∗ < 0 in the 𝑙 𝑡ℎ interval. For this case, (𝑢 𝑙 ) 𝑖𝑛 = 𝐾 − 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ and (𝑤 𝑙 ) 𝑗 ∗ 𝑛 = (𝑦 𝑙 ) 𝑗 ∗ 𝑛 = (𝑧) 𝑗 ∗ 𝑛 = −𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ based on (13)-(16). Fig 1b shows that there exists a possibility for (𝑥 𝑚 ) 𝑖 corresponding to a point in (0, ( 𝑥 𝑙 ) 𝑀 ) . For this case, (𝑢 𝑙 ) 𝑗 ∗ 𝑛 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 have to be evaluated using an interpolation based on (𝑢 𝑙 ) 𝑗𝑛 , (𝑤 𝑙 ) 𝑗𝑛 , (𝑦 𝑙 ) 𝑗𝑛 and (𝑧 𝑙 ) 𝑗𝑛 . Here, we employ the Hermite interpolation (Burden et al., 2016) to ensure higher-order accuracy. Fig. 1c shows that there exists a possibility for (𝑥 𝑚 ) 𝑖 corresponding to (𝑥 𝑙 ) 𝑗 ∗ >(𝑥 𝑙 ) 𝑀 . For this case, we set (𝑢 𝑙 ) 𝑗 ∗ 𝑛 = (𝑤 𝑙 ) 𝑗 ∗ 𝑛 = (𝑦 𝑙 ) 𝑗 ∗ 𝑛 = (𝑧 𝑙 ) 𝑗 ∗ 𝑛 = 0 . In overall, we can evaluate (𝑢 𝑙 ) 𝑗 ∗ 𝑛 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 based on the following formulas as: (𝑢 𝑙 ) 𝑗 ∗ 𝑛 = { 𝐾 − 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑐 (𝑢 𝑙 ) 𝑗𝑛 + 𝑏 𝑐 (𝑢 𝑙 ) 𝑗+1𝑛 + 𝑐 𝑐 (𝑤 𝑙 ) 𝑗𝑛 + 𝑑 𝑐 (𝑤 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − 𝑙𝑛 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (51𝑎) (𝑤 𝑙 ) 𝑗 ∗ 𝑛 = { −𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑒 𝑐 (𝑢 𝑙 ) 𝑗𝑛 + 𝑓 𝑐 (𝑢 𝑙 ) 𝑗+1𝑛 + 𝑔 𝑐 (𝑤 𝑙 ) 𝑗𝑛 + 𝑜 𝑐 (𝑤 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − 𝑙𝑛 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (51𝑏) (𝑦 𝑙 ) 𝑗 ∗ 𝑛 = { −𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑐 (𝑦 𝑙 ) 𝑗𝑛 + 𝑏 𝑐 (𝑦 𝑙 ) 𝑗+1𝑛 + 𝑐 𝑐 (𝑧 𝑙 ) 𝑗𝑛 + 𝑑 𝑐 (𝑧 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − 𝑙𝑛 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (51𝑐) (𝑦 𝑙 ) 𝑗 ∗ 𝑛 = { −𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑐 (𝑧 𝑙 ) 𝑗𝑛 + 𝑏 𝑐 (𝑧 𝑙 ) 𝑗+1𝑛 + 𝑐 𝑐 (𝑧̇) 𝑗𝑛 + 𝑑 𝑐 (𝑧̇) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − 𝑙𝑛 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (51𝑑) where the coefficients are given based on the cubic Hermite Interpolation as follows: 𝑎 𝑐 = 1ℎ [1 + 2[(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ]ℎ ] [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ] , (52𝑎) 𝑏 𝑐 = 1ℎ [1 − 2[(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ]ℎ ] [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] , (52𝑏) 𝑐 𝑐 = 1ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ][(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ] , 𝑑 𝑐 = 1ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ][(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] , (52𝑐) 𝑒 𝑐 = 2ℎ [1 + 2[(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ]ℎ ] [(𝑥 𝑙 ) 𝑗 − (𝑥 𝑙 ) 𝑗+1 ] + 2ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ] , (52𝑑) 𝑓 𝑐 = 2ℎ [1 − 2[(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ]ℎ ] [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] + − 2ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] , (52𝑒) 𝑔 𝑐 = 2ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ][(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ] + 1ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ] , (52𝑓) 𝑜 𝑐 = 2ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ][(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] + 1ℎ [(𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ] . (52𝑔) Alternatively, one may use the higher accurate quintic Hermite interpolation as (𝑢 𝑙 ) 𝑗 ∗ 𝑛 = { 𝐾 − 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑞 (𝑢 𝑙 ) 𝑗−1𝑛 + 𝑏 𝑞 (𝑢 𝑙 ) 𝑗𝑛 + 𝑐 𝑞 (𝑢 𝑙 ) 𝑗+1𝑛 + 𝑑 𝑞 (𝑤 𝑙 ) 𝑗−1𝑛 + 𝑒 𝑞 (𝑤 𝑙 ) 𝑗𝑛 +𝑓 𝑞 (𝑤 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝑡 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (53𝑎) (𝑤 𝑙 ) 𝑗 ∗ 𝑛 = { −𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑔 𝑞 (𝑤 𝑙 ) 𝑗−1𝑛 + 𝑜 𝑞 (𝑤 𝑙 ) 𝑗𝑛 + 𝑝 𝑞 (𝑤 𝑙 ) 𝑗+1𝑛 + 𝑞 𝑞 (𝑦 𝑙 ) 𝑗−1𝑛 + 𝑟 𝑞 (𝑦 𝑙 ) 𝑗𝑛 +𝑠 𝑞 (𝑦 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (53𝑏) (𝑦 𝑙 ) 𝑗 ∗ 𝑛 = { 𝐾 − 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑞 (𝑦 𝑙 ) 𝑗−1𝑛 + 𝑏 𝑞 (𝑦 𝑙 ) 𝑗𝑛 + 𝑐 𝑞 (𝑧 𝑙 ) 𝑗+1𝑛 + 𝑑 𝑞 (𝑧 𝑙 ) 𝑗−1𝑛 + 𝑒 𝑞 (𝑧 𝑙 ) 𝑗𝑛 +𝑓 𝑞 (𝑧 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (53𝑐) (𝑧 𝑙 ) 𝑗 ∗ 𝑛 = { 𝐾 − 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑒 (𝑥 𝑙 ) 𝑗∗ , (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ 0; 𝑎 𝑞 (𝑧 𝑙 ) 𝑗−1𝑛 + 𝑏 𝑞 (𝑧 𝑙 ) 𝑗𝑛 + 𝑐 𝑞 (𝑧 𝑙 ) 𝑗+1𝑛 + 𝑑 𝑞 (𝑧̇ 𝑙 ) 𝑗−1𝑛 + 𝑒 𝑞 (𝑧̇ 𝑙 ) 𝑗𝑛 +𝑓 𝑞 (𝑧̇ 𝑙 ) 𝑗+1𝑛 , (𝑥 𝑙 ) 𝑗 ≤ (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) ≤ (𝑥 𝑙 ) 𝑗+1 ; 0, (𝑥 𝑚 ) 𝑖 − ln 𝑠 𝑓(𝑙) (𝜏 𝑛 )𝑠 𝑓(𝑚) (𝜏 𝑛 ) > (𝑥 𝑙 ) 𝑓𝑏 , (53𝑑) with the coefficients given as 𝑎 𝑞 = [1 + 3ℎ ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )] [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] , (54𝑎) 𝑏 𝑞 = [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] ℎ , (54𝑏) 𝑐 𝑞 = [1 − 3ℎ ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )] , (54𝑐) 𝑑 𝑞 = ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 ) [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] , (54𝑑) 𝑒 𝑞 = ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 ) [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] ℎ , (54𝑒) 𝑓 𝑞 = ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 ) [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )] , (54𝑓) 𝑔 𝑞 = [1 + 3ℎ ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )] [2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗 + (𝑥 𝑙 ) 𝑗+1 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )]2ℎ + 3[((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] , (54𝑔) 𝑜 𝑞 = 2[2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗−1 + (𝑥 𝑙 ) 𝑗+1 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )]ℎ , (54ℎ) 𝑝 𝑞 = [1 − 3ℎ ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] [2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗−1 + (𝑥 𝑙 ) 𝑗 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )]2ℎ − 3[((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )] , (54𝑖) 𝑞 𝑞 = ((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )[2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗 + (𝑥 𝑙 ) 𝑗+1 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )]2ℎ + [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] (54𝑗) 𝑟 𝑞 = 2((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )[2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗−1 + (𝑥 𝑙 ) 𝑗+1 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )]ℎ + [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )] ℎ (54𝑘) 𝑠 𝑞 = 2((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗+1 )[2(𝑥 𝑙 ) 𝑗 ∗ − ((𝑥 𝑙 ) 𝑗−1 + (𝑥 𝑙 ) 𝑗 )][((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )]2ℎ + [((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗−1 )((𝑥 𝑙 ) 𝑗 ∗ − (𝑥 𝑙 ) 𝑗 )] . (54𝑙) The subscripts “c” and “q” denote the cubic and quintic Hermite interpolations, respectively. It should be pointed out that we have compared with other high-order interpolations when estimating these (𝑢 𝑙 ) 𝑗 ∗ 𝑛 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 . Hermite interpolation proves to be accurate, more efficient in handling large state space and very fast in computation. Moreover, it is worth noting that the derivative of (𝑧 𝑙 ) 𝑗 ∗ 𝑛 , (𝑧̇) 𝑗 ∗ 𝑛 is employed in the cubic and quintic Hermite interpolations. To evaluate (𝑧̇) 𝑗 ∗ 𝑛 with fourth-order accuracy, we obtain it from the speed option PDE by further taking derivative. To carry out an extensive 8 analysis, we further investigate the performance of both cubic and quintic Hermite interpolations in the numerical example section using both Gauss-Seidel and Newton iterative methods. Stability Analysis
The stability analysis of our numerical schemes is carried out using the matrix form of von Neumann method (see Hirsch (2001) and Liao and Khaliq (2009)). Due to the complex system of the present method, we ignore the coupled regimes (𝑢 𝑙 ) 𝑖𝑛 and (𝑤 𝑙 ) 𝑖𝑛 , (𝑦 𝑙 ) 𝑖𝑛 and (𝑧 𝑙 ) 𝑖𝑛 and let (𝑢 𝑚 ) 𝑖𝑛 = 𝜆 𝑚𝑛 𝑒 Ι𝛽𝑖ℎ , (𝑤 𝑚 ) 𝑖𝑛 = Φ 𝑚𝑛 𝑒 Ι𝛽𝑖ℎ , (𝑦 𝑚 ) 𝑖𝑛 = Υ 𝑚𝑛 𝑒 Ι𝛽𝑖ℎ , (𝑧 𝑚 ) 𝑖𝑛 = Ψ 𝑚𝑛 𝑒 Ι𝛽𝑖ℎ , Ι = √−1. (55)
Denote 𝜇 = 𝜎 𝑚2 𝑘4ℎ , 𝜅 = (𝑟 𝑚 − 𝑞 𝑚𝑚 )𝑘, 𝜔 = ( ( 𝑠 𝑓 ( 𝑚 ) 𝑛+1 − 𝑠 𝑓 ( 𝑚 ) 𝑛 )( 𝑠 𝑓 ( 𝑚 ) 𝑛+1 + 𝑠 𝑓 ( 𝑚 ) 𝑛 ) 𝑘 + 𝑟 𝑚 − 𝜎 ) 𝑘. (56) Substituting (55), (56) into (46a)-(46d), we obtain 𝜆 𝑚𝑛+1 [1 − 13 sin (𝛽ℎ2 ) + 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )] −𝜆 𝑚𝑛 [1 − 13 sin (𝛽ℎ2 ) − 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− 𝜔 [Φ 𝑚𝑛+1 (12 − 16 sin (𝛽ℎ2 )) + Φ 𝑚𝑛 (12 − 16 sin (𝛽ℎ2 ))] = 0, (57𝑎) Φ 𝑚𝑛+1 [1 − 13 sin (𝛽ℎ2 ) + 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 ) 7]− Φ 𝑚𝑛 [1 − 13 sin (𝛽ℎ2 ) − 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− 𝜔 [−4𝜆 𝑚𝑛+1 sin (𝛽ℎ2 ) −4𝜆 𝑚𝑛 sin (𝛽ℎ2 )] = 0, (57𝑏) Υ 𝑚𝑛+1 [1 − 13 sin (𝛽ℎ2 ) + 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− Υ 𝑚𝑛 [1 − 13 sin (𝛽ℎ2 ) − 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− 𝜔 [−4Φ 𝑚𝑛+1 sin (𝛽ℎ2 ) −4Φ 𝑚𝑛 sin (𝛽ℎ2 )] = 0, (57𝑐) Ψ 𝑚𝑛+1 [1 − 13 sin (𝛽ℎ2 ) + 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− Ψ 𝑚𝑛 [1 − 13 sin (𝛽ℎ2 ) − 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 )]− 𝜔 [−4Υ 𝑚𝑛+1 sin (𝛽ℎ2 ) −4Υ 𝑚𝑛 sin (𝛽ℎ2 )] = 0, (57𝑑) which can be simplified to 𝑝𝜆 𝑚𝑛+1 − 𝑞𝜆 𝑚𝑛 = 𝑟Φ 𝑚𝑛+1 + 𝑟Φ 𝑚𝑛 , 𝑝Φ 𝑚𝑛+1 − 𝑞Φ 𝑚𝑛 = 𝑠𝜆 𝑚𝑛+1 + 𝑠𝜆 𝑚𝑛 , (58𝑎) 𝑝Υ 𝑚𝑛+1 − 𝑞Υ 𝑚𝑛 = 𝑠Φ 𝑚𝑛+1 + 𝑠Φ 𝑚𝑛 , 𝑝Ψ 𝑚𝑛+1 − 𝑞Ψ 𝑚𝑛 = 𝑠Υ 𝑚𝑛+1 + 𝑠Υ 𝑚𝑛 , (58𝑏) where 𝑝 = 1 − 13 sin (𝛽ℎ2 ) + 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 ) , 𝑟 = − 2𝜔ℎ sin (𝛽ℎ2 ), (59𝑎) 𝑠 = 𝜔 [12 − 16 sin (𝛽ℎ2 )] , 𝑞 = 1 − 13 sin (𝛽ℎ2 ) − 4 𝜇 sin (𝛽ℎ2 ) + 𝜅2 − 𝜅2 sin (𝛽ℎ2 ). (59𝑏) We then obtain a system of equations from (58) and present it in matrix-vector form as [ 𝜆 𝑚𝑛+1 Φ 𝑚𝑛+1 Υ 𝑚𝑛+1 Ψ 𝑚𝑛+1 ] = [𝑝 − 𝑟 0 0−𝑠 𝑝 0 00 − 𝑠 𝑝 00 0 − 𝑠 𝑝] −1 [ 𝑞 𝑟 0 0𝑠 𝑞 0 00 𝑠 𝑞 00 0 𝑠 𝑞 ] [ 𝜆 𝑚 𝑛 Φ 𝑚𝑛 Υ 𝑚𝑛 Ψ 𝑚𝑛 ] = [ 𝑝𝑝 − 𝑠𝑟 𝑟𝑝 − 𝑠𝑟
0 0 𝑠𝑝 − 𝑠𝑟 𝑝𝑝 − 𝑠𝑟
0 0 𝑠 𝑝(𝑝 − 𝑠𝑟) 𝑠𝑝 − 𝑠𝑟
1𝑝 0 𝑠 𝑝 (𝑝 − 𝑠𝑟) 𝑠 𝑝(𝑝 − 𝑠𝑟) 𝑠𝑝 ] [ 𝑞 𝑟 0 0𝑠 𝑞 0 00 𝑠 𝑞 00 0 𝑠 𝑞 ] [ 𝜆 𝑚 𝑛 Φ 𝑚𝑛 Υ 𝑚𝑛 Ψ 𝑚𝑛 ]= [ 𝑝𝑞 + 𝑟𝑠𝑝 − 𝑠𝑟 𝑝𝑟 + 𝑞𝑟𝑝 − 𝑠𝑟
0 0 𝑝𝑠 + 𝑞𝑠𝑝 − 𝑠𝑟 𝑝𝑞 + 𝑟𝑠𝑝 − 𝑠𝑟
0 0 𝑝𝑠 + 𝑞𝑠 𝑝(𝑝 − 𝑠𝑟) 𝑝𝑞𝑠 + 𝑝 𝑠𝑝(𝑝 − 𝑠𝑟) 𝑞𝑝 0 𝑝𝑠 + 𝑞𝑠 𝑝 (𝑝 − 𝑠𝑟) 𝑝𝑞𝑠 + 𝑝 𝑠 𝑝 (𝑝 − 𝑠𝑟) 𝑝𝑠 + 𝑞𝑠𝑝 𝑞𝑝 ] [ 𝜆 𝑚 𝑛 Φ 𝑚𝑛 Υ 𝑚𝑛 Ψ 𝑚𝑛 ] = 𝐴 [ 𝜆 𝑚 𝑛 Φ 𝑚𝑛 Υ 𝑚𝑛 Ψ 𝑚𝑛 ]. (60) Here, 𝐴 represents the amplification matrix. To show our numerical method to be unconditionally stable, we need to confirm that the modulus of the eigenvalue of the matrix A is less than or equal to 1 (see 0 Hirsch (2001) and Liao and Khaliq (2009)). Denoting the eigenvalue of the matrix 𝐴 as 𝜑 , we obtain the equation below [𝜑 − 2𝜑 𝑝𝑞 + 𝑟𝑠𝑝 − 𝑠𝑟 + ( 𝑝𝑞 + 𝑟𝑠 ) ( 𝑝 − 𝑠𝑟 ) − ( 𝑝𝑟 + 𝑞𝑟 ) (𝑝𝑠 + 𝑞𝑠) ( 𝑝 − 𝑠𝑟 ) ] [( 𝑞𝑝 − 𝜑 ) ( 𝑞𝑝 − 𝜑 )] = 0. (61) Note that 𝑟𝑠 = − 2𝜔 ℎ sin (𝛽ℎ2 ) [12 − 16 sin (𝛽ℎ2 )] ≤ 0, 𝑝 ≥ 𝑞. (62) Since 𝜇 > 0, we obtain 𝑝 > 𝑞 and hence (𝑞𝑝 − 𝜑) (𝑞𝑝 − 𝜑) = 0, | 𝜑 | = |𝑞𝑝| < 1. (63) Furthermore, we need to obtain 𝜑 by solving 𝜑 − 2𝜑 𝑝𝑞 + 𝑟𝑠𝑝 − 𝑠𝑟 + ( 𝑝𝑞 + 𝑟𝑠 ) ( 𝑝 − 𝑠𝑟 ) − ( 𝑝𝑟 + 𝑞𝑟 )( 𝑝𝑠 + 𝑞𝑠 )( 𝑝 − 𝑠𝑟 ) = 0, (64) which gives 𝜑 = ( 𝑝𝑞 + 𝑟𝑠 ) ± ( 𝑝 + 𝑞 )√ 𝑟𝑠𝑝 − 𝑟𝑠 . (65) Noticing 𝜛 ≡ −𝑟𝑠 ≥ 0 , we obtain the complex conjugate values of the eigenvalues as 𝜑 = ( 𝑝𝑞 − 𝜛 ) ± ( 𝑝 + 𝑞 ) 𝐼 √ 𝜛𝑝 + 𝜛 . (66) Thus, we have |𝜑 | = ( 𝑝𝑞 − 𝜛 ) + ( 𝑝 + 𝑞 ) 𝜛 ( 𝑝 + 𝜛 ) = ( 𝑝 + 𝜛 )( 𝑞 + 𝜛 )( 𝑝 + 𝜛 ) = ( 𝑞 + 𝜛 )( 𝑝 + 𝜛 ) ≤ 1. (67) Based on the von Neumann analysis, we have proved that our numerical schemes are unconditionally stable.
Computational Procedure with Gauss-Seidel Iteration
The system in (44), (46)-(54) must be solved iteratively. Here, we first present an iterative procedure based on the Gauss-Seidel (GS) iterative method (Kwok, 2011; Chapra, 2012). We initialize 𝑠 𝑓(𝑚) 𝑛 , (𝑢 𝑚 ) 𝑖𝑛 , (𝑤 𝑚 ) 𝑖𝑛 , (𝑦 𝑚 ) 𝑖𝑛 , (𝑧 𝑚 ) 𝑖𝑛 , (Θ 𝑚 ) 𝑖𝑛 , and (Κ 𝑚 ) 𝑖𝑛 where (𝑢 𝑙 ) 𝑗 ∗ 𝑛 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 are calculated based on (51)-(54). We assume that (𝑢 𝑚 ) 𝑖𝑛+1(It=0) = (𝑢 𝑚 ) 𝑖𝑛 , (𝑤 𝑚 ) 𝑖𝑛+1(It=0) = (𝑤 𝑚 ) 𝑖𝑛 , (𝑦 𝑚 ) 𝑖𝑛+1(It=0) = (𝑦 𝑚 ) 𝑖𝑛 , and (𝑧 𝑚 ) 𝑖𝑛+1(It=0) = (𝑧 𝑚 ) 𝑖𝑛 , where “It ” is the iteration counter. Next, (𝑢 𝑚 ) 𝑖𝑛+1(It=1) is computed and 𝑠 𝑓(𝑚) 𝑛+1(It=1) is obtained from (𝑢 𝑚 ) . Subsequently, we compute (𝑤 𝑚 ) 𝑖𝑛+1(It=1) , (𝑦 𝑚 ) 𝑖𝑛+1(It=1) , (𝑧 𝑚 ) 𝑖𝑛+1(It=1) , (Θ 𝑚 ) 𝑖𝑛+1(It=1) , (Κ 𝑚 ) 𝑖𝑛+1(It=1) , and (Γ 𝑚 ) 𝑖𝑛+1(It=1) . T he iterative process continues until the convergence criterion of both max 𝑚 |𝑠 𝑓(𝑚)𝑛+1(It+1) − 𝑠 𝑓(𝑚)𝑛+1(It) | <𝜀 and max 𝑚,𝑖 |(𝑢 𝑚 ) 𝑖𝑛+1(It+1) − (𝑢 𝑚 ) 𝑖𝑛+1(It) | < 𝜀 is satisfied . An algorithm for obtaining the numerical solutions of the optimal exercise boundary, asset option and the option Greeks in each regime using the GS method is described below.
Algorithm 1.
An algorithm based on the Gauss-Seidel Iteration Initialize 𝑠 𝑓(𝑚)𝑛 , (𝑢 𝑚 ) 𝑖𝑛 , (𝑤 𝑚 ) 𝑖𝑛 , (𝑦 𝑚 ) 𝑖𝑛 , (𝑧 𝑚 ) 𝑖𝑛 , (Θ 𝑚 ) 𝑖𝑛 , (Κ 𝑚 ) 𝑖𝑛 , and (Γ 𝑚 ) 𝑖𝑛 for 𝑖 = 0,1, … , 𝑀 and 𝑚 =1,2, … , 𝐼 𝐟𝐨𝐫 𝐧 = 𝟏 𝐭𝐨 𝐍 Compute (𝑢 𝑙 ) 𝑗 ∗ 𝑛 and (𝑤 𝑙 ) 𝑗 ∗ 𝑛 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛 for 𝑙 = 1,2, … , 𝐼 and 𝑙 ≠ 𝑚 based on (51)-(54) Set 𝑠 𝑓(𝑚)𝑛+1(It=0) = 𝑠 𝑓(𝑚)𝑛 , ( 𝑢 𝑚 ) 𝑖𝑛+1(It=0) = ( 𝑢 𝑚 ) 𝑖𝑛 , ( 𝑤 𝑚 ) 𝑖𝑛+1(It=0) = ( 𝑤 𝑚 ) 𝑖𝑛 , ( 𝑦 𝑚 ) 𝑖𝑛+1(It=0) = ( 𝑦 𝑚 ) 𝑖𝑛 , ( 𝑧 𝑚 ) 𝑖𝑛+1(It=0) = ( 𝑧 𝑚 ) 𝑖𝑛 𝐰𝐡𝐢𝐥𝐞 𝐭𝐫𝐮𝐞 Compute (𝑢 𝑙 ) 𝑗 ∗ 𝑛+1 , (𝑤 𝑙 ) 𝑗 ∗ 𝑛+1 , (𝑦 𝑙 ) 𝑗 ∗ 𝑛+1 and (𝑧 𝑙 ) 𝑗 ∗ 𝑛+1 for 𝑙 = 1,2, … , 𝐼 and 𝑙 ≠ 𝑚 based on (51)-(54) 𝐟𝐨𝐫 𝐦 = 𝟏 𝐭𝐨 𝐈 Compute (𝑢 𝑚 ) 𝑖𝑛+1(It+1) and evaluate 𝑠 𝑓(𝑚) 𝑛+1(It+1) based on (44), (46a) and (47a) Evaluate (𝑤 𝑚 ) 𝑖𝑛+1(It+1) , (𝑦 𝑚 ) 𝑖𝑛+1(It+1) , (𝑧 𝑚 ) 𝑖𝑛+1(It+1) based on (46b)-(46d) 𝐞𝐧𝐝 𝐢𝐟 max 𝑚 |𝑠 𝑓(𝑚)𝑛+1(It+1) − 𝑠 𝑓(𝑚)𝑛+1(It) | < 𝜀 and max 𝑚,𝑖 |(𝑢 𝑚 ) 𝑖𝑛+1(It+1) − (𝑢 𝑚 ) 𝑖𝑛+1(It) | < 𝜀 Calculate (Θ 𝑚 ) 𝑖𝑛+1(It+1) , and (Κ 𝑚 ) 𝑖𝑛+1(It+1) , and (Γ 𝑚 ) 𝑖𝑛+1(It+1) based on (48), (50) Set 𝑠 𝑓 ( 𝑚 ) 𝑛 = 𝑠 𝑓 ( 𝑚 ) 𝑛+1 , ( 𝑢 𝑚 ) 𝑖𝑛 = ( 𝑢 𝑚 ) 𝑖𝑛+1 , (𝑤 𝑚 ) 𝑖𝑛 = (𝑤 𝑚 ) 𝑖𝑛+1 , (𝑦 𝑚 ) 𝑖𝑛 = (𝑦) 𝑖𝑛+1 , and (𝑧 𝑚 ) 𝑖𝑛 = (𝑧 𝑚 ) 𝑖𝑛+1(It+1) 𝐞𝐥𝐬𝐞 Set 𝑠 𝑓(𝑚) 𝑛+1(It) = 𝑠 𝑓(𝑚) 𝑛+1(It+1) , (𝑢 𝑚 ) 𝑖𝑛+1(It) = (𝑢 𝑚 ) 𝑖𝑛+1(It+1) , (𝑤 𝑚 ) 𝑖𝑛+1(It) = (𝑤 𝑚 ) 𝑖𝑛+1(It+1) , (𝑦 𝑚 ) 𝑖𝑛+1(It) = (𝑦 𝑚 ) 𝑖𝑛+1(It+1) , 𝑎𝑛𝑑 (𝑧 𝑚 ) 𝑖𝑛+1(It) = (𝑧 𝑚 ) 𝑖𝑛+1(It+1) 𝐞𝐧𝐝 𝐞𝐧𝐝 𝐞𝐧𝐝 Computational Procedure with Newton Iterative Method
The Newton method is known to provide quadratic convergence to the solution
𝐹(𝒙) = 0, and solving our numerical scheme with this method presents an alternative and good choice. Based on (44), (46), we start our iteration in the form 2
𝐹 (𝒖 𝑚n+1(It=0) ) = 𝐴 𝑚𝑢 𝒖 𝑚n+1(It=0) − (𝒃 𝑚𝑢 ) 𝑛 ; 𝐹 (𝒘 𝑚n+1(It=0) ) = 𝐴 𝑚 𝒘 𝑚n+1(It=0) − (𝒃 𝑚𝑤 ) 𝑛 , (68𝑎) 𝐹 (𝒚 𝑚n+1(It=0) ) = 𝐴 𝑚 𝒚 𝑚n+1(It=0) − (𝒃 𝑚𝑦 ) 𝑛 ; 𝐹 (𝒛 𝑚n+1(It=0) ) = 𝐴 𝑚 𝒛 𝑚n+1(It=0) − (𝒃 𝑚𝑧 ) 𝑛 . (68𝑏) Matrix 𝐴 𝑚𝑢 is symmetric, sparse and tridiagonal with constant coefficients likewise 𝐴 𝑚 . The former differs from the latter because of the boundary treatment in (44). Next, we evaluate 𝐽 (𝒖 𝑚n+1(It=0) ) ∆𝒖 𝑚n+1(It=0) = 𝐹 (𝒖 𝑚n+1(It=0) ) ; 𝐽 (𝒘 𝑚n+1(It=0) ) ∆𝒘 𝑚n+1(It=0) = 𝐹 (𝒘 𝑚n+1(It=0) ), (69𝑎)
𝐽 (𝒚 𝑚n+1(It=0) ) ∆𝒚 𝑚n+1(It=0) = 𝐹 (𝒚 𝑚n+1(It=0) ) ; 𝐽 (𝒛 𝑚n+1(It=0) ) ∆𝒛 𝑚n+1(It=0) = 𝐹 (𝒛 𝑚n+1(It=0) ). (69𝑏)
The advantage of our model is that the generated discrete Jacobian matrix is symmetric, sparse and tridiagonal with constant coefficients as one can easily observe from (44), (46). More precisely,
𝐽(𝒖 𝑚n ) ≡ 𝐴 𝑚𝑢 ; 𝐽(𝒘 𝑚n ) = 𝐽 (𝒚 𝑚n) ) = 𝐽(𝒛 𝑚n ) ≡ 𝐴 𝑚 for all 𝑛 where 𝑚 = 1,2, … , 𝐼. (70) It presents some nice properties that reduce the cost of computing the Jacobian matrix and enable the use of the Thomas Algorithm for solving ∆𝒖 𝑚n+1(It=0) , ∆𝒖 𝑚n+1(It=0) , ∆𝒖 𝑚n+1(It=0) , and ∆𝒖 𝑚n+1(It=0) . The next iteration is obtained as follows: 𝒖 𝑚n+1(It=1) = 𝐺 (𝒖 𝑚n+1(It=0) ) = 𝒖 𝑚n+1(It=0) − ∆𝒖 𝑚n+1(It=0) ; (71𝑎) 𝑠 𝑓(𝑚)𝑛+1(It=1) = 𝐾 − ( 𝑢 𝑚n+1 ( It=1 ) ) 𝑖=0 ; (71𝑏) 𝒘 𝑚n+1(It=1) = 𝐺 (𝒘 𝑚n+1(It=0) ) = 𝒘 𝑚n+1(It=0) − ∆𝒘 𝑚n+1(It=0) ; (71𝑐) 𝒚 𝑚n+1(It=1) = 𝐺 (𝒚 𝑚n+1(It=0) ) = 𝒚 𝑚n+1(It=0) − ∆𝒚 𝑚n+1(It=0) ; (71𝑑) 𝒛 𝑚n+1(It=1) = 𝐺 (𝒛 𝑚n+1(It=0) ) = 𝒛 𝑚n+1(It=0) − ∆𝒛 𝑚n+1(It=0) , for 𝑚 = 1,2, … , 𝐼. (71𝑒) The iterative process continues until the convergence criterion of both max 𝑚 |𝒔 𝑓(𝑚)𝑛+1(It+1) − 𝒔 𝑓(𝑚)𝑛+1(It) | < 𝜀 and max 𝑚,𝑖 |( 𝒖 𝑚 ) 𝑖𝑛+1 ( It+1 ) − ( 𝒖 𝑚 ) 𝑖𝑛+1 ( It ) | < 𝜀 is satisfied . It should be pointed out that to facilitate computation using the Newton method, we adopt the procedure used in the work of Egorova et al. (2016) by treating the coupled regime in the set of the system of PDEs explicitly. Moreover, in the work of Khaliq and Liu (2009), the linear implicit approach was adopted in treating the coupled regime. An algorithm for obtaining the numerical solutions of the optimal exercise boundary, asset option and the option Greeks in each regime using the Newton method is described below. Algorithm 2.
An algorithm based on the Newton iteration Initialize 𝑠 𝑓(𝑚)𝑛 , 𝒖 𝑚n , 𝒘 𝑚n , 𝒚 𝑚n , 𝒛 𝑚n , 𝚯 𝑚n , 𝚱 𝑚n , and 𝚪 𝑚n for 𝑚 = 1,2, … , 𝐼 𝐟𝐨𝐫 𝐧 = 𝟏 𝐭𝐨 𝐍 Compute 𝒖 𝑙n , 𝒘 𝑙n , 𝒚 𝑙n , 𝒛 𝑙n for 𝑙 = 1,2, … , 𝐼 and 𝑙 ≠ 𝑚 based on (51)-(54) Set 𝑠 𝑓(𝑚)𝑛+1(It=0) = 𝑠 𝑓(𝑚)𝑛 , 𝒖 𝑚n+1(It=0) = 𝒖 𝑚n , 𝒘 𝑚n+1(It=0) = 𝒘 𝑚n , 𝒚 𝑚n+1(It=0) = 𝒚 𝑚n , 𝒛 𝑚n+1(It=0) = 𝒛 𝑚n 𝐰𝐡𝐢𝐥𝐞 𝐭𝐫𝐮𝐞 𝐟𝐨𝐫 𝐦 = 𝟏 𝐭𝐨 𝐈 Compute 𝐹 ( 𝒖 𝑚n+1 ( It=0 ) ) . Obtain ∆𝒖 𝑚n+1(It=0) using the Thomas Algorithm with 𝐽(𝒖 𝑚n ) ≡ 𝐴 𝑚 𝑢 based on (70). Compute 𝒖 𝑚n+1(It=1) based on (71) and evaluate 𝑠 𝑓(𝑚)𝑛+1(It=1) from (𝑢 𝑚n+1(It=1) ) 𝑖=0 based on (47a) 𝐞𝐧𝐝 𝐢𝐟 max 𝑚 |𝑠 𝑓(𝑚)𝑛+1(It+1) − 𝑠 𝑓(𝑚)𝑛+1(It) | < 𝜀 and max 𝑚,𝑖 | 𝒖 𝑚n+1 ( It+1 ) − 𝒖 𝑚n+1 ( It ) | < 𝜀 Calculate 𝚯 𝑚n+1 based on (48)-(50) Set 𝑠 𝑓(𝑚)𝑛 = 𝑠 𝑓(𝑚)𝑛+1 , 𝒖 𝑚n = 𝒖 𝑚n+1 𝐞𝐥𝐬𝐞 Set 𝑠 𝑓(𝑚) 𝑛+1(It) = 𝑠 𝑓(𝑚) 𝑛+1(It+1) , 𝒖 𝑚n+1(It) = 𝒖 𝑚n+1(It+1) 𝐞𝐧𝐝 𝐞𝐧𝐝 𝐰𝐡𝐢𝐥𝐞 𝐭𝐫𝐮𝐞 Compute 𝒘 𝑚n+1(It+1) in the same manner as 𝒖 𝑚n+1(It+1) based on (46), (70), (71) Calculate 𝚱 𝑚n+1 based on (48)-(50) 𝐞𝐧𝐝 𝐰𝐡𝐢𝐥𝐞 𝐭𝐫𝐮𝐞 Compute 𝒚 𝑚n+1(It+1) in the same manner as 𝒖 𝑚n+1(It+1) based on (46), (70), (71) Calculate 𝚪 𝑚n+1 based on (48)-(50) 𝐞𝐧𝐝 𝐰𝐡𝐢𝐥𝐞 𝐭𝐫𝐮𝐞 Compute 𝒛 𝑚n+1(It+1) in the same manner as 𝒖 𝑚n+1(It+1) based on (46), (70), (71) 𝐞𝐧𝐝 𝐞𝐧𝐝 Numerical Experiments
To test the accuracy and applicability of the present scheme, we consider the American put options pricing problems with two regimes, four regimes, eight regimes, and sixteen regimes, respectively. The numerical code was written with MATLAB 2019a on Intel Core i5-3317U CPU 1.70GHz 64-bit ASUS Laptop. The numerical procedures were carried out on the mesh with a uniform grid size. 4
Numerical Examples: Two Regimes
We first consider the American put options with two regimes. We label our present method as “FF-CS1, FF-CS2, FF-CS3 and FF-CS4 which we denote as the front fixing-compact scheme with cubic Hermite interpolation and GS iteration, with quintic Hermite interpolation and GS iteration, with cubic Hermite interpolation and Newton Iteration, and with quintic Hermite interpolation and Newton Iteration, respectively. We further compare them with MTree (Liu, 2010), IMS1, IMS2 (Khaliq and Liu, 2009), MOL (Chiarella et al., 2016), RBF-FD (Li et al., 2018), FF-expl (Egorova et al., 2016), ETD-CN (Khaliq et al., 2013), Iterated Optimal Stopping and Local Optimal Iteration (Babbin et al., 2011) as listed in Tables 1-4, and 7. The option Greeks results were also listed in Tables 5 and 6.
Example 1:
We consider a switching regime problem with the strike price chosen to be
𝐾 = 9 at the expiration time
𝑇 = 1.
In our computation, we chose the interval 𝑚 ≤ 3 with the grid size ℎ = 0.1,0.05, and for FF-CS1 and FF-CS2 and ℎ = 0.01 for FF-CS3 and FF-CS4. The time step 𝑘 was determined using 𝑘 = ℎ . The parameters were given as 𝑄 = [−6 69 − 9] , 𝒓 = [0.100.05] , 𝝈 = [0.800.30] , 𝜀 = 10 −8 (72) Figs. 2 and 3 show the profiles of the option prices, Greek parameters, and optimal exercise boundaries for the two-regime case. From Tables 1-3, one can easily observe that the data obtained based on FF-CS1 and FF-CS2, FF-CS3 and FF-CS4, respectively, when ℎ = 0.01 are the same as those obtained from MOL, MTree and RBF-FD up to 5 digits in most cases. In particular, FF-CS3 and FF-CS4 are more than five times faster than FF-CS1 and FF-CS2. Moreover, Chiarella et al. (2016) pointed out that data obtained from MTree data was used as the benchmark in the work of Khaliq and Liu (2009). Generally, our data slightly decreases in direct proportion with ℎ . Example 2:
We investigate the performance of our method as compared with MOL when there is no jump between regimes (Chiarella et al., 2016; Meyer and van der Hoek, 1997 ) . We use the same data provided in the first example. The grid size and generator matrix were chosen to be ℎ = 0.01 and 𝑄 = [0 00 0], (73) respectively. The obtained result was listed in Table 8. It can be seen from Table 8 that the data obtained based on FF-CS1 and FF-CS2 are virtually the same. This is because there was no jump between different states. Fig. 2.
Option Greeks for the two-regime case when 𝜏 = 𝑇 (example 1). Fig. 3.
Asset options and optimal exercise boundaries for the two-regime case when 𝜏 = 𝑇 (example 1).
Example 3 : We compare our results with Iterated Optimal Stopping and Local Optimal Iteration (Babbin et al., 2011). The strike price was chosen to be
𝐾 = 10 at the expiration time
𝑇 = 1.
The grid size was chosen to be ℎ = 0.01.
The parameters were given as
𝑄 = [ −3 32 − 2] , 𝒓 = [0.050.05] , 𝝈 = [0.30.4]. (74)
Like other given examples, the convergence criterion 𝜀 = 10 −8 was chosen and the time step 𝑘 was determined using 𝑘 = ℎ . In Table 7, we compared the result with RBF-FD, IOS, and LOP methods. The gamma and speed plots from this example were shown in Fig. 4. To check the accuracy of our present method, we calculated the convergence rate from the asset option in regime 1. To obtain the convergence rate of our numerical scheme, we defined the maximum error using the notation
Table 1.
Comparison of American put option price in regime 1 for example 1. S MTree IMS1 IMS2 MOL FF-CS1 FF-CS2 FF-CS3 FF-CS4 ℎ = 0.1 0.05 0.01 0.1 0.05 0.01 0.01 0.01
Table 2.
Comparison of American put option price in regime 2 for example 1. S MTree IMS1 IMS2 MOL FF-CS1 FF-CS2 FF-CS3 FF-CS4 ℎ = 0.1 0.05 0.01 0.1 0.05 0.01 0.01 0.01
Table 3.
Further
Comparison of American put option price for example 1. S RBF-FD ETD-CN FF-expl FF-CS2 (ℎ = 0.01)
Regime 1 Regime 2 Regime 1 Regime 2 Regime 1 Regime 2 Regime 1 Regime 2 9.0 1.9718 1.8825 1.9756 1.8859 1.9713 1.8817 1.9720 1.8825 10.5 1.5185 1.4274 1.5213 1.4301 1.5177 1.4265 1.5185 1.4273 12.0 1.1803 1.0924 1.1825 1.0945 1.1796 1.0915 1.1803 1.0923
Table 4.
Comparing FF-CS3, and FF-CS4 up to sixteen digits at strike price for example 1.
Strike price FF-CS3 FF-CS4 Regime 1 Regime 2 Regime 1 Regime 2 9.0 1.971738757801249 1.882203676543793 1.971733636374602 1.882198321043946
Table 5.
American put option Greeks for the two-regime case for example 1 with FF-CS2.
Delta Gamma Speed S Regime 1 Regime 2 Regime 1 Regime 2 Regime 1 Regime 2 3.5 -1.0000 -1.0000 0.0000 0.0000 0.0000 0.0000 4.0 -0.9652 -1.0000 0.0164 0.0000 0.0171 0.0000 4.5 -0.8749 -0.9171 0.0508 0.0497 0.0438 0.0470 6.0 -0.6426 -0.6571 0.0851 0.0905 0.0381 0.0462 9.5 -0.3165 -0.3181 0.0560 0.0594 0.0015 0.0013 12.0 -0.1945 -0.1913 0.0347 0.0361 -0.0025 -0.0031
Table 6.
American put option Greeks for the two-regime case for example 1 with FF-CS2.
Theta Delta-Decay Color S Regime 1 Regime 2 Regime 1 Regime 2 Regime 1 Regime 2 3.5 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 4.0 -0.0300 0.0000 -0.0211 0.0000 -0.0086 0.0000 4.5 -0.1200 -0.0850 -0.0690 -0.0722 -0.0229 -0.0295 6.0 -0.4083 -0.4279 -0.1160 -0.1358 -0.0125 -0.0206 9.5 -0.7904 -0.8467 -0.0310 -0.0317 0.0108 0.0132 12.0 0.8248 -0.8700 0.0169 0.0240 0.0061 0.0069 8
Table 7.
Comparison of American put option price with IOS and LOI in Regime 1 for example 3. S IOS LOI RBF-FD FF-CS2 FF-CS4 Maximum Refinement ℎ = 0.01 . Table 8.
Comparison of American put option price with no jump between regimes for example 2. S MOL FF-CS1 FF-CS2 Regime 1 Regime 2 Regime 1 Regime 2 Regime 1 Regime 2 6.00 3.666762424 3.000000000 3.666746420 3.000000000 3.666746420 3.000000000 9.00 2.375385605 0.888311178 2.375408073 0.888393716 2.375408073 0.888393716 12.00 1.604853957 0.203543056 1.604912489 0.204583983 1.604912489 0.203637945
Table 9.
The maximum errors and convergence rates for regime 1 for example 1. ℎ maximum error convergence rate (FF-CS1) (FF-CS2) (FF-CS1) (FF-CS2) −1 −1 −2 −2 −2 −3 −3 −2 −4 −4 Table 10.
Average CPU time(s) per each time step for the two-regime example. ℎ CPU Time(s) FF-CS1 FF-CS2 FF-CS3 FF-CS4 0.1 0.182 0.209 0.058 0.061 0.05 0.311 0.316 0.084 0.078 0.01 3.697 4.167 0.768 0.771
Fig. 4.
Gamma and speed options for the two-regime case when 𝜏 = 𝑇 (example 3). 9
𝐸(ℎ, 𝑘) = max |(𝑢 ) 𝑖𝑛 (ℎ, 𝑘) − (𝑢 ) 𝑖𝑛 (ℎ/2, 𝑘/4)|, (75𝑎) 𝐸(ℎ/2, 𝑘/4) = max |(𝑢 ) 𝑖𝑛 (ℎ/2, 𝑘/4) − (𝑢 ) 𝑖𝑛 (ℎ/4, 𝑘/16)|, (75𝑏) where 𝑘 = ℎ . (𝑢 ) 𝑖𝑛 (ℎ, 𝑘), (𝑢 ) 𝑖𝑛 (ℎ/2, 𝑘/4), and (𝑢 ) 𝑖𝑛 (ℎ/4, 𝑘/16) are the numerical solutions from regime 1 obtained based on ℎ and 𝑘 , ℎ/2 and 𝑘/4 , and ℎ/4 and 𝑘/16 , respectively. As such, the convergence rate was evaluated using the following equation as:
𝑅𝑎𝑡𝑒 = log 𝐸(ℎ, 𝑘)𝐸(ℎ/2, 𝑘/4) . (76)
Table 9 lists the maximum errors and convergence rates of FF-CS1 and FF-CS2 obtained based on ℎ =0.2, 0.1, 0.05, 0.025, 0.0125.
It can be seen from Table 9 that the convergence rate is above 3.0, indicating that our present method provides a more accurate solution than the existing methods do.
Besides , t he computational speed of FF-CS1, FF-CS2, FF-CS3, and FF-CS4 is very fast as seen from Table 10.
Numerical Examples beyond Two Regimes
Commonly, previous works of literature have limited the regime-switching analysis to two and four regimes. Moreover, Chiarella et al. (2016) and Khaliq and Liu (2008) pointed out that the method proposed by Buffington and Elliot (2002) cannot be extended beyond two regimes. To show that our method can compute a large finite state space, we wrote a sequence of MATLAB function files and used it to write a few lines of code that can take any number of finite state spaces. We then considered the American put options pricing problems with four, eight, and sixteen regimes, respectively. The strike price and expiration time were chosen to be
𝐾 = 9 and
𝑇 = 1, respectively.
In our computation, we chose the interval 𝑚 ≤ 3 where the grid size ℎ = 10 −2 and 𝑘 = 10 −4 . The four-regime example was computed with the convergent criterion of 𝜀 = 10 −8 while eight- and sixteen-regime examples were computed with 𝜀 = 10 −7 and 𝜀 = 10 −8 for FF-CS2 and FF-CS4, respectively. The parameters are given in (77)-(79), respectively. Four-regime example : 𝑄 = [−1 1/3 1/3 1/31/3 − 1 1/3 1/31/3 1/3 − 1 1/31/3 1/3 1/3 − 1 ] , 𝒓 = [0.020.100.060.15] , 𝝈 = [0.900.500.700.20]. (77) Eight-regime example : 𝑄 = [ −1 0.2 0.2 0.2 0.1 0.1 0.1 0.10.2 − 1 0.1 0.1 0.1 0.2 0.2 0.10.2 0.1 − 1 0.1 0.2 0.1 0.1 0.20.2 0.1 0.2 − 1 0.2 0.1 0.1 0.10.1 0.2 0.1 0.1 − 1 0.2 0.1 0.20.2 0.2 0.2 0.1 0.1 − 1 0.1 0.10.1 0.1 0.2 0.2 0.2 0.1 − 1 0.10.1 0.1 0.1 0.2 0.1 0.2 0.2 − 1] , 𝒓 = [ 0.030.150.200.090.050.120.150.18] , 𝝈 = [ 0.800.400.500.700.450.380.300.25] . (78)
Sixteen-regime example: 𝒓 = [ 0.04 0.15 0.03 0.30 0.13 0.12 0.10 0.18 0.08 0.25 0.06 0.20 0.21 0.07 0.12 0.19], 𝝈 = [ 0.07 0.30 0.90 0.80 0.25 0.15 0.12 0.28 0.85 0.35 0.39 0.72 0.45 0.18 0.20 0.25],
𝑄 = [ −3 0.2 0.2 0.2 0.2 ⋯ 0.2 0.2 0.2 0.2 0.20.2 − 3 0.2 0.2 0.2 ⋯ 0.2 0.2 0.2 0.2 0.20.2 0.2 − 3 0.2 0.2 ⋯ 0.2 0.2 0.2 0.2 0.20.2 0.2 0.2 − 3 0.2 ⋯ 0.2 0.2 0.2 0.2 0.20.2 0.2 0.2 0.2 − 3 ⋯ 0.2 0.2 0.2 0.2 0.2⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱ ⋱0.2 0.2 0.2 0.2 0.2 ⋯ 0.2 − 3 0.2 0.2 0.20.2 0.2 0.2 0.2 0.2 ⋯ 0.2 0.2 − 3 0.2 0.20.2 0.2 0.2 0.2 0.2 ⋯ 0.2 0.2 0.2 − 3 0.20.2 0.2 0.2 0.2 0.2 ⋯ 0.2 0.2 0.2 0.2 − 3 ] . (79)
Figs. 5-8 plot the profiles of the option prices, Greek parameters, and optimal exercise boundaries for the four, eight, and sixteen regimes. Tables 11-14 list the option prices and Greeks of the four, eight and sixteen regimes using the asset values in the interval of . Fig. 5.
Asset options and optimal exercise boundaries for the four-regime case when 𝜏 = 𝑇 . Fig. 6.
Option Greeks for the four-regime case when 𝜏 = 𝑇 . Fig. 7.
Optimal exercise boundaries, asset options, and option Greeks for the eight-regime case when 𝜏 = 𝑇 . Fig. 8.
Optimal exercise boundaries, asset options, and option Greeks for the sixteen-regime case when 𝜏 = 𝑇 . Table 11.
Comparison of American put options price for the four-regime example.
MTree RBF-FD FF-expl S Reg 1 Reg 2 Reg 3 Reg 4 Reg 1 Reg 2 Reg 3 Reg 4 Reg 1 Reg 2 Reg 3 Reg 4 7.5 3.1433 2.2319 2.6746 1.6574 3.1424 2.2320 2.6744 1.6576 3.1421 2.2313 2.6739 1.6573 9.0 2.5576 1.5834 2.0568 0.9855 2.5564 1.5835 2.0566 0.9857 2.5563 1.5827 2.0559 0.9850 10.5 2.1064 1.1417 1.6014 0.6533 2.1052 1.1415 1.6013 0.6554 2.1047 1.1406 1.6004 0.6546 12.0 1.7545 0.8377 1.2625 0.4708 1.7527 0.8377 1.2625 0.4708 1.7524 0.8368 1.2614 0.4700 ETD-CN FF-CS2 7.5 3.1513 2.2384 2.6813 1.6664 3.1418 2.2319 2.6746 1.6578 9.0 2.5641 1.5884 2.0623 0.9903 2.5545 1.5835 2.0567 0.9858 10.5 2.1113 1.1451 1.6057 0.6580 2.1015 1.1414 1.6012 0.6553 12.0 1.7578 0.8377 1.2658 0.4725 1.7525 0.8374 1.2621 0.4706 5
Table 12.
American put options price for the eight- and sixteen-regime examples using FF-CS2.
Eight regimes Sixteen regimes S Reg 1 Reg 2 Reg 4 Reg 6 Reg 8 Reg 1 Reg 2 Reg 4 Reg 6 Reg 8 Reg 12 Reg 16 3.5 5.5551 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 4.0 5.1238 5.0000 5.0006 5.0000 5.0000 5.0074 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 4.5 4.7190 4.5000 4.5319 4.5000 4.5000 4.5385 4.5000 4.5000 4.5000 4.5000 4.5000 4.5000
Table 13.
American put options price for the eight- and sixteen-regime examples using FF-CS4.
Eight regimes Sixteen regimes S Reg 1 Reg 2 Reg 4 Reg 6 Reg 8 Reg 1 Reg 2 Reg 4 Reg 6 Reg 8 Reg 12 Reg 16 3.5 5.5555 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 5.5000 4.0 5.1244 5.0000 5.0006 5.0000 5.0000 5.0075 5.0000 5.0000 5.0000 5.0000 5.0000 5.0000 4.5 4.7197 4.5000 4.5320 4.5000 4.5000 4.5387 4.5000 4.5000 4.5000 4.5000 4.5000 4.5000
Table 14.
American put option Greeks for the four-regime example with FF-CS2.
Delta Gamma Speed S Reg 1 Reg 2 Reg 4 Reg 1 Reg 2 Reg 4 Reg 1 Reg 2 Reg 4 3.5 -0.8246 -1.0000 -1.0000 0.0515 0.0000 0.0000 0.0949 0.0000 0.0000 6.0 -0.5739 -0.7442 -1.0000 0.0638 0.0723 0.0000 0.0204 0.0262 0.0000 9.0 -0.3401 -0.3546 -0.3026 0.0488 0.0727 0.1086 0.0005 0.0011 -0.0027 12.0 -0.2055 -0.1679 -0.0958 0.0289 0.0370 0.0269 -0.0023 -0.0048 -0.0081 At the money option, volatility has a negligible impact on the delta option for all the regimes. Hence, the plot for each regime intersects at the strike price. For long put options, as we move deep in the money and out of the money, delta converges to -1 and 0, respectively.
Gamma is maximum when at the money. Ignoring the sign convention, the theta of ATM is maximum. Delta decay and color options measure the rate at which delta and gamma options decay, respectively. 6 Conclusion
We have developed an accurate numerical method for solving American put options with regime-switching. Through the front-fixing transformation, we were able to map the optimal exercise boundary for each regime to a fixed interval. The derivative transformation enables us to employ the higher-order compact finite difference method coupled with the Hermite interpolation for solving the system of the asset, delta, gamma, and speed options while capturing the optimal exercise boundary and theta, delta decay and color options. Moreover, our method has a substantial advantage because it simultaneously and accurately calculates asset, delta, gamma, speed, theta, delta decay, and color options, as well as optimal exercise boundary during iteration. Greek parameters are difficult to estimate correctly as can be seen from previous works of literature. However, by formulating a set of systems of PDEs that consist of the asset option and its derivatives for each regime, we were able to estimate those parameters with higher-order accuracy. Our numerical discretization also presents a system where the coefficient matrix is tridiagonal and positive definite with constant-elements, which enables us to implement both Gauss-Seidel and Newton iterations (with Thomas algorithm) with simple computation. The present scheme has been tested in two-, four-, eight-, and sixteen-regime problems and with cubic and quintic Hermite interpolations. The results show that the method provides an accurate solution and is fast in computation as compared with the existing methods. Future research will include applying this method to non-constant volatility and/or interest rate cases.
References Adam, Y. (1975). A Hermitian finite difference method for the solution of parabolic equation. Computers and Mathematics with Applications, 1, 393-406. 2.
Adam, Y. (1976). Highly accurate compact implicit methods and boundary condition. Journal of Computational Physics, 24, 10-22. 3.
Babbin, J., Forsyth, P. A., and Labahn, G. (2014). A comparison of iterated optimal stopping and local policy iteration for American options under regime switching. Journal of Scientific Computing, 58, 409-430. 4.
Bidégaray-Fesquet, B. (2006). Von-Neumann stability analysis of FD-TD methods in complex media. 5.
Bingham, N. H. (2006). Financial modeling with jump processes. Journal of the American Statistical Association, 101 (475), 1315-1316 (https://doi:10.1198/jasa.2006.s130). 7 6.
Blackwell, B. F., and Hogan, R. E. (1994). One-dimensional ablation using Landau transformation and finite control volume procedure. Journal of Thermophysics and Heat Transfer, 8 (2), 282-287 (https://doi:10.2514/3.535). 7.
Boyarchenko, S. I., and Levendorskii, S. Z. (2008). Pricing American options in regime-switching models: FFT Realization. SSRN Electronic Journal (https://doi:10.2139/ssrn.1127562). 8.
Bravo, E. G., and Mcgraw, T. (2015). Visualizing Aldo Giorgini’s ideal’s flow. Springer International Publishing. Doi:10.1007/978-3-319-27863-6_72. 9.
Buffington, J., and Elliott, R. J. (2002). American options with regime-switching. International Journal of Theoretical and Applied Finance, 5 (05), 497-514 (https://doi:10.1142/s0219024902001523). 10.
Burden, R. L., Faires, D. J., and Burden, A. M. (2016). Numerical Analysis. Cengage Learning, Boston, MA. 11.
Cao, H. H., Liu, L. B., Zhang, Y., and Fu, S. M. (2011). A fourth-order method of the convection-diffusion equation with Neumann boundary conditions. Applied Mathematics and Computation, 217, 9133-9141. 12.
Chapra, S. C. (2012). Applied Numerical Methods with MATLAB for Engineers and Scientist. MC Graw-Hill, New York. 13.
Chen, Y., Xiao, A., and Wang, W. (2019). An IMEX-BDF2 compact scheme for pricing options under regime-switching jump-diffusion models. Mathematical Methods in the Applied Sciences, (https://doi:10.1002/mma.5539). 14.
Chiarella, C., Nikitopoulos Sklibosios, C., Schlogl, E., and Yang, H. (2016). Pricing American options under regime switching using method of lines. SSRN Electronic Journal, (https://doi:10.2139/ssrn.2731087). 15.
Chockalingam, A., and Muthuraman, K. (2011). American options under stochastic volatility. Operations Research, 59(4), 793-809 (https://doi:10.1287/opre.1110.0945). 16.
Company, R., Egorova, V.N., and Jódar, L. (2016). A positive, stable and consistent front-fixing numerical scheme for American Options. Springer International Publishing AG 2016 (https://doi 10.1007/978-3-319-23413-7_10). 17.
Company, R., Egorova, V., Jódar, L., and Vázquez, C. (2016). Computing American option price under regime switching with rationality parameter. Computers & Mathematics with Applications, 72 (3), 741-754 (https://doi:10.1016/j.camwa.2016.05.026). 8 18.
Company, R., Egorova, V. N., and Jódar, L. (2016). Constructing positive reliable numerical solution for American call options: A new front-fixing approach. Journal of Computational and Applied Mathematics, 291, 422-431 (https://doi:10.1016/j.cam.2014.09.013). 19.
Company, R., Egorova, V., Jódar, L., and Vázquez, C. (2016). Finite difference methods for pricing American put option with rationality parameter: Numerical analysis and computing. Journal of Computational and Applied Mathematics, 304, 1-17 (https://doi:10.1016/j.cam.2016.03.001). 20.
Cont., R., and Tankov, P. (2004). Financial Modelling with Jump Diffusion. Chapman & Hall, London. 21.
Crank, J. (1984). Free and moving boundary problem, Oxford University Press, London. 22.
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. 23.
Düring, B., and Fournié, M. (2012). High-order compact finite difference scheme for option pricing in stochastic volatility models. Journal of Computational and Applied Mathematics, 236 (17), 4462-4473 (https://doi:10.1016/j.cam.2012.04.017). 24.
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 (1), 224–237 (https://doi:10.1016/j.camwa.2015.11.019). 25.
Elliott, R. J., Siu, T. K., Chan, L., and Lau, J. W. (2007). Pricing options under a generalized markov-modulated jump-diffusion model. Stochastic Analysis and Applications, 25 (4), 821–843.doi:10.1080/07362990701420118 . 26.
Gao, G. and Sun, Z. (2013). Compact difference schemes for heat equation with Neumann boundary condition (II), Numerical Methods for Partial Differential Equations, 29 (5), 1459-1486. 27.
Garnier, J., and Sølna, K. (2017). Correction to Black--Scholes formula due to fractional stochastic volatility. SIAM Journal on Financial Mathematics, 8 (1), 560–588 (https://doi:10.1137/ 15m1036749 Stochastic Volatility Correction to Black-Scholes). 28.
Guoqing, Y., and Hanson, F. B. (2006). Option pricing for a stochastic-volatility jump-diffusion model with log-uniform jump-amplitudes. American Control Conference, (https://doi:10.1109/ acc.2006.1657175). 29.
Hamilton, J. D. (1989). A new approach to the economic analysis of nonstationary time series and the business cycle. Econometrica, 57 (2), 357-384 (https://doi:10.2307/1912559). 9 30.
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, 2016, 1-14 (https://doi:10.1155/2016/2474305). 31.
Hirsch, C. (2001). Numerical Computation of Internal and External Flows: Fundamental of Numerical Discretization. Wiley, New York. 32.
Hirsh, R. S. (1975). Higher order accurate difference solutions of fluid mechanics by a compact differencing technique. Journal of Computational Physics, 19, 90-109. 33.
Huang, Y., Forsyth, P. A., and Labahn, G. (2011). Methods for pricing american options under regime switching. SIAM Journal on Scientific Computing, 33 (5), 2144-2168 (https:// doi:10.1137/110820920). 34.
Hull, J., and White, A. (1987). The pricing of options on assets with stochastic volatilities. The Journal of Finance, 42 (2), 281-300. 35.
Ikonen, S., and Toivanen, J. (2007). Efficient numerical methods for pricing American options under stochastic volatility. Numerical Methods for Partial Differential Equations, 24 (1), 104-126 (https://doi:10.1002/num.20239). 36.
Kangro, R., and Nicolaides, R. (2000). Far field boundary conditions for Black--Scholes equations. SIAM Journal on Numerical Analysis, 38 (4), 1357-1368 (https://doi:10.1137/ s0036142999355921). 37.
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 (03), 319-340 (https://doi:10.1142/s0219024909005245). 38.
Khaliq A. Q. M., Kleefeld, B., and Liu, R. H. (2013). Solving complex PDE systems for pricing American options with regime-switching by efficient exponential time differencing schemes. Numerical Methods Partial Differential Equation, 29 (1), 320-336. 39.
Kou, S. G. (2002). A jump-diffusion model for option pricing. Management Science, 48 (8), 1086-1101 (https://doi:10.1287/mnsc.48.8.1086.166). 40.
Kwok Y. K. (2008). Mathematical Models of Financial Derivatives. Springer, Berlin. 41.
Landau, H. G. (1950). Heat conduction in a melting solid. Quarterly of Applied Mathematics, 8 (1), 81-94 (https://doi:10.1090/qam/33441). 42.
Li, H., Mollapourasi, R., and Haghi, M. (2018). A local radial basis function method for pricing options under the regime switching model. Journal of Scientific Computing, 79, 517-541. 43.
Liao, W., Zhu, J., and Khaliq, A. Q. M. (2002). An efficient high-order algorithm for solving systems of reaction-diffusion equations. Numerical Methods for Partial Differential Equations, 18 (3), 340-354. 0 44.
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 (6), 1009-1023 (https://doi:10.1080/00207160802609829). 45.
Liu, R. H. (2010). Regime-switching recombining tree for option pricing. International Journal of Theoretical and Applied Finance, 13, 479-499. 46.
Liu, R. H., Zhang, Q., and Yin, G. (2006). Option pricing in a regime-switching model using the fast Fourier transform. Journal of Applied Mathematics and Stochastic Analysis, 2006, 1-22 (https://doi:10.1155/jamsa/2006/18109). 47.
Mamon, R. S., and Rodrigo, M. R. (2005). Explicit solutions to European options in a regime-switching economy. Operations Research Letters, 33 (6), 581-586 (https://doi:10.1016/ j.orl.2004.12.003). 48.
Meyer, G. H., and van der Hoek, J. (1997). The evaluation of American options with the method of lines. Advances in Futures and Options Research, 9 (9), 265-285. 49.
Mitchell, S. L., and Vynnycky, M. (2009). Finite-difference methods with increased accuracy and correct initialization for one-dimensional Stefan problems. Applied Mathematics and Computation, 215 (4), 1609-1621 (https://doi:10.1016/j.amc.2009.07.054). 50.
Mitchell, S. L., and Vynnycky, M. (2012). An accurate finite-difference method for ablation-type Stefan problems. Journal of Computational and Applied Mathematics, 236 (17), 4181-4192 (https://doi:10.1016/j.cam.2012.05.011). 51.
Morton, W.K., and Mayer, D. (2005). Numerical Solutions of Partial Differential Equations. Cambridge University Press, London. 52.
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 (4), 69-97. 53.
Norris, J. R. (1998). Markov Chain. Cambridge University Press, London. 54.
Shang, Q., and Bryne, B., (2019). An efficient lattice search algorithm for the optimal exercise boundary in American options. SSRN Electronic Journal. 55.
Tauryawati, M. L., Imron, C., and Putri, E. R. (2018). Finite volume method for pricing European call option with regime-switching volatility. Journal of Physics: Conference Series, 974, (https://doi:10.1088/1742-6596/974/1/012024). 56.
Toivanen, J. (2010). Finite difference methods for early exercise options. Encyclopedia of Quantitative Finance, (https://doi:10.1002/9780470061602.eqf12002). 57.
Ugur, O. (2009). Introduction to Computational Finance. Imperial College Press, London. (https://doi.org/10.1142/p556). 1 58.
Wu, L., and Kwok, Y.K. (1997). A front-fixing method for the valuation of American option, Journal of Financial Engineering, 6 (2), 83-97. 59.
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. 60.
Zhang, K., Teo, K. L., and Swartz, M. (2013). A Robust numerical scheme for pricing American options under regime switching based on penalty method. Computational Economics, 43 (4), 463-483 (https://doi:10.1007/s10614-013-9361-3). 61.