Analysing the Crank-Nicolson Method¶
Click here for the interactive version of this notebook.
Here we aim to study the convergence properties of a Crank-Nicolson method solving the European put option problem in finance. We know from the theory that the truncation errors of the method are O((Δt)2,(ΔS)2), so what happens in reality? There are many sources of error in our algorithm, in particular
- finite difference approximations, which are O((Δt)2,(ΔS)2)
- nonlinearity errors (caused by discontinuity in payoff and any barriers or early exercise boundaries)
- truncation errors in the grid (choice of
SMax) - interpolation errors (since S0 is not necessarily a grid point)
- SOR matrix solve (bounded by tolerance) The key point to note here is that any algorithm is only as good as its worst performing error, so improving one or two errors and not improving any of the others is useless. The most efficient implementation of the scheme will make sure each of the errors is roughly of a similar magnitude.
We solve the problem ∂P∂t+12σ2S2∂2P∂S2+rS∂P∂S−rP=0 with P(S,T)=max and S_0 = 1.973 \text{ , } X=2 \text{ , } T=1 \text{ , } r=0.05 \text{ , } \sigma=0.4
To implement our analysis we wish to put our algorithm into a function, so that we can easily call it for several numerical parameter combinations. We pass in all option contract and model parameters as well as all of the numerical parameters in the scheme.
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;
Now we include our function:
/*
ON INPUT:
S0 -- initial stock price
X -- exercise (strike) price
T -- Time to expiry (years)
r -- interest rate (per annum)
sigma -- volatility (per annum^12)
iMax -- number of time steps
jMax -- number of space steps
SMax -- Maximum value of S
omega -- is the relaxation parameter
tol -- is the tolerance level
iterMax -- is maximum iterations
ON OUTPUT:
return -- the value of a European put option at S=S0, t=0
*/
double europeanPut_CN(double S0,double X,double T,double r,double sigma,int iMax,int jMax,double SMax,double omega,double tol,int iterMax)
{
double dS=SMax/jMax;
double dt=T/iMax;
// create storage for the stock price and option price (old and new)
std::vector<double> S(jMax+1),vOld(jMax+1),vNew(jMax+1);
// setup and initialise the stock price
for(int j=0;j<=jMax;j++)
{
S[j] = j*dS;
}
// reset initial condition
for(int j=0;j<=jMax;j++)
{
vOld[j] = std::max(X-S[j],0.);
vNew[j] = std::max(X-S[j],0.);
}
// run through timesteps
for(int i=iMax-1;i>=0;i--)
{
// declare vectors for matrix equations
std::vector<double> a(jMax+1),b(jMax+1),c(jMax+1),d(jMax+1);
// set up matrix equations a[j]=
a[0] = 0.;b[0] = 1.;c[0] = 0.;
d[0] = X*exp(-r*(iMax-i)*dt);
for(int j=1;j<jMax;j++)
{
a[j] = 0.25*(sigma*sigma*j*j-r*j);
b[j] =-0.5*sigma*sigma*j*j - 0.5*r - 1./dt;
c[j] = 0.25*(sigma*sigma*j*j+r*j);
d[j] =-a[j]*vOld[j-1] - (b[j]+2./dt)*vOld[j] - c[j]*vOld[j+1];
}
a[jMax] = 0.;b[jMax] = 1.;c[jMax] = 0.;
d[jMax] = 0.;
// solve with SOR method
int sor;
for(sor=0;sor<iterMax;sor++)
{
double error=0.;
// implement sor in here
{
double y = (d[0] - c[0]*vNew[1])/b[0];
vNew[0] = vNew[0] + omega*(y-vNew[0]);
}
for(int j=1;j<jMax;j++)
{
double y = (d[j] - a[j]*vNew[j-1] - c[j]*vNew[j+1])/b[j];
vNew[j] = vNew[j] + omega*(y-vNew[j]);
}
{
double y = (d[jMax] - a[jMax]*vNew[jMax-1])/b[jMax];
vNew[jMax] = vNew[jMax] + omega*(y-vNew[jMax]);
}
// calculate residual norm ||r|| as sum of absolute values
error += fabs(d[0] - b[0]*vNew[0] - c[0]*vNew[1]);
for(int j=1;j<jMax;j++)
error += fabs(d[j] - a[j]*vNew[j-1] - b[j]*vNew[j] - c[j]*vNew[j+1]);
error += fabs(d[jMax] - a[jMax]*vNew[jMax-1] - b[jMax]*vNew[jMax]);
// make an exit condition when solution found
if(error<tol)
break;
}
if(sor>=iterMax)
{
std::cout << " Error NOT converging within required iterations\n";
std::cout.flush();
throw;
}
// set old=new
vOld=vNew;
}
int jStar=S0/dS;
double sum=0.;
sum+=(S0 - S[jStar])/dS * vNew[jStar+1];
sum+=(S[jStar+1] - S0)/dS * vNew[jStar];
return sum;
}
In order to complete our analysis we will need the analytic solutions of the put option to compare against. So we put these in as functions as well:
double normalDistribution(double x)
{
return 0.5*erfc(-x/sqrt(2.));
}
// return the value of a put option using the black scholes formula
double europeanPut_exact(double S,double X,double T,double r,double sigma)
{
if(fabs(T)<1.e-14) // check if we are at maturity
{
if(S<X)return X-S;
else return 0;
}
if((T)<=-1.e-14)return 0.; // option expired
if(X<1.e-14*S)return 0.; // check if strike << asset then exercise with certainty
if(S<1.e-14*X)return X*exp(-r*(T)) - S; // check if asset << strike then worthless
if(sigma*sigma*(T)<1.e-14) // check if variance very small then no diffusion
{
if(S<X*exp(-r*(T)))return X*exp(-r*(T)) - S;
else return 0.;
}
// calculate option price
double d1=(log(S/X) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
double d2=(log(S/X) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
return normalDistribution(-d2)*X*exp(-r*(T)) - normalDistribution(-d1)*S ;
}
Now declare our Black-Scholes parameters and the put option contract parameters, as well as the grid parameters;
// declare and initialise Black Scholes parameters
double S0,X,T,r,sigma;
// declare and initialise grid paramaters
int iMax,jMax;
double SMax;
and initialise them:
// initialise Black Scholes parameters
S0=1.973;X=2.;T=1.;r=0.05;sigma=0.4;
// initialise grid paramaters
iMax=40;jMax=40;SMax=2*X;
First let us check that the code is running correctly and the values look sensible, we expect the result from our Crank-Nicolson implementation to be broadly similar to the exact value.
cout << europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,SMax,1,1.e-6,1000);
cout << " " << europeanPut_exact(S0,X,T,r,sigma);
Now lets run for some different values of iMax and jMax...
cout << " n | V_cn | V_exact | error \n";
for(int k=1;k<=4;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value_exact = europeanPut_exact(S0,X,T,r,sigma);
double value_cn = europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,SMax,1,1.e-6,1000);
cout << n << " | " << value_cn;
cout << " | " << value_exact << " | " ;
cout << value_cn - value_exact << "\n";
}
If you are running this in the notebook you will notice that it takes quite a while to run, which is to be expected as the code is compiled on the fly. Running this code as a compiled executable (in VS or commandline) will be 100s or even 1000s of times faster.
We notice that the last result with n=160 is actually worse than the result with n=80, which does not seem to make sense. In fact this is caused by an inappropriate boundary condition, since S_\text{max}=2X is much too small it cause large truncation errors. So we try running the code again this time with S_\text{max}=5X.
cout << " n | V_cn | V_exact | error \n";
for(int k=1;k<=4;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value_exact = europeanPut_exact(S0,X,T,r,sigma);
double value_cn = europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,5*X,1.2,1.e-6,1000);
cout << n << " | " << value_cn;
cout << " | " << value_exact << " | " ;
cout << value_cn - value_exact << "\n";
}
This gives slightly better results than before, but we still see problems in the solution. One thing to note here is that the errors with n=160 are similar magnitude to the tolerance in our scheme. We should not expect the solution will be more accurate than the tolerance, so we increase the tolerance before increasing n to further investigate the solution.
cout << " n | V_cn | V_exact | error \n";
for(int k=1;k<=5;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value_exact = europeanPut_exact(S0,X,T,r,sigma);
double value_cn = europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,5*X,1.4,1.e-8,1000);
cout << n << " | " << value_cn;
cout << " | " << value_exact << " | " ;
cout << value_cn - value_exact << "\n";
}
So this doesn't look good -- the results haven't actually improved at all. So if it isn't S_\text{max} or the tolerance which is the problem, what can it be? There shouldn't be any non-linearity errors (grid is aligned to the strike price) so the only one left is the interpolation. To see if this is the problem, we can return the value on a grid node and see how that behaves. So set S_0=X and run again:
S0 = X;
cout << " n | V_cn | V_exact | error \n";
for(int k=1;k<=5;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value_exact = europeanPut_exact(S0,X,T,r,sigma);
double value_cn = europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,5*X,1.4,1.e-8,1000);
cout << n << " | " << value_cn;
cout << " | " << value_exact << " | " ;
cout << value_cn - value_exact << "\n";
}
In fact this gives very good results and we can very easily extrpolate these values. However, it also implies that the linear interpolation we use in the function is not good enough for extrapolation of values. We need our interpolation to several orders more accurate than the accuracy of the scheme. I recommend using at least cubic interpolation (left as an exercise). To examine how good this method can be if we have accurate results and stable convergence, consider the results above. Let a new estimate of the value be found using the formula (from Richardson extrapolation) V_\text{extrap} = \frac{4V_{2n} - V_n}{3}
To implement the Richardson extrapolation we introduce some extra storage to keep hold of the results from previous steps
double valueOld=0;
double valueOlder=0;
Now run the code again and this time include the extrapolated result.
S0 = X;
cout << " n | V_cn | V_exact | error | V_extrap | error \n";
for(int k=1;k<=5;k++)
{
int n=10*pow(2,k);
iMax = n;
jMax = n;
double value_exact = europeanPut_exact(S0,X,T,r,sigma);
double value_cn = europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,5*X,1.2,1.e-8,1000);
double value_extrap = (4.*value_cn - valueOld)/3.;
cout << n << " | " << value_cn;
cout << " | " << value_exact << " | " ;
cout << value_cn - value_exact << " | ";
cout << value_extrap << " | ";
cout << value_extrap - value_exact << "\n";
valueOld = value_cn;
}
So we can see that using just n=80 and n=160 we can get a solution which is accurate to six digits. In fact we can even extrapolate the extrapolated results to an even better solution. The formula would be V_{\text{extrap}^2} = \frac{16V_{4n} -8V_{2n} + V_n}{9} . This works for S=X, but for values that are not on grid nodes you would need an extremely accurate interpolation to make it work. The final solution is also on the web click here.