/*     A procedure that generates a random sample path for an SDE.
       Written by Jonathan Goodman for assignment 8 in the class 
	   http://www.math.nyu.edu/faculty/goodman/teaching/DerivSec09/index.html
	   
	   It uses the native C random number generator rand(), which may not be that good.
	   Feel free to replace this with a better one, but remember to set the seed.                */

#include <iostream>
#include <math.h>

using namespace std;

int randn( double Z[], int nZ);

//    Put the definition of the model here, as macros to be used in the SDE below.
//    Specify the model parameters and give functional forms for the drift and noise.
//    a(x) is the drift, and b(x) is the noise.

#define mu    .1
#define sigma .3
#define a(x) (mu*(x))      //      See http://en.wikipedia.org/wiki/C_preprocessor to ...
#define b(x) (sigma*(x))   //  ... learn about parentheses here.



int path(             // Return 0 if it worked, return 1 if not.
   double X[],        // The computed sample path.  The array of length n+1 must be allocated already.
   int    n,          // The number of steps in the path.  X[0] will be xStart.  X[n] is the last value (n+1 in all)
   double tStart,     // The starting time of the path.
   double T,          // The ending time of the path.
   double xStart) {   // The starting x value of the path.
   
   if ( n <= 0 ) {    // Error checking. Should do more, but am lazy.
      cout << "Path quitting because n = " << n << endl;
	  return 1; }

   
   double* Z;           // Allocate storage for the standard normals that drive the path.
   Z = new double[n];   // Should error check this, but lazy.
   randn( Z, n);        // Get the normals.
   
   double dt = (T - tStart)/n;  // Calculate the time step needed to fit n time steps in the interval [tStart,T]
   double sdt = sqrt(dt);       // Avoid evaluating this over and over.
   double t;                    // The present value of t.
   double XX;                   // The value of X at the present time.
   
   t    = tStart;               // Initializations
   XX   = xStart;
   X[0] = XX;
   
   for ( int i = 1; i <= n; i++ ) {   // Note that i starts at i=1 and ends at n rather than i=0 and n-1.
      XX  += a(XX)*dt + b(XX)*sdt*Z[i-1];
	  X[i] = XX;
	}


   delete[] Z;    // Plug those memory leaks.
   return 0;
   
}