A Simple Numerical Method for Pricing an American Put Option

Abstract

We present a simple numerical method to find the optimal exercise boundary in an American put option. We formulate an intermediate function with the fixed free boundary that has Lipschitz character near optimal exercise boundary. Employing it, we can easily determine the optimal exercise boundary by solving a quadratic equation in time-recursive way. We also present several numerical results which illustrate a comparison to other methods.

1. Introduction

The owner of a put (call) option has the right but no obligation to sell (buy) an underlying asset at the exercise price. European options can be exercised only on the expiry date, while American options can be exercised at any time until the expiry date. Closed-form solutions for the European options are derived in papers by Black and Scholes [1] and Merton [2]. In the case of American options, because of the early exercise possibility, the pricing problem leads to complications for analytic calculation. 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, mostly called an optimal exercise boundary. Hence, finance researchers have studied methods to quickly and accurately find the optimal exercise boundary. These methods are basically of two types, that is, analytical approximations such as those developed by Geske and Johnson [5], MacMillan [6], Barone-Adesi and Whaley [7], and Ju [8] and numerical methods such as those of Brennan and Schwartz [9], Hull and White [10], and Longstaff and Schwartz [11]. Zhu [12] finds an exact and explicit solution of the Black-Scholes equation for the valuation of American put options using Taylor series with infinitely many terms. His work is an excellent result for the valuation of American put options; however, it seems difficult to perform his solution numerically. The infinite sum is likely to yield many computation errors. Zhao and Wong [13] study an extension of Zhu’s work [12] to price American options under general diffusion processes.

The majority of numerical methods for pricing American options, such as the finite difference method of Brennan and Schwartz [9], the binomial method of Cox et al. [14], the Monte Carlo simulation method of Grant et al. [15], the least squares method of Longstaff and Schwartz [11], the integral-equation method of Ševčovič [16], and the Laplace transform method of Zhu [17], are time-recursive ways. Their idea is to discretize the lifetime of an option and find its optimal exercise boundary backward in time. Since time-recursive ways yield repeated calculations for every time step, they require fast computation times and small pricing errors. Also, front-fixing methods developed by Wu and Kwok [18] and Nielsen et al. [19] apply a nonlinear transformation to fix the boundary and solve the resulting nonlinear problem. A secant method developed by Zhu et al. [20] needs to solve a nonlinear problem, and a moving boundary approach developed by Muthuraman [21] converts the arising linear free boundary partial differential equation (PDE) problem into a sequence of linear fixed-boundary PDE problems. More recently, Zhu and Zhang [22] introduced a new predictor-corrector scheme to price American put options under the Black-Scholes model, and then Zhu and Chen [23] proposed an extension of Zhu and Zhang’s work [22] to solve for the valuation of American put options with stochastic volatility model.

The main contribution of this paper is the development of a simple numerical method to find optimal exercise boundary in a time-recursive way. Our result is motivated by the necessity for better understanding of the solution surface near optimal exercise boundary. We adopt the front-fixing transformation [18] to change the unknown free boundary to a known and fixed boundary. We exploit an intermediate function with the fixed free boundary that has Lipschitz character which avoids the degeneracy of the solution surface near optimal exercise boundary as in Kim et al. [24]. Indeed, our function from the Black-Scholes equation and the boundary conditions transforms the surface above the exercise region onto a new Lipschitz surface which forms a sufficiently large angle with the hyperplane corresponding to the exercise region, thereby making the borderline more easily distinguishable (see Figure 2). We use implicit scheme in the continuation region and apply extrapolation near optimal exercise boundary. Thus we can determine the optimal exercise boundary by solving a quadratic equation in a time-recursive way. Our method also provides fast and accurate results for calculating the optimal exercise boundary and pricing American put options.

The structure of the paper is as follows. Section 2 presents the model formulation. The intermediate function with the fixed free boundary to calculate the optimal exercise boundary is presented in Section 3. Numerical results and comparative studies are presented in Section 4. Section 5 summarizes the paper.

