Finite-Difference Schemes for the Black-Scholes Equation with Non-smooth Payoff Initial Conditions

Issues with the Greeks in simple schemes

Introduction

[This lecture will be supplemented with C++ code in a lab session provided separately.]
In previous lectures I gave some examples which show how various commonly used finite-difference schemes behave with simple smooth "payoffs". With smooth payoffs the Douglas finite-difference scheme gives much better results than Crank-Nicolson. In this session we will look at how these methods behave when applied to real-world option-pricing problems. This will highlight the potentially nasty behaviour of both Crank-Nicolson and Douglas finite-difference schemes when applied to simple option-pricing problems. The discontinuous nature of the payoff or its derivative in the neighbourhood of the strike induces slowly decaying oscillations into the solution of the finite-difference equations, when we use a larger time-step (which is the main point of introducing FD schemes in the first place.)  These introduce small errors into the valuation itself, and undermine attempts to compute δ, Γ and other "derivative" quantities in the neighbourhood of the strike. Note that, even with the oscillatory components small in the valuation error, the effect on the slope of the function is larger (generating larger errors in δ), and there are even larger errors in Γ.
In this session I also look at the three-time-level version of the Douglas scheme. This chapter shows how the problem of larger errors in the Greeks can be cured, in that the oscillations can be removed, at least for our test payoffs, and the method ("three-time-level Douglas") will be used initially as the basis of our programme to define benchmark numerical algorithms. Later we shall look at some recent work by Evans and Chawla, and Giles and Carter.
The problems with the Crank-Nicolson scheme, and its resolution, are exemplified by a test problem with α = 8. Although the error in the valuation is small, we get an error in Γ near the strike of 12% of its exact value, and an error in θ of 200%. In the corresponding three-time-level Douglas solution, the error in Γ near the strike is reduced to 0.2%, and that in θ to 2.6%.

Finite-Difference Schemes with Three Time-Levels

Various schemes have been proposed for treating problems introduced by considering discontinuous or non-smooth boundary/initial conditions. These schemes aim to damp fast oscillations more effectively, by adjusting the spectrum of eigenvalues of the difference matrix. Richtmyer and Morton (1957) give a list of 14 difference schemes for the diffusion equation, and recommend schemes (their numbering) 9, 11 and 13 for non-smooth initial data. Schemes 9 and 13 are also recommended for these purposes by Smith (1985). These two schemes may be regarded as the three-time-level versions of the Crank-Nicolson scheme and the Douglas scheme, since they have truncation errors of a similar character to their two-time-level counterparts. We already know that the truncation error characteristics of the Douglas scheme make it preferable, so we shall use this. As in previous lectures, we let m denote the time-step, and n denote the x-step. The three time-level Douglas scheme is then given by

((1/8)-α)((Subscript[u, n-1]^(m+1))+(Subscript[u, n+1]^(m+1)))+((5/4)+2α)(Subscript[u, n]^(m+1))
=(1/6)((Subscript[u, n-1]^m)+(Subscript[u, n+1]^m)+10(Subscript[u, n]^m))-(1/24)((Subscript[u, n-1]^(m-1))+(Subscript[u, n+1]^(m-1))+10(Subscript[u, n]^(m-1)))

This type of process requires a kick-off procedure, since initially we only know (u^1). We use the ordinary Douglas two-time-level scheme:

(1 - 6α)((Subscript[u, n-1]^(m+1))+(Subscript[u, n+1]^(m+1))) + (10 + 12α)(Subscript[u, n]^(m+1))=(1 + 6α)((Subscript[u, n-1]^m)+(Subscript[u, n+1]^m)) +(10-12α)(Subscript[u, n]^(m ))

once with α → α/4 , then the three-time-level scheme once with α/4 and then again with α/2.

This gives us our vector pair of vectors to allow the three-time-level iteration to proceed normally thereafter.

Some Required Functions from previous lectures

I am going to define many functions with similar sounding names so I will ask Mathematica to not give spelling warnings (only do this if you have sure there are no mistakes in function names!)

Off[General::spell1]

ExplicitSolver =
	Compile[{{initial,_Real,1}, {lower,_Real,1}, {upper,_Real,1},alpha},
		Module[{wold=initial, wnew=initial, wvold=initial,m,k, tsize = Length[lower],xsize=Length[initial]},
			For[m=2, m <= tsize, m++,
				(wvold = wold; wold=wnew;
					For[k=2, k < xsize, k++,
					(wnew[[k]] = alpha (wold[[k-1]]+wold[[k+1]])+(1-2alpha) wold[[k]])];
						wnew[[1]] = lower[[m]];
						wnew[[xsize]] = upper[[m]])];
			{wvold,wold,wnew}]
		];

CompTridiagSolve = 			Compile[{{a, _Real, 1},{b,_Real, 1},{c, _Real, 1},{r, _Real, 1}},
    Module[{len=Length[r], solution=r, aux = 1/(b[[1]]), aux1=r, a1 = Prepend[a, 0.0], iter}, 	solution[[1]] = aux *r[[1]];
Do[aux1[[iter]] = c[[iter-1]] aux;
			  aux = 1/(b[[iter]] - a1[[iter]]*aux1[[iter]]);
				solution[[iter]] = (r[[iter]]-a1[[iter]] solution[[iter-1]]) aux,
	    {iter, 2, len}]; 				
Do[solution[[iter]] -= aux1[[iter+1]] solution[[iter+1]],
	    {iter, len-1, 1, -1}];
	    solution]];

