## 1.INTRODUCTION

Financial derivatives are becoming more and more prevalent all over the world. Among them, an option is a contract which gives the buyer the right, but not the obligation, to buy or sell an underlying asset or instrument at a specified strike price on or before a specified date. American style option, which allows early exercise before the expiration date, is getting popular since it provides an investor with a greater degree of flexibility than a European style option.

However, there are no closed-form solutions for valuing American options, but only a choice of models to approximate the price is available. There have been mainly two approaches on numerical approximations: one is to compute the solution of the Black-Scholes partial differential equations based on finite difference method, projected successive over relaxation method (PSOR), policy or penalty methods. See Higham (2004), Duffy (2006), Kwok (2008). The other is to compute the expectation under risk-neutral measure such as binomial tree method (Cox *et al*., 1979) or Monte Carlo method (Glasserman, 2004; Kwok, 2008). The early exercise feature of American options can be more efficiently handled by using the former approach than the latter method. Although the binomial model is an excellent tool for pricing simple American options, it is vulnerable to extensions to, for example exotic options such as path dependent options or high dimensions. Therefore, we focus on the PDE approach in this work.

The price of American options can be formulated as a linear complementarity problem (LCP), so called the free boundary problem, see Higham (2004). The finite difference method such as *θ*-method with projection for the boundary conditions of American options requires less computational cost than other methods, but it has weakness in accuracy. The PSOR in Gilli *et al*. (2011) is easy to implement but its convergence rate depends on the choice of the relaxation parameter. Its computational cost also increases exponentially as the number of mesh points increases. Ikonen and Toivanen (2004) numerically showed that the operator splitting method, *θ*-method with projec- tion for one dimensional problem, performed better than the PSOR for one dimensional Black-Scholes problems, later Jeong and Kim (2013) extended to higher dimensional problems. The penalty method is another efficient method since the free boundary can be removed by adding a small penalty term and is solved on a fixed domain, see Khaliq *et al*. (2006), Nielsen *et al*. (2002), Zvan *et al*. (1998), Forsyth and Vetzal (2002), Zhang *et al*. (2010). It is being noted that the power penalty method improves the rate of convergence by considering the nonlinear term as power of the usual penalty method in Wang *et al*. (2006) and Sun *et al*. (2015). The policy method has been proposed in Reisinger and Witte (2012), Babbin *et al*. (2014), Forsyth and Labahn (2007) which is a Newtonbased penalty method. Recently Chockalingam and Muthuraman (2015) used a moving boundary method to convert the free boundary problem into a sequence of fixed boundary problems.

A novel idea in this work is to combine the penalty method and the *θ* -method with projection to compute the price of American options efficiently and accurately. We combine the two methods in such a way that the penalty method is used near the expiration date to handle the nonlinear error from the free boundary, where careful computations are required, and then the *θ*-method with projection is used in the remaining contract period to reduce the computational cost. Numerical tests in Section 4 show that for a given error tolerance, the proposed hybrid method approximates the solution with less computational costs than other commonly used methods such as the *θ*- method with projection, PSOR, or policy method.

The outline of the paper is as follows. In Section 2 the problem is first set up to price American options. A new numerical scheme is proposed and compared with other existing numerical methods in Section 3. Section 4 shows the results from the numerical experiments and the conclusions are drawn in Section 5.

## 2.AMERICAN OPTIONS

In this section, we describe the underlying asset model and the linear complementarity problem (LCP). Let *s*(*t*) be the price of the underlying asset at time *t*(0 ≤ *t* ≤ *T*) with a given expiry date *T* following a geometric Brownian motion with a constant interest rate *r* > 0 and a constant volatility *σ* > 0.Then the value *V* (*s*, *t*) of European options under classical Black-Scholes model can be computed by solving the following *Black-Scholes* partial differential equations (PDE)

with the final condition *V* (*s*, *T* ) = Λ(*s*) . Here Λ(*s*) is a payoff function, for instance, Λ(*s*) = maxmax(*K* − *s*, 0) for the American vanilla put option, where *K* is a pre-defined exercise price.

For American options, the early exercise feature makes the Black-Scholes PDE in (1) into the following partial differential inequality, known as a *linear complementarity problem* (LCP)

with conditions

