// OptionTree 1.0
// Calculate the price for plain vanilla European options
// on non-dividend paying equity using a binomial tree

#include <iostream>
#include <cmath>
#include <algorithm>

using namespace std;

const char PROGNAME[] = "OptionTree";
const char PROGVERSION[] = "1.0";
const char PROGAUTHOR[] = "J. Goodwin";
const char PROGDATE[] = "18 Nov 2004";
const char HRULE[] = "-------------------------------\n";

// --------------------------------
// MAIN program
// No arguments are used
// --------------------------------

int main(int argc, char *argv[])
{
    // --------------------------------
    // Variables
    // --------------------------------
    char   optionClass; // should be C (Call) or P (Put)
    double k;           // strike price
    double s0;          // initial stock price
    double sigma;       // stock volatility
    double r;           // risk-free interest rate (constant)
    double t;           // time to maturity (years)
    int    n;           // number of time steps
    double u;           // stock price up-step
    double d;           // stock price down-step
    double p;           // probability of up-step
    double dt;          // time step size
    double payoff;      // payoff calculated for each node

    int    i,j;         // counters
    
    // --------------------------------
    // Display Header
    // --------------------------------

    cout << PROGNAME << " v" << PROGVERSION << endl << PROGAUTHOR << " -- " << PROGDATE << endl;
    cout << HRULE;

    // --------------------------------    
    // Request Program Inputs
    // Limited validation
    // --------------------------------
    
    // Option Class
    for(;;)
    {
        cout << "(C)all or (P)ut option?";
        cin >> optionClass;
        optionClass = toupper(optionClass);
        if ('C' == optionClass || 'P' == optionClass) break;
    }

    // Strike Price
    cout << "What is the strike price?";
    cin >> k;
    
    // Initial stock price
    cout << "What is the initial stock price?";
    cin >> s0;
    
    // Stock volatility
    cout << "What is the stock volatility (sigma)?";
    cin >> sigma;
    
    // Risk-free interest rate
    cout << "What is the risk-free interest rate?";
    cin >> r;
    
    // Time to maturity
    cout << "What is the time to maturity (in years)?";
    cin >> t;
    
    // Number of time steps - model parameter
    cout << "How many time steps should be used?";
    cin >> n;

    // --------------------------------
    // Deduce binomial tree values
    // --------------------------------

    dt = t/n;
    u = exp(sigma*sqrt(dt));
    d = 1.0/u;
    p = (exp(r*dt)-d)/(u-d);

    // --------------------------------
    // Build Stock Price Tree
    // --------------------------------
    
    // Create a pointer array, each element will represent a time step
    double *s[n+1];   // stock price array

    // Populate first node in first time step
    // In this prototype we are not trapping memory allocation errors
    s[0] = new double;
    *s[0] = s0;

    // Create tree going forward, generating stock price values
    // Nodes are numbered [i][j]
    //   i represents the time step (0 to n)
    //   j is the node number within the time step, from highest stock value (0) to lowest (i)
    // This means that the [i][j] node will step forward to a stock
    // price increase at [i+1][j] and a decrease at [i+1][j+1]
    for(i=1; i<n+1; i++)
    {
        // i+1 nodes will be in the ith column (time step)
        s[i] = new double[i+1];
        // populate stock prices
        s[i][0] = s[i-1][0]*u;  // topmost node first
        for(j=1; j<i; j++)
        {
            s[i][j] = s[i-1][j-1]*d;
        }        
    }
    
    // --------------------------------
    // Build Option Value Tree
    // --------------------------------
    
    double *v[n+1];   // option value array

    // Populate the final column of option price array
    // Is equal to the payoff from the option which depends
    // on where the stock price is relative to the strike price
    v[n] = new double[n+1];
    for(i=0; i<n; i++)
    {
        if ('C' == optionClass)
        {
            v[n][i] = max(0.0, s[n][i]-k);
        }
        else
        {
            v[n][i] = max(0.0, k-s[n][i]);
        }
    }
    
    // Now traverse the tree backwards, populating option value
    for(i=n-1; i>=0; i--)
    {
        v[i] = new double[i+1]; // option values for ith time step
        for(j=0; j<=i; j++)
        {
            v[i][j] = (p*v[i+1][j]+(1-p)*v[i+1][j+1])*exp(-r*dt);
        }
    }

    // --------------------------------
    // Finish up
    // --------------------------------
    
    // Output the predicted option price
    cout.precision(10);
    cout << "Time step was:" << dt << endl;
    cout << "Probability of an upward stock movement was: " << p << endl;
    cout << "Upward stock factor was: " << u << endl;
    cout << "Downward stock factor was: " << d << endl;
    cout << HRULE;
    cout << "Calculated option price is:" << v[0][0] << endl;
    cout << HRULE;

    // Release the memory we had allocated
    for(i=0; i<n+1; i++)
    {
        delete s[i];
        delete v[i];
    }
    
    system("PAUSE");	
    return 0;
}