• Editorial Board +
• For Contributors +
• Journal Search +
Journal Search Engine
ISSN : 1598-7248 (Print)
ISSN : 2234-6473 (Online)
Industrial Engineering & Management Systems Vol.16 No.2 pp.224-233
DOI : https://doi.org/10.7232/iems.2017.16.2.224

# An Efficient Hybrid Penalty Method for Pricing American Options

Hongjoong Kim, Taeyoung Oh, Kyoung-Sook Moon*
Department of Mathematics, Korea University, Korea
Department of Mathematical Finance, Gachon University, Korea
Corresponding Author, ksmoon@gachon.ac.kr
May 17, 2016 December 31, 2016 March 27, 2017

## ABSTRACT

We propose a hybrid numerical method for computing the prices of American options. In order to solve efficiently and accurately the linear complementarity problem arising in the valuation of American options, the proposed method initially applies the penalty method to annihilate the nonlinear error from the free boundary, then performs the θ- method with projection to solve the rest of the problem quickly. Numerical computations show that the proposed hybrid method is more efficient than other existing methods for a given level of accuracy.

## 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 ≤ tT) 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)

$∂ V ∂ t + r s ∂ V ∂ s + 1 2 σ 2 s 2 ∂ 2 V ∂ s 2 − r V = 0$
(1)

with the final condition V (s, T ) = Λ(s) . Here Λ(s) is a payoff function, for instance, Λ(s) = maxmax(Ks, 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)

$∂ V ∂ t + r s ∂ V ∂ s + 1 2 σ 2 s 2 ∂ 2 V ∂ s 2 − r V ≤ 0$
(2)

with conditions

$V ( s , t ) ≥ Λ ( s ) , for all 0 ≤ t ≤ T , s ≥ 0$
(3)

for the payoff function Λ(s) and

$o n e o f ( 2 ) a n d ( 3 ) i s a t e q u a l i t y f o r a l l 0 ≤ t ≤ T , s ≥ 0$
(4)

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 ≤ tT. The LCP (2) – (4) can be rewritten as a forward problem by changing the parameters τ = Tt and v(s, τ ) = V (s, t)

