Lecture 12: Numerical Methods for Partial Differential Equations

Objectives

Introduction
- Second Order PDE Classification
- Initial and Boundary Conditions
- Black-Scholes Equation

Finite Difference Methods for Solving PDEs
- Discretization
- Finite Difference Methods
- The Explicit Methods
- The Implicit Methods
- The Crank-Nicholson Method
- Consistency, Stability and Convergence

Numerical Examples
- Vanilla European Call
- Up-And-Out Barrier Call

Introduction

In finance, most partial differential equations (PDE) we encounter are of the form

$$ \renewcommand{PDut}{\frac{\partial u}{\partial t}} \renewcommand{PDux}{\frac{\partial u}{\partial x}} \renewcommand{PDutt}{\frac{\partial ^2u}{\partial t^2}} \renewcommand{PDuxx}{\frac{\partial ^2u}{\partial x^2}} \renewcommand{PDuyy}{\frac{\partial ^2u}{\partial y^2}} \large{ a(x,t)\PDutt + 2b(x,t)\frac{\partial ^2u}{\partial t\partial x} + c(x,t)\PDuxx + d(x,t)\PDut + e(x,t)\PDux = f(t,x,u) } $$

Classified into: Hyperbolic, Parabolic and Elliptic

In general, the different types of equations make a big difference in how the solutions behave and

On how we can solve them more effectively

Classification

Hyperbolic

$b^2(x,t) - a(x,t)c(x,t) > 0 $

$\;$
The canonical form:

$$ \large{ \PDutt = c^2 \PDuxx } $$

Mostly referred to as the Wave equation

Parabolic

$b^2(x,t) - a(x,t)c(x,t) = 0 $

$\;$
The canonical form:

$$ \large{ \PDut = d \;\PDuxx } $$

Mostly referred to as the Heat equation

Elliptic

$b^2(x,y) - a(x,y)c(x,y) < 0 $

$\;$
The canonical form:

$$ \large{ \PDuxx + \PDuyy = 0 } $$

Mostly referred to as the Laplace or Poisson equation. The solutions are Harmonic functions, describing e.g. the gravitational, electric potentials.

Initial and Boundary Conditions

Depending on the problem, we usually have initial conditions (IC)

$$ \large{ x(x, t=0) = u_0(x) } $$

and/or boundary conditions (BC)

$$ \large{ a \; u + b\; \PDux = c, \;\forall \; x \in \partial Q, \;\forall \; t } $$

and the BC is called Dirichlet if $b = 0$, Neumann if $a = 0$, and Robin otherwise .

Black-Scholes Equation

$$ \renewcommand{PDuS}{\frac{\partial u}{\partial S}} \renewcommand{PDuSS}{\frac{\partial ^2u}{\partial S^2}} \PDut + \frac{1}{2}\sigma^2S^2\PDuSS + rS\PDuS - ru = 0 $$

It belongs to the parabolic type, and can be converted into the standard Heat equation form by a change-of-variables
For European Call option, the "initial" (the payoff function) and boundary conditions are

\begin{aligned} &u(S,T) = \max(S-K, 0); \\ &u(0,t) = 0; \\ &u(S,t) = S - Ke^{-r(T-t)}, \;\mbox{as } \; S \rightarrow \infty \end{aligned}

For European Put option, the "initial" (the payoff function) and boundary conditions are

\begin{aligned} \renewcommand{FDut}{\frac{u_{i,k+1}-u_{i,k}}{\triangle t}} \renewcommand{FDutb}{\frac{u_{i,k}-u_{i,k-1}}{\triangle t}} \renewcommand{FDutc}{\frac{u_{i,k+1}-u_{i,k-1}}{2\triangle t}} \renewcommand{FDutt}{\frac{u_{i,k+1}-2u_{i,k}+u_{i,k-1}}{\triangle t^2}} \renewcommand{FDux}{\frac{u_{i+1,k}-u_{i,k}}{\triangle x}} \renewcommand{FDuxb}{\frac{u_{i,k}-u_{i-1,k}}{\triangle x}} \renewcommand{FDuxc}{\frac{u_{i+1,k}-u_{i-1,k}}{2\triangle x}} \renewcommand{FDuxx}{\frac{u_{i+1,k}-2u_{i,k}+u_{i-1,k}}{\triangle x^2}} &u(S,T) = \max(K-S, 0); \\ &u(0,t) = Ke^{-r(T-t)}; \\ &u(S,t) = 0, \;\mbox{as } \; S \rightarrow \infty \end{aligned}

