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

Authors: R. Company; L. Jódar; M. Fakharany; M.-C. Casabán
Journal: Abstract and Applied Analysis
 (2013)

---

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


---

## Body

## 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)dS(t)=pS(S(t),v(t),t)dt+qS(S(t),v(t),t)dW1,dv(t)=pv(S(t),v(t),t)dt+qv(S(t),v(t),t)dW2,
where the tradable security S(t) and its variance v(t) are correlated; that is, the Wiener stochastic processes satisfy dW1dW2=ρdt, and the functional form of pS,qS,pv, and qv 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)dS(t)=μS(t)dt+v(t)s(t)dW1,dv(t)=κ(θ-v(t))dt+σv(t)dW2,dW1dW2=ρdt,
where W1, W2 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 U=U(S,v,t) of a contingent claim:
(3)∂U∂t+12vS2∂2U∂S2+ρσvS∂2U∂S∂v+12σ2v∂2U∂v2+rS∂U∂S+κ(θ-v)∂U∂v-rU=0,(4)0<S<∞,0<v<∞,0≤t<T.A European vanilla call option with strike priceE and maturity T satisfies (3) together with the payoff final condition
(5)U(S,v,T)=max(S-E,0),
and the boundary conditions
(6)U(0,v,t)=0,(7)∂U∂S(∞,v,t)=1,(8)rS∂U∂S(S,0,t)+κ¯θ¯∂U∂v(S,0,t)-rU(S,0,t)+∂U∂t(S,0,t)=0,(9)U(S,∞,t)=S,
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 matrixA=(aij)m×n in ℝm×n, we denote by ∥A∥=max{|aij|;1≤i≤m,1≤j≤n} and by ∥A∥1 the maximum absolute column sum norm
(10)∥A∥1=max1≤j≤n{∑i=1m‍|aij|};
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:
(11)U(Si,vj,τn)≈ui,jn;∂U∂τ(Si,vj,τn)≈ui,jn+1-ui,jnk;∂U∂S(Si,vj,τn)≈ui+1,jn-ui-1,jn2h1;∂2U∂S2(Si,vj,τn)≈ui-1,jn-2ui,jn+ui+1,jnh12;∂U∂v(Si,vj,τn)≈ui,j+1n-ui,j-1n2h2;∂2U∂v2(Si,vj,τn)≈ui,j-1n-2ui,jn+ui,j+1nh22;∂2U∂S∂v(Si,vj,τn)≈ui+1,j+1n+ui-1,j-1n-ui-1,j+1n-ui+1,j-1n4h1h2,
where τ=T-t,k=Δτ,h1=ΔS,h2=Δv,Si=ih1,vj=jh2, and τn=nk.By discretizing (3) one achieves the following scheme:
(12)ui,jn+1=(1-rk-kh12vjSi2-kh22σ2vj)ui,jn+k2Si(vjSih12+rh1)ui+1,jn+k2Si(vjSih12-rh1)ui-1,jn+k2(σ2vjh22+κ(θ-vj)h2)ui,j+1n+k2(σ2vjh22-κ(θ-vj)h2)ui,j-1n+k4h1h2ρσvjSi(ui+1,j+1n+ui-1,j-1n-ui-1,j+1n-ui+1,j-1n).Note that ifρ≠0, that is, when variables S and v 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 σ, vj, and Si 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 termrU by means of the substitution
(13)V=erτU,τ=T-t,
obtaining
(14)∂V∂τ=12vS2∂2V∂S2+ρσvS∂2V∂S∂v+12σ2v∂2V∂v2+rS∂V∂S+κ(θ-v)∂V∂v.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
(15)Δ(S,v)=B2-4AC=σ2v2S2(ρ2-1),
where
(16)B=ρσvS;A=12vS2;C=12σ2v.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:
(17)dvdS=σ(ρ+i1-ρ2)S=σ(ρ+iρ~)S,
where
(18)ρ~=1-ρ2.Solving (17) one gets the following:
(19)v+Z0=σ(ρ+iρ~)lnS,
where the integration constant Z0 is related to the new variables by
(20)Z0=y+ix.
From (19) and (20), it follows the expression of the new variables
(21)x=ρ~σlnS;y=ρσlnS-v.By denotingW(x,y,τ)=V(S,v,τ), (14) takes the following equivalent elliptic form:
(22)∂W∂τ=12ρ~2σ2v(∂2W∂x2+∂2W∂y2)+ρ~σ(r-12v)∂W∂x+(ρσ(r-12v)-κ(θ-v))∂W∂y,
where v from (21) takes the expression in terms of x and y, given by the following:
(23)v=ρρ~x-y.
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:
(24)-∞<x<+∞,0<τ≤T,
and v means the variance of an underlying asset that must be positive, from (23) it follows that (ρ/ρ~)x-y>0; thus,
(25)-∞<y<ρρ~x=mx;m=ρρ~.
Thus, the domain of the transformed problem is as follows:
(26)D={(x,y,τ);x∈ℝ;y<mx;0<τ≤T}.
In the following we will assume that m=(ρ/ρ~)≠0, because when m=0, ρ=0 and variables S and v 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:
(27)W(x,y,0)=max(ex/σρ~-E,0);x∈ℝ,y<mx.The boundary conditions (6), (9) of problem (3) are transformed into the following Dirichlet conditions:
(28)W(x,y,τ)→0,asx→-∞,(29)W(x,y,τ)→erτ+(x/σρ~),asmx-y→+∞,
and the boundary condition (7) is transformed into the following Neumann condition:
(30)σρ~∂W∂x(x,y,τ)+σρ∂W∂y(x,y,τ)→erτ+x/σρ~,asx→+∞.Boundary condition (8) can be replaced with a more practicable homogeneous Neumann condition:
(31)∂U∂v(S,0,t)=0;
see, for instance, [11, 13]. This boundary condition is transformed into the condition
(32)∂W∂y(x,y,τ)→0,asmx-y→0+.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 [S1,S2]×[v1,v2] is transformed into the rhomboid ABCD; see Figure 1, where the sides are described by the following:
(33)AD¯={(x,y)∈ℝ2;x=a=ρ~σlnS1;y=ma-v;v1≤v≤v2},AB¯={(x,y)∈ℝ2;a≤x≤b=ρ~σlnS2;y=mx-v2},BC¯={(x,y)∈ℝ2;x=b;y=mb-v;v1≤v≤v2},CD¯={(x,y)∈ℝ2;a≤x≤b;y=mx-v1}.Figure 1
RhomboidABCD.Remark 2.
For the case|ρ|=1, fully correlated variables S and v, from (15), the discriminant Δ=0, and (14) become a parabolic PDE. Following the techniques for reduction to canonical form, an appropriate substitution is x=S; y=v-ρσlnS and the transformed equation takes the following form:
(34)∂W∂τ=12vx2∂2W∂x2+rx∂W∂x+(ρσ(12v-r)+κ(θ-v))∂W∂y,
where |ρ|=1 and v=y+ρσlnx.


