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.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 variance = 0.00609616
From simple statistics we know that each estimate of the solution V N is
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
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;
}
}
No comments:
Post a Comment