#include <iostream>
#include <cmath>
#include <math.h>
#include <limits>
#include <algorithm>
#include <Eigen/Dense>
#include <Eigen/Geometry>
#include <random>
#include <vector>
#include <time.h>

using namespace Eigen;
using namespace std;


double FM(double T, double r, double q, double sigma, double S0, double K, int M, int N);
MatrixXd generateGaussianNoise(int M, int N);        // Generates Normally distributed random numbers
double Black_Scholes(double T, double K, double S0, double r, double sigma);
double phi(long double x);
VectorXd time_vector(double min, double max, int N);
MatrixXd call_payoff(MatrixXd S, double K);

int main(){
    double r = 0.03;        // Riskless interest rate
    double q = 0.0;         // Divident yield
    double sigma = 0.15;    // Volatility of stock
    int T = 1;              // Time (expiry)
    int N = 100;            // Number of time steps
    double K = 100;         // Strike price
    double S0 = 102;        // Initial stock price
    int M = 100;            // Number of paths              // Current issue

    FM(T,r,q,sigma,S0,K,M,N);

    return 0;
}

MatrixXd generateGaussianNoise(int M, int N){
    MatrixXd Z(N,M);
    random_device rd;
    mt19937 e2(time(0));
    normal_distribution<double> dist(0.0, 1.0);
    for(int i = 0; i < M; i++){
        for(int j = 0; j < N; j++){
            Z(i,j) = dist(e2);
        }
    }
    return Z;
}

double phi(double x){
    return 0.5 * erfc(-x * M_SQRT1_2);
}

double Black_Scholes(double T, double K, double S0, double r, double sigma){
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double d_2 = (log(S0/K) + (r-sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C = S0*phi(d_1) - phi(d_2)*K*exp(-r*T);
    return C;
}

VectorXd time_vector(double min, double max, int N){
    VectorXd m(N + 1);
     double delta = (max-min)/N;
     for(int i = 0; i <= N; i++){
             m(i) = min + i*delta;
     }
    return m;
}

MatrixXd call_payoff(MatrixXd S, double K){
    MatrixXd result(S.rows(),S.cols());
    for(int i = 0; i < S.rows(); i++){
        for(int j = 0; j < S.cols(); j++){
            if(S(i,j) - K > 0){
                result(i,j) = S(i,j) - K;
            }else{
                result(i,j) = 0.0;
            }
        }
    }
    return result;
}

double FM(double T, double r, double q, double sigma, double S0, double K, int M, int N){
    MatrixXd Z = generateGaussianNoise(M,N);
    double dt = T/N;
    VectorXd t = time_vector(0.0,T,N);

    // Generate M paths of stock prices
    MatrixXd S(M,N+1);
    for(int i = 0; i < M; i++){
        S(i,0) = S0;
        for(int j = 1; j <= N; j++){
            S(i,j) = S(i,j-1)*exp((double) (r - q - pow(sigma,2.0))*dt + sigma*sqrt(dt)*(double)Z(i,j-1));
        }
    }

    //
    // If path i is "alive" at time index j - 1 < N, generate the price for time index j, denoted as S = S_ij
    // Case for call option:
    // If j = N, the option is expired with value V = exp(-rT)(S-K)^+ and path i is finished
    // If j < N, calculate S_c = f_C(S)
    // If S > S_c, the option is exercised with value V_i = exp(-rT)(S-K)^+ and path i is stopped. Otherwise,
    // the option is held and path continues to live to the next step j+1
    //
    // Case for put option:
    // If j = N, the option is expired with value V = exp(-rT)(K-S)^+ and path i is finished
    // If j < N, calculate S_p = f_p(S)
    // if S < S_p, the option is exercised with value V_i and path i is stopped. Otherwise,
    // the option is held and path continues to live to the next step j+1.

    // Compute S_c parameters and S_p
    double m = 2*r/(pow(sigma,2.0));
    double n = 2*(r-q)/(pow(sigma,2.0));
    VectorXd k(t.size());
    for(int i = 0; i < k.size(); i++){
        k(i) = 1.0 - exp((double) -r*(double)(T - t(i)));                       // Note the t vector (not sure if this is correct)
    }
    VectorXd Q_2(k.size());
    VectorXd Q_1(k.size());
    for(int i = 0; i < Q_2.size(); i++){
        Q_1(i) = (-1*(n-1) + sqrt((double)(n-1)*(n-1) + (double)4*m/(double)(k(i))))/2.0;                       // Q_1 < 0
        Q_2(i) = (-1*(n-1) + sqrt((double)(n-1)*(n-1) + (double)4*m/(double)(k(i))))/2.0;                       // Q_2 > 0
    }
    double d_1 = (log(S0/K) + (r+sigma*sigma/2.)*(T))/(sigma*sqrt(T));
    double C_e = Black_Scholes(T, K, S0, r, sigma);                 // C_e(S) is the European call option price calculated by Black-Scholes
    double Delta = exp(-q*T)*phi(d_1);
    MatrixXd V(M,N+1);
    VectorXd S_c(Q_2.size());
    MatrixXd call_fun = call_payoff(S,K);
    for(int j = 0; j < N + 1; j++){
        for(int i = 0; i < M; i++){
            if(j == N){
                V(i,j) = exp(-r*T)*call_fun(i,j); //////////////
                //cout << "The option is expired with value " << V(i) << " and path " << i << " is finished" << endl;
            }
            else if(j < N){
                S_c(j) = Q_2(j)*(C_e + K)/(Q_2(j) - (1 - Delta));
            }
            else if (S(i,j) > S_c(j)){
                V(i,j) = exp(-r*T)*call_fun(i,j); ///////////////
                //cout << "The option is expired with value " << V(i) << " and path " << i << " is finished" << endl;
            }
        }
    }
    double Value = 0.0;
    for(int i = 0; i < V.rows(); i++){
        for(int j = 0; j < V.cols(); j++){
            Value += V(i,j);
        }
    }
    Value = 1.0/M * Value;
    cout << C_e << endl;
    cout << endl;
    cout << Value << endl;

}