Sunday, February 22, 2015

C++ Coding - Black Scholes Option Pricing - Monte Carlo

3.2 Monte Carlo: Black Scholes European Call Option

Now we are going to value an European call option using Monte-Carlo. The setup is very simple, we just need to sum up the payoffs from a bunch of sample paths and then take the average. First start with an empty program except for the random number generator, as follows

#include <iostream>
#include <random>
#include <cmath>
using namespace std;

int main()
{
  // declare the random number generator
  mt19937 rng;
  
  
  
}  

where we are including the random library, as well as cmath for access to mathematical functions required later. Next declare the variable which contains the methods to get random draws from a normal distribution

  // declare the distribution
  normal_distribution<> ND(0.,1.);

Now we have everything we need to start, declare your parameters for the option at the top of your main function alongside the number of simulations we are going to run (a counting number)

  // parameters
  double stock0=9.576,strikePrice=10.,interestRate=0.05,sigma=0.4,maturity=0.75;
  // number of paths
  int N=1000;
  

We put these at the top to indicate that they will be input variables when we move this out to its own function later on. Now we can write the Monte-Carlo algorithm in less than 10 lines, start by initialising a sum variable with zero and then run through each path, adding in the payoff of a call for each particular stock price at maturity

  // initialise sum to zero
  double sum=0.;
  // run through each simulation
  for(int i=0;i<N;i++)
  {
    double phi=ND(rng); // get a random draw from N(0,1)
    double ST=S0 * exp( (interestRate - 0.5*sigma*sigma)*maturity + phi*sigma*sqrt(maturity) ); // get a random draw for the value of stock price at maturity
    sum = sum + max( ST - strikePrice , 0. ); // add in the payoff of the option in that case
  }
  cout << sum/N*exp(-interestRate*maturity) << endl; // output the average value over all paths

Running this code should return a value of 1.2671, which we can check is reasonably close to the analytic value we would expect.

In order to do any sort of meaningful analysis, we must now put this algorithm into a function. Use the following definition of the function

double monteCarlo(double S0,double strikePrice,double interestRate,double sigma,double maturity,int N)
{
  
}

and copy/paste everything from main into this function (except the parameter declarations). Our full program at this stage should now look like

#include <iostream>
#include <random>
#include <cmath>
using namespace std;

double monteCarlo(double S0,double strikePrice,double interestRate,double sigma,double maturity,int N)
{
  // declare the random number generator
  mt19937 rng;
  // declare the distribution
  normal_distribution<> ND(0.,1.);
  ND(rng);
  // initialise sum
  double sum=0.;
  for(int i=0;i<N;i++)
  {
    double phi=ND(rng);
    double ST=S0 * exp( (interestRate - 0.5*sigma*sigma)*maturity + phi*sigma*sqrt(maturity) );
    sum = sum + max( ST - strikePrice , 0. );
  }
  return sum/N*exp(-interestRate*maturity);
}

int main()
{
  cout << monteCarlo(9.576,10.,0.05,0.4,0.75,1000) << endl;
}

Run it and check you get the same result as last time. Now to do some analysis, we are going to want to run this algorithm lots of times with a different number of paths. So try running the function several times in main (copy/paste the single line of code)

  cout << monteCarlo(9.576,10.,0.05,0.4,0.75,1000) << endl;
  cout << monteCarlo(9.576,10.,0.05,0.4,0.75,1000) << endl;
  cout << monteCarlo(9.576,10.,0.05,0.4,0.75,1000) << endl;
  cout << monteCarlo(9.576,10.,0.05,0.4,0.75,1000) << endl;
  // output is
  // 1.2671
  // 1.2671
  // 1.2671
  // 1.2671

The output is the same each time, but why? We are supposed to be generating random numbers. Unfortunately because of the way we declare the random number generator, each time we call this function a new version of the random number generator is created which gets reset back to the first number in the sequence. To avoid this we need to make sure that one and only one random number generator is ever created for this function, so we use the keyword static to tell the compiler this. The appropriate line in your code is changed to

  static mt19937 rng;

and when we run the program again we get

1.2671  
1.33243  
1.22289  
1.31114

Ok now we are confident that we can generate random samples of the solution for a given N, we want to see what happens expected distribution of the results as we increase the number of paths. We are going to do some analysis now which requires the stl storage container vector which allows for dynamic allocation of memory. They can perform all of the same things a simple array can do and much more. In the main function, we declare vector<double> samples(M) which is an array of type double intialised with M values. First check you can run the results by outputting to screen

  // we want to run M calculations
  int M=100;
  // now store all the results
  vector<double> samples(M);
  // number of paths in each calculation
  int N=1000;
  // run some calculations
  for(int i=0;i<M;i++)
  {
    cout << monteCarlo(9.576,10.,0.05,0.4,0.75,N) << endl;
  }

