Abstract and Applied Analysis

Volume 2013 (2013), Article ID 246724, 11 pages

http://dx.doi.org/10.1155/2013/246724

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

Instituto Universitario de Matemática Multidisciplinar, Universitat Politècnica de València, Camino de Vera s/n, 46022 Valencia, Spain

Received 23 March 2013; Revised 18 May 2013; Accepted 20 May 2013

Academic Editor: T. Diagana

Copyright © 2013 R. Company et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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: 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: 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:

A European vanilla call option with strike price and maturity satisfies (3) together with the payoff final condition and the boundary conditions 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 and by the maximum absolute column sum norm 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: where , and .

By discretizing (3) one achieves the following scheme:

Note that if , 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 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 obtaining

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 where

Under the assumption of correlated variables with , (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: where

Solving (17) one gets the following: where the integration constant is related to the new variables by From (19) and (20), it follows the expression of the new variables

By denoting , (14) takes the following equivalent elliptic form: where from (21) takes the expression in terms of and , given by the following: 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: and means the variance of an underlying asset that must be positive, from (23) it follows that ; thus, Thus, the domain of the transformed problem is as follows: In the following we will assume that , because when , 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:

The boundary conditions (6), (9) of problem (3) are transformed into the following Dirichlet conditions: and the boundary condition (7) is transformed into the following Neumann condition:

Boundary condition (8) can be replaced with a more practicable homogeneous Neumann condition: see, for instance, [11, 13]. This boundary condition is transformed into the condition

*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:

*Remark 2. *For the case , fully correlated variables and , from (15), the discriminant , and (14) become a parabolic PDE. Following the techniques for reduction to canonical form, an appropriate substitution is ; and the transformed equation takes the following form:
where and .

#### 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 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 , . Let and of the form

Once is fixed, the admissible values of so that belongs to the rhomboid are given by the following: Note that with this discretization the boundary sides of the rhomboid are partitioned in the following way:

Thus, the rhomboid spatial domain together with its boundary contains mesh points. See Figure 1.

For the time variable , we consider the standard discretization

Hence, the numerical domain is defined by the following:

Note that because of (23) the value of at the spatial variables is as follows:

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: where is the approximation of the exact solution of (22) at the mesh point and operators , and are defined by the following: where and are the stepsize spatial variables related by the following: 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):

where

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: The Dirichlet conditions (28)-(29) are translated to the numerical boundary as follows along the side of the rhomboid, and 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: or

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: 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:

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 be the approximating difference equation (41). In accordance with [25, page 100], the scheme (41) is consistent with (22) if where denotes the theoretical solution of (22) evaluated at and is the operator:

Assuming that is four times continuously differentiable with respect to and and twice with respect to , and using Taylor expansion about it follows that where With respect to the approximations of the spatial partial derivatives, let us denote where Analogously where

Otherwise,

where

Let us write

where

From (56), (59), (62), (65), and (68), the local truncation error takes the following form:

From (58), (61), (64), (67), (70), and (71), one gets the following:

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
*

#### 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 is equivalent to the following condition:

Now, from (45) and taking into account that , it follows that Hence, (74) holds if Note that the simultaneous nonnegativity of coefficients and is equivalent to the condition From (46), condition (77) is equivalent to Note that if , condition (78) holds true. If , then let us write (78) in the following form: As the minimum of real function in the interval is given by then from (80), condition (79) holds true if Note that the simultaneous nonnegativity of coefficients and is equivalent to the condition For the sake of clarity, let us introduce the following: From (46) and (83), condition (82) can be written in the following form: Note that if , then (84) holds true. If , then (84) is equivalent to As the continuous positive function is bounded, in the compact , thus From (85) and (87), condition (85) holds true if Summarizing, the following result has been established.

Lemma 4. *Let , be defined by (81) and (88), respectively, and let . 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 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 and let be the pentaset
**Then under conditions of Lemma 4, it follows that
**for .*

*Proof. *From (45), Lemma 4, and using that 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 defined by 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
for all such that , with .

Note that under conditions of Theorem 5, the entries of matrix are nonnegative and let us denote

Now, we show using induction principle that there exists a number sequence so that Let us denote and for , let us take ; then from the initial condition (47) one gets the following: If , for the internal mesh points, from Lemma 6, one gets the following: From (48)–(51), and (53) and Theorem 5, since for , it follows that Hence, taking , it follows that Let and assume the induction hypothesis From Lemma 6, Theorem 5, and (48)–(51), and (53), one gets the following:

Thus, (100)-(101) holds true for . Note that , and using (100)-(101) one gets the following:

Let be the positive number, , defined by . Note that from (47) and definition of one gets the following: and from (103),

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 . 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
*

#### 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).

; year; ; ; ; ; , , and .

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.

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 . 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 one gets convergence time rate of 0.97. For the rate is 1.

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 with , , not verifying the stability condition (76). Figure 4 shows that spurious oscillations of the numerical solution appear when the stability condition is broken.

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.

; year; ; ; ; ; and , . Figure 5 shows the numerical solution in 3D.

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: ; year; ; ; and , with time-dependent parameters ; . Figure 6 shows the numerical results in this time-dependent case for .

