Using an explicit scheme for an application of finite difference methods may lead to stability issues. If one wants to increase the accuracy by raising the number of spatial grid points, the number of time intervals have to be increased to a certain extent in order to sustain a converging behavior.

As for quite accurate results ridiculously many grid points in time are needed, the practical use of the explicit scheme is rather limited due to high computational effort. Implicit methods for finite difference methods are designed to overcome these stability limitations imposed by the already mentioned convergence restrictions. Since such methods are unconditionally stable, both accuracy and limited computational effort can be combined.

This text offers an introductory treatment of Finite Difference Methods employing an implicit scheme. It includes a theoretical derivation of the implicit scheme and the Crank-Nicolson scheme, a numerical application to European puts as well as a theoretical discussion and comparison of the truncation error for both schemes. Finally, Richard-Extrapolation is introduced as a nice tool for lowering the truncation error.

Excerpt

Contents

List of Tables

List of Figures

Theoretical Aspects of Implicit Finite Difference Methods

1.1

Derivation of a Fully Implicit Scheme

. . . . . . . . . . . . . . . . . . . . . .

1.1.1

Exemplary Boundary Conditions for a European Put . . . . . . . . . .

1.1.2

The Fully Implicit Scheme in Matrix Notation

. . . . . . . . . . . . .

1.2

Solving Tridiagonal System of Linear Equations . . . . . . . . . . . . . . . . .

1.3

Derivation of the Crank-Nicolson Scheme

. . . . . . . . . . . . . . . . . . . .

1.3.1

Exemplary Boundary Conditions for a European Put . . . . . . . . . .

1.3.2

The Crank-Nicolson Scheme in Matrix Notation

. . . . . . . . . . . .

Numerical Application for a European Put

2.1

Results and Interpretations . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Technical Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Truncation Error and Rate of Convergence

3.1

Richardson Extrapolation for the Fully Implicit Scheme . . . . . . . . . . . .

3.2

Richardson Extrapolation for Crank-Nicolson Scheme

. . . . . . . . . . . . .

3.3

Preliminary Efficiency Analysis . . . . . . . . . . . . . . . . . . . . . . . . . .

3.4

Technical Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Appendix

Spurious oscillations for Crank-Nicolson . . . . . . . . . . . . . . . . . . . . . . . .

Problems of Crank-Nicolson for non-smooth Boundaries . . . . . . . . . . . . . . .

Continued Efficiency Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

References

List of Tables

Results for Richardson Extrapolation Implicit Scheme . . . . . . . . . . . . .

Results for Richardson Extrapolation CN Scheme . . . . . . . . . . . . . . . .

List of Figures

Errors of numerical solutions (U

- u

BSM

(S, = T )) depending on S . . . . .

3-D Plot for Error of Numerical Solutions depending on S and

. . . . . . .

Errors of numerical solutions (U

- u

BSM

(S, = T )) depending on S (zoom) .

Approximation Error for Implicit Scheme and Richardson Extrapolation . . .

Comparison of Approximation Error for Implicit Scheme . . . . . . . . . . . .

Approximation Error for CN Scheme and Richardson Extrapolation

. . . . .

Comparison of Approximation Error for CN Scheme . . . . . . . . . . . . . .

Comparison of Approximation Error for both Schemes with RE . . . . . . . .

Influence of Total Number of grid points using relation N =

. . . . . .

Maximum Absolute Error as a Function of Computational Time with N =

Influence of Total Number of grid points using relation N =

1.2

. . . . . .

Maximum Absolute Error as a Function of Computational Time with N =

1.2

M 21

Comparison of Maximum Absolute Error as a Function of Computational Time 22

Comparison of the Effect of varying Grid Points on Approximation Error

. .

Comparison of Maximum Absolute for a fixed M = 120 . . . . . . . . . . . .

Spurious Oscillation of CN method for M = 2000 and N = 50 . . . . . . . . .

Spurious Oscillation of CN method for M = 2000 and N = 50 3D plot . . . .

Righthand-Boundary Problem of CN method M = 500 and N = 1000 3D plot

Approximation error for fixed N and varying M

. . . . . . . . . . . . . . . .

Convergence of the Approximation Error for N = 100 and varying M

. . . .

Approximation error for fixed M and varying N

. . . . . . . . . . . . . . . .

Theoretical Aspects of Implicit Finite Difference Methods

Using an explicit scheme for an ap-plication of finite difference methods may lead to stability

issues.

