Removing the Correlation Term in Option Pricing Heston Model: Numerical Analysis and Computing

Abstract

This paper deals with the numerical solution of option pricing stochastic volatility model described by a time-dependent, two-dimensional convection-diffusion reaction equation. Firstly, the mixed spatial derivative of the partial differential equation (PDE) is removed by means of the classical technique for reduction of second-order linear partial differential equations to canonical form. An explicit difference scheme with positive coefficients and only five-point computational stencil is constructed. The boundary conditions are adapted to the boundaries of the rhomboid transformed numerical domain. Consistency of the scheme with the PDE is shown and stepsize discretization conditions in order to guarantee stability are established. Illustrative numerical examples are included.

1. Introduction

It is well recognized that the Black-Scholes model, where the prices depend only on the variance of the stock returns, leads to unreliable prices associated with the hypothesis of lognormal distribution of the asset returns and constant volatility. Since 1993, [1] it is known that volatility is a function of both the strike and the expiry date of the derivative security. Empirical evidence suggests that asset price volatility is not constant by variable and stochastic [2].

There are two prominent ways of working around this problem, namely, local volatility models [1] and stochastic volatility models [2–4].

Stochastic volatility models have the following general pattern: 𝑑𝑆(𝑡)=𝑝_𝑆(𝑆(𝑡),𝑣(𝑡),𝑡)𝑑𝑡+𝑞_𝑆(𝑆(𝑡),𝑣(𝑡),𝑡)𝑑𝑊₁,𝑑𝑣(𝑡)=𝑝_𝑣(𝑆(𝑡),𝑣(𝑡),𝑡)𝑑𝑡+𝑞_𝑣(𝑆(𝑡),𝑣(𝑡),𝑡)𝑑𝑊₂, (1) where the tradable security 𝑆(𝑡) and its variance 𝑣(𝑡) are correlated; that is, the Wiener stochastic processes satisfy 𝑑𝑊₁𝑑𝑊₂=𝜌𝑑𝑡, and the functional form of 𝑝_𝑆,𝑞_𝑆,𝑝_𝑣, and 𝑞_𝑣 are determined by the model.

The model proposed by Heston, see [3], [5, chapter 10], takes into account nonlognormal distribution of the assets returns, leverage effect, and important mean-reverting property of volatility and has a closed formula when the parameters are constant [3] or piecewise constant [6].

According to (1), the model of Heston is specified as follows: 𝑑𝑆(𝑡)=𝜇𝑆(𝑡)𝑑𝑡+𝑣(𝑡)𝑠(𝑡)𝑑𝑊₁,𝑑𝑣(𝑡)=𝜅(𝜃−𝑣(𝑡))𝑑𝑡+𝜎𝑣(𝑡)𝑑𝑊₂,𝑑𝑊₁𝑑𝑊₂=𝜌𝑑𝑡, (2) where 𝑊₁, 𝑊₂ are standard Brownian motions, 𝜇 represents the deterministic drift, 𝜅 is the mean reversion rate, 𝜃 is the long-run variance, 𝜎 is the volatility of the variance, and 𝜌 is the correlation parameter.

Applying the Itô lemma and standard arbitrage arguments leads to the partial differential equation [3, page 329–335] for the price 𝑈=𝑈(𝑆,𝑣,𝑡) of a contingent claim: 𝜕𝑈+1𝜕𝑡2𝑣𝑆²𝜕²𝑈𝜕𝑆²𝜕+𝜌𝜎𝑣𝑆²𝑈+1𝜕𝑆𝜕𝑣2𝜎²𝑣𝜕²𝑈𝜕𝑣²+𝑟𝑆𝜕𝑈()𝜕𝑆+𝜅𝜃−𝑣𝜕𝑈𝜕𝑣−𝑟𝑈=0,0<𝑆<∞,0<𝑣<∞,0≤𝑡<𝑇. (3)(4)

A European vanilla call option with strike price 𝐸 and maturity 𝑇 satisfies (3) together with the payoff final condition 𝑈(𝑆,𝑣,𝑇)=max(𝑆−𝐸,0), (5) and the boundary conditions 𝑈(0,𝑣,𝑡)=0,𝜕𝑈()𝜕𝑆∞,𝑣,𝑡=1,𝑟𝑆𝜕𝑈()+𝜕𝑆𝑆,0,𝑡𝜅𝜃𝜕𝑈()()+𝜕𝑣𝑆,0,𝑡−𝑟𝑈𝑆,0,𝑡𝜕𝑈()𝑈()𝜕𝑡𝑆,0,𝑡=0,𝑆,∞,𝑡=𝑆, (6)(7)(8)(9) where 𝜅,𝜃 are related to 𝜅,𝜃 by means of the expression 𝜅=𝜅+𝜆;  𝜃=𝜅𝜃/(𝜅+𝜆), and 𝜆 is the market price of the volatility risk.

For constant parameters, Heston [3] uses the method of characteristic functions to derive a closed-form solution involving infinite integrals. When the parameters are piecewise constant in time, one can still derive a recursive closed formula using a PDE method [6] or a Markov argument in combination with affine models [7].

For the case where coefficients are time dependent, the authors in [8] use a small volatility of volatility expansion and Malliavin calculus techniques, to derive approximate analytical solutions. Reliable numerical methods for solving problems (3)–(9) are suitable for both situations, the general time-dependent coefficient case as well as for the constant parameter case, where the closed form proposed solutions require a further numerical treatment.

There are some approaches for the numerical treatment of stochastic volatility models such as sparse wavelet [9], spectral methods [10], and finite-difference methods [11–15].

A feature of the time dependent, two-dimensional convection-diffusion-reaction equation (3) is the presence of a mixed spatial derivative. Dealing with finite-difference methods, this fact involves the existence of negative coefficient terms into the numerical scheme and deteriorates the quality of the numerical solution; see the introduction of  [16].

Furthermore, finite difference schemes in the presence of a mixed spatial derivative produces four terms more in the numerical scheme with the corresponding additional computational cost and possible rounding accumulation error.

Both papers [13, 15] construct difference schemes involving the mixed spatial derivative with associated drawbacks. Reference [13] derives a compact finite difference scheme using a nine-point computational stencil. [15] propose three splitting schemes of the alternating direction implicit (ADI) type.