#### 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.

#### References

- B. Dupire, “Arbitrage pricing with stochastic volatility,” Tech. Rep., Banque Paribas Swaps and Options Research Team Monograph, 1993. View at Google Scholar
- J. Hull and A. White, “The pricing of options with stochastic volatilities,”
*Journal of Finance*, vol. 42, pp. 281–300, 1987. View at Google Scholar - S. L. Heston, “A closed-form solution for options with stochastic volatility with applications to bond and currency options,”
*The Review of Financial Studies*, vol. 6, pp. 327–343, 1993. View at Google Scholar - P. Hagan, D. Kumar, and A. S. Lesniewski, “Managing smile risk,”
*Wilmott Magazine*, vol. 15, pp. 84–108, 2002. View at Google Scholar - A. Pascucci,
*PDE and Martingale Methods in Option Pricing*, vol. 2 of*Bocconi & Springer Series*, Springer, Milan, Italy, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. Mikhailov and U. Nögel, “Heston's stochastic volatility model implementation, calibration and some extensions,”
*Wilmott Magazine*, vol. 4, pp. 74–79, 2003. View at Google Scholar - A. Elices, “Models with time-dependent parameters using transform methods: application to Heston's model,” http://arxiv.org/abs/0708.2020.
- E. Benhamou, E. Gobet, and M. Miri, “Time dependent Heston model,”
*SIAM Journal on Financial Mathematics*, vol. 1, pp. 289–325, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - N. Hilber, A. M. Matache, and C. Schwab, “Sparse wavelet methods for option pricing under stochastic volatility,”
*Journal of Computational Finance*, vol. 8, no. 4, pp. 1–42, 2005. View at Google Scholar - W. Zhu and D. A. Kopriva, “A spectral element approximation to price European options with one asset and stochastic volatility,”
*Journal of Scientific Computing*, vol. 42, no. 3, pp. 426–446, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - N. Clarke and K. Parrott, “Multigrid for american option pricing with stochastic volatility,”
*Applied Mathematical Finance*, vol. 6, no. 3, pp. 177–195, 1999. View at Google Scholar - R. Sheppard,
*Pricing equity derivatives under stochastic volatility: a partial differential equation approach [Ph.D. thesis]*, Faculty of Science, University of the Witwatersrand, 2007. - B. Düring and M. Fournié, “High-order compact finite difference scheme for option pricing in stochastic volatility models,”
*Journal of Computational and Applied Mathematics*, vol. 236, no. 17, pp. 4462–4473, 2012. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - D. J. Duffy,
*Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach*, vol. 312, John Wiley & Sons, Chichester, UK, 2006. View at MathSciNet - K. J. Hout and S. Foulon, “ADI finite difference schemes for option pricing in the Heston model with correlation,”
*International Journal of Numerical Analysis and Modeling*, vol. 7, no. 2, pp. 303–320, 2010. View at Google Scholar · View at MathSciNet - R. Zvan, P. A. Forsyth, and K. R. Vetzal, “Negative coefficients in two-factor 25 option pricing models,”
*Journal of Computational Finance*, vol. 7, no. 1, pp. 37–73, 2003. View at Google Scholar - G. H. Golub and C. F. Van Loan,
*Matrix Computations*, Johns Hopkins University Press, London, UK, 3rd edition, 1996. View at MathSciNet - R. Company, L. Jódar, and J.-R. Pintos, “Consistent stable difference schemes for nonlinear Black-Scholes equations modelling option pricing with transaction costs,”
*Mathematical Modelling and Numerical Analysis*, vol. 43, no. 6, pp. 1045–1061, 2009. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - R. Company, L. Jódar, and J.-R. Pintos, “Numerical analysis and computing for option pricing models in illiquid markets,”
*Mathematical and Computer Modelling*, vol. 52, no. 7-8, pp. 1066–1073, 2010. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M.-C. Casabán, R. Company, L. Jódar, and J.-R. Pintos, “Numerical analysis and computing of a non-arbitrage liquidity model with observable parameters for derivatives,”
*Computers & Mathematics with Applications*, vol. 61, no. 8, pp. 1951–1956, 2011. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - P. R. Garabedian,
*Partial Differential Equations*, AMS Chelsea, Providence, RI, USA, 1998. View at MathSciNet - R. Ehrhardt,
*Discrete artificial boundary conditions [Ph.D. thesis]*, Technische Universitat Berlin, 2001. - R. Kangro and R. Nicolaides, “Far field boundary conditions for Black-Scholes equations,”
*SIAM Journal on Numerical Analysis*, vol. 38, no. 4, pp. 1357–1368, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M. Ehrhardt and R. E. Mickens, “A fast, stable and accurate numerical method for the Black-Scholes equation of American options,”
*International Journal of Theoretical and Applied Finance*, vol. 11, no. 5, pp. 471–501, 2008. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - G. D. Smith,
*Numerical Solution of Partial Differential Equations: Finite Difference Methods*, Clarendon Press, Oxford, UK, 3rd edition, 1985. View at MathSciNet