#### 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>usingnamespacestd; intmain() {// declare the random number generatormt19937 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 distributionnormal_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)

// parametersdouble stock0=9.576,strikePrice=10.,interestRate=0.05,sigma=0.4,maturity=0.75;// number of pathsint 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 zerodouble sum=0.;// run through each simulationfor(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 maturitysum = 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

doublemonteCarlo(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>usingnamespacestd; doublemonteCarlo(double S0,double strikePrice,double interestRate,double sigma,double maturity,int N) {// declare the random number generatormt19937 rng;// declare the distributionnormal_distribution<>ND(0.,1.);ND(rng);// initialise sumdouble 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. ); }returnsum/N*exp(-interestRate*maturity); } intmain() { 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

staticmt19937 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 calculationsint M=100;// now store all the resultsvector<double>samples(M);// number of paths in each calculationint N=1000;// run some calculationsfor(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 calculationsfor(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

^{∗}is the true solution and σ

^{2}is the variance. The if we choose to take a sample mean from that distribution with M samples we have

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>usingnamespacestd;// payoff from the European call optiondoublepayoff(double S,double strikePrice) {returnmax( S - strikePrice , 0. );// change this line here to solve for different European options} doublemonteCarlo(double S0,double strikePrice,double interestRate,double sigma,double maturity,int N) {// declare the random number generatorstaticmt19937 rng;// declare the distributionnormal_distribution<>ND(0.,1.);ND(rng);// initialise sumdouble sum=0.;for(int i=0;i<N;i++) { double phi=ND(rng);// calculate stock price at Tdouble ST=S0 *exp( (interestRate - 0.5*sigma*sigma)*maturity + phi*sigma*sqrt(maturity) );// add in payoffsum = sum +payoff(ST,strikePrice); }// return discounted valuereturnsum/N*exp(-interestRate*maturity); } intmain() {// run for differentfor(int M=100;M<=100000;M*=10) {// now store all the resultsvector<double>samples(M);// number of paths in each calculationint N=1000; cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl; cout << " Run results with M="<<M<<" samples from V_N, where N="<<N<<"."<<endl;// run some calculationsfor(int i=0;i<M;i++) { samples[i] =monteCarlo(9.576,10.,0.05,0.4,0.75,N); }// estimate the mean from the sampledouble 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 sampledouble 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 meandouble 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