## 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 m=(ρ/ρ~)≠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 S is available and generally accepted. In an analogous way, considering an admissible range of the variance v, we can identify a convenient-bounded numerical domain [S1,S2]×[v1,v2] in the Sv plane.In accordance with (33), let xi∈[a,b], with xi=a+ih, 0≤i≤Nx,Δx=h=(b-a)/Nx. Let Ny=(v2-v1)/mh and yj∈[ma-v2,mb-v1] of the form
(35)y0=ma-v2;yj=y0+jmh;Δy=mh=mΔx.Oncei is fixed, the admissible values of j so that (xi,yi) belongs to the rhomboid ABCD are given by the following:
(36)i≤j≤Ny+i;0≤i≤Nx.
Note that with this discretization the boundary sides of the rhomboid are partitioned in the following way:
(37)P(AD¯)={(x0,yj);0≤j≤Ny},P(AB¯)={(xi,yi);0≤i≤Nx},P(BC¯)={(xNx,yj);Nx≤j≤Nx+Ny},P(CD¯)={(xi,yi+Ny);0≤i≤Nx}.Thus, the rhomboid spatial domain together with its boundary contains(Nx+1)(Ny+1) mesh points. See Figure 1.For the time variableτ, we consider the standard discretization
(38)τn=kn;0≤n≤Nτ;Δτ=k=TNτ.Hence, the numerical domainΩ is defined by the following:
(39)Ω={(xi,yj,τn)∣0≤i≤Nx;i≤j≤Ny+i;0≤n≤Nτ}.Note that because of (23) the value of v at the spatial variables (xi,yj) is as follows:
(40)vi,j=mxi-yj;0≤i≤Nx;i≤j≤Ny+i.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:
(41)F(wi,jn)=wi,jn+1-wi,jnk-12vi,jσ2ρ~2((Δx)i,jn+(Δy)i,jn)-σρ~(r-12vi,j)(∇x)i,jn-[σρ(r-12vi,j)-κn(θn-vi,j)](∇y)i,jn=0,
where
(42)1≤i≤Nx-1;i+1≤j≤Ny+i-1;0≤n≤Nτ-1,wi,jn is the approximation of the exact solution W(xi,yj,τn) of (22) at the mesh point (xi,yj,τn) and operators Δx,Δy,∇x, and ∇y are defined by the following:
(43)(Δx)i,jn=wi+1,jn-2wi,jn+wi-1,jn2h12;(Δy)i,jn=wi,j+1n-2wi,jn+wi,j-1n2h22;(∇x)i,jn=wi+1,jn-wi-1,jn2h1;(∇y)i,jn=wi,j+1n-wi,j-1n2h2,
where h1 and h2 are the stepsize spatial variables related by the following:
(44)h1=h=Δx;h2=mh1=mΔx,
in our numerical domain Ω. In (41) time-dependent parameters are described by κn=κ(τn) and θn=θ(τn). Thus, for the internal mesh points, the scheme (41) can be written in an explicit five points stencil scheme, (see Figure 2):
(45)wi,jn+1=(1-kh2αi,j)wi,jn+kh(ρ22hαi,j+βi,j)wi+1,jn+kh(ρ22hαi,j-βi,j)wi-1,jn+kh(ρ~22hαi,j+βi,j-γi,jn)wi,j+1n+kh(ρ~22hαi,j-βi,j+γi,jn)wi,j-1n,1≤i≤Nx-1;i+1≤j≤Ny+i-1;0≤n≤Nτ-1,Figure 2
Five-point versus nine-point stencils.where(46)αi,j=σ2m2vi,j;βi,j=σρ~(r-(1/2)vi,j)2;γi,jn=κn(θn-vi,j)2m.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:
(47)wi,j0=max(exi/σρ~-E,0),0≤i≤Nx;i≤j≤Ny+i.
The Dirichlet conditions (28)-(29) are translated to the numerical boundary as follows
(48)w0,jn=0,0≤j≤Ny-1,0≤n≤Nτ,
along the side AD of the rhomboid, and
(49)wi,in=erτn+xi/σρ~,1≤i≤Nx,
along the side AB 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:
(50)wi,Ny+in-wi,Ny+i-1nmh=0,
or
(51)wi,Ny+in=wi,Ny+i-1n,0≤i≤Nx.The Neumann boundary condition (30) related to the numerical boundary BC of the rhomboid is regarded assuming that option price has a linear behaviour for large values of S. Hence, second partial derivatives terms vanish in (14) and taking into account the leading first partial derivative term of (14) for large values of S, one gets the following:
(52)∂V∂τ(∞,v,τ)=rS∂V∂S(∞,v,τ).
From (52), the change of variables (21), and (30), it follows that the boundary numerical condition for the BC side of the rhomboid is as follows:
(53)wNx,jn+1-wNx,jnk=renrk+b/σρ~,Nx+1≤j≤Nx+Ny.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

