A high order compact method for nonlinear Black-Scholes option pricing equations with transaction costs

(1)

A high-order compact method for nonlinear Black-Scholes option pricing equations with

transaction costs

Master’s Thesis in Financial Mathematics Ekaterina Dremkova

School of Information Science, Computer and Electrical Engineering Halmstad University

(2)

A high-order compact method for nonlinear Black-Scholes option pricing

equations with transaction costs

Ekaterina Dremkova

Halmstad University Project Report IDE0915

Master’s thesis in Financial Mathematics, 15 ECTS credits Supervisor: Priv.-Doz. Dr. Matthias Ehrhardt

Examiner: Prof. Ljudmila A. Bordag External referees: Prof. Daniel ˇSevˇcoviˇc

June 2, 2009

Department of Mathematics, Physics and Electrical Engineering School of Information Science, Computer and Electrical Engineering

Halmstad University

i

(3)

(4)

Contents

1 Introduction 1

2 Compact schemes 5

2.1 The valuation of options . . . . 5

2.2 High-order compact schemes . . . . 6

2.3 Classical finite difference schemes . . . 10

2.4 The compact schemes of higher order . . . 13

2.5 The fixed domain transformation . . . 18

2.6 The method of Liao and Khaliq . . . 21

3 The numerical solutions 25 3.1 Formulation of initial-boundary problem for American options in case of Black-Scholes model . . . 25

3.2 The grid for American options . . . 26

3.3 The free boundary problem . . . 26

3.4 The method of Liao and Khaliq for American Options . . . 27

3.5 The algorithm description . . . 29

4 Conclusions 31

Bibliography 33

iii

(5)

(6)

Chapter 1 Introduction

In the recent several years stock option was one of the most popular financial derivatives. There are many types of options on the market including Eu- ropean call (put) options, American call (put) options, Exotic options, etc.

But it was difficult to accurately price options until 1973 when Fisher Black and Myron Scholes published their famous Black-Scholes model in [4]. In an idealized financial market, the price of an European option can be obtained by analytically solving the Black-Scholes equation

0 = Vt+ 1

2σ²S²VSS + rSVS− rV, (1.1) where t - current date, σ - volatility, S - the price of the underlying asset, r - the risk-free interest rate, V - the option price, VS, VSS - the first and second derivatives of V with respect to S, Vt - the first derivative of V with respect to t. This is not very useful in practice, as mentioned in [9], because the Black-Scholes model had been derived under very restrictive assumptions, such as frictionless, liquid and complete market. In the real financial market the traders actually work in a different environment: transaction cost arising [3] and [5], the market is incomplete, etc.

Although the Black-Scholes model has been used in practice, it has also caused some criticism, for example because the volatility is not observable.

It is often deduced by calculating the so-called implied volatility out of the considered option prices by inverting the Black-Scholes formula. A widely observed unique property - the so-called smile/skew effect - is that, in the con- tradiction to the Black-Scholes model assumptions, these computed volatil- ities are not constant. This leads to a natural generalization of the Black- Scholes model replacing the constant volatility σ in the model by a local volatility function σ = σ(E, T ), where E denotes the exercise price, T - the maturity.

1

(7)

In practice, transaction costs arise when trading securities. Although they are generally small for institutional investors, they lead to a notable increase in the option price. In the past years, different models have been proposed to weaken unrealistic assumptions of the Black-Scholes model. These models result in fully nonlinear Black-Scholes equations.

In 1992, Boyle and Vorst [5] derived from a binomial model an option price taking into account transaction costs and that is equal to a Black- Scholes price but with a modified volatility of the form

σ = σ0(1 + cA)¹², A = µ σ0√

△T, c = 1. (1.2)

Here, µ is the proportional transaction cost, △T - the transaction period, and σ0 - is the original volatility constant. Leland [18] computed c = ²_π¹₂

. A more complex model has been proposed by Barles and Soner [3]. In their model the nonlinear volatility reads

σ² = σ₀²(1 + Ψ[exp(r(T − t)a²S²VSS)]), (1.3) where r is the risk-free interest rate, T - the maturity and a = µ√

γN in which γ is the risk aversion factor and N is the number of options to be sold.

The function Ψ is the solution of the nonlinear singular initial-value problem Ψ^′(A) = Ψ(A) + 1