2. Problem Formulation

In this section, we present a mathematical formula for pricing an American put option.

Consider an American put option on an underlying asset (stock) with exercise price 𝐾 and expiration 𝑇. In risk-neutral probability, an underlying asset price 𝑆(𝑡) is governed by the following stochastic differentiable equation: 𝑑𝑆(𝑡)=(𝑟−𝛿)𝑆(𝑡)𝑑𝑡+𝜎𝑆(𝑡)𝑑𝑊(𝑡), (1) where 𝑟>0 represents the risk-free interest rate, 𝛿≥0 represents the continuous dividend yield, 𝜎>0 represents the volatility of the underlying asset price, and 𝑊(𝑡) is the standard Brownian motion. The payoff function of the put option at 𝑇 is defined as (𝐾−𝑆(𝑇))⁺=max{𝐾−𝑆(𝑇),0}. (2) The valuation of an American put option is denoted by 𝑃(𝜏,𝑆), where 𝜏(∶=𝑇−𝑡) is the time to expiration for 𝜏∈[0,𝑇] and 𝑆 is the underlying asset price for 𝑆∈[0,∞).

As seen in the previous article by McKean [3], the valuation of an American put option is considered the solution to a free boundary problem with a parabolic PDE. We suppose that the optimal exercise boundary 𝛽(𝜏) is continuously nonincreasing with 𝛽(0)=min{𝐾,𝐾𝑟/𝛿}. The region where it is optimal to hold, generally called the continuation region, is defined as Ω_𝑐=[0,𝑇]×(𝛽(𝜏),∞), and the region where it is optimal to exercise, generally called the exercise region, is defined as Ω_𝑒=[0,𝑇]×[0,𝛽(𝜏)]. Then, 𝑃(⋅,⋅) and 𝛽(⋅) uniquely solve ℒ𝑃−𝑟𝑃−𝑃_𝜏=0inΩ_𝑐,𝑃(0,𝑆)=(𝐾−𝑆)⁺,𝑃𝜏,𝛽(𝜏)=𝐾−𝛽(𝜏),lim_𝑆↑∞max_{0≤𝜏≤𝑇}|𝑃(𝜏,𝑆)|=0,lim_{(𝜏)𝑆↓𝛽}𝑃_𝑠()𝑃()𝜏,𝑆=−1,𝜏,𝑆=𝐾−𝑆inΩ_𝑒,𝑃(𝜏,𝑆)≥(𝐾−𝑆)⁺, (3) where ℒ𝑃=(1/2)𝜎²𝑆²𝑃_𝑠𝑠+(𝑟−𝛿)𝑆𝑃_𝑠 and 𝑃_𝑠,𝑃_𝑠𝑠, and 𝑃_𝜏 are defined by the infinitesimal generator and partial derivatives, respectively. Here, we assume 𝑟>𝛿 for using the Black-Scholes equation at 𝜏=0. See, for example, Karatzas and Shreve [25] for general reference.

Figure 1 shows an illustration of an optimal exercise boundary with 𝑃. The two regions are separated by the optimal exercise boundary. From Figure 1 it is necessary that 𝑃 must satisfy max_𝜏,𝑆{ℒ𝑃−𝑟𝑃−𝑃_𝜏,(𝐾−𝑆)⁺−𝑃}=0. This condition is known as the Hamilton-Jacobi Bellman (HJB) equation. As Chockalingam and Muthuraman [26] point out, the continuation and exercise regions are determined by which term in the HJB equation is tight. Their method requires iterations till convergence of the boundaries. However, we emphasize that there is no iteration in our method using a square root transformation. Refer to Chockalingam and Muthuraman [26], and Pham [27].

Figure 1 

Boundary conditions.

(a)

(b)

Figure 2 

Transformation from 𝒫(𝜏,𝑢) to 𝑄(𝜏,𝑢).

