/*
      A simple code for option valuation using the explicit forward Euler method
   for the class Derivative Securities, fall 2010
	  http://www.math.nyu.edu/faculty/goodman/teaching/DerivSec10/index.html
	  Written for this purpose by Jonathan Goodman, instructor.
   
      Assignment 8
	  
*/


#include <iostream>
#include <fstream>
#include <math.h>
#define NSPOTS  100   /* The number of spot prices computed  */

/*   A program to compute a simple binomial tree price for a European style put option  */

using namespace std;


//    The pricer, main is at the bottom of the file



void FE(                   // Solve a pricing PDE using the forward Euler method
   
   double T, double sigma, double r, double K,     // The standard option parameters
   double Smin, double Smax,                       // The min and max prices to return
   int nPrices,                                    // The number of prices to compute between Smin and Smax, 
                                                   // Determines the accuracy and the cost of the computation
   double prices[],                                // An array of option prices to be returned.  
   double intrinsic[],                             // The intrinsic value at the same prices
   double spots[],                                 // The corresponding spot prices, computed here for convenience.
                                                   // Both arrays must be allocated in the calling procedure
   double *SEarly ) {                              // The early exercise boundary
                                                   
//      Setup for the computation, compute computational parameters and allocate the memory

   double xMin = log(Smin);                                 //  Work in the log variable
   double xMax = log(Smax);
   double dx   = ( xMax - xMin )/ ( (double( nPrices - 1 ) ) );  // The number of gaps is one less than the number of prices
   double CFL  = .8;                                        // The time step ratio
   double dt   = CFL*dx*dx/sigma;                           // The forward Euler time step size, to be adjusted slightly
   int    nTimeSteps = (int) (T/dt);                        // The number of time steps, rounded down to the nearest integer
          nTimeSteps++;                                     // Now rounded up
          dt   = T / ( (double) nTimeSteps );               // Adjust the time step to land precisely at T in n steps.
   int    nx   = nPrices + 2*nTimeSteps;                    // The number of prices at the final time, growing by 2 ...
                                                            // ... each time step
          xMin = xMin - nTimeSteps*dx;                      // The x values now start here
   double *fOld;                                            // The values of the pricing function at the old time
           fOld   = new double [nx];                        // Allocated using old style C++ for simplicity
   double *fNew;                                            // The values of the pricing function at the new time
           fNew   = new double [nx];
   double *V;                                               // The intrinsic value = the final condition
           V      = new double [nx];
           
//       Get the final conditions and the early exercise values 

   double x;      // The log variable
   double S;      // A stock price = exp(x)
   int    j;
   for ( j = 0; j < nx; j++ ) {
      x = xMin + j*dx;
      S = exp(x);
      if ( S < K ) V[j] = K-S;      // A put struck at K
      else         V[j] = 0;
      fOld[j] = V[j];               // The final condition is the intrinsic value
    }
    
//      The time stepping loop

   double alpha, beta, gamma;                   // The coefficients in the finite difference formula
   alpha = beta = gamma = .333333333333;        //   XXXXXXXXXXXXXXXXXXXXXXXXXXX
   int jMin = 1;                                // The smallest and largest j ... 
   int jMax = nx - 1;                           // ... for which f is updated.  Skip 1 on each end the first time.
   int jEarly    ;                              // The last index of early exercise
   for ( int k = nTimeSteps; k > 0; k-- ) {     // This is, after all, a backward equation
      jEarly = 0;                               // re-initialize the early exercise pointer
      for ( j = jMin; j < jMax; j++ ) {                                // Compute the new values
         x = xMin + j*dx;                                              // In case the coefficients depend on S
         S = exp(x);
         fNew[j] = alpha*fOld[j-1] + beta*fOld[j] + gamma*fOld[j+1];   // Compute the continuation value
         if ( fNew[j] < V[j] ) {
            fNew[j] = V[j];                                            // Take the max with the intrinsic value
            jEarly  = j;                                               // and record the largest early exercise index
          }
       }
      for ( j = jMin; j < jMax; j++ )                                  // Copy the new values back into the old array
         fOld[j] = fNew[j];
      jMin++;                                                          // Move the boundaries in by one
      jMax--;
    }
    
//     Copy the computed solution into the desired place
   
   jMin--;      // The last decrement and increment were mistakes
   jMax++;
   int i = 0;                             // The index into the output array
   for ( j = jMin; j < jMax; j++ ) {      // Now the range of j should match the output array
      x = xMin + j*dx;                                              
      S = exp(x);
      prices[i]    = fOld[j];
      intrinsic[i] = V[j];
      spots[i++]   = S;                     // Increment i after all copy operations
    }
   double xEarly = xMin + jEarly*dx;
         *SEarly = exp(xEarly);             // Pass back the computed early exercise boundary
   
   delete fNew;       // Be a good citizen and free the memory when you're done.
   delete fOld;
   delete V;
   return;
 }
       
int main() {

   cout << "Hello " << endl;
   
   ofstream csvFile;              // The file for output, will be csv format for Excel.
   csvFile.open ("PutPrice.csv");
   
   double sigma = .3;
   double r     = .003;
   double T     = .5;
   double K     = 100;
   double Smin  = 60;
   double Smax  = 180;
   double prices[NSPOTS];
   double intrinsic[NSPOTS];
   double spots[ NSPOTS];
   double SEarly;
   
   FE( T, sigma, r,  K, Smin, Smax, NSPOTS, prices, intrinsic, spots, &SEarly ); 
   
   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the spot prices for plotting
      csvFile << spots[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", "; //   Don't put a comma after the last value
    }
   csvFile << endl;
   
   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the intrinsic prices for plotting
      csvFile << intrinsic[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", "; //   Don't put a comma after the last value
    }
   csvFile << endl;

   
   for ( int j = 0; j < NSPOTS; j++ ) {        //   Write out the computed option prices
      csvFile << prices[j];
      if ( j < (NSPOTS - 1) ) csvFile << ", ";
    }
   csvFile << endl;
   
   csvFile << "Critical price," << SEarly << endl;
   
   csvFile << "T ," << T << endl;
   csvFile << "r ," << r << endl;
   csvFile << "sigma ," << sigma << endl;
   csvFile << "strike ," << K << endl;
   
   return 0 ;
   
 }