C++ Coding - Black Scholes Option Pricing - Explicit Finite Difference

The example question for these solutions can be found on my website (click here).

6.1 Explicit Finite Difference For Option Pricing

In this example we are going to price a European call option with explicit finite difference.

Initial Setup

For the explicit method we shall need:

• All parameters for the option, such as X and S₀ etc.
• The number of divisions in stock, jMax, and divisions in time iMax
• The size of the divisions ΔS = S_max∕jMax and Δt = T∕iMax
• A vector to store the stock prices, and two to store the option values at the current and previous time level.

It will make things easier if you just declare these at the top of the program before you start doing anything with them.

The program structure will look something like this:

Use the code below to as a template to get started.

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace std;

/* Template code for the Explicit Finite Difference
 */
int main()
{
  // declare and initialise Black Scholes parameters
  
  // declare and initialise grid paramaters 
  
  // declare and initialise local variables (ds,dt)
  
  // create storage for the stock price and option price (old and new)
  
  // setup and initialise the stock price 
  
  // setup and initialise the final conditions on the option price 
  
  // loop through time levels, setting the option price at each grid point, and also on the boundaries
  
  // output the estimated option price
  
}

download

Inside the template the comments list the rough layout of your program. First declare the parameters for the problem, local grid variables and the vectors required. Set the values of the parameters to the following: σ = 0.4, r = 0.05, X = 2, T = 1, and iMax = jMax = 4. Next we need to note that in option pricing problems we must for S over the semi inifinite domain, therefore numerically we need to choose an appropriate S_max. In general, you should set this according to mathematical reasoning if you want your code to be generic and work in all cases. Often we just choose S_max = K × X where K is some constant (say 5) but it should really be linked to the probability distribution of the underlying asset. For simplicity here we choose S_max = 2X so that small grids still work, for any decent level of accuracy S_max will need to set much higher. Next assuming we have set the correct level of storage to your vector, calculate the values of ΔS and Δt for your grid and using them to assign initial values to the stock price vector, and the option value vectors.

It is a good idea at this stage to output the stock values and option values to check that the results are as you would expect them to be. The lines added to our template might look something like this

  // declare and initialise Black Scholes parameters
  double S0=100.,X=100.,T=1.,r=0.06,sigma=0.2;
  // 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(S[j]-X,0.);
    vNew[j] = max(S[j]-X,0.);
    cout << iMax << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl;
  }
/* OUTPUT
4 0 0 0 0
4 1 1 0 0
4 2 2 0 0
4 3 3 1 1
4 4 4 2 2
*/  

download

Timestep Calculation

Next we must setup a loop to count backwards through time, and at each timestep another loop to go through all stock price values (except the boundaries). We wish to use two time levels of storage for V , which we shall call vNew and vOld. Here let V ⁱ be represented by the vector vNew, and V ⁱ⁺¹ be represented by the vector vOld. Then at the end of each timestep we must overwrite the old values with the new ones, please refer to the solution at the end to see this in action. This allows us to move recursively through time with only two levels of storage. Since we can often have timestep constraints, it can be extremely inefficient to store all levels of time.

Firstly add in your code the timestep loop over i = iMax − 1 down through to i = 0 to step backwards through time. Now inside the loop, you must first implement the boundary condition at j = 0. For example, for a call we may write

and remember that vNew= V ⁱ. Next add in the condition at j = jMax, for a call we may write

Check at this stage that you are implementing the boundary conditions correctly, so output them to screen. Your code for the loop may look like

  // loop through time levels, setting the option price at each grid point, and also on the boundaries
  for(int i=iMax-1;i>=0;i--)
  {
    // apply boundary condition at S=0
    vNew[0] = 0.;
    cout << i << " " << 0 << " " << S[0] << " " << vNew[0] << " " << vOld[0] << endl;
    // apply boundary condition at S=S_max
    vNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt));
    cout << i << " " << jMax << " " << S[jMax] << " " << vNew[jMax] << " " << vOld[jMax] << endl;
  }
  

download

Now because we are only using five nodes in space and five nodes in time, this will allow us to check our calculations by hand. Make sure you are checking the values outputed from your code at every step.

Next we can put a loop in between the boundary conditions to calculate for the values 1 ≤ j ≤ jMax − 1. Inside the loop, to make things easier follow the notation from the lectures (click here) and declare variables to store the coefficients A, B and C. Then simply calculate the values of A, B and C for the current value of j, and write

Again output the values of the code and check against the solution in figure 1. You must remember to update the old values to new by putting a line vOld=vNew; at the end of the time loop.

Once all code is checked and working you should run grid checks (change the grid size and check the effect on the solution) to satisfy yourself that the code is correct. You will notice that if you increase jMax by a factor of 10 you will need to increase iMax by a factor of 100 to maintain stability. Please refer to the course notes and lectures for more details on this (click here).

---

Figure 1: Solution for given parameters.

