Processing math: 21%

Analysing the Crank-Nicolson Method

Click here for the interactive version of this notebook. Binder

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 Pt+12σ2S22PS2+rSPSrP=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.

In [1]:
#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
using namespace std;

Now we include our function:

In [3]:
/* 
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:

In [4]:
double normalDistribution(double x)
{
  return 0.5*erfc(-x/sqrt(2.));
}
In [5]:
// 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;

In [6]:
// 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:

In [7]:
// 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.

In [8]:
cout << europeanPut_CN(S0,X,T,r,sigma,iMax,jMax,SMax,1,1.e-6,1000);
cout << " " << europeanPut_exact(S0,X,T,r,sigma);
0.273017 0.273152

Now lets run for some different values of iMax and jMax...

In [9]:
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";
}
 n    | V_cn     | V_exact  | error   
20 | 0.271973 | 0.273152 | -0.00117863
40 | 0.273017 | 0.273152 | -0.000134771
80 | 0.273116 | 0.273152 | -3.58068e-05
160 | 0.273085 | 0.273152 | -6.66597e-05

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.

In [10]:
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";
}
 n    | V_cn     | V_exact  | error   
20 | 0.26121 | 0.273152 | -0.0119412
40 | 0.271037 | 0.273152 | -0.00211423
80 | 0.2729 | 0.273152 | -0.000251823
160 | 0.273158 | 0.273152 | 6.52624e-06

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.

In [13]:
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";
}
 n    | V_cn     | V_exact  | error   
20 | 0.26121 | 0.273152 | -0.0119412
40 | 0.271037 | 0.273152 | -0.00211422
80 | 0.2729 | 0.273152 | -0.000251822
160 | 0.273158 | 0.273152 | 6.52653e-06
320 | 0.273123 | 0.273152 | -2.89602e-05

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:

In [14]:
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";
}
 n    | V_cn     | V_exact  | error   
20 | 0.247269 | 0.262918 | -0.0156485
40 | 0.259262 | 0.262918 | -0.00365595
80 | 0.262016 | 0.262918 | -0.000902232
160 | 0.262693 | 0.262918 | -0.00022486
320 | 0.262862 | 0.262918 | -5.61718e-05

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

In [15]:
double valueOld=0;
double valueOlder=0;

Now run the code again and this time include the extrapolated result.

In [16]:
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;
}
 n    | V_cn     | V_exact  | error  | V_extrap | error 
20 | 0.247269 | 0.262918 | -0.0156485 | 0.329693 | 0.0667747
40 | 0.259262 | 0.262918 | -0.00365595 | 0.263259 | 0.000341548
80 | 0.262016 | 0.262918 | -0.000902232 | 0.262934 | 1.56755e-05
160 | 0.262693 | 0.262918 | -0.00022486 | 0.262919 | 9.31133e-07
320 | 0.262862 | 0.262918 | -5.61718e-05 | 0.262918 | 5.74938e-08

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.