2pAΨ(A) − A, A 6= 0, Ψ(0) = 0. (1.4) In the mathematical literature, only a few results can be found on the numerical discretization of Black-Scholes equation, mainly for linear Black- Scholes equations. The numerical discretization of the Black-Scholes equations with the nonlinear volatility (1.3) has been performed using explicit finite-difference schemes [3]. However, explicit schemes have the disadvantage that restrictive conditions on the discretization parameters (for instance, the ratio of the time and the space step) are needed to obtain stable and con- vergent schemes [24]. Moreover, the convergence order is only one in time and two in space. D¨uring et al [8,9] combined high-order compact difference schemes derived by Rigal [21] and techniques to construct numerical solutions with frozen values of the nonlinear coefficient of the nonlinear Black-Scholes equation

0 = Vt+ 1

2σ(VSS)²S²VSS+ rSVS− rV, (1.5) to linearize formulation.

As one can see, analytical solutions to nonlinear Black-Scholes equations are seldom available. We have to rely on numerical approaches such as binomial approximations, Monte-Carlo methods, finite element methods and

(8)

A compact method for nonlinear Black-Scholes equations 3

finite difference methods to get accurate option prices. Further we consider so called compact schemes for European and American options and focus on the transaction cost model of Barles and Soner (1.3). Instead of solving the singular differential equation (1.4) we propose to use some properties of Ψ = Ψ(A) described recently in [6]. As compact schemes cannot be directly applied to American type options and multi-dimensional problems, we will try to employ them using a fixed domain transformation proposed by Sevˇcoviˇc [2]. Finally, we consider the method of Liao and Khaliq [19] forˇ solving nonlinear Black-Scholes equation (1.5) with transaction costs.

(9)

(10)

Chapter 2

Compact schemes

2.1 The valuation of options

The financial derivatives used in trading are stocks, stock indices, foreign currency, debt instruments, goods and future contracts. There exists two types of vanilla options. A call option is a contract that gives the right to the holder to buy the underlying asset on a particular date at a specified value. A put option is a contract that gives the right to the holder to sell the underlying asset on a particular date at a specified value. The price in the contract is called exercise price (or stike price). The date in the contract is called expiration date (or exercise date, maturity date). Options are also divided into two types according to the expiration date. American options can be exercised at any time before maturity date. European options can only be exercised at the maturity date.

Remark 1. The terms ”American” and ”European” do not mean the place of the contracting but solely the type of the option.

Most options traded are American type options. Whereas European options are easier to analyze, some of the properties of American type options are carried over from the properties of European options.

Remark 2. It should be emphasized that options give the right to the holder to do something. But it does not mean that the holder has to do it.

This fact separates options from futures and forwards, where the holder has to buy or sell the underlying asset.

There exists two sides in every option contract. From one side there is an investor who has the long position, i.e. he has bought the option or he is the holder. From the other side there is an investor who has the short position, i.e. he has sold or he is the writer. The holder of an option receives the money immediately, but has the potential liabilities later. His gains or

5

(11)

losses are controversial to those of the writer.

Often it is beneficial to characterize the positions of European options through the profit of the holder. If E is the exercise price, ST - the final value of the underlying asset, then in European call option the long position will have the profit max(ST − E, 0). This shows that the option will be exercised if ST > E, and will not be exercised if ST ≤ K. The writer of the option will have the opposite value, i.e. in this case - max(ST − E, 0) = min(K − S^T). The profit of the long position in the European put option is equal to max(E − S^T, 0).

2.2 High-order compact schemes

In [9] Düring, Fournié and Jüngel presented an approach to use high-order compact schemes which needs a stencil of three points in space only. Consider the equation (1.5), where the nonlinear volatility is given by (1.3) with the terminal and boundary conditions

V (S, T ) = V0(S), S ≥ 0,

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

V (S, t) ∼ S − Ee^{r(t−T )}, S → ∞.

The last asymptotic condition can be rewritten in the form

S→∞lim

V (S, t)

S − Ee^{r(t−T )} = 1, (2.2)

uniformly for 0 ≤ t ≤ T.

In the Black-Scholes equation (1.5) we apply a time reversal substitution τ = T − t and obtain the following equation

0 = Uτ − S²

2 σ²USS − rSU^S+ rU, 0 < S < ∞, 0 < τ ≤ T. (2.3) The terminal condition in (2.1) takes the form

U(S, 0) = max(0, S − E), (2.4)

