Finite Difference Methods in Option Pricing
Finite-Difference Methods for Option Pricing¶
$$ \frac{\partial V}{\partial t} + rS\frac{\partial V}{\partial S} + \frac{1}{2}\sigma^2 S^2\frac{\partial^2 V}{\partial S^2} = rV $$
We will do two transformations
$$ \tau = T - t $$$$ X = ln (S) $$Derivatives of SDE
$$ \Theta = \frac{\partial V}{\partial t} $$$$ \Delta = \frac{\partial V}{\partial S} $$$$ \Gamma = \frac{\partial^2 V}{\partial S^2} $$$$ \Theta + rS\Delta + \frac{1}{2}\sigma^2 S^2\Gamma = rV $$FTCS¶
Theta¶
We are going to approximate $\Theta$ using a forward Euler method.
Using Taylor
$$ V(S, t - \delta t) = V(S, t) + ((t - \delta t) - t)\frac{\partial V}{\partial t}(S,t) + \ldots $$$$ V(S, t - \delta t) = V(S, t) -\delta t\frac{\partial V}{\partial t}(S,t) + \mathcal{O}(\delta t^2) $$Rearrange for Theta
$$ \delta t\frac{\partial V}{\partial t}(S,t) = V(S, t) - V(S, t - \delta t) + \mathcal{O}(\delta t^2) $$$$ \frac{\partial V}{\partial t}(S,t) = \frac{V(S, t) - V(S, t - \delta t)}{\delta t} + \mathcal{O}(\delta t) $$In our discrete grid we have an axis over $\tau = T - t$ so we have to take care choosing the correct indices:
$$ \frac{\partial V}{\partial t}(S,t) = \frac{V^{n}_{i} - V^{n+1}_{i}}{\delta t} + \mathcal{O}(\delta t) $$Delta¶
For Delta we use a Centered approximation. We have:
$$ V(S - \delta S, t) = V(S, t) - \delta S\frac{\partial V}{\partial S}(S,t) + \frac{1}{2}\delta S^2\frac{\partial^2 V}{\partial S^2}(S,t) + \mathcal{O}(\delta S^3) $$and
$$ V(S + \delta S, t) = V(S, t) + \delta S\frac{\partial V}{\partial S}(S,t) + \frac{1}{2}\delta S^2\frac{\partial^2 V}{\partial S^2}(S,t) + \mathcal{O}(\delta S^3) $$The odd numbered terms on the right side cancel leaving us with a with a $\mathcal{O}(\delta S^3)$ error term. We then divide by $2\delta S$. (Note that we can just move the error term to the other side since the upper bound notation absors multiplication by some negative constant).
$$ \frac{\partial V}{\partial S}(S,t) = \frac{V(S + \delta S, t) - V(S - \delta S, t)}{2\delta S} + \mathcal{O}(\delta S^2) $$Since the other axis in our grid runs in the direction of $S$ we can straightforwardly formulate our discretisation:
$$ \frac{\partial V}{\partial S}(S,t) = \frac{V^{n}_{i+1} - V^{n}_{i-1} }{2\delta S} + \mathcal{O}(\delta S^2) $$Gamma¶
The Gamma is approximated by taking the difference between the forward Delta and the backward Delta. We have already shown the Taylor expansions for both the forward and the backward in the previous section.
Forward $$ \frac{V(S + \delta S, t) - V(S, t)}{\delta S} = \frac{\delta S \frac{\partial V}{\partial S}(S,t)+ \frac{1}{2}\delta S^2\frac{\partial^2 V}{\partial S^2}(S,t) + \mathcal{O}(\delta S^3)}{\delta S} $$
Backward $$ \frac{V(S, t) - V(S - \delta S, t)}{\delta S} = \frac{\delta S \frac{\partial V}{\partial S}(S,t) - \frac{1}{2}\delta S^2\frac{\partial^2 V}{\partial S^2}(S,t) + \mathcal{O}(\delta S^3)}{\delta S} $$
Since the forward is sort of a central-difference at $S + \frac{1}{2}\delta S$ (and similarly the backward at $S - \frac{1}{2}\delta S$), we can now take the difference between these two (forward-backward)
$$ \frac{\partial^2 V}{\partial S^2}(S,t) = \frac{V(S + \delta S, t) - 2V(S, t) + V(S - \delta S, t)}{\delta S^2} + \mathcal{O}(\delta S^2) $$In our discrete grid this is
$$ \frac{\partial^2 V}{\partial S^2}(S,t) = \frac{V^{n}_{i+1} - 2V^{n}_{i} + V^{n}_{i-1} }{\delta S^2} + \mathcal{O}(\delta S^2) $$Approximation¶
we insert all the terms and lose the higher order error terms, accepting them as approximation errors
$$ \frac{V^{n}_{i} - V^{n+1}_{i}}{\delta t} + (r-\frac{1}{2}\sigma^2) \frac{V^{n}_{i+1} - V^{n}_{i-1} }{2\delta S} + \frac{1}{2}\sigma^2\frac{V^{n}_{i+1} - 2V^{n}_{i} + V^{n}_{i-1} }{\delta S^2} = rV^{n}_{i} $$rewriting for $V^{n+1}_i$
$$ V^{n}_{i} + \delta t(r-\frac{1}{2}\sigma^2) \frac{V^{n}_{i+1} - V^{n}_{i-1} }{2\delta S} + \delta t\frac{1}{2}\sigma^2\frac{V^{n}_{i+1} - 2V^{n}_{i} + V^{n}_{i-1} }{\delta S^2} - \delta trV^{n}_{i} = V^{n+1}_{i} $$As required by the assignment. The error of this approximation is $\mathcal{O}(\delta t, \delta S^2)$
Remember we use $S = i \delta S$
$$ A_i = \frac{1}{2}(\sigma^2 i^2 - ri)\delta t $$$$ B_i = 1 - (\sigma^2 i^2 + r)\delta t $$$$ C_i = \frac{1}{2}(\sigma^2 i^2 + ri)\delta t $$In order to aid readability, in our code we will refer to these matrix constants as:
In [1]:
def A_FTCS(i, volatility, r, dt):
return 0.5*(volatility * volatility * i*i - r*i)*dt
def B_FTCS(i, volatility, r, dt):
return 1 - (volatility * volatility * i*i + r)*dt
def C_FTCS(i, volatility, r, dt):
return 0.5*(volatility * volatility * i*i + r*i)*dt
Crank-Nicolson¶
The Crank-Nicolson method can be seen as a combination of the explicit and implicit method.
$$ -A^{n+1}_{i} V^{n+1}_{i-1} + (1 - B^{n+1}_{i}) V^{n+1}_{1} - C^{n+1}_{i}V^{n+1}_{i+1} = A^{n}_{i} V^{n}_{i-1} + (1 - B^{n}_{i}) V^{n}_{1} + C^{n}_{i}V^{n}_{i+1} $$Note that here we have a matrix inversion to solve (as opposed to the FTCS scheme where there were no constraints on other $V^{n+1}$ values).
$$ A_i = \frac{1}{4}(\sigma^2 i^2 - ri)\delta t $$$$ B_i = -\frac{1}{2} (\sigma^2 i^2 + r)\delta t $$$$ C_i = \frac{1}{4}(\sigma^2 i^2 + ri)\delta t $$In [2]:
def A_CN(i, volatility, r, dt):
return (dt/4.)*(volatility * volatility*i*i - r*i)
def B_CN(i, volatility, r, dt):
return -(dt/2.)*(volatility * volatility*i*i + r)
def C_CN(i, volatility, r, dt):
return (dt/4.)*(volatility * volatility*i*i + r*i)
We find in literature that the Crank-Nicolson scheme has a better error $\mathcal{O}(\delta t^2, \delta S^2)$
For the CN scheme, we use a L-U decomposition for the matrix inversion step. This is achieved using numpy's lu_solve. We can go from our tridiagonal matrix to $L, U$ matrices as follows:
Where $A\_, B\_, C\_$ are the values on the three diagonals in our matrix.
We can use bilinear interpolation on the resulting values (see source code) in order to find the value for a choice of parameters that is not exactly an existing grid point.
| S | Analytical | Explicit FTCS | Crank-Nicolson |
|---|---|---|---|
| 100 | 9.625357828843697 | 9.1839796120981561 | 9.221853409581783 |
| 110 | 15.128591111967928 | 14.593425290034666 | 14.604914128121614 |
| 120 | 21.788808338829327 | 21.299628855180373 | 21.263916199436547 |
[Table: 20 asset steps]
Boundary Conditions¶
Value¶
When $S=0$, the option value must be $0$. Therefore, $V[0, n]= 0$ for all $n$. When $S\to\infty$ (in the grid, when S goes near the maximum value) the diffusion term disappears in the limit and only $r$ influences the option price. Therefore we can set:
$E = max(S-K, 0)$.
$V[S_{MAX}, n] = E e^{(r(n \delta t))}$
Delta¶
As $S\to0$, the diffusion and drift go to zero as well. Therefore the boundary condition in our grid (where S is actually 0) is that the delta becomes zero here.
$$ V[0, k] = (V[0, k-1] (1 - r \delta t)) $$As $S\to\infty$, the delta goes to $1$. When we approach our maximum value $S_{MAX}$ we would like the dela to become close to 1. My takeaway from this exercise is that it is probably hard to enfore this when we have a small $S_{max}$ (see the last part of the assignment).
Gamma¶
For large values of $S$, e.g. at $S_{max}$ the payoff becomes almost linear in the underlying therefore the second derivative becomes 0. The rule we derive from this is:
$$ V[N_{max}, k] = 2 V[N_{{max}-1}, k] - V[N_{{max}-2}, k]) $$\problem{Check whether your computer program is inline with your expectations from the theoretical analysis of exercise A. Determine your own optimal mesh size.}
In the papers we find the constraint for stability:
\begin{equation} \delta t \leq \frac{\delta S^2}{2a} = \frac{\delta S^2}{\sigma^2 S^2}=\frac{1}{\sigma^2} \left( \frac{\delta S}{S} \right )^2 \end{equation}\label{c1} Let $N$ be the number of asset steps. If we choose $\delta t = 0.9 \frac{N^2}{\sigma^2}$ we are sure that equation \ref{c1} holds. We now experiment with increasing number of asset steps $N$ to see the convergence:
In [3]:
import numpy
import math
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib.colors
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import matplotlib.pyplot as pyplot
pyplot.style.use('ggplot')
%matplotlib inline
from pylab import rcParams
rcParams['figure.figsize'] = 15, 5
import sys
from scipy import linalg
In [4]:
def Payoff(S, K):
return max(S-K, 0.0)
In [5]:
def Visualise(NAS, T=1., S=100., K=110.0, volatility=0.1, r=-0.1):
sarr = numpy.zeros(NAS)
dS = 2*(S / NAS)
dt = 0.9 / (volatility * volatility) / (NAS * NAS)
NTS = int(1.0 / dt) + 1
dt = T / NTS
V = numpy.zeros((NAS, NTS))
for i in range(NAS):
sarr[i] = i*dS
V[i,0] = max(sarr[i] - K, 0.0)
for k in range(1, NTS):
V[0, k] = V[0, k - 1] * (1. - r * dt)
for i in range(1, NAS-1):
Delta = (V[i + 1, k - 1] - V[i - 1, k - 1]) / 2. / dS
Gamma = (V[i + 1, k - 1] - 2. * V[i, k - 1] + V[i - 1, k - 1]) / dS / dS
Theta = -0.5 * volatility * volatility * sarr[i] * sarr[i] * Gamma - r * sarr[i] * Delta + r * V[i, k - 1]
V[i, k] = max(sys.float_info.epsilon, V[i, k - 1] - dt * Theta)
V[NAS-1, k] = max(sys.float_info.epsilon, 2. * V[NAS - 2, k] - V[NAS - 3, k])
price = V[int(NAS/2), -1]
if NAS & 0x1:
price = (V[int(NAS/2), -1] + V[int(NAS/2 + 1), -1]) / 2
figure = pyplot.figure()
ax = figure.gca(projection='3d')
X = numpy.arange(0.0, 1.0, 1.0 / NAS)
Y = numpy.arange(0.0, 1.0, 1.0 / NTS)
X,Y=numpy.meshgrid(Y,X)
surf = ax.plot_surface( X, Y, V
, rstride=1
, cstride=1
, cmap=cm.coolwarm
#, norm=matplotlib.colors.LogNorm(vmin=max(0.0001, V.min()), vmax=V.max() * 1)
, linewidth=0, )
In [6]:
Visualise(64, 1.)
In [7]:
def MakeGrids(NAS, NTS0=None, S=100., K=110.0, volatility=0.30, r=0.04):
sarr = numpy.zeros(NAS)
dS = ((2*S)/NAS)
dt = 0.9/(volatility*volatility)/(NAS*NAS)
NTS = int(1.0 / dt) + 1 if NTS0 is None else NTS0
dt = 1.0 / NTS
V = numpy.zeros((NAS, NTS))
V2 = numpy.zeros((NAS, NTS))
V3 = numpy.zeros((NAS, NTS))
for i in range(NAS):
sarr[i] = (i*dS)
V[i,0] = max(sarr[i] - K, 0.0)
V2[i,0] = max(sarr[i] - K, 0.0)
V3[i,0] = max(sarr[i] - K, 0.0)
for n in range(1, NTS):
V2[0, n] = 0
V2[NAS-1, n] = (sarr[i] - K) * math.exp(r*(n*dt))
V3[0, n] = 0
V3[NAS-1, n] = (sarr[i] - K) * math.exp(r*(n*dt))
Mat = numpy.zeros((NAS-2, NAS-2))
Mat2 = numpy.zeros((NAS-2, NAS-2))
for i in range(1, NAS-1):
Mat[i-1, i-1] = 1 + B_CN(i, volatility, r, dt)
Mat2[i-1, i-1] = 1 - B_CN(i, volatility, r, dt)
if i < NAS - 2:
Mat[i, i-1] = A_CN(i, volatility, r, dt)
Mat2[i, i-1] = -A_CN(i, volatility, r, dt)
Mat[i-1, i] = C_CN(i, volatility, r, dt)
Mat2[i-1, i] = -C_CN(i, volatility, r, dt)
P, L, U = linalg.lu(Mat2)
for k in range(1, NTS):
operand = V[:, k-1]
V[0, k] = V[0, k - 1] * (1. - r * dt)
for i in range(1, NAS-1):
Delta = (V[i + 1, k - 1] - V[i - 1, k - 1]) / 2. / dS
Gamma = (V[i + 1, k - 1] - 2. * V[i, k - 1] + V[i - 1, k - 1]) / dS / dS
Theta = -0.5 * volatility * volatility * sarr[i] * sarr[i] * Gamma - r * sarr[i] * Delta + r * V[i, k - 1]
V[i, k] = V[i, k - 1] - dt * Theta
V[NAS-1, k] = 2. * V[NAS - 2, k] - V[NAS - 3, k]
first =(V[0, k]* (1. - r * dt))
last = (2. * V[NAS - 2, k] - V[NAS - 3, k])
vec6=Mat.dot(operand[1:-1])
vec6[0] =vec6[0] + A_CN(NAS-2, volatility, r, dt) * V[0, k-1]
vec6[-1] =vec6[-1] + C_CN(NAS-2, volatility, r, dt) * V[NAS-1, k-1]
V2[0, k] = first
V2[1:NAS-1,k] = vec6
V2[NAS-1, k] = last
vec7 = Mat.dot(operand[1:-1])
offset = numpy.zeros(len(vec7))
offset[0] = A_CN(1, volatility, r, dt)*(V3[0,k]+V3[0,k-1])
offset[-1] = offset[-1] + C_CN(NAS-2, volatility, r, dt)*(V3[-1,k]+V3[-1,k-1])
V3[0, k] = first
V3[1:NAS-1,k] = linalg.inv(Mat2).dot(vec7 + offset)
V3[NAS-1, k] = last
price = V2[int(NAS / 2), -1]
if NAS & 0x1:
price = (V2[int(NAS / 2), -1] + V2[int(NAS / 2) + 1, -1])/2
price3 = V3[int(NAS / 2), -1]
if NAS & 0x1:
price3 = (V3[int(NAS / 2), -1] + V3[int(NAS / 2) + 1, -1])/2
return price, price3, V
In [8]:
x = range(5,41,1)
y0 =[]
y1 =[]
y2 =[]
for i in x:
r = MakeGrids(i)
y0.append(9.625357828843697)
y1.append(r[0])
y2.append(r[1])
In [9]:
pyplot.plot(x, y0, label="Analytical result")
pyplot.plot(x, y1, label="Forward-Time Central-Space")
pyplot.plot(x, y2, label="Crank-Nicolson")
pyplot.xlabel("Asset steps")
pyplot.ylabel("Option price")
pyplot.legend(bbox_to_anchor=(1, 0.27))
Out[9]:
