Enhancing finite difference approximations for double barrier options: mesh optimization and repeated Richardson extrapolation

We show that the performances of the finite difference method for double barrier option pricing can be strongly enhanced by applying both a repeated Richardson extrapolation technique and a mesh optimization procedure. In particular, first we construct a space mesh that is uniform and aligned with the discontinuity points of the solution being sought. This is accomplished by means of a suitable transformation of coordinates, which involves some parameters that are implicitly defined and whose existence and uniqueness is theoretically established. Then, a finite difference scheme employing repeated Richardson extrapolation in both space and time is developed. The overall approach exhibits high efficacy: barrier option prices can be computed with accuracy close to the machine precision in less than one second. The numerical simulations also reveal that the improvement over existing methods is due to the combination of the mesh optimization and the repeated Richardson extrapolation.


Introduction
The most common approach for pricing double barrier options is based on solving a partial differential equation of Black-Scholes type (see Wilmott 1998). However, such an equation does not have a closed-form solution and thus requires numerical approximation. To this aim, some scholars (see, for example, Boyle and Lau 1994;Cheuk and Vorst 1996;Ritchken 1996;Ahn and Gao 1999;Boyle and Tian 1999;Figlewski and Gao 1999;Tian 1999) have proposed the use of binomial and trinomial B Luca Vincenzo Ballestra luca.ballestra@unibo.it 1 Department of Statistical Sciences "Paolo Fortunati", University of Bologna, Via delle Belle Arti 41, 40126 Bologna, Italy lattices, according to which the time variable and the price of the underlying asset (hereafter referred to as the space variable) are discretized on a suitable tree of nodes. This approach, albeit relatively simple to implement, has the disadvantage of not being very flexible, since the number of time-steps and the number of space discretization nodes in the tree cannot be chosen independently. As a consequence, binomial and trinomial lattices have progressively left the scene to a more modern and versatile numerical technique, namely the finite difference method, which has been employed, for instance, by Zvan et al. (2000), Duffy (2006), Wade et al. (2007), Ndogmo and Ntwiga (2011), and Milev and Tagliani (2013). In the present paper, we show that the performances of the finite difference method for double barrier option pricing can be considerably improved by applying a suitable mesh optimization and a Richardson extrapolation technique.
When pricing double barrier options by finite difference approximation (as well as by binomial/trinomial trees), some difficulties arise. In fact, the solution to be computed is not smooth at the strike price and at the barriers, and thus the space mesh, i.e. the set of nodes that are used to discretize the underlying asset price, must be chosen appropriately. First of all, to avoid losses of accuracy, a node of the mesh must be located exactly at the strike price, and other two nodes must be located exactly at the barriers (see Tian 1999;Pooley et al. 2003;Duffy 2006, Chapter 13,). Furthermore, if we want to improve the accuracy of the computed solution by applying some Richardson extrapolation procedure (see, for example, Arciniega and Allen 2004;Feng and Linetsky 2008;Chan et al. 2009;Ballestra 2014), we need a smooth pattern of error reduction, which requires us to use a mesh with equally spaced nodes. It is also worth noticing that, when using finite difference methods, a uniform mesh does also guarantee the highest rates of truncation error reduction (see, for example, Hirsch 1988, Chapter 4;Hyman et al. 2000;Duffy 2004). In the following, for the sake of brevity, a space mesh that satisfies the above two requirements will be said aligned and uniform. Now, an aligned and uniform mesh is straightforward to obtain only if the strike price and the barriers are the terms of some arithmetic progression. This is the case that is usually considered when proposing finite difference schemes for double barrier option pricing, see, e.g., Duffy (2004), Wade et al. (2007), Ndogmo and Ntwiga (2011), and Milev and Tagliani (2013). Nevertheless, it can also happen that the strike price and the barriers are not in arithmetic progression. Anyway, even if they are, in the frequently encountered case where the barriers are discretely monitored, another issue arises. In fact, when pricing double barrier options with discretely monitored barriers, we have to discretize not only the space interval between the two barriers, but also the space interval below the lower barrier and the space interval above the upper barrier (by contrast, in the case of continuously monitored barriers, only the space interval between the barriers needs to be considered). Then, clearly, if we want to use an aligned and uniform mesh, the number of space discretization nodes below the lower barrier, denote it N L , the number of space discretization nodes between the lower and the upper barriers, denote it N M , and the number of space discretization nodes above the upper barriers, denote it N U , cannot be chosen arbitrarily. In the following, a mesh that is not only aligned and uniform but also such that N L , N M and N U are freely and a-priori decided by the user will be called optimal. As we may easily understand (and as will be clear in Sect. 3) in the general case an optimal mesh is impossible to have, at least in the original space domain of the partial differential problem that we have to solve.
In the present paper we propose a definitive and also very effective remedy to this issue. Specifically, we map the original space domain to a new computational one in which an optimal space mesh can always be constructed. This is accomplished by means of a suitable transformation of coordinates, which involves some parameters that are implicitly defined and whose existence and uniqueness is theoretically established. Such a change of coordinates is thoroughly analyzed from both the theoretical and the computational standpoint, and it is proven to be monotone, infinitely regular, and also simple to implement, so that it can be used without introducing any kind of singularity or incurring in any computational issue.
After applying the transformation of coordinates, a new partial differential problem is obtained which we solve by means of a finite difference scheme with repeated Richardson extrapolation in both space and time. The overall approach exhibits high efficacy. In particular, barrier option prices are computed with an error (in the maximum norm) close to the machine precision (i.e. of the order 10 −11 or 10 −12 ) in less than 0.6 s. As shown by numerical simulations, the improvement over existing methods is due to the combination of the mesh optimization and the repeated Richardson extrapolation.
It is also worth to remark that our mesh optimization approach does not depend on the partial differential equation being solved, and thus it could also be applied to option pricing models other than the Black-Scholes model, such as the standard CEV model (Cox 1996), the CEV model with time dependent parameters (Lo et al. 2000(Lo et al. , 2009 or the fractional Brownian motion (Mandelbrot and Van Ness 1968).
The remainder of the paper is organized as follows: in Sect. 2 the partial differential problem that yields the price of a double barrier option under the Black-Scholes model is briefly presented. In particular, our attention is focused on the case of discretely monitored barriers, but continuously monitored barriers could be considered as well; in Sect. 3 the procedure for constructing an optimal mesh is developed and theoretically analyzed; in Sect. 4 it is shown how to compute the optimal mesh by means of suitable bisection/Newton algorithms whose convergence is a-priory theoretically guaranteed; in Sect. 5 the finite difference scheme is briefly sketched; in Sect. 6 some numerical results are presented and discussed; finally, in Sect. 7 some conclusions are drawn.

The mathematical problem
In this section, the partial differential problem that is satisfied by the price of a double barrier option is briefly recalled. In doing that, we only consider the case of the popular Black-Scholes model, even if the proposed mesh optimization approach can also be used in conjunction with other models (for example the CEV model, see Cox 1996, or the fractional Brownian motion, see Mandelbrot and Van Ness 1968).
Let the price of an asset satisfy the following stochastic differential equation (under the risk-neutral measure): where r and σ are the (constant) interest rate and volatility, respectively, W is a Wiener process.
Let us consider a double barrier call option written on the above asset, with maturity T , strike price E, lower barrier L and upper barrier U . Moreover, let the barriers be applied discretely, at the I b times t b,1 , t b,2 , . . ., t b,I b . Following a common approach, the barrier dates are assumed to be equally spaced in denotes the time between two consecutive barriers. Let us observe that, as will be perfectly clear later on, the mesh optimization technique which we are going to present is not at all restricted to the pricing of the aforementioned kind of double barrier option. By contrast, the case of a put option, the case of continuously monitored barriers, and the case of unequally spaced barrier dates could be dealt with as well. Finally, we will focus our attention onto the case L < E < U (the strike lies between the barriers), which is the one that is usually encountered in the financial practice. Nevertheless, the reader will easily figure it out that the very unfrequent case E < L < U could also be considered.
The barrier option price, denote it V (S, t), can be obtained by recursively solving, for i = I b , I b − 1, . . . , 1, the following partial differential problems: with boundary conditions and final condition where and for i = 1, 2, . . . , I b − 1 we have: That is, first of all we have to solve problem (2)-(3) on the time interval The partial differential problem (2)-(3) must be solved by numerical approximations. In particular, in this paper we will use a finite difference scheme enhanced by repeated Richardson extrapolation in both space and time.

Domain localization of problems (2)-(3)
To solve problem (2)-(3) by finite difference approximation, first of all, domain localization is needed, i.e. the unlimited space domain of problem (2)-(3) must be substituted with a bounded one: where S min and S max have to be chosen such that the error due to the truncation of the domain is negligible. In this paper, since the domain localization is not the focus of our analysis, S min and S max are chosen according to a simple and common procedure, which is described below. First of all, let S b,i denote the price of the underlying asset at the barrier date t b,i , i = 0, . . . , I b . Using the closed-form solution of (1) we have (see, for example, Wilmott 1998): Now, based on the probability distribution of the increments of the Wiener process, we have that for any (given) t Then, we choose S min and S max as follows: The above relations yield suitable values of S min and S max . In fact, if at the barrier the option expires), then at any (given) time t ∈ ]t b,i , t b,i+1 [ the probability that S(t) < S min or S(t) > S max is almost null (smaller than 2.6 × 10 −12 ). Therefore, values of S that lie outside the interval [S min , S max ] are very extreme, and thus they can be reasonably neglected (they contribute with an almost null probability to the option price).
Thus, we can safely replace the boundary conditions (3) with the following ones:

The optimal space mesh
Let us think to solve problems (2), (4), (11) by a finite difference method and to enhance its accuracy by Richardson extrapolation. Then, in order to discretize the derivatives with respect to S in (2), a finite set of nodes (space mesh) in [S min , S max ] has to be employed. However, due to the fact that the solution V (S, t) is not smooth at S = E, S = L and S = U , the space mesh must be chosen appropriately. In particular, to obtain a smooth pattern of error reduction, which is necessary for the Richardson extrapolation to be effective, one node must be located exactly at S = E and other two nodes must be located exactly at S = L and S = U (see Ritchken 1996;Hyman et al. 2000;Pooley et al. 2003). Moreover, to achieve the highest rate of spatial consistency (i.e., the highest rate of reduction of the truncation error due the space discretization), the nodes must be equally spaced (see, for example, Hirsch 1988, Chapter 4;Duffy 2004). In this paper, for the sake of brevity, a space mesh that satisfies the above two requirements is said aligned and uniform.
In addition, we would also like to freely decide how many discretization nodes to use in each one of the intervals it is almost flat and null. Therefore, one should be able to arbitrarily choose both the number of space discretization intervals between S min and L, denote it N L , the number of space discretization intervals between L and U , denote it N M , and the number of space discretization intervals between U and S max , denote it N U .
In substance, we would like to freely specify three integers N L , N M , N U , and, setting N = N L + N M + N U , we would like to have a mesh with N + 1 equally spaced nodes S 1 , S 2 , . . ., S N +1 such that S 1 = S min , S N L +1 = L, S j * = E for some integer j * ∈ {N L + 2, N L + 3, . . . , N L + N M }, S N L +N M +1 = U , S N +1 = S max . Throughout this paper, a space discretization mesh that is aligned and uniform and also such that N L , N M and N U can be decided arbitrarily is referred to as optimal.

The mesh optimization procedure
It is clear that in the space domain S an optimal mesh is impossible to obtain. More precisely, we cannot even construct an aligned and uniform mesh, unless S min , S max , L, U , E, N L , N M and N U do not satisfy very particular (and also restrictive) requirements. For example, let us assume that the nodes S 1 , S 2 , . . ., S N +1 are such that S N L +1 = L and S N L +N M +1 = U . Therefore, if we want this mesh to be uniform, the distance between two consecutive nodes must be: Then, clearly, we can have one of the nodes equal to the strike price E only in the very special circumstance where is an integer value. By contrast, if l E is not an integer, the strike price E is internal to the interval S N L + j E +1 , S N L + j E +2 , where j E denotes the integer part of l E .
Nevertheless, as we are going to show, it is possible to find a function φ : S → R that satisfies the following conditions: then we have: That is, if we denote: the function φ which we are going to determine maps the original space domain S = [S min , S max ] to the domain z , where we can choose N + 1 equally spaced points z 1 , z 2 , . . ., z N +1 such that z 1 = z min , z N +1 = z max and, in addition, z N L +1 , z j E and z N L +N M +1 coincide with the transformed lower barrier, strike price and upper barrier, respectively. Therefore, if we apply the change of coordinates to problem (2), (4), (11), we end up with a partial differential problem that can be discretized by using an optimal mesh (actually, we will have an optimal mesh in the computational domain z ).
Remark 1 Condition (C1) ensures that the function φ is invertible in S , such that we can actually apply the change of variables (16) to solve the partial differential equation (2) (if the function φ was not invertible, then we could not go back from the computational domain z to the true domain S ). Furthermore, the fact that φ is C ∞ is crucial as well, since it guarantees that the partial differential equation that we are going to obtain after applying the transformation of coordinates is a partial differential equation with smooth coefficients (non-smooth coefficients would inevitably spoil the accuracy of the space discretization scheme).
The function φ will be obtained as the composition product of other two functions, which we denote T 1 and T 2 . Precisely, we are going to construct the function φ : S → z as follows: where the functions T 1 and T 2 will be chosen appropriately. In particular, T 1 will be used to transform the interval [L, U ] into a new interval [T 1 (L), T 1 (U )] where we can choose N M + 1 equally spaced points such that one of them coincide with T 1 (E). That is, first of all, by means of a suitable function T 1 , we will accommodate the interval between the two barriers, such to align the space discretization mesh with the strike price. Then, the function T 2 will be used to transform the whole set of points [T 1 (S min ), T 1 (S max )] into the interval z where we can choose nodes z 1 , z 2 , . . ., z N +1 that satisfy condition (C2). In particular, T 2 will accommodate the intervals , such to obtain an optimal mesh in the whole computational domain [z min , z max ].

Constructing T 1
In order to construct the function T 1 the following result is needed. (12), (13), and let j E denote the integer part of l E . Moreover, assume that l E is greater than one and is not an integer. Then, the number of positive solutions of the following transcendental equation in the β variable

Proposition 1 Let l E be defined according to
is exactly equal to one. Moreover, the positive solution of (18), which we denote β * , is such that: Before proving Proposition 1 let us make some useful comments on it.

Remark 2
The assumption that l E is not an integer is not restrictive at all. In fact, if l E is an integer, then Proposition 1 is not needed at all, because the strike price E does already coincide with one of the equally spaced meshpoints S N L +1 , S N L +2 , . . ., S N L +N M +1 and the change of coordinates T 1 is not applied (better saying, in such a case we would choose T 1 to be the identity map, see (34)).

Remark 3
The assumption that l E is greater than one is a sort of "minimal requirement" that one has to satisfy when solving problem (2), (4), (11) by a lattice-based scheme. In fact, if it was l E < 1, then, according to (13), both the lower barrier L and the strike E would be located within the same discretization interval of length S M . This means that the space mesh which we are trying to use is unreasonably too coarse, as it does not even allows us to distinguish between the strike price and the lower barrier. Obviously, if we had l E < 1, then we would simply have to select a larger value of N M , such to have L and E separated by at least one space discretization node.

Proof of Proposition 1 Let us consider the function
We shall prove that f has a unique positive root. To this aim, we note that and, due to the fact that U > E > L, Furthermore, f is of class C ∞ and which yields Since j E is the integer part of l E , and owing to (12) and (13) we have From relations (24) and (25) we easily obtain Relations (21), (22), (26) (and the continuity of f ), imply that f has at least one root β * > 0.
Let us show that β * is the unique positive solution of (18). To this aim, let us assume instead that the function f has two positive roots, say β * 1 and β * 2 , with β * 2 > β * 1 . Then, by using (21), (22), (26), we immediately obtain that there exist two values of β, say β 1 and β 2 , such that Let us consider the following function: Note that, according to (20), (23), the function f 1 takes the same values of f at all those points β such that f (β) = 0. In fact, from f (β) = 0 and (23) it follows that which, if substituted in (20), yields f 1 (β) = f (β). Therefore, according to (27), we must have: These two relations and the continuity of f 1 imply that f 1 must have at least one root in the interval ]β 1 , β 2 ]. Moreover, from (28) we have: which, together with (25), implies By setting (28) equal to zero, we obtain a very simple exponential equation from solving which we can see that f 1 has a unique and positive real root. Now, according to (32) and to the first of (30), such a unique root should be located between 0 and β 1 . Nevertheless, this is in contradiction with the fact that f 1 must have a root in the interval ]β 1 , β 2 ], and so, by absurd, f cannot have more than one positive solution.
Finally, relation (19) easily follows from (21), (26), from the continuity of f and from the fact that We are now in the position to state the first change of coordinates to be applied. The function T 1 : S → R + is chosen as follows: if l e is an integer, where β * is the positive solution of (18), which exists unique according to Proposition (1). Note that T 1 reduces to the identity function if l e is an integer, because in such a case the node S j E does already coincide with the strike price E, and thus, actually, in the price interval [L, U ] no transformation of variables is needed (see also Remark 2). Let us define: Then, let us discretize the interval [y L , y U ] using N M + 1 equally spaced points y L , y L + y, y L + 2 y, . . ., y U , where Then, the following Proposition holds.
and is strictly increasing. Moreover, if y L , y E are defined as in (35) and j E is the integer part of l E given by (13), we have: Proof of Proposition 2 The first part of Proposition 2 follows immediately from (34). Moreover, (37) is a direct consequence of (34), (35), (36) and of the fact that β * is a solution of (18).
Remark 4 According to Proposition 2, one of the N M + 1 equally spaced points that we have chosen in [y L , y U ], i.e. the point y j E , exactly coincides with the (transformed) strike price y E .

Constructing T 2
In the general case, the transformation of coordinates (35) where erf(·) is the error function: Then the equation has a unique positive root, which, in order to stress out its dependence on q, we denote by a(q).

Proof of Lemma 3
Let us observe that h(0) = 1 and lim x→+∞ h(x) = 0. Then, since −1 < q < 0, and h is a continuous function in [0, +∞), equation (40) has at least one positive root. Moreover, we have: so that h is a strictly decreasing function in [0, +∞[, and the uniqueness of a(q) follows.
Let us define: and let the functions g L : [y min , y L [→ R and g U : ]y U , y max ] → R be defined as follows: where g(y, q, y 1 , and a(q) denotes the (unique) positive solution of equation (40).

Lemma 4
The functions g L and g U are infinitely smooth in [y min , y L [ and ]y U , y max ], respectively and are such that lim y→y − L g L (y) = 1, lim U (y) = 0, n = 2, 3, 4, 5, . . . , g L (y) > 0 ∀y ∈]y min , y L [, g U (y) > 0 ∀y ∈]y U , y max [, Proof of Lemma 4 The fact that g L and g U are infinitely smooth in their domains, as well as relations (46) are trivial. Moreover, (47)-(53) can be immediately obtained by taking derivatives in (43)-(45). In particular, we can easily compute the following relations, which will be useful for later purposes: Finally, let us focus our attention on the first of relations (54) (the proof of the second one is analogous). If q L ≥ 0, then the first of relations (54) is rather trivial (it can be obtained by setting y = y min in (43)). Instead, if 0 < q L < 1, the first of (53) follows from setting y = y min in (43) and from the fact that a(q L ) satisfies equation (40).
We are now in the position to state the change of coordinates T 2 : [y min , y max ] → R: and we have the following result.

Constructing
Once the transformations of coordinates T 1 and T 2 are determined according to (34) and (59), respectively, we can finally construct the function φ : S → z using (17), and we have the following proposition.

Proposition 6 The function φ satisfies conditions (C1) and (C2).
Proof Proposition 6 is a direct consequence of Lemma 5 and Proposition 2. Therefore, let us consider the transformation of coordinates: which maps the original domain S to the domain z , where Note that, since the function φ satisfies condition (C1), relation (61) is invertible: Then let us define: By using the transformation of coordinates (61) and (64), the partial differential problem (2), (4), (11) can be rewritten as follows: where Note that the coefficients α 1 and α 2 in (68) need to be computed only one time, at the beginning of the numerical simulation. Moreover, once the inverse function φ −1 is obtained (see the next section), the functions φ and φ in (71) and (72) can be evaluated very easily and quickly on any standard computer.

Applying the transformation of coordinates
To use the change of coordinates (17) we proceed as follows. First of all, once N L , N U , N M are chosen (i.e. decided by the user), we obtain S min and S max according to (10), we calculate j E as the integer part of l E defined as in (13), and thus we compute the unique positive solution β * of equation (18) (this task is accomplished by applying the bisection method, whose convergence is theoretically established below). To this point, the function T 1 is known (see (34)), and so we can evaluate y min , y L , y E , y U , y max (according to (35)), y (according to (36)), q L and q U (according to (42)). Then, we compute a(q L ) and a(q U ) as the unique positive solutions of equation (40) with q = q L and q = q U , respectively (this is done by applying the bisection method, whose convergence is theoretically established below).
In summary, the overall mesh optimization procedure requires us to solve equations (18), (40), (75), (76). Now, these equations do not have exact closed-form solutions. Nevertheless, we can solve them (with machine precision) by using bisection and Newton methods that are very simple to implement and whose convergence can be proven theoretically. In fact, we have the following results: Lemma 7 Let the bisection method be used to compute the (unique) positive solution β * of equation (18), starting from the interval 0, . Then, such an algorithm does always converge (with its usual dichotomic convergence rate).

Lemma 8 Let q be any real number in
Moreover, let the functions g L : [y min , y L ] → R and g U : [y U , y max ] → R be defined as the continuations of the functions g L and g U (see (43) and (44)), respectively, that is Note that, from Lemma 4, it follows that the g L and g U are infinitely smooth in [y min , y L ] and [y U , y max ], respectively.
The solutions of equations (75), (76) are obtained based on the following results.
Lemma 9 Let z be any real number in ]g L (y min ), y L [ and let the Newton method be used to approximate the (unique) solution of equation: starting from the initial guess y = y L if q L < 0 or y = y min if q L ≥ 0. Then, such an algorithm does always converge (with its usual quadratic convergence rate).

Proof of Lemma 9
For the sake of brevity, let us only consider the case q L < 0 (the case q L ≥ 0 is perfectly analogous). According to Lemma 4, g L is an infinitely smooth increasing and convex function. Then, the thesis follows.
Lemma 10 Let z be any real number in ]y U , g U (y max )[ and let the Newton method be used to approximate the (unique) solution of the equation starting from the initial guess y = y U if q U < 0 or y = y max if q U ≥ 0. Then, such an algorithm does always converge (with its usual quadratic convergence rate).

Proof of Lemma 10
The proof is identical to the proof of Lemma 9 and therefore is omitted.
Thus, based on all the above Lemmas, in order to obtain β * , we solve equation (18) by using the bisection method, with starting interval 0, ; in order to obtain a(q L ) and a(q U ), we solve equation (40), with q = q L and q = q U , respectively, by using the bisection method with starting interval [0, − ln(−q)]; in order to obtain φ −1 (z i ), i = 2, 3, . . . , N L , we solve equation (79) by using the Newton method, with initial guess y = y L if q L < 0 or y = y min if q L ≥ 0; in order to obtain φ −1 (z i ), i = N L + N M + 2, N L + N M + 3, . . . , N , we solve equation (80) by using the Newton method, with initial guess y = y U if q U < 0 or y = y max if q U ≥ 0.
All the above root finding algorithms, whose convergence has been theoretically established, can be easily implemented on any computer. Therefore, we can compute almost exact (up to the machine precision) values of β * , a(q L ), a(q U ) and φ −1 (z i ) for i = 2, 3, . . . , N L and i = N L + N M + 2, N L + N M + 3, . . . , N , in an extremely small time (actually, much smaller than the time that is required for solving the partial differential problem (65)-(67)).

The finite difference scheme with repeated Richardson extrapolation in space and time
We can now solve problem (65)-(67) by numerical approximation. To this aim, we shall also discretize the time derivative in (65), which is done by using a mesh of equally spaced M + 1 time levels t 0 , t 1 , . . . , t M , such that t k = k t, k = 0, 1, . . . , M, where t = T M . Then, we apply a finite difference scheme enhanced by repeated Richardson extrapolation in both space and time. Note that, as already mentioned, the Richardson extrapolation will be effective since the space mesh is aligned and uniform, and thus the finite difference scheme achieves a smooth convergence pattern (see, for example, Hairer et al. 1993;Chan et al. 2009). To provide evidence of this fact, in Sect. 6 we will also present a numerical simulation in which the Richardson extrapolation technique is applied without any mesh optimization procedure (so that the mesh is not aligned). The results will be considerably less accurate than that obtained by using the optimized mesh.
The numerical method that we are going to apply is analogous to that proposed in Ballestra (2014) for the pricing of vanilla options (with no barriers) and is based on a finite difference scheme that is second-order accurate in space and first-order accurate in time. Actually, this is a very standard three-point finite difference approximation, and thus, for the sake of brevity, we will only give a brief description of it (for more details the interested reader is referred, for example, to Quarteroni et al. (2007), Tavella and Randall (2000) and Ballestra (2014)). In the partial differential equation (65), the time derivative is evaluated (at z = z i and t = t k ) by the implicit (first-order) Euler scheme: whereas the space derivatives are discretized by using the following (second-order) approximation: where z is given by (14). By substituting (81) and (82) in equation (66), and by imposing the boundary conditions V k 1 = 0 and V k N +1 = 0, we obtain, for each value of k ∈ {0, 1, . . . , M − 1}, a tridiagonal system of linear equations that allows us to obtain the option price at time t k . These systems must be solved recursively for k = M − 1, M − 2, . . ., 0 starting from the knowledge of V M i (since, according to (67), we have V M i = ψ I b (z i ), i = 1, 2, . . . , N + 1). In addition, following (Pooley et al. 2003;Ballestra 2014), at every barrier date t b,i , right after the barrier constraints are imposed, the numerical solution (which is discontinuous at the barriers) is projected onto the space of piecewise linear functions.
The accuracy of the resulting approximation is enhanced by Richardson extrapolation, which is repeated four times in both the space and time variables. Precisely, problem (65)-(67) is solved by using, first of all, a mesh with N + 1 space points and M + 1 time levels, then a mesh with 2N + 1 space points and 4M + 1 time levels, then a mesh with 3N + 1 space points and 9M + 1 time levels, then a mesh with 4N + 1 space points and 16M + 1 time levels, and finally a mesh with 5N + 1 space points and 25M + 1 time levels. That is, the mesh is progressively refined, such that the space length and the time step pass from ( z, t) to ( z 2 , t 4 ), then to ( z 3 , t 9 ), then to ( z 4 , t 16 ), and finally to ( z 5 , t 25 ). Note that the time step is reduced by a factor that is the square of the factor by which the space length is reduced. This is due to the fact that, as already mentioned, the error of the space discretiztion is O( z 2 ), whereas the error of the time discretization is O( t). By doing that, we end up with a sequence of five finite difference approximations, from which an enhanced numerical solution is obtained by standard (repeated) Richardson extrapolation (see, for example, Ballestra 2014; Hairer et al. 1993).

Numerical results
Let us show the computational performances that the mesh optimization approach and the repeated Richardson extrapolation proposed in this paper allow us to reach. To this aim, let V ap (S, t) denote the approximate value of the price of the double barrier option obtained by using the numerical procedure described in Sects. 3 and 4. The error on the option price (at the current time t = 0) is measured using the (discrete) maximum norm: Note that the exact option price V needed in (83) is not available and thus we can only compute it by numerical approximation. Then, to obtain an "exact", or, better saying, reference value of the option price, we employ the approach proposed in Sullivan (2000), Andricopoulos et al. (2003), and Andricopoulos et al. (2007). Now, such a numerical method, being based on a suitable integral formulation of the barrier option price, does not require us to discretize differential operators, but only to compute certain integrals that involve the probability density function of the price of the option's underlying asset. Therefore, if these integrals are approximated by some high-order quadrature rule, the price of the double barrier option can be evaluated very quickly and with extreme accuracy (close to the machine precision). Nevertheless, the approach  (2000), Andricopoulos et al. (2003), and Andricopoulos et al. (2007), albeit mathematically ingenious and very efficient from the computational standpoint, is not very versatile, since it can be used only if the transition probability density function of the price of the stochastic process (1), or at least, the characteristic function associated to it (see, e.g., Fang and Oosterlee 2011), is explicitly known. By contrast, a partial differential approach (such as the one developed in the present manuscript) is much more flexible and can also be applied to those models that do not admit an explicit representation of the transition probability density function (among which, let us recall the CEV model with time-varying parameters, see Lo et al. 2000Lo et al. , 2009. The numerical experiments that follow are performed on a computer AMD Ryzen 3 3200U CPU 2500Ghz 8,00 GB and the software programs are written in Fortran. TEST CASE 1. Let us consider a discrete double barrier option with maturity T = 0.5 years, lower barrier L = 20, strike price E = 25.96, upper barrier U = 30. Note that L, E and U are chosen such that the ratio (13) is not an integer, so that the first case in (34) applies. Finally, we set I b = 100 (number of equally spaced barrier dates), σ = 0.3 and r = 0.03. Based on these data and on relations (10) we compute the lower and the upper bound of the domain S : S min = 17.238837791978 and S max = 34.799910178911.
As is usually done, the mesh parameters N L , N M , N U (as well as the number of time steps M) are varied, even though the error obtained (in the discrete maximum norm) is already of the order 10 −11 when the mesh size parameter N (which is equal to N L + N M + N U ) is not much bigger than one hundredth (see Table 1). The errors and computer times experienced are shown in Table 1, where, for the sake of completeness, we also report the values of β * , q L , q U , a(q L ), a(q U ), y min , y L , y E , y U , y max , i.e. the values of the various quantities involved by the mesh optimization procedure (see Sects. 3 and 4). As we may observe, the mesh optimization procedure allows us to obtain excellent results. In fact, if N L = 5, N M = 40, N U = 10 (so that the overall number of space discretization intervals N is equal to 56) and the number of time discretization steps M is equal to 400, the error Err Max is already of the order 10 −9 .
Moreover, the approach developed in this paper is also extremely efficient from the computational standpoint, since the barrier option price can be computed with an error Err Max of the order 10 −11 in only 0.27 s.

TEST CASE 2.
As a second test case, let us consider a discrete double barrier option with maturity T = 1 year, lower barrier L = 5, strike price E = 7.85, upper barrier U = 10. Moreover, we set I b = 50, σ = 0.2 and r = 0.05. These data and relations (10) yield: S min = 4.104352572684 and S max = 12.196816703217. The errors and computer times obtained are shown in Table 2 (again we also report the values of β * , q L , q U , a(q L ), a(q U ), y min , y L , y E , y U , y max ).
Again, the proposed approach is extraordinarily efficient from the computational standpoint, as an error Err Max of the order 10 −11 can be obtained in only 0.16 s.

TEST CASE 3.
As a third test case, let us consider a discrete double barrier option with maturity T = 1 year, lower barrier L = 45, strike price E = 52.38, upper barrier U = 60. Moreover, we set I b = 400, σ = 0.35 and r = 0.03. These data and relations (10) yield: S min = 39.808655549752 and S max = 67.813849040528. The errors and computer times are shown in Table 3.
As we may observe, the proposed numerical method confirms to be extraordinarily efficient, as an error Err Max of the order 10 −12 can be obtained in only 0.59 s.

Further tests and comparisons
The proposed numerical approach relies on both the mesh optimization procedure and the space-time Richardson extrapolation. Then, it is interesting to investigate To this aim, first we apply the finite difference scheme with space-time Richardson extrapolation (as described in Sect. 5) directly to problem (2)-(4), i.e., without using the mesh optimization procedure. More precisely, we solve the Black-Scholes problem on its original space domain (in the S variable) using equally spaced nodes S 1 = S min , S 2 , S 3 , . . ., S N +1 = S max . Note that none of these points coincides with the strike price of with the barriers (S min and S max are chosen as in the previous section). For the sake of brevity, we only report the results obtained for Test Case 1, since the results experienced in Test Case 2 and Test Case 3 are substantially analogous. As we may observe in Table 4, the finite difference scheme with Richardson extrapolation yields satisfactory levels of computational efficiency, since, for example, the solution is computed with the error 9.25 × 10 −7 in 0.27 s. Nevertheless, as shown in Table 1, if the finite difference scheme with Richardson extrapolation is used in conjunction with the mesh optimization procedure, the results obtained are considerably more accurate (in 0.27 s the barrier option price is computed with the error 3.47 × 10 −11 ). Therefore, the use of the mesh optimization procedure is crucial to obtain the excellent computational performances reported in the previous subsection.
Finally, we employ again the mesh optimization procedure, but we solve problem (65)-(67) without applying the Richardson extrapolation. In particular, we focus on The error due to the time discretization is negligible the space approximation, which we perform by using both the second-order accurate finite difference scheme (82) and the following discretization The above finite difference approximation is employed, for example, by Company et al. (2008), and, by Taylor polynomial expansion, it can be shown to be fourthorder accurate when the solution being computed is regular enough (see Company et al. 2008). Therefore, in the following the scheme (84) will be referred to as fourthorder consistent. Note that for i = 2 and i = N the approximation (84) is not defined and thus at the nodes z 2 and z N we shall use the second-order method (82). In this last numerical experiment, we are not concerned with the error due to the time discretization. Therefore, the time derivative in (65) is computed by applying the (implicit) Euler scheme (81) with a very large number of time discretization steps. In particular, we choose M = 100,000, so that the error due to the time discretization is negligible with respect to the error due to the space discretization (we empirically checked that). Therefore, in summary, we are focusing on the space approximation of problem (65)-(67) alone, which we perform by using both a second-order and a fourth-order consistent finite difference method. We apply the mesh optimization procedure, but not the Richardson extrapolation.
The results obtained are reported in table 5 (again we only consider Test Case 1). As we may observe, the fourth-order consistent scheme fails to reach fourth-order accuracy. In particular, when the mesh size parameters N L , N M and N U are doubled, the error is reduced by a factor that is approximately equal to four, showing that the fourth-order consistent algorithm is only second-order accurate. This is clearly due to the fact that the solution being approximated is not regular at the strike price and at the barriers, and thus the scheme (84), which is based on a fourth-order Taylor polynomial expansion, does not reach its optimal (fourth-order) convergence. In this respect, it is worth observing that, unlike the space discretization scheme (84), the Richardson extrapolation turns out to be very effective to reduce the error of the space approximation even in the presence of a non-smooth solution. To this fact we can give the following explanation: the discretization (84) is applied at every barrier date, when the solution is not a regular function, and thus it fails to achieve fourth-order accuracy. By contrast, the Richardson extrapolation is effective since it is used only at time t = 0, when the solution (obtained by finite difference approximation) is smooth. Finally, the error obtained by using the second-order scheme is (approximately) from 1.3 to 2 times greater than the error achieved by applying the fourth-order consistent scheme.

Conclusions
When pricing double barrier options by lattice-based numerical methods, losses of convergence can occur due to the fact that the solution being computed is not a regular function (see Boyle and Lau 1994;Ndogmo and Ntwiga 2011, as well as the numerical simulations that we performed in Sect. 6.1). In the present paper, we show that extremely high levels of computational efficiency can be reach by employing a mesh optimization procedure coupled with Richardson extrapolation. In particular, first of all an aligned and uniform mesh is constructed by applying a suitable transformation of coordinates, which involves some implicitly defined parameters whose existence and uniqueness is theoretically established. Such a transformation of variables is thoroughly analyzed, both from the theoretical and the computational standpoint, and it is shown to be monotone, infinitely smooth, and also simple to compute. Then, we employ a finite difference scheme enhanced by repeated Richardson extrapolation in both space of time. The overall approach exhibits high efficacy: the price of double barrier options can be computed with an error (in the maximum norm) of the order 10 −11 or 10 −12 in less than 0.6 s. The numerical simulations reveal that the improvement over existing methods is due to the combination of the mesh optimization and the repeated Richardson extrapolation.
Funding Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.