If one wants to increase the accuracy by raising the number of spatial grid points, the

number

of time intervals have to be increased to a certain extent in order to sustain a

converging

behavior. As for quite accurate results ridiculously many grid points in time are

needed,

the practical use of the explicit scheme is rather limited due to high computational

effort.

Implicit methods for fi-nite difference methods are designed to overcome these stability

limitations

imposed by the already mentioned convergence restrictions. Since such methods

are

unconditionally stable, both accuracy and limited computational effort can be combined.

1.1

Derivation of a Fully Implicit Scheme

As a starting point, assume that the price of an asset follows a geometric Brownian motion:

= (r

- q)S

dt + S

(1)

According to Hirsa (2013, p. 115), the value of a derivative on this financial asset specified

in (1) must satisfy the Black-Scholes-Merton PDE subject to a given terminal condition at

every point in time:

(S, t) + (r

- q)SF

(S, t)

- rF (S, t) + F

(S, t) = 0

s.t. F (S, T ) = f (S, T ) .

(2)

In the PDE, is the so-called diffusion coefficient of the underlying Brownian motion, S the

price of the underlying asset, i.e. a stock, r the risk-free interest rate and q the continuously-

paid dividend rate.

In order to obtain a numerical solution for the above PDE, the initial variables should be

transformed. Furthermore, the continuous pricing problem has to be converted into a discrete

one by dividing the feasible time and space domain into a rectangular grid and substituting the

partial derivatives with differential quotients. Following this idea, a new set of independent

variables is introduced:

x = S = ih

i [0, M]

= T

- t = nk n [0, N]

where i and n are index numbers for the discretized grid points, and h and k are step sizes

in space and time direction, respectively. Consequently, the rectangular grid reads:

[0, x

max

]

with x

max

= M

· h

[0, T ] with T = N · k

Subsequently, the financial derivative and the terminal condition are also transformed

F (S, t) = u(x, )

f (S, T ) = f (x, 0)

By the definition of , the terminal condition becomes an initial condition.

As a final step, the PDE as given in (2) is discretized with the newly introduced variables

x = ih and = nk

a(i, n)u

(i, n) + b(i, n)u

(i, n) + c(i, n)u(i, n)

- u

(i, n) = 0

s.t. u(x, 0) = f (x, 0)

(3)

By comparing (2) and (3) one may simply define

a(i, n) = a(i) = 0.5 i

b(i, n) = b(i) = (r

- q)ih c(i, n) = c = -r

(4)

Since the functional form of the underlying partial derivatives in (3) is not known, they are

approximated by differences using Taylor series expansions. For the implementation of the

fully implicit scheme, the following approximations around x = ih and = (n + 1)k are used:

· u

(i, n + 1) =

n+1

i+1

- 2U

n+1

+ U

n+1

-1

) + O(h

)

· u

(i, n + 1) =

n+1

i+1

- U

n+1

-1

) + O(h

)

· u

(i, n + 1) =

n+1

- U

) + O(k)

· u(i, n + 1) = U

n+1

+ O(k)

In contrast to the explicit scheme, the fully implicit scheme employs a backward approxima-

tion for the partial time derivative. Plugging these approximations into (3) and dropping the

error terms leads to

(

n+1

i+1

- 2U

n+1

+ U

n+1

-1

)

(

n+1

i+1

- U

n+1

-1

)

+ cU

n+1

(

n+1

- U

)

= 0 ,

which can be rearranged for an implicit solution of the values at the next time level:

(i, n + 1)U

n+1

-1

+ d

(i, n + 1)U

n+1

+ d

(i, n + 1)U

n+1

i+1

= d

(i, n + 1) = U

(5)

where

(i, n + 1) = d

(i) =

-a

+ b

k +

- q)ik

(i, n + 1) = d

(i) = a

- ck + 1 = 1 + rk +

(i, n + 1) = d

(i) =

-a

- b

- q)ik

In contrast to an explicit scheme, it is not possible to compute U

n+1

directly from the already

known values because three unknown values are linked to only one known value. However,

since (5) has to hold for all i = 1, ..., M

-1, a system of M -1 equations emerges for the M +1

unknown values. Appropriate boundary conditions yield the two missing values for each time

step and the terminal conditions give the values in the first time layer. Thus, after adding

two additional equations (von-Neumann boundaries) or fixing two values in an appropriate

manner (Dirichlet boundary), the N

- 1 linear equation systems can be solved recursively in