LetF(wi,jn)=0 be the approximating difference equation (41). In accordance with [25, page 100], the scheme (41) is consistent with (22) if
(54)Ti,jn(W)=F(Wi,jn)-L(Wi,jn)→0,ash1=Δx→0,h2=Δy→0,k=Δτ→0,
where Wi,jn denotes the theoretical solution of (22) evaluated at (xi,yj,τn) and L is the operator:
(55)L(W)=∂W∂τ-12ρ~2σ2v(∂2W∂x2+∂2W∂y2)-ρ~σ(r-12v)∂W∂x-(ρσ(r-12v)-κ(θ-v))∂W∂y.Assuming thatW is four times continuously differentiable with respect to x and y and twice with respect to τ, and using Taylor expansion about (xi,yj,τn) it follows that
(56)Wi,jn+1-Wi,jnk=∂W∂τ(xi,yj,τn)+kEi,jn(1),
where
(57)Ei,jn(1)=12∂2W∂τ2(xi,yj,τn),nk<δ<(n+1)k,(58)|Ei,jn(1)|≤12|Di,jn(1)|max≡12max{|∂2W∂τ2(xi,yj,τ)|;ink≤τ≤(n+1)k|∂2W∂τ2(xi,yj,τ)|}.
With respect to the approximations of the spatial partial derivatives, let us denote
(59)(Δx)i,jn=∂2W∂x2(xi,yj,τn)+h12Ei,jn(2),
where
(60)Ei,jn(2)=112∂4W∂x4(ξ1,yj,τn);xi-h1<ξ1<xi+h1,(61)|Ei,jn(2)|≤112|Djn(2)|max≡112max{|∂4W∂x4(x,yj,τn)|;a≤x≤b}.
Analogously
(62)(∇x)i,jn=∂W∂x(xi,yj,τn)+h12Ei,jn(3),
where
(63)Ei,jn(3)=16∂3W∂x3(ξ2,yj,τn);xi-h1<ξ2<xi+h1,(64)|Ei,jn(3)|≤16|Djn(3)|max≡16max{|∂3W∂x3(x,yj,τn)|;a≤x≤b}.Otherwise,(65)(Δy)i,jn=∂2W∂y2(xi,yj,τn)+h22Ei,jn(4),where(66)Ei,jn(4)=112∂4W∂y4(xi,ξ3,τn);yj-h2<ξ3<yj+h2,(67)|Ei,jn(4)|≤112|Din(4)|max≡112max{|∂4W∂y4(xi,y,τn)|;-v2+mxi<y<-v1+mxi|∂4W∂y4(xi,y,τn)|}.Let us write(68)(∇y)i,jn=∂W∂y(xi,yj,τn)+h22Ei,jn(5),where(69)Ei,jn(5)=16∂3W∂y3(xi,ξ4,τn),yj-h2<ξ4<yj+h2,(70)|Ei,jn(5)|≤16|Din(5)|max≡16max{|∂3W∂y3(xi,y,τn)|;-v2+mxi<y<-v1+mxi|∂3W∂y3(xi,y,τn)|}.From (56), (59), (62), (65), and (68), the local truncation error Ti,jn(W) takes the following form:
(71)Ti,jn(W)=F(Wi,jn)-L(Wi,jn)=kEi,jn(1)-12vi,jσ2ρ~2(h12Ei,jn(2)+h22Ei,jn(4))-(r-12vi,j)σρ~(h12Ei,jn(3))-{(r-12vi,j)σρ-κn(θn-vi,j)}h2Ei,jn(5).From (58), (61), (64), (67), (70), and (71), one gets the following:
(72)|Ti,jn(W)|≤k2|Di,jn(1)|max+h126{14v2σ2ρ~2|Djn(2)|max+|r-12vi,j|σρ~|Djn(3)|max}+h226{14v2σ2ρ~2|Din(4)|max+(|r-12vi,j|σρ~+κn|θn-vi,j|)|Din(5)|max}=O(k)+O(h12)+O(h22).Summarizing, the following result has been established.Theorem 3.
Assuming that the exact solutionW of (22) admits four times continuous partial derivatives with respect to the spatial variables x and y and twice continuous partial derivatives with respect to τ, the scheme (41) is consistent with (22) and the local truncation error behaves
(73)Ti,jn(W)=O(k)+O(h12)+O(h22),asΔτ=k→0,Δx=h1→0,Δy=h2→0.


