#!/usr/bin/env python
# coding: utf-8

# In[1]:


get_ipython().run_line_magic('pylab', 'inline')
lecture = 12


# #Lecture 12: Numerical Methods for Partial Differential Equations
# 
# ## Topics
# 
# * I. Introduction
#     * Second Order PDE Classification
#     * Initial and Boundary Conditions
#     * Black-Scholes Equation

# * II. Finite Difference Methods for Solving PDEs
#     * Discretization
#     * Finite Difference Methods
#     * The Explicit Methods
#     * Consistency, Stability and Convergence
#     * The Implicit Methods
#     * The Crank-Nicholson Method
#     
# * III. Numerical Examples
# 
# 

# # I. Introduction
# 
# 
# * In finance, most partial differential equations (PDE) we encounter are of the second order (the highest derivative)
# 
# $$
# \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}}
# 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
# }
# $$
# 
# * No steady state and no diffusion, mostly referred to as the **Wave equation** or **Advection equation**
# 

# ### Parabolic 
# 
# 
# * $b^2(x,t) - a(x,t)c(x,t) = 0 $
#      $\;$
# 
# * The canonical form:
# $$
# \large{
# \PDut = d \;\PDuxx
# }
# $$
# 
# 
# * Evolving to a steady state due to diffusion, referred to as the **Heat equation**
# 

# ### Elliptic 
# 
# 
# * $b^2(x,y) - a(x,y)c(x,y) < 0 $
#      $\;$
# 
# * The canonical form:
# $$
# \large{
# \PDuxx + \PDuyy = f(x,y)
# }
# $$
# 
# 
# * Steady state (time independent), mostly referred to as the **Laplace** ($f(x,y) = 0$) or **Poisson** equation.

# ## Initial and Boundary Conditions 
# 
# 
# * Depending on the problem, we usually have **initial conditions (IC)**
# 
# $$
# \large{
# u(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 $a\cdot b \neq 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}
# 
# 
# 

# # II. 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$
# 
# $$
# \small
# \PDut - \frac{1}{2}\sigma^2S^2\PDuSS - rS\PDuS + ru = 0
# $$
# 
# * Approximate the Black-Scholes Eqn by
# 
# $$
# \small
# \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}
# \small
# 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 accurate 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 satisfied 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 }
# $$
# 
# 
# * This is the CFL condition for the explicit method, which says, in order for the explicit method to be stable, the time step size needs to be really 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
# 
# $$
# \small
# \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}
# \small
# 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}
# \small
# \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)
# $$
# 
# 

# # III. Numerical Examples 
# 
# * Will go through the Black Scholes example in next lecture.
# 
# 
# 
# 

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

# # Reference
# 
# * L. Andersen and V. Piterbarg, Chapter 4 of Interest Rate Modeling, Volumn I, Atlantic Financial Press, 2010.
# 
# * D. Duffy, Finite Difference Method in Financial Engineering, John Wiley & Sons, 2006.