Implicit Finite Difference Method (FDM) — Part 2 [Financial Engineering Programming #20]

Last time we picked up some serious firepower!!!! For the systems here — with nonzero pivots — we’ve got a function that solves Ax = b!!!!! Like, we are LOADED. Cavalry’s here lol

So now let’s actually price a European call with Implicit FDM! For the code right now, let’s keep it simple: it’s the 3x3 matrix equation we derived back when we split the grid into 5 points — yeah, that one:

The 3 × 3 matrix equation from last time, for the three unknowns f3,3, f3,2, f3,1

Let’s code this thing up. Wait — not “code it ourselves” exactly. More like, read and understand the code that’s given. Let’s do that, heh. Here’s the code:

VBA program ImplicitFDMEuropean: m = n = 4, the coefficients a, b, c for j = 3, 2, 1, matA and matb, Solve_SE, and the three option prices in the Immediate window

I tried it in Python too! heh heh heh

The same program in Python, first part: parameters, boundary values and the coefficients for j = 3, 2, 1

The Python program, second part: matA, matb, Solve_SE and the printed option prices 25.3248, 0.6927, 0.006

And here’s how I worked through it on my own:

Handwritten notes on the VBA program, top half: the declarations and the boundary-condition loops, the same as for Explicit

Handwritten notes, bottom half: the j = 2 and j = 1 coefficients are copy-paste, the matrix equation from before, and Solve_SE, the firepower we picked up

OK, now let’s study the code that generalizes — the one that gives us all the $f_{i,j}$ values. But before that, there’s one thing to go over one more time.

Quick recap. Implicit FDM was: approximate the Black-Scholes PDE with central difference for the first-order derivative in $S$, and forward difference for the first-order derivative in $t$. And when we did that:

this is what fell out — and the number of unknowns depended on how many pieces we cut $S$ into. The unknowns here aren’t $f_{i+1, j}$ — they’re $f_{i, j+1}$, $f_{i, j}$, $f_{i, j-1}$.

With 5 grid points, those 3 things were the unknowns. 3 unknowns → we needed 3 equations…

Ahahaha, OK, so with 7 grid points we get 5 unknowns, right? And we’d need 5 equations to pin all 5 of them down??

Let’s quickly see where those 5 come from and move on.

Just write it out and plug stuff in. The boundaries will be $j=1$, $j=7$ and the time step $i+1$~~~ let’s just say that (I’m not going to be precise about $i$). Writing out the 5 equations that pull in boundary values:

The relation written out for j = 6, 5, 4, 3, 2

So we’ve got these 5 equations. And if we pick out just the unknowns from these:

The same five equations with the unknowns f(i,·) in red

These are the unknowns, right??? So the unknown $x$-vector is that, and we can write the matrix equation like:

The 5 × 5 matrix equation for the unknowns f(i,6) … f(i,2), with the boundary terms a6 f(i,7) and c2 f(i,1) moved to the right-hand side

Yeah, something like that… Now that it’s written out, I can kind of see the pattern.

And through this process we find $f_{i, \,\sim}$, then through the next pass we find $f_{i-1, \,\sim}$, and if we keep stepping backward all the way down to $i=0$… boom, we’ve filled out every $f_{i,j}$ on the grid, right??

Oh — and one thing to watch carefully: the $a_j$, $b_j$, $c_j$ we need in the matrix aren’t $j=1$ through $7$, they’re $j=2$ through $6$!!!! Not $j=1$ through $j=7$ — from 2 through 6!!!!!! We need everything from just above the boundary to just below the boundary~~~

Also, the shape of the matrix is worth a careful look.

The matrix: the b’s on the diagonal, the c’s just above it, the a’s just below it

It’s like $a$, $b$, $c$ laid out along the diagonals like this!!!~~

Before the code, a naming heads-up: the seven-point schematic labels the $f_{i,j+1}$ coefficient $a_j$ and the $f_{i,j-1}$ coefficient $c_j$. In the code below those names are swapped: $A$ multiplies $f_{i,j-1}$ and $C$ multiplies $f_{i,j+1}$. The matrix order follows that code convention.

OK then — let’s study the code that finds all the $f_{i,j}$! heh heh heh

VBA program ImplicitFDMEuropeancall2 for the whole grid, m = n = 10: boundary values, the coefficient arrays A, B, C, the loops that build matA and matb

The end of the VBA program: Solve_SE, writing the solution back into f, and the table of stock price and option price in the Immediate window

Did this one in Python too! heh heh heh

The general program in Python, first part

The general program in Python, second part, with the printed table

This is exactly the code!!! And here’s how I read it:

Handwritten notes on the general program, top: the declarations, why matdim = m − 2, and the boundary values

Handwritten notes, middle: the coefficient loop, how the loops build matrix A row by row while mm counts down from matdim, and how b gets its top, bottom and middle entries

Handwritten notes, end: Solve_SE gives x, x is written back into f (jj and ii run in opposite directions), and the print loop shows the present option values f(0,j)

We can also pull out the Greeks pretty easily here, and yeah, we covered the Greeks-extraction part in class too.

Once you know the principle, it’s a piece of cake — there’s no deep interpretation or anything lol More like — “oh nice, we can also check the Greeks this way~~~~” — that level??

For the Greeks in detail, see The Greeks (Derivatives I Studied #20).

Handwritten VBA program for the Greeks, first part: the declarations, including delta, gamma and theta, and the same grid setup

The Greeks program, second part: the same loops that build matA and matb and solve each time step

The Greeks from the grid: delta from neighbouring f’s in S over ds, gamma from three neighbours over ds², theta from neighbouring f’s in time over dt

Python code

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))

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 = 0.5 * dt * (r*j - (v**2)*(j**2))
b3 = 1 + ((v*j)**2)*dt + r*dt
c3 = 0.5 * dt * (-r*j - (v**2)*(j**2))

j=2
a2 = 0.5 * dt * (r*j - (v**2)*(j**2))
b2 = 1 + ((v*j)**2)*dt + r*dt
c2 = 0.5 * dt * (-r*j - (v**2)*(j**2))

j=1
a1 = 0.5 * dt * (r*j - (v**2)*(j**2))
b1 = 1 + ((v*j)**2)*dt + r*dt
c1 = 0.5 * dt * (-r*j - (v**2)*(j**2))

matA = np.zeros((3, 3))
matA[0, 0] = b3; matA[0, 1] = a3; matA[0, 2] = 0
matA[1, 0] = c2; matA[1, 1] = b2; matA[1, 2] = a2
matA[2, 0] = 0; matA[2, 1] = c1; matA[2, 2] = b1

matb = np.zeros((3, 1))
matb[0] = f[4, 3] - c3*f[3, 4]
matb[1] = f[4, 2]
matb[2] = f[4, 1] - a1*f[3, 0]

matx = Solve_SE(matA, matb)

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

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((len(matb), 1))
    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((len(modB), 1))
    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=10; n=10

ds = smax / m
dt = t/n

matdim=m-2

f = np.zeros((n+1, m+1))

matA = np.zeros((matdim+1, matdim+1))
matb = np.zeros((matdim+1, 1))

A = np.zeros((matdim+1, 1))
B = np.zeros((matdim+1, 1))
C = np.zeros((matdim+1, 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)

for j in range(m-1, 1-1, -1):
    A[j-1] = 0.5 * dt * (r * j - (v ** 2) * (j ** 2))
    B[j-1] = 1 + ((v * j) ** 2) * dt + r * dt
    C[j-1] = 0.5 * dt * (-r * j - (v ** 2) * (j ** 2))

for i in range(n-1, 0-1, -1):
    mm = matdim
    for ii in range(0, matdim+1):
        for jj in range(0, matdim+1):
            if ii == jj + 1:
                matA[ii, jj] = C[mm, 0]
            elif ii == jj:
                matA[ii, jj] = B[mm, 0]
            elif ii+1 == jj:
                matA[ii, jj] = A[mm, 0]
            else:
                matA[ii, jj] = 0
        mm=mm-1

    matb[0] = f[i+1, m-1] - C[matdim]*f[i,m]
    matb[matdim] = f[i+1, 1] - A[0]*f[i,0]

    for jj in range(matdim-1, 1-1, -1):
        matb[jj] = f[i+1, m-1-jj]

    matx = Solve_SE(matA, matb)

    jj=m-1
    for ii in range(len(matx)):
        f[i, jj]=matx[ii, 0]
        jj=jj-1

for j in range(m+1):
    print("stock price: ", j*ds, "option price: ", round(f[0, j], 2))

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

In the general VBA program (typed and handwritten), the right-hand side fills row jj with the wrong grid row: it should be matb(jj) = f(i + 1, m - 1 - jj), because row 0 of the matrix is the top grid row (the handwritten note’s numbering of b is reversed for the same reason). The code in the post is corrected. With it the prices for S = 0, 10, …, 100 are 0.00, 0.00, 0.00, 0.07, 0.79, 4.18, 12.04, 21.45, 31.30, 41.25, 51.23, not the 0.65, 13.98, 4.19, … printed in the screenshots; For comparison, Black-Scholes gives 0.88 at 40, 4.82 at 50 and 12.23 at 60. This coarse grid still has a noticeable error, especially near the strike.
In the 3 x 3 example (typed and handwritten), the last entry on the right is f4,1 - a1 f3,0, as the code has it, not f4,1 - c1 f3,0 (f3,0 is 0, so the numbers are unaffected); in the row equations for j = 3 and j = 2 the middle coefficients are b3 and b2, not b4 and b3.
On the matrix size: with m = 10 there are 11 grid points in S and 9 unknowns, so the system is (m - 1) x (m - 1), and matdim = m - 2 = 8 is its last index (arrays start at 0). The code is right.
In the five-equation schematic, the j = 2 middle coefficient is b2, not b3. The handwritten Greeks program inherits the same reversed middle-RHS indexing as the general VBA program: use f(i + 1, m - 1 - jj). The first three-price example only solves the first backward time step, t = 0.375, not time zero.

I tried it in Python too! heh heh heh

Did this one in Python too! heh heh heh

Python code

Comments