#include <vector>
#include <gsl/gsl_cdf.h>
#include "var_alea.hpp"
#include <algorithm>
#include <iostream>
#include "low_discrepancy.hpp"


template <typename L>
std::vector<double> monte_carlo(int n, L X)
{
	std::vector<double> result(3,0);
	double x;
	for (int j = 0; j < n; j++) {
		x = X();
		result[0] += x;
		result[1] += x*x;
	}
	result[0] /= (double) n;
	result[1] = (result[1] - n*result[0]*result[0])/(double)(n-1);
	result[2] = 1.96*sqrt(result[1]/(double) n);
	return result;
}

double pos (double x){
    return x>0?x:0;
}

double frac_part(double x){
    return x - floor(x);
}

struct asian_option : public var_alea<double>
{
	asian_option(double r, double T, double S0, double V, int d, double K)
        : r(r), T(T), S0(S0), V(V), d(d), K(K), G(0,1) {};
        
	double operator()() {
        std::vector<double> S(d);
        S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*G());
        for(int i=1;i<d;i++){
            S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*G());
        }
        double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
        return exp(-r*T)*pos(temp-K);
    };
    
	private:
        double r;
        double T;
        double S0;
        double V;
        int d;
        double K;
        gaussian G;
};

struct asian_option_qmc : public var_alea<double>
{
	asian_option_qmc(double r, double T, double S0, double V, int d, double K)
        : r(r), T(T), S0(S0), V(V), d(d), K(K), G(d), U(0,1) {};
        
	double operator()() {
        std::vector<double> S(d);
        std::vector<double> sob(d); 
        sob = G();
        S[0]= S0*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(U()+sob[0]), 1));
        for(int i=1;i<d;i++){
            S[i]=S[i-1]*exp((r-V*V/2)*(T/d)+V*sqrt(T/d)*gsl_cdf_gaussian_Pinv(frac_part(U()+sob[i]), 1));
        }
        double temp = std::accumulate(S.begin(), S.end(), 0.)/d;
        return exp(-r*T)*pos(temp-K);
    };
    
	private:
        double r;
        double T;
        double S0;
        double V;
        int d;
        double K;
        sobol G;
        uniform U;
};


int main(){
    init_alea(1);
    asian_option A(0.05, 1.0, 50.0, 0.1, 16, 45);
    asian_option_qmc B(0.05, 1.0, 50.0, 0.1, 16, 45);
    int M= 1000000;
    int N=10000;
    int I=100;
    
    std::vector<double> meanvar = monte_carlo(M, A);
    std::cout<<"espérance "<<meanvar[0] <<" IC "<<meanvar[2]<<std::endl;
    std::vector<double> temp;
    double m = 0;
    double s = 0;
    for(int i=0;i<I;i++){
        temp = monte_carlo(N,B);
        m+=temp[0];
        s+=temp[0]*temp[0];
        }
    m = m/I;
    s = s/I - m*m;
    std::cout<<"espérance "<<m <<" IC "<<sqrt(s)*1.96/10<<std::endl;

    
    return 0;
}