Finite Difference Methods for Solving PDEs

Very rare that the PDEs have closed form solutions, in general, can only be solved numerically.

Will focus on the finite difference method (FDM): other numerical methods exist and can be more appropriate---very much depends on the actual problems at hand.

Discretization

First step in solving PDEs using FDM is to represent the solution $u(x,t)$ as a discrete collection of values at a well distributed grid points in space and time in the proper domain.

For example, for the rectangular domain in space and time with $t \in [0, T]$ and $x\in [x_{min}, x_{max}]$, we can represent the points simply as $t_0 = 0, t_1, t_2, \cdots, t_M = T$ in time and $x_0 = x_{min}, x_1, x_2, \cdots, x_N = x_{max}$ in space. The grid is uniform, i.e. $t_{k+1} = t_k + \triangle t, \;\triangle t = \frac{T}{M}$ and $x_{i+1} = x_i + \triangle x, \;\triangle x = \frac{x_{max}-x_{min}}{N}$.

The discretized version of the solution will be represented by the values of

$$ u_{i,k} = u(x_i, t_k), \;\mbox{for } i = 0,1,\cdots,N, \mbox{ and } k = 0,1,\cdots,M $$

The values of $u(x,t)$ at any other desired points can be approximated via interpolation.

A uniform mesh may not necessary be the most efficient form to work with, in fact, it rarely is.

A greatly simplified rule-of-thumb: the mesh needs to be refined around the region where the function varies a lot

On the other hand, the mesh can be relatively coarse if the function is smooth and changing slowly

How do I know where to put more points and where to put less before knowing the solution: concepts of adaptive mesh

In finance, the time grid is well advised to keep the important dates as grid points: such as cash flow dates, contractual schedule dates, etc. The spacing of these dates is most likely not uniform.

Finite Difference Methods(FDM)

Toolkit for finite difference approximation

Partial Derivative	Finite Difference	Type	Order
$\PDux$	$\FDux$	forward	1st in $x$
$\PDux$	$\FDuxb$	backward	1st in $x$
$\PDux$	$\FDuxc$	central	2nd in $x$
$\PDuxx$	$\FDuxx$	symmetric	2nd in $x$
$\PDut$	$\FDut$	forward	1st in $t$
$\PDut$	$\FDutb$	backward	1st in $t$
$\PDut$	$\FDutc$	central	2nd in $t$
$\PDutt$	$\FDutt$	symmetric	2nd in $t$

The Explicit Method

For convenience, reverse the time for Black-Scholes notation: $ t \leftarrow T - t$

$$ \PDut - \frac{1}{2}\sigma^2S^2\PDuSS - rS\PDuS + ru = 0 $$

Approximate the Black-Scholes Eqn by

$$ \FDut - \frac{1}{2} \sigma^2 x_i^2 \FDuxx - r x_i \FDuxc + r u_{i,k} = 0 $$

This can be rewritten simply as

\begin{aligned} u_{i,k+1} & = \left\{ \frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} - r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i-1,k} \\ & + \left\{ 1 - \sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} u_{i,k} \\ & + \left\{ \frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i+1,k} \end{aligned}

The Implicit Method

Approximate the Black-Scholes Eqn by

$$ \FDutb - \frac{1}{2} \sigma^2 x_i^2 \FDuxx - r x_i \FDuxc = - r u_{i,k} $$

Which can be rewritten as

\begin{aligned} u_{i,k-1} & = \left\{ -\frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i-1,k} \\ & + \left\{ 1 + \sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} + r \triangle t \right\} u_{i,k} \\ & + \left\{ -\frac{1}{2}\sigma^2 x_i^2 \frac{\triangle t}{\triangle x^2} - r x_i \frac{\triangle t}{2\triangle x} \right\} u_{i+1,k} \end{aligned}

The Crank-Nicholson Method

Crank-Nicholson is an average of the explicit and implicit scheme:

\begin{aligned} \FDut & = \frac{1}{4} \sigma^2 x_i^2 \left\{ \FDuxx + \frac{u_{i+1,k+1}-2u_{i,k+1}+u_{i-1,k+1}}{\triangle x^2} \right\} \\ & + \frac{1}{2} r x_i \left\{ \FDuxc + \frac{u_{i+1,k+1}-u_{i-1,k+1}}{2\triangle x} \right\} \\ & - \frac{1}{2} r \left\{ u_{i,k} + u_{i,k+1} \right\} \end{aligned}

Consistency, Stability and Convergence

For any FDM that are used to solve practical problems, we should ask
- 1) How acurate is the method?
- 2) Does it converge?
- 3) What is the best choice of step sizes?

Consistency

For the explicit method above, the local truncation error is

$$ T(x,t) = O(\triangle x^2 + \triangle t) $$

Taylor series expansion can verify this (or one can simply read this off of the Toolkit Table).

So the explicit method above is consistent of order 1 in time and order 2 in spatial variable.

Fourier (Von Neumann) Stability Analysis

The Fourier method of analyzing the stability of FDM boils down to substitute $u_{j,k}$ with

$$ \epsilon_{j,k} = e^{a t_k}e^{il_mx_j} $$

The logic behind this: assuming $v_{j,k}$ is the solution of the FDM that did not suffer from finite precision arithmetic and is without the round off error, then the error

$$ \epsilon_{j,k} = v_{j,k} - u_{j,k} $$

$\hspace{0.2in}$ again satifies the FDM due to linearity.

Looking for $\epsilon_{j,k} $ in the Fourier expansion representation,

$$ \epsilon(x,t) = \sum_{m= 1}^{M} e^{at}e^{il_mx} $$

where $e^{at}$ is a special form of the amplitude, and $l_m$ is the wavelength: $l_m = \frac{\pi m}{L}, m = 1, \cdots, M, \mbox{ and } M = \frac{L}{\triangle x}$

Due to linearity again, the error analysis can be performed (without any loss of generality) by looking at a single term of the expansion

$$ e^{at_k}e^{il_mx_j} $$

In fact, we are really interested in the amplification of the error,

$$ \frac{\epsilon_{j,k+1}}{\epsilon_{j,k}} = e^{a\triangle t} $$

The stability analysis in the Fourier method is simply to analyze the form of $e^{a\triangle t}$ so that

$$ ||e^{a\triangle t}|| < 1 $$

For the explicit method,

\begin{aligned} e^{a (t_k + \triangle t)}e^{il_mx_j} & = \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} - r x_j \frac{\triangle t}{2\triangle x} \right\} e^{a t_k}e^{il_mx_{j-1}} \\ & + \left\{ 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} e^{a t_k}e^{il_mx_j} \\ & + \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} + r x_j \frac{\triangle t}{2\triangle x} \right\} e^{a t_k}e^{il_mx_{j+1}} \end{aligned}

Simplify it into,

\begin{aligned} e^{a\triangle t} & = \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} - r x_j \frac{\triangle t}{2\triangle x} \right\} e^{-il_m\triangle x} \\ & + \left\{ 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} -r \triangle t \right\} \\ & + \left\{ \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} + r x_j \frac{\triangle t}{2\triangle x} \right\} e^{il_m\triangle x} \end{aligned}

Looking at the simple case that $r = 0$,

\begin{aligned} e^{a\triangle t} & = \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} e^{-il_m\triangle x} \\ & + 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \\ & + \frac{1}{2}\sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} e^{il_m\triangle x} \end{aligned}

which is

\begin{aligned} e^{a\triangle t} &= \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cos(l_m\triangle x) + 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \\ & = 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2) \end{aligned}

The stability condition will be saisfied if

$$ || 1 - \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2) || < 1 $$

Or if

$$ \triangle t < \frac{\triangle x^2}{ \sigma^2 x_{max}^2 } $$

Which is the CFL condition for the explicit method. This says for the explicit method to be stable, the time step needs to be pretty small.
$r\neq 0$ case?

Lax Theorem: Consistency + Stability = Convergence

The Lax equivalence theorem holds for PDE as well.

Takeaway: when implementing FDM for PDEs, study the CFL (Courant–Friedrichs–Lewy) stability condition before building your discretization grid.

For the Implicit method, the local truncation error

$$ T(x,t) = O(\triangle x^2 + \triangle t) $$

And the amplification factor for $r=0$