A front-fixing method, proposed in Wu and Kwok [18], uses a change in variables to transform the free boundary problem into a nonlinear problem on a fixed domain. The following transformation of state variable serves for such a purpose: 𝑢=ln𝑆−ln𝛽(𝜏). (4) They derive the equation and the boundary conditions with respect to 𝑢 as follows: ℒ𝒫−𝑟𝒫−𝒫_𝜏(),𝒫()(),𝒫()(𝜏),𝒫=0,𝑢∈0,∞0,𝑢=0,𝑢∈0,∞𝜏,0=𝐾−𝛽_𝑢()(𝜏),𝒫()𝜏,0=−𝛽𝜏,∞=0, (5) where ℒ𝒫=(1/2)𝜎²𝒫_𝑢𝑢+(𝑟−𝛿−(𝜎²/2)+(𝛽^′(𝜏)/𝛽(𝜏)))𝒫_𝑢 and 𝒫_𝑢 and 𝒫_𝑢𝑢 are defined by infinitesimal generator and partial derivatives, respectively. Note that 𝛽(𝜏) is a monotonically decreasing function of 𝜏 with a nontrivial asymptotic limit as follows: 𝛽(∞)=𝜃𝜃−1𝐾, (6) where −()𝜎𝜃=𝑟−𝛿−1/2²−𝑟−𝛿−(1/2)𝜎²²+2𝜎²𝑟𝜎². (7) Namely, the optimal exercise boundary does not change with time. Especially, plugging 𝛿=0 into (6), we have the asymptotically optimal exercise boundary as follows: 𝛽(∞)=𝛾𝛾+1𝐾, (8) where 𝛾=2𝑟/𝜎². Note that transformation (4) is valid because 𝛽(𝜏)>0 holds for all 𝜏≥0. Refer to Kim [28].

3. Intermediate Function with the Fixed Free Boundary

In this section, we present an intermediate function with the fixed free boundary and can determine the optimal exercise boundary by solving a quadratic equation in a time-recursive way.

Under the assumption of the Black-Scholes model, the time for optimal exercise can be shown to be the first hitting time of a boundary, the optimal exercise boundary, in the plane consisting of pairs of the underlying asset price and the time to expiration. Namely, the price curve of an American put option touches the line representing the intrinsic value tangentially. With a careful examination of the solution surface near optimal exercise boundary, we find a Lipschitz surface which avoids the degeneracy of the solution surface near optimal exercise boundary. To find the optimal exercise boundary, we present an intermediate function with the fixed free boundary that has Lipschitz character near optimal exercise boundary as follows: 𝑄(𝜏,𝑢)∶=𝒫(𝜏,𝑢)−𝐾+𝑒^𝑢𝛽(𝜏). (9) The transformed function 𝑄(𝜏,𝑢) provides that the solution surface in Ω_𝑒 is a horizontal plane, and it is an inclined plain in Ω_𝑐. Namely, this function forms a sufficiently large angle with the hyperplane corresponding to the exercise region, thereby making the borderline more easily distinguishable. 𝑄(𝜏,𝑢) also has a Lipschitz character with nonsingularity in (0,∞) and a nondegeneracy property near optimal exercise boundary. Hence, we have 𝑄()(∞)(𝜏),𝑄()().𝜏,𝑢=0if𝑢∈ln𝛽−ln𝛽,0𝜏,𝑢>0if𝑢∈0,∞ (10)Figure 2(a) shows that 𝒫(𝜏,𝑢) is transformed to 𝑄(𝜏,𝑢), and Figure 2(b) is a magnified view of the optimal exercise boundary.