Adjustment for non-zero boundary conditions

Our implicit schemes involve the solution of a matrix problem A.x = r where A has the tridiagonal form :

A = {{(Subscript[b, 0]),(Subscript[c, 0]),0,…,…,0,0},{(Subscript[a, 1]), (Subscript[b, 1]), (Subscript[c, 1]),0,…,…,0},{0,(Subscript[a, 2]),(Subscript[b, 2]),(Subscript[c, 2]),…,…,0},{0,…,⋱,⋱,⋱,⋱,0},{0,…,…,0,(Subscript[a, "N-2"]),(Subscript[b, "N-2"]),(Subscript[c, "N-2"])},{0,…,…,…,0,(Subscript[a, "N-1"]),(Subscript[b, "N-1"])}};

MatrixForm[A]

In the particular case, for example, of a fully implicit scheme, we have the matrix whose diagonal and off-diagonal terms are

A = {{1+2α,-α,0,…,…,0,0},{-α, 1+2α, -α,0,…,…,0},{0,-α,1+2α,-α,…,…,0},{0,…,⋱,⋱,⋱,⋱,0},{0,…,…,0,-α,1+2α,-α},{0,…,…,…,0,-α,1+2α}};MatrixForm[A]

This came from the difference equation

(1+2 α ) (Subscript[u, n]^(m+1))-α  ((Subscript[u, n-1]^(m+1))+(Subscript[u, n+1]^(m+1)))= (Subscript[u, n]^m)

When we are at the edge of the grid we need to allow for the possibly non-zero value of the boundary condition at the top at bottom end. In the fully implicit case we correct the first and last entries on the right hand side by +α (Subscript[u, p]^(m+1))where p=1 or (Subscript[N, max]). Similarly with other schemes. In the general θ-method case the discrete problem is

(1+2 α θ) (Subscript[u, n]^(m+1))-α θ ((Subscript[u, n-1]^(m+1))+(Subscript[u, n+1]^(m+1)))=(1-2 α (1-θ)) (Subscript[u, n]^m)+α (1-θ) ((Subscript[u, n-1]^m)+(Subscript[u, n+1]^m))