In this paper, we construct explicit finite difference schemes with positive coefficients for solving the Heston model (3)–(9) for the continuously time-dependent coefficient case of the mean reversion rate and long-run variance, after removing the mixed spatial derivative. The organization of the paper is as follows. The problems coming out from negative coefficients arising from the discretization of the diffusion term are not easily remedied [16]; thus, Section 2 addresses the removing of the mixed spatial derivative in (3) by means of the classical technique for the reduction of second-order linear partial differential equations in two variables to canonical form. We also determine the rhomboid nonrectangular numerical domain where the problem is discretized after the transformation. In Section 3, we construct an explicit difference scheme with only a five-points computational stencil. The boundary conditions are adapted to the boundaries of the numerical domain. Section 4 deals with the consistency of the scheme with (3). In Section 5, we study stepsize discretization conditions in order to guarantee the positivity and stability of the numerical scheme. In Section 6, illustrative numerical examples are included.

For a matrix 𝐴=(𝑎_𝑖𝑗)_𝑚×𝑛 in ℝ^𝑚×𝑛, we denote by ‖𝐴‖=max{|𝑎_𝑖𝑗|;1≤𝑖≤𝑚,1≤𝑗≤𝑛} and by ‖𝐴‖₁ the maximum absolute column sum norm ‖𝐴‖₁=max_{1≤𝑗≤𝑛}^𝑚_𝑖=1𝑎_𝑖𝑗; (10) see [17, page 56].

2. On the Transformation Problem and Its Motivation

Recently, the authors used space-centered forward in time explicit finite difference schemes for the computation and numerical analysis of several one-dimensional option pricing problems [18–20].

Following these ideas for the two-dimensional problem (3) one gets the following: 𝑈𝑆_𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_𝑖,𝑗;𝜕𝑈𝑆𝜕𝜏^𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛+1_𝑖,𝑗−𝑢^𝑛_𝑖,𝑗𝑘;𝜕𝑈𝑆𝜕𝑆^𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_𝑖+1,𝑗−𝑢^𝑛_{𝑖−1,𝑗}2ℎ₁;𝜕²𝑈𝜕𝑆²𝑆_𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_{𝑖−1,𝑗}−2𝑢^𝑛_𝑖,𝑗+𝑢^𝑛_𝑖+1,𝑗ℎ²₁;𝜕𝑈𝑆𝜕𝑣^𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_𝑖,𝑗+1−𝑢^𝑛_{𝑖,𝑗−1}2ℎ₂;𝜕²𝑈𝜕𝑣²𝑆_𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_{𝑖,𝑗−1}−2𝑢^𝑛_𝑖,𝑗+𝑢^𝑛_𝑖,𝑗+1ℎ²₂;𝜕²𝑈𝑆𝜕𝑆𝜕𝑣^𝑖,𝑣_𝑗,𝜏^𝑛≈𝑢^𝑛_{𝑖+1,𝑗+1}+𝑢^𝑛_{𝑖−1,𝑗−1}−𝑢^𝑛_{𝑖−1,𝑗+1}−𝑢^𝑛_{𝑖+1,𝑗−1}4ℎ₁ℎ₂, (11) where 𝜏=𝑇−𝑡,𝑘=Δ𝜏,ℎ₁=Δ𝑆,ℎ₂=Δ𝑣,𝑆_𝑖=𝑖ℎ₁,𝑣_𝑗=𝑗ℎ₂, and 𝜏^𝑛=𝑛𝑘.

By discretizing (3) one achieves the following scheme: 𝑢^𝑛+1_𝑖,𝑗=𝑘1−𝑟𝑘−ℎ²₁𝑣_𝑗𝑆²_𝑖−𝑘ℎ²₂𝜎²𝑣_𝑗𝑢^𝑛_𝑖,𝑗+𝑘2𝑆_𝑖𝑣_𝑗𝑆_𝑖ℎ²₁+𝑟ℎ₁𝑢^𝑛_𝑖+1,𝑗+𝑘2𝑆_𝑖𝑣_𝑗𝑆_𝑖ℎ²₁−𝑟ℎ₁𝑢^𝑛_{𝑖−1,𝑗}+𝑘2𝜎²𝑣_𝑗ℎ²₂+𝜅𝜃−𝑣_𝑗ℎ₂𝑢^𝑛_𝑖,𝑗+1+𝑘2𝜎²𝑣_𝑗ℎ²₂−𝜅𝜃−𝑣_𝑗ℎ₂𝑢^𝑛_{𝑖,𝑗−1}+𝑘4ℎ₁ℎ₂𝜌𝜎𝑣_𝑗𝑆_𝑖𝑢^𝑛_{𝑖+1,𝑗+1}+𝑢^𝑛_{𝑖−1,𝑗−1}−𝑢^𝑛_{𝑖−1,𝑗+1}−𝑢^𝑛_{𝑖+1,𝑗−1}. (12)

Note that if 𝜌≠0, that is, when variables 𝑆 and 𝑣 are correlated, in the last term of the right-hand side of the scheme, (12) involves two terms with negative coefficients because 𝜌∈[−1,1] and 𝜎, 𝑣_𝑗, and 𝑆_𝑖 are positive. The existence of these terms with negative coefficients does not allow the technique developed in [18–20]. This fact motivates the transformation of problem (3) into an equivalent one where the mixed spatial derivative term disappears.

Firstly, we eliminate the reaction term 𝑟𝑈 by means of the substitution 𝑉=𝑒^𝑟𝜏𝑈,𝜏=𝑇−𝑡, (13) obtaining 𝜕𝑉=1𝜕𝜏2𝑣𝑆²𝜕²𝑉𝜕𝑆²𝜕+𝜌𝜎𝑣𝑆²𝑉+1𝜕𝑆𝜕𝑣2𝜎²𝑣𝜕²𝑉𝜕𝑣²+𝑟𝑆𝜕𝑉()𝜕𝑆+𝜅𝜃−𝑣𝜕𝑉.𝜕𝑣 (14)

Now, following the classical techniques for reduction of second-order linear partial differential equations in two independent variables to canonical form, see for instance chapter 3 of [21], we proceed to classify the right-hand side of (14) by means of the sign of the discriminant Δ(𝑆,𝑣)=𝐵²−4𝐴𝐶=𝜎²𝑣²𝑆²𝜌²−1, (15) where 1𝐵=𝜌𝜎𝑣𝑆;𝐴=2𝑣𝑆²1;𝐶=2𝜎²𝑣. (16)

Under the assumption of correlated variables with −1<𝜌<1, (14) becomes of elliptic type and the suitable substitution for eliminating the mixed spatial derivative term is given by solving the following ordinary differential equation: 𝑑𝑣=𝜎𝑑𝑆𝜌+𝑖1−𝜌²𝑆=𝜎𝜌+𝑖𝜌𝑆, (17) where 𝜌=1−𝜌². (18)

Solving (17) one gets the following: 𝑣+𝑍₀=𝜎𝜌+𝑖𝜌ln𝑆, (19) where the integration constant 𝑍₀ is related to the new variables by 𝑍₀=𝑦+𝑖𝑥. (20) From (19) and (20), it follows the expression of the new variables 𝑥=𝜌𝜎ln𝑆;𝑦=𝜌𝜎ln𝑆−𝑣. (21)

