{'fraction_non_alphanumeric': 0.07652917308568577, 'fraction_numerical': 0.023759411547724953, 'mean_word_length': 3.347104319939634, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 42, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 22, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A hybrid approach for the implementation of the Heston model

We propose a hybrid tree-finite difference method in order to approximate the Heston model. We prove the convergence by embedding the procedure in a bivariate Markov chain and we study the convergence of European and American option prices. We finally provide numerical experiments that give accurate option prices in the Heston model, showing the reliability and the efficiency of the algorithm.


Introduction
Black-Scholes model was the most popular model for derivative pricing and hedging, although it has shown several problems with capturing dramatic moves in financial markets. In fact, the assumption of a constant volatility in the Black-Scholes model over the lifetime of the derivative is not realistic. As an alternative to the Black-Scholes model, stochastic volatility models emerged. The Heston model [17] is perhaps the most popular stochastic volatility model, allowing one to obtain closed-formulae in the European case using Fourier transform. In the American option pricing case, the main algorithms turn out to be tree methods, Fouriercosine methods and finite difference methods. Approximating trees for the Heston model have been considered in different papers, see e.g. [22], [13] , [14], [18], [15]. The tree approach of Vellekoop and Nieuwenhuis [25] actually provides at our knowledge the best tree procedure in the literature. They use an approach which is based on a modification of an explicitly defined stock price tree where the number of nodes grows quadratically in the number of time steps. Fang and Oosterlee [11] use a Fourier-cosine series expansion approach for pricing Bermudan options under the Heston model. As for finite difference methods for solving the parabolic partial differential equation associated to the option pricing problems, they can be based on implicit, explicit or alternating direction implicit schemes. The implicit scheme requires to solve a sparse system at each time step. Clarke and Parrott [8] and Oosterlee [26] formulate the American put pricing problem as a linear complementarity problem (LCP) and use an implicit finite difference scheme combined with a multigrid procedure, whereas Forsyth, Vetzal and Zvan [33] use a penalty method. The explicit scheme is a quick approach although it requires small time steps to retain the stability. This request brings to a large number of time steps and is not economic in computation. The ADI schemes are good alternative methods. For example, Hout and Foulon [20] investigate four splitting schemes of the ADI type for solving the PDE Heston equation: the Douglas scheme, the Craigh-Sneyd scheme, the Modified Craigh-Sneyd scheme and the Hundsdorfer-Verwer scheme. Ikonen and Tovainen [19] propose a componentwise splitting method for pricing American options in the Heston model. The linear complementarity problem associated to the American option problem is decomposed into a sequence of five one-dimensional LCP's problems at each time step. The advantage is that LCP's need the use of tridiagonal matrices. In Haentjens, Hout and Foulon [20], the splitting method of Ikonen and Tovainen is combined with ADI schemes in order to obtain more efficient numerical results. In this paper we propose a new approach based both on tree and finite difference methods. Roughly speaking, our method approximates the CIR type volatility process through a tree approach already studied in Appolloni, Caramellino and Zanette [3], which turns out to be very robust and reliable. And at each step, we make use of a suitable transformation of the asset price process allowing one to take care of a new diffusion process with null correlation w.r.t. the volatility process. Then, by taking into account the conditional behavior with respect to the evolution of the volatility process, we consider a finite difference method to deal with the evolution of the (transformed) underlying asset price process. The paper is organized as follows. In Section 2, we introduce the model, we study in details the partial differential equation associated to the pricing problem (Section 2.1) and then we describe the hybrid-finite difference scheme (Section 2.2). The convergence of the approximating algorithm is studied in Section 3. Section 4 is devoted to numerical results and comparisons. the continuous dividend rate. We recall that the dynamics of V follows a CIR process with mean reversion rate κ and long run variance θ. The parameter σ is called the volatility of the volatility. From now on we setρ = 1 − ρ 2 and Z V = W, Z S = ρW +ρZ, in which (W, Z) denotes a standard 2-dimensional Brownian motion. So, the dynamics can be written as dS(t) S(t) = (r − δ)dt + V (t) (ρdW (t) +ρdZ(t)), (2.1) We consider the diffusion pair (Y, V ), where One has

4)
dV (t) = κ(θ − V (t))dt + σ V (t) dW (t), (2.5) (recall that W and Z are independent Brownian motions), with In the following, we define µ Y and µ V to be the drift coefficient of Y t and V t respectively, i.e.
This means that any functional of the pair (S t , V t ) can be written as a suitable functional of the pair (Y t , V t ) by using the transformation (2.3), so (Y t , V t ) will be our underlying process of interest.

