Numerics For Black Scholes

We have the PDE 
	({{(∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V, =, 0}})
and the Boundary Conditions
	({{V[0,t], =, 0                                 }, {V[s,t], (Underscript[⟶, s->∞]), s                                  }, {v[s, T], =, max[s-k,0]}})
We are going to fake the asymptotic boundary condition by picking a maximum possible stock price which we will call S.  So out boundary conditions are
	({{V[0,t], =, 0                                 }, {V[S,t], =, S                                  }, {v[s, T], =, max[s-k,0]}})

Finite Differences (FD) are Straight forward

Pick n and m and then define
	({{Δt, =, T/n}, {Δs, =, S/m}})
and define values (Subscript[𝕍, i,j])≈V[i Δs, j Δt].  We are going to define the values (Subscript[V, i,j]) by the linear equations we get by plugging finite difference approximations for the derivatives into the PDE.  Note that the BC gives us values for
	({{(Subscript[𝕍, 0,j]), =, 0}, {(Subscript[𝕍, m,j]), =, S}, {(Subscript[𝕍, i,n]), =, max[i Δs-K,0]}})
We can visualize the point approximations we are trying to compute quite easily.

{m,n}={12,8};
{S,T}={100.0,64};
{Δs,Δt}={S,T}/{m,n};
Graphics[{
PointSize[0.03],
{Pink,Rectangle[ {0,0},{S,T}]},
{Black, Table[ Point[{i Δs, j Δt}],
{i, 1, m-1},
{j, 0, n-1}]},
{Green, Table[Point[{S,j Δt}],{j,0,n-1}]},
{Cyan, Table[Point[{0,j Δt}],{j,0,n-1}]},
{Red, Table[Point[{i Δs,T}],{i,1,m-1}]}
},
Frame->True,
FrameLabel->{"s","t"}
]

Note: the grid points are nicely lined up and it is not clear what values to pick at the Green/Red and Cyan/Red corners. At each of the interior black points we simply replace the partial derivatives with appropriate finite difference approximations. It is worthwhile noting that at each black point except the ones on the bottom we have top and bottom and left and right neighbors.  At the bottom blck dots we only have top, left, and right neighbors.

The first thing we could do is use centered differences in s and a backwards difference in t. These approximations are 
	({{(∂_s)V(Subscript[|, {i,j}]), ≈, (((Subscript[𝕍, i+1,j])-(Subscript[𝕍, i-1,j]))/(2 Δs))}, {(∂_ss)V(Subscript[|, {i,j}]), ≈, (((Subscript[𝕍, i+1,j])-2 (Subscript[𝕍, i,j])+(Subscript[𝕍, i-1,j]))/((Δs)^2))}, {(∂_t)V(Subscript[|, {i,j}]), ≈, (((Subscript[𝕍, i,j])-(Subscript[𝕍, i,j-1]))/(Δt))}})
Making these replacements in the PDE
	({{(∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V, =, 0}})
gives (at each of the interior i, j value i.e black)
	({{(((Subscript[𝕍, i,j])-(Subscript[𝕍, i,j-1]))/(Δt))+ 0.5 (((σ^2)(i Δs))^2) ((((Subscript[𝕍, i+1,j])-2 (Subscript[𝕍, i,j])+(Subscript[𝕍, i-1,j]))/((Δs)^2)))+r (i Δs) ((((Subscript[𝕍, i+1,j])-(Subscript[𝕍, i-1,j]))/(2 Δs)))- r (Subscript[𝕍, i,j]), =, 0}})
Along with the know values on the Red/Blue/Cyan boundaries we have the same number of linear equations as unknowns in a linear system.   The easiest way to organize solving these equations would be to solve for (Subscript[𝕍, i,j-1]) and march down from the top boundary. 
	({{(Subscript[𝕍, i,j-1]), =, (Subscript[𝕍, i,j]), +, Δt (0.5 (σ^2)(((i Δs))^2) ((((Subscript[𝕍, i+1,j])-2 (Subscript[𝕍, i,j])+(Subscript[𝕍, i-1,j]))/((Δs)^2)))+r i Δs ((((Subscript[𝕍, i+1,j])-(Subscript[𝕍, i-1,j]))/(2 Δs)))- r (Subscript[𝕍, i,j]))}})
This is an explicit Finite Difference (FD) scheme because we can simply compute the next layer in j from the known j layers. It is convenient to define matrices that do the differencing in the i direction.

SecondDerMat[m_] :=SparseArray[{
Band[{1,1}]->-2,
Band[{1,2}]->-1,
Band[{2,1}]->-1
},
{m-1,m-1}]
M=SecondDerMat[m];
MatrixForm[M]

You need to make sure that you take small steps in the time direction to make this explicit  computation stable.

FD Semi-Discretization is also fairly simple.

All I do now is discretize in the “s” variable using the same differences. I do not do anything with the time derivatives. The result is a set of coupled Ordinary Differential Equations (ODE) which can be solved with an ODE solver.  

This looks like 
The first thing we could do is use centered differences in s and a backwards difference in t. These approximations are 
	({{(∂_s)V(Subscript[|, i]), ≈, (((Subscript[𝕍, i+1])-(Subscript[𝕍, i-1]))/(2 Δs))}, {(∂_ss)V(Subscript[|, i]), ≈, (((Subscript[𝕍, i+1])-2 (Subscript[𝕍, i])+(Subscript[𝕍, i-1]))/((Δs)^2))}})
Making these replacements in the PDE
	({{(∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V, =, 0}})
gives (at each of the interior i values) 
	({{((Subscript[d𝕍, i])/dt)+ 0.5 (σ^2)((i )^2) ((Subscript[𝕍, i+1])-2 (Subscript[𝕍, i])+(Subscript[𝕍, i-1]))+0.5r i  ((Subscript[𝕍, i+1])-(Subscript[𝕍, i-1]))- r (Subscript[𝕍, i]), =, 0}})

Clear[V,σ,r,K,T,S]
m=24;

{S,T}={100.0,64};
Δs=S/m;
{σ,r,K}={0.1,0.01,23};
V[0][t_]=0;
V[m][t_]=S;
ODESys= Table[{
V[i]'[t]+0.5 (σ^2) (i^2)(V[i+1][t]-2 V[i][t]+V[i-1][t]) + 0.5 r i (V[i+1][t]-V[i-1][t])- r V[i][t]==0,
V[i][T]==Max[i Δs -K, 0]
},
{i,1,m-1}];

SubThing = NDSolve[
ODESys,
Table[ V[i],{i,1,m-1}],
{t, 0, T}][[1]];

Sols=Table[
{t,i Δs,V[i][t]},
{i, 0, m}]/.SubThing;
ParametricPlot3D[Sols,{t, 0, T},
AxesLabel->{"t","s","V"}]

One can put them all together as single function

Clear[VSD,s,t]
VSD[s_,t_] =Interpolation[ 
Table[
{i Δs,V[i][t]}/.SubThing,
{i, 0, m}],
InterpolationOrder->1
][s];
ContourPlot[
VSD[s,t],
{s, 0, S},
{t, 0,T},
ContourLabels->True,
FrameLabel->{"s","t"}
]
Plot3D[
VSD[s,t],
{s, 0, S},
{t, 0,T},
AxesLabel->{"s","t","V"}
]

It is pretty easy to plot these guys.

Finite Elements Semi-Discretization is also fairly simple

Now you start by choosing a set of shape functions (traditionally written (Subscript[{(Subscript[ϕ, i])[s]}, i=1]^m))  in the S direction.  Our solution is going to be a time varying linear combination of these shape functions i.e.
	V=(∑_(i=1)^m)(Subscript[α, i])[t](Subscript[ϕ, i])[s]
We are going to plug these into the PDE
	({{(∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V, =, 0}})
but we need to work out how to get m ODEs out (one for each of the αs) after we have plugged in.

One standard way to do this is to multiply the PDE by a test function w(s) (which in the end will look like one of our shape functions) and integrate in s.  
	({{w((∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V), =, 0}, {(∫_(s=0)^S)w((∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V)ds, =, 0}, {(∫_(s=0)^S)w((∂_t)V+r s (∂_s)V- r V)+0.5 (σ^2)((s^2)w (∂_ss)V)ds, =, 0}})
One last trick is to do an integration by parts on the pink term
	({{(∫_(s=0)^S)(s^2)w (∂_ss)V ds, =, (s^2)w (∂_s)V}})(Subscript[|, s=0]^S) - (∫_(s=0)^S)(∂_s)((s^2)w) (∂_s)V ds
or expanding the derivatives and evaluating the boundary terms 
	({{(∫_(s=0)^S)(s^2)w (∂_ss)V ds, =, (S^2)w[S] (∂_s)V[S,t]}}) - (∫_(s=0)^S)((s^2)(∂_s)w+ 2s w)(∂_s)V ds
If you plug in the form V=(∑_(i=1)^m)(Subscript[α, i])[t](Subscript[ϕ, i])[s]  and choose w=(Subscript[ϕ, j]) you get an ODE.  If you look at the system of ODEs you get by choosing (Subscript[ϕ, 1]), (Subscript[ϕ, 2]),(Subscript[ϕ, 3]),…,(Subscript[ϕ, m]) then you get a system of m 1st order ODEs for the  coefficients (Subscript[α, i])[t].  The only thing that we have not thought through are how to incorporate the boundary conditions for in the s direction
	 ({{V[0,t], =, 0                                 }, {V[S,t], =, S                                  }})
It is easy to implement the first condition:  just choose all your shape functions to be zero at s=0.  If the boundary condition at s=S was was homogeneous you could do the same thing there.   We are going to do a simple trick (not sure if you remember but we talked about this before) to make the boundary condition homogeneous. 

Put W[s,t] = V[s,t] - s.  Then 
	({{W[0,t], =, V[0,t]-0, =, 0}, {W[S,t], =, V[S,t]-S, =, 0}})
The terminal condition needs to be modified as does the PDE
	({{(∂_t)V+ 0.5 (σ^2)(s^2) (∂_ss)V+r s (∂_s)V- r V, =, 0}})
To find the new terminal condition and PDE for W just plug
	V[s,t]=W[s,t]+s 
into the PDE and compute terms:  You will wind up with terms on the RHS.
The terminal conditions are simply
	({{W[s,T], =, V[s,T]-s=Max[s-K,0]-s}})