and the boundary conditions can be rewritten in the following way U(0, τ ) = 0, lim

S→∞

U(S, τ )

S − Ee^−rτ = 1. (2.5)

(12)

The derivatives of U with respect to S can be approximated at (S, τ ) by

∂U

∂S(S, τ )

= U(S − 2h, τ) − 8U(S − h, τ) + 8U(S + h, τ) − U(S + 2h, τ)

12h + O(h⁴),

(2.6)

∂²U

∂S²(S, τ )

= −U(S − 2h, τ) + 16U(S − h, τ) − 30U(S, τ) + 16U(S + h, τ) − U(S + 2h, τ) 12h²

+O(h⁴),

see [26]. We semidiscretize the nonlinear Black-Scholes equation (1.5) with respect to the variable S on a bounded interval [E − L, E + L], where E is the strike price and L is a suitable value such that E − L < 0. Grid lines are chosen as Si = E − L + ih, 0 ≤ i ≤ N, where Nh = 2L divides interval [E − L, E + L] into N equal subintervals. Let us rewrite expressions (2.6) at every mesh point Si. Doing so, we get a system of ordinary differential equations

du(τ )

dτ = Bu(τ ) + ω, (2.7)

where u(τ ) = [u1, u2, . . . , uN −1]^T ∈ Rⁿ⁻¹ with

B =







γ₁ δ₁ ξ₁ 0 0 0 0 . . . 0

β2 γ2 δ2 ξ2 0 0 0 . . . 0

α3 β3 γ3 δ3 ξ3 0 0 . . . 0 0 α₄ β₄ γ₄ δ₄ ξ₄ 0 . . . 0 . . . . . . . . . . . . . . . . . . .

0 . . . 0 0 αN −3 βN −3 γN −3 δN −3 ξN −3

0 . . . 0 0 0 α_{N −2} β_{N −2} γ_{N −2} δ_{N −2} 0 . . . 0 0 0 0 αN −1 βN −1 γN −1







, (2.8)

and ω is a column vector

ω =







α1u⁻1+ β1u0

α₂u₀ 0... 0 ξN −2uN

δ_{N −1}uN + ξ_{N −1}u_N+1







, (2.9)

(13)

where

αi = αi(τ ) = −σ_i²S_i²

24h² + rSi

12h, βi = βi(τ ) = 2σ²_iS_i²

3h² − 2rSi

3h , γi = γi(τ ) = −15σ²_iS_i²

12h² − r, (2.10)

δi = γi(τ ) = 2σ_i²S_i²

3h² +2rSi

3h , ξi = ξi(τ ) = −σ²_iS_i²

24h² − rSi

12h, where σi = σ(USS(Si, τ )) according to (1.3).

Remark 3. We should note that from (1.4) one gets σi = σ0 for the linear case, i.e. no transaction costs, a = 0, Ψ(0) = 0.

We need to search for an exact solution of the nonlinear volatility correction function Ψ - a unique solution of (1.4), see Theorem 3.1 of [3]. In order to solve (1.4), we distinguish two subdomains to the right and the left of the origin A = 0, where the initial condition Ψ(0) = 0 is given.

All the values u⁻1, u0, uN and uN+1 appearing in (2.9) must be specified because they correspond to approximations of the numerical solution at the values of the underlying variable

S−1 = E − L − h, S⁰ = E − L, S^N = E + L, SN+1 = E + L + h, (2.11) respectively. Unlike to the approach of [9] where the boundary conditions (2.5) are translated to the boundary of the bounded discretization domain, here using the semidiscretization technique we do not need the boundary values but they are obtained by interpolating the values of the numerical solution at the neighbour’s mesh internal nodes. Here, we use fourth-order Lagrange interpolating polynomial according to the order of approximation of scheme (2.6). Thus we have

u−1 = 10u1− 20u²+ 15u3− 4u⁴,

u₀ = 4u₁− 6u2+ 4u₃− u4, (2.12) uN = 4uN −1− 6u^{N −2}+ 4uN −3− u^{N −4},

uN+1 = 10uN −1− 20u^{N −2}+ 15uN −3− 4u^{N −4}.

Taking into account (2.7)-(2.9) and (2.12), the semidiscretized problem takes the form

du

dτ = M(τ )u(τ ), 0 < τ ≤ T, (2.13)

(14)

together with the initial condition