The associated pricing PDE in a small time interval
Let f = f (y, v) be a function of the time and the space-variable pair (y, v). For h small, we need to compute (an estimate for) the quantity u(t, y, v) defined through in which (Y t,y,v , V t,v ) denotes the solution to (2.4) and (2.5) with the starting condition (Y t , V t ) = (y, v). In our mind, the time instant t plays the role of a dicretization instant in [0, T ], that is t = nh, so t + h = (n + 1)h stands for the next discretizing time.
We first notice that where F W t+h = σ(W u : u ≤ t + h). But conditionally to F W t+h , the volatility process V can be considered deterministic and the process Y turns out to have constant coefficients. More precisely, for g ∈ L 2 ([t, t + h]) with g ≥ 0 a.e. and g t = v, set Let us now discuss the quantityū(t, y; g), for g fixed as required above, that is g : [t, t + h] → R + is continuous and g t = v. (2.10) Set U (g),s,y t+h as the solution U (g) at time t + h, with starting condition U (g) s = y of the following stochastic differential equation with deterministic (although path dependent) coefficients: Recall that the associated infinitesimal generator is given by So, we getū (t, y; g) =ū(s, y; g)| s=t , withū(s, y; g) = E f (U (g),s,y t+h ; g t+h ) . Now, from the Feynman-Kac formula, the function (s, y) →ū(s, y; g) solves the parabolic PDE Cauchy problem ∂ sū (s, y; g) + L (g) sū (s, y; g) = 0 y ∈ R, s ∈ (t, t + h), u(t + h, y; g) = f (y, g t+h ) y ∈ R. (2.13) Once the problem (2.13) is solved, we can proceed to compute u(t, y, v) by using (2.9). We stress that the fixed path g plays the role of a parameter and the solution to (2.13) depends in general on the whole trajectory of g.
We consider now the case h ≃ 0, so that, by (2.10), This numerically brings to replace (2.13) with a PDE problem with constant coefficients. More precisely, we consider the approximationû h (s, y; v, g t+h ) forū(t, y; g) given by the solution to with yy . Let us remark that the solution to (2.14) actually depends on g only through v = g t (appearing in the coefficients of the second order operator) and g t+h (appearing in the Cauchy condition), that is why we used the notationû h (s, y; v, g t+h ). In contrast, the function solving (2.13) depends in principle on the whole trajectory g over the time interval [t, t + h]. Now, problem (2.14) can be easily solved by using a finite difference numerical method. Numerical reasonings suggest the use of an implicit approximation (in time) if v is "far enough" from zero, otherwise an explicit method should be considered -details are given in Section 3.1.1 and 3.1.2. This means that one fixes a space-step ∆y h and a space-grid Y h = {y j = Y 0 + j∆y h } j∈Z splitting the real line and approximates the solutionû h (s, y; v, g t+h ) to (2.14) on the grid Y h by means of a linear operator Now, recalling (2.9) and the fact thatū h ≃û h , on the grid Y the function u is approximated through We stress that the expectation on the r.h.s. above is now written in terms of the process V only, and this is the key point of our story because we can now use the tree method in [3]. But we will examine in depth this point in a moment. In practice, one cannot solve the PDE problem over the whole real line. So, one takes a positive integer M h > 0 such that M h ∆y h → +∞ as h → 0 and considers a discretization of h is finite and the approximation ofû h n (nh, y; v, g t+h ) is done by adding to (2.14) suitable boundary conditions. By calling again Π h (v) the matrix (now, finite dimensional) giving the solution from the finite difference approach, we still obtain

The hybrid tree-finite difference approach
We describe the main ideas of our approximating algorithm by means of an example, which was our starting point.
Consider an American option with maturity T and payoff function (Φ(S t )) t∈[0,T ] . First of all, by using (2.3) we replace the pair (S, V ) with the pair (Y, V ), so the obstacle will be given by The price at time 0 of such an option is then approximated by a backward dynamic programming algorithm, working as follows. First, consider a discretization of the time interval [0, T ] into N subintervals of length h = T /N: Then the price P (0, Y 0 , V 0 ) of such an American option is numerically approximated through the quantity P h (0, Y 0 , V 0 ) which is iteratively defined as follows: for (y, v) ∈ R × R + , P h (T, y, v) = Ψ(y, v) and as n = N − 1, . . . , 0 From the financial point of view, this means to allow the exercise at the fixed times nh, n = 0, . . . , N. Now, what we are going to set up is a contamination of a tree method for the process V with a finite difference method to handle the noise in Y (which is independent of the noise driving V ). In fact, the expectations appearing in the backward induction can be written as expectations of functions of the process V only, such functions being solution to parabolic PDE's. So, we proceed as described in the previous section: we fix a grid on the . . , M h }, and we approximate the above conditional expectations for y = y i ∈ Y h M h by using the matrix Π h (v) from the finite difference method. So, as already seen in (2.16), we write (2.17) By resuming, the price P (0, y, v) in y ∈ Y h M h can be numerically computed from the function P h (0, y, v) defined on the grid Y h M h as follows: (2.18) We stress that the backward induction (2.18) is now written in terms of the process V only, and here the binomial tree method in [3] comes on, we see how in a moment. First, let us briefly recall the binomial tree procedure in [3]. For n = 0, 1, . . . , N, consider the lattice (notice that v 0,0 = V 0 ) and for each fixed v n,k ∈ V h n , we define k h d (n, k) = max{k * : 0 ≤ k * ≤ k and v n,

