Adaptive algorithm for solution of early exercise boundary problem of American put option implemented in Mathematica

The paper is focused on American option pricing problem. Assuming non-dividend paying American put option leads to two disjunctive regions, a continuation one and a stopping one, which are separated by an early exercise boundary. We present variational formulation of American option problem with special attention to early exercise action effect. Next, we discuss financially motivated additive decomposition of American option price into a European option price and another part due to the extra premium required by early exercising the option contract. As the optimal exercise boundary is a free boundary, its determination is coupled with the determination of the option price. Therefore, a closed-form expression of the free boundary is not attainable in general. We discuss in detail a derivation of an asymptotic expression of the early exercise boundary. Finally, we present some numerical results of determination of free boundary based upon this approach. All computations are performed by sw Mathematica, and suitable numerical procedure is discussed in detail, as well.


Introduction
American option pricing problems given with different levels of abstraction are presented in [3], [4], [6], and [7].We know that in contrast to European options, American options can be exercised at any time between the writing and the expiration of the option contract.Exercising flexibility gives American option more profiteering opportunity than European option offers, and it means that payoff functions of American options take the following general form, where a subscript C, and P denotes an option type, i.e. call, and put, respectively V c (S,t) max(S -K, 0) V p (S,t) max(K -S, 0) here S is an underlying risky asset price, e.g. stock price, 0 t T, T is the option expiration date, and K is the option exercise price, called strike price, too.Now, let us focus on non-dividend paying American put option.In fact, when S falls below a certain point, one should exercise the option immediately to avoid loss.For example, if at time t, the current stock price S t is S t < K(1 -e -r(T-t) ), then the option holder should exercise the option immediately.Here, r denotes a riskfree interest rate.In fact, the payoff at the option expiration data T will never exceed K in any case.Though, if the option is exercised at the time t, the immediate gain is K -S t > K -K(1 -e -r(T-t) ) = Ke -r(T-t) , and by depositing this gain in a saving account, assuming a regular compounding process applies, the total payoff will exceed K at t = T, clearly.
The optimal exercise boundary is reasonably defined by mapping ī:[0,T]ĺī(t), thus specifying the underlying asset price at ī, denoted S Ș =ī(t e ), precisely.Here, we adopt the subscript Ș instead of usual t in order to point out the specific asset price at the optimal exercise time t e , in particular.
The optimal exercise boundary is not given a priori, and has to be determined together with option pricing function V(S,t) by solving the corresponding option pricing problem.Here, we already drop the subscript P for simplicity, as we are handling the American put option in the paper excusively.Using optimal exercise boundary ī(t), the subregions Ȉ 1 and Ȉ 2 are defined in following way where Ȉ 1 is called the continuation subregion, since the payoff is zero when S > S Ș , and the holder should continue to keep the option, whereas Ȉ 2 is called the stoping subregion.strike price.Noting, the continuation subregion Ȉ 1 is sometimes called retained subregion alternatively, whereas the stopping subregion Ȉ 2 is called selling subregion, or exercise subregion, as well.In Figure 1, there are depicted both subregions Ȉ 1 , Ȉ 2 , the expiration date T, and the option exercise price, all schematically.The horizontal axis represents an underlying asset price S, and the vertical axis carries the time t elapsed from entering the option contract into power.Following classical procedure based on construction of self-financing portfolio, ¨-hedging principle, and It‫'۾‬s formula, which is common and typical for derivation of almost all models within the field of financial option pricing problems, we can infer that function V(S,t) satisfies the Black-Scholes PDE in the subregion Ȉ 1 provided V(S,t) belongs to C 2,1 (Ȉ 1 ) class of functions, in classical formulation framework, see for example [3] and [6].
The boundary conditions on Ȉ 1 are specified in following way.First, on the optimal exercise boundary ī, it holds Further, the terminal condition at t = T, and the condition when S ĺ +, are following Concluding, we see that American put option pricing problem leads to finding a function pair {V(S,t), ī(t)} in subregion Ȉ 1 satisfying the PDE (5) and boundaryterminal conditions (6) and (7).

Decomposition formula for American put option
Early exercise premium is discussed from several viewpoints in [1], and [2].Following [4], Chapter 6.3, we present an interesting decomposition formula for American put option price V(S,t) which gives a link to European put option price V E (S,t) in following way here V E (S,t) is the well-known European put option price on the same underlying asset, which is computed by closed Black-Scholes formula, see for example [2], [3], and [6], as well, and e(S,t) is the early exercise premium given as follows The function G(S,t;ȟ,Ș) is fundamental solution, sometimes called Green function, too, of the well-known Black-Scholes PDE, being expressed explicitly in relation (5).Inspecting (13) thoroughly, we reckon the upper bound S Ș of the inner integral to cause direct dependence of early exercise premium e(S,t) upon the optimal exercise boundary ī, evidently.Now, we refer to [4], Chapter 6.3, and Theorems 6.3 and 6.4 therein, for more technical details regarding the function G(S,t;ȟ,Ș).However, we ought to present here at least the function G(S,t;ȟ,T) in explicit form which represents a solution of the Black-Scholes PDE with special terminal condition expressing concentrated effect at an arbitrary point ȟ where į(.) is the well-known Dirac function being given by its properties į(x) = + LI x = 0, otherwise į(x) = 0, and ‫‬ ‫ݔ݀)ݔ(ߜ‬ = 1 +λ െλ .

Early exercise boundary equation
As the optimal early exercise boundary ī is not known a priori but has to be determined together with pricing function V(S,t), therefore, in general, the problem is called free-boundary value problem for parabolic PDE.Albeit we have already formulated the problem by expressions ( 5)-( 7) explicitly, the main difficulty therein is that one needs to solve both for unknown function V(S,t) and unknown boundary ī(t), simultaneously.Since any American option has early exercise term in the contract, the determination of ī(t) is of special importance to the holder of American options, put ones in particular.In [4], Chapters 6.3-6.5, and [8], there are several interesting results related to both qualitative and quantitative properties of the optimal exercise boundary.
Except other ones, as for example position of ī(T), the monotonicity of ī(t), upper and lower bounds of ī(t) and convexity of ī(t), the most important are following two ones.First, asymptotic expression of ī(t) near t = T, and at second, the particular equation for determination of ī(t), which is given in [4], Chapter 6.3, Theorem 6.5.
Solution of that equation by successive approximations realized by an adaptive iteration algorithm implemented in sw package Mathematica, there is the core of our paper.
Following [4], Theorem 6.5, the equation for determination of optimal exercise boundary ī(t) for American options is following here q is dividend rate, Ɏ(x) is cumulative distribution function of standard normal distribution N(0,1), parameters ȕ 1 , ȕ 2 , are defined as follows and functions ݃ ݅ ‫;ݐ(‬ ‫,ܭ‬ ܶ, ߪ, ߚ ݅ ), ݄ ݅ ‫;ݐ(‬ Ȟ(ߟ), ߟ, ߪ, ߚ ݅ ), i = 1, 2, are given by following expressions However, it is a nonlinear Volterra integral equation of the second kind.In general, solving such type of equation is a rather challenging problem for numerical mathematics.Regarding (16), we remark that it is applicable in any case of American options when early exercise is interesting, i.e. in case of American put options, and in case of American call options paying dividends.

Solution algorithm
The first step is time discretization.Let t i , i = 0, …, n, such that t 0 = 0, t n = T, denote regular mesh of time points with a mesh size ¨t = T/n.As the ī(t) represents a uncountable set of points (t,S t ) laying on the optimal early exercise boundary ī, then adopting time discretization turns the problem to solve equation ( 16) directly to looking for a finite set of points {t i ,S i }, i = 0, …, n, satisfying a discrete version of ( 16), where we just introduce shorter notation S i instead of cumbersome S t | t=ti .
Time discretization being applied upon (16) yields the following set of fixed-point looking problems here symbol Ȍ expresses formally the right-hand side of (16) after discretization thus denoting its proper dependence upon variable S i , but also upon parameters t i and ī(Ș) thereby making the numerical solution of (19) rather involved.In order to overcome numerically a complication caused by function ī(Ș) appearing in (19), which is itself a solution of equation ( 16) being seeked, we propose the following algorithm based upon successive approximation of ī(Ș) starting with some initial guess.

Adaptive iterative algorithm -steps:
1. Set initial approximation ī 0 (Ș), Ș ‫א‬ [t i ,T], and set index k = 0, 2. Increment current k by 1, for loop Our first attempt to solve (10) for a given t i , which is a core of the step 2 within the proposed algorithm, was simply lured by computational power of sw Mathematica to be able to solve it just by two following commands However, it failed with a Mathematica process control message SystemException["MemoryAllocationFailure"] signalizing an extreme demand of internal Mathematica procedures upon memory allocation when executing such complex command.It was simply over the memory capacity of our Dell Latitude E6510 computer with 4GB RAM.Hence, we turned to break the given algorithm into a sequence of less complex steps and control manually their performance and convergence in adaptive way, too.Following [4], Chapter 6.5, we also know some advanced information about optimal exercise boundary.First, the location point at t = T ī(T) = min(rK/q, K), (20) and second, an asymptotic expression of ī(t) near t = T, in case of q = 0, i.e. for non-dividend paying American put option, As already stated above, we know that ī : t ĺ ī(t) is a convex and monotone increasing function on (0,T).

Numerical example:
For numerical determination of the optimal early exercise boundary ī(t), we select the American put option with following data: strike price K = 40, risk-free rate r = 0.05, underlying asset price S ‫א‬ (0,65], with volatility ı = 0.3, option expiration date T = 1 [year], and dividend rate q = 0, i.e. non-dividend paying case.In Figures 2 and 3, we show two different initial variants of ī 0 (t), constructed either by using asymptotic expression (21) with a proper prolongation, or by ad-hoc way using a cubic polynomial with proper interpolation conditions maintaining ī(T) = K as required by (20), and both qualitative properties, i.e. monotonicity and convexity of ī 0 (t), on (0,T), as well.
In Figure 2, we present an idea of construction of ī 0 (t) using (21) with proper prolongation.A curve given Fig. 2. ī 0 (t) -initial early exercise boundary asymptoticderived with proper prolongation (full line); ī(t) |asymp -curve given by ( 21) on (0,T) (dashed line).Fig. 3. ī 0 (t) -variants: 1) ad-hoc guess cubic curve (full line); 2) asymptotic-derived with proper prolongation one (dashed line).by ( 21) is convex on (0,T), but it violates the monotonicity.Hence, to fix that defect, we propose the following procedure.First, we find a global minimum point t min of function ( 21) on (0,T) issuing the minimum value ī min of that function.Second, we make a smooth prolongation from t = 0 to t = t min by a horizontal line keeping the calculated minimal value ī min , as it is depicted in Figure 2. Finally, we bond both branches together to get a feasible guess for ī 0 (t t T. Assumption q = 0, i.e. the case of non-dividend paying option, simplifies the equation ( 16) for determination of optimal early exercise boundary ī(t) so that the first term is modified and the last one is dropped completely.Hence, we obtain another set of fixed-point looking problems, written formally in the same form (19), however, with a slightly different form of Ȍ.It is defined as the right-hand side of the following equation here parameters ȕ 1 , ȕ 2 , are defined as follows and functions ݃ ݅ ‫;ݐ(‬ ‫,ܭ‬ ܶ, ߪ, ߚ ݅ ), i = 1, 2, are given by (17), and function ݄ 1 ‫;ݐ(‬ Ȟ(ߟ), ߟ, ߪ, ߚ 1 ) by (18), correspondingly.
Classical trapeziodal rule is used for numerical approximation of integral ‫‬ ߱(ߟ)݀ߟ ܶ ‫ݐ‬ , which appears in equation ( 22), as the third term on its right-hand side.Because of the time discretization, we point out that the lower bound t = t i , in particular case, whereas the upper bound T is fixed in any case.
Using function ݄ 1 ‫;ݐ(‬ Ȟ(ߟ), ߟ, ߪ, ߚ 1 ) DQG Ɏ(x), we can define the function Ȧ(Ș) in following way For trapezoidal rule, let m is a number if intervals covering [t i ,T] with a step ¨Ș = (T -t i ) /m, thus providing Now, we formulate the inner-most problem.For given set of points {t i , S i,0 }, find {t i , S i,k }, k 1, such that equation (25) holds, with given tolerance appearing in convergence criterion (26) Now, we present some snippets of our Mathematica code we developed for numerical solution of the problem to find optimal early exercise boundary for American put option.
First, it is a construction of the initial early exercise boundary ī 0 (t) based upon asymptotic expression of ī(t) near t =T, as given by ( 21) Results of our calculations are depicted in Figures 4  and 5. Setting İ 1 = 10 -3 , we computed two variants of optimal early exercise boundary using from two initial functions ī 0 (t), as discussed above and shown in Figure 3, respectively.Fig. 4. Local comparison of computed optimal early exercise boundaries ī(t) computed from two initial curves ī 0 (t) given in Figure 3 (full line and dashed one), and location of asymptoticderived with proper prolongation curve ī 0 (t) (doted line), with set of nodal points {t i ,S i }.Fig. 5. Global comparison of computed optimal early exercise boundaries ī(t) computed from two initial curves ī 0 (t) given in Figure 3 (full line and dashed one), and location of asymptoticderived with proper prolongation curve ī 0 (t) (doted line), with set of nodal points {t i ,S i }.
In Figures 4 and 5, we show also the selected set of points for time discretization, i.e. {t i }, i = 0, ..., n, which are used for basic time discretization of nonlinear Volterra integral equation ( 16), in general.Inspecting the results presented, we may conclude that different initial functions ī 0 (t) cause almost negligible differences between our computed approximations of the optimal early exercise boundary ī(t)., which sounds rather