#### 4.1 Time To First Exit

Calculate the expected hitting time E[t^{∗}] for a Brownian motion X(t), where the process must
hit either X(t^{∗}) = 0 or X(t^{∗}) = 1 for t^{∗} < T. If neither boundary is hit within the time T then
t^{∗} = T.

Assume that the process X follows the SDE

Now we need to consider the problem, which is to track the path of Brownian motion through time, and check at each point in time whether or not the process has crossed either one or zero. So there are obvious stages in which we can build up our code, namely

- Stage 1: Generate random numbers from a normal distribution
- Stage 2: Generate a single path in X
- Stage 3: Stop the path when X < 0 or X > 1 and output the current time
- Stage 5: Calculate expected value
- Stage 6: Move into function returning expected value for N runs
- Stage 7: Test and validate the code

Stage 1 has already been dealt with in a previous example. You can use this sample code to start.

#include<iostream>#include<random>#include<cmath>usingnamespacestd; intmain() {// declare the random number generatormt19937 rng;// declare the distributionnormal_distribution<>ND(0.,1.); }

Now we need to think about variables and parameters for this problem. The drift, standard deviation, maximum time T and the initial value of X are the input parameters, as well as possibly the position of the bounds. For a given path we also need to specify the discretisation size dt. Declare all of these variables and assign them with some trivial values to begin with.

// declare parametersdouble X0=0.56,mu=0.01,sigma=0.75,T=1.;// grid paramsint K=100;// local paramsdouble dt=T/K;

Compile and run the program to check everything is ok so far. Next we need to keep track of the current value of X(t) independently of storing the initial value which will become clearer later on. Declare another variable Xt to store X(t), initialise it with an appropriate value. Output to screen with the command cout << Xt << endl;, you should get 0.56. Now we want to calculate the path using a recurrence relation for X. We choose to use a for loop for this, which must run between k = 0 and k = K − 1. Input the equation for the next value of Xt overwriting the previous value. It should look something like this

// first generate a path in X(t) and print to screen]// initialise value of X at time 0double t=0.; double Xt=X0; cout << " t \t| X(t) "<<endl; cout << "-------------------- "<<endl; cout <<t<< " \t| " << Xt << endl;for(int k=0;k<K;k++) { Xt = Xt + mu*dt +sqrt(dt)*sigma*ND(rng); t = (k+1)*dt; cout <<t<< " \t| " << Xt << endl; }

The final value should t = 1 and X(1) = 0.412219. Inspecting the path we can see that it hits a value of X(t) = 1.06307 at t = 0.54.

Now we want to implement that inspection of the path within the code. To make the code more flexible for future modifications and make out code look more like a standard Monte Carlo algorithm we are going to first reimplement the path generation to store all the values in an array. We choose to use a stl vector to store the values, again for flexibility. First step is to delete the declaration for Xt and redeclare it vector<double> Xt(K+1);, which is a vector with K + 1 elements. Now we initialise the first value in the vector to Xt[0]=X0; and then use appropriate values in the recurrence relation. This part of your code should now look like

// now generate the path in a vector X(K+1) and print to screenvector<double>Xt(K+1);// initialise value of X at time 0double t=0.; Xt[0]=X0; cout << " t \t| X(t) "<<endl; cout << "-------------------- "<<endl; cout <<t<< " \t| " << Xt[0] << endl;for(int k=0;k<K;k++) { Xt[k+1] = Xt[k] + mu*dt +sqrt(dt)*sigma*ND(rng); t = (k+1)*dt; cout <<t<< " \t| " << Xt[k+1] << endl; }

Now we want to implement the payoff, defined as the first hitting time. To do this we seperate the generation of the path from evaluating the payoff. So once the path has been generated, we move back through the path and calculate the payoff. For each element in the vector Xt, we check if it is above or below the bounds. If either condition is satisfied we break from the loop and return the current time. If we get to the final element just return T. My implementation of this looks like this

// now calculate first exit time from the pathdouble exitTime=0.;for(int k=0;k<=K;k++) { t=k*dt;// output to screen to check what is going oncout <<t<< " \t| " << Xt[k] << endl;//if(k==K)exitTime = t;if(Xt[k]<0.)// set exittime and stop{ exitTime = t;break; }if(Xt[k]>1.)// set exittime and stop{ exitTime = t;break; } }// finished getting exit timecout << " Exit time " << exitTime << endl;