20)
k h u (n, k) = min{k * : k + 1 ≤ k * ≤ n + 1 and v n, The transition probabilities are defined as follows: starting from the node (n, k) the probability that the process jumps to k h u (n, k) at time-step n + 1 is set as And of course, the jump to (n + 1, k h d (n, k)) happens with probability p h . This gives rise to a Markov chain (V h n ) n=0,...,N that weakly converges, as h → 0, to the diffusion process (V t ) t∈[0,T ] and turns out to be a robust tree approximation for the CIR process V . This means that we can approximate the expectation of suitable functionals of the diffusion V with the same expectation evaluated on the Markov chainV h . In particular, for a function g we write . So, at step n, we can numerically compute the expectation in the backward induction (2.18) on the lattice V h n as We can finally write the backward induction giving our approximating algorithm: for n = 0, 1, . . . , N, where the sum above is done for k * ∈ {k h u (n, k), k h d (n, k)} and j ∈ J M h . Notice that, at time step n, for every fixed i ∈ J M h and k = 0, . . . , n the sum in (2.23) can be seen as an integral w.r.t. the measure for every Borel sets A and B, δ {a} denoting the Dirac mass in a, so that µ h (y i , v n,k ; ·) is a discrete measure on Y h M h × V h n+1 . Now, in next Section 3 we shall be able to prove that, for small values of h, Π h (v) is a stochastic matrix. This gives that µ h (y i , v n,k ; ·) is actually a probability measure, that can be interpreted as a transition probability measure. Thus, we are doing expectations on a Markov chain (Ȳ h n ,V h n ) n=0,1,...,N , whose state-space, at time step n, is given by Y h M h × V h n and whose transition probability measure at time step n is given by µ h (y i , v n,k ; ·) in (2.24). Moreover, we shall prove that, under appropriate conditions on ∆y h and M h such that, as h → 0, ∆y h → 0 and M h ∆y h → ∞ (see next (3.21) and (3.22)), the family of Markov chains (Ȳ h ,V h ) h weakly converges to the diffusion process (Y, V ) . And this gives the convergence of our hybrid tree-finite difference algorithm approximating the Heston model. We shall finally discuss the convergence of the backward induction (2.23) to the price of the associated American option written on the Heston model.

The convergence of the algorithm
We first set up the finite difference method we take into account. Then, in Section 3.2, we formally define the approximating Markov chain and prove the weak convergence to the Heston model in the path space. 3.1 The finite difference scheme for the PDE problem (2.14) As described in Section 2.2, at each time step n we need to numerically solve (2.14) for t = t n = nh. So, we briefly describe the finite difference method we apply to problem (2.14), outlining some key properties of the associated operator allowing us to prove the convergence. For further information on finite difference methods for partial differential equations we refer for instance to [28]. Let t = nh, v and v * = g nh+h be fixed and let us set u n j =ũ h n (nh, y j ) the discrete solution of (2.14) at time nh on the point y j of the grid Y M -for simplicity of notations, we do not stress on u n j the dependence on v (from the coefficients of the PDE), v * (from the Cauchy problem) and h. It is well known that the behavior of the solution to (2.14) changes with respect the magnitude of the rate between the diffusion coefficient (ρ 2 v/2) and the advection term (µ Y (v)). To deal with these cases, we fix a small real threshold ǫ > 0 and in the following we shall describe how to solve both the case v < ǫ and v > ǫ by applying an explicit in time and an implicit in time approximation respectively. It is well known that for a big enough diffusion coefficient, to avoid over-restrictive conditions on the grid steps, it is suggested to apply implicit finite differences to problem (2.14). In this case, the discrete solution {u n j } j∈Z at time nh will then be computed in terms of the solution {u n+1 j } j∈Z at time (n + 1)h by solving the following discrete problem: where ∆y = y j − y j−1 , ∀j ∈ Z.
On the other hand, when the diffusion coefficient is small compared with the reaction one, it is suggested to apply an explicit in time approximation coupled with a forward or backward finite differences for the first order term u x depending on the sign of the reaction coefficient.
Specifically, for v close to 0, that is v < ǫ, we solve the problem by the following approximation schemes: As previously mentioned at the end of the Section 2.1, for the numerical tests one has to deal with a finite grid Y h M = {y j } j∈J M (for simplicity, here we simply write M instead of M h ) and problems (3.1) and (3.2) have to be coupled with suitable numerical boundary conditions. Here, we assume that the two-off the domain values are defined by the following relations (Neumann-type boundary conditions): in the implicit case we set whereas in the explicit case we set Other conditions can surely be selected, for example the two boundary values u n+1 −M and u n+1 M may be a priori fixed by a known constant (this typically appears in financial problems). All arguments that follow apply as well. So, hereafter we set By applying implicit finite differences (3.1) coupled with boundary conditions (3.4), we get the solution u n = (u n −M , . . . , u n M ) T by solving the following linear system where A is the (2M + 1) × (2M + 1) tridiagonal real matrix given by We immediately prove that the solution u n to (3.7) actually exists at least when β > |α| (we will see later that this is not a restrictive condition).
In next Section 3.2 we need to deal with A −1 ψ i l (y) for l ≤ 4 and i ∈ J M . So, we study such objects. By Proposition 3.1, for β > |α| one has that A is invertible and we may then compute In the following, the symbol [·] will stand for the floor function and we use the understanding that is a n = 2β if n is even and a n = 2α if n is odd, and the coefficients b ±M can be bounded as follows: On the other hand, following the same reasonings for k = ±M we easily obtain and by using (3.12) we get the full form (3.10). Finally, the estimate in (3.13) immediately follows from the Newton's binomial formula.
We are now ready to characterize A −1 ψ i l (y), for every polynomial ψ i l (y) as in (3.9).
Proposition 3.3 Suppose that β > |α| and for l ≥ 1 let γ l,k , k = 0, 1, . . . , l − 1, be iteratively (backwardly) defined as follows: where the coefficients a n are given in (3.11). Then the b ±M j,i 's being given in (3.12). Moreover, if ∆y ≤ 1, M∆y ≥ 1 and l2 l+2 (β∆y 2 + |α|∆y) ≤ 1, the following estimate holds for T ±M l,i : for every i ∈ J M , We call ( * ) the r.h.s. above. By using Lemma 3.2 and setting By the definition of the γ j,k 's and φ i l,M (y)'s, each coefficients in the above (first) sum is null and the sum of the last three terms is zero, so that ( * ) = ψ i l (y). Let us discuss (3.16). By using (3.13) and the fact that |y ±M − y i | ≤ 2M∆y, since ∆y(1 + 2M) ≥ 1 we can write It remains to prove that, under our constraints, l−1 j=0 |γ l,j | ≤ 1. For every k = 0, 1, . . . , l − 1 and j = k + 1, . . . , l − 1 we consider the estimates By inserting in (3.14) we can write We also notice that, for m ≥ 1 and ∆y ≤ 1, |a m |∆y m ≤ 2(β∆y 2 + |α|∆y) so that we get and by summing over k, from which it follows that l−1 k=0 |γ l,k | ≤ 1 and the statement holds.
Remark 3.4 For further use, we write down explicitly the vector A −1 ψ i l (y) for l = 1, 2, 4. Straightforward computations give the following: In particular, since (ψ i l (y)) i = 0 for every l ≥ 1, the ith entry of the above sequences are given by By replacing the α and β expressions (3.6), we get the formulas Furthermore, to deal with the numerical boundary conditions, as those given in (3.4), we need to study the behavior of the ith component of the boundary term φ i l,M (y) in (3.15) as i is "far from the boundary" and l = 1, 2, 4. Here, we use a quite general result (allowing one to set up different boundary conditions) whose precise statement and proof are postponed in Appendix A.1.

The case v < ǫ
Here we go through the already described procedure for the explicit in time approximation. We recall that, for v < ǫ, we consider the forward finite differences (3.2) or the backward finite differences (3.3) for the first order term depending on the sign of the reaction coefficient: µ Y (v) > 0 or µ Y (v) < 0 respectively, and from (3.6) this means that α > 0 or α < 0 respectively. So, by considering also the case α = 0 and by coupling with the boundary conditions (3.5), straightforward computations give that the solution u n of the explicit in time scheme satisfies the condition

18)
α and β being given in (3.6) and 1l denoting the indicator function (1l Γ = 1 on Γ and 1l Γ = 0 otherwise). We remark that C is a stochastic matrix if and only if 2β + 2|α| ≤ 1. (3.19) We also notice that if ǫ is small enough such a condition is not restrictive, but we will discuss deeper this point later. In next Section 3.2 we need to deal with Cψ i l (y) for l ≤ 4 and i ∈ J M , where the function ψ i l (y) ∈ R 2M +1 are defined in (3.9). So, we get Lemma 3.6 Let ψ i l (y) be defined in (3.9). Then for every l ∈ N and i ∈ J M one has that is d n = 2(β + |α|) if n is even and d n = 2α if n is odd, and c ±M l,i are given by where A + and A − denotes the set of the even and the odd numbers respectively.
The proof is straightforward (it suffices to follow the same arguments developed for Lemma 3.2), so we omit it.
Remark 3.7 In the special case l = 1, 2, 4, Lemma 3.6 gives In particular, since (ψ i l (y)) i = 0 for every l ≥ 1 and assuming |i| < M, the ith entries of the above sequences are given by (3.20) in which we have inserted the formulas for α and β in (3.6),