---

The following code generates the results found in figure 1 and at the end you have estimates for the option price at time t = 0 for a discrete set of values in S. The value estimates might not though correspond to the value of S₀ you need the option price at, so we will need to look at interpolation in the next section.

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace std;

/* Template code for the Explicit Finite Difference
 */
int main()
{
  // declare and initialise Black Scholes parameters
  double S0=1.,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(S[j]-X,0.);
    vNew[j] = max(S[j]-X,0.);
    cout << iMax << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl;
  }
  // loop through time levels, setting the option price at each grid point, and also on the boundaries
  for(int i=iMax-1;i>=0;i--)
  {
    // apply boundary condition at S=0
    vNew[0] = 0.;
    cout << i << " " << 0 << " " << S[0] << " " << vNew[0] << " " << vOld[0] << endl;
    for(int j=1;j<=jMax-1;j++)
    {
      double A,B,C;
      A=0.5*sigma*sigma*j*j*dt+0.5*r*j*dt;
      B=1.-sigma*sigma*j*j*dt;
      C=0.5*sigma*sigma*j*j*dt-0.5*r*j*dt;
      vNew[j] = 1./(1.+r*dt)*(A*vOld[j+1]+B*vOld[j]+C*vOld[j-1]);
      cout << i << " " << j << " " << S[j] << " " << vNew[j] << " " << vOld[j] << endl;
    }
    // apply boundary condition at S=S_max
    vNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt));
    cout << i << " " << jMax << " " << S[jMax] << " " << vNew[jMax] << " " << vOld[jMax] << endl;
    // set old values to new
    vOld=vNew;    
  }
}

download

Interpolation

Example - Linear Interpolation

We show here a simple example of how to take the code detailed above and output a value for the option at a given stock price using linear interpolation. We use a style here which allows for a more generic higher order Lagrange polynomial approximation to be derived. We assume that we have the set of points S_j and vNew_j and we wish to find the value of the interpolated function V (S). The linear approximation says that for points j∗ and j ∗ +1 we should choose

There are three stages:

• Given the value S₀, find j^∗ such that S ₀ ∈ [j ∗ ΔS, (j ∗ +1)ΔS]
  

  int jstar;
  jstar = S0/dS;

  download
• Evaluate the Lagrange polynomial at S = S₀
  

    double sum=0.;
    sum = sum + (S0 - S[jstar+1])/(S[jstar]-S[jstar+1])*vNew[jstar];
    sum = sum + (S0 - S[jstar])/(S[jstar+1]-S[jstar])*vNew[jstar+1];

  download
• Test and debug the code!!! You can check that the polynomial above is the equation of a line passing through the points (S_j^∗,vNew_j^∗) and (S_j^∗+1,vNew_j^∗+1). Try different values of S₀ and check they look ok.

Once you are convinced the code is working move everything out into a function. You will now be able to pass in S₀ as an argument to the function, simply return the value from the interpolation.

The final code might look like this

#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
using namespace std;

/* Template code for the Explicit Finite Difference
 */
double explicitCallOption(double S0,double X,double T,double r,double sigma,int iMax,int jMax)
{
  // 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(S[j]-X,0.);
    vNew[j] = max(S[j]-X,0.);
  }
  // loop through time levels, setting the option price at each grid point, and also on the boundaries
  for(int i=iMax-1;i>=0;i--)
  {
    // apply boundary condition at S=0
    vNew[0] = 0.;
    for(int j=1;j<=jMax-1;j++)
    {
      double A,B,C;
      A=0.5*sigma*sigma*j*j*dt+0.5*r*j*dt;
      B=1.-sigma*sigma*j*j*dt;
      C=0.5*sigma*sigma*j*j*dt-0.5*r*j*dt;
      vNew[j] = 1./(1.+r*dt)*(A*vOld[j+1]+B*vOld[j]+C*vOld[j-1]);
    }
    // apply boundary condition at S=S_max
    vNew[jMax] = S[jMax] - X*exp(-r*(T-i*dt));
    // set old values to new
    vOld=vNew;    
  }
  // get j* such that S_0 \in [ j*dS , (j*+1)dS ]
  int jstar;
  jstar = S0/dS;
  double sum=0.;
  // run 2 point Lagrange polynomial interpolation
  sum = sum + (S0 - S[jstar+1])/(S[jstar]-S[jstar+1])*vNew[jstar];
  sum = sum + (S0 - S[jstar])/(S[jstar+1]-S[jstar])*vNew[jstar+1];
  return sum;
}

int main()
{
  // declare and initialise Black Scholes parameters
  double S0=1.639,X=2.,T=1.,r=0.05,sigma=0.4;
  // declare and initialise grid paramaters 
  int iMax=4,jMax=4;
  cout << explicitCallOption(S0,X,T,r,sigma,iMax,jMax) << endl;
  /* OUTPUT
0.194858
   */
}

download