By denoting 𝑊(𝑥,𝑦,𝜏)=𝑉(𝑆,𝑣,𝜏), (14) takes the following equivalent elliptic form: 𝜕𝑊=1𝜕𝜏2𝜌²𝜎²𝑣𝜕²𝑊𝜕𝑥²+𝜕²𝑊𝜕𝑦²+1𝜌𝜎𝑟−2𝑣𝜕𝑊+1𝜕𝑥𝜌𝜎𝑟−2𝑣−𝜅(𝜃−𝑣)𝜕𝑊,𝜕𝑦 (22) where 𝑣 from (21) takes the expression in terms of 𝑥 and 𝑦, given by the following: 𝜌𝑣=𝜌𝑥−𝑦. (23) It is important to remark that the previous substitution has also computational advantages because the elimination of the mixed derivative saves four terms using the same discretization scheme; see the last term of the right-hand side of (12).

Note that from (21), (4) the domain of the transformed problem (22) one gets the following: −∞<𝑥<+∞,0<𝜏≤𝑇, (24) and 𝑣 means the variance of an underlying asset that must be positive, from (23) it follows that (𝜌/𝜌)𝑥−𝑦>0; thus, 𝜌−∞<𝑦<𝜌𝜌𝑥=𝑚𝑥;𝑚=𝜌. (25) Thus, the domain of the transformed problem is as follows: 𝐷=𝑥,𝑦,𝜏;𝑥∈ℝ;𝑦<𝑚𝑥;0<𝜏≤𝑇. (26) In the following we will assume that 𝑚=(𝜌/𝜌)≠0, because when 𝑚=0, 𝜌=0 and variables 𝑆 and 𝑣 are uncorrelated and the cross-derivative term disappears in (3).

Note that after substitution of (21) the final condition (5) of problem (3) takes the following form: 𝑊𝑥,𝑦,0=max𝑒^𝑥/𝜎𝜌−𝐸,0;𝑥∈ℝ,𝑦<𝑚𝑥. (27)