order to calculate the corresponding derivative price at = T .

In the functional formulation, the step sizes h and k will be omitted as they are the same for all values.

The index numbers i and n determine the underlying dynamics.

Note that all newly introduced coefficients are independent of the time index n.

1.1.1

Exemplary Boundary Conditions for a European Put

The choice of boundary conditions can be crucial for ensuring the accuracy of the pricing

algorithm given in (5).

This is true as any error on a boundary is propagated through

the ongoing finite difference scheme. This makes the entire scheme quite sensitive to the

boundary conditions (Hirsa, 2013, p. 121). Hence, a sensible and complete formulation

of these conditions at i = 0 and i = M is the key to a successful implementation of the

underlying pricing problem.

If a European put is the underlying derivative for the recursive pricing algorithm given in

(5), one can simply specify the starting point for the implicit scheme by taking the terminal

values for a plain-vanilla European put:

f (S, 0) = u(x, 0) = max(E

- S, 0)

= max(E

- ih, 0) i [0, M] .

(6)

For the left-sided boundary, i.e. the stock price is 0, a Dirichlet boundary can be used since

it can be assumed that the option value just equals the discounted strike price:

= Ee

-r(nk)

n [1, N] .

(7)

Hence for i = 1, the equation system in (5) reads

(1, n + 1)U

n+1

+ d

(1, n + 1)U

n+1

= d

(1, n + 1)

- d

(1, n + 1)Ee

-r(n+1)k

(1,n+1)

(8)

As there is no "natural" bound on the right side of the grid, an artificial one needs to be

introduced. It can be assumed that for really high stock prices, the option value does not

change if the price of the underlying stock changes slightly, i.e. the Delta of the European

put is equal to 0 if i = M . Thus, the central approximation for the first derivative around

i = M can be used and set equal to zero

n+1

(

n+1

M +1

- U

n+1

-1

)

= 0 .

(9)

Additionally, (5) is also valid at i = M

(M, n + 1)U

n+1

-1

+ d

(M, n + 1)U

n+1

+ d

(M, n + 1)U

n+1

M +1

= d

(M, n + 1) .

(10)

Solving (9) for U

n+1

M +1

= U

n+1

-1

gives the possibility to eliminate U

n+1

M +1

in (10) in order to

obtain

(M, n + 1) + d

(M, n + 1))

(M,n+1)

n+1

-1

+ d

(M, n + 1)U

n+1

= U

= d

(M, n + 1)

(11)

The expression in (11) adds an additional equation to the linear system of equations such

that there are now M equations for M unknowns at each time layer n

. At first glance,

it seems that the implicit scheme requires much more computational effort compared to an

explicit scheme as a system of equations needs to be solved for every n. However, if the

implicit scheme is written in matrix notation, it exhibits a quite nice structure. As shown

below, the matrix is tridiagonal which simplifies the solution of the system dramatically.

Note that the value U

n+1

is fixed via a Dirichlet boundary and is therefore not an unknown anymore.

1.1.2

The Fully Implicit Scheme in Matrix Notation

The M equations given in (5), (8) and (11) can be written in compact form using matrix

notation:

×M

× U

n+1

×1

(12)

(1)

(2)

. ..

- 1) d

- 1)

(M )

n+1

-1

n+1

(1)

(2)

- 1)

(M )

The straightforward method to solve (12) for U

n+1

involves taking the inverse of A. However,

in practice, there are far more efficient solution techniques than matrix inversion. The matrix

A has the property that it is tridiagonal, which means that only the diagonal, the super-

diagonal and the sub-diagonal elements are non-zero. For such special systems of equations,

two efficient methods are outlined in the next Section.

1.2

Solving Tridiagonal System of Linear Equations

In order to solve (12) for U

n+1

, the tridiagonal structure of the underlying system can be

exploited. The first proposed method involves a matrix decomposition into an upper and

lower matrix which makes it possible to subdivide the given problem into two minor problems.

Assume A

× U

n+1

= d

is decomposed into

× U

n+1

= (L

× R)U

n+1

= L(R

× U

n+1

) = L

× y = d

(13)

This means that the first step employs the so-called Cholesky-decomposition according to

A = L

× R

. ..

-1

. ..

-1

. (14)

From (14) it follows for the first values of and r

(1) =

= d

(1)

(15)

(1) =

(1)

In order to keep it well-arranged, the spatial indices are omitted. However, it should be kept in mind

that the values for d

, d

and d