u(0) = [u1(0), . . . , uN −1(0)]^T, ui(0) = max(Si−E, 0), 1 ≤ i ≤ N −1, (2.14) where M(τ ) is given by







a1,1 a1,2 a1,3 a1,4 0 . . . . . . 0

a_2,1 a_2,2 a_2,3 a_2,4 0 . .. . .. 0

α₃ β₃ γ₃ δ₃ ξ₃ 0 . .. ...

0 α4 β4 γ4 δ4 ξ4 . .. ...

... ... ... ... . .. . .. . .. 0

... . . . 0 αN −3 βN −3 γN −3 δN −3 ξN −3

... . . . 0 aN −2,N −4 aN −2,N −3 aN −2,N −2 aN −2,N −1

0 . . . . . . 0 aN −1,N −4 aN −1,N −3 aN −1,N −2 aN −1,N −1





 ,

(2.15) where entry values are given by (2.10) and

a1,1 = γ1+ 10α1+ 4β1, a1,2 = δ1− 20α¹− 6β¹, a1,3 = ξ1+ 15α1+ 4β1, a1,4 = −4α¹− β¹,

a2,1 = β2+ 4α2, a2,2 = γ2− 6α², a_2,3 = δ₂+ 4α₂, a_2,4 = ξ₂− α2,

aN −2,N −4 = αN −2− ξ^{N −2}, aN −2,N −3 = βN −2+ 4ξN −2, aN −2,N −2 = γN −2− 6ξ^{N −2}, aN −2,N −1 = δN −2+ 4ξN −2,

aN −1,N −4 = −δ^{N −1}− 4ξ^{N −1}, aN −1,N −3 = αN −1+ 4δN −1+ 15ξN −1, aN −1,N −2 = βN −1− 6δ^{N −1}− 20ξ^{N −1}, aN −1,N −1 = γN −1+ 4δN −1 + 10ξN −1, (2.16) now are nonconstant and involve numerical approximations of USS, see [9].

The solution of (2.13) is given by

u(τ ) = e^Aτu(0). (2.17)

Using the Euler method the numerical solution of (2.17) and (2.13) takes the form

u(τ ) =

" _m=0 Y

m=ℓ−1

(I + kM(mk))

#

u(0), (2.18)

where k = △τ, lk = τ, M(τ) is given by (2.10) and (2.13) and

σ_i² = σ²₀(1 + Ψi(τ )), (2.19)

(15)

and

Ψi(mk) = Ψ

e^mkra²S_i² −ui−2+ 16u_i−1+ 16u_i+1− ui+2

12h²

, (2.20) uj evaluated at τ = mk [9].

Here instead of solving numerically the singular ODE Ψ^′(A) = Ψ(A) + 1

2pAΨ(A) − A, A 6= 0, Ψ(0) = 0, (2.21) we propose to use the following theorem presented recently in [6].

Theorem 1. The nonlinear volatility correction function Ψ, unique solution of (2.21) satisfies the following properties:

(i) Ψ is implicitly defined by

A = −arcsinhp(Ψ)

√Ψ + 1 +√ Ψ

!2

, if Ψ > 0, (2.22)

A = − arcsinp(−Ψ)

√Ψ + 1 −√

−Ψ

!2

, if 0 > Ψ > −1. (2.23)