${ L v = ∂ v ∂ τ − 1 2 σ 2 s 2 ∂ 2 v ∂ s 2 − r s ∂ v ∂ s + r v ≥ 0 v ( s , τ ) ≥ Λ ( s ) ( L v ) ⋅ ( v ( s , τ ) − Λ ( s ) ) = 0$
(5)

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(Ks, 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, smax]×[0, T] where smax is sufficiently large value in order to approximate the boundary v(smax , t) = 0 for put options. Then let us discretize the domain D into Nτ uniform intervals for time variable τ and Ns uniform intervals for space varia ble s and denote 0 = τ0 < τ1 <…<τNτ = T and $0 = s 0 < s 1 < ⋯ < s N s = s max$ with $Δ t = T / N τ$ and $Δ s = s max / N s$. Let us denote the approximation of the value of option v(s,τ ) in (5) $b y v i k = v ( s i , t k )$ for i = 0,1, 2, …, Ns and k = 0, 1, 2, …, Nτ.

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

$v i k − v i k − 1 Δ t = θ ( 1 2 σ 2 s i 2 v i + 1 k − 2 v i k + v i − 1 k Δ s 2 + r s i v i + 1 k − v i − 1 k 2 Δ s − r v i k ) + ( 1 − θ ) ( 1 2 σ 2 s i 2 v i + 1 k − 1 − 2 v i k − 1 + v i − 1 k − 1 Δ s 2 + r s i v i + 1 k − 1 − v i − 1 k − 1 2 Δ s − r v i k − 1 )$
(6)

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

$a i k v i − 1 k + b i k v i k + c i k v i + 1 k = A i k v i = 1 k − 1 + B i k v i k − 1 + C i k v i + 1 k − 1$
(7)

where(8)

${ a i k = − 1 2 θ Δ t ( σ 2 s i 2 1 Δ s 2 − r s i 1 Δ s ) b i k = 1 + θ Δ t ( σ 2 s i 2 1 Δ s 2 + r ) c i k = − 1 2 θ Δ t ( σ 2 s i 2 1 Δ s 2 + r s i 1 Δ s )$
(8)

and(9)

${ A i k = 1 2 ( 1 − θ ) Δ t ( σ 2 s i 2 1 Δ s 2 − r s i 1 Δ s ) B i k = 1 − ( 1 − θ ) Δ t ( σ 2 s i 2 1 Δ s 2 + r ) C i k = 1 2 ( 1 − θ ) Δ t ( σ 2 s i 2 1 Δ s 2 + r s i 1 Δ s )$
(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)

$M 1 v k = M 2 v k − 1$
(10)

where

$M 1 = [ b 0 c 0 0 ⋯ 0 a 1 b 1 c 1 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 ⋯ a N s − 1 b N s − 1 c N s − 1 0 ⋯ 0 a N s b N s ]$

$M 2 = [ B 0 C 0 0 ⋯ 0 A 1 B 1 C 1 ⋯ 0 0 ⋱ ⋱ ⋱ 0 0 ⋯ A N s − 1 B N s − 1 c C − 1 0 ⋯ 0 A N s B N s ]$

and

$v k = [ v 0 k , v 1 k , ⋯ , v N s k ] T , k = 0 , 1 , 2 , ⋯ , N τ$

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

$b N s = 1 , B N s = 0 , a N s = 0 , A N s = 0$
(11)

and the LCP (5) becomes(12)

${ M 1 v k − M 2 v k − 1 ≥ 0 v k ≥ Λ ( M 1 v k − M 2 v k − 1 ) ⋅ ( v k − Λ ) = 0$
(12)

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

### 3.2.Projected SOR Method

Let us consider(13)

$A v = f$
(13)

for vRN +1, where $A = ( a i j ) ∈ R ( N + 1 ) × ( N + 1 )$ and fRN +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 vk from the iteration (D + L)vk + 1 = −Uvk + f or vk+1 = -(D + L)-1Uvk + (D + L)-1f The Successive Over Relaxation (SOR) method updates vk by(14)

$v k + 1 = v k + ω ( v k + 1 − v k ) = v k + ω ( − D − 1 L v k + 1 − D − 1 U v k + D − 1 f − v k ) = v k + ω ( − 1 ω ∑ j = 1 i − 1 a i j v j k + 1 − 1 a i i ∑ j = i + 1 N + 1 a i j v j k + 1 a i i f i )$
(14)

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

$min { M 1 v k − M 2 v k − 1 , v k − Λ } = 0 min { v − A − 1 f , v − Λ } = 0 ⇔ where A = M 1 f = M 2 v k − 1 v = v k ⇔ v = max { A − 1 f , Λ }$
(15)

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

### 3.3.Policy Method

Given the LCP (5)} above, let us define $( ϕ n ) i$ for 0 ≤ iN such that

$( ϕ n ) i ∈ argmin ϕ ∈ ( 0 , 1 ) { ϕ ( M 1 v k − 1 − f ) i + ( 1 − ϕ ) ( v k − 1 − Λ ) i }$

and

$A i = ( ϕ n ) i ( M 1 ) i + ( 1 − ( ϕ n ) i ) ( I N s ) i f i = ( ϕ n ) i ( f ) i + ( 1 − ( ϕ n ) i ) ( Λ ) i$

such that

$( A v ( k − 1 ) − f ) i = = min ϕ ∈ ( 0 , 1 ) { ϕ ( M 1 v k − 1 − f ) i + ( 1 − ϕ ) ( v k − 1 − Λ ) i }$

Then the policy iteration method (Reisinger and Witte, 2012; Forsyth and Labahn, 2007) computes vk from

$A v k = f$

Following stopping criteria is used for vk for i = 0,1, …, N as in Reisinger and Witte (2012)(16)

$[ M 1 v k − b ‖ b ‖ ∞ ≥ − ∈ ] ∧ [ v k − Λ ‖ Λ ‖ ∞ ≥ − ∈ ] ∧ [ ( M 1 v k − b ) i ‖ b ‖ ∞ ≤ ∈ ∧ ( v k − Λ ) ‖ Λ ‖ ∞ ≤ ∈ ]$
(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)

$∂ v ∂ τ = 1 2 σ 2 s 2 ∂ 2 v ∂ s 2 + r s ∂ v ∂ s − r v + ρ max ( Λ − v , 0 )$
(17)

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

$M 1 v k = M 2 v k − 1 + P ¯ ( v k ) ( g − v k )$
(18)

or(19)

$[ M 1 + P ¯ ( v k ) ] v k = M 2 v k − 1 + P ¯ ( v k ) Λ$
(19)

where(20)

$P ¯ ( v k ) i j = { ρ if v i k < Λ i and i = j 0 otherwise$
(20)

The penalty method applies generalized Newton iteration method for vk :

$( v k ) m + 1 = ( v k ) m − [ F ( v k ) − 1 ] [ F ( v k ) ] = ( v k ) m − [ M 1 + P ¯ ( v k ) ] − 1 × [ ( M 1 + P ¯ ( v k ) ) v k − M 2 v k − 1 − P ¯ ( v k ) Λ ] − 1 = [ M 1 + P ¯ ( v k ) ] − 1 [ M 2 v k − 1 + P ¯ ( 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)

(21)

where τ = Tt 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τ = 216 time steps and Ns = 216 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 Ns and the number of time stepNτ increase. The CPU times for each Ns and 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 Ns and the number of time stepNτ 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.

## ACKNOWLEDGMENT

This work was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by Ministry of education (NRF-2014R1A1A2058271); and by the Gachon University research fund of 2016.

## Figure

θ -method with projection

PSOR method for the American option

Policy method for the American option

Penalty method for the American option

Penalty method for the American option

The errors of the option prices vs CPU time for the American vanilla put option with r = 0.05, σ = 0.25, K = 2, T =1 or the θ-method with projection (solid), PSOR (star), policy (triangle), penalty (square) and hybrid penalty (circle) methods.

The errors of the option prices vs CPU time for the American down-and-out barrier put option with r = 0.05, σ = 0.25, K = 2, T =1, B =1 for the hybrid penalty method with α = 7 / 8.

## Table

The option prices and delta and corresponding errors from the θ-method with projection for the American vanilla put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the PSOR method for the American vanilla put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the policy method for the American vanilla put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the penalty method for the American vanilla put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the hybrid penalty method with α = 7 / 8 for the American vanilla put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the θ-method with projection for the American downand- out put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the PSOR method for the American down-and-out put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the policy method for the American down-and-out put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the penalty method for the American down-and-out put option with r = 0.05,σ = 0.25, K = 2, T =1

The option prices and delta and corresponding errors from the hybrid penalty method with α = 7 / 8 for the American down-and-out put option with r = 0.05,σ = 0.25, K = 2, T =1

## REFERENCES

1. Babbin J , Forsyth P A , Labahn G (2014) A comparison of iterated optimal stopping and local policy iteration for American options under regime switching , Journal of Scientific Computing, Vol.58 (2) ; pp.409-430
2. Chockalingam A , Muthuraman K (2015) An approximate moving boundary method for American option pricing , European Journal of Operational Research, Vol.240 (2) ; pp.431-438
3. Cox J C , Ross S A , Rubinstein M (1979) Option pricing: A simplified approach , Journal of Financial Economics, Vol.7 (3) ; pp.229-263
4. Duffy D (2006) Finite Difference Methods in Financial Engineering: A Partial Differential Equation Approach, John Wiley & Sons,
5. Forsyth P A , Labahn G (2007) Numerical methods for controlled Hamilton-Jacobi Bellman PDEs in finance , Journal of Computational Financ, Vol.11 ; pp.1-44
6. Forsyth P A , Vetzal K R (2002) Quadratic convergence for valuing American options using a penalty method , SIAM Journal on Scientific Computing, Vol.23 (6) ; pp.2095-2122
7. Glasserman P (2004) Monte Carlo Methods in Financial Engineering, Springer,
8. Gilli M , Maringen D , Schumann E (2011) Numerical Methods and Optimization in Finance, Academic Press,
9. Higham D J (2004) An Introduction to Financial Option Valuation, Cambridge University Press,
10. Ikonen S , Toivanen J (2004) Operator splitting methods for American option pricing , Applied Mathematics Letters, Vol.17 (7) ; pp.809-914
11. Jeong D , Kim J (2013) A comparison study of ADI and operator splitting methods on option pricing models , Journal of Computational and AppliedMathematics, Vol.247 (1) ; pp.162-171
12. Khaliq A Q M , Voss D A , Kazmi S H K (2006) A linearly implicit predictor-corrector schemefor pricing American options using penalty method approach , Journal of Banking and Finance, Vol.30 (2) ; pp.489-502
13. Kwok Y K (2008) Mathematical Models of Financial Derivatives, Springer-Verlag,
14. Nielsen B F , Skavhaug O , Tveito A (2002) Penalty and front-fixing methods for the numerical solution of American option problems , Journal of Computational Finance, Vol.5 ; pp.65-98
15. Reisinger C , Witte J H (2012) On the use of policy iteration as an easy way of pricing American options , SIAM Journal on Financial Mathematic , Vol.3 (1) ; pp.459-478
16. Sun Z , Liu Z , Yang X (2015) On power penalty methods for linear complementarity problems arising from American option pricing , Journal of GlobalOptimization, Vol.63 (1) ; pp.165-180
17. Wang S , Yang X Q , Teo K L (2006) Power penalty method for a linear complementarity problem arising from American option valuation , Journal of Optimization Theory and Applications, Vol.129 (2) ; pp.227-254
18. Zhang K , Yang X Q , Wang S , Teo K L (2010) Numerical performance of penalty method for Americanoption pricing , Optimization Methods and Software, Vol.25 (5) ; pp.737-752
19. Zvan R , Forsyth P A , Vetzal K R (1998) Penalty methods for American options with stochastic volatility , Journal Computational and Applied Mathematics, Vol.91 (2) ; pp.199-218