$$ e^{a\triangle t} = \frac{1}{1 + \sigma^2 x_j^2 \frac{\triangle t}{\triangle x^2} \cdot 2\sin^2(l_m\triangle x/2)} $$

So the implicit method is unconditionally stable.

For the Crank-Nicholson method, the local truncation error is (Note: expanding the Taylor series at $(i, k+\frac{1}{2})$)

$$ T(x,t) = O(\triangle x^2 + \triangle t^2) $$

The stability anaylsis will be one of your homeworks.

Numerical Examples

Black Scholes formulae and FDM methods for vanilla European call and up-and-out Barrier Call.

In [101]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
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 = np.clip(Ss-K, 0.0, 1e600)
#print "payoff = ", payoff

fig = plt.figure(figsize=(15,4))
ax = fig.add_subplot(1,3,1)
ax.plot(Ss, payoff)
plt.xlabel('Stock', fontsize=14);
plt.title('Payoff' , fontsize=16);

ax = fig.add_subplot(1,3,2)
ax.plot(Ss, px)
plt.xlabel('Stock', fontsize=14);
plt.title('Analytical Solution' , fontsize=16);

ax = fig.add_subplot(1,3,3)
ax.set_axis_off()
K = 100.0
r = 0.05
sigma = 0.35
T = 1
putCall ='C'
ax.text(0, .2, "   Interest Rate = 0.05 ", size="x-large");
ax.text(0, .3, "   Strike = 100.0 ", size="x-large");
ax.text(0, .4, "   $\sigma$ = 0.35 ", size="x-large");
ax.text(0, .5, "   T = 1.0 ", size="x-large");
ax.text(0, .6, "   putCall =  'C' ", size="x-large");
ax.text(0, .75, "Black Scholes Equation with:", size="xx-large");

plt.show()

Elapsed Time: 0.00399994850159

In [72]:

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 = np.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

print "ESTIMATE THE TIME STEPS THROUGH CFL Condition: \n"
dx = (Smax - Smin)/(ns-1)
print "\tns(number of spatial steps) =", ns
print "\tdx (spatial size) =", dx
print "\tvolatility =", sigma
print "\tSmax =", Smax
dt_max = dx*dx/(sigma*sigma*Smax*Smax)
print "\tby CFL, dt <", dt_max
mt_min = int(T/dt_max)+1
print "\tTime Steps ~= ", mt_min

ESTIMATE THE TIME STEPS THROUGH CFL Condition: 

	ns(number of spatial steps) = 201
	dx (spatial size) = 0.6
	volatility = 0.35
	Smax = 120
	by CFL, dt < 0.000204081632653
	Time Steps ~=  4901

In [97]:

print "SOLVING for the Vanilla Call Using the Explicit FDM method: "
n = ns-2
X = np.linspace(0.0, Smax, n+2)
X = X[1:-1]

payoff = np.clip(X-K, 0.0, 1e600)
#print "payoff = ", payoff
  
m = 4555 
Fl = np.zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))
    
t = time.time()
bs1 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs1.solve()
elapsed=time.time()-t
print "\tElapsed Time1:", elapsed

##print  px_fd_mat.shape
nrow = len(px_fd_mat[:,1])
#print px_fd_mat[nrow-1,:]

fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(1,3,1)
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 4560
Fl = np.zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))

t = time.time()
bs2 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs2.solve()
elapsed=time.time()-t
print "\tElapsed Time2:", elapsed

ax = fig.add_subplot(1, 3, 2)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 4565
Fl = np.zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))
    
t = time.time()
bs3 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs3.solve()
elapsed=time.time()-t
print "\tElapsed Time3:", elapsed

ax = fig.add_subplot(1, 3, 3)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

plt.show()

SOLVING for the Vanilla Call Using the Explicit FDM method: 
	Elapsed Time1: 0.0469999313354
	Elapsed Time2: 0.0460000038147
	Elapsed Time3: 0.0439999103546

In [98]:

print "SOLVING for the UP-AND-OUT Barrier Call Using the Explicit FDM method: "
Smax = 120
n = ns-2
X = np.linspace(0.0, Smax, n+2)
X = X[1:-1]

payoff = np.clip(X-K, 0.0, 1e600)
#print "payoff = ", payoff
  
m = 4570
Fl = np.zeros((m+1,))
Fu = np.zeros((m+1,))
    
t = time.time()
bs1 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs1.solve()
elapsed=time.time()-t
print "\tElapsed Time1:", elapsed