and we will add a term +α θ (Subscript[u, p]^(m+1))+α(1- θ) (Subscript[u, p]^m) at the edges of the grid.NEEDS CARE - possible source of error. You need to think about boundary conditions. Usually such thought pays off. Sometimes a misplaced love of trees can arise from wanting to avoid thinking about boundary conditions. (But then you get bitten by exotics - see Boyle's work on barriers.)

In our calculations we shall just use some simple lists

FullyImpCMatrix[alpha_, nminus_, nplus_] :=
Sequence[Table[-alpha, {nplus+nminus-2}], 
Table[1+2*alpha, {nplus+nminus-1}],
Table[-alpha, {nplus+nminus-2}]]

FullyImpCMatrix[α,2,2]

Case Study: the Vanilla European Put Option

We now begin a detailed study of the first of two examples that we have picked for detailed investigation. You might wonder why we are considering such a trivial case for which there is a known analytic solution. Our goal here is to try out various difference schemes and find out what works well, by testing them on a case for which the solution is known and where the errors can be precisely described. In this way we can see what is happening without all the complications of other real-world effects.

Factors Common to All Our FD Schemes

The Black-Scholes differential equation

(1/2) (S^2) (((∂^2)V)/(∂(S^2))) (σ^2)-r V+(r-q) S ((∂V)/(∂S))+((∂V)/(∂t))=0

with constant coefficients (in particular r is a constant) can be transformed into the diffusion equation (1) by a standard change of variables, as given by Wilmott et al (Chapter 17, equations 17.1, 17.2). We can therefore implement our four standard schemes on the problem, and use solutions of (3) with known payoffs. In the following example we use the vanilla European Put.

Standardization of Variables

With a strike K and constant parameters r, q, σ, we make the changes of variables

τ=(((σ^2) (T-t))/2)       (Subscript[k, 1])=((2 r)/(σ^2))       (Subscript[k, 2])=((2 (r-q))/(σ^2))
V(S, t)=K (e^(-(1/2) ((Subscript[k, 2])-1) x-((1/4) ((((Subscript[k, 2])-1))^2)+(Subscript[k, 1])) τ)) u(x, τ)

The Mathematica implementation of this requires the following functions.

NonDimExpiry[T_,σ_]:=(((σ^2) T)/2);
kone[r_,σ_]:=((2r)/(σ^2));    ktwo[r_,q_,sd_]:=((2 (r-q))/(sd^2));
ValuationMultiplier[strike_,r_,q_,x_,tau_,sd_]:=strike Exp[-(1/2) (ktwo[r,q,sd]-1) x-((1/4) (((ktwo[r,q,sd]-1))^2)+kone[r,sd]) tau]

Initial (Expiry) Conditions

V(S, T)=K (e^(-(1/2) ((Subscript[k, 2])-1) x)) u(x, 0)
C(S, T) = Max(S-K, 0)
P(S, T) = Max(K-S, 0)

So, e.g. for a Put

(Subscript[u, P])(x, 0)=(e^(+(1/2) ((Subscript[k, 2])-1) x)) 1/K Max(K-S,0) =(e^(+(1/2) ((Subscript[k, 2])-1) x)) Max(1-S/K,0) =(e^(+(1/2) ((Subscript[k, 2])-1) x)) Max(1-(e^x),0)

and for a Call

(Subscript[u, C])(x, 0)=(e^(+(1/2) ((Subscript[k, 2])-1) x)) 1/K Max(S-K,0) =(e^(+(1/2) ((Subscript[k, 2])-1) x)) Max(S/K-1,0) =(e^(+(1/2) ((Subscript[k, 2])-1) x)) Max((e^x)-1,0)

The Mathematica functions to do this are:

CallExercise[x_,r_,q_,sd_]:=Max[Exp[(1/2) (ktwo[r,q,sd]-1) x] (Exp[x]-1),0];
PutExercise[x_,r_,q_,sd_]:=Max[Exp[(1/2) (ktwo[r,q,sd]-1) x] (1-Exp[x]),0];

Black-Scholes Model for Verification

From last time - nothing new here:

Ncdf[(z_)?NumberQ] := N[0.5*Erf[z/Sqrt[2]] + 0.5]; 
Ncdf[x_] := (1 + Erf[x/Sqrt[2]])/2;
done[s_, σ_, k_, t_, r_, q_] := 
((r - q)*t + Log[s/k])/(σ*Sqrt[t]) + (σ*Sqrt[t])/2; 
dtwo[s_, σ_, k_, t_, r_, q_] := 
((r - q)*t + Log[s/k])/(σ*Sqrt[t]) - (σ*Sqrt[t])/2;

BlackScholesCall[s_, k_, σ_, r_, q_, t_] := 
s*Exp[-q*t]*Ncdf[done[s, σ, k, t, r, q]] - k*Exp[-r*t]*Ncdf[dtwo[s, σ, k, t, r, q]]; 
BlackScholesPut[s_, k_, σ_, r_, q_, t_] := 
k*Exp[-r*t]*Ncdf[-dtwo[s, σ, k, t, r, q]] - s*Exp[-q*t]*Ncdf[-done[s, σ, k, t, r, q]]

Just check the function calls are working:

{BlackScholesCall[8, 10, 0.2, 0.05, 0, 3], BlackScholesPut[8, 10, 0.2, 0.05, 0, 3]}

Build Greeks automatically for verification

In Mathematica we can just calculcate the derivatives automatically

BlackScholesCallDelta[s_,k_, v_, r_, q_, t_]= 
Evaluate[Simplify[D[BlackScholesCall[s,k, v, r, q, t], s]]];

BlackScholesPutDelta[s_,k_, v_, r_, q_, t_]= 
Evaluate[Simplify[D[BlackScholesPut[s,k, v, r, q, t], s]]];

BlackScholesCallGamma[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesCall[s,k, v, r, q, t], {s, 2}]];

BlackScholesPutGamma[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesPut[s,k, v, r, q, t], {s, 2}]];

BlackScholesCallTheta[s_,k_, v_, r_, q_, t_]= 
-Evaluate[D[BlackScholesCall[s,k, v, r, q, t], t]];

BlackScholesPutTheta[s_,k_, v_, r_, q_, t_]= 
-Evaluate[D[BlackScholesPut[s,k, v, r, q, t], t]];

BlackScholesCallRho[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesCall[s,k, v, r, q, t], r]];

BlackScholesPutRho[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesPut[s,k, v, r, q, t], r]];

BlackScholesCallVega[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesCall[s,k, v, r, q, t], v]];

BlackScholesPutVega[s_,k_, v_, r_, q_, t_]= 
Evaluate[D[BlackScholesPut[s,k, v, r, q, t], v]];

Put Boundary Conditions

For our boundary conditions the upper boundary condition is to set the function to zero, while the lower takes what we expect to be the Put value as the stock price approaches zero. In general the determination of suitable boundary conditions can take quite a bit of thought. It is also a common source of error, though sometimes a bad choice might be washed out by a payoff constraint in the interior of the solution. (Not a good idea to ever repy on that!) So what is the Put value as S becomes small, or equivalently as x becomes large and negative? If S is very small compared to K the Put will be exercised at maturity with a probability that tends to 1 as S tends to zero, so that we can take the present value of K-(Subscript[S, T]). This present value (or solve the BS equation) is given by

(e^(-r (T-t))) K-(e^(-q (T-t))) S

and if we make the transformation to u (I messed this up in my book - the bits below do satisfy the diffusion equation!) you get

(e^((1/4) τ ((((Subscript[k, 2])-1))^2)+(1/2) x ((Subscript[k, 2])-1)))-(e^((1/4) τ ((((Subscript[k, 2])+1))^2)+(1/2) x ((Subscript[k, 2])+1)));

