**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>usingnamespacestd; // 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();returncos(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 timefor(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 breakMy program now prints out a final lineif( Xt<a || Xt>b ){ cout << " exit time t=" << i*dt << endl;break; }

`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>usingnamespacestd; // 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();returncos(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 timefor(int i=1;i<=maxTime/dt;i++) { Xt = Xt + sigma*sqrt(dt)*normalRandom(); // if X<a or X>b then output t and breakif( Xt<a || Xt>b )returni*dt; }returnmaxTime; } 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); }returnsum/N; } int main() { cout << " Average time to exit " << averageTimeToExit(1000,0.56,0.4,0,1,0.1,20.) << "\n"; }

## No comments:

## Post a Comment