Now we have generated a single path we check and run the code. For the default seed (none
supplied), we get t^{∗}=0.54 as expected. We can run the code for another path by resetting the
seed (place rng.seed(1) directly underneath where it is declared). By changing the value of the
seed you should see a range of values are possible. Once we are happy that this is behaving
correctly add a Monte Carlo loop to calculate the average exit time over N samples. Underneath
where you have declared the normal distribution declare N, your sum variable and start the
Monte Carlo loop

// declare the distributionnormal_distribution<>ND(0.,1.);// add in the monte carlo loopint N=1000; double sum=0.;for(int i=0;i<N;i++) {// now generate the path in a vector X(M+1) and print to screen

close the loop towards the bottom of your program and output the average exitTime

}// finished getting exit time//cout << " Exit time " << exitTime << endl;sum = sum + exitTime; }// finish Monte Carlo loopcout << " Expected exit time " << sum/N <<endl;

Again we should test now with different values of N and K to ensure that things are working properly. After I have removed unnecessary screen outputs my code looks something like this

#include<iostream>#include<random>#include<cmath>usingnamespacestd; intmain() {// declare parametersdouble X0=0.56,mu=0.01,sigma=0.75,T=1.;// grid paramsint M=100;// local paramsdouble dt=T/K;// declare the random number generatormt19937 rng;// declare the distributionnormal_distribution<>ND(0.,1.);// monte carlo loopint N=1000; double sum=0.;for(int i=0;i<N;i++) {// now generate the path in a vector X(K+1) and print to screenvector<double>Xt(K+1);// initialise value of X at time 0double t=0.; Xt[0]=X0;for(int k=0;k<K;k++) { Xt[k+1] = Xt[k] + mu*dt +sqrt(dt)*sigma*ND(rng); }// now calculate first exit time from the pathdouble exitTime=0.;for(int k=0;k<=K;k++) {// get current time this loopt=k*dt;if(k==K)exitTime = t;if(Xt[k]<0.)// set exittime and stop{ exitTime = t;break; }if(Xt[k]>1.)// set exittime and stop{ exitTime = t;break; } }// finished getting exit timesum = sum + exitTime; }// finish Monte Carlo loopcout << " Expected exit time " << sum/N <<endl; }

Once we are convinced we can put the Monte Carlo part in its own function. Remember that once it goes inside the function we need to make the random number generator static so that the sequence is not reset every time the function is called. We can now analyse the results, taking different value of N and K and trying to get confidence interval just like we did in a previous example. Finally your code might look something like:

#include<iostream>#include<random>#include<cmath>usingnamespacestd; doublemonteCarlo(double X0,double mu,double sigma,double T,int K,int N) {// local paramsdouble dt=T/K;// declare the random number generatorstaticmt19937 rng;// declare the distributionnormal_distribution<>ND(0.,1.); double sum=0.;for(int i=0;i<N;i++) {// now generate the path in a vector X(K+1) and print to screenvector<double>Xt(K+1);// initialise value of X at time 0double t=0.; Xt[0]=X0;for(int k=0;k<K;k++) { Xt[k+1] = Xt[k] + mu*dt +sqrt(dt)*sigma*ND(rng); }// now calculate first exit time from the pathdouble exitTime=0.;for(int k=0;k<=K;k++) {// get current time this loopt=k*dt;if(k==K)exitTime = t;if(Xt[k]<0.)// set exittime and stop{ exitTime = t;break; }if(Xt[k]>1.)// set exittime and stop{ exitTime = t;break; } }// finished getting exit timesum = sum + exitTime; }// finish Monte Carlo loopreturnsum/N; } intmain() {// declare parametersdouble X0=0.56,mu=0.01,sigma=0.75,T=1.;// grid paramsint K=1000; cout << " Solve for the first exit time with parameters " << endl; cout << " X0="<< X0 << ", mu="<< mu <<", sigma="<< sigma << ", T="<< T << " and K="<< K << "." << endl; cout << " ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;// run for differentfor(int M=10;M<=1000;M*=10) {// now store all the resultsvector<double>samples(M);// monte carlo loopint 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(X0,mu,sigma,T,K,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