// Black-Scholes by simulation.
// TestCandidat.com


#include <iostream>
#include <iomanip>
#include <vector>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <algorithm>
#include "NormalFunctions.h"
using namespace std;

//BSSimul compute the price of an european option using simulation
    // S Spot Price
    // K Strike Price
    // T Maturity in Years
    // rf Interest Rate
    // q Dividend  yeild
    // v Volatility
    // PutCall 'P' for put and 'C' for call
    // Nsims Number of simulations

double BSSimul(double S,double K,double v,double T,double rf,double q,char PutCall) {

double u1,u2,Z;
double pi = 3.141592653589793;
int Nsims = 1e5; // Number of simulations
vector<double> ST(Nsims, 0.0);// Initialize terminal prices S(T)
vector<double> ST_K(Nsims, 0.0);// Initialize call payoff
vector<double> K_ST(Nsims, 0.0);// Initialize put payoff
for (int i=0; i<=Nsims-1; i++) {
// Independent uniform random variables
u1 = ((double)rand() / ((double)(RAND_MAX)+(double)(1)) );
u2 = ((double)rand() / ((double)(RAND_MAX)+(double)(1)) );
// Floor u1 to avoid errors with log function
u1 = max(u1,1.0e-10);
// Z ~ N(0,1) by Box-Muller transofmration
Z = sqrt(-2.0*log(u1)) * sin(2*pi*u2);
ST[i] = S*exp((rf - q - 0.5*v*v)*T + v*sqrt(T)*Z);// Simulated terminal price S(T)
ST_K[i] = max(ST[i] - K, 0.0);// Call payoff
K_ST[i] = max(K - ST[i], 0.0);// Put payoff
}
// Simulated prices as discounted average of terminal prices
double BSCallSim = exp(-rf*T)*VecMean(ST_K);
double BSPutSim  = exp(-rf*T)*VecMean(K_ST);
if(PutCall =='C')
    return BSCallSim;
else
    return BSPutSim;
}