The boundary conditions (6), (9) of problem (3) are transformed into the following Dirichlet conditions: 𝑊𝑊𝑥,𝑦,𝜏→0,as𝑥→−∞,𝑥,𝑦,𝜏→𝑒^{𝑟𝜏+(𝑥/𝜎}𝜌₎,as𝑚𝑥−𝑦→+∞, (28)(29) and the boundary condition (7) is transformed into the following Neumann condition: 𝜎𝜌𝜕𝑊𝜕𝑥𝑥,𝑦,𝜏+𝜎𝜌𝜕𝑊𝜕𝑦𝑥,𝑦,𝜏→𝑒^{𝑟𝜏+𝑥/𝜎}𝜌,as𝑥→+∞. (30)

Boundary condition (8) can be replaced with a more practicable homogeneous Neumann condition: 𝜕𝑈()𝜕𝑣𝑆,0,𝑡=0; (31) see, for instance, [11, 13]. This boundary condition is transformed into the condition 𝜕𝑊𝜕𝑦𝑥,𝑦,𝜏→0,as𝑚𝑥−𝑦→0⁺. (32)

Remark 1. In the next section, in order to introduce a numerical scheme, we will need to determine a bounded numerical domain.
Because of the transformation of the spatial variables (21), a rectangle [𝑆₁,𝑆₂]×[𝑣₁,𝑣₂] is transformed into the rhomboid 𝐴𝐵𝐶𝐷; see Figure 1, where the sides are described by the following: 𝐴𝐷=𝑥,𝑦∈ℝ²;𝑥=𝑎=𝜌𝜎ln𝑆₁;𝑦=𝑚𝑎−𝑣;𝑣₁≤𝑣≤𝑣₂,𝐴𝐵=𝑥,𝑦∈ℝ²;𝑎≤𝑥≤𝑏=𝜌𝜎ln𝑆₂;𝑦=𝑚𝑥−𝑣₂,𝐵𝐶=𝑥,𝑦∈ℝ²;𝑥=𝑏;𝑦=𝑚𝑏−𝑣;𝑣₁≤𝑣≤𝑣₂,𝐶𝐷=𝑥,𝑦∈ℝ²;𝑎≤𝑥≤𝑏;𝑦=𝑚𝑥−𝑣₁. (33)

Figure 1 

Rhomboid ABCD.

Remark 2. For the case |𝜌|=1, fully correlated variables 𝑆 and 𝑣, from (15), the discriminant Δ=0, and (14) become a parabolic PDE. Following the techniques for reduction to canonical form, an appropriate substitution is 𝑥=𝑆; 𝑦=𝑣−𝜌𝜎ln𝑆 and the transformed equation takes the following form: 𝜕𝑊=1𝜕𝜏2𝑣𝑥²𝜕²𝑊𝜕𝑥²+𝑟𝑥𝜕𝑊+1𝜕𝑥𝜌𝜎2𝑣−𝑟+𝜅(𝜃−𝑣)𝜕𝑊,𝜕𝑦 (34) where |𝜌|=1 and 𝑣=𝑦+𝜌𝜎ln𝑥.

3. Numerical Scheme and Boundary Numerical Domain Considerations

Dealing with numerical solutions using finite difference schemes requires the selection of our numerical domain and the transfer of the analytic boundary conditions of the problem to the boundary conditions of the numerical domain avoiding artificial behaviour of the numerical solution [22]. Unlike the classical stepsize discretizations, due to the transformation (21) and Remark 1, we use a discretization of the numerical domain where the stepsize discretizations are related by the slope 𝑚=(𝜌/𝜌)≠0 appearing in (25). Thus, we guarantee that the rhomboid boundary of our numerical domain includes meshpoints of the discretization.

From [23, 24], a suitable bound for the underlying asset variable 𝑆 is available and generally accepted. In an analogous way, considering an admissible range of the variance 𝑣, we can identify a convenient-bounded numerical domain [𝑆₁,𝑆₂]×[𝑣₁,𝑣₂] in the 𝑆𝑣 plane.

In accordance with (33), let 𝑥_𝑖∈[𝑎,𝑏], with 𝑥_𝑖=𝑎+𝑖ℎ, 0≤𝑖≤𝑁_𝑥,Δ𝑥=ℎ=(𝑏−𝑎)/𝑁_𝑥. Let 𝑁_𝑦=(𝑣₂−𝑣₁)/𝑚ℎ and 𝑦_𝑗∈[𝑚𝑎−𝑣₂,𝑚𝑏−𝑣₁] of the form 𝑦₀=𝑚𝑎−𝑣₂;𝑦_𝑗=𝑦₀+𝑗𝑚ℎ;Δ𝑦=𝑚ℎ=𝑚Δ𝑥. (35)

Once 𝑖 is fixed, the admissible values of 𝑗 so that (𝑥_𝑖,𝑦_𝑖) belongs to the rhomboid 𝐴𝐵𝐶𝐷 are given by the following: 𝑖≤𝑗≤𝑁_𝑦+𝑖;0≤𝑖≤𝑁_𝑥. (36) Note that with this discretization the boundary sides of the rhomboid are partitioned in the following way: 𝑃𝐴𝐷=𝑥₀,𝑦_𝑗;0≤𝑗≤𝑁_𝑦,𝑃𝐴𝐵=𝑥_𝑖,𝑦_𝑖;0≤𝑖≤𝑁_𝑥,𝑃𝐵𝐶=𝑥_𝑁𝑥,𝑦_𝑗;𝑁_𝑥≤𝑗≤𝑁_𝑥+𝑁_𝑦,𝑃𝐶𝐷=𝑥_𝑖,𝑦_{𝑖+𝑁𝑦};0≤𝑖≤𝑁_𝑥. (37)

Thus, the rhomboid spatial domain together with its boundary contains (𝑁_𝑥+1)(𝑁_𝑦+1) mesh points. See Figure 1.

For the time variable 𝜏, we consider the standard discretization 𝜏^𝑛=𝑘𝑛;0≤𝑛≤𝑁_𝜏𝑇;Δ𝜏=𝑘=𝑁_𝜏. (38)

Hence, the numerical domain Ω is defined by the following: Ω=𝑥_𝑖,𝑦_𝑗,𝜏^𝑛∣0≤𝑖≤𝑁_𝑥;𝑖≤𝑗≤𝑁_𝑦+𝑖;0≤𝑛≤𝑁_𝜏. (39)

Note that because of (23) the value of 𝑣 at the spatial variables (𝑥_𝑖,𝑦_𝑗) is as follows: 𝑣_𝑖,𝑗=𝑚𝑥_𝑖−𝑦_𝑗;0≤𝑖≤𝑁_𝑥;𝑖≤𝑗≤𝑁_𝑦+𝑖. (40)

Using centered differences approximations for the spatial partial derivatives of (22) and forward in time partial derivatives, the numerical scheme for the internal mesh points of the numerical domain takes the following form: 𝐹𝑤^𝑛_𝑖,𝑗=𝑤^𝑛+1_𝑖,𝑗−𝑤^𝑛_𝑖,𝑗𝑘−12𝑣_𝑖,𝑗𝜎²𝜌²Δ_𝑥^𝑛_𝑖,𝑗+Δ_𝑦^𝑛_𝑖,𝑗𝜌1−𝜎𝑟−2𝑣_𝑖,𝑗∇_𝑥^𝑛_𝑖,𝑗−1𝜎𝜌𝑟−2𝑣_𝑖,𝑗−𝜅^𝑛𝜃^𝑛−𝑣_𝑖,𝑗∇_𝑦^𝑛_𝑖,𝑗=0, (41) where 1≤𝑖≤𝑁_𝑥−1;𝑖+1≤𝑗≤𝑁_𝑦+𝑖−1;0≤𝑛≤𝑁_𝜏−1, (42)𝑤^𝑛_𝑖,𝑗 is the approximation of the exact solution 𝑊(𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) of (22) at the mesh point (𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) and operators Δ_𝑥,Δ_𝑦,∇_𝑥, and ∇_𝑦 are defined by the following: Δ_𝑥^𝑛_𝑖,𝑗=𝑤^𝑛_𝑖+1,𝑗−2𝑤^𝑛_𝑖,𝑗+𝑤^𝑛_{𝑖−1,𝑗}2ℎ²₁;Δ_𝑦^𝑛_𝑖,𝑗=𝑤^𝑛_𝑖,𝑗+1−2𝑤^𝑛_𝑖,𝑗+𝑤^𝑛_{𝑖,𝑗−1}2ℎ²₂;∇_𝑥^𝑛_𝑖,𝑗=𝑤^𝑛_𝑖+1,𝑗−𝑤^𝑛_{𝑖−1,𝑗}2ℎ₁;∇_𝑦^𝑛_𝑖,𝑗=𝑤^𝑛_𝑖,𝑗+1−𝑤^𝑛_{𝑖,𝑗−1}2ℎ₂, (43) where ℎ₁ and ℎ₂ are the stepsize spatial variables related by the following: ℎ₁=ℎ=Δ𝑥;ℎ₂=𝑚ℎ₁=𝑚Δ𝑥, (44) in our numerical domain Ω. In (41) time-dependent parameters are described by 𝜅^𝑛=𝜅(𝜏^𝑛) and 𝜃^𝑛=𝜃(𝜏^𝑛). Thus, for the internal mesh points, the scheme (41) can be written in an explicit five points stencil scheme, (see Figure 2): 𝑤^𝑛+1_𝑖,𝑗=𝑘1−ℎ²𝛼_𝑖,𝑗𝑤^𝑛_𝑖,𝑗+𝑘ℎ𝜌²𝛼2ℎ^𝑖,𝑗+𝛽_𝑖,𝑗𝑤^𝑛_𝑖+1,𝑗+𝑘ℎ𝜌²𝛼2ℎ^𝑖,𝑗−𝛽_𝑖,𝑗𝑤^𝑛_{𝑖−1,𝑗}+𝑘ℎ𝜌²𝛼2ℎ^𝑖,𝑗+𝛽_𝑖,𝑗−𝛾^𝑛_𝑖,𝑗𝑤^𝑛_𝑖,𝑗+1+𝑘ℎ𝜌²𝛼2ℎ^𝑖,𝑗−𝛽_𝑖,𝑗+𝛾^𝑛_𝑖,𝑗𝑤^𝑛_{𝑖,𝑗−1},1≤𝑖≤𝑁_𝑥−1;𝑖+1≤𝑗≤𝑁_𝑦+𝑖−1;0≤𝑛≤𝑁_𝜏−1, (45)

Figure 2 

Five-point versus nine-point stencils.

where 𝛼_𝑖,𝑗=𝜎²𝑚²𝑣_𝑖,𝑗;𝛽_𝑖,𝑗=𝜎𝜌𝑟−(1/2)𝑣_𝑖,𝑗2;𝛾^𝑛_𝑖,𝑗=𝜅^𝑛𝜃^𝑛−𝑣_𝑖,𝑗.2𝑚 (46)

Finally, it remains translating the initial condition (27) and the boundary conditions (28)–(32) of the transformed problem (22), to the numerical domain. The initial condition gives the following: 𝑤⁰_𝑖,𝑗=max𝑒^{𝑥𝑖/𝜎}𝜌−𝐸,0,0≤𝑖≤𝑁_𝑥;𝑖≤𝑗≤𝑁_𝑦+𝑖. (47) The Dirichlet conditions (28)-(29) are translated to the numerical boundary as follows 𝑤^𝑛_0,𝑗=0,0≤𝑗≤𝑁_𝑦−1,0≤𝑛≤𝑁_𝜏, (48) along the side 𝐴𝐷 of the rhomboid, and 𝑤^𝑛_𝑖,𝑖=𝑒^{𝑟𝜏𝑛+𝑥𝑖/𝜎}𝜌,1≤𝑖≤𝑁_𝑥, (49) along the side 𝐴𝐵 of the rhomboid. The Neumann boundary conditions are translated by replacing the partial derivatives values for the corresponding finite differences approximations values.

So, condition (32) takes the discretized backward form in order to involve mesh points of the numerical domain: 𝑤^𝑛_{𝑖,𝑁𝑦+𝑖}−𝑤^𝑛_{𝑖,𝑁𝑦+𝑖−1}𝑚ℎ=0, (50) or 𝑤^𝑛_{𝑖,𝑁𝑦+𝑖}=𝑤^𝑛_{𝑖,𝑁𝑦+𝑖−1},0≤𝑖≤𝑁_𝑥. (51)

The Neumann boundary condition (30) related to the numerical boundary 𝐵𝐶 of the rhomboid is regarded assuming that option price has a linear behaviour for large values of 𝑆. Hence, second partial derivatives terms vanish in (14) and taking into account the leading first partial derivative term of (14) for large values of 𝑆, one gets the following: 𝜕𝑉()𝜕𝜏∞,𝑣,𝜏=𝑟𝑆𝜕𝑉().𝜕𝑆∞,𝑣,𝜏 (52) From (52), the change of variables (21), and (30), it follows that the boundary numerical condition for the 𝐵𝐶 side of the rhomboid is as follows: 𝑤_𝑁^𝑛+1_𝑥,𝑗−𝑤^𝑛_{𝑁𝑥,𝑗}𝑘=𝑟𝑒^{𝑛𝑟𝑘+𝑏/𝜎}𝜌,𝑁_𝑥+1≤𝑗≤𝑁_𝑥+𝑁_𝑦. (53)

Summarizing, the numerical scheme is expressed by (45) for the internal mesh points, together with the initial condition (47) and the boundary conditions (48)–(51) and (53).

4. Consistency

Let 𝐹(𝑤^𝑛_𝑖,𝑗)=0 be the approximating difference equation (41). In accordance with [25, page 100], the scheme (41) is consistent with (22) if 𝑇^𝑛_𝑖,𝑗(𝑊)=𝐹𝑊^𝑛_𝑖,𝑗−𝐿𝑊^𝑛_𝑖,𝑗→0,asℎ₁=Δ𝑥→0,ℎ₂=Δ𝑦→0,𝑘=Δ𝜏→0, (54) where 𝑊^𝑛_𝑖,𝑗 denotes the theoretical solution of (22) evaluated at (𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) and 𝐿 is the operator: 𝐿(𝑊)=𝜕𝑊−1𝜕𝜏2𝜌²𝜎²𝑣𝜕²𝑊𝜕𝑥²+𝜕²𝑊𝜕𝑦²−1𝜌𝜎𝑟−2𝑣𝜕𝑊−1𝜕𝑥𝜌𝜎𝑟−2𝑣−𝜅(𝜃−𝑣)𝜕𝑊.𝜕𝑦 (55)

Assuming that 𝑊 is four times continuously differentiable with respect to 𝑥 and 𝑦 and twice with respect to 𝜏, and using Taylor expansion about (𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) it follows that 𝑊^𝑛+1_𝑖,𝑗−𝑊^𝑛_𝑖,𝑗𝑘=𝜕𝑊𝑥𝜕𝜏^𝑖,𝑦_𝑗,𝜏^𝑛+𝑘𝐸^𝑛_𝑖,𝑗(1), (56) where 𝐸^𝑛_𝑖,𝑗(1)=12𝜕²𝑊𝜕𝜏²𝑥_𝑖,𝑦_𝑗,𝜏^𝑛()𝐸,𝑛𝑘<𝛿<𝑛+1𝑘,_{𝑛𝑖,𝑗}(1≤1)2𝐷^𝑛_𝑖,𝑗(1)_max≡12𝜕max²𝑊𝜕𝜏²𝑥_𝑖,𝑦_𝑗;()𝑘.,𝜏𝑛𝑘≤𝜏≤𝑛+1 (57)(58) With respect to the approximations of the spatial partial derivatives, let us denote Δ_𝑥^𝑛_𝑖,𝑗=𝜕²𝑊𝜕𝑥²𝑥_𝑖,𝑦_𝑗,𝜏^𝑛+ℎ²₁𝐸^𝑛_𝑖,𝑗(2), (59) where 𝐸^𝑛_𝑖,𝑗(2)=1𝜕12⁴𝑊𝜕𝑥⁴𝜉₁,𝑦_𝑗,𝜏^𝑛;𝑥_𝑖−ℎ₁<𝜉₁<𝑥_𝑖+ℎ₁,𝐸^𝑛_𝑖,𝑗(2≤1)𝐷12^𝑛𝑗(2)_max≡1𝜕12max⁴𝑊𝜕𝑥⁴𝑥,𝑦_𝑗,𝜏^𝑛;𝑎≤𝑥≤𝑏. (60)(61) Analogously ∇_𝑥^𝑛_𝑖,𝑗=𝜕𝑊𝑥𝜕𝑥^𝑖,𝑦_𝑗,𝜏^𝑛+ℎ²₁𝐸^𝑛_𝑖,𝑗(3), (62) where 𝐸^𝑛_𝑖,𝑗(3)=16𝜕³𝑊𝜕𝑥³𝜉₂,𝑦_𝑗,𝜏^𝑛;𝑥_𝑖−ℎ₁<𝜉₂<𝑥_𝑖+ℎ₁,𝐸^𝑛_𝑖,𝑗(3≤1)6𝐷^𝑛_𝑗(3)_max≡16𝜕max³𝑊𝜕𝑥³𝑥,𝑦_𝑗,𝜏^𝑛;𝑎≤𝑥≤𝑏. (63)(64)

Otherwise, Δ_𝑦^𝑛_𝑖,𝑗=𝜕²𝑊𝜕𝑦²𝑥_𝑖,𝑦_𝑗,𝜏^𝑛+ℎ²₂𝐸^𝑛_𝑖,𝑗(4), (65)

where 𝐸^𝑛_𝑖,𝑗(4)=1𝜕12⁴𝑊𝜕𝑦⁴𝑥_𝑖,𝜉₃,𝜏^𝑛;𝑦_𝑗−ℎ₂<𝜉₃<𝑦_𝑗+ℎ₂,𝐸^𝑛_𝑖,𝑗(4≤1)𝐷12^𝑛𝑖(4)_max≡1𝜕12max⁴𝑊𝜕𝑦⁴𝑥_𝑖,𝑦,𝜏^𝑛;−𝑣₂+𝑚𝑥_𝑖<𝑦<−𝑣₁+𝑚𝑥_𝑖. (66)(67)

Let us write ∇_𝑦^𝑛_𝑖,𝑗=𝜕𝑊𝑥𝜕𝑦^𝑖,𝑦_𝑗,𝜏^𝑛+ℎ²₂𝐸^𝑛_𝑖,𝑗(5), (68)

where 𝐸^𝑛_𝑖,𝑗(5)=16𝜕³𝑊𝜕𝑦³𝑥_𝑖,𝜉₄,𝜏^𝑛,𝑦_𝑗−ℎ₂<𝜉₄<𝑦_𝑗+ℎ₂,𝐸^𝑛_𝑖,𝑗(5≤1)6𝐷^𝑛_𝑖(5)_max≡16𝜕max³𝑊𝜕𝑦³𝑥_𝑖,𝑦,𝜏^𝑛;−𝑣₂+𝑚𝑥_𝑖<𝑦<−𝑣₁+𝑚𝑥_𝑖. (69)(70)

From (56), (59), (62), (65), and (68), the local truncation error 𝑇^𝑛_𝑖,𝑗(𝑊) takes the following form: 𝑇^𝑛_𝑖,𝑗(𝑊)=𝐹𝑊^𝑛_𝑖,𝑗−𝐿𝑊^𝑛_𝑖,𝑗=𝑘𝐸^𝑛_𝑖,𝑗(1)−12𝑣_𝑖,𝑗𝜎²𝜌²ℎ²₁𝐸^𝑛_𝑖,𝑗(2)+ℎ²₂𝐸^𝑛_𝑖,𝑗(4)−1𝑟−2𝑣_𝑖,𝑗𝜎𝜌ℎ²₁𝐸^𝑛_𝑖,𝑗(3)−1𝑟−2𝑣_𝑖,𝑗𝜎𝜌−𝜅^𝑛𝜃^𝑛−𝑣_𝑖,𝑗ℎ²𝐸^𝑛_𝑖,𝑗(5). (71)

From (58), (61), (64), (67), (70), and (71), one gets the following: 𝑇^𝑛_𝑖,𝑗(𝑊≤𝑘)2𝐷^𝑛_𝑖,𝑗(1)_max+ℎ²₁614𝑣₂𝜎²𝜌²𝐷^𝑛_𝑗(2)_max+1𝑟−2𝑣_𝑖,𝑗𝜎𝜌𝐷^𝑛_𝑗(3)_max+ℎ²₂614𝑣₂𝜎²𝜌²𝐷^𝑛_𝑖(4)_max+1𝑟−2𝑣_𝑖,𝑗𝜎𝜌+𝜅^𝑛𝜃^𝑛−𝑣_𝑖,𝑗𝐷^𝑛_𝑖(5)_max=𝑂(𝑘)+𝑂ℎ²₁+𝑂ℎ²₂. (72)

Summarizing, the following result has been established.

Theorem 3. Assuming that the exact solution 𝑊 of (22) admits four times continuous partial derivatives with respect to the spatial variables 𝑥 and 𝑦 and twice continuous partial derivatives with respect to 𝜏, the scheme (41) is consistent with (22) and the local truncation error behaves 𝑇^𝑛_𝑖,𝑗(𝑊)=𝑂(𝑘)+𝑂ℎ²₁+𝑂ℎ²₂,𝑎𝑠Δ𝜏=𝑘→0,Δ𝑥=ℎ₁→0,Δ𝑦=ℎ₂→0. (73)

5. Positivity and Stability

We begin this section showing that under suitable relationship between the stepsize discretizations, the coefficients of scheme (45) are nonnegative. This fact together with the positivity of the boundary values will guarantee the positivity of the values of the solution in all time steps at every spatial mesh point.

Let us start considering the first coefficient of the right-hand side of (45). Note that 1−(𝑘/ℎ²)𝛼_𝑖,𝑗≥0 is equivalent to the following condition: 𝛼_𝑖,𝑗≤ℎ²𝑘. (74)

Now, from (45) and taking into account that 𝑣₁≤𝑣_𝑖,𝑗≤𝑣₂, it follows that 𝑘0≤ℎ²𝛼_𝑖,𝑗≤𝑘ℎ²𝑣₂𝜎²𝑚²;𝛼_𝑖,𝑗≤𝑣₂𝜎²𝑚². (75) Hence, (74) holds if 𝑘ℎ²≤𝑚²𝑣₂𝜎². (76) Note that the simultaneous nonnegativity of coefficients (𝜌²/2ℎ)𝛼_𝑖,𝑗+𝛽_𝑖,𝑗 and (𝜌²/2ℎ)𝛼_𝑖,𝑗−𝛽_𝑖,𝑗 is equivalent to the condition 𝛽_𝑖,𝑗≤𝜌²𝛼2ℎ^𝑖,𝑗. (77) From (46), condition (77) is equivalent to 1𝑟−2𝑣_𝑖,𝑗≤𝜎𝜌ℎ𝑣_𝑖,𝑗. (78) Note that if 𝑟=(1/2)𝑣_𝑖,𝑗, condition (78) holds true. If 𝑟≠(1/2)𝑣_𝑖,𝑗, then let us write (78) in the following form: 𝜎ℎ≤𝜌𝑣_𝑖,𝑗𝑟−(1/2)𝑣_𝑖,𝑗. (79) As the minimum of real function 𝑔(𝑧)=𝜎𝜌𝑧/|𝑟−(1/2)𝑧| in the interval 0<𝑧₁≤𝑧≤𝑧₂ is given by min𝑔(𝑧)∶𝑧₁≤𝑧≤𝑧₂𝜎=min𝜌𝑧_𝑖𝑟−(1/2)𝑧_𝑖;𝑖=1,2 (80) then from (80), condition (79) holds true if ℎ≤ℎ_𝑎=min_𝑠=1,22𝜎𝜌𝑣_𝑠2𝑟−𝑣_𝑠. (81) Note that the simultaneous nonnegativity of coefficients (𝜌²/2ℎ)𝛼_𝑖,𝑗+𝛽_𝑖,𝑗−𝛾^𝑛_𝑖,𝑗 and (𝜌²/2ℎ)𝛼_𝑖,𝑗−𝛽_𝑖,𝑗+𝛾^𝑛_𝑖,𝑗 is equivalent to the condition 𝛽_𝑖,𝑗−𝛾^𝑛_𝑖,𝑗≤𝜌²𝛼2ℎ^𝑖,𝑗. (82) For the sake of clarity, let us introduce the following: 𝐴𝜏^𝑛𝜅=𝜎𝑟−^𝑛𝜃^𝑛𝜌;𝐵𝜏^𝑛=𝜎2−𝜅^𝑛𝜌. (83) From (46) and (83), condition (82) can be written in the following form: 𝐴𝜏^𝑛−𝐵𝜏^𝑛𝑣_𝑖,𝑗≤𝜎²𝜌𝑣_𝑖,𝑗ℎ𝑚². (84) Note that if 𝐴(𝜏^𝑛)=𝐵(𝜏^𝑛)𝑣_𝑖,𝑗, then (84) holds true. If 𝐴(𝜏^𝑛)≠𝐵(𝜏^𝑛)𝑣_𝑖,𝑗, then (84) is equivalent to 𝜎ℎ≤²𝜌𝑣_𝑖,𝑗𝐴(𝜏^𝑛)−𝐵(𝜏^𝑛)𝑣_𝑖,𝑗𝑚². (85) As the continuous positive function 𝐺()=|𝐴(𝜏)(𝜏)𝑧|=𝜅(𝜏)𝜃(𝜏)𝜏,𝑧−𝐵𝜎𝑟−𝜌−𝜎2−𝜅(𝜏)𝜌𝑧 (86) is bounded, in the compact [0,𝑇]×[𝑣₁,𝑣₂], thus max𝐺(𝜏,𝑧);(𝜏,𝑧)∈[0,𝑇]×𝑣₁,𝑣₂=𝐶>0. (87) From (85) and (87), condition (85) holds true if ℎ≤ℎ_𝑏=𝜎²𝜌𝑣₁𝑚²𝐶. (88) Summarizing, the following result has been established.

Lemma 4. Let ℎ_𝑎, ℎ_𝑏 be defined by (81) and (88), respectively, and let ℎ₀=min{ℎ_𝑎,ℎ_𝑏}. Then, for ℎ<ℎ₀ and 𝑘 satisfying (76), all the coefficients of scheme (45) are nonnegative.

The nonnegativity of the coefficients of scheme (45) is not sufficient to guarantee the positivity of the constructed numerical solution due to the influence on the values at the boundaries of the numerical domain. From (47)–(51), (53), and Lemma 4 it is easy to show the positivity of 𝑤^𝑛_𝑖,𝑗 in all the boundaries of the rhomboid for all 𝑛. Hence, the following result has been established.

Theorem 5. Let (𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) be a mesh point in the domain described by Ω (39) with 𝑎≤𝑥_𝑖≤𝑏,𝑣₁≤𝑚𝑥_𝑖−𝑦_𝑗≤𝑣₂. Let 𝑤^𝑛_𝑖𝑗 be the numerical approximation of 𝑊(𝑥_𝑖,𝑦_𝑗,𝜏^𝑛) given by (45) together with the initial condition (47) and the boundary conditions (48)–(51) and (53). Assume that ℎ<ℎ₀, where ℎ₀ is given by Lemma 4, and 𝑘 satisfying (76), then 𝑤^𝑛_𝑖,𝑗≥0 for all 𝑛,𝑖,𝑗.

The next result is a local discrete maximum principle that guarantees the boundedness of the numerical solution at the internal mesh points of the domain and that will be used below to prove the stability of the scheme.

Lemma 6. Let (𝑖,𝑗) be spatial internal mesh point with 1≤𝑖≤𝑁_𝑥−1;𝑖+1≤𝑗≤𝑁_𝑦+𝑖−1 and let 𝑆^𝑛_𝑖,𝑗 be the pentaset 𝑆^𝑛_𝑖,𝑗=𝑤^𝑛_𝑖,𝑗;𝑤^𝑛_{𝑖−1,𝑗};𝑤^𝑛_𝑖+1,𝑗;𝑤^𝑛_{𝑖,𝑗−1};𝑤^𝑛_𝑖,𝑗+1. (89)
Then under conditions of Lemma 4, it follows that min𝑆^𝑛_𝑖,𝑗≤𝑤^𝑛+1_𝑖,𝑗≤max𝑆^𝑛_𝑖,𝑗, (90)
for 0≤𝑛≤𝑁_𝜏−1.

Proof. From (45), Lemma 4, and using that 𝜌²+𝜌²=1 and that the sum of all the five coefficients of the right-hand side of (45) is one, the result is established.

We conclude this section by introducing a concept of stability for the numerical solution of scheme (45), (47)–(51), (53). Let {𝑤^𝑛_𝑖,𝑗} be the numerical solution of scheme (45), (47)–(51), (53) and consider the matrix 𝑤^𝑛 in ℝ^{(𝑁𝑦+1)×(𝑁𝑥+1)} defined by 𝑤^𝑛=𝑤^𝑛_{𝑖,𝑗𝑖≤𝑗≤𝑁𝑦+𝑖,0≤𝑖≤𝑁𝑥},0≤𝑛≤𝑁_𝜏, (91) whose entries are the values of the numerical solution of the scheme for a fixed temporal step 𝑛.

Definition 7. With previous notation, we say that the solution {𝑤^𝑛_𝑖,𝑗} of the scheme (45), (47)–(51), (53) is ‖⋅‖₁-stable, if for small enough values of ℎ₁=ℎ=Δ𝑥, ℎ₂=𝑚Δ𝑥 and 𝑘=Δ𝜏, there is a constant 𝑀, independent of 𝑛, ℎ, and 𝑘, such that 𝑤^𝑛₁≤𝑀𝑤⁰₁, (92) for all 𝑛 such that 0≤𝑛≤𝑁_𝜏, with 𝑘𝑁_𝜏=𝑇.

Note that under conditions of Theorem 5, the entries of matrix 𝑤^𝑛 are nonnegative and let us denote 𝑤^𝑛=max𝑤^𝑛_𝑖,𝑗0≤𝑖≤𝑁_𝑥,𝑖≤𝑗≤𝑁_𝑦+𝑖. (93)

Now, we show using induction principle that there exists a number sequence {𝑄^𝑛}_{0≤𝑛≤𝑁𝜏} so that 𝑤^𝑛≤𝑄^𝑛. (94) Let us denote 𝛽=𝑒^𝑏/𝜎𝜌, (95) and for 𝑛=0, let us take 𝑄⁰=𝛽; then from the initial condition (47) one gets the following: 𝑤⁰≤𝑄⁰. (96) If 𝑛=1, for the internal mesh points, from Lemma 6, one gets the following: 0≤𝑤¹_𝑖,𝑗≤𝑄⁰;1≤𝑖≤𝑁_𝑥−1,𝑖+1≤𝑗≤𝑁_𝑦+𝑖−1. (97) From (48)–(51), and (53) and Theorem 5, since 𝑒^𝑟𝑘>1+𝑘𝑟 for 𝑘𝑟>0, it follows that 0≤𝑤¹_𝑖,𝑗≤max𝛽+𝑘𝑟𝛽,𝛽𝑒^𝑟𝑘=𝛽𝑒^𝑟𝑘,0≤𝑖≤𝑁_𝑥,𝑖≤𝑗≤𝑁_𝑦+𝑖. (98) Hence, taking 𝑄¹=𝛽𝑒^𝑟𝑘, it follows that 𝑤¹≤𝑄¹. (99) Let 𝑄^𝑛=𝛽𝑒^𝑛𝑟𝑘,𝑛≥1, (100) and assume the induction hypothesis 𝑤^𝑛≤𝑄^𝑛. (101) From Lemma 6, Theorem 5, and (48)–(51), and (53), one gets the following: 0≤𝑤^𝑛+1_𝑖,𝑗≤max𝛽𝑒^𝑛𝑟𝑘+𝑘𝑟𝛽𝑒^𝑛𝑟𝑘,𝛽𝑒^{(𝑛+1)𝑟𝑘}=𝛽𝑒^𝑛𝑟𝑘max1+𝑘𝑟,𝑒^𝑟𝑘=𝛽𝑒^{(𝑛+1)𝑟𝑘}=𝑄^𝑛+1. (102)

Thus, (100)-(101) holds true for 𝑛+1. Note that 0≤𝑛≤𝑁_𝜏,𝑘=𝑇/𝑁_𝜏, and using (100)-(101) one gets the following: 𝑤^𝑛≤𝛽𝑒^𝑛𝑟𝑘≤𝛽𝑒^{𝑁𝜏𝑟𝑘}=𝛽𝑒^𝑟𝑇. (103)

Let 𝑞 be the positive number, 0<𝑞<1, defined by 𝑞=1−𝐸𝑒^−𝑏/𝜎𝜌. Note that from (47) and definition of ‖⋅‖₁ one gets the following: 𝑤⁰₁=𝑁_𝑦+1𝑒^−𝑏/𝜎𝜌−𝐸=𝑞𝛽 (104) and from (103), 𝑤^𝑛₁≤𝑁_𝑦+1𝛽𝑒^𝑟𝑇. (105)

Finally, from (104)-(105) 𝑤^𝑛₁≤𝑒^𝑟𝑇𝑞𝑤⁰₁ (106)

Summarizing, the following result has been established.

Theorem 8. Consider the schemes (45), (47)–(51), (53) in the numerical domain Ω defined by (39) where 𝑥₀=𝑎,𝑥_𝑁𝑥=𝑏, and 𝑦_𝑗∈[𝑚𝑎−𝑣₂,𝑚𝑏−𝑣₁],𝜏^𝑛∈[0,𝑇]. Let ℎ₀ be defined by Lemma 4 and suppose that ℎ=Δ𝑥<ℎ₀ and 𝑘=Δ𝜏 is small enough so that condition (76) is satisfied. Then the scheme is ‖⋅‖₁-stable and 𝑤^𝑛₁≤𝑒^𝑟𝑇𝑞𝑤⁰₁,1≤𝑛≤𝑁_𝜏. (107)

6. Examples and Simulations

In this section, we illustrate the good properties of the proposed numerical scheme. Firstly, we consider an example where the closed form solution is available.

Example 1. Consider the European call option for the Heston model (3).
𝐸=100; 𝑇=1 year; 𝑟=0.01; 𝜅=2; 𝜃=0.01; 𝜎=0.1; 𝜌=0.5, 𝑣₁=0.4, and 𝑣₂=0.9.
The difference between the closed form solution 𝑈 and the numerical solution 𝑢 is shown in Figure 3, by respecting the stability condition (76) and Lemma 4.

Figure 3 

Error of numerical solution.

The following example shows the variation of the absolute and relative error of the numerical solution in light of the stability and positivity conditions hold at the strike for a given variance for different values of the stepsize discretizations ℎ and 𝑘, respectively.

Example 2. Consider the European call option for the Heston mode (3) with data of Example 1 at 𝑆=𝐸=100,𝑣=0.5. Table 1 shows the behaviour of the error with parameters ℎ and 𝑘. Notice that the numerical solution exhibits the expected second order convergence rate in space. Analogously, for fixed ℎ=0.1 one gets convergence time rate of 0.97. For ℎ=0.05 the rate is 1.

Table 1 

Errors and convergence rates.

In the next example, we illustrate the fact that if we do not respect the stability condition, the numerical results are bad and unreliable.

Example 3. Consider the problem of Example 1 at 𝑣=0.9 with ℎ=0.0082, 𝑘=0.0065, not verifying the stability condition (76). Figure 4 shows that spurious oscillations of the numerical solution appear when the stability condition is broken.

Figure 4 

Breaking positivity and stability conditions.

The next example shows the dependence of the option price on both the underlying asset and the variance.

Example 4. Consider the European call option for the Heston model with data.
𝐸=100; 𝑇=1 year;  𝑟=0.01;  𝜅=2;  𝜃=0.01;  𝜎=0.1;  𝜌=0.4 and 𝑣₁=0.2, 𝑣₂=0.9. Figure 5 shows the numerical solution in 3D.

Figure 5 

Numerical solution in 3𝐷.

Finally, in the next example, we show the numerical result of a situation with time dependent mean reversion rate 𝜅(𝜏) and long run variance 𝜃(𝜏).

Example 5. Consider the European call option for the Heston model with data: 𝐸=100; 𝑇=1 year; 𝑟=0.01; 𝜎=0.1; 𝜌=0.4 and 𝑣₁=0.4, 𝑣₂=0.9 with time-dependent parameters 𝜅(𝜏)=1.5+𝜏; 𝜃(𝜏)=0.005(1+2𝜏). Figure 6 shows the numerical results in this time-dependent case for 𝑣=0.9.

Figure 6 

A time-dependent case.

Acknowledgments

This work has been partially supported by the European Union in the FP7-PEOPLE-2012-ITN Program under Grant Agreement no. 304617 (FP7 Marie Curie Action, Project Multi-ITN STRIKE-Novel Methods in Computational Finance) and by the Spanish MEYC Grant DPI2010-20891-C02-01.