for the payoff function Λ(*s*) and

For American vanilla put option, the boundary conditions are *V* (0, *t*) = *K* from (3) and *V* (*s*, *t*)→0 as *s* → *∞* for all 0 ≤ *t* ≤ *T*. The LCP (2) – (4) can be rewritten as a forward problem by changing the parameters *τ* = *T* − *t* and *v*(*s*, *τ* ) = *V* (*s*, *t*)

For barrier options, which have a payoff that switches on or off depending on whether the asset crosses a predefined barrier level, the same equation (1) for European or (5) for American is satisfied unless the barrier is crossed. For instance, a down-and-out put option has a payoff that is zero if the asset crosses some pre-defined barrier level *B* at some time in [0, *T*] . If the barrier is not crossed then the payoff becomes that of the usual put, Λ(*s*) = max(*K* − *s*, 0) .

## 3.HYBRID PENALTY METHOD

In this section, we review numerical schemes for solving LCP (5) and propose the hybrid penalty method.

### 3.1.*θ*-Method with Projection

The *θ*-method with projection for one dimensional problem was also known as operator splitting method in Ikonen and Toivanen (2004). Let us first consider the computational domain *D* = [0, *s*_{max}]×[0, *T*] where *s*_{max} is sufficiently large value in order to approximate the boundary *v*(*s*_{max} , *t*) = 0 for put options. Then let us discretize the domain *D* into *N _{τ}* uniform intervals for time variable

*τ*and

*N*uniform intervals for space varia ble

_{s}*s*and denote 0 =

*τ*

_{0}<

*τ*

_{1}<…<

*τ*

*N*=

_{τ}*T*and $0={s}_{0}<{s}_{1}<\cdots <{s}_{{N}_{s}}={s}_{\text{max}}$ with $\Delta t=T/{N}_{\tau}$ and $\Delta s={s}_{\text{max}}/{N}_{s}$. Let us denote the approximation of the value of option

*v*(

*s*,

*τ*) in (5) $by{v}_{i}^{k}=v\left({s}_{i},\text{\hspace{0.17em}}{t}_{k}\right)$ for

*i*= 0,1, 2, …,

*N*and

_{s}*k*= 0, 1, 2, …,

*N*.

_{τ}Let us consider the finite difference discretization for the Black-Scholes PDE in (1) to get the difference equation

The formula for *θ*-method in (6) becomes the explicit scheme for *θ* = 0 and the fully implicit scheme if *θ* =1 . The parameter *θ* = 0.5 gives the Crank-Nicolson method, which is second-order accurate in both Δ*t* and Δ*s* . The equation (6) then can be written as

where(8)

and(9)

Since *v*(*s*, *τ* ) = 0 as *s* → *∞* , the boundary condition becomes ${v}_{{N}_{s}}^{k}={v}_{{N}_{s}}^{k-1}=0$ so that (7) can be converted to a matrix form(10)

where

and

From the boundary conditions, the values of parameters at the boundary can be set as(11)

The algorithm of the *θ* -method with projection can be summarized as follows.

### 3.2.Projected SOR Method

Let us consider(13)

for *v*∈*R*^{N +1}, where $A=\left({a}_{ij}\right)\in {R}^{\left(N+1\right)\times \left(N+1\right)}$ and *f* ∈ *R*_{N +1} . Given the decomposition of *A* = *L* + *D* + *U* with a lower triangular matrix, *L* , a diagonal matrix, *D* , and a upper triangular matrix, *U*, the *Gauss-Seidel* method computes *v ^{k}* from the iteration (

*D*+

*L*)

*v*

^{k + 1}= −

*Uv*+

^{k}*f*or

*v*

^{k}

^{+1}= -(

*D*+

*L*)

^{-1}

*Uv*+ (

^{k}*D*+

*L*)

^{-1}

*f*The

*Successive Over Relaxation (SOR)*method updates

*v*by(14)

^{k}

where 0 < *ω* < 2 is a parameter for SOR. Then the LCP (5) is equivalent to(15)

This problem can be solved by *Projected SOR (PSOR)*: algorithm

### 3.3.Policy Method

Given the LCP (5)} above, let us define ${\left({\varphi}^{n}\right)}_{i}$ for 0 ≤ *i* ≤ *N* such that