The associated Markov chain and the convergence of the hybrid algorithm
We denote, as in Section 2.2, (V h n ) n=0,1,...,N , with Nh = T , the Markov chain approximating the volatility process V introduced in [3]. We recall that the state-space at step n is given by V h n defined in (2.19). We define now the Y -component of our Markov chain. First, given the time-step h, we set up the dependence on h for the space-step ∆y, the number M giving the points of the grid Y M = {y i = Y 0 + i∆y ; i = −M, . . . , M} and the threshold ǫ that allows us to use the explicit or the implicit finite difference method. So, we assume that where c M > 0 and the constants c y , c ǫ , p, q > 0 are chosen as follows in which the parameters κ, θ, r, δ, σ and ρ come from our model, see Section 2. Let us stress that the last constraint in (3.22) can be really satisfied, for example by choosing c y > 4|r − δ − ρ σ κθ|. We also notice that (3.21) and (3.22) give for every h small enough. Now, Proposition 3.8 ensures that there exists h 0 > 0 such that for every h < h 0 and for every v ∈ ∪ N n=0 V h n , the matrix A −1 discussed in Section 3.1.1 and the matrix C discussed in Section 3.1.2 are both well defined and stochastic matrices. So, for h small, that is h < h 0 with h 0 as in Proposition 3.8, we define Π h (v) as follows: As a consequence, we can assert that for every h is a stochastic matrix. We now define the following transition probability law: at time-step n ∈ {0, 1, . . . , N}, for (y i , v n,k ) ∈ Y h M h × V h n we set µ h (y i , v n,k ; ·) the law in R 2 as in (2.24), that is So, we callX h = (X h n ) n=0,1,...,N the 2-dimensional Markov chain having (3.23) as its transition probability law at time-step n ∈ {0, 1, . . . , N}, that is Since j Π(v) ij = 1, one gets that the second component ofX h is a Markov chain itself and it equals, in law, toV h . So, we writē X h = (Ȳ h ,V h ) and for every function f : (3.24) So, coming back to the discussion in Section 2.2, by (3.24) we can assert that our algorithm is actually given by approximating the diffusion pair X = (Y, V ) with the Markov chainX h = (Ȳ h ,V h ). So, we set X h = (Y h , V h ) as the piecewise constant and càdlàg interpolation in time ofX h , that is X h t =X h n , t ∈ [nh, (n + 1)h), n = 0, 1, . . . , N − 1. (3.25) We set D([0, T ]; R 2 ) the space of the R 2 -valued and càdlàg functions on the interval [0, T ], that we assume to be endowed with the Skorohod topology (see e.g. Billingsley [5]). Our main result is the following: Proof. The idea of the proof is standard, see e.g. Nelson and Ramaswamy [24] or also classical books such as Billingsley [5], Ethier and Kurtz [10] or Stroock and Varadhan [29].
To simplify the notations, let us set It is clear that M Y n,i,k (h; l) is the local moment of order l at time nh related to Y , as well as M V n,k (h; l) is related to the component V , and M Y,V n,i,k (h) is the local cross-moment of the two components at the generic time step n. So, the proof reduces in checking that, for fixed r * , v * > 0 and setting Λ * = {(n, i, k) : v n,k ≤ v * , |y i | ≤ r * }, then the following properties i), ii) and iii) hold: iii) (fast convergence to 0 of the fourth order local moments) We recall that the V -component of the 2-dimensional Markov chain is a Markov chain itself and we have The limits in i), ii) and iii) containing M V n,k (h; l) for l = 1, 2, 4, have been already proved in [3] (see Theorem 7 therein), so we prove the validity of the remaining limits. We set ψ i l (y) the vector in R 2M h +1 given by (ψ i l (y)) j = (y j − y i ) l , j ∈ J M h . From (3.24) we get and we notice that the above quantity has been already discussed in the previous sections. We set Λ * = Λ * ,1,h ∪ Λ * ,2,h , with For (n, i, k) ∈ Λ * ,1,h , Π h (v n,k ) is the inverse of the matrix A in (3.8). So, by using (3.17), we have being given in (3.15). In Lemma 3.10 below, we prove that for l ≤ 4, And since (µ Y ) n,k is bounded on Λ * , the limits in i), ii) and iii) associated to M Y n,i,k (h; l), l = 1, 2, 4, hold uniformly in Λ * ,1,h . We prove the same on the set Λ * ,2,h . For (n, i, k) ∈ Λ * ,2,h , the matrix Π h (v n,k ) to be taken into account is given by the matrix C in (3.18). Moreover, for h small enough, we have that if (n, i, k) ∈ Λ * ,2,h then |i| < M h . Therefore, by (3.20) we obtain and again the limits in i), ii) and iii) concerning M Y n,i,k (h; l), l = 1, 2, 4, hold uniformly in Λ * ,2,h . It remains to study the cross-moment. By using (3.24), it is given by and the convergence as in ii) immediately follows from the already proved limits in i).
In order to conclude, we only need to prove next Lemma 3.10 Assume that (3.21) and (3.22) both hold. Let v * , r * > 0 and set Proof. We use Proposition 3.5. Here, C denotes a positive constant that may vary from line to line. Under (3.21) and (3.22), for (n, i, k) ∈ Λ * ,1,h we have already observed that β h > |α h |, α h , β h being given in (3.6), and the constraints ∆y h ≤ 1 and M h ∆y h ≥ 1 are trivially satisfied for h small. Moreover, on the set Λ * ,1,h , there exists C > 0 such that for every l ≤ 4, for h small enough. And since β h + |α h | ≤ Ch 1−2p , we also have Then, by applying Proposition 3.5, we can write Now, on the set Λ * ,1,h we have |y i | ≤ r * , so that |i| ≤ r * +Y 0 ∆y h ≤ Ch −p . And by recalling that q > p, we get M h − i ≥ C(h −q − h −p ) ≥ ch −q for some c > 0. So, we have proved that there exist C, c > 0 such that for every h close to 0 Now, straightforward computations give that, as h → 0, the r.h.s. above tends to 0 if q > p and p ≤ 1, and this completes the proof.
Remark 3.11 Theorem 3.9 proves the convergence of the algorithm in the case we introduce suitable boundary conditions in the finite difference component of our procedure. We stress that our Neumann-type conditions (3.4) and (3.5) may be replaced by other types of boundary conditions, that can be handled by using Proposition A.1 (see Appendix A.1). We also recall that the use of a boundary is a numerical requirement which is necessary to solve the problem in practice. However, one could prove the convergence result even in the whole grid Y h = {y i = Y 0 + i∆y h } i∈Z , that is by considering the solutions to (3.1) and (3.2) without linking to these equations any boundary condition. In fact, under the requirements of Theorem 3.9, the inverse of the implicit difference linear operator A associated to (3.1) remains a Markovian transition function, as well as the explicit difference linear operator C related to (3.2). Moreover, formulas similar to (3.17) and (3.20) hold, with a simplification from the fact that boundary contributions do not exist. For the sake of completeness, we give these proofs in Appendix A.2. This means that the arguments used to prove Theorem 3.9 can be applied also on the infinite grid Thanks to Theorem 3.9 we can deal with the convergence of the price evaluated on our Markov chain to the one computed on the Heston model. The problem of pricing European options is simple when the payoff is not too complicated. In fact, let T denote the maturity date and f be the payoff function, that is f : D([0, T ]; R) → R + . First, we write the payoff in terms of the transformation (2.3), so we get the transformed payoff-function g(y, v) = f (e y+ ρ σ v ), (y, v) ∈ D([0, T ]; R 2 ), and the associated option prices on the continuous and the discrete model as seen at time 0 are given by and respectively,g denoting the discounted payoff, i.e.g = e −rT g (it is clear that the writing E for both prices is an abuse of notations, since in principle one should use the notations E P and E P h related to the measures P and P h of the probability space where the processes (Y, V ) and (Y h , V h ) are defined). Now, the weak convergence in Theorem 3.9 ensures the convergence P h Eu → P Eu of the European price when the discounted payoff-function fulfills the following requests: (y, v) →g(y, v) is continuous and there exists a > 0 and h * > 0 such that This a consequence of standard results on the convergence of the expectations for sequences of random variables which are weakly convergent and satisfy uniform integrability properties. As for American style options, even for simple payoffs things are more difficult because of the presence of optimal stopping times. However, due to the results in Amin and Khanna [2], we can deduce the convergence of the prices for suitable payoffs. In fact, let f : [0, T ] × D([0, T ]; R) → R + denote a payoff function. By passing to the pair (Y, V ) and by considering the resulting discounted payoff functiong(t, y, v) = e −rt f (t, e y+ ρ σ v ), the associated option prices on the continuous and the discrete model as seen at time 0 are given by Let ρ h denote the optimal stopping time related to the discrete problem, that is Then by using the arguments in Amin and Khanna [2], properties i)-iii) in the proof of Theorem 3.9 ensure the existence of ρ ∈ G 0,T such that the triple (ρ h , Y h , V h ) converges in law to (ρ, Y, V ) -this is a very roughly speaking: in the words by Amin and Khanna, ρ is "in some appropriate sense a legitimate stopping time w.r.t. the filtration generated by (Y, V )", and to be precise one should refer to a further probability space (for technical details, see the discussion at pages 299-300 and in particular Theorem 4.1 in [2]). We stress that in [2] this result is based on two starting assumptions: the (global) Lipschitz continuity and the sublinearity property for both the drift and the diffusion coefficient of the pair (Y, V ). Here, the Lipschitz continuity property does not hold because of the presence of the square-root function, so the diffusion coefficient is only Hölder continuous. Nevertheless, they use such a property just to ensure the existence of the strong solution, a condition that here holds, and use specifically the sublinearity property of the drift coefficients, that here holds. So, their arguments apply as well. Once the weak convergence of the triple is achieved, the plan to prove the convergence of the prices is the following. Under (H1) one gets that {g(ρ h , Y h , V h )} h converges in law, as h → 0, tõ g(ρ, Y, V ). Moreover, (H2) implies that the {g(ρ h , Y h , V h )} h is a uniformly integrable family of random variables, and this suffices to get the convergence of the expectations. Finally, one gets P Am = E(g(ρ, Y, V )), and the convergence P h Am → P Am of the American prices is achieved. As an immediate consequence, both European and American put options can be numerically evaluated by means of the approximating algorithm described in Section 2.2 (and, due to the call-put parity formula, also European call option can be numerically priced with our method). However, in next Section 4 we apply our procedure to numerically price European/American barrier options written on the Heston model, that is options which are much more sophisticated and, for example, do not fulfill (H1). Nevertheless, the numerical experiments give a real evidence of the goodness (that is, convergence) of our algorithm also in this case.