## 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-(k/h2)αi,j≥0 is equivalent to the following condition:
(74)αi,j≤h2k.Now, from (45) and taking into account that v1≤vi,j≤v2, it follows that
(75)0≤kh2αi,j≤kh2v2σ2m2;αi,j≤v2σ2m2.
Hence, (74) holds if
(76)kh2≤m2v2σ2.
Note that the simultaneous nonnegativity of coefficients (ρ2/2h)αi,j+βi,j and (ρ2/2h)αi,j-βi,j is equivalent to the condition
(77)|βi,j|≤ρ22hαi,j.
From (46), condition (77) is equivalent to
(78)|r-12vi,j|≤σρ~hvi,j.
Note that if r=(1/2)vi,j, condition (78) holds true. If r≠(1/2)vi,j, then let us write (78) in the following form:
(79)h≤σρ~vi,j|r-(1/2)vi,j|.
As the minimum of real function g(z)=σρ~z/|r-(1/2)z| in the interval 0<z1≤z≤z2 is given by
(80)min{g(z):z1≤z≤z2}=min{σρ~zi|r-(1/2)zi|;i=1,2}
then from (80), condition (79) holds true if
(81)h≤ha=mins=1,2{2σρ~vs|2r-vs|}.
Note that the simultaneous nonnegativity of coefficients (ρ~2/2h)αi,j+βi,j-γi,jn and (ρ~2/2h)αi,j-βi,j+γi,jn is equivalent to the condition
(82)|βi,j-γi,jn|≤ρ~22hαi,j.
For the sake of clarity, let us introduce the following:
(83)A(τn)=σr-κnθnρ;B(τn)=σ2-κnρ.
From (46) and (83), condition (82) can be written in the following form:
(84)|A(τn)-B(τn)vi,j|≤σ2ρ~vi,jhm2.
Note that if A(τn)=B(τn)vi,j, then (84) holds true. If A(τn)≠B(τn)vi,j, then (84) is equivalent to
(85)h≤σ2ρ~vi,j|A(τn)-B(τn)vi,j|m2.
As the continuous positive function
(86)G(τ,z)=|A(τ)-B(τ)z|=|σr-κ(τ)θ(τ)ρ-(σ2-κ(τ)ρ)z|
is bounded, in the compact [0,T]×[v1,v2], thus
(87)max{G(τ,z);(τ,z)∈[0,T]×[v1,v2]}=C>0.
From (85) and (87), condition (85) holds true if
(88)h≤hb=σ2ρ~v1m2C.
Summarizing, the following result has been established.Lemma 4.
Letha, hb be defined by (81) and (88), respectively, and let h0=min{ha,hb}. Then, for h<h0 and k 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 wi,jn in all the boundaries of the rhomboid for all n. Hence, the following result has been established.Theorem 5.
Let(xi,yj,τn) be a mesh point in the domain described by Ω (39) with a≤xi≤b,v1≤mxi-yj≤v2. Let wijn be the numerical approximation of W(xi,yj,τn) given by (45) together with the initial condition (47) and the boundary conditions (48)–(51) and (53). Assume that h<h0, where h0 is given by Lemma 4, and k satisfying (76), then wi,jn≥0 for all n,i,j.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(i,j) be spatial internal mesh point with 1≤i≤Nx-1;i+1≤j≤Ny+i-1 and let Si,jn be the pentaset
(89)Si,jn={wi,jn;wi-1,jn;wi+1,jn;wi,j-1n;wi,j+1n}.
Then under conditions of Lemma4, it follows that
(90)minSi,jn≤wi,jn+1≤maxSi,jn,
for0≤n≤Nτ-1.Proof.
From (45), Lemma 4, and using that ρ2+ρ~2=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 {wi,jn} be the numerical solution of scheme (45), (47)–(51), (53) and consider the matrix wn in ℝ(Ny+1)×(Nx+1) defined by
(91)wn=[wi,jn]i≤j≤Ny+i,0≤i≤Nx,0≤n≤Nτ,
whose entries are the values of the numerical solution of the scheme for a fixed temporal step n.Definition 7.
With previous notation, we say that the solution{wi,jn} of the scheme (45), (47)–(51), (53) is ∥·∥1-stable, if for small enough values of h1=h=Δx, h2=mΔx and k=Δτ, there is a constant M, independent of n, h, and k, such that
(92)∥wn∥1≤M∥w0∥1,
for all n such that 0≤n≤Nτ, with kNτ=T.Note that under conditions of Theorem5, the entries of matrix wn are nonnegative and let us denote
(93)∥wn∥=max{wi,jn0≤i≤Nx,i≤j≤Ny+i}.Now, we show using induction principle that there exists a number sequence{Qn}0≤n≤Nτ so that
(94)∥wn∥≤Qn.
Let us denote
(95)β=eb/σρ~,
and for n=0, let us take Q0=β; then from the initial condition (47) one gets the following:
(96)∥w0∥≤Q0.
If n=1, for the internal mesh points, from Lemma 6, one gets the following:
(97)0≤wi,j1≤Q0;1≤i≤Nx-1,i+1≤j≤Ny+i-1.
From (48)–(51), and (53) and Theorem 5, since erk>1+kr for kr>0, it follows that
(98)0≤wi,j1≤max{β+krβ,βerk}=βerk,0≤i≤Nx,i≤j≤Ny+i.
Hence, taking Q1=βerk, it follows that
(99)∥w1∥≤Q1.
Let
(100)Qn=βenrk,n≥1,
and assume the induction hypothesis
(101)∥wn∥≤Qn.
From Lemma 6, Theorem 5, and (48)–(51), and (53), one gets the following:
(102)0≤wi,jn+1≤max{βenrk+krβenrk,βe(n+1)rk}=βenrkmax{1+kr,erk}=βe(n+1)rk=Qn+1.Thus, (100)-(101) holds true for n+1. Note that 0≤n≤Nτ,k=T/Nτ, and using (100)-(101) one gets the following:
(103)∥wn∥≤βenrk≤βeNτrk=βerT.Letq be the positive number, 0<q<1, defined by q=1-Ee-b/σρ~. Note that from (47) and definition of ∥·∥1 one gets the following:
(104)∥w0∥1=(Ny+1)(e-b/σρ~-E)=qβ
and from (103),
(105)∥wn∥1≤(Ny+1)βerT.Finally, from (104)-(105)
(106)∥wn∥1≤erTq∥w0∥1Summarizing, the following result has been established.Theorem 8.
Consider the schemes (45), (47)–(51), (53) in the numerical domain Ω defined by (39) where x0=a,xNx=b, and yj∈[ma-v2,mb-v1],τn∈[0,T]. Let h0 be defined by Lemma 4 and suppose that h=Δx<h0 and k=Δτ is small enough so that condition (76) is satisfied. Then the scheme is ∥·∥1-stable and
(107)∥wn∥1≤erTq∥w0∥1,1≤n≤Nτ.


## 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).
E
=
100; T=1 year; r=0.01; κ=2; θ=0.01; σ=0.1; ρ=0.5, v1=0.4, and v2=0.9.
The difference between the closed form solutionU and the numerical solution u 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 discretizationsh and k, respectively.Example 2.
Consider the European call option for the Heston mode (3) with data of Example 1 at S=E=100,v=0.5. Table 1 shows the behaviour of the error with parameters h and k. Notice that the numerical solution exhibits the expected second order convergence rate in space. Analogously, for fixed h=0.1 one gets convergence time rate of 0.97. For h=0.05 the rate is 1.Table 1
Errors and convergence rates.
k
=
0.015
k
=
0.0075
h
0.1
0.05
0.1
0.05
|
U
-
u
|
0.06056
0.01594
0.03079
0.00797
Convergence rate
1.93
1.95
Relative error
2.3
×
10
-
3
1.7907
×
10
-
4
1.2
×
10
-
3
8.9537
×
10
-
5In 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 Example1 at v=0.9 with h=0.0082, k=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.
E
=
100; T=1 year;  r=0.01;  κ=2;  θ=0.01;  σ=0.1;  ρ=0.4 and v1=0.2, v2=0.9. Figure 5 shows the numerical solution in 3D.Figure 5
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:E=100; T=1 year; r=0.01; σ=0.1; ρ=0.4 and v1=0.4, v2=0.9 with time-dependent parameters κ(τ)=1.5+τ; θ(τ)=0.005(1+2τ). Figure 6 shows the numerical results in this time-dependent case for v=0.9.Figure 6
A time-dependent case.



---

*Source: 246724-2013-06-12.xml*