We find the intermediate function with the fixed free boundary to decide the optimal exercise boundary by the Taylor series. From 𝒫(𝜏,𝑢)=𝑄²(𝜏,𝑢)+𝐾−𝑒^𝑢𝛽(𝜏), we obtain the following relations near optimal exercise boundary (𝑢=0): 𝒫_𝑢𝑢=2𝑄²_𝑢−𝛽,𝒫_𝑢=−𝛽,𝒫=𝐾−𝛽,𝒫_𝜏=−𝛽^′. (11) Plugging (11) into (5), we obtain 12𝜎²2𝑄²_𝑢+1−𝛽𝑟−𝛿−2𝜎²+𝛽^′𝛽−𝛽−𝑟𝐾−𝛽+𝛽^′=0, (12) and then we get 𝑄²_𝑢=(𝑟𝐾−𝛿𝛽)/𝜎². More precisely, we have 𝑄_𝑢()=√𝜏,0𝑟𝐾𝜎,√𝑟𝜎√𝐾<𝑄_𝑢()<√𝜏,0𝑟𝐾𝜎𝑒^𝑢∗ (13) with 𝛿=0, where 𝑢^∗=ln𝛽(∞)−ln𝛽(𝜏)<0. Then, we have |𝒫−𝐾+𝑒^𝑢𝛽(𝜏)|≤𝑐|𝑢−0|² for some constant 𝑐 because 𝑄(𝜏,𝑢) is Lipschitz and a natural candidate for computation in (0,∞). We obtain an angle between exercise surface (𝑄=0) and 𝑄 surface (𝑄>0) such that 0<𝜆₀<𝜕𝑄/𝜕𝑢<𝜆₁ for some constants 𝜆₀ and 𝜆₁. We also calculate the partial derivative with respect to 𝑢 in (5) as follows: 12𝜎²𝒫_𝑢𝑢𝑢+1𝑟−𝛿−2𝜎²+𝛽^′𝛽𝒫_𝑢𝑢−𝑟𝒫_𝑢−𝒫_𝜏𝑢=0. (14) From 𝒫(𝜏,𝑢)=𝑄²(𝜏,𝑢)+𝐾−𝑒^𝑢𝛽(𝜏), we obtain the following relations near optimal exercise boundary (𝑢=0): 𝒫_𝑢𝑢𝑢=6𝑄_𝑢𝑄_𝑢𝑢−𝛽,𝒫_𝜏𝑢=−𝛽^′. (15) Hence, plugging (15) into (14), we get 𝑄_𝑢𝑢=−2𝜉𝜂3𝜎³−𝛿𝛽,3𝜎𝜂 (16) where 𝜉=𝑟−𝛿−𝜎²/2+𝛽^′/𝛽 and 𝜂=𝑟𝐾−𝛿𝛽. From (6) we easily show that 𝜉 and 𝜂 are bounded parameters such that 𝑐₁≤𝜉≤𝑐₂ for some negative constants 𝑐₁ and 𝑐₂ and (𝑟−𝛿)𝐾≤𝜂≤√𝑟𝐾, respectively.

Using the similar arguments, we can obtain the following equations near optimal exercise boundary (𝑢=0) as follows: 12𝜎²𝒫_{𝑢𝑢𝑢𝑢}+1𝑟−𝛿−2𝜎²+𝛽^′𝛽𝒫_𝑢𝑢𝑢−𝑟𝒫_𝑢𝑢−𝒫_𝜏𝑢𝑢𝒫=0,_{𝑢𝑢𝑢𝑢}=6𝑄²_𝑢𝑢+8𝑄_𝑢𝑄_𝑢𝑢𝑢−𝛽,𝒫_𝜏𝑢𝑢=−2𝛿𝜎²+1𝛽^′. (17)(18) Hence, plugging (18) into (17), we get 𝑄_𝑢𝑢𝑢=2𝜉²𝜂3𝜎⁵+𝜉𝛿𝛽6𝜎³𝜂−𝛿²𝛽²12𝜎𝜂³−𝛿𝛽+4𝜎𝜂𝑟𝜂2𝜎³−𝛿𝛽^′2𝜎³𝜂. (19)

