Crank-Nicolson Finite Difference Method (FDM)

OK, this is the last FDM method we’re covering. It sits somewhere in the middle between the Explicit and Implicit methods we did before, and it can improve time accuracy without the small-step restriction of Explicit FDM. That does not make it the fastest or best choice for every grid, though.

Anyway. Let me throw the Crank-Nicolson FDM method at you first. Take the Black-Scholes PDE:

For each term in this partial differential equation that’s a derivative, we put in approximations like these.

For $t$, we use the difference between the two time levels. Think of it at their midpoint. For the first and second derivatives with respect to $S$, we’re going to put in something a little weird.

For $\partial f / \partial S$: at $i$ and at $i+1$ we each take the central difference centered at $j$, and put in the average of the two.

Same thing for $\partial^2 f / \partial S^2$: at $i$ and at $i+1$ we each approximate the second derivative with central differences centered at $j$, and put in the average of the two.

And one more thing: the $-rf$ term has to be averaged too, as $-r(f_{i,j}+f_{i+1,j})/2$. The original notes leave it all at time $i$; the code block below corrects that.

Ahahahaha — to somehow convince myself that approximating it this way is actually “reasonable”, I thought about it like this.

Say, on some arbitrary grid —

Thinking about the midpoint made this click for me! With all the terms averaged, Crank-Nicolson is second order in time for smooth solutions, while the earlier time-stepping methods are first order. A coarse grid still needs care, especially at the payoff kink.

Alright. Now let’s plug all of that into the Black-Scholes PDE, the way we just saw, and tidy up the equation.

Breaking the equation above into tiny pieces and carrying the algebra forward in as much detail as I can —

— we can bring it all the way down to here and call it done.

This was less math and more, like, arithmetic, hmm…

From here, the way it gets computed is the same as Implicit FDM. Like we did with Implicit, let’s say we built a grid with 5 points each way.

With 5 points each way, our boundary conditions now look like this!!!! In our Black-Scholes approximation, plugging in $i=3$, $j=2$ gives an equation in those 3 unknowns.

OK, OK, OK. But with just this one equation, we’ve got 1 equation and 3 unknowns — nowhere near enough to pin down each of the 3 unknowns. So let’s also write out the equations for $j=1$ and $j=3$, and line all three up one after another.

LOL. The way this looks is genuinely gross, lmao. It only looks insanely complicated;;;

So now all we do is define $a_{1,2,3}$, $b_{1,2,3}$, $c_{1,2,3}$, $d_{1,2,3}$ and the boundary conditions and feed them to our cavalry, the Solve_SE() function procedure.

And out come $f_{3,3}$ / $f_{3,2}$ / $f_{3,1}$, right?!?!?

Right now!!!! Let’s study how exactly this situation gets coded~

It’s this!!!! And here’s how I studied it!!!!

Python version!

import numpy as np
import math

def Solve_SE(matA, matb):
    rowA = len(matA)
    colA = len(matA.T)
    rowB = len(matb)
    if (rowA != colA) | (rowB != colA):
        raise NotImplementedError
    else:
        pass
    L = np.zeros((rowA, colA))
    U = np.zeros((rowA, colA))
    for j in range(0, len(L.T)):
        for i in range(0, len(L)):
            if j==0:
                L[i, j] = matA[i, j]
                U[j, i] = matA[j, i] / L[0, 0]
            else:
                if i < j:
                    L[i, j] = 0
                else:
                    Lsum = 0
                    for k in range(0, j):
                        Lsum = Lsum + L[i, k]*U[k, j]
                    L[i, j] = matA[i, j] - Lsum
                if i < j:
                    U[j, i] = 0
                elif i == j:
                    U[j, i] = 1
                else:
                    Usum = 0
                    for k in range(0, i):
                        Usum = Usum + L[j, k]*U[k, i]
                    U[j, i] = (matA[j, i] - Usum) / L[j, j]
    modB = np.zeros((3, len(matb.T)))
    modB[0] = matb[0] / L[0, 0]
    for i in range(0, len(modB)):
        sumofB = 0
        for k in range(0, i):
            sumofB = sumofB + L[i, k]*modB[k]
        modB[i] = (matb[i] - sumofB) / L[i, i]
    x = np.zeros((3, len(modB.T)))
    for no in range(len(x)-1, 0-1, -1):
        sum=0
        for i in range(0, len(x)):
            sum = sum + U[no, i] * x[i]
        x[no] = modB[no] - sum
    return x

k=50; r=0.05; t=0.5; v=0.3; smax=100
m=4; n=4
ds = smax / m
dt = t/n
f = np.zeros((n+1, m+1))
matA = np.zeros((3, 3))
matb = np.zeros((3, 1))
matx = np.zeros((3, 1))

for j in range(1, m):
    f[n, j] = max(j*ds - k, 0)
for i in range(0, n+1):
    f[i, 0] = max(0, 0)
for i in range(0, n+1):
    f[i, m] = max(m*ds - k*math.exp(-r*(n-i)*dt), 0)

j=3
a3 = (-dt/(4*ds)) * (r*(j*ds) + (v**2) * (j*ds)**2 / ds)
b3 = 1 + (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) + r*dt/2
c3 = -dt / (4*ds) * (-r*(j*ds) + (v**2 * (j*ds)**2)/ds)
d3 = 1 - (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) - r*dt/2

j=2
a2 = (-dt/(4*ds)) * (r*(j*ds) + (v**2) * (j*ds)**2 / ds)
b2 = 1 + (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) + r*dt/2
c2 = -dt / (4*ds) * (-r*(j*ds) + (v**2 * (j*ds)**2)/ds)
d2 = 1 - (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) - r*dt/2

j=1
a1 = (-dt/(4*ds)) * (r*(j*ds) + (v**2) * (j*ds)**2 / ds)
b1 = 1 + (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) + r*dt/2
c1 = -dt / (4*ds) * (-r*(j*ds) + (v**2 * (j*ds)**2)/ds)
d1 = 1 - (v**2) * ((j*ds)**2) * (dt/(2*ds**2)) - r*dt/2

matA[0, 0] = b3; matA[0, 1]=c3; matA[0, 2]=0
matA[1, 0] = a2; matA[1, 1]=b2; matA[1, 2]=c2
matA[2, 0] = 0;  matA[2, 1]=a1; matA[2, 2]=b1

matb[0] = -a3 * f[4, 4] + d3*f[4, 3] - c3*f[4, 2] - a3*f[3, 4]
matb[1] = -a2 * f[4, 3] + d2*f[4, 2] - c2*f[4, 1]
matb[2] = -a1 * f[4, 2] + d1*f[4, 1] - c1*f[4, 0] - c1*f[3, 0]

matx = Solve_SE(matA, matb)

for i in range(len(matx)):
    print(round(float(matx[i, 0]), 4))

Heads-up. The figures here are the original notes, slips and all. A few things to watch for:

• The VBA and Python screenshots use a mistaken leading minus sign in b. Changing it to plus gives 25.3182, 0.7031 and 0.0030, but there is one more correction for full Crank-Nicolson: average the discount term too. Then $b_j=1+(\sigma^2j^2+r)\Delta t/2$ and $d_j=1-(\sigma^2j^2+r)\Delta t/2$. The corrected code block gives 25.3192, 0.7053 and 0.0031 for the first backward step, at time 0.375, not at time zero.
• In the derivation pictures, the averaged second derivative takes one term at time $i$ and one at time $i+1$ (the picture repeats the time-$i$ term), and the last term of the long line needs $f_{i+1,j-1}$; the off-diagonal coefficients $a_j$ and $c_j$ are right, but $b_j$ and $d_j$ also need the discount-term correction above.
• In the 3 x 3 system, the $j = 2$ row uses $a_2$, $d_2$, $c_2$ (the picture writes $a_3$, $d_3$, $c_3$), and the last right-hand entry ends in $-c_1 f_{3,0}$ rather than $-a_1 f_{3,0}$; since $f_{3,0} = 0$, the numbers are unaffected.