are different in every row of A.

and then in general

-1

- 1)

-1

i = 2, ..., M

= d

(i)

- d

(i)r

-1

i = 2, ..., M

(16)

= d

(i)

i = 2, ..., M

Due to the fact that all values in A are time-invariant, this decomposition only needs to be

done once as it is the same for every time step

. Subsequently, L

× y = d

can be solved by

forward iteration. As d

is time-dependent, this step has to be done for every n.

. ..

-1

(1)

(2)

- 1)

(M )

(17)

According to (17), it follows directly that

= d

(1)

(18)

(i)

- d

(i)y

-1

)

i = 2, ..., M

(19)

which can be used for a backward iterative solution of R

× U

n+1

= y

. ..

-1

n+1

-1

n+1

-1

(20)

Using (20), one can finally solve for the unknown values U

n+1

= y

(21)

n+1

= y

- r

n+1

i+1

= y

(i)

n+1

i+1

i = M

- 1, ..., 1

(22)

Although this workaround seems quite time-consuming, solving a tridiogonal system of equa-

tion with the above algorithm requires only a little bit more computational effort compared to

an explicit scheme, to be exact twice as many operations per time step (Wilmott et al., 1995,

p. 144). Yet it can be shown that the implicit scheme is unconditionally stable which means

that it is not necessary to increase N to really large numbers for a finer spatial grid

. This

more than compensates the loss in time for solving the above linear system. Thus, implicit

methods can be considered superior to explicit methods.

It should be noted that actually only (15) and (16) are needed for later steps. Hence, only these values

must be stored somewhere.

A proof of the unconditional stability of the implicit scheme is given by Brandimarte (2006, p. 311).

For the sake of completeness, an algorithm for a Gauss elimination is given as a second

method in the following. As for the first method, there are some calculations which do not

involve d

and are therefore time-invariant:

(i) :=

(i)

- 1)

i = 2, ..., M

(23)

(i) := d

(i)

- d

(i)d

- 1)

i = 2, ..., M .

(24)

As a next step, one needs to do a forward iteration according to

(i) := d

(i)

- d

(i)d

- 1)

i = 2, ..., M

(25)

and finally a backward iteration which results in

n+1

(M )

(26)

n+1

(i)

- d

(i)U

n+1

i+1

(i)

i = M

- 1, ..., 1 .

(27)

As for the first method, the forward and backward iterations have to be done on each time

layer n. The implementation of both methods using MATLAB and a comparison with respect

to computational time are given below.

1.3

Derivation of the Crank-Nicolson Scheme

So far, the explicit and the implicit scheme were derived using a forward or a backward

approximation in time, respectively. As a central approximation for the first derivative in

space was proven to be more accurate than a forward or backward iteration, one could follow

the same line of argumentation and simply combine them to a central approximation with

respect to time. However, such approximations with step size k in each direction usually

result in very unstable solutions which make them defective for practical use (Wilmott et al.,

1995, p. 138).

This gave rise to the Crank-Nicolson (CN) method by Crank & Nicolson (1947) who managed

to combine a forward and backward approximation in a feasible way in order to make use of a

more accurate central approximation in time. The basic idea is to consider a grid point which

is actually not on the underlying grid but positioned halfway between two time layers

. This

means that (2) is discretized around x = ih and = (n +

)k. For a discretization around

this point, there is no other possibility than using a central approximation since a forward or

backward approximation would involve a differential quotient with a grid point outside the

actual grid. For a derivation of this central approximation, one needs to consider a Taylor

series expansion in time direction using an appropriate time step. Obviously the appropriate

time step is

(

)

= u(x, )

(

)

(

)

+ . . . .

(28)

For the original CN method the grid point is halfway between n and n + 1. However, there are also other

methods which vary the positioning between two time layers.

Excerpt out of 32 pages - scroll top

Details

Title: Finite difference methods with an implicit scheme
College: University of Tubingen
Grade: 1,0
Author: Pascal Sturm (Author)
Publication Year: 2016
Pages: 32
Catalog Number: V373359
ISBN (eBook): 9783668526068
ISBN (Book): 9783668526075
Language: English
Tags: Finite Difference Methods Implicit Scheme Crank-Nicolson-Scheme

Quote paper: Pascal Sturm (Author), 2016, Finite difference methods with an implicit scheme, Munich, GRIN Verlag, https://www.grin.com/document/373359

Finite difference methods with an implicit scheme

Excerpt

Details