Furthermore, we recognize that 𝒫 is analytic up to the optimal exercise boundary and 𝑢 is locally analytic. Hence, the approximation for 𝑄(𝜏,𝑢) at 𝛽(𝜏) can be written as 𝑄(𝜏,𝑢)=𝑄(𝜏,0)+𝑄_𝑢()1𝜏,0𝑢+𝑄2!^𝑢𝑢(𝜏,0)𝑢²+1𝑄3!^𝑢𝑢𝑢(𝜏,0)𝑢³+𝒪𝑢⁴. (20) We introduce the equilibrium parameter 𝑢>0 which enables us to adjust the location of optimal exercise boundary in a mesh size. So, plugging 𝑄_𝑢(𝜏,0), 𝑄_𝑢𝑢(𝜏,0), and 𝑄_𝑢𝑢𝑢(𝜏,0) into (20), we obtain 𝑄(𝜏,𝑢) as follows: 𝑄(𝑢)=𝜂𝜏,𝜎1𝑢−22𝜉𝜂3𝜎³+𝛿𝛽𝑢3𝜎𝜂²+162𝜉²𝜂3𝜎⁵+𝜉𝛿𝛽6𝜎³𝜂−𝛿²𝛽²12𝜎𝜂³−𝛿𝛽+4𝜎𝜂𝑟𝜂2𝜎³−𝛿𝛽^′2𝜎³𝜂𝑢³. (21) We rewrite (21) with respect to 𝜉 as follows: 𝑎𝛽^′2+𝑏̃𝛽^′+𝑐̃=0, (22) where 𝜂𝑢𝑎=³9𝜎⁵𝛽²,𝑏̃=−𝜂𝑢²3𝜎³𝛽+2𝜂𝜈𝑢³9𝜎⁵𝛽+𝛿𝑢³36𝜎³𝜂−𝛿𝑢³12𝜎³𝜂,̃(𝑢)+𝜂𝑢𝑐=−𝑄𝜏,𝜎−13𝜂𝜈𝜎³+𝛿𝛽𝑢2𝜎𝜂²+162𝜂𝜈²3𝜎⁵+𝛿𝜈𝛽6𝜎³𝜂−𝛿²𝛽²12𝜎𝜂³−𝛿𝛽+4𝜎𝜂𝑟𝜂2𝜎³𝑢³,𝛽𝜉=𝜈+^′𝛽,𝜎𝜈=𝑟−𝛿−²2. (23) Combining (22) with 𝛽^′<0, we have 𝛽^′=−𝑏̃−𝑏̃²−4𝑎𝑐̃2𝑎. (24) For discretization (Δ𝜏,Δ𝑢), we introduce a two-dimensional mesh in the first quadrant of the 𝜏−𝑢 plane. From (24) we have 𝛽_𝑛+1−𝛽_𝑛=−̂Δ𝜏𝑏−𝑏̂²−4𝑎𝑐̂2𝑎, (25) where ,𝑄𝑛=0,1,2,…,𝑁,𝑚=0,1,2,…,𝑀,𝑢=𝜌Δ𝑢𝜌>0_𝑛𝑧(),𝛽isthenumericalapproximationto𝑄𝑛Δ𝜏,𝑧Δ𝑢𝑧=𝜌,𝑚_𝑛(),𝜂𝑢isthenumericalapproximationto𝛽𝑛Δ𝜏𝑎=³9𝜎⁵𝛽²,𝑏̂=−𝜂𝑢²3𝜎³𝛽+2𝜂𝜈𝑢³9𝜎⁵𝛽+𝛿𝑢³36𝜎³𝜂−𝛿𝑢³12𝜎³𝜂,̂(𝑢)+𝜂𝑢𝑐=−𝑄𝜏,𝜎−13𝜂𝜈𝜎³+𝛿𝛽𝑢2𝜎𝜂²+162𝜂𝜈²3𝜎⁵+𝛿𝜈𝛽6𝜎³𝜂−𝛿²𝛽²12𝜎𝜂³−𝛿𝛽+4𝜎𝜂𝑟𝜂2𝜎³𝑢³, (26) respectively. Hence, we rewrite (25) with respect to 𝛽_𝑛 as follows: 𝛽_𝑛+1=𝛽_𝑛+−𝑏̂−𝑏̂²−4𝑎𝑐̂2𝑎Δ𝜏. (27) When the initial values are given by 𝑄⁰_𝜌 (transformed price of the American put option) and 𝛽₀ (optimal exercise boundary) at 𝜏=0, we can determine 𝛽₁ (optimal exercise boundary) at Δ𝜏 using (27). More importantly, for updating the optimal exercise boundary our method dose not include iteration until sufficient accuracy is obtained. So, we repeat the previously mentioned process until 𝑁Δ𝜏 and obtain the optimal exercise boundary in a time-recursive way.

