Implicit Finite Difference Method (FDM) — Part 2
Walking through the Implicit FDM code for pricing a European call — from a 3×3 matrix to the fully generalized grid, step by step and hehe all the way!
We got a thousand-troop reinforcement in the previous post!!!!
No matter what Ax = b matrix equation is given to us,
we’ve got a function procedure that can crack them all wide open!!!!!
This is seriously full-on thousand-troop-reinforcement energy for real lolololololololololol
So now let’s price a European call with Implicit FDM!
For now, to keep the coding simple — that equation we derived earlier by splitting the grid into 5 pieces, the one that became a 3x3 matrix equation, you know that one~
Let’s try coding this up~
Wait no, not coding it — understanding the given code!
Let’s do that hehe
The code goes like this.
I did it in Python too! hehehe
My own studying of this code went like this.
This time, let’s study the code that generalizes so all f_i,j get derived
And for that, there’s one thing I want to nail down one more time
The Implicit FDM method
approximated the Black-Scholes PDE
with a central difference for the first derivative in S
and a forward difference for the first derivative in t
And when I worked it out
it came out like this,
and depending on how many pieces you chopped S into, the number of unknowns changed
Here the unknowns aren’t f_i+1, j but rather f_i,j+1, f_i,j, f_i,j-1
When we cut it into five pieces, these three were the unknowns
and we needed 3 equations to find those 3 unknowns… yeah
ahhhh haha so then if we cut it into 7, 5 unknowns show up?
And to find all 5 unknowns we’d need 5 equations?????
Let’s lightly look at how those five come about and move on
Just write it down and substitute, that should do it
The boundaries — j=0, j=7, i=i+1 will be the boundaries~~~ is what we’ll say (I won’t say exactly what i is)
If we write out the 5 equations that have boundary values
it’s these 5 equations,
and if we check just the unknowns here
these here are the unknowns, y’see??? So the unknown x-vector gets written with these, and we can write the matrix equation somethin’ like this~
We can write it somethin’ like this~…
Writing it like this, I kinda get the rule~
And through this process we find f_i,~~,
and through the next process we find f_i-1,~~,
and goiiiing back and back until i=0 is all done,
wouldn’t we end up finding all the f_i,j on the grid????
Ah…. and one thing to watch carefully — a_j, b_j, c_j needed in the matrix are not from 1 to 7 but from 2 to 6!!!! hehe
That is, not j=1, j=7 but from 2 to 6!!!!!!
From just above the boundary to just below the boundary is what’s needed~~~
Also, the shape of the matrix is somewhat worth paying attention to.
a, b, c get written in this kinda diagonal shape!!!~~
Alright then, let’s study the code that finds all f_i,j! hehehe
This one too — I did it in Python! hehehe
This right here is the code~!!!!
This is how I looked at this code!!!
From here we can also simply pull out the Greeks
We went over pulling out these Greeks in class too
Once you get the principle it’s ridiculously simple
so there’s no particular interpretation for this one lololol
Just “oh, so we can check the Greeks like this too~~~~”
something like that??????????????????
For detailed Greeks, here’s the link http://gdpresent.blog.me/220887482357
Black-Scholes formula practice problems [ Derivatives I studied #17 ]
So now we’re going to solve the problems in Chapter 13, and there’s one thing to point out before solving them.Black and …
gdpresent.blog.me
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[k, i]*U[j, k]
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]), 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[k, i]*U[j, k]
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]
elif ii == jj:
matA[ii, jj] = B[mm]
elif ii+1 == jj:
matA[ii, jj] = A[mm]
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, jj+1]
matx = Solve_SE(matA, matb)
jj=m-1
for ii in range(len(matx)):
f[i, jj]=matx[ii]
jj=jj-1
for j in range(m+1):
print("stock price: ", j*ds, "option price: ", round(f[0, j], 2))
Originally written in Korean on my Naver blog (2016-12). Translated to English for gdpark.blog.