Note that as x gets large and negative the second term is exponentially small compared to the first, but best not to use just the K piece as you have not yet thought about where to truncate the grid.

g[x_,tau_,r_,q_,sd_]:=0;
f[x_,tau_,r_,q_,sd_]:=Exp[(1/2) (ktwo[r,q,sd]-1) x+(1/4) (((ktwo[r,q,sd]-1))^2) tau]-Exp[(1/2) (ktwo[r,q,sd]+1) x+(1/4) (((ktwo[r,q,sd]+1))^2) tau];

Explicit Scheme for Put

dx=0.025; dtau=0.00025; alpha=(dtau/(dx^2));

M=400; nminus=160; nplus=160;

Setting the Initial (Expiry) Condition and Boundary Conditions

initial=Table[PutExercise[(k-1-nminus) dx,0.05,0,0.2], {k,nminus+nplus+1}];
lower=Table[f[-nminus dx,(m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
upper=Table[g[+nplus dx,  (m-1) dtau,0.05,0,0.2], {m, 1, M+1}];

Note that the payoff looks slightly odd in these coordinates - the main point to notice is the discontinuity in slope at the strike.

ListPlot[initial]

Solving the PDE (Explicit Method)

We use the same functions as we used for our smooth test problem in Chapter 14 - load the function ExplicitSolver now if you are using the electronic form:

soln = ExplicitSolver[initial,lower,upper, alpha ];

Interpolating to Supply a Continuous Function

interpoldata=Table[{(k-nminus-1) dx,soln[[3,k]]},{k,1,nminus+nplus+1}];

ufunc=Interpolation[interpoldata,InterpolationOrder->3];

Valuation[strike_,r_,q_,S_,T_,sd_]:=ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] ufunc[Log[(S/strike)]]

Error Plot

Plot[Valuation[10,0.05`,0,S,5,0.2`]-BlackScholesPut[S,10,0.2`,0.05`,0,5],{S,1,20},PlotPoints->50,PlotRange->All]

samples=TableForm[Join[{{"S","Explicit FD", "Exact", "Error"}},
	Table[(PaddedForm[N[#1],{5,5}]&)/@
		{S,Valuation[10,0.05,0,S,5,0.2],
BlackScholesPut[S,10,0.2,0.05,0,5],
						Valuation[10,0.05,0,S,5,0.2]-BlackScholesPut[S,10,0.2,0.05,0,5]},
		{S,2,16,1}
]]]

Fully Implicit Scheme for Put

These algorithms are by now self-explanatory - first the initialization:

dx=0.025; dtau=0.00025; alpha=(dtau/(dx^2));
M=400; nminus=160; nplus=160;

initial=Table[PutExercise[(k-1-nminus) dx,0.05,0,0.2], {k,nminus+nplus+1}];
lower=Table[f[-nminus dx,(m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
upper=Table[g[+nplus dx,  (m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
wold = initial;
wvold = wold; wnew = wold;

CMat = FullyImpCMatrix[alpha,nminus,nplus];

Evolving the solution (this may take some time):

For[m=2, m<=M+1, m++,
(wvold = wold;
wold = wnew;
(* Adjust the rhs for non-zero BCs *)
rhs = Take[wold, {2, -2}]+
Table[
Which[
k==1, alpha*lower[[m]],
k== nplus + nminus-1,  alpha*upper[[m]],
True, 0],
{k, 1, nplus + nminus-1}];
temp = CompTridiagSolve[CMat, rhs];
wnew = Join[{lower[[m]]}, temp, {upper[[m]]}])
]

Interpolation and construction of valuation function:

interpoldata = Table[{(k - nminus - 1)*dx, wnew[[k]]}, 
    {k, 1, nminus + nplus + 1}];

ufuncb = Interpolation[interpoldata, InterpolationOrder -> 3];

Valuation[strike_, r_, q_, S_, T_, sd_] := 
  ValuationMultiplier[strike, r, q, Log[S/strike], (sd^2*T)/2, sd]*
   ufuncb[Log[S/strike]]

Error Plot

Plot[Valuation[10,0.05`,0,S,5,0.2`]-BlackScholesPut[S,10,0.2`,0.05`,0,5],{S,1,20},PlotPoints->50,PlotRange->All]

samples=TableForm[Join[{{"S","Implicit FD", "Exact", "Error"}},
	Table[(PaddedForm[N[#1],{5,5}]&)/@
		{S, Valuation[10,0.05,0,S,5,0.2],
						BlackScholesPut[S,10,0.2,0.05,0,5],
						Valuation[10,0.05,0,S,5,0.2]-BlackScholesPut[S,10,0.2,0.05,0,5]},	{S,2,16,1}]]]

Note that we obtain no improvement in accuracy over the explicit scheme. The only advantage of the fully implicit scheme over the explicit scheme is the fact that we can, if we wish, increase α, that is, the time-step for a given price-step, without the system going unstable.

Crank-Nicolson

We now increase the time-step by a factor of 20, so that α=8. Otherwise all this proceeds as before:

Initialization:

M=20; nminus = 160; nplus = 160;
dx = 0.025; dtau = 0.005; alpha = dtau/dx^2

initial=Table[PutExercise[(k-1-nminus) dx,0.05,0,0.2], {k,nminus+nplus+1}];
lower=Table[f[-nminus dx,(m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
upper=Table[g[+nplus dx,  (m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
wold = initial;
wvold = wold;
wnew = wold;

CNCMatrix[alpha_, nminus_, nplus_] :=
Sequence[Table[-alpha/2, {nplus+nminus-2}],
Table[1+alpha, {nplus+nminus-1}],
Table[-alpha/2, {nplus+nminus-2}]];

CNDMatrix[alpha_, vec_List] := Module[{temp},
temp = (1 - alpha)*vec + (alpha/2)*(RotateRight[vec] + RotateLeft[vec]);
temp[[1]] = Simplify[First[temp] - 
alpha*Last[vec]/2];
temp[[-1]] = Simplify[Last[temp] - alpha*First[vec]/2];
temp];

CMat = CNCMatrix[alpha,nminus,nplus];

Evolution:

For[m=2, m<=M+1, m++,
(wvold = wold;
wold = wnew;
rhs = CNDMatrix[alpha, Take[wold, {2, -2}]]+
Table[
Which[
k==1, alpha*(lower[[m-1]] + lower[[m]])/2,
k== nplus + nminus-1,  alpha*(upper[[m-1]] + upper[[m]])/2,
True, 0],
{k, 1, nplus + nminus-1}];
temp = CompTridiagSolve[CMat, rhs];
wnew = Join[{lower[[m]]}, temp, {upper[[m]]}])
]

Interpolation

interpoldatab = 
Table[{(k - nminus - 1)*dx, wnew[[k]]}, {k, 1, nminus+nplus+1}];

ufuncb = Interpolation[interpoldatab, InterpolationOrder -> 3];

Valuation[strike_, r_, q_, S_, T_, sd_] :=
ValuationMultiplier[strike, r, q, Log[S/strike], (sd^2*T)/2, sd]*ufuncb[Log[S/strike]]

Error Plot

Plot[Valuation[10,0.05`,0,S,5,0.2`]-BlackScholesPut[S,10,0.2`,0.05`,0,5],{S,1,20},PlotPoints->50,PlotRange->All]

samples=TableForm[Join[{{"S","Crank-Nic", "Exact", "Error"}},
	Table[(PaddedForm[N[#1],{5,5}]&)/@
		{S, Valuation[10,0.05,0,S,5,0.2],
						BlackScholesPut[S,10,0.2,0.05,0,5],
						Valuation[10,0.05,0,S,5,0.2]-BlackScholesPut[S,10,0.2,0.05,0,5]},{S,2,16,1}]]]

Note that we obtain comparable accuracy to the explicit method with 1/20 the number of time-steps, but near the strike things are becoming "interesting". We obtain a sharply oscillatory error in the neighbourhood of the strike. Note that the size of the error is quite small, but it is steeply sloped. Bearing in mind the warnings given in our early lecture this is where we need to start thinking about getting the Greeks numerically..

Construction and Verification of Interpolated Valuation and Greeks

We begin by writing down a function that does numerical differentiation of a list. This function uses a simple central difference algorithm for points in the interior of the list. The end points are treated using a special difference algorithm.

listd[data_, step_] :=
Module[{dleft,dright,len},
len = Length[data];
dleft = (4*data[[2]]-3*data[[1]]-data[[3]])/(2*step);
dright = (3*data[[len]]-4*data[[len-1]]+data[[len-2]])/(2*step);
Join[{dleft}, Take[RotateLeft[data]-RotateRight[data], {2, -2}]/(2*step), {dright}]
]

Next we fix the values of the k-parameters:

kt = ktwo[0.05, 0, 0.2];
ko = kone[0.05, 0.2];

deltadata=listd[wnew,dx]-(1/2) (kt-1) wnew;

gammadata=listd[deltadata,dx]-(1/2) (kt+1) deltadata;

thetadata=((3wnew-4wold+wvold)/(2dtau))-((1/4) (((kt-1))^2)+ko)wnew;

points = Table[(k - nminus - 1)*dx, {k, 1, nminus + nplus + 1}];

deltainterpoldata=Transpose[{points,deltadata}];
gammainterpoldata=Transpose[{points,gammadata}];
thetainterpoldata=Transpose[{points,thetadata}];

dfunc=Interpolation[deltainterpoldata,InterpolationOrder->3];
gfunc=Interpolation[gammainterpoldata,InterpolationOrder->3];
tfunc=Interpolation[thetainterpoldata,InterpolationOrder->3];

CNDelta[strike_,r_,q_,S_,T_,sd_]:=(1/S)(ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] dfunc[Log[(S/strike)]])

CNGamma[strike_,r_,q_,S_,T_,sd_]:=(1/(S^2))(ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] gfunc[Log[(S/strike)]])

CNTheta[strike_,r_,q_,S_,T_,sd_]:=-(1/2) (sd^2) ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] tfunc[Log[(S/strike)]]

Delta Analysis

To investigate the errors in Delta we make a comparison with the exact solution, which we have already loaded as a function:

?BlackScholesPutDelta

We make a table of values of the percentage errors in delta, given in the last column of the following table:

deltasamples = TableForm[
Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, CNDelta[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutDelta[S,10,0.2,0.05,0, 5],
100*(CNDelta[10,0.05,0,S,5,0.2]/BlackScholesPutDelta[S,10,0.2,0.05,0, 5]-1)}],
{S, 2, 16, 1}]]

The error in Delta is at most 0.2 % of the exact value, near the strike. This is still not too bad, and we can plot the error in Delta.

Plot of Computed and Exact Delta

If we overlay the computed and exact Delta, the wobbles are starting to show:

Plot[{CNDelta[10,0.05`,0,S,5,0.2`],BlackScholesPutDelta[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

Plotting the difference makes the problem clear:

Plot[CNDelta[10,0.05`,0,S,5,0.2`]-BlackScholesPutDelta[S,10,0.2`,0.05`,0,5],{S,8,12},PlotPoints->50,PlotRange->All]

Gamma Analysis

Here we tabulate the exact and numerical solution, with the percentage error in Gamma in the last column:

gammasamples = TableForm[Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, CNGamma[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutGamma[S,10,0.2,0.05,0, 5],
100*(CNGamma[10,0.05,0,S,5,0.2]/BlackScholesPutGamma[S,10,0.2,0.05,0,5]-1)}],{S, 2, 16, 1}]]

The error in Γ is 12.4% of the exact value at the strike.

Gamma Plot

The wobbles in Γ are now manifest in a plot of the computed and exact version - the error is of a similar scale to the value of Γ:

Plot[{CNGamma[10,0.05`,0,S,5,0.2`],BlackScholesPutGamma[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

Theta Analysis

Now we make a table of percentage errors with the Crank-Nicolson scheme:

thetasamples = TableForm[
Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, CNTheta[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutTheta[S,10,0.2,0.05,0, 5],
100*(CNTheta[10,0.05,0,S,5,0.2]/BlackScholesPutTheta[S,10,0.2,0.05,0,5]-1)}], {S, 2, 16, 1}]]

Near the strike the error in Theta peaks at 197% of the exact value.

Theta Plots

Plot[{CNTheta[10,0.05`,0,S,5,0.2`],BlackScholesPutTheta[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

Remark on Two-Time-Level Douglas

This time, there is no benefit in going to the simple Douglas scheme as described in Chapter 14. The oscillations in the neighbourhood of the strike are just as bad. For example, here is the error in valuation with α=8, for comparison. There are corresponding problems in the Greeks.

Douglas Three-Time-Level Solution

We have to work a little harder to develop a scheme that eliminates the oscillation problem. What we do is to implement the three-time-level extension of the Douglas scheme discussed previously. We shall not go into a detailed theoretical discussion of this algorithm - we shall content ourselves with an explicit demonstration that it works.

Necessary Functions

DougCMatrix[alpha_, nminus_, nplus_] :=
Sequence[Table[1-6*alpha, {nplus+nminus-2}], Table[10+12*alpha, {nplus+nminus-1}],
Table[1-6*alpha, {nplus+nminus-2}]]

DougCCMatrix[alpha_, nminus_, nplus_] :=
Sequence[Table[1/8-alpha, {nplus+nminus-2}], 
Table[5/4+2*alpha, {nplus+nminus-1}],
Table[1/8-alpha, {nplus+nminus-2}]]

DougDMatrix[alpha_, vec_List] := Module[{temp},
temp = (10 - 12*alpha)*vec + (1+6*alpha)*(RotateRight[vec] + RotateLeft[vec]);
temp[[1]] = Simplify[First[temp] - (1+6*alpha)*Last[vec]];
temp[[-1]] = Simplify[Last[temp] - (1 + 6*alpha)*First[vec]];
temp]

DougDDMatrix[vec_List] := Module[{temp},
temp = (10*vec + RotateRight[vec] + RotateLeft[vec]);
temp[[1]] = Simplify[First[temp] - Last[vec]];
temp[[-1]] = Simplify[Last[temp] - First[vec]];
temp/6]

Douglas Three Time-Level Solution Evolution

The initialization consists of defining the grid parameters, setting boundary and initial conditions and defining the various matrices.

dx = 0.025;
dtau = 0.005;
alpha = dtau/dx^2;
M=20;
nminus = 160;
nplus = 160;
alpha

initial=Table[PutExercise[(k-1-nminus) dx,0.05,0,0.2], {k,nminus+nplus+1}];
lower=Table[f[-nminus dx,(m-1) dtau,0.05,0,0.2], {m, 1, M+1}];
upper=Table[g[+nplus dx,  (m-1) dtau,0.05,0,0.2], {m, 1, M+1}];

w = Table[0, {m, 1, 3},   {k, 1, nminus+nplus+1}];

vold = Take[initial, {2, -2}];

w[[1,1]] = f[-nminus*dx, dtau/4, 0.05, 0, 0.2];
w[[1, nminus+nplus+1]] = g[nplus*dx, dtau/4, 0.05, 0, 0.2];
w[[2,1]] = f[-nminus*dx, dtau/2, 0.05, 0, 0.2];
w[[2, nminus+nplus+1]] = g[nplus*dx, dtau/2, 0.05, 0, 0.2];
w[[3,1]] = f[-nminus*dx, dtau, 0.05, 0, 0.2];
w[[3, nminus+nplus+1]] = g[nplus*dx, dtau, 0.05, 0, 0.2];

CMat = DougCMatrix[alpha,nminus,nplus];
CCMat = DougCCMatrix[alpha,nminus,nplus];

CMatQ = DougCMatrix[alpha/4,nminus,nplus];
CMatH = DougCMatrix[alpha/2,nminus,nplus];

CCMatQ = DougCCMatrix[alpha/4,nminus,nplus];
CCMatH = DougCCMatrix[alpha/2,nminus,nplus];

Kick-Off Phase

This begins with the simple Douglas scheme with two time-levels and 1/4 the basic time-step.

rhs = DougDMatrix[alpha/4, vold]+
Table[
Which[
k==1, (6*alpha/4+1)*lower[[1]] + (6*alpha/4-1)*w[[1, 1]],
k== nplus + nminus-1,  (6*alpha/4+1)*upper[[1]] +
(6*alpha/4-1)*w[[1, nplus+nminus+1]],
True, 0],
{k, 1, nplus + nminus-1}];
vnew = CompTridiagSolve[CMatQ, rhs];
w[[1]] = Join[{w[[1,1]]}, vnew, {w[[1,nplus+nminus+1]]}];
vvold = vold;
vold = vnew;

Now we have two iterations of the three-time-level Douglas system, doubling the time-step at each stage:

rhs = DougDDMatrix[vold] - DougDDMatrix[vvold]/4 + 
Table[
Which[
k==1, (alpha/4-1/8)*w[[2, 1]] + w[[1,1]]/6 - lower[[1]]/24,
k==nplus + nminus-1,  (alpha/4-1/8)*w[[2, nplus+nminus+1]] + 
w[[1,nplus+nminus+1]]/6 - upper[[1]]/24,
True, 0],
{k, 1, nplus + nminus-1}];
vnew = CompTridiagSolve[CCMatQ, rhs];
w[[2]] = Join[{w[[2,1]]}, vnew, {w[[2,nplus+nminus+1]]}];
vold = vnew;

rhs = DougDDMatrix[vold] - DougDDMatrix[vvold]/4 + 
Table[
Which[
k==1, (alpha/2-1/8)*w[[3, 1]] + w[[2,1]]/6 - lower[[1]]/24,
k==nplus + nminus-1,  (alpha/2-1/8)*w[[3, nplus+nminus+1]] + 
w[[2,nplus+nminus+1]]/6 - upper[[1]]/24,
True, 0],
{k, 1, nplus + nminus-1}];
vnew = CompTridiagSolve[CCMatH, rhs];
w[[3]] = Join[{w[[3,1]]}, vnew, {w[[3,nplus+nminus+1]]}];
wold = initial;
wvold = initial;
wnew = w[[3]];

Main Evolution Phase and Interpolation of Solution

For[m=3, m<=M+1, m++,
(wvold = wold;
wold = wnew;
rhs = DougDDMatrix[Take[wold, {2, -2}]] -
DougDDMatrix[Take[wvold, {2, -2}]]/4 + 
Table[
Which[
k==1, 
(alpha-1/8)*lower[[m]] + lower[[m-1]]/6 - lower[[m-2]]/24,
k==nplus + nminus-1, 
(alpha-1/8)*upper[[m]] + upper[[m-1]]/6 - upper[[m-2]]/24,
True, 0],
{k, 1, nplus + nminus-1}];
temp = CompTridiagSolve[CCMat, rhs];
wnew = Join[{lower[[m]]}, temp, {upper[[m]]}])
]

points = Table[(k-nminus-1)*dx,  
{k, 1, nplus+nminus+1}];

finalstep = wnew;

prevstep = wold;

pprevstep = wvold;

interpoldata = Transpose[{points,finalstep}];

ufunc = Interpolation[interpoldata, InterpolationOrder -> 3];

Valuation[strike_, r_, q_, S_, T_, sd_] :=
ValuationMultiplier[strike, r, q, Log[S/strike], sd^2*T/2, sd]*
ufunc[Log[S/strike]]

Valuation Analysis

Plot[Valuation[10,0.05`,0,S,5,0.2`]-BlackScholesPut[S,10,0.2`,0.05`,0,5],{S,1,20},PlotPoints->50,PlotRange->All]

samples = TableForm[
Table[Map[PaddedForm[N[#], {5, 4}]&, {S, Valuation[10, 0.05, 0, S, 5, 0.2],
BlackScholesPut[S,10,0.2,0.05,0,5],
Valuation[10, 0.05, 0, S, 5, 0.2]-
BlackScholesPut[S,10,0.2,0.05,0,5]}],
{S, 2, 16, 1}]]

Construction and Verification of Interpolated Valuation and Greeks

In this section we define some functions to compute Delta, Gamma and Theta directly from the finite-difference grid.

listd[data_, step_] :=
Module[{dleft,dright,len},
len = Length[data];
dleft = (4*data[[2]]-3*data[[1]]-data[[3]])/(2*step);
dright = (3*data[[len]]-4*data[[len-1]]+data[[len-2]])/(2*step);
Join[{dleft}, Take[RotateLeft[data]-RotateRight[data], {2, -2}]/(2*step), {dright}]
]

kt = ktwo[0.05, 0, 0.2];
ko = kone[0.05, 0.2];

deltadata=listd[finalstep,dx]-(1/2) (kt-1) finalstep;

gammadata=listd[deltadata,dx]-(1/2) (kt+1) deltadata;

thetadata=((3finalstep-4prevstep+pprevstep)/(2dtau))-((1/4) (((kt-1))^2)+ko)finalstep;

deltainterpoldata=Transpose[{points,deltadata}];
gammainterpoldata=Transpose[{points,gammadata}];
thetainterpoldata=Transpose[{points,thetadata}];

dfunc=Interpolation[deltainterpoldata,InterpolationOrder->3];
gfunc=Interpolation[gammainterpoldata,InterpolationOrder->3];
tfunc=Interpolation[thetainterpoldata,InterpolationOrder->3];

DougDelta[strike_,r_,q_,S_,T_,sd_]:=(1/S)(ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] dfunc[Log[(S/strike)]])

DougGamma[strike_,r_,q_,S_,T_,sd_]:=(1/(S^2))(ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] gfunc[Log[(S/strike)]])

DougTheta[strike_,r_,q_,S_,T_,sd_]:=-(1/2) (sd^2) ValuationMultiplier[strike,r,q,Log[(S/strike)],(((sd^2) T)/2),sd] tfunc[Log[(S/strike)]]

Delta Analysis

The percentage error in Delta is given in the last column of the following table.

deltasamples = TableForm[
Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, DougDelta[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutDelta[S,10,0.2,0.05,0,5],
100*(DougDelta[10,0.05,0,S,5,0.2]/BlackScholesPutDelta[S,10,0.2,0.05,0,5]-1)}],
{S, 2, 16, 1}]]

The error in Delta is at most 0.1% of the exact value.

Delta Plots

Plot[{DougDelta[10,0.05`,0,S,5,0.2`],BlackScholesPutDelta[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

Gamma Analysis

The percentage error in Gamma is given by the last column of the following table:

gammasamples = TableForm[Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, DougGamma[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutGamma[S,10,0.2,0.05,0,5],
100*(DougGamma[10,0.05,0,S,5,0.2]/BlackScholesPutGamma[S,10,0.2,0.05,0,5]-1)}],
{S,2,16,1}]]

The error is at most 0.7% of the exact value - near the strike it is now only 0.2%.

Plot of Computed and Exact Gamma

Plot[{DougGamma[10,0.05`,0,S,5,0.2`],BlackScholesPutGamma[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

Theta Analysis

The percentage error in Theta is given by the last column of the following table:

thetasamples = TableForm[
Table[Map[PaddedForm[N[#], {5, 4}]&, 
{S, DougTheta[10, 0.05, 0, S, 5, 0.2],
BlackScholesPutTheta[S,10,0.2,0.05,0,5],
100*(DougTheta[10,0.05,0,S,5,0.2]/BlackScholesPutTheta[S,10,0.2,0.05,0,5]-1)}],
{S, 2, 16, 1}]]

Near the strike the error in θ peaks at 2.6% of the exact value. Note that with the two-time-level Crank-Nicolson scheme the error was almost 200%.

Theta Plots

Plot[{DougTheta[10,0.05`,0,S,5,0.2`],BlackScholesPutTheta[S,10,0.2`,0.05`,0,5]},{S,8,12},PlotPoints->50,PlotRange->All]

As a reminder, here is the corresponding plot for the Crank-Nicolson scheme.

Summary

The introduction of implicit schemes is motivated by the desire to get accurate and stable solutions for larger values of α. When we need to compute the values of both a function and its first and second derivatives, popular schemes such as the ordinary two-time-level Crank-Nicolson scheme are unsuitable for option-pricing problems, unless small to moderate values of α are used. This is because of the fact that non-smoothness in initial data (almost always present in option payoffs) propagates through the solution causing small oscillatory errors. When amplified by the process of differentiation, these may cause substantial errors in the Greeks. The errors in the valuation may be deceptively small - to quote Wilmott et al (1993) when discussing a verification example with the Crank-Nicolson scheme, they remark that "Even with α=10, the numerical and exact results differ only marginally." This is absolutely right, at least when we just look at the valuation. When we inspect the Greeks a rather more disturbing picture emerges. We have given an example where estimates of θ from a computation with α=8 are in error by about 200%. 

The oscillation problem and the resulting corruption of the Greeks can be cured by the use of a more suitable difference scheme. When the initial data were smooth, we saw previously that the two-time-level Douglas scheme was sufficient, but in the presence of non-smooth data the corresponding three-time-level Douglas scheme cures the problems quite dramatically. 

The transition to such a scheme for practitioners should not be a great leap. We have already seen that the two-time-level Douglas scheme is the natural implicit generalization of the trinomial model, when carried out on a grid rather than a tree. The next step to a three-time-level scheme allows larger time steps still to be taken without a loss of accuracy or corruption of the Greeks.