Lecture 12: Numerical Methods for Partial Differential Equations

Topics

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
- Consistency, Stability and Convergence
- The Implicit Methods
- The Crank-Nicholson Method

Numerical Examples

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}} \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) } $$

It can be classified into three categories: 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,t) - a(x,t)c(x,t) < 0 $

$\;$
The canonical form:

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

Mostly referred to as the Laplace or Poisson equation

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$, or Robin if $c = 0$ .

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

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}

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.

Von Neumann Stability Analysis

The Von Neumann 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 Von Neumann 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.

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

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}

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

Numerical Examples

Will go through the Black Scholes example in next lecture.

Higher Dimension Methods

What to do when there are more than one spatial variables (Heston Stochastic Vol Model)? - ADI methods.

Homework

Derive the CFL condition for the Crank-Nicholson method for the Black-Scholes PDE, assuming $r=0$.