## Option Pricing Using The Crank-Nicolson Finite Difference Method

This tutorial discusses the specifics of the Crank-Nicolson finite difference method as it is applied to option pricing. Example code implementing the Crank-Nicolson method in MATLAB and used to price a simple option is given in the Crank-Nicolson Method - A MATLAB Implementation tutorial.

The Finite Difference Methods tutorial covers general mathematical concepts behind finite diffence methods and should be read before this tutorial. Alternative finite difference methods, namely the implicit method and the explicit method, are covered in companion tutorials.

### Discretizing the Black-Scholes-Merton PDE

The Crank-Nicolson finite difference method represents an average of the
implicit method and the
explicit method.
Consider the grid of points shown in Figure 1.
This represent a small portion of the general
pricing grid used in
finite difference methods.
Indices *i* and *j* represent nodes on the
pricing grid.

Whereas the explicit method
prices the node *ƒ _{i-1, j}* (the central node in the left
hand side) based on the values of

*ƒ*,

_{i, j+1}*ƒ*and

_{i, j}*ƒ*(the nodes on the right hand side), and the implicit method nominally prices the nodes

_{i, j-1}*ƒ*,

_{i-1, j+1}*ƒ*and

_{i-1, j}*ƒ*(the three nodes on the left hand side) based on the value of

_{i-1, j-1}*ƒ*(the central node in the right hand side), the Crank-Nicolson method prices all three of the left side nodes based on the values of all three of the right side nodes. Both the implicit method and the Crank-Nicolson method achieve this by solving a set of

_{i, j}*M-1*simultaneous equations.

To derive the Crank-Nicolson difference equations consider the
node *ƒ _{i-1/2, j}* which lies at the centre
of Figure 1.
(Note that a price will not be calculated for this node.
Rather it is used as a mathematical convenience that will not appear in the final
equations.)

The goal is to discretize the Black-Scholes-Merton partial differential equation,

To achieve this (at node *ƒ _{i-1/2, j}*) the following approximations
are used.

- use a central approximation
for
*∂ƒ/∂t*(Compare this with the explicit method where the backward approximation is used, and the implicit method where the backward approximation is used) - use a central approximation
for
*∂ƒ/∂S* - use a standard approximation
for
*∂*^{2}ƒ/∂S^{2}

where the indices *i* and *j* represent nodes on the
pricing grid.

Substituting these approximations into the Black-Scholes-Merton PDE and collecting like terms this reduces to

where

When Equation 1 is written for all values of *i* and *j*
it leads to a set of *M-1* equations in *M-1* unknowns.
Hence the value for *ƒ* at each node can be calculated uniquely by
solving this set of simultaneous equations.
(This can arguably be seen easier using the matrix notation shown in the following
A Matrix Formulation section.)

In the option pricing framework, given the option payoff at expiry nodes then the
prices *Δt* before expiry
can be calulated, then from those prices the value *2Δt* before expiry can be
calculated, and working iteratively backwards through time until the option
price at grid nodes for *t=0* (i.e. today) can be calulated.

### A Matrix Formulation

The formulation for the Crank-Nicolson method given in Equation 1 may be written in the matrix notation

where

and

This matrix notation is used in the Crank-Nicolson Method - A MATLAB Implementation tutorial.

### Stability and Convergence

Two important questions to ask about any numerical algorithm are *when is it stable?* and
if it's stable then *how fast does it converge?*
(An iterative algorithm that is unstable will lead to the calculation of ever increasing numbers
that will at some point approach infinity.
On the other hand, a stable algorithm will converge to a finite solution.
Typically the faster that finite solution is reached the better the algorithm.

From standard results in matrix algebra it is known that a matrix equation of the form given in Equation 3 is stable if and only if

Equation 4 shows the *infinity norm* of the product of the matrices *C ^{-1}D*.
Heuristically, if the infinity norm of

*C*is less than 1 then successive values of

^{-1}D*F*in Equation 3 get smaller and smaller, and hence the algorithm converges, or is stable. (Alternatively if the infinity norm of

_{i}*C*is greater than 1 then successive values of

^{-1}D*F*get larger and larger and hence diverge.)

_{i}
It can be shown that the infinity norm of *C ^{-1}D* is less than 1 for all
values of

*ρ*,

*σ*and

*δt*. Hence the Crank-Nicolson finite difference method is always stable. (Compare this with the explicit method which can be unstable if

*δt*is chosen incorrectly, and the implicit method which is also guaranteed to be stable.)

The rate of convergence of the algorithm is directly related to the truncation error
introduced when approximating the partial derivatives.
Hence the Crank-Nicolson method converges at the rates of *Ο(δt ^{2})*
and

*Ο(δS*. This is a faster rate of convergence than either the explicit method, or the implicit method.

^{2})
A disadvantage of the Crank-Nicolson method (over the
explicit method)
is that it requires the inverse of a matrix
(i.e. *C ^{-1}*) to be calculated, and the inverse of a matrix is
(computationally) an expense operation to perform.
Fortunately, for tri-diagonal matrices such as

*C*(i.e. matrices that only have non-zero elements down their diagonal and the terms directly above and below their diagonal) fast inversion algorithms are available.

Another issue with the Crank-Nicolson method is that it is known to be sensitive
to non-smooth boundary conditions.
This is particularly the case when calculating *the Greeks*.
A modified approach, using the
explicit method
for the first few time steps, then reverting to the Crank-Nicolson method for the remaining steps
often overcomes this problem.

### Pricing American Style Options

When pricing options that include the possibility of early exercise special care must
be taken when solving Equation 3 for
*F _{i}*.
Taking the inverse of

*C*to calculate a value for

*F*then comparing the calculated values to the intrinsic value of the option and taking the larger value (on an element by element basis) results in incorrect option values. This is because the modified

_{i}*F*will no longer satisfy Equation 3 as correct values must do.

_{i}However an iterative approach to calculating the inverse (such as the Guass-Seidel matrix inversion method) may be used successfully.