For the Crank-Nicolson method we shall need:
Again, the easiest thing to do is just to declare these at the top of the program before you start doing anything with them.
A sample program structure is given in the example code click here and the final solution is click here. I have initialised the variables and put the loop in to save time. At each timestep we must solve the matrix equation AV=b. Since the matrix is tridiagonal we need only use three vectors to store all values in the matrix A.
Now, assuming that our vectors are indexed from 0 to iMax (and therefore have iMax+1 elements) I will go through an example for a European put with σ=0.4, r=0.05, X=2, dS=1, dt=0.25, and jMax=4. Note here that these calculations could have been done by hand, and trying to replicate an example of this scale with a code is a good way to start, and also to check for bugs/errors.
#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace std;
// declare and initialise Black Scholes parameters
double S0=1.974,X=2.,T=1.,r=0.05,sigma=0.4;
// declare and initialise grid paramaters
int iMax=4,jMax=4;
// declare and initialise local variables (ds,dt)
double S_max=2*X;
double dS=S_max/jMax;
double dt=T/iMax;
// create storage for the stock price and option price (old and new)
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;
}
// setup and initialise the final conditions on the option price
for(int j=0;j<=jMax;j++)
{
vOld[j] = max(X-S[j],0.);
vNew[j] = max(X-S[j],0.);
cout << iMax << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl;
}
// start looping through time levels
for(int i=iMax-1;i>=0;i--)
{
// declare vectors for matrix equations
// set up matrix equations a[j]=
// solve matrix equations with SOR
// output results
for(int j=0;j<=jMax;j++)
{
cout << i << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl;
}
// set old=new
}
// finish looping through time levels
// output the estimated option price
I tend to use a, b and c to represent the tridagonal matrix A, and
d for the right hand side of the equation. Just remember that it makes
things easier if you are consistent with your notation.
a, b, c and d, and initialise them to be jMax+1 in size.for(int i=iMax-1;i>=0;i--)
{
// declare vectors for matrix equations
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.;
// output matrix
for(int j=0;j<=jMax;j++)
{
cout << j << " " << a[j] << " " << b[j] << " " << c[j] << " " << d[j] << endl;
}
// just check output at first step
break;
}
Still inside the time loop, we now need to solve the matrix equation. This can be done using either iterative methods (SOR) or direct methods (Thomas's solver). We shall go through the SOR method here as it is more easily adapted to American options.
The method is as follows:-
Create a loop in your code after the matrix has been setup, to iterate with SOR until a solution is found.
// store the value of sor OUTSIDE the loop to check on convergence
int sor;
for(sor=0;sor<iterMax;sor++)
{
//implement sor in here...
}
where iterMax is the maximum number of iterations you permit (may need to be as high as 10000). Like the simple example from before, you need to write an equation to update the value at the boundary j=0, then a loop of equations to update all values in between 1≤j<jMax and then another at the boundary j=jMax.
// reset initial condition
for(int j=0;j<=jMax;j++)
{
vOld[j] = max(X-S[j],0.);
vNew[j] = max(X-S[j],0.);
}
// run through timesteps
for(int i=iMax-1;i>=0;i--)
{
cout << "## Solving at timestep "<< i << "\n";
// declare vectors for matrix equations
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,iterMax=10;
double tol = 1.e-8,omega=1.;
for(sor=0;sor<iterMax;sor++)
{
double error=0.;
cout << sor << " == \n";
// implement sor in here
{
double y = (d[0] - c[0]*vNew[1])/b[0];
vNew[0] = vNew[0] + omega*(y-vNew[0]);
cout << " ( " << 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]);
cout << vNew[j] << " , ";
}
{
double y = (d[jMax] - a[jMax]*vNew[jMax-1])/b[jMax];
vNew[jMax] = vNew[jMax] + omega*(y-vNew[jMax]);
cout << vNew[jMax] << " )\n";
}
// 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;
}
cout << " Solved after " << sor << " iterations.\n";
// set old=new
vOld=vNew;
}
Once we have a solution in terms of S and V at particular points, we can get the value of the option at any value of S using interpolation. So to get the value at S0=1.974 we do:
double optionValue;
{
int jStar=S0/dS;
double sum=0.;
sum+=(S0 - S[jStar])/dS * vNew[jStar+1];
sum+=(S[jStar+1] - S0)/dS * vNew[jStar];
optionValue = sum;
}
cout << "Value of the option V(S="<<S0<<") = " << optionValue << "\n";
and if we want the value at another value (say S0=2.141) we have
S0=2.141;
{
int jStar=S0/dS;
double sum=0.;
sum+=(S0 - S[jStar])/dS * vNew[jStar+1];
sum+=(S[jStar+1] - S0)/dS * vNew[jStar];
optionValue = sum;
}
cout << "Value of the option V(S="<<S0<<") = " << optionValue << "\n";
Finally, you want to wrap all this code up into a function so that it can be used again and again.