Lecture 13: Numerical Methods for Partial Differential Equations
Topics
- More on the stability of FDM for Black-Scholes
- Multi spatial variables: ADI method
- Numerical examples
More on the stability of FDM for Black-Scholes:
- When r≠0:
ea△t={12σ2x2j△t△x2−rxj△t2△x}e−ilm△x+{1−σ2x2j△t△x2−r△t}+{12σ2x2j△t△x2+rxj△t2△x}eilm△x
- which is
ea△t=1−r△t+σ2x2j△t△x2⋅2sin2(lm△x/2)+irxj△t△xsin(lm△x)
- To satyisfiy |ea△t|<1, we require
r>0,△t<△x2σ2x2max,△t<σ2r2.
The first condition is trivial (until the FED decrees all rates to be negative!)
The last condition does not have much impact in practice unless the volatility is very small.
So the important one is the second condition, which is what we went through last lecture.
The Case of more than one spatial variables
- Most numerical methods, FDM included, suffer the "dimentionality" curse, i.e. the size, complexity of the problem grow exponentially with the dimension.
- Usually, Monte Carlo method is the only practical option in dimension ≥3.
- But for two dimensional problems, it is worthwhile to explore the FDM further.
- Finance examples that you will need many sptial variables: stochastic vol model, convertible bonds, credit risky bonds, variable annuities, etc.
Heston Stochastic Volatility Model
- Heston Stochastic Vol Model
dSt=rStdt+√νtStdW1tdνt=κ(θ−νt)+ξ√νtdW2t⟨dW1t,dW2t⟩=ρdt
PDE for Heston Stochastic Vol Model:
∂u∂t+rS∂u∂S+[κ(θ−ν)−λν]∂u∂ν+12νS2∂2u∂S2+ρξνS∂2u∂S∂ν+12ξ2ν∂2u∂ν2−ru=0
Alternating Direction Implicit Method
- Write the Crank-Nicolson method as
ui,k+1−ui,k△t=14σ2x2i{ui+1,k−2ui,k+ui−1,k△x2+ui+1,k+1−2ui,k+1+ui−1,k+1△x2}+12rxi{ui+1,k−ui−1,k2△x+ui+1,k+1−ui−1,k+12△x}−12r{ui,k+ui,k+1}=12D⋅(ui,k+ui,k+1)
where D⋅ui,k can be considered as the Crank-Nicolson finite difference operator on ui,k.
Numerical Examples
- Black Scholes formulae and FDM methods for vanilla European call and up-and-out Barrier Call.
In [242]:
import numpy as np
from scipy.stats import norm
import time
#Black and Scholes
def BlackScholesFormula(type, S0, K, r, sigma, T):
dtmp1 = np.log(S0 / K)
dtmp2 = 1.0/(sigma * np.sqrt(T))
sigsq = 0.5 * sigma * sigma
d1 = dtmp2 * (dtmp1 + T * (r + sigsq))
d2 = dtmp2 * (dtmp1 + T * (r - sigsq))
if type=="C":
return S0 * norm.cdf(d1) - K * np.exp(-r * T) * norm.cdf(d2)
else:
return K * np.exp(-r * T) * stats.cdf(-d2) - S0 * norm.cdf(-d1)
K = 100.0
r = 0.05
sigma = 0.35
T = 1
putCall ='C'
Smin = 0.0
Smax = 200.0
ns = 201
Ss = np.linspace(Smin, Smax, ns, endpoint=True)
Ss = Ss[1:-1]
#print Ss
t=time.time()
px = BlackScholesFormula(putCall, Ss, K, r, sigma, T)
elapsed=time.time()-t
print px
print "Elapsed Time:", elapsed
#idx = 200-1
#print Ss[idx]
#print px[idx]
payoff = clip(Ss-K, 0.0, 1e600)
#print "payoff = ", payoff
figure(figsize=[10, 4])
subplot(1, 2, 1)
plot(Ss, payoff)
xlabel('Stock', fontsize=14);
title('Payoff' , fontsize=16);
subplot(1, 2, 2)
plot(Ss, px)
xlabel('Stock', fontsize=14);
title('Analytical Solution' , fontsize=16);
In [243]:
from scipy import sparse
import scipy.sparse.linalg.dsolve as linsolve
class BS_FDM_explicit:
def __init__(self,
r,
sigma,
maturity,
Smin,
Smax,
Fl,
Fu,
payoff,
nt,
ns):
self.r = r
self.sigma = sigma
self.maturity = maturity
self.Smin = Smin
self.Smax = Smax
self.Fl = Fl
self.Fu = Fu
self.nt = nt
self.ns = ns
self.dt = float(maturity)/nt
self.dx = float(Smax-Smin)/(ns+1)
self.xs = Smin/self.dx
self.u = empty((nt + 1, ns))
self.u[0,:] = payoff
## Building Coefficient matrix:
A = sparse.lil_matrix((self.ns, self.ns))
for j in xrange(0, self.ns):
xd = j + 1 + self.xs
sx = self.sigma * xd
sxsq = sx * sx
dtmp1 = self.dt * sxsq
dtmp2 = self.dt * self.r
A[j,j] = 1.0 - dtmp1 - dtmp2
dtmp1 = 0.5 * dtmp1
dtmp2 = 0.5 * dtmp2 * xd
if j > 0:
A[j,j-1] = dtmp1 - dtmp2
if j < self.ns - 1:
A[j,j+1] = dtmp1 + dtmp2
self.A = A.tocsr()
### Building bc_coef:
nxl = 1 + self.xs
sxl = self.sigma * nxl
nxu = self.ns + self.xs
sxu = self.sigma * nxu
self.blcoef = 0.5 * self.dt * (sxl * sxl - self.r * nxl)
self.bucoef = 0.5 * self.dt * (sxu * sxu + self.r * nxu)
def solve(self):
for i in xrange(0, m):
self.u[i+1,:] = self.A * self.u[i,:]
self.u[i+1,0] += self.blcoef * self.Fl[i]
self.u[i+1,self.ns-1] += self.bucoef * self.Fu[i]
return self.u
dx = (Smax - Smin)/(ns-1)
print "ns =", ns
print "dx =", dx
print "sigma =", sigma
dt_max = dx*dx/(sigma*sigma*Smax*Smax)
print "by CFL, dt <", dt_max
mt_min = int(T/dt_max)+1
print "mt_min ~= ", mt_min
Homework
- Implement the FDM for Black-Scholes PDE using the Crank-Nicolson Scheme and test your code using the same example here for both the vanilla European Call and the Up-and-out Barrier Call.