Document [original]

Fixed Domain Transformations and Split–Step Finite

Difference Schemes for Nonlinear Black–Scholes

Equations for American Options

Julia Ankudinova and Matthias Ehrhardt

Institut für Mathematik, TU Berlin, Strasse des 17. Juni 136, 10623 Berlin, Germany

Abstract. Due to transaction costs, illiquid markets, large investors or risks from an unprotected

portfolio the assumptions in the classical Black–Scholes model become unrealistic and the model

results in strongly or fully nonlinear, possibly degenerate, parabolic diffusion–convectionequations,

where the stock price, volatility, trend and option price may depend on the time, the stock price or

the option price itself.

In this chapter we will be concerned with several models from the most relevant class of nonlin-

ear Black–Scholes equations for American options with a volatility depending on different factors,

such as the stock price, the time, the option price and its derivatives.

We will analytically approach the option price by following the ideas proposed by Ševˇ

coviˇ

c and

transforming the free boundary problem into a fully nonlinear nonlocal parabolic equation defined

on a fixed, but unboundeddomain. Finally, we will present the results of a split–step finite difference

scheme for various volatility models including the Leland model, the Barles and Soner model and

the Risk adjusted pricing methodology model.

1 Introduction

The strong interest in pricing financial derivatives – among them in pricing options – arises

from the fact that financial derivatives, also called contingent claims, can be used to mini-

mize losses caused by price fluctuations of the underlying assets. This process of protection

is called hedging. There is a variety of financial products on the market, such as futures,

forwards, swaps and options. In this chapter we will focus on American Call options.

We recall that an American Call option is a contract where at any time before a pre-

scribed time in the future, known as the expiry date T, the owner of the option, known as

the holder, may purchase a prescribed asset, known as the underlying asset S(t), for a pre-

scribed amount, known as the exercise or strike price K. The opposite party, or the writer,

has the obligation to sell the asset if the holder chooses to excercise his right. The value of

the American Call option at the time of execution, known as the pay-off function, is

V(S, t) = (S−K)+.

Option pricing theory has made a great leap forward since the development of the

Black–Scholes option pricing model by Black and Scholes [6] in 1973 and previously by

2 Julia Ankudinova and Matthias Ehrhardt

Merton [44]. The solution of the famous (linear) Black–Scholes equation

0 = Vt+1

2σ2S2VSS +rSVS−rV, 0< S < Sf(t), t ∈(0, T),(1)

where Vdenotes the value of the option and rthe riskless interest rate, provides both an

option pricing formula for an American Call option and a hedging portfolio that replicates

the contingent claim. This is true under the assumption that the market is complete, which

means that any derivative and any asset can be replicated or hedged with a portfolio of other

assets in the market, cf. [61].

However, this assumption of a complete market is never fulfilled in reality. Due to

transaction costs [4,8,41], large investor preferences [23,24,50] and incomplete markets

[55] the classical model results in strongly or fully nonlinear, possibly degenerate, parabolic

convection–diffusion equations, where both the volatility σand the drift µcan depend on

the time t, the stock price Sor the derivatives of the option price Vitself. Here, we will

be concerned with several transaction cost models from the most relevant class of nonlinear

Black–Scholes equations for American options with a constant drift µand a nonconstant

modified volatility function

eσ2:= eσ2(t, S, VS, VSS).

Under these circumstances (1) becomes the following nonlinear Black–Scholes equa-

tion:

0 = Vt+1

2eσ2(t, S, VS, VSS)S2VSS +rSVS−rV, S > 0, t ∈(0, T)(2)

Studying (1) for an American Call option would be redundant, since the value of an

American Call option equals the value of a European Call option if no dividends are paid

and the volatility is constant. In order to make the model more realistic, we will consider a

modification of (2) for American options, where Spays out a continuous dividend qSdt in

a time step dt:

0 = Vt+1

2eσ2(t, S, VS, VSS)S2VSS + (r−q)SVS−rV, S > 0, t ∈(0, T),(3)

where the dividend yield qis constant.

Remark 1.1 Most dividend payments on an index — such as the Dow Jones Industrial Av-

erage (DJIA) or the Standard and Poor’s 500 (S&P500) — are so frequent that they can

be modeled as a continuous payment, which is the case in (3). However, if companies

only make two or four dividend payments per year, then one has to treat the dividend pay-

ments discretely and the question of how to incorporate discrete dividend payments into the

Black–Scholes equation arises.

Even though in this work we will focus on the case of continuous dividend payments,

we briefly review the results for discrete dividend payments from [63] in the sequel.

We assume that there is only one dividend payment of the dividend yield qduring the

lifetime of the option at the dividend date tq. Neglecting other factors, such as taxes, the

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 3

asset price Smust decrease exactly by the amount of the dividend payment qat time tq.

Thus we have the jump condition

S(t+

q) = (1 −q)S(t−

q),

where t−

q,t+

qdenote the moments just before and after the dividend date tq. This leads to

the following effect on the option price:

V(S, t−

q) = V((1 −q)S, t+

q),(4)

i.e. the value of the option at Sand time t−

qis the same as the value immediately after the

dividend date tqbut at the asset value (1 −q)S. In order to calculate the value of a Call

option with one dividend payment we solve the Black–Scholes equation from expiry t=T

until t=t+

qand use the relation (4) to compute the values at t=t−

q. Finally, we continue

to solve the Black–Scholes equation backwards starting at t=t−

qusing these values as the

initial data. The boundary conditions, that are discussed in the next section, do not need to

be modified for this case.

In the mathematical sense equations (2) and (3) are called convection–diffusion equa-

tions. The second-order term 1

2eσ2(t, S, VS, VSS)S2VSS is responsible for the diffusion, the

first-order term rSVSor (r−q)SVSis the convection term and −rV can be interpreted as

the reaction term (cf. [53,62]).

In the financial sense, the partial derivatives indicate the sensitivity of the option price

Vto the corresponding parameter and are called Greeks. The option delta is denoted by

∆ = VS, the option gamma by Γ = VSS and the option theta by θ=Vt[33].

Since American options can be exercised at any time before expiry, we need to find

the optimal time tof exercise, known as the optimal exercise time. At this time, which

mathematically is a stopping time, the asset price reaches the optimal exercise price or

optimal exercise boundary Sf(t).

This leads to the formulation of the problem for American options by dividing the do-

main [0,∞[×[0, T]of (3) into two parts along the curve Sf(t)and analyzing each of them

(see Fig. 1(a)). Since Sf(t)is not known in advance but has to be determined in the process

of the solution, the problem is called free boundary value problem (see e.g. [67]).

For different numerical approaches, the free boundary problem for American options

can be reformulated into a linear complementary problem, a variational inequality and a

minimization problem [26]. Here, we will only consider the formulation as a free boundary

problem.

For the American Call option the spatial domain is divided into two regions by the free

boundary Sf(t), the stopping region Sf(t)< S < ∞,0≤t≤T, where the option is

exercised or dead with V(S, t) = S−Kand the continuation region 0≤S≤Sf(t),

0≤t≤T, where the option is held or stays alive and equation (3) is valid under the

4 Julia Ankudinova and Matthias Ehrhardt

hold exercise

Sf(0)Sf(T)

Sf(t)

(a) Exercising and holding regions.

V(S, t)

Sf(0)

Sf(T)

Sf(t)

(b) Schematical value V(S, t).

Figure 1: American Call.

following terminal and boundary conditions:

V(S, T) = (S−K)+for 0≤S≤Sf(T)

V(0, t) = 0 for 0≤t≤T

V(Sf(t), t) = Sf(t)−Kfor 0≤t≤T(5)

VS(Sf(t), t) = 1 for 0≤t≤T

Sf(T) = max(K, rK/q).

For the sake of simplicity we will assume r > q in this chapter, and therefore we have

Sf(T) = rK/q for the American Call.

The structure of the value of an American Call can be seen Fig. 1(b), where we notice

that the free boundary Sf(t)determines the position of the exercise.

For American options, in general, analytic valuation formulae are not available, except

for a few special types, which we are not going to address in this chapter. Those types are

Calls on an asset that pays discrete dividends and perpetual Calls — meaning Calls with

an infinite time to expiry [40]. For the other types, there are various kinds of analytical and

numerical approximations that will be discussed in this chapter.

The structure of this chapter is as follows. In Section 2 several nonconstant volatility

models that lead to the nonlinearity of the Black–Scholes equation will be introduced. The

focus of this chapter is the solution of the resulting nonlinear problems for American Call

options. Since in general, a closed–form solution to the nonlinear Black–Scholes equation

for American options does not exist (even in the linear case), we have to solve the prob-

lems numerically. The numerical solution and the comparison study for American options

will be achieved by the transformation of the free boundary problem (3) subject to (5) into a

forward-in-time parabolic equation defined on a fixed (but unbounded) spatial domain (Sec-

tion 3). This new problem will be numerically solved by the method of finite differences

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 5

using an operator splitting technique (Section 4). It will then be evaluated and concisely

discussed in Section 5 thereafter.

2 Volatility Models

The essential parameter of the standard Black–Scholes model, that isnot directly observable

and is assumed to be constant, is the volatility σ. There have been many approaches to

improve the model by treating the volatility in different ways and using a modified volatility

function eσ(·)to model the effects of transaction costs, illiquid markets and large traders,

which is the reason for the nonlinearity of (2) and (3). In this section we will first give

a brief overview of several volatility models and then focus on the volatility models of

transaction costs.

•The constant volatility σin the standard Black–Scholes model can be replaced by

the estimated volatility from the former values of the underlying. This volatility is

known as the historical volatility [26].

•If the price of the option and the other parameters are known, which is e.g. the case

for the European Call and Put options, then the implied volatility can be calculated

from those Black–Scholes formulae. The implied volatility is the value σ, for which

the Black–Scholes equation. is true compared to the real market data. It can be

calculated implicitly via the difference between the observed option price V(from

the market data) and the Black–Scholes formulae, where all the parameters - except

for the implied volatility σ- are taken from the market data (the stock price S, the

time t, the expiration date T, the strike price K, the interest rate rthe dividend rate

q).

Considering options with different strike prices Kbut otherwise identical parameters,

we see that the implicit volatility changes depending on the strike price. If the implicit

volatility for a certain strike price Kis less than the implicit volatility for both the

strike price greater and less than K, this effect is called volatility smile (see e.g. [39]).

•Replacing the constant volatility with the observed implicit volatilities at each stock

price and time leads to the term of the local volatility eσ:= eσ(S, t). Dupire [16]

examines the dependencies and expresses the local volatility as a function of implicit

volatilities.

•Hull and White [32] and Heston [28] develop a model, in which the volatility follows

the dynamics of a stochastic process. This is known as the stochastic volatility.

•The assumption, that each security is available at any time and any size, or that in-

dividual trading will not influence the price, is not always true. Therefore, illiquid

markets and large trader effects have been modeled by several authors. In [23] Frey

6 Julia Ankudinova and Matthias Ehrhardt

and Stremme and later Frey and Patie [24] consider these effects on the price and

come up with the result

eσ=σ

1−ρλ(S)SVSS

,(6)

where σis the historical volatility, ρis constant, λ(S)is a strictly convex function

and λ(S)≥1. The function λ(S)depends on the pay-off function of the financial

derivative. For the European Call option, Frey and Patie show that λ(S)is a smooth,

slightly increasing function for S≥K. Bordag and Chmakova [7] assume that λ(S)

is constant and solve the problem (2) with the modified volatility (6) explicitly using

Lie-group theory (see also [12]).

The main scope of this chapter is the numerical solution of the nonlinear Black–Scholes

equation for the American Call option, where the nonlinearity results from transaction costs.

Therefore, after this general overview, we devote our attention to a more detailed description

of several transaction cost models.

2.1 Transaction Costs

The Black–Scholes model requires a continuous portfolio adjustment in order to hedge the

position without any risk. In the presence of transaction costs it is likely that this adjustment

easily becomes expensive, since an infinite number of transactions is needed [40]. Thus, the

hedger needs to find the balance between the transaction costs that are required to rebalance

the portfolio and the implied costs of hedging errors. As a result to this "imperfect" hedging,

the option might be over- or underpriced up to the extent where the riskless profit obtained

by the arbitrageur is offset by the transaction costs, so that there is no single equilibrium

price but a range of feasible prices. It has been shown that in a market with transaction costs

there is no replicating portfolio for the European Call option and the portfolio is required to

dominate rather than replicate the value of the option (see [4]). Soner, Shreve and Cvitaniˇ

prove in [54] that the minimal hedging portfolio that dominates a European Call is the trivial

one (hence holding one share of the stock that the Call is written on), so that efforts have

been made to find an alternate relaxation of the hedging conditions to better replicate the

pay-offs of derivative securities.

2.1.1 The model of Leland

Leland’s idea of relaxing the hedging conditions is to trade at discrete times [41], which

promises to reduce the expenses of the portfolio adjustment. He assumes that the transac-

tion cost κ|∆|S/2, where κdenotes the round trip transaction cost per unit dollar of the

transaction and ∆the number of assets bought (∆>0) or sold (∆<0) at price S, is

proportional to the monetary value of the assets bought or sold. Leland derives the relation

rBδt −κ

2|δ∆|S= (Vt+σ2

2S2VSS)δt, (7)

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 7

where Bis the bond and rthe riskless interest rate, and shows that

2|δ∆|S=σ2

2Le S2|VSS|δt. (8)

Here, Le denotes the Leland number, which is given by

Le = r2

πκ

σ√δt,(9)

with δt being the transaction frequency (interval between successive revisions of the port-

folio) and κthe round trip transaction cost per unit dollar of the transaction. Plugging (8)

and B= Π −∆S=V−SVSinto the equation (7) becomes

rV −rSVS−σ2

2LeS2|VSS|=Vt+σ2

2S2VSS.(10)

Therefore, Leland deduces that the option price is the solution of the nonlinear Black–

Scholes equation

0 = Vt+1

2eσ2S2VSS +rSVS−rV,

with the modified volatility

eσ2=σ21 + Le sign(VSS),(11)

where σrepresents the historical volatility and Le the Leland number. It follows from the

definition of the Leland number (9) that the more frequent the rebalancing (δt smaller), the

higher the transaction cost and the greater the value of V.

Leland’s model has played a significant role in financial mathematics, even though it has

been partly criticized by e.g. Kabanov and Safarian in [37], who prove that Leland’s result

has a hedging error. The restriction of his model is the convexity of the resulting option

price V(hence VSS >0) and the possibility to only consider one option in the portfolio.

Hoggard, Whalley and Wilmott study equation (2) with the modified volatility (11) for

several underlyings in [30]. An extension to this approach to general pay–offs is obtained

by Avellaneda and Parás in [3].

2.1.2 The model of Barles and Soner

In [4] Barles and Soner derive a more complicated model by following the utility function

approach of Hodges and Neuberger [29].

Supposing that the proportional transaction cost κis equal to a√εfor some constant

a > 0, they prove that as εand κgo to 0,Vis the unique (viscosity) solution of the

nonlinear Black–Scholes equation

0 = Vt+1

2eσ2S2V2

SS +rSVS−rV,

8 Julia Ankudinova and Matthias Ehrhardt

where

eσ2=σ21 + Ψ(er(T−t)a2S2VSS).(12)

Here σdenotes the historical volatility, a=κ/√εand Ψ(x)is the solution to the following

nonlinear ordinary differential equation (ODE)

Ψ′(x) = Ψ(x) + 1

2pxΨ(x)−x, x 6= 0,(13a)

with the initial condition

Ψ(0) = 0.(13b)

The analysis of this ODE (13) by Barles and Soner in [4] implies that

lim

x→∞

Ψ(x)

x= 1 and lim

x→−∞ Ψ(x) = −1.(14)

The property (14) encourages to treat the function Ψ(·)as the identity for large arguments

and therefore to simplify the calculations. In this case the volatility becomes

eσ2=σ2(1 + er(T−t)a2S2VSS).(15)

The existence of a viscosity solution to (2) for European options with the volatility given

by (12) is proved by Barles and Soner in [4] and their numerical results indicate an eco-

nomically significant price difference between the standard Black–Scholes model and the

nonlinear model with transaction costs.

2.1.3 Risk Adjusted Pricing Methodology

In this model, proposed by Kratka in [39] and improved by Jandaˇ

cka and Ševˇ

coviˇ

c in [35],

the optimal time-lag δt between the transactions is found to minimize the sum of the rate

of the transaction costs and the rate of the risk from an unprotected portfolio. That way the

portfolio is still well protected with the Risk Adjusted Pricing Methodology (RAPM) and

the modified volatility is now of the form

eσ2=σ21 + 3C2M

2πSVSS1

3,(16)

where M≥0is the transaction cost measure and C≥0the risk premium measure.

It is worth mentioning that these nonlinear transaction cost models that are described

above are all consistent with the linear model if the additional parameters for transaction

costs are equal to zero and vanish (Le,Ψ(·),M). We will study these models – more

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 9

precisely equations (2) and (3) where the volatility is given by the equations (11), (12), (15)

and (16) – for American Call options.

In general, an exact analytical solution leading to a closed expression is not known for

American options in a market with transaction costs. In the next section we will analytically

approach the solution of (3) by a transformation, that facilitates the numerical solution in

the section thereafter. We will compare and evaluate the results in the section thereafter.

3 The fixed Domain Transformation

The equation (3) subject to (5) is a backward-in-time free boundary problem. In order to

ease the numerical solution of (3) (5) for American Call options, we transform the prob-

lem into a problem posed on a a fixed (unbounded) domain additionally to the forward

transformation in time. Hence, the domain does not depend on the free boundary Sf(t)

anymore and we simply calculate an algebraic constraint equation for the position of the

free boundary. Following the idea of Ševˇ

coviˇ

c [52] we change the variables to:

τ=T−t, x = ln (τ)

S⇔S=e−x(τ), (τ) = Sf(T−τ),

such that x∈R+and τ∈[0, T].

Then, we construct a portfolio

Π(x, τ) = V(S, t)−SVS(S, t)

by buying ∆ = VSshares Sand selling an option V. Differentiating Πwith respect to x

and τgives us

Πx=VSSx−SxVS−SVSSSx=S2VSS

and

Πτ=VSSτ+Vttτ−SτVS−S(VSSSτ+VSttτ)

=−Vt−′(τ)

(τ)S2VSS +SVSt

=−Vt−′(τ)

(τ)Πx−S∂S(−Vt).

(17)

Substituting

−Vt=eσ2

2S2VSS −r(V−SVS)−qSVS=eσ2

2Πx−rΠ−qSVS

from (3) into (17) and using the fact that −S∂S=∂x, we get

Πτ=eσ2

2Πx−rΠ−qSVS−′(τ)

(τ)Πx+∂xeσ2

2Πx−rΠ+S∂S(qSVS)

=eσ2

2−′(τ)

(τ)−r+qΠx−rΠ + 1

2∂x(eσ2Πx).

10 Julia Ankudinova and Matthias Ehrhardt

We therefore obtain

0 = Πτ+b(τ)−eσ2

2Πx−1

2∂x(eσ2Πx) + rΠ,(18)

defined on x∈R+,0≤τ≤T, where the coefficient b(τ)is

b(τ) = ′

(τ) + r−q.

The terminal condition from (5) in the original variables (S, T )becomes the initial condi-

tion in the new variables (x, 0):

Π(x, 0) = V(S, T)−SVS(S, T)

=(−Kfor S > K ⇔x < ln (0)

0otherwise .(19a)

and the boundary conditions from (5) transform to

Π(x, τ) = 0 as x→ ∞,0≤τ≤T, (19b)

Π(0, τ) = −Kfor 0≤τ≤T. (19c)

To complete the system of equations that enables the computation of the portfolio Πwe

need to use the last two conditions of (5) to obtain an expression at the free boundary

position (τ). Differentiating and evaluating V(Sf(t), t) = Sf(t)−Kat the free boundary

gives us

VS(Sf(t), t)S′

f(t) + Vt(Sf(t), t) = S′

f(t).

Using (5), we conclude that

Vt(Sf(t), t) = 0 for 0≤τ≤T.

Computing (3) at the point (Sf(t), t)or at (0, τ)in the transformed variables yields:

0 = Vt(Sf(t), t) + 1

2eσ2Πx(0, τ) + (r−q)Sf(t)VS(Sf(t), t)−rV (Sf(t), t)

2eσ2Πx(0, τ) + rK −q(τ).

We remind the reader that we have assumed r≥qand therefore we obtain the last condi-

tion:

(τ) = 1

2qeσ2Πx(0, τ) + rK

qwith (0) = rK

q,(19d)

where 0≤τ≤Tand eσ2depends on the volatility model we choose. The volatility (11)

from Leland’s model becomes

eσ2=σ21 + Le sign(Πx),(20a)

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 11

for Barles’ and Soner’s model (12) we get

eσ2=σ2(1 + Ψ(erτ a2Πx)),(20b)

for the identity model (15) we obtain

eσ2=σ2(1 + erτ a2Πx),(20c)

and for the Risk Adjusted Princing Methodology (16) we derive

eσ2=σ21 + 3C2M

2πΠx(τ)e−x1

3.(20d)

This transformed problem (18) subject to (19) with the corresponding volatilities (20) is

solved by the split-step finite-difference method proposed by Ševˇ

coviˇ

c [52]

Once we have numerically solved the transformed problem by calculating the solution

to our portfolio Π(x, τ)and the free boundary (τ), we calculate the value of the American

Call V(S, t)option by transforming

Π(x, τ) = V(S, t)−SVS(S, t)

back to the original variables. Since we know that

Π(x, τ)

S2=V(S, t)

S2−VS(S, t)

S=∂S−V(S, t)

S,

we integrate the above equation from Sto Sf(t), take into account the boundary condition

V(Sf(t), t) = Sf(t)−Kand obtain:

ZSf(t)

Π(ln((τ)/S), τ)

S2dS =ZSf(t)

∂S−V(S, t)

SdS

Zln (τ)

Sf(t)

ln (τ)

Π(x, τ)

S2(−S)dx =−V(Sf(t), t)

Sf(t)+V(S, t)

SZln (τ)

Π(x, τ)

e−x(τ)dx =−S(τ)−K

(τ)+V(S, t)

V(S, T −τ) = S

(τ)(τ)−K+Zln (τ)

exΠ(x, τ)dx.(21)

Therefore, (21) yields the price of the American Call option V(S, t)in the presence (and

absence) of transaction costs.

12 Julia Ankudinova and Matthias Ehrhardt

4 Numerical Solution

Due to the lack of general closed–form solutions to the Black–Scholes equations, there are

various numerical methods for solving Black–Scholes equations for American options.

For European Call and Put options, the Black–Scholes formulae provide the correct

answer, but for more complicated contracts in more general settings analytical formulae

are seldom available and numerical methods have to be used to solve the problem. These

vary from lattice methods (including binomial and trinomial approximations [14]), Monte-

Carlo methods using the least-square techniques [34], analytical approximations [5,11,42],

finite-element discretizations [26] to finite-difference methods [2,9,13].

There are numerous other methods for pricing American options including the method

of lines [45], front-tracking algorithms [64], penalty methods [68] and many others. One

of the standard approaches for solving the Black–Scholes equation for American options

consists of the transformation of the original equation into the heat equation posed on a

semi–unbounded domain with a free boundary Sf(t)[53,63]. For a new alternative direct

method using the Mellin transformation we refer to [36,47].

Up to now, an exact analytical formula for the free boundary profile Sf(t)in (3) subject

to (5) is not known, but several authors derived approximate expressions to evaluate Amer-

ican Call and Put options in the linear case [25]. Recently, in a promising approach [51],

Ševˇ

coviˇ

c obtained a semi–explicit formula for an American Call in the case of r > q.

By transforming the linear Black–Scholes equation for the American Call option into a

nonlinear parabolic equation on a fixed domain and applying Fourier sine and cosine trans-

formations, he derives a nonlinear singular integral equation determining the shape of the

free boundary. This integral equation can be solved effectively by the means of successive

iterations.

Another standard method consists of the reformulation of the free boundary problem

into a linear complementary problem (LCP) and the solution by the Projected Successive

Over Relaxation (PSOR) method of Cryer [15]. Alternatively, penalty and front–fixing

methods are developed (e.g. in [22,46]). A disadvantage of these methods is the change of

the underlying model.

A different approach [31] is based on a recursive calculation of the early exercise bound-

ary, estimating the boundary only at some points and then approximating the whole bound-

ary by Richardson extrapolation. Explicit boundary tracking algorithms are e.g. a finite-

difference bisection scheme [38] or the front–tracking strategy of Han and Wu [27].

This emphasis of this chapter is on finite-difference schemes, thus other methods will

not be further elaborated on here. For more information on numerical methods we refer the

reader to [48,49,66] and the references therein.

4.1 American Call option

Now we want to solve the transformed problem from the previous section.

0 = Πτ+b(τ)−eσ2

2Πx−1

2∂x(eσ2Πx) + rΠ, x ∈R+,0≤τ≤T(22)

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 13

with the corresponding volatilities (20) subject to the conditions

Π(x, 0) = −Kfor x < ln (0)

0otherwise

Π(x, τ) = 0 as x→ ∞,0≤τ≤T,

Π(0, τ) = −Kfor 0≤τ≤T,

(23)

and the constraint

(τ) = 1

2qeσ2Πx(0, τ) + rK

qwith (0) = rK

q.(24)

We therefore first describe the solution of (22) subject to (23) and (24) with the correspond-

ing volatilities (20) by finite-difference schemes and then present the numerical results.

4.2 Finite-Difference Schemes

There have been many approaches to calculate the value of an American option numeri-

cally by compact finite-difference schemes in the absence of transaction costs. Recently,

Tangman et al. [59, 60] introduced a compact scheme of order (2,4). Two other com-

pact schemes, known as the Numerov-type (see [58,66]) and the Crandall-Douglas scheme

(see [43]), are analyzed for linear Black–Scholes equations. However, these schemes are

not directly transferable to the model in the presence transaction costs.

In order to find a solution for the nonlinear Black–Scholes equation (22) subject to (23)

with the corresponding volatilities (20) and the constraint (24), Ševˇ

coviˇ

c suggests to com-

bine two approaches that solve the problem for the American Call with a constant volatility

numerically [52]. One of them is the transformation of the problem into a variational in-

equality and its solution by the PSOR algorithm [26,53]. The other one is the derivation of

a nonlinear integral equation for the position of the free boundary without the knowledge

of the price itself [40,64].

Even though these methods are not directly applicable, since they require a constant

volatility σ, this approach is successful when it is combined with an operator splitting tech-

nique. The idea is to discretize (22) in time, to split the equation into a convective and a

diffusive part and to find an approximation for the solution pair (Π, )at each time level.

The detailed derivation is given in the sequel.

4.2.1 Grid

We discretize the problem (22) subject to the conditions (23) with the corresponding volatil-

ities (20) by confining the unbounded domain x∈R+and τ∈[0, T]to x∈(0, R)with

R > 0sufficiently large (see [52]). For the calculation Ševˇ

coviˇ

c chooses to take R= 3,

since this is equivalent to S∈(Sf(t)e−R, Sf(t)) and yields a good approximation for

S∈(0, Sf(t)) (as the transformation was S=Sf(t)e−x). In the sequel we refer to h > 0

14 Julia Ankudinova and Matthias Ehrhardt

time

space

xixi+1

τn

τn+1

R=Nh

T=Mk

Figure 2: Uniform grid for an American Call option.

as the spatial step and to k > 0as the time step, xi=ih,i∈[0, N],R=Nh and τn=nk,

n∈[0, M],T=Mk (see Fig. 2).

The approximate solution of (22) in xiat time τnis denoted by Πn

i:= Π(xi, τn), the

value of the free boundary at time τnby n:= (τn)and the value of the coefficient b(τ)

at τnby bn:= b(τn).

We treat the initial and boundary conditions (23) in the following way:

Π0

i= Π(xi,0) = (−Kfor xi<ln (0)

K= ln r

0otherwise ,

Πn

0=−K,

Πn

N= 0.

(25)

4.2.2 Difference Quotients

We denote the forward difference quotient with respect to the spatial variable in xiat time

τnwith the spatial step size hby:

D+

hΠn

i:= Πn

i+1 −Πn

h≈Πx(xi, τn),

the backward difference quotient by:

D−

hΠn

i:= Πn

i−Πn

i−1

h≈Πx(xi, τn)

and the central difference quotient by

hΠn

i:= Πn

i+1 −Πn

i−1

2h≈Πx(xi, τn),

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 15

omitting the truncation error O(h),O(h)and O(h2), respectively. For the second spatial

derivative we introduce the standard difference quotient

hΠn

i:= Πn

i+1 −2Πn

i+ Πn

i−1

h2≈Πxx(xi, τn),

with the error term O(h2).

4.3 Volatility Functions

The volatilities (20) can all be written in the form

eσn

i2=σ2(1 + sn

i),

where sn

idenotes the volatility correction in xiat time τn. We choose forward differences

to approximate Πxin the volatility formulae, so that for Leland’s model with the volatility

(20a) our volatility correction becomes

i=Le sign(D+

hΠn

i),(26a)

for the volatility correction in Barles’ and Soner’s model with the volatility (20b) we get

i= Ψerτna2D+

hΠn

i,(26b)

for the volatility correction in case of treating Ψ(·)as the identity with the original volatility

(20c) we obtain

i=erτna2D+

hΠn

i,(26c)

and for the volatility (20d) in the Risk Adjusted Pricing Methodology (RAPM) the volatility

correction is

i= 3C2M

2πD+

hΠn

ine−xi1

3.(26d)

4.4 The Treatment of the Free Boundary

We discretize the free boundary (24) by approximating the spatial derivative at the origin

x= 0 by forward differences and obtain:

n=1

2qσ2(1 + sn

0)D+

hΠn

0+rK

qwith 0=rK

q,(27)

where sn

0denotes (26) at x= 0 depending on the volatility model.

Note, that in case of the RAPM, where the volatility correction is given by equation

(26d), sn

0depends on nand therefore nin (27) is expressed by a fixed point equation.

16 Julia Ankudinova and Matthias Ehrhardt

Remark 4.1 For the American Call option (in contrast to the American Put option) it is

possible to derive a series for the location of the optimal exercise boundary close to expiry

using standard asymptotic analysis [1,63]. This local analysis of the free boundary Sf(t)

yields

Sf(t)∼Sf(T) 1 + ξ0r1

2σ2(T−t) + ...!,as t→T, (28)

where ξ0= 0.9034 ... is a universal constant of Call option pricing. Equation (28) can be

rewritten as

(τ)∼(0) 1 + ξ0r1

2σ2(τ) + ...!,as τ→0.(29)

With only very few terms we get a fairly accurate result for the free boundary and thus

equation (29) will serve us as a check for the case of a constant volatility eσ2=σ2(see

Fig. 3). Note that this result is especially useful in the first time levels of a numerical

calculation where rapid changes in (τ)influence the whole solution region.

(τ)

Figure 3: Asymptotic solution for the free boundary (τ)with T= 1,K= 10,σ= 0.2,

r= 0.1,q= 0.05.

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 17

4.5 The Splitting in Time Method

We approximate the time derivative of (22) by backward differences D−

kΠn

i, the first and

second spatial derivatives by central differences D0

hΠn

iand D2

hΠn

i. Then, (22) becomes:

0 = D−

kΠn

i+bn−σ2

2(1 + sn

i)D0

hΠn

i−1

2∂xσ2(1 + sn

i)D0

hΠn

i+rΠn

i(30)

subject to the Dirichlet conditions (25). We introduce an intermediate step at time τn−1

such that

D−

kΠn

i=Πn

i−Πn−1

k=Πn

i−Πn−1

i+ Πn−1

i−Πn−1

and then split the problem (30) into a convective part with the linear first-order term

bnD0

hΠn

0 = Πn−1

i−Πn−1

k+bnD0

hΠn

i(31)

and a diffusive part with the nonlinear first- and second-order terms σ2/2(1 + sn

i)D0

hΠn

and −∂xσ2/2(1 + sn

i)D0

hΠn

i:

0 = Πn

i−Πn−1

k−σ2

2(1 + sn

i)D0

hΠn

i−1

2∂xσ2(1 + sn

i)D0

hΠn

i+rΠn

i.(32)

Assuming that D0

hΠn

i≈D0

hΠn−1

i, which is reasonable for small time steps k, we can

approximate the convective part (31) as

0 = Πn−1

i−Πn−1

k+bnD0

hΠn−1

i.(33)

Now the solution to (32)-(33) gives a good approximation to the solution of (30) (see [52]).

This decomposition of the problem is called Lie-Splitting and is a spitting of order 1in time.

4.5.1 Convective part

First, we solve the convective part (33), which can be approximated by an explicit solution

to the transport equation

Πτ+b(τ)Πx= 0,(34)

for (x, τ)∈R×[0, T ], subject to the boundary and initial conditions

Π(0, τ) = −K,

Π(x, 0) = (−Kfor x < ln r

0otherwise = Π0(x).(35)

18 Julia Ankudinova and Matthias Ehrhardt

We then know by the theory of partial differential equations (see e.g. [20]) that the solution

for this problem (34)–(35) is

Π(x, τ) = Π(x−Zτ

b(s)ds, 0) = Π0(x−Zτ

b(s)ds)(36)

with the primitive function Rb(s)ds =B(τ) + c= ln (τ) + (r−q)τ+c. Hence,

considering the problem (34) for (xi, τj)∈R×[τn−1, τn]subject to the boundary and

initial conditions

Π(0, τj) = −K,

Π(xi, τn−1) = Πn−1(xi),(37)

we know that the solution is given by

Π(xi, τj) = Πxi−Zτj

τn−1

b(s)ds, τn−1

=(Π(ξj

i, τn−1)for ξj

i>0

−Kotherwise, .

(38)

where we set ξj

i=xi−B(τj) + B(τn−1) = xi−ln j

n−1−(τj−τn−1)(r−q). Then we

can write

Πn−1

i=(Π(ξn

i, τn−1)ξn

i=xi−ln n

n−1−k(r−q)>0

−Kotherwise. (39)

Here, we use a linear approximation between the discrete values Π(xi, τn−1),i∈Nin

order to compute the value of Π(ξn

i, τn−1).

Hence, (39) is the solution to the convective part (33) of the problem (30).

4.5.2 Diffusive part

We solve the diffusive part (32) of the problem (30) by the finite-difference method. We

approximate the second spatial derivative by central differences D2

hΠn

iand the first spatial

derivative by both central D0

hΠn

iand backward differences D−

hΠn

i. Then, (32) becomes:

0 = Πn

i−Πn−1

k−σ2

2(1 + sn

i)Πn

i+1 −Πn

i−1

2h+rΠn

−σ2

2(1 + sn

i)Πn

i+1 −2Πn

i+ Πn

i−1

h2+(1 + sn

i)−(1 + sn

i−1)

Πn

i−Πn

i−1

h

=Πn

i−Πn−1

k−σ2

2(1 + sn

i)Πn

i+1 −Πn

i−1

2h+rΠn

−σ2

2(1 + sn

i)Πn

i+1 −Πn

h2−(1 + sn

i−1)Πn

i−Πn

i−1

h2.

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 19

Rearranging leads to a tridiagonal system of equations

Πn−1

i=an

iΠn

i−1−bn

iΠn

i+cn

iΠn

i+1,(40)

with the coefficients

i=σ2

2(1 + sn

i)k

2h−σ2

2(1 + sn

i−1)k

h2,

i= 1 + kr +σ2

2(1 + sn

i)k

h2+σ2

2(1 + sn

i−1)k

h2,

i=−σ2

2(1 + sn

i)k

2h−σ2

2(1 + sn

i)k

h2.

Equation (40) can be written in the form of matrices:

Πn−1

2=AnΠn+dn,(41)

where

Πn=Πn

1,···,Πn

N−1⊤∈RN−1,

An=







1cn

10··· 0

2bn

.......

0.........0

.......bn

N−2cn

N−2

0··· 0an

N−1bn

N−1







∈R(N−1)×(N−1),

and

dn=an

1Πn

0,0,···,0, cn

N−1Πn

N⊤∈RN−1.

Therefore, (41) solves the diffusive part (32) of the problem (30).

Now, we have a set of nonlinear equations (26), (27), (39) and (41) that delivers the

solution to our portfolio Π(x, τ)and to the free boundary (τ), from which we can calculate

the value of the American Call option V(S, t)with equation (21).

In order to see the dependencies of the equations, we rewrite them in the following

abstract form:

sn=D(Πn, n),

n=FΠn, sn=FΠn, n,

Πn−1

2=GΠn−1, n, n−1=GΠn−1, n,

AsnΠn=AΠn, nΠn= Πn−1

2−d(sn),

(42)

where

sn=sn

0,···, sn

N⊤∈RN+1,

20 Julia Ankudinova and Matthias Ehrhardt

D(·)denotes the right-hand side of (26), F(·)is the right-hand side of (27), G(·)is the

right-hand side of the transport equation (39), A(·)is the tridiagonal matrix and d(sn)the

vector as defined in (41).

As we can see by this notation (42), both nand Πnare given in terms of themselves,

hence each is given in terms of nand Πn. This problem can be approximately solved by a

successive fixed point iteration over p∈Nat each time level n.

Following Ševˇ

coviˇ

c [52] we define for n≥1:Πn,0= Πn−1,n,0=n−1and sn,0=

sn−1. Then the (p+1)-th approximation of Πn,nand snis obtained as the solution of the

system:

sn,p+1 =D(Πn,p, n,p),

n,p+1 =FΠn,p, sn,p+1,

Πn−1

2,p+1 =GΠn−1,p, n,p+1,

Asn,p+1Πn,p+1 = Πn−1

2,p+1 −d(sn,p+1).

(43)

Both the volatility correction sn,p+1

i, the free boundary n,p+1 and the solution Πn−1

2,p+1 to

the convective part (31) can be directly computed from (26), (27) and (39) respectively. The

solution Πn,p+1 to the diffusive part (32) has to be calculated from the system of equations

(41).

Assuming that the system (43) converges to some limiting values sn,pmax ,n,pmax ,

Πn−1

2,pmax and Πn,pmax at each time level n[52], we can calculate V(Si, tn) =

V(e−xin, T −τn)with these values and proceed to the next time level n+ 1.

From (21) we then know that:

V(Si, tn) = e−xin−K+Ii,(44)

where

Ii=

i−1

j=0 Ik+Zxi

xi−1

exΠ(x, τ)dx

i−1

j=0 Ik+xi−xi−1

2exi−1Πn

i−1+exiΠn

i.

Here, we use the trapezoidal rule in order to approximate the integral in equation (21).

We summarize the calculation of the price V(S, t)for the American Call option in the

presence or absence of transaction costs by the Algorithm 1 given in the appendix.

5 Comparison Study

Based on the iterative algorithm described in the previous section (Algorithm 1), we solve

the transformed Black–Scholes equation (22) with the corresponding volatilities (20) for

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 21

the American Call option and finally transform Π(x, τ)back to the original option price

V(S, t).

The main purpose of this section is to compare the resulting option value V(S, t)and

the free boundary Sf(T−t) = (τ)for the four different transaction cost models (26) to

the linear model (σconstant) and to each other.

We choose pmax = 5 for the successive iteration over pin our algorithm in order to

solve the system (42) with the precision of 10−7[52]. We use the following parameters to

calculate Π(x, τ)and (τ):

r= 0.1, σ = 0.2, K = 10, T = 1 (one year), R = 3.

We start by comparing the free boundary (τ)computed with Algorithm 1 to the asymptotic

solution (29) from Remark 4.1 for the linear case (sn

i= 0). In Fig. 4 we observe that for

smaller spatial steps h→0the freeboundary computed bythe iterative algorithm converges

monotonically towards the asymptotic solution (29) from below.

0 1

20.5

21.5

22.5

asymptotic solution

h=0.0086

h=0.01

h=0.012

h=0.015

h=0.03

(τ)

Figure 4: Free boundary for various spatial steps hwith a constant time step k= 0.0008

and a constant volatility σ2computed by Algorithm 1 vs. the asymptotic solution of (29).

We keep the time step k= 0.0008 constant and see that for h= 0.0086 (purple line)

the free boundary at Tis computed by our algorithm as (T)≈22.2201. The asymptotic

solution at Tis (T)≈22.5552, which means a relative error of 1.49%. The free boundary

values for the other spatial steps can be seen in Table 1.

22 Julia Ankudinova and Matthias Ehrhardt

h0.03 0.015 0.012 0.01 0.0086

(T) 21.8764 22.1111 22.1619 22.1955 22.2201

Table 1: Values of the free boundary position for various spatial steps hwith a constant

time step k= 0.0008 and a constant volatility σ2.

Since the asymptotic solution of (29) is only an approximation, we are satisfied by our

results and take the free boundary (T)≈22.1111 for k= 0.0008,h= 0.015 (blue line

in Fig. 4) as our reference solution in the absence of transaction costs for the sake of the

computational time.

Fig. 5 shows the structure of the price for the American Call option V(S, t)without

transaction costs with k= 0.0008 and h= 0.015. It is computed with the iterative algo-

rithm described in the previous sections and the parameters above.

0.5

V(S, t)

Figure 5: Value of an American Call option V(S, t)in the absence transaction costs com-

puted with Algorithm 1 determined by the free boundary (red line).

The corresponding synthetic portfolio Π(x, τ)in the absence of transaction costs is

illustrated in Fig. 6. Note, that we include rounding and discretization errors when trans-

forming Π(x, τ)back into V(S, t), since equation (44) involves an integral approximation.

However, the analysis of V(S, t)is more interesting for us and we therefore assume that

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 23

these errors are sufficiently small due to the chosen mesh.

0.5

−10

Π(x, τ)

xτ

(a) 3-D profile.

0 1 2 3

−10

Π(x, τ)

τ= 0

τ= 0.2

τ= 0.4

τ= 0.6

τ= 0.8

τ= 1

(b) Profile at different time points.

Figure 6: Value of the synthetic portfolio Π(x, τ)in the absence of transaction costs com-

puted with Algorithm 1.

We now compare the price V(S, 0) computed by Algorithm 1 to the price VP SOR(S, 0)

computed by the PSOR algorithm in the linear case sn

i= 0. Fig. 7 shows that with the

given mesh size k= 0.0008 and h= 0.015 the price computed by our algorithm (Fig. 7(a))

only slightly differs from the price computed by the PSOR algorithm (Fig. 7(b)).

0 10 20 25

V(S, 0)

(a) Computed with Algorithm 1.

0 10 20 25

V(S, 0)

(b) Computed with PSOR.

Figure 7: Price of an American Call option V(S, 0) in the absence of transaction costs and

the pay-off V(S, T )(red dotted line).

We calculate the error of accuracy of our computation one year to expiry at t= 0,

denoted by the ℓ2-error

err2(0) = h

i=0 |VPSOR(Si,0) −V0

i|21

24 Julia Ankudinova and Matthias Ehrhardt

where VPSOR(Si,0) denotes the solution computed by the PSOR algorithm at Si=

e−ih(T)and (T)depends on the step size h. For this purpose, we interpolate the so-

lution computed by the PSOR algorithm by the MATLAB routines spline and ppval.

For V0

iwe use our corresponding solution, where k= 0.0008. The error can be seen in

Table 2, which reveals that it is reasonable to assume the accuracy O(h).

h0.03 0.015 0.012 0.01 0.0086

ℓ2-error 0.0365 0.0162 0.0257 0.0084 0.0167

Table 2: ℓ2-error of accuracy of Algorithm 1 compared to the PSOR algorithm in the ab-

sence of transaction costs.

We further compute the free boundary profiles for the four different transaction cost

models (26) by Algorithm 1 and compare them to the profile of the free boundary in the

absence transaction costs. For our computations we take k= 0.0008 and h= 0.015.

0 1

20.5

21.5

22.5

23.5

no transaction costs

Identity

Barles and Soner

Leland

RAPM

(τ)

Figure 8: Free boundary positions for various transaction cost models vs. the free boundary

profile in the absence of transaction costs.

As expected, we see that for all the transaction cost models the free boundary values

are greater than in the case without transaction costs (Fig. 8). With the given parameters

the free boundary in the absence of transaction costs is (T)≈22.11, followed by the

identity model with a= 0.02 ((T)≈22.16), Barles’ and Soner’s model with a= 0.02

((T)≈22.34), Leland’s model with δt = 0.1,κ= 0.02 ((T)≈22.44) and finally the

RAPM with C= 0.01,R= 30 ((T)≈23.39).

Furthermore, we compute the corresponding values V(S, t)for the American Call op-

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 25

tion by Algorithm 1 and check the price difference between the American Call option with

transaction costs and the American Call option without transaction costs

Vnonlinear(S, t)−Vlinear(S, t).

The influence of transaction costs for the four models can be seen in Fig. 9. We notice that

the difference is maximal one year to expiry at t= 0 and S≈9.5. The difference is not

symmetric, but decreases towards the expiry. This seems plausible, since towards expiry the

portfolio can not be adjusted as often at it could be adjusted before. Hence, the transaction

costs and the value of the American Call option with transaction costs decrease towards

t= 1.

0.5

010 20 25

0.15

V(S, t)

(a) Barles’ and Soner’s model (a= 0.02) vs. linear

model.

0.5

010 20 25

5x 10−3

V(S, t)

(b) Ψ(x) := xchosen as the identity (a= 0.02) vs.

linear model.

0.5

010 20 25

0.1

V(S, t)

model.

0.5

010 20 25

0.2

0.4

V(S, t)

(d) RAPM model (M= 0.01,C= 30) vs. linear

model.

Figure 9: The influence of transaction costs Vnonlinear(S, t)−Vlinear(S, t).

The corresponding prices V(S, 0) in the presence of transaction costs can be seen in

Fig. 10. At S≈9.5with the parameters as indicated above and k= 0.0008,h= 0.015

the price of the American Call option evaluated with the RAPM transaction cost model is

the highest (≈1.06). It is followed by Barles’ and Soner’s model (≈0.82), Leland’s model

26 Julia Ankudinova and Matthias Ehrhardt

(≈0.78), the identity model (≈0.74) and finally the model in the absence of transaction

costs (≈0.71). As already shown in Table 2, the linear price computed by our algorithm

(light blue solid line in Fig. 10) only slightly deviates from the price computed by the PSOR

algorithm (black dotted line in Fig. 10).

9 9.5 10 10.5 11

0.5

1.5

Pay−off V(S,T)

no transaction costs Algorithm 2

no transaction costs PSOR

Identity

Leland

Barles and Soner

RAPM

V(S, 0)

Figure 10: Price of an American Call option V(S, 0) for different transaction cost models

vs. the price without transaction costs.

For other numerical experiments in the future is recommendable to use rather Cor C++

in order to reduce the computational time which is relatively high in MATLAB.

Conclusion

In this chapter we solved the nonlinear Black–Scholes equation for American options in the

presence of transaction costs. Summing up, our numerical results showed a considerable

price difference between linear and nonlinear prices for American Call options.

While we focused in this chapter on standard options (known as plain–vanilla options)

of American type, our future work will deal with extensions: forward and future contracts,

options on futures, more general pay–off functions (e.g. ‘cash–or–nothing call’) with trans-

action costs and instalment options.

Moreover, we will consider a higher-order splitting in time, e.g. the well–known

Strang–Splitting [56] and combine this with modern compact finite difference of high spa-

tial order, like the Crandall–Douglas Scheme [43] which is fourth-order accurate in ‘space’

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 27

(i.e. asset price) or the high–order compact methods proposed in [59], [60], [66]. Espe-

cially, the method of [60] is promising, since it is already an improvement of the Han and

Wu method [27] with a higher order interior scheme and more accurate tracking of the free

boundary.

Acknowledgements

The authors acknowledge very stimulating discussions with Daniel Ševˇ

coviˇ

c in the frame-

work of a bilateral German–Slovakian DAAD project Fin-Diff-Fin: Finite differences for

Financial derivative models.

28 Julia Ankudinova and Matthias Ehrhardt

Algorithm 1 Computation of the price V(S, t)for the American Call

Appendix

Require: R,T,h,k,M,N,r,K,D,σ,Le,a,C,M

1: solve the ODE (13) required for the volatility model of Barles and Soner and interpolate

the solution

2: initialize Π0

3: initialize the free boundary 0=rK/q

4: transform Π0into V0

5: set π= Π0and v=V0

6: set Π1,0= Π0and 1,0=0

7: calculate Πn,nat each time level

8: for n= 1 : Mdo

9: calculate sn,p,n,p,Πn−1/2,p and Πn,p in the successive loop over p

10: for p= 1 : pmax do

11: calculate the volatility correction sn,p depending on the volatility model using

Πn,p−1and n,p−1(in the case of Barles’ and Soner’s model use the interpo-

lated solution of (13), in the case without transaction costs sn,p = (0,···,0)⊤∈

RN+1)

12: calculate n,p using Πn,p−1and sn,p

13: calculate Πn−1/2,p using Πn−1and n,p

14: fill the matrix An,p and the vector dn,p with the corresponding coefficients using

sn,p

15: L-R-decompose An,p =Ln,pRn,p

16: solve Ln,pyn,p = Πn−1/2,p −dn,p for yn,p

17: solve Rn,pΠn,p =yn,p for Πn,p

18: start over with the loop over p

19: end for

20: set Πn= Πn,p and n=n,p

21: transform Πninto Vn

22: save the solution in the transformed variables in the array

π=π[−K; Πn; 0]

23: save the solution in the original variables in the array

v=v[n−K;Vn; 0]

24: start over with the loop over n

25: end for

26: plot vat each time level and each stock price, plot at each time level

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 29

References

[1] G. Alobaidi, American options and their strategies, Ph.D. thesis, University of Western

Ontario, London, Canada, 2000.

[2] J. Ankudinova and M. Ehrhardt, The numerical solution of nonlinear Black–Scholes

equations, in press: Comput. Math. Appl., 2008.

[3] M. Avellaneda and A. Parás, Dynamic hedging portfolios for derivative securities in

the presence of large transaction costs, Appl. Math. Finance 1(1994), 165–193.

[4] G. Barles and H.M. Soner, Option pricing with transaction costs and a nonlinear

Black–Scholes equation, Finance Stoch. textbf2 (1998), 369–397.

[5] G. Barone-Adesi and R. Whaley, Efficient Analytic Approximation of American Option

Values, J. Finance 42 (1987), 302–320.

[6] F. Black and M. Scholes, The pricing of options and corporate liabilities, J. Polit.

Econ. 81 (1973), 637–659.

[7] L. Bordag and A. Chmakova, Explicit solutions for a nonlinear model of financial

derivatives, Int. J. Theor. Appl. Finance 10 (2007), 1–21.

[8] P. Boyle and T. Vorst, Option Replication in Discrete Time with Transaction Costs,

J. Finance 47 (1992), 271–293.

[9] M. Brennan and E. Schwartz, The valuation of American put options, J. Finance 32

(1977), 449–462.

[10] M. Broadie and J. Detemple, American Option valuation: New bounds, approxima-

tions, and a comparison of existing methods, Rev. Fin. Stud. 9(1996), 1211–1250.

[11] P. Carr, R. Jarrow and R. Myneni, Alternative Characterizations of American Put

Options, Math. Finance 2(1992), 87–105.

[12] A. Chmakova, Symmetriereduktionen und explizite Lösungen für ein nichtlineares

Modell eines Preisbildungsprozesses in illiquiden Märkten, Dissertation, Technische

Universität Cottbus, 2005 (in german).

[13] G. Courtadon, A More Accurate Finite Difference Approximation for the Valuation of

Options, J. Fin. Quant. Anal. 15 (1982), 697–703.

[14] J. Cox and S. Ross and M. Rubinstein, Option pricing: A simplified approach, J. Fin.

Econ. 7(1979), 229–263.

[15] E. Cryer, The solution of a quadratic programming problem using systematic overre-

laxation, SIAM J. Control 9(1971), 385–392.

30 Julia Ankudinova and Matthias Ehrhardt

[16] B. Dupire, Pricing with a Smile, Risk 7(1994), 18–20.

[17] M. Ehrhardt, Discrete Transparent Boundary Conditions for Parabolic Equations, Z.

Angew. Math. Mech. 77 (1997) S2, 543–544.

[18] M. Ehrhardt, Discrete Artificial Boundary Conditions, Ph.D. Thesis, Technische Uni-

versität Berlin, 2001.

[19] M. Ehrhardt and R.E. Mickens, A fast, stable and accurate numerical method for the

Black–Scholes equation of American options, accepted for publication: Int. J. Theoret.

Appl. Finance, 2008.

[20] L.C. Evans, Partial Differential Equations, American Mathematical Society, Provi-

dence, 1998.

[21] G. Ferreyra, The Mathematics Behind the 1997 Nobel Prize in Economics, e-MATH,

What’s New in Mathematics, (electr. magazine of the Amer. Math. Soc.), December

1997.

[22] P.A. Forsyth and K.R. Vetzal, Quadratic Convergence for Valuing American Options

using a Penalty Method, SIAM J. Sci. Comput. 23 (2002), 2095–2122.

[23] R. Frey and A. Stremme, Market Volatility and Feedback Effects from Dynamic Hedg-

ing, Math. Finance 4(1997), 351–374.

[24] R. Frey and P. Patie, Risk Management for Derivatives in Illiquid Markets: A Simula-

tion Study, Springer, Berlin, 2002.

[25] R. Geske and R. Roll, On valuing American call options with the Black–Scholes Eu-

ropean formula, J. Finance 89 (1984), 443–455.

[26] M. Günther and A. Jüngel, Finanzderivate mit MATLAB, Vieweg, 2003 (in german).

[27] H. Han and X. Wu, A fast numerical method for the Black–Scholes equation of Amer-

ican options, SIAM J. Numer. Anal. 41 (2003), 2081–2095.

[28] S. Heston, A Closed–Form Solution for Options with Stochastic Volatility with Appli-

cations to Bond and Currency Options, Rev. Fin. Studies 6(1993), 327–343.

[29] S. Hodges and A. Neuberger, Optimal replication of contingent claims under transac-

tion costs, Review of Futures Markets 8(1989), 222–239.

[30] T. Hoggard and A. E. Whalley and P. Wilmott, Hedging option portfolios in the pres-

ence of transaction costs, Adv. Fut. Opt. Res. 7(1994), 21–35.

[31] J.-Z. Huang, M.G. Subrahmanyam and G.G. Yu, Pricing and Hedging American Op-

tions: A Recursive Integration Method, Rev. Fin. Studies 9(1996), 277–300.

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 31

[32] J. Hull and A. White, The Pricing of Options on Assets with Stochastic Volatilities,

J. Finance 42 (1987), 281–300.

[33] J. Hull, Options, Futures and Other Derivatives, 5th Edition, Prentice Hall, 2002.

[34] P. Jäckel, Monte Carlo methods in finance, John Wiley & Sons, Inc., London, 2002.

[35] M. Jandaˇ

cka and D. Ševˇ

coviˇ

c, On the risk-adjusted pricing-methodology-based valu-

ation of vanilla options and explanation of the volatility smile, J. Appl. Math. 3(2005),

235–258.

[36] L. Jódar, P. Sevilla–Peris, J.C. Cortés, and R. Sala, A new direct method for solving

the Black–Scholes equation, Appl. Math. Lett. 18 (2005), 29–32.

[37] Y. Kabanov and M. Safarian, On Leland’s Strategy of Option Pricing with Transac-

tions Costs, Finance Stoch. 1(1997), 239–250.

[38] T. Kim and Y. Kwon, Finite–Difference Bisection Algorithms for free boundaries of

American Call and Put Options, Preprint Pohang University of Science and Technol-

ogy, 2002.

[39] M. Kratka, No Mystery Behind the Smile, Risk 9(1998), 67–71.

[40] Y.-K. Kwok, Mathematical Models of Financial Derivatives, Springer Finance, 1998.

[41] H.E. Leland, Option pricing and replication with transactions costs, J. Finance 40

(1985), 1283–1301.

[42] L. MacMillan, Analytic Approximation for the American Put Option, Adv. Fut. Opt.

Res. 1(1986), 119–139.

[43] B.J. McCartin and S.M. Labadie, Accurate and efficient pricing of vanilla stock op-

tions via the Crandall–Douglas Scheme, Appl. Math. Comp. 143 (2003), 39–60.

[44] R.C. Merton, Theory of rational option pricing, Bell J. Econ. Manag. Sci. 4(1973),

141–183.

[45] G. Meyer and J. van der Hoek, The valuation of American options with the method of

lines, Adv. Fut. Opt. Res. 9(1997), 265–286.

[46] B.F. Nielsen, O. Skavhaug and A. Tveito, Penalty and front–fixing methods for the

numerical solution of American option problems, J. Comp. Finance 5(2002), 69–97.

[47] R. Panini and R. P. Srivastav, Option Pricing with Mellin Transforms, Math. Comput.

Modelling 40 (2004), 43–56.

[48] O. Pauly, Numerical simulation of American Options, Diploma thesis, Universität

Ulm, 2004.

32 Julia Ankudinova and Matthias Ehrhardt

[49] J. Persson, Accurate Finite Difference Methods for Option Pricing, Dissertation, Up-

psala University, 2006.

[50] P. Schönbucher and P. Wilmott, The feedback–effect of hedging in illiquid markets,

SIAM J. Appl.Math. 61 (2000), 232–272.

[51] D. Ševˇ

coviˇ

c, Analysis of the free boundary for the pricing of an American call option,

Euro J. Appl. Math 12 (2001), 25–37.

[52] D. Ševˇ

coviˇ

c, An iterative algorithm for evaluating approximations to the optimal exer-

cise boundary for a nonlinear Black–Scholes equation, Canad. Appl. Math. Quarterly

15 (2007), 77–79.

[53] R. Seydel, Tools for Computational Finance, Second ed., Springer, 2004.

[54] H.M. Soner, S.E. Shreve and J. Cvitaniˇ

c, There is No Nontrivial Hedging Portfolio for

Option Pricing with Transaction Costs, Ann. Appl. Probab. 5(1995), 327–355.

[55] H.M. Soner and N. Touzi, Superreplication Under Gamma Constraints, SIAM

J. Contr. Optim. 39 (2001), 73–96.

[56] G. Strang, On the Construction and Comparison of Difference Schemes, SIAM J.

Numer. Anal. 5(1968), 506–517.

[57] J.C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wads-

worth & Brooks/Cole, 1989.

[58] H. Sun and J. Zhang, A high order compact boundary value method for solving one

dimensional heat equations, Tech. Rep. No. 333–02, University of Kentucky, 2002.

[59] D.Y. Tangman, A. Gopaul and M. Bhuruth, Numerical pricing of options using high–

order compact finite difference schemes, in press: J. Comput. Appl. Math. (2007).

[60] D.Y. Tangman, A. Gopaul and M. Bhuruth, A Fast High–Order Finite Difference Al-

gorithm for Pricing American Options, submitted to: J. Comput. Appl. Math.

[61] D. Tavella and C. Randall, Pricing financial instruments: The finite difference method,

John Wiley & Sons, 2000.

[62] P. Wilmott, J. Dewynne and S. Howison, Option Pricing: Mathematical Models and

Computation, Oxford, Financial Press, 1993.

[63] P. Wilmott, S. Howison, and J. Dewynne, The Mathematics of Financial Derivatives,

A Student Introduction, Cambridge University Press, 2002.

[64] L. Wu and Y. Kwok, A front-fixing finite difference method for the valuation of Amer-

ican options, J. Fin. Engin. 6(1997), 83–97.

Split–Step Schemes for Nonlinear Black–Scholes equations for American Options 33

[65] X. Wu and Z.-Z. Sun, Convergence of difference scheme for heat equation in un-

bounded domains using artificial boundary conditions, Appl. Numer.Math. 50 (2004),

261–277.

[66] J. Zhao, M. Davison and R.M. Corless, Compact finite difference method for American

option pricing, J. Comput. Appl. Math. 206 (2007), 306–321.

[67] Y. Zhu, X. Wu and I. Chern, Derivative securities and difference methods, Springer

Finance, New York, 2004.

[68] R. Zvan, P. Forsyth and K. Vetzal, A penalty method for American options with

stochastic volatility, J. Comp. Appl. Math. 91 (1998), 199–218.