Now rather than output to screen write them into the vector samples using the array syntax like this

// run some calculations
  for(int i=0;i<M;i++)
  {
    // cout << monteCarlo(9.576,10.,0.05,0.4,0.75,N) << endl;
    samples[i] =  monteCarlo(9.576,10.,0.05,0.4,0.75,N);
  }

Calculate the mean of the vector

  double sum=0.;
  for(int i=0;i<M;i++)
  {
    sum+=samples[i];
  }
  double mean = sum/M;
  cout << " sample mean = " << mean << endl;

and variance

  double sumvar=0.;
  for(int i=0;i<M;i++)
  {
    sumvar+=(samples[i]-mean)*(samples[i]-mean);
  }
  double variance = sumvar/(M-1);
  cout << " sample variance = " << variance << endl; 
  

and the output is

 sample mean = 1.28323  
 sample variance = 0.00609616

From simple statistics we know that each estimate of the solution V N is

V   ∼ N (V ∗,σ2)
 N
where V is the true solution and σ2 is the variance. The if we choose to take a sample mean from that distribution V with M samples we have
        (       )
 ¯         ∗ σ2-
V  ∼ N   V  ,M
So to get a confidence interval for our result we use this result to output the following

  double sd = sqrt(variance/M);
  cout << " 95% confident result is in ["<<mean-2.*sd << "," << mean+2.*sd << "] with " << N*M << " total paths. " << endl;

and the output is

 95% confident result is in [1.26761,1.29884] with 100000 total paths.

We can now run this for a range of values for N and M. Note first that the sample mean estimate uses NM total paths to get the confidence interval. We could set N = 1 and M = 100000 and you should get a similar result to the one above. This is because even with a single path the European put/call option will be approximately normally distributed. For other more complex options this might not be the case and you will need to ensure both N >> 1 and M >> 1. Finally we present some code to run different values of M for fixed N on a European call option. Note here that I have included a payoff function, so that if you need to solve for a different option you should change the payoff here and it should work fine.

#include <iostream>
#include <random>
#include <cmath>
using namespace std;

// payoff from the European call option
double payoff(double S,double strikePrice)
{
  return max( S - strikePrice , 0. ); // change this line here to solve for different European options
}

double monteCarlo(double S0,double strikePrice,double interestRate,double sigma,double maturity,int N)
{
  // declare the random number generator
  static mt19937 rng;
  // declare the distribution
  normal_distribution<> ND(0.,1.);
  ND(rng);
  // initialise sum
  double sum=0.;
  for(int i=0;i<N;i++)
  {
    double phi=ND(rng);
    // calculate stock price at T
    double ST=S0 * exp( (interestRate - 0.5*sigma*sigma)*maturity + phi*sigma*sqrt(maturity) );
    // add in payoff
    sum = sum + payoff(ST,strikePrice);
  }
  // return discounted value
  return sum/N*exp(-interestRate*maturity);
}

int main()
{
  // run for different 
  for(int M=100;M<=100000;M*=10)
  {
  // now store all the results
  vector<double> samples(M);
  // number of paths in each calculation
  int N=1000;

  cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
  cout << " Run results with M="<<M<<" samples from V_N, where N="<<N<<"."<<endl;
  
  // run some calculations
  for(int i=0;i<M;i++)
  {
    samples[i] = monteCarlo(9.576,10.,0.05,0.4,0.75,N);
  }
  // estimate the mean from the sample
  double sum=0.;
  for(int i=0;i<M;i++)
  {
    sum+=samples[i];
  }
  double mean = sum/M;
  cout << " mean = " << mean << endl;

  // estimate the variance from the sample
  double sumvar=0.;
  for(int i=0;i<M;i++)
  {
    sumvar+=(samples[i]-mean)*(samples[i]-mean);
  }
  double variance = sumvar/(M-1);
  cout << " variance = " << variance << endl; 
  
  // get the standard deviation of the sample mean
  double sd = sqrt(variance/M);
  cout << " 95% confident result is in ["<<mean-2.*sd << "," << mean+2.*sd << "] with "<< N*M << " total paths." << endl;
  }
}

1 comment:

  1. Varieties of languages offered - Being more convenient with use of the binary trading software program, it is essential how the language of the program should be your own mother language. Binary Options There are some expert binary options broker platforms on the market that offer varieties associated with languages from English, Persia, Spanish, Russian, Turkish, and so on.

    ReplyDelete