and

such that

Then the policy iteration method (Reisinger and Witte, 2012; Forsyth and Labahn, 2007) computes *v ^{k}* from

Following stopping criteria is used for *v ^{k}* for

*i*= 0,1, …,

*N*as in Reisinger and Witte (2012)(16)

for a tolerance *ε* and the algorithm for the policy iteration algorithm can be summarized as follows:

### 3.4.Penalty Method

The penalty method Forsyth and Labahn (2007) adds a penalty term to the Black-Scholes PDE(17)

where *ρ* is a penalty parameter. It can be written as(18)

or(19)

where(20)

The penalty method applies generalized Newton iteration method for *v ^{k}* :

and the penalty algorithm can be summarized as follows.

### 3.5.Hybrid Penalty Method

When the numerical methods above are applied to various options, for example American barrier options, it has been observed that the policy and penalty methods show superior accuracy when the number of grids is fixed, while the *θ*-method with projection has advantages in the computational speeds. It is known that the accuracy of option pricing depends on the accuracy of the computation near the expiry. Thus, a hybrid model of the penalty and *θ*-method with projection has been derived(21)

where *τ* = *T* − *t* represents the time to expiry and 0 < *α* <1 is a parameter for the hybrid method. The initial application of the penalty method promptly reduces the nonlinear error from the free boundary, which determines the accuracy of the whole computation. Then the explicit quickly solves the problem for the remaining of the time period. As seen in Section 4, the proposed hybrid method solves efficiently and accurately the linear complementarity problem for the valuation of the American options. The algorithm for the hybrid penalty method is as follows:

## 4.NUMERICAL COMPUTATIONS

### 4.1.American Option

In this section, numerical results of the linear complementarity problem (5) are performed for American vanilla and down-and-out barrier put options. The interest rate *r* = 0.05, the volatility *σ* = 0.25, the strike price *K* = 2 and the time to expiry *T* = 1 are used for the computations. The down-and-out barrier *B* = 1 is used for the barrier option with zero rebate. The solution from the penalty method with *N _{τ}* = 2

^{16}time steps and

*N*= 2

_{s}^{16}price nodes is considered the exact solution for the corresponding American option.

*α*is set to

*α*= 7 / 8 for the hybrid penalty method in this paper.

The Table 1-Table 5 show the values of the American vanilla at-the-money (*s* = *K*) put option prices and delta and corresponding errors from the *θ* -method with projection, PSOR, policy, penalty and hybrid penalty methods, respectively, as the number of nodes *N _{s}* and the number of time step

*N*increase. The CPU times for each

_{τ}*N*and

_{s}*N*are also shown. The computational cost measured by CPU time is minimized with the

_{τ}*θ*-method with projection, but its convergence rate of the error is lower than those of other methods.Table 2Table 3Table 4

### 4.2.American Barrier Option

The Table 6-Table 10 show the values of the American down-and-out barrier at-the-money (*s* = *K*) put option prices and delta and corresponding errors from the *θ*-method with projection, PSOR, policy, penalty and hybrid penalty methods, respectively, as the number of nodes *N _{s}* and the number of time step

*N*increasesTable 7Table 8Table 9

_{τ}Figure 2 compares the errors of the American down-and-out barrier at-the-money put option prices with respect to the CPU times required for the computations. Similarly to the Figure 1, the convergence rate for the *θ*-method with projection (solid) is smaller than those of the others, and given a fixed tolerance rate, the hybrid penalty method (circle) requires less computational time than the *θ*-method with projection (solid), PSOR (star), policy (triangle), or penalty (square) methods.

## 5.CONCLUSIONS

The proposed hybrid penalty method initially applies the penalty method to reduce the nonlinear errors from the linear complementarity problem, then uses the *θ*- method with projection to gain computational speed up. Numerical experiments for American options show the superiority of the method over well-known finite difference schemes for the free boundary problems.

The proposed hybrid method can be also seen as a technique to compute other problems in much more general context. The hybrid method can be extended in the direction of path-dependent options, higher dimensional American options and portfolio optimization problems as well. Also based on real financial market data, we may illustrate the performance of the hybrid method for future research. More in-depth analysis of the convergence rate of the hybrid method will be another research direction.