Solving American Option Pricing Models by the Front Fixing Method: Numerical Analysis and Computing

and Applied Analysis 3 The approximate value of p(x, τ) at the point x j and time τ n is denoted by p j ≈ p(x j , τ n ). Then a forward-time centralspace explicit scheme is constructed for internal spacial nodes as follows:


Introduction
American options are contracts allowing the holder the right to sell (buy) an asset at a certain price at any time until a prespecified future date.The pricing of American options plays an important role both in theory and in real derivative markets.
The American option pricing problem can be posed either as a linear complementarity problem (LCP) or a free boundary value problem.These two different formulations have led to different methods for solving American options.
The most algebraic approach of LCPs for American option pricing can be found in [1,2] and the references therein.
In this paper we follow the free boundary value approach.McKean [3] and van Moerbeke [4] show that the valuation of American options constitutes a free boundary problem looking for a boundary changing in time to maturity, known as the optimal exercise boundary.The problem of finding optimal exercise boundary can be treated analytically or numerically.With respect to the analytical approach, Geske and Johnson [5] obtained a valuation formula for American puts expressed in terms of a series of compound-option functions.See also Barone-Adesi and Whaley [6], MacMillan [7], and Ju [8].
Free boundary value problems are challenging because one has to find the solution of a partial differential equation that satisfies auxiliary initial conditions and boundary conditions on a fixed boundary as well as on an unknown free boundary.This complexity is reduced by transforming the problem into a new nonlinear partial differential equation where the free boundary appears as a new variable of the PDE problem.
This technique which originated in physics problems is the so called front fixing method based on Landau transform [18] to fix the optimal exercise boundary on a vertical axis.The front fixing method has been applied successfully to a wide range of problems arising in physics (see [19]) and references therein.In 1997, Wu and Kwok [20] proposed the front fixing technique for solving free boundary problem to the field of option pricing.Other relevant papers related to the fixed domain transformation are Nielsen et al. [21], Šev č ovi č [22], and Zhang and Zhu [23].All these papers use finite difference schemes after the corresponding fixed domain transformations.Nielsen et al. [21] use both explicit and implicit schemes while Šev č ovi č [22] uses difference schemes of operator splitting type.

2
Abstract and Applied Analysis Wu and Kwok [20] transform the original problem into one more manageable equation because the coefficients of the transformed equation do not depend on the spatial variable.In [20], one uses a tree-level scheme and performs numerical tests versus the binomial method.The drawback of the initialization and iteration of the three-level method is overcome by Zhang and Zhu [23] using a predictor-corrector approach.
The front fixing method combined with the use of an explicit finite difference scheme avoids the drawbacks of alternative algebraic approaches since it avoids the use of iterative methods and underlying difficulties such as how to initiate the algorithm, when to stop it, and which is the error after the stopping.
As the best model can be wasted with a disregarded numerical analysis, in this paper, we focus on the numerical issues of explicit finite difference scheme for American put option problem after using the fixed domain transformation used in [20].Consistency of the numerical scheme with the PDE problem is stated.Conditions on the size of discretization steps in order to guarantee stability and positivity of the numerical solution are given.Furthermore, the numerical solution preserves properties of the theoretical solution presented in [24], such as the positivity and the monotonicity of the optimal exercise boundary.
This paper is organized as follows.In Section 2, dimensionless fixed domain transformation of the American put option problem is introduced as well as the discretization of the continuous problem.Section 3 is addressed to guarantee the positivity of the numerical solution as well as the nonincreasing time behaviour of the optimal exercise boundary.In Section 4, stability and consistency are treated.Finally, in Section 5, numerical experiments and comparison with other methods are illustrated.

Fixed Domain Transformation and Discretization
It is well-known that the American put option price model is given by [1] the moving free boundary PDE as follows: together with the initial and boundary conditions: where  =  −  denotes the time to maturity ,  is the asset's price, (, ) is the option price, () is the unknown early exercise boundary,  is a volatility of the asset,  is the riskfree interest rate, and  is the strike price.Note that if asset price  ≤ (), the optimal strategy is to exercise the option, while if  > (), the optimal strategy is to hold the option.
Let us consider the dimensionless change of the following variables [20]: that transforms the problem (1)-( 2) into the normalized form as follows: where    denotes the derivative of   with respect to .The new boundary and initial conditions are  (, 0) = 0,  ≥ 0, lim The Equation ( 4) is a nonlinear differential equation on the domain (0, ∞) × (0, ].In order to solve numerically problem (4)-( 9), one has to consider a bounded numerical domain.Let us introduce  max large enough to translate the boundary condition (8); that is, ( max , ) = 0. Then the problem (4)-( 9) can be studied on the fixed domain [0,  max ] × (0, ].The value  max is chosen following the criterion pointed out in [25].
The approximate value of (, ) at the point   and time   is denoted by    ≈ (  ,   ).Then a forward-time centralspace explicit scheme is constructed for internal spacial nodes as follows: By denoting  = /ℎ 2 , the scheme ( 11) can be rewritten in the following form: where From the boundary conditions ( 6) and ( 7), we can obtain where  −1 = −ℎ is an auxiliary point out of the domain.By considering (4) at the point  0 = 0,  > 0 (see [26], p. 341), and replacing of the boundary conditions ( 6) and ( 7) into ( 4) at (0 + , ), one gets the new boundary condition as follows: (see [20,23,26]), and its central difference discretization is From ( 14) and ( 16), the value of   −1 can be eliminated obtaining the relationship between the free boundary approximation    and   1 , where By using the scheme (12) for  = 1 and evaluating (17) at the (+1) st time level, the free boundary  +1  can be expressed as where After expression (19), the value  +1  can be replaced in ( 12), (17), and ( 14) to obtain values  +1  , 0 ≤  ≤ .Then the numerical scheme for the problem ( 4)-( 9) can be rewritten for any  = 0, . . .,  − 1 in the following algorithmic form: where with the initial conditions

Positivity and Monotonicity
In this section, we will show the free boundary nonincreasing monotonicity as well as the positivity and nonincreasing spacial monotonicity of the numerical option price.Let us start this section by showing that constant coefficients , , and  of scheme (12) for the transformed problem (4)-( 9) are positive under appropriate conditions on the stepsize discretization ℎ and .
Lemma 1. Assuming that the stepsizes ℎ and  satisfy the following conditions: then the coefficients of the scheme , , and  are nonnegative.If  =  2 /2, then under the condition C2, coefficients , , and  are nonnegative.
The following lemma prepares the study of positivity of the numerical solution    as well as the monotonicity of the free boundary sequence    , that will be established in a further result.Lemma 2. Let {   ,    } be the numerical solution of scheme (21)-( 26) for a transformed American put option problem (4) and let   be defined by (20).Then under hypothesis of the Lemma 1, for small enough ℎ, one verifies the following.
Let us assume the induction hypothesis that conclusions hold true for index  − 1; that is, Now we are going to prove that conclusions hold true for index .Let us consider value of   for  ≥ 1.By denoting then from (20), For  > 2, using Taylor's expansion, the approximation of the involved derivatives, and boundary conditions ( 6), (7), and (15), From ( 37) and ( 34), numerator   takes the form and verifies   > 0 since  < ℎ/(ℎ +  2 ) under condition C2 of lemma and for ℎ < 1.
In order to prove the positivity of { +1  }, it will be useful to show the nonnegativity of coefficients ã and c appearing in (21).Coefficient ã is positive since In order to show the positivity of the coefficient c , note that from ( 13) and (39), the sign of c is the same as the sign of (2ℎ  +   −   ) and from ( 13), (38), (39), and    ≤ 1, one gets the following for small enough values of ℎ: Under hypotheses of induction (33) together with positivity of coefficients ã and c , the positivity of { +1  } is proved.Moreover, { +1  } is nonincreasing with respect to index  from (24), since Summarizing, the following result has been established.

Stability and Consistency
In this section, we study the stability and consistency properties of the scheme ( 21)- (26).For the sake of clarity in the presentation knowing that several different concepts of the stability are used in the literature, we begin the section with the following definition.
Proof of Theorem 5. Since for each fixed , {   } is a nonincreasing sequence with respect to , then according to the boundary condition (22) and positivity of    since (28), one gets Thus, the scheme is ‖ ⋅ ‖ ∞ -stable.
Classical consistency of a numerical scheme with respect to a partial differential equation means that the exact theoretical solution of the PDE approximates well the solution of the difference scheme as the discretization stepsizes tend to zero [29].However, in our case, we have to harmonize the behaviour not only of the PDE (4) with the scheme (24), but also the boundary conditions ( 6), (7), and (15) with the numerical scheme for the free boundaries ( 21) and (20).With respect to the first part of consistency, let us write the numerical scheme (24) in the following form: Let us denote by p  = (  ,   ) the exact theoretical solution value of the PDE at the mesh point (  ,   ), and let S  =   (  ) be the exact solution of the free boundary at time   .The scheme ( 46) is said to be consistent with Assuming the existence of the continuous partial derivatives up to order two in time and up to order four in space, using Taylor's expansion about (  ,   ), one gets where Equations ( 50) and (51) show the local truncation error of the numerical scheme (24) with respect to the PDE (4).Note that from (50) and (51), In order to complete the consistency of the solution of the free boundary problem with the scheme ( 21)-( 26), it is necessary to rewrite the boundary conditions ( 6), (7), and (15) in the following form: Furthermore, we have finite difference approximation for the following boundary conditions: It is not difficult to check using Taylor's expansion that the local truncation error satisfies The truncation error for the boundary condition behaves as ℎ 2 .Theorem 6.Assuming that the solution of the PDE problem (4)-( 9) admits two times continuous partial derivative with respect to time and up to order four with respect to space, the numerical solution computed by the scheme (21)-( 26) is consistent with (4) and boundary conditions of order two in space and order one in time.

Numerical Experiments
In this section we illustrate the previous theoretical results with numerical experiments showing the potential advantages of the proposed method such as preserving the qualitative properties of the theoretical solution [24] such as positivity and monotonicity.A comparison with other approaches is also presented in this section.

Confirmations of Theoretical Results
. In this experiment we check that the stability condition can be broken showing that in such case the results could become unreliable.Let us consider an American put option pricing problem as in [21] with the following parameters: We consider fixed space step ℎ = 0.01 that satisfies condition C1 of Lemma 1, and change time step (Figure 1).Numerical tests show that the condition C2 is critical for positivity of coefficient  and, as a result, for the stability of the scheme.The monotonicity of the free boundary is numerically shown by numerical tests (Figure 1(a)).In Figure 1(b), one can see that when condition C2 is not satisfied, spurious oscillations appear so the monotonicity is broken.It was theoretically proved in Section 4 that the scheme has order of approximation (ℎ 2 ) + ().The results of the numerical experiments for the problem with the following parameters: are presented in Table 1.For "True" values, we use the reference values offered in [30].We consider the different space steps for fixed time step  = 0.002 to guarantee the stability.The root-mean-square error (RMSE) is used to measure the accuracy of the scheme.The last row presents the CPU-time in seconds for each experiment.The plot of the RMSE against computational time is presented in Figure 2.
To check the order of approximation in time, the space step ℎ is fixed (ℎ = 10 −2 ) and time step  is variable.From Table 2 by analogous to (58), formula one gets (5 ⋅ 10 −4 , 10 −3 ) = 0.80, that is close to 1.It proves the second order of approximation in space and the first order in time.
(59) By the numerical tests, it is shown that the considered scheme is stable for  ≤ 24 while the condition for a stability of the scheme in [21] is  ≤ 6.The results are represented in Table 3.The front fixing method with the logarithmic transformation has a good advantage: weaker condition for the stability.It reduces the computational costs.
Let us compare front fixing method with another approach [31], based on Mellin's transform.The parameters of the problem are  = 0.0488,  = 0.3, and  = 0.5833.
To compare results of the explicit front fixing method with Mellin's transform [31], we have to multiply our dimensionless value on  = 45.Then, for Mellin's transform method   () = 32.77,while for front fixing method   () = 32.7655.
Close to maturity, exercise boundary for the American put can be analytically approximated [32] as follows: This approximation can be used only near expiry date.In Figure 3, the value of the free boundary for both approaches is presented.Near initial point, front fixing method and analytical approximation (60) give close results.
We compare explicit front fixing method (FF) with ℎ 1 = 0.001 and ℎ 2 = 0.002 and  = 5 for American put with other numerical methods shown in [15] in Table 4 for the following problem parameters: There are several considered methods as follows: (i) penalty method (PM) is considered in [21,33]; (ii) Ikonen and Toivanen [34] proposed an operator splitting technique (OS) for solving the linear complementarity problem;  (iii) Han-Wu algorithm (HW) transforms the Black-Scholes equation into a heat equation on an infinite domain [35]; (iv) the front fixing method proposed by Wu and Kwok (WK) [20].
Apart from the previous comparisons showing that the method is competitive, we remark that the proposed scheme here guarantees the positivity of the solution as well as the monotonicity of both solution and the free boundary.

Conclusion
The front fixing method for American put option is considered.We found that the proposed explicit numerical scheme is stable under appropriate conditions (see Lemma 1) and consistent with the differential equation with the second order in space and first order in time.Moreover, the monotonicity and positivity of the solution and the free boundary under that condition are proved in agreement with Kim's results in [24] about the properties of the theoretical solution.The theoretical study is confirmed by numerical tests.It is shown that the condition (27) is critical; that is, violation of the condition leads to destabilization of the scheme.
By the comparison with other approaches, advantages of the method are discovered.The logarithmic transformation requires weaker stability condition than Landau's transformation, presented in [21].The explicit calculations with quite weak stability conditions avoid iterations to solve the nonlinear partial differential equation.It makes the method simpler and reduces computational costs.

Table 1 :
Convergence in space for fixed time step  = 2 ⋅ 10 −3 for the problem with parameters (56).

Table 2 :
Convergence in time for fixed space step ℎ = 0.01 for the problem with parameters (56).

Table 4 :
Comparison with other methods: put option value for different asset prices.