C++ Coding - Time to first exit

Question Calculate the expected time to first exit from the interval [0,1] for a Brownian motion X with the following SDE

where dW is a standard Wiener process and X(t=0)=0.56.

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 4: Move into function returning exit time for one path
• 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 post, which you can find here. Now we need to think about variables and parameters for this problem. The standard deviation and the initial value of X are the only 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.

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double uniformRandom()
{
  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}
 // return a normally distributed random number
double normalRandom()
{
  double u1=uniformRandom();
  double u2=uniformRandom();
  return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)); 
}

int main()
{
  // declare input parameters
  double X0=0.56,sigma=0.4;
  // bounds
  double a=0.,b=1.;
  // discretisation variable
  double dt=0.1;
}

Compile and run the program to check everything is ok so far. Next we need to keep track of the current value of X independently of storing the initial value which will become clearer later on. Declare another variable to store X and output the first step in the path to screen. You should get t=0 X_t=0.56 on the screen.

  // store the current value of X
  double Xt=X0;
  // plot out the first value in the path
  cout << "t=" << 0 << " X_t=" << Xt << endl;

Now we want to calculate the path using a recurrence for X. We choose to use a for loop for this and set a maximum number of iterations, in this case 20.

  // now plot out a path in time
  for(int i=1;i<=20;i++)
  {
    Xt = Xt + sigma*sqrt(dt)*normalRandom();
    cout << "t=" << i*dt << " X_t=" << Xt << endl;
  }

The final value should be something like t=1 X_t=-0.245419 (here I have used the in built rand() function on a linux 64bit machine). Obviously given this path X has gone out of bounds before t=1, so let us know include a conditional break in the loop.

    // if X<a or X>b then output t and break
    if( Xt<a || Xt>b ){
      cout << " exit time t=" << i*dt << endl;
      break;
    }

My program now prints out a final line exit time t=0.35. We can verify that this is the first time the path exits the bounds. You now want to stress test your code so that it works under a variety of conditions. In particular the choice of dt and maximum number of steps need to either left to the user or set as something sensible given the other input variables (X0, a, b and sigma). To be brief we leave them as inputs and move the algorithm into a function. The function definition could look something like this

double timeToExit(double X0,double sigma,double a,double b,
		  double dt,double maxTime)

Inside the function you will need a local variable Xt to keep track of the path, and calculate the maximum number of steps given the values for dt and maxTime. You can then call the function in your main code.

  cout << " Time to exit " << timeToExit(0.56,0.4,0,1,0.1,20.) << "\n";

Now you can check this works by running the algorithm for multiple paths and check that each value returned looks sensible. Now we are going to calculate the expectation so we will need to store the number of paths and keep track of a running sum. So add them in your main file

  // number of paths
  int N=1000;
  // running sum
  double sum=0;

and now run a loop to add up the first exit time from each path. At the end you can output the average exit time over the number of runs

  cout << " Average time to exit " << sum/N << "\n";

Again we should test now with different values of N and dt to ensure that things are working properly. Once we are convinced we can put the Monte Carlo part in its own function. The definition might look like

double averageTimeToExit(int N,double X0,double sigma,double a,double b,double dt,double maxTime)

Finally your code might look something like:

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double uniformRandom()
{
  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}
 // return a normally distributed random number
double normalRandom()
{
  double u1=uniformRandom();
  double u2=uniformRandom();
  return cos(8.*atan(1.)*u2)*sqrt(-2.*log(u1)); 
}

double timeToExit(double X0,double sigma,double a,double b,
		  double dt,double maxTime)
{
  // store the current value of X
  double Xt=X0;
  // now plot out a path in time
  for(int i=1;i<=maxTime/dt;i++)
  {
    Xt = Xt + sigma*sqrt(dt)*normalRandom();
    // if X<a or X>b then output t and break
    if( Xt<a || Xt>b )return i*dt;
  }
  return maxTime;
}

double averageTimeToExit(int N,double X0,double sigma,double a,double b,
		  double dt,double maxTime)
{
  // running sum
  double sum=0;
  for(int i=0;i<N;i++)
  {
    sum = sum + timeToExit(X0,sigma,a,b,dt,maxTime);
  }
  return sum/N;
}

int main()
{
  cout << " Average time to exit " << averageTimeToExit(1000,0.56,0.4,0,1,0.1,20.) << "\n";
}