(ii) Ψ is an increasing function mapping the real line onto the interval ] − 1, +∞[.

2.3 Classical finite difference schemes

We consider the nonlinear Black-Scholes equation (1.5) with the volatility (1.3) as proposed in [8]. In order to transform problem (1.5) into a convection- diffusion problem, we use the following transformation:

x(S) = ln S

E, t(r) = 1

2σ₀²(T − t), u = e^−xV E. Then (1.5) may be rewritten in the following form

ut− 1 + Φ[e^(Kt+x)a²E(uxx+ ux)] (uxx− u^x) − Ku^x= 0, (2.24) where x ∈ (−∞, ∞), 0 ≤ t ≤ T = ^σ²⁰₂^T, K = _σ^2r2

0. The transformed problem (2.24) has the following boundary and initial conditions

u(x, 0) = u0(x) = max(1 − e^−x, 0),

(16)

u(x, t) = 0 (x → −∞), (2.25)

u(x, t) ∼ 1 (x → +∞).

Here we will consider both standard difference schemes and compact schemes, developed by Rigal [22]. All considered difference schemes have two time levels. Let Aⁿ and Bⁿ be discretization matrices

Aⁿ= [a−1, a0, a1], Bⁿ= [b−2, b−1, b0, b1, b2],

where ai, bi - matrices’ main diagonals, superdiagonals and subdiagonals.

Then the schemes may be rewritten in the following form

AⁿUⁿ⁺¹ = BⁿUⁿ. (2.26)

The matrix Aⁿ is tridiagonal and the obtained linear systems may be solved using the Thomas algorithm [27]. Assume that

1

X

i=−1

ai =

2

X

i=−2

bi = 1.

Let us define the volatility correction in the following way σi = Ψ

exp (Knk + xi)a²E U_i−2ⁿ − 2Uiⁿ+ U_i+2ⁿ

4h² + U_i+1ⁿ − Ui−1ⁿ

2h

. (2.27) This formula gives an explicit discretization of nonlinearity and uses a special stencil for the second derivative (step 2h instead of h).

Another problem lies in the initial condition for u0, as it is nondifferentiable in the point x = 0. Oosterlee et al. [20] solved this problem of reduced accuracy and proposed a grid stretching technique, which is based on an idea of placing more points in the neighborhood of the nondifferentiable payment condition. We use interpolation of high order to smooth the initial data but only with approximation (2.27) useful results can be obtained.

The Figure 2.1 shows solutions of equations (2.15), (2.16) (on the left), and of the ODE (2.14) (on the right) and their spline interpolation. It is easily seen that figures show indistinguishable results. Results are obtained using the standard Matlab routine fsolve.

For further investigation and comparison of the schemes we need to introduce some values and give some definitions.

Remark 4. Stability of a finite difference scheme means that in a relatively small time period approximate and precise values have a small difference.

Let us also introduce the notations:

(17)

Figure 2.1: Solutions of equations (2.15), (2.16) (left figure) and ODE (2.14) (right figure, for A > 0), where T = 1, K = 10, r = 0.1, q = 0.05.

1. λ = −(1 + K) - the linear part of the coefficient of the convection term in (2.24);

2. α = ^λh₂ - the cell Reynolds number;

3. r = _h^k2 - the parabolic mesh ratio;

4. µ = _h^k - the hyperbolic mesh ratio.

The Forward-Time Central-Space explicit scheme (FTCS).

This scheme is given by

a⁻1 = 0, a0 = 1, a1 = 0,

b−1 = r − µ

2(σi− λ), b⁰ = 1 − 2r −r

2σi, b1 = r + µ

2(σi− λ), b−2 = b2 = r

4σi.

It is of order 1 in space, 2 in time, with a very strict stability condition r ≤ 1

2. (2.25)

The condition

|α| ≤ 1, (2.26)

where α is a cell Reynold’s number, must be satisfied to avoid oscillations.

(18)

The Backward-Time Central-Space semi-explicit scheme (BTCS).

This scheme is a scheme with an explicit treatment of the nonlinearity and it is given by

a⁻1 = λ

2µ − r, a0 = 1 + 2r, a1 = −λ 2µ − r, b⁻2 = r

4σi, b⁻1 = −1

2µσi, b0 = 1 − r

2σi, b1 = 1

2µσi, b2 = r 4σi, where λ, µ, r are values explained in the notation above.

This scheme is of order 1 in space and 2 in time. It is unconditionally stable and if (2.26) is satisfied, then it is non-oscillatory.

The Crank-Nicolson scheme (CN).

This scheme, with an explicit treatment of the nonlinearity, is given by a−1 =

−r 2+ µ

4

σi− r 2 − λ

4µ, b−1 =r 2 − µ

4

σi +r 2 +λ

4µ, a0 = 1 + r(1 + σi), b0 = 1 − r(1 + σⁱ),

a1 =

−r 2− µ

4

σi− r 2+ λ

4µ, b1 =r 2 +µ

4

σi+ r 2− λ

4µ,

and b−2 = b₂ = 0. It is of the order 2 both in time and space and unconditionally stable.

2.4 The compact schemes of higher order

In [22] Rigal introduced several finite difference schemes (FDS) for linear convection-diffusion problems. We consider only two of them - R3A and R3B - and apply them to the problem (2.24). These schemes are both compact two-level schemes of order 2 in space and 4 in time in the linear case.

In R3 methods several propositions were made. The class of the two-level three-point schemes of the order (2,4) is defined in the following way

a−1vⁿ⁺¹_i−1 + a₀v_iⁿ⁺¹+ a₁v_i+1ⁿ⁺¹= b−1v_i−1ⁿ + b₀v_iⁿ+ b₁v_i+1ⁿ . (2.27) Matrices A and B are positive and all their entries are positive too. For the FDS (2.27) to be positive matrix A⁻¹B should be positive too. A and B are matrices with diagonals a0 and b0, superdiagonals a1 and b1, subdiagonals a⁻1 and b⁻1 respectively.

Definition 1. A positive matrix is a matrix in which all the elements are greater than zero.

(19)

For the general construction and description of properties of the fourth- order schemes we should refer to two lemmas that correspond with two-level three-point schemes (2.27), [22]. Consider the model diffusion-convection problem P

∂tu = ∂_x²u − λ∂^xu + f = Au + f in ]0, 1[×[0, T [.

We assume that

1

X

i=−1

ai =

1

X

i=−1

bi = 1,

which can always be obtained after a possible normalization of the coefficients in the schemes.

Lemma 1 [21]. The FDS (2.27) is stable iff the coefficients ai, bi fulfill (a₁− a⁻1)²− (b1− b⁻1)² > a₁+ a−1− b1− b⁻1, (2.28) (a1− a⁻¹)²− (b¹− b⁻¹)² > a1+ a−1− b¹− b⁻¹. (2.29) Lemma 2 [21]. The FDS (2.27) is non-oscillatory if the coefficients ai, bi fulfill the condition

(a1− b¹)(a⁻1− b⁻¹) ≥ 0. (2.30)

The general two-level three-point scheme (P_h) is defined in the following way

(1 + C)Dtv_iⁿ= 1 2+ A1

D+D−v_iⁿ+ 1 2+ A2

D+D−v_iⁿ⁺¹ (2.31)

−λ 1 2 + B1

D0v_iⁿ− λ 1 2+ B2

D0vⁿ⁺¹_i ,

where A1, A2, B1, B2, C are real constants chosen in the way to eliminate the lower terms in the truncation error. Dt, D0, D+, D⁻ are the standard difference operators

Dtv_iⁿ= v_iⁿ⁺¹− vⁿi

△t , D0v_iⁿ= v_i+1ⁿ − vⁿi−1

2h ,

D₊v_iⁿ= v_i+1ⁿ − vⁿi

h , D−vⁿ_i = vⁿ_i − vi−1ⁿ

h .

(2.32)

Remark 5. The truncation error or the local truncation error is the error made by numerical algorithms that arises from taking finite number of steps in

(20)

the computation. It is present even with infinite-precision arithmetic, because it is caused by truncation of the infinite Taylor series to form the algorithm.

In order to get the truncation error, we need to apply (Ph) to the homo- geneous problem u, associated with (P) :

Eu(△t, h) = (1 + C)D^tu(xi, tn) (2.33)

− 1 2 + A2

D+D⁻u(xi, tn+1)

− 1 2 + A1

D+D−u(xi, tn)

−λ 1 2 + B₂

D₀u(xi, t_n+1)

−λ 1 2 + B1

D0u(xi, tn).

Here u(x, t) satisfies

∂tu + λ∂xu = ∂_x²u. (2.34) By decomposition each term in (2.29) we get

Eu(△t, h) =

6

X

i=1

ei∂_xⁱu + HOD (higher order derivatives), (2.35)

where

e1 = λ(B1+ B2− C), (2.36)

e2 = C − A¹− A²− λ²△t

B2−C 2

, (2.37)

e3 = λ

△t(A²+ B2− C) + h²

6 (1 + B1+ B2) (2.38) +λ²△t²

6

1

2 + 3B2− C

,

e4 = △t C 2 − A²

−(1+A¹+A2)h²

12+λ²△t² 2

C − 1

2− A²− 2B²

(2.39)

− 1 2+ B₂

λ²h²△t

6 + λ⁴△t³

24 (C − 4B2− 1), e5 = λ △t²

2

2A2+ B2+ 1 2− C

+ h⁴

120(1 + B1+ B2) (2.40)

(21)

+h²△t 12

3

2 + A2+ 2B2

,

e₆ = △t² 2

C − 3A2− 1 2

− (1 + A1+ A₂) h⁴

360 −h²△t 12

1 2+ A₂

. (2.41) The class of the R3 schemes is defined by

e1 = e2 = e3 = 0, (2.42)

ER3 = O(△t²+ h⁴), (2.43)

that means that e₁, e₂, e₃ and the terms of order less than (2,4) in e₄ should dissappear.

As mentioned before Ai, Bi, C are real constants chosen in the way to eliminate the lower terms in the truncation error. We use C = 0 and considering (2.38), we express A₁, A₂, B₁ as functions of B₂

B1 = −B², (2.44)

A2 = −λ²△t 12 − 1

6r − B²

1 + λ²△t 2

, (2.45)

A1 = λ²△t 12 + 1

6r + B2

1 − λ²△t 2

. (2.46)

The value of B2 must be chosen in such way that

e4 = h²

12− λ²△t² 6 + B2

12[−λ²h²△t + 12△t + λ⁴△t³] (2.47) is of order (2,4).

The R3A scheme.

In this scheme the value of B2 is chosen as B2 = − 1

12r, (2.48)

which eliminates in e4 the only terms depending on △t and h². Hence,

e4 = −λ²△t²

6 +λ²h⁴

144 − λ⁴h²△t² 144 .

(22)

Replacing B2 in expressions of A1, A2 (2.45), (2.46) the coefficients will be written in the following way

a−1 = 1 12− r

2

(1 + α) − α²r

6 + α²r² 3 , a₀ = 5

6+ r + α²r

3 − 2α²r² 3 , a1 = 1

12− r 2

(1 − α) −α²r

6 +α²r² 3 , b⁻2 = r

4σi, b−1 = 1

12+ r 2

(1 + α) + α²r

6 +α²r² 6 − 1

2µσi, b0 = 5

6− r − α²r

3 − 2α²r² 3 − r

2σi, b₁ = 1

12+ r 2

(1 − α) − α²r

6 + α²r² 3 + 1

2µσi, b2 = r

4σi.

It is stable in the linear case σi = 0 if r ≤ 1

√2|α|, (2.49)

cf. [22]. If α is arbitrary, then this scheme is non-oscillatory, cf. Rigal [22].

The R3B scheme.

The value B₂ for this scheme is defined in the following form

B₂ = − 1

12r − λ²△t

12 . (2.50)

From (2.44)-(2.46) we have

B1 = −B², A₁ = 1

12r +λ²h²

24 +λ⁴△t²

24 , (2.51)

A2 = − 1

12r +λ²h²

24 + λ⁴△t² 24 .

(23)

The coefficients for this scheme defined by (2.51) are the following:

a−1 = 1 12− r

2

(1 + α) − α²r

6 + α²r²

3 − 2α⁴r³ 3 , a0 = 5

6+ r + α²r

3 + 4α⁴r³ 3 , a1 = 1

12− r 2

(1 − α) − α²r

6 − α³r²

3 − 2α⁴r³ 3 , b⁻2 = r

4σi, b−1 = 1

12+ r 2

(1 + α) + α²r

6 +α³r²

3 +2α⁴r³ 3 − r

4+1 2µ

σi, b0 = 5

6− r − α²r

3 −4α⁴r³

3 − 2rσⁱ, b1 = 1

12− r 2

(1 + α) + α²r

6 − α³r²

3 + 2α⁴r³ 3 − r

4 −1 2µ

σi, b₂ = r

4σi.

It is unconditionally stable and non-oscillatory in the linear case σi = 0 [22].

2.5 The fixed domain transformation

Compact schemes which many authors applied to the Black-Scholes equation with transaction costs have one disadvantage: these schemes cannot be generalized to multi-dimensional problems, and are (directly) applicable to European type options only. However, with the fixed domain transformation introduced by ˇSevˇcoviˇc [2], [23] we overcome this shortcoming.

We consider the nonlinear Black-Scholes equation with dividend yield 0 = Vt+ 1

2σ˜²(t, S, VS, VSS)S²VSS+ (r − q)SV^S− rV, (2.52) where dividend yield q is constant, S > 0, t ∈ (0, T ). This equation is supplied with the following terminal and boundary conditions

V (S, T ) = (S − E)⁺ for 0 ≤ S ≤ S^f(T ),

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

V (Sf(t), t) = Sf(t) − E for 0 ≤ t ≤ T, VS(Sf(t), t) = 1 for 0 ≤ t ≤ T, Sf(T ) = max(K, rK/q).

(2.53)