Numerical results
In this section we provide numerical results in order to asses the efficiency and the robustness of our hybrid tree-finite difference method in the case of plain vanilla options and barrier options.

European and American vanilla options
We compare here the performance of the hybrid tree-finite difference algorithm (called HTFD) introduced in Section 2.2 with the tree method of Vellekoop Nieuwenhuis (called VN) in [25] for the computation of European and American options in the Heston model. All the numerical results (except for our method) are obtained by using the Premia software [27]. In the European and American option contracts we are dealing with, we consider the following set of parameters: initial price S 0 = 100, strike price K = 100, maturity T = 1, interest rate r = log(1.1), dividend rate δ = 0, initial volatility V 0 = 0.1, long-mean θ = 0.1, speed of mean-reversion κ = 2, correlation ρ = −0.5. In order to study the numerical robustness of the algorithms, we choose three different values for σ: we set σ = 0.04, 0.5, 1. We first consider the case σ = 0.04, that is σ close to zero (which implies that the Heston PDE is convection-dominated in the V -direction). Moreover, for σ = 1, we stress that the Feller condition 2κθ ≥ σ 2 is not satisfied. In the (pure) tree method VN, we fix the number of points in the V coordinate as N V = 50, with varying number of time and space steps: N t = N S = 50, 100, 200, 400. The numerical study of the hybrid tree-finite difference method HTFD is split in two cases: HTFD1 refers to the (fixed) number of time steps N t = 100 and varying number of space steps N S = 50, 100, 200, 400; we add the situation HTFD2 where the number of time steps is equal to the number of space steps N t = N S = 50, 100, 200, 400. Table 1 reports European put option prices. Comparisons are given with a benchmark value obtained using the Heston closed formula CF in [17]. In Table 2 we provide results for American put option prices. In this case we use a benchmark from the Monte Carlo Longstaff-Schwartz algorithm, called MC-LS, as in [23], with a huge number of Monte Carlo simulation (1 million iterations) which are done by means of the accurate Alfonsi [1] discretization scheme for the CIR process with M = 100 discretization time steps and bermudan exercise dates. We recall that the Alfonsi method provides a Monte Carlo weak second-order scheme for the CIR process, without any restriction on its parameters.  Table 1: Prices of European put options. σ = 0.04, 0.5, 1. S 0 = 100, K = 100, T = 1, r = log(1.1),  Table 2: Prices of American put options. σ = 0.04, 0.5, 1. S 0 = 100, K = 100, T = 1, r = log(1.1), The numerical results show that the hybrid tree-finite difference method is very accurate, reliable and efficient.

European and American barrier options
We study here the continuously monitered barrier options case and we compare our hybrid tree-finite difference algorithm with the numerical results of the method of lines provided in Chiarella et al. [7]. So, we consider European and American up-and-out call options with the following set of parameters: K = 100, T = 0.5, r = 0.03, δ = 0.05, V 0 = 0.1, θ = 0.1, κ = 2, ρ = −0.5. The up barrier is H = 130. We choose different values for S 0 : S 0 = 80, 100, 120.
We also compare with a benchmark value obtained by using the method of lines, called MOL, with mesh parameters 100, 200, 6400 (see Chiarella et al. [7]).  Table 3: Computational times (in seconds) for European put options for σ = 0.5. Table 4 and Table 6 report European and American Up-and-Out option prices respectively, while

A Appendix
A.1 Boundary sensitivity for the implicit finite difference operator We study here the behavior of the solution x = (x 1 , . . . , x N ) T of the two following linear systems where v i , i = 1, . . . , N, denotes the standard orthonormal basis in R N , i.e. (v i ) k = 0 for k = i and (v i ) i = 1, i = 1, . . . , N, and where A has the following general tridiagonal form The result we are going to present is due for matrices A as in (A.3) and that satisfy the hypotheses (P2)-(P3) in the proof of Proposition 3.1, ensuring that they are invertible Mmatrices (see for instance [4]).
Assume moreover the following stability conditions on the "boundary" values a 1 , a N , c 1 and b N : where z + = (a + a 2 − 4|bc|)/2. Then the solution x of (A.1) is defined by a sequence {x k } k=1,...,N of positive terms and there exists a positive valueγ * > |b| such that, for k = 2, . . . , N − 1 Similarly, for the solution x of (A.2) it holds x k > 0, for all k = 1, . . . , N and there exists a positive value γ * > |c| such that for k = 2, . . . , N − 1, Proof. Let us start by estimating first the solution x of system (A.2). By applying the Thomas algorithm, also known as tridiagonal matrix algorithm [30], the solution is given by back substitutions: where the coefficients γ k are recursively defined by It is easy to verify that under assumptions (A.4), for k = 3, . . . , N − 1, the sequence {γ k } has two strictly positive fixed points z ± = (a ± a 2 − 4|bc|)/2. Moreover, z − is an unstable fixed point while z + is stable. Since γ 2 = a−|bc 1 |/a 1 and by (A.5), we have that γ 2 > z − . So, starting from γ 2 the sequence converges to z + and we have that for γ * = min Furthermore, by (A.4) we have that γ N > 0. Going back to the sequence {x k } i=1,...,N , we first notice that since γ N > 0 then x N > 0 and accordingly x k > 0 for all k = 1, ..., N. Moreover, from condition (A.9) we obtain (A.7). In fact, for k = 2, . . . , N − 1 we have and thus It remains to verify that γ * > |c|. We have already seen that γ 2 > |c|; so, it remains to verify the case γ * = z + which is obtained as a result of the strictly diagonal dominant condition on A.
To obtain the estimate (A.6) we introduce the N × N matrix U satisfying Uv i = v N +1−i , i = 1, . . . , N, so Since Uv N = v 1 and UU = I (i.e. U −1 = U), to compute (A.6) we use that So, following the same reasoning as above, we getγ 2 = a − |b N c| a N andγ * = min (γ 2 , z + ) > |b| such thatx In fact, as for the second inequality, (A.7) gives |( Similarly, one hasγ N ≥ 1, and the first inequality holds as well.

A.2 The finite difference operators on the infinite grid
For fixed ∆y > 0 and y 0 ∈ R, we consider the grid Y = {y k } k∈Z on R by setting Any real function ϕ on Y can be seen as a sequence on R. So, as usual we set R Z = {x = (x k ) k∈Z : x k ∈ R ∀k},and we think to x ∈ R Z as "column-type" vectors. We also set ϕ(y), y ∈ Y, as the point in R Z defined as (ϕ(y)) k = ϕ(y k ), k ∈ Z.
We split our discussion in two parts: the case v > ǫ, giving the implicit finite difference operator, and v ≤ ǫ, related to the explicit finite difference operator

A.2.1 The case v > ǫ
For the implicit in time approximation, equation (3.1) can be written as u n = Au n+1 , in which u n = (u n j ) j∈Z is the unknown, u n+1 = (u n+1 j ) j∈Z is given and A = (A kj ) k,j∈Z is the infinite dimensional matrix given by and A kj = 0 for |j − k| > 1.
A is actually a linear operator on R Z : by using the standard row/column product, we get A : R Z → R Z and y = Ax if and only if where I denotes the identity map over R Z and T : R Z → R Z is the linear operator defined through the infinite dimensional matrix Throughout this section, we assume β > |α| ≥ 0.
This gives that T ≥ 0, in the sense that T ij ≥ 0 for every i, j ∈ Z. First, we want to show that A is an invertible operator, at least over a nice subset of R Z . Secondly, we want to show that the inverse is a stochastic infinite dimensional matrix, that is it defines a Markov transition function over Z. So, for a > 0, we set ℓ 2 a = x ∈ R Z : k∈Z a 2|k| |x k | 2 < ∞ .
Remark A.3 If a < 1 then ℓ 2 a contains all the sequences with polynomial growth: for any fixed l ∈ N, the sequence x = (x k ) k such that |x k | ≤ C(1 + |k| l ), k ∈ Z, belongs to ℓ 2 a . Moreover, ℓ 2 a 1 ⊂ ℓ 2 a 2 if a 1 ≤ a 2 ≤ 1. Finally, we notice that, for a < 1, ℓ 2 a contains also special sequences with exponential growth, namely points x ∈ R Z such that |x k | ≤ C(1 + e γ|k| ) for every γ < −2 log a.
On ℓ 2 a , we define the inner product and the norm by x, y a = k∈Z a 2|k| x k y k and |x| 2 a = k∈Z a 2|k| x 2 k respectively. When a = 1 one has which is the standard ell-2 space of sequences. For a > 0 the map D a : ℓ 2 a → ℓ 2 , (D a x) k = a |k| x k , k ∈ Z, (A.10) is well defined and invertible, with D −1 a = D a −1 . Moreover, since x, y a = D a x, D a y , x, y ∈ ℓ 2 a , (A.11) ·, · denoting the standard scalar product in ℓ 2 , D a is an isometry between ℓ 2 a and ℓ 2 . Therefore, the Hilbert structure of ℓ 2 can be transferred to ℓ 2 a through D a , so that the space ℓ 2 a , endowed with the scalar product ·, · a , is a Hilbert space. Since the objects of our interest will belong to ℓ 2 a with a < 1, we study the invertibility of A when it is restricted on ℓ 2 a .
Then A : ℓ 2 a → ℓ 2 a . Moreover, A is invertible and the inverse A −1 can be identified with an infinite dimensional matrix (Π i,j ) i,j∈Z that does not depend on the choice of a: for every x ∈ ℓ 2 a , Finally, Π defines a Markov transition function, that is Π kj ≥ 0 and j∈Z Π kj = 1 for every k, j ∈ Z.
Proof. For x ∈ ℓ 2 a , we set (I ± x) k = x k±1 , k ∈ Z. So, we have Straightforward computations give |I − x| 2 a = 1 a 2 k≤0 a 2|k| x 2 k + a 2 k≥1 a 2|k| x 2 k , so that |I − x| a ≤ 1 a |x| a . And similarly, one gets |I + x| a ≤ 1 a |x| a . Therefore, |T x| a ≤ (β − α)|I − x| a + (β + α)|I + x| a ≤ 2β a |x| a and this gives T : ℓ 2 a → ℓ 2 a , so that A : ℓ 2 a → ℓ 2 a . Moreover, one has T a ≤ 2β a , · a denoting the standard operator norm on ℓ 2 a . As a consequence one gets T n a ≤ 2β a n for every n ≥ 0.
We prove now that A is invertible on ℓ 2 a . For N ∈ N, we set S N : ℓ 2 a → ℓ 2 a through S N = 1 1 + 2β N n=0 1 (1 + 2β) n T n .

We notice that
Since 2β a(1+2β) < 1, the series converges and the limit (linear) operator S : ℓ 2 a → ℓ 2 a exists. Moreover, we can write Let us prove that S = A −1 . For every N, one has S N A = S N (1 + 2β)I − T = I − 1 (1 + 2β) n+1 T n+1 and by taking the limit as N → ∞ one gets SA = I. Similarly, one proves that AS = I, so that S = A −1 on ℓ 2 a . Now, let {e i } i∈Z ⊂ ℓ 2 denote the standard orthonormal basis: (e i ) k = 0 for k = i and (e i ) i = 1. Since ℓ 2 = ℓ 2 a with a = 1 > 2β 1+2β , we can set Π kj = e k , Se j , k, j ∈ Z.
We prove now that Π is actually the matrix we are looking for: if x ∈ ℓ 2 a then (Sx) k = j∈Z Π kj x j , k ∈ Z. Let D a denote the linear operator in (A.10). Since D −1 a = D a −1 , (A.11) immediately gives that {D a −1 e i } i is an orthonormal basis in ℓ 2 a . Moreover, for x ∈ ℓ 2 a one has x = j∈Z x j e j = j∈Z a |j| x j a −|j| e j = j∈Z (D a x) j D a −1 e j and since Sx ∈ ℓ 2 a , we can also write D a −1 e k , Sx a = a |k| (Sx) k .
So, (Sx) k = a −|k| D a −1 e k , Sx a = a −|k| j∈Z (D a x) j D a −1 e k , SD a −1 e j a = a −|k| j∈Z a |j| x j D a −1 e k , SD a −1 e j a = j∈Z a |j|−|k| D a −1 e k , SD a −1 e j a x j .
Therefore, we get (Sx) k = j Π kj x j with Π kj = a |j|−|k| D a −1 e k , SD a −1 e j a .
We show now that actually Π is independent of a and moreover Π kj = e k , Ae j . Since D a −1 e j = a −|j| e j and e j ∈ ℓ 2 a , using (A.11) and the fact that D a x ∈ ℓ 2 if x ∈ ℓ 2 a , we can write Π kj = a −|k| D a −1 e k , Se j a = a −|k| D a −1 e k , D a −1 D a Se j a = a −|k| e k , D a Se j = a −|k| (D a Se j ) k = a −|k| a |k| (Se j ) k = e k , Se j .
In order to show that Π is a Markov transition function on Z, we proceed as follows. One has Π kj = e k , Se j = lim N →∞ e k , S N e j = lim N →∞ 1 1 + 2β N n=0 1 (1 + 2β) n (T n ) i,j ≥ 0 because the entries of T are all ≥ 0, and so it is for T n , for every n. Finally, let 1 ∈ R Z be the unit vector, that is (1) k = 1 for every k ∈ Z. Notice that 1 ∈ ℓ 2 a for every a < 1 and T 1 = 2β1. Then S1 is well defined and S1 = 1 1 + 2β ∞ n=0 2β 1 + 2β n 1 = 1 and this automatically gives j∈Z Π kj = 1 for every k ∈ Z.
For each l ∈ N, we consider the polynomial (y − y i ) l and we call ψ i l (y) ∈ R Z the associated (vector) function of y ∈ Y: ψ i l (y) k = (y k − y i ) l = ∆y l (k − i) l , k ∈ Z. (A.12) By Remark A.3, one has ψ i l (y) ∈ ℓ 2 a for every fixed a < 1, i ∈ Z and l ∈ N, so the quantity A −1 ψ i l (y) makes sense. We also notice that ψ i 0 (y) = 1, so that A −1 ψ i 0 (y) = A −1 1 = 1. As already developed in Section 3.1.1, we need to deal with A −1 ψ i l (y) for l ≤ 4 and i ∈ Z. First of all we have where a n = (β − α)(−1) n + (β + α), n ∈ N, (A.13) that is a n = 2β if n is even and a n = 2α if n is odd.
The proof is identical to the one of Proposition 3.2, so we omit it. We only notice that, due to the use of the grid on the whole real line, we do not have any contribution from the boundary -see the terms b ±M l,i e ±M in (3.10). This gives that the inverse of A on any polynomial ψ i l (y) is much simpler, and in fact we have Proposition A.6 For l ≥ 1 let γ l,k , k = 0, 1, . . . , l, be iteratively (backwardly) defined as follows: γ l,k = l k a l−k ∆y l−k + l−1 j=k+1 γ l,j j k a j−k ∆y j−k , k = l − 1, . . . , 0, where the coefficients a n are given in (A.13). Then A −1 ψ i l (y) = ψ i l (y) + l−1 j=0 γ l,j ψ i j (y).
For the proof, we again refer to the finite case, that is Proposition 3.3. As a consequence, for l = 1, 2, 4 and i ∈ Z, the formulas for (A −1 ψ i l (y)) i are identical to the ones for the finite grid (see (3.17)) except for the boundary terms, which are null here. In fact, by inserting the formula (3.6) for β and α, we get (A.14) for every i ∈ Z.

A.2.2 The case v ≤ ǫ
The operator associated to the the explicit in time approximation is given by and C kj = 0 for |j − k| > 1.
α and β being given in (3.6), and the solution u n ∈ R Z is given by u n = Cu n+1 . The infinite dimensional matrix C defines a Markov transition function if and only if 2β + 2|α| ≤ 1.
Remark A.7 By using (A.14) and (A.15), we can prove the convergence result in Theorem 3.9 also in the case of the infinite grid. The proof is identical, even simpler because we do not have here contributions from the boundary.