4. Numerical Examples

In this section, we provide numerical examples to illustrate our method. We also make runtimes and computation errors compared with the results obtained by other numerical methods such as the binomial method (Binomial) developed by Cox et al. [14], the front-fixing method (Front-fixing) developed by Wu and Kwok [18], and the finite difference implementation of the moving boundary method (MBM-FDM) developed by Muthuraman [21].

All implementations are carried out using a C++ implementation with the a 2.66 Ghz Intel 4 Core CPU with 3 GB RAM. A finite difference method with Crank-Nicolson scheme is proposed for our method. The benchmark results are obtained using the Binomial with 10,000 time steps, and we consider these results to be the exact values of the American put options. Here, root mean squared error (RMSE) is calculated by the values of the Binomial.

The parameter values used to calculate the optimal exercise boundary and values of the American put options are 𝑟=0.05,  𝜎=0.20,  𝛿=0,  𝐾=100,  𝑇=1,  𝑢∈[0,2], and a discrete mesh of 2000×300 nodes.

In Figure 3, we find a numerical optimization 𝜌=5.5702. Table 1 also shows the results of the optimal exercise boundary and the values of American put options with safety parameters. One can see from Table 1 that the optimal exercise boundary monotonically increases as the value of 𝜌 increases, but the value of American put options monotonically decreases when 𝜌 increases. They are so gradual that they are not very susceptible to change in the value of 𝜌.

Table 1 

Valuation of the American put options.

Figure 3 

RMSE of the safety parameter 𝜌.

Table 2 reports the values of the American put options for the specific parameter set associated with the table. In Figure 4 and Table 3, we take the parameter values used in Figure 3 except for the discrete mesh and plot runtimes and computational errors compared with various methods. Note that the discrete meshes of 125×25,  250×50,  500×100,  1000×200, and 2000×400 nodes are plotted in Figure 4 and Table 3.

Table 2 

Comparison of the values of the American put options.

Table 3 

Comparison of RMSE and runtime.

Figure 4 

Comparison of RMSE and runtime.

As is shown in Table 2, Figure 4, and Table 3, although four different methods have similar values of the American put option, our method is computationally faster and more accurate than other methods. Especially, Figure 4 and Table 3 show the numerical convergence of our method. So our method is superior to the others in accuracy and computational efficiency.

5. Final Remarks

The front-fixing method suggested by Wu and Kwok [18] shows a degeneracy near optimal exercise boundary, while our method adopts a square root function to avoid the quadratic behavior of solution surface that causes degeneracy. Our method employing an intermediate function with the fixed free boundary solves a nonlinear problem on a fixed domain derived from a free boundary problem. Since the computation process depends on Lipschitz surface, we need to focus on the motion of the solution surface which would be simple to see the minute behavior of solution surface. The moving boundary approach developed by Mutheraman [21] requires iterations till the convergence of the boundaries. However, we emphasize that there is no iteration in our method. In such a rapidly changing environment, our straightforward method is a very powerful tool to understand financial market. Numerical study also shows that overall speed and accuracy comparisons have demonstrated the superiority of our method over other methods. Our method can be easily extended to other models under stochastic volatility and jump diffusion processes. These remain as topics for future research.

Acknowledgment

This work was supported by the National Research Foundation (NRF-2011-0028951) to the first author and third author.