##print  px_fd_mat.shape
nrow = len(px_fd_mat[:,1])
#print px_fd_mat[nrow-1,:]


fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(1,3,1)
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 4573
Fl = np.zeros((m+1,))
Fu = np.zeros((m+1,))

t = time.time()
bs2 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs2.solve()
elapsed=time.time()-t
print "\tElapsed Time2:", elapsed

ax = fig.add_subplot(1, 3, 2)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 4580
Fl = np.zeros((m+1,))
Fu = np.zeros((m+1,))
    
t = time.time()
bs3 = BS_FDM_explicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs3.solve()
elapsed=time.time()-t
print "\tElapsed Time3:", elapsed

ax = fig.add_subplot(1, 3, 3)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

SOLVING for the UP-AND-OUT Barrier Call Using the Explicit FDM method: 
	Elapsed Time1: 0.047000169754
	Elapsed Time2: 0.0480000972748
	Elapsed Time3: 0.0480000972748

In [99]:

print "SOLVING for the Vanilla Call Using the Implicit FDM method: "
Smax = 200
n = ns-2
X = np.linspace(0.0, Smax, n+2)
X = X[1:-1]

payoff = np.clip(X-K, 0.0, 1e600)
  
m = 4555 
Fl = np.zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))
    
t = time.time()
bs1 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs1.solve()
elapsed=time.time()-t
print "\tElapsed Time1:", elapsed

#print  px_fd_mat.shape

fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(1, 3, 1)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 1000
Fl =np. zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))

t = time.time()
bs2 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs2.solve()
elapsed=time.time()-t
print "\tElapsed Time2:", elapsed

ax = fig.add_subplot(1, 3, 2)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 100
Fl = np.zeros((m+1,))
Fu = Smax - K*np.exp(-r * np.linspace(0.0, T, m+1))
    
t = time.time()
bs3 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs3.solve()
elapsed=time.time()-t
print "\tElapsed Time3:", elapsed

ax = fig.add_subplot(1, 3, 3)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);


#print px_fd_mat[nrow-1,:]

SOLVING for the Vanilla Call Using the Implicit FDM method: 
	Elapsed Time1: 0.061999797821
	Elapsed Time2: 0.0260000228882
	Elapsed Time3: 0.00899982452393

Out[99]:

<matplotlib.text.Text at 0x27c476a0>

In [100]:

print "SOLVING for the UP-AND-OUT Barrier Call Using the Implicit FDM method: "
Smax = 120

n = ns-2
X = np.linspace(0.0, Smax, n+2)
X = X[1:-1]

payoff =np.clip(X-K, 0.0, 1e600)
  
m = 4555 
Fl = np.zeros((m+1,))
Fu = np.zeros((m+1,))
    
t = time.time()
bs1 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs1.solve()
elapsed=time.time()-t
print "\tElapsed Time1:", elapsed

#print  px_fd_mat.shape

fig = plt.figure(figsize=(16,5))
ax = fig.add_subplot(1, 3, 1)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 1000
Fl = np.zeros((m+1,))
Fu = np.zeros((m+1,))

t = time.time()
bs2 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs2.solve()
elapsed=time.time()-t
print "\tElapsed Time2:", elapsed

ax = fig.add_subplot(1, 3, 2)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);

m = 100
Fl = zeros((m+1,))
Fu = zeros((m+1,))
    
t = time.time()
bs3 = BS_FDM_implicit(r, sigma, T, Smin, Smax, Fl, Fu, payoff, m, n)
px_fd_mat = bs3.solve()
elapsed=time.time()-t
print "\tElapsed Time3:", elapsed

ax = fig.add_subplot(1, 3, 3)
nrow = len(px_fd_mat[:,1])
ax.plot(X, px_fd_mat[nrow-1,:])
plt.xlabel('Stock', fontsize=14);
plt.title('Time Steps: %.f, CPU Time:%.3g'%(m, elapsed), fontsize=14);


#print px_fd_mat[nrow-1,:]

SOLVING for the UP-AND-OUT Barrier Call Using the Implicit FDM method: 
	Elapsed Time1: 0.0620000362396
	Elapsed Time2: 0.0190000534058
	Elapsed Time3: 0.00800013542175

Out[100]:

<matplotlib.text.Text at 0x27b26a90>