A Front Fixing Implicit Finite Difference Method for the American Put Options Model

In this paper, we present an implicit finite difference method for the numerical solution of the Black-Scholes model of American put options without dividend payments. We combine the proposed numerical method by using a front fixing approach where the option price and the early exercise boundary are computed simultaneously. Consistency and stability properties of the method are studied. We choose to improve the accuracy of the computed solution via a mesh refinement based on Richardson's extrapolation. Comparisons with some proposed methods for the American options problem are carried out to validate the obtained numerical results and to show the efficiency of the proposed numerical methods. Finally, by \textit{a posteriori} error estimator, we find a suitable computational grid requiring that the computed solution verifies a prefixed tolerance.


Introduction
American options are contracts allowing the holder to exercise the option, to sell or to buy an asset, at a certain price at any time prior to and including its maturity date. The pricing of American options plays an important role both in theory and in the real derivative markets. The option pricing model, developed by Black and Scholes [2] and extended by Merton [17], gives rise to a partial differential equation governing the value of an option. Schwartz [24], Brennan and Schwartz [3,4] were the first to apply a finite difference method to price American options. The accuracy of their finite difference method was proved by Jaillet et al. [13]. Other papers present finite difference methods, (see, for example, Hull and White [12], Duffy [6], Wilmott et al. [27]). An alternative approach is based on fronttracking methods that keep track of the boundary and discretize the problem in changing domain; they have been considered in [10], [19], [25].
When we price an American option we also need to determine the optimal exercise moment as a function of the value of the underlying asset. This leads to formulating free boundary of the non-linear problem for the price of the American option looking for a boundary changing in time to maturity, known as the optimal exercise boundary. In particular, the American call option problem is a free boundary problem defined on a finite domain. On the other hand, the American put option problem is a free boundary problem defined on a semi-infinite domain so that it is a non-linear model complicated by a boundary condition at infinity. The difficulty associated with a free boundary problem can be reduced by using a front-fixing method, an approach relying on a change of variable to map the changing domain into a fixed domain.
The front-fixing approach has been considered in several papers. In Holmes and Yang [11], a front-fixing finite element method was used. Tangman et al. [25] introduced a fourth-order accurate finite difference scheme. In [28], the original problem is transformed into one more manageable equation with coefficients not depending on the spatial variable and an explicit update for the location of the free boundary at each time step was used, and in Zhu et al. [30], the secant method was employed to solve the non-linear problems. Moreover, in Tangman et al. [25] and [30], the difference between the prices of American options and European options was computed. Sevčovič [26] studied and approximated the early exer-cise boundary for a class of nonlinear Black-Scholes equations, applying a fixeddomain transformation and using an operator splitting iterative numerical technique. Lauko and Ševčovič [16] introduced a local iterative numerical scheme and compared analytical and numerical solution for the early exercise boundary position computation of the American put option. In [29], Zhang and Zhu presented a predictor-corrector method after the fixed-domain transformation. Nielsen et al. [18] used both explicit and implicit schemes for solving Black-Scholes model of American options. While Company et al. [5] presented an explicit finite difference method for the free boundary value problem under logarithmic front-fixing transformation. In [7,8], Fazio proposed one a posteriori error estimator for the numerical solution for the American put option obtained by an explicit finite difference scheme.
In this paper, we present an implicit finite difference scheme combined with the use of the front-fixing method in order to solve the American put option problem without dividend payments. Preliminary results have been presented to the International Conference ICNAAM2015 [9]. We use a front-fixing variable transformation to reformulate the variable domain problem into a non-linear problem on a fixed rectangular domain. For the non-linear problem, we introduce a suitable choice of truncated boundary that allows to impose the asymptotic boundary condition. Then, the original problem is transformed into a new non-linear partial differential equation where the free boundary appears as a new variable and computed as part of the solution. We develop an implicit finite difference method, we investigate the consistency and the stability. Finally, we choose to validate the obtained numerical results via a mesh refinement and the Richardson's extrapolation and we report the comparison with numerical methods available in the literature.

The American put options model
Let us suppose that at time t the price of a given underlying asset is S. We consider here the following mathematical model for the value P = P(S, τ) of an American put option to sell the asset where τ = T − t denotes the time to maturity T , S * (τ) is a free boundary, that is the unknown early exercise boundary, σ , r and E are given constant parameters representing the volatility of the underlying asset, the interest rate and the exercise price of the option, respectively. Equation (1) is known as the Black-Scholes-Merton equation and it has been developed by the three economists Fischer Black, Myron Scholes and Robert Merton in 1973, [2] and [17].

The front-fixing method
The front-fixing method is widely employed for solving free boundary problems. The basic idea of the front fixing method is to use a variable change in order to remove the free boundary and, then, to transform the original equation into a new non-linear partial differential equation on a bounded domain, where the free boundary appears as a new unknown of the problem. The main advantage of the front-fixing method is that the free boundary is computed directly. Then, we first transform the Black-Scholes equation into a new parabolic equation over a bounded domain, we introduce a truncated boundary and, then, we use finite difference schemes for the new approximate problem over a bounded domain.
According to Wu and Kwok [28], we consider the dimensionless new variables it follows that S f (τ) is mapped on the fixed line x = 0, 0 ≤ p(x, τ) ≤ 1 and 0 ≤ S f (τ) ≤ 1. Then, the put option problem (1) can be rewritten as follows on the domain defined by 0 < τ ≤ T and 0 < x < ∞.
In order to solve numerically the obtained problem (3-6), we have to consider a finite computational domain. Then, we introduce a truncated boundary value x ∞ , which is a suitable large value where it would be convenient to impose the asymptotic boundary condition. In other words, we replace the asymptotic boundary condition (5) with the side condition For the choice of x ∞ and the accuracy of the related numerical solution, we can refer to the study by Kangro and Nicolaides [14]. By setting an integer J and a positive value µ, we define the step-sizes where ⌈·⌉ : IR + → IN is the ceiling function which maps a real number to the least integer that is greater than or equal to that number. Therefore, µ is the grid-ratio Therefore, within the finite domain, we can introduce a mesh of grid-points for j = 0, 1, . . ., J and n = 0, 1, . . ., N. We would like to define a numerical scheme that allows us to compute the grid values for j = 0, 1, . . ., J, n = 0, 1, . . ., N and the free boundary values
The initial conditions (4) are Moreover, by considering the governing differential equation (3) at x 0 = 0, τ > 0 and taking into account the side conditions (6), one gets a new boundary condition σ 2 2 (see Wu and Kwok [28], Zhang and Zhu [29] or Kwok [15, p. 341]), and its central finite difference discretization From the two boundary conditions (6), using a central finite difference formula at time n + 1 and considering (10), we obtain where x −1 = −∆x is a fictitious point out of the computational domain. Considering µ = ∆τ/(∆x) 2 and rearranging the (8), our implicit numerical scheme can be written as followsā for j = 1, 2, . . ., J − 1 and n = 0, 1, . . ., N − 1, wherē Considering (11) and putting j = 1 in (12), we get Putting For j = 3, 4, · · · , J − 1 we have the equations Then, at each time step we obtain a system J equations, given by (14)- (17), in J unknowns, p n+1 2 , p n+1 3 , · · · , p n+1 J and S n+1 f . The system (14)- (17) can be written in the following compact form The system (18) can now be written as a non-linear problem in the form The implicit discretization leads to a system of non-linear equations (19) for the price and the location of the free boundary at each time step.

Consistency and stability
In this section, we discuss the consistency and stability of the implicit finite difference scheme (8).

Consistency
We write the PDE (3) and the numerical scheme (8) as follows In order to study the consistency, let us take an arbitrary point (x, τ) in the domain where we denote withp n+1 j = p(x j , τ n+1 ) the exact solution value of the PDE and withS n+1 f = S f (τ n+1 ) the exact solution of the free boundary. Using Taylor's expansion, assuming the existence of the continuous partial derivatives up to order two in time and up to order four in space, we obtain Then, the local truncation error is In our case, we have also to consider the boundary conditions (6) and (9) The numerical scheme for the boundary conditions is Using Taylor's expansion, the local truncation error for the boundary conditions is of O(∆x) 2 So, assuming the existence of the continuous partial derivatives up to order two in time and up to order four in space of the solution of the problem (3-6), the implicit finite difference method defined by (8)

Stability
Now, we perform the Von Neumann analysis to investigate the stability of the implicit method. The Von Neumann analysis is only valid for linear problems with constant coefficients, it is not possible to apply it to non-linear problems or to problems with variable coefficients. In order to apply the method of Von Neumann, it is necessary to linearize the model and freeze the coefficients, considering the problem locally [1]. Then, the von Neumann analysis can be applied, and a stability condition can be derived, this condition will depend on the frozen coefficients involved.
In our context, the non-linear nature of the model (3) is due to the presence of the term with S f unknown of the problem. Then, in the stability analysis, we decide to take it into account and to derive the stability condition in relation to the value of S f and of its first derivative. Moreover, note that we approximate the term (22) with that assumes nonpositive values because S n f is positive and not increasing. Using the Fourier analysis, we set p n+1 j = λ p n j p n+1 j±1 = λ e ±ik∆x p n j where λ = λ (k∆x) is called the amplification factor. Substituting these expressions in the numerical method (12) and dividing by p n j , we obtain λ = 1 a n+1 e −ik∆x +b +c n+1 e ik∆x .
By applying the Euler's formulas and by (13), after some manipulations, we derive . Now if we compute the magnitude of amplification factor λ , then we have Because A > 0, (1 + A) 2 > 1 and B 2 > 0, thus |λ | < 1. So, we can conclude that the implicit finite difference method defined by (8) for the fixed domain problem (3-6) is unconditionally stable.
In figure (1), we show the amplification factor module for implicit method at different time steps τ n . In fact, in the stability analysis, we have decided to evaluate the amplification factor λ depending on the value of (23). We observe that the variation of (23) does not influence the stability of the implicit method. In figure (2) we show the amplification factor module for different values of µ at final time τ = T . We observe that the unconditionally stable implicit finite difference scheme produces results that are stable for any value of grid-ratio µ.
We can compare the result reported in figure (2) with the ones reported in figure (3), where we show the amplification factor module λ , for different values of µ, of the explicit method [5] In order to evaluate λ , we use the same assumptions given for the implicit method. We observe that the explicit finite difference scheme produces results that are stable only for values of grid-ratio µ < 26.
The main limitation of the explicit method is the restriction on the choice of ∆τ and ∆x that have to be chosen in relation at the values of the parameters of the model σ and r, the volatility of the underlying asset and the interest rate, respectively. On the other hand, the new implicit method is unconditionally stable, no restriction on step sizes is required.

A posteriori error estimator
In this section, we describe a posteriori error estimator for the American put options problem based on Richardson's extrapolation. In this way, we can find a more suitable computational grid requiring that the numerical solution obtained by the implicit method verifies a prefixed error tolerance.
For a scalar U of interest, either a value of the solution p n j or a free boundary value S n f , the numerical error e can be defined by where u is the exact, usually unknown, value. When the numerical error is caused prevalently by the discretization error and in the case of smooth enough solutions the global error can be split into a sum of powers of the inverse of N where C 0 , C 1 , C 2 , . . . are coefficients that depend on u and its derivatives, but are independent of N, and q 0 , q 1 , q 2 , . . . are the true orders of the discretization error, see Schneider and Marchi [23] and the references quoted therein. Each q k is  usually a positive integer with q 0 < q 1 < q 2 < · · · and all together constitute an arithmetic progression of ratio q 1 − q 0 . The value of q 0 is called the asymptotic order or the order of accuracy of the method or of the numerical value U N . By replacing into equation (25) N = N g and N = N g+1 and subtracting, to the second obtained equation the first times (1/s) q 0 , s = N g+1 /N g , we get the first extrapolation formula that has a leading order of accuracy equal to q 1 . This type of extrapolation is due to Richardson [20,21]. Taking into account equation (26), we can conclude that the error estimator by a first Richardson's extrapolation is given by where q 0 is the order of the numerical method used to compute the numerical solutions. Hence, (27) gives the real value of the numerical solution error without knowledge of the exact solution. In comparison with (27), a safer error estimator can be defined by Of course, q 0 can be estimated with the formula where u is again the exact solution (or, if the exact solution is unknown, a reference solution computed with a suitable large value of N), and both u and U g+1 are evaluated at the same grid-points of U g . Within the above framework, in order to improve the numerical accuracy by using only a small number of grid-nodes, we can generalize (26) introducing the following repeated extrapolation formula where g ∈ {0, 1, 2, . . ., G − 1}, k ∈ {0, 1, 2, . . ., G − 1}, s = N g+1 /N g is the grid refinement ratio, and q k is the true order of the discretization error. The formula (30) is asymptotically exact in the limit as N 0 goes to infinity if we use uniform grids. We notice that to obtain each value of U g+1,k+1 we need two computed solutions U in two adjacent grids, namely g + 1 and g at the extrapolation level k. For any g, the level k = 0 represents the numerical solution of U without any extrapolation. We recall that the theoretical orders of accuracy of the numerical values U g,k , with N = N g and k extrapolations, verify the relation valid for k ∈ {0, 1, 2, . . ., G − 1}.

Numerical results
In this section, we show numerical results obtained by using the implicit finite difference scheme. Comparisons with some numerical methods available in the literature for the American put options problem are carried out. The parameters considered for solving the American put option problem (3-6) are the following In order to choose the truncated boundary value x ∞ , we computed the numerical solution S N f for three different values of x ∞ . As we can see in table 1, the change of these values does not affect much the numerical solution; for this reason, we decided to set x ∞ = 1.  In figure 4, we show the plots of p N j versus x j and of S n f versus τ n ; these results are obtained by the implicit finite difference scheme, considering the parameters (32) with N = 320 and µ = 20. In figure (5), the 3D plot of numerical solution p n j is shown.
In order to solve the non-linear system, we use the MATLAB routine "fsolve". Initially, we used the classical Newton iteration method, but the obtained numerical results showed less accurate when compared with ones obtained by Matlab routine fsolve. By our numerical experiments, we conclude that the MATLAB routine, principally regarding the choice of the initial iterate, reveals more robust than the Newton iteration method and then more suited for our application.
Looking at the results listed in Table 1, we realize also that, fixing a value of the truncated boundary x ∞ , the computed values of S N f for different values of the grid-steps are in agreement only up to the first decimal place. For this reason, we decide to improve the numerical accuracy by performing a mesh refinement applying repeated Richardson's extrapolation. The implicit difference scheme is accurate of first-order in time and second-order in space, then, we can perform a mesh refinement keeping constant the grid-ratio µ, so that we end up with secondorder truncation error T n j = O(∆τ 2 ) in time. Therefore, the global error is of firstorder, which is the value of q 0 = 1. In our case the sequence of s q k , for k = 0, 1, . . ., with s = 4 and q k = k + 1 is given by 4, 16, 64, 256, 1024  by Fazio [7], with S N f ≈ 0.86222 computed by Nielsen et al. [18] and with S N f ≈ 0.8628 found by the explicit method of Company et al. [5]. In Table 3, we compare the option price P(S, T ) obtained by different methods with the following parameters We report the "true value" as the reference offered in [22], the penalty method (PM) of Nielsen et al. given in [22], the explicit method (EM) of Company et al. [5] with ∆x = 0.02, the explicit method with the Richardson's extrapolation (EMR) of Fazio [7] with µ = 5 and our implicit method (IM) without extrapolation, setting ∆x = 0.02 and µ = 5, and with Richardson's extrapolation (IMR) with µ = 5. In order to find the value of the option price P(S, T ) in correspondence to each different asset price S we use a piecewise cubic spline interpolation. Our goal is to solve the American put option problem with a given tolerance ε, where 0 < ε ≪ 1; then we use the error estimator defined by equation (27), or alternatively by equation (28). To this end, we should solve the given problem twice, for two grids defined with given values of J g = J and J g+1 = 2J of space intervals but for the same value of the grid-ratio µ. The corresponding time intervals N g and N g+1 verify the relation s = N g+1 /N g . Hence, we can apply (componentwise) to p n and to S n f the error estimator formula (27), or (28). Then, we can verify whether, for n = 1, 2, . . ., N,   If the two inequalities (34) hold true, for n = 1, 2, . . ., N, then we can accept the numerical solution computed on the grid defined by J g+1 and N g+1 , otherwise, we have to increase these two integers and repeat the computation. Fig. 6 shows the error estimator results computed by setting ε = 0.005. We set µ = 20 and start with J 0 = 5 and J 1 = 10 repeating the computation by doubling the number of spatial grid-intervals if the required accuracy is not achieved. Our algorithm stops when J 5 = 160 that for µ = 20 means N 5 = 1280. From the numerical results shown in Figure 6 we can easily realize that the greatest errors are found within a few time steps for both numerical methods. The error of the implicit method is almost twice of the error of the explicit scheme only in the first time steps, but it becomes smaller in the next time steps, so we can conclude that the numerical results obtained by the implicit method are more accurate than those of the explicit one. We suggest that in order to obtain better accuracy also in the first time steps it can be developed an adaptive version of the both finite difference schemes.

Concluding Remarks
In this paper, we consider an American put option model, a free boundary problem defined on an infinite domain. We overcome the numerical difficulty of solving a free boundary problem using a front fixing approach combined with finite difference schemes. We propose an unconditionally stable implicit finite difference scheme and we prove the consistency. In order to improve the accuracy solution, we use Richardson's extrapolation. Comparisons with some methods available in the literature are carried out to validate the obtained numerical results. Finally, by a posteriori error estimator, we find a suitable computational grid requiring that the computed solution verifies a prefixed tolerance.