/*
 Copyright (C) 2003 Niels Elken Sønderby

 This file is part of QuantLib for Mathematica, a Mathematica extension for
 QuantLib, a free-software/open-source financial C++ library
 http://www.nielses.dk/quantlib/mma
 http://quantlib.org/

 QuantLib for Mathematica is free software: you can redistribute it and/or
 modify it under the terms of the QuantLib license.  You should have received
 a copy of the license along with this program; if not, please email
 ferdinando@ametrano.net The license is also available online at
 http://quantlib.org/html/license.html

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include "mmamontecarlo.hpp"

using namespace QuantLib::RandomNumbers;
using namespace QuantLib::VolTermStructures;
using namespace QuantLib::TermStructures;
using namespace QuantLib::Calendars;
using namespace QuantLib::DayCounters;

void qlMonteCarloPathBlackScholes(double s0, double riskFreeRate,
                                  double dividendYield, double volatility,
                                  double length, int timeSteps, int anti)
{
    try {
        using namespace QuantLib::MonteCarlo;

        // generate the risk neutral path
        static long seed = 0;

        if (seed == 0)
            seed = QL_TIME(0);
        else
            seed++;

        RandomNumbers::GaussianRandomSequenceGenerator gen =
            MonteCarlo::PseudoRandom::make_sequence_generator(timeSteps, seed);

        RelinkableHandle<BlackVolTermStructure> volTS = RelinkableHandle<BlackVolTermStructure>(makeFlatVolatility(volatility));

        RelinkableHandle<TermStructure> riskFreeTS = RelinkableHandle<TermStructure>(Handle<TermStructure>(new ConstantTS(riskFreeRate)));
        RelinkableHandle<TermStructure> dividendTS = RelinkableHandle<TermStructure>(Handle<TermStructure>(new ConstantTS(dividendYield)));
//    Handle<TermStructure> rfCurve = makeFlatCurve(rRate);

        Handle<DiffusionProcess> bs(new
            BlackScholesProcess(riskFreeTS, dividendTS, volTS, s0));

        // Exercise dates
        TimeGrid grid(length, timeSteps);

        Handle<MonteCarlo::GaussianPathGenerator> pathGenerator =
            Handle<MonteCarlo::GaussianPathGenerator> (
                new MonteCarlo::GaussianPathGenerator(bs, grid, gen));

        MonteCarlo::GaussianPathGenerator::sample_type 
            pathSample = pathGenerator->next();

        MonteCarlo::Path path = pathSample.value;

        if (anti) {
            MLPutFunction(stdlink, "List", 2);
        }

        double *ret = new double[timeSteps];

        for (int i=0; i<=anti; i++) { // if antithetic run twice
            if (i==1) { 
                pathSample = pathGenerator->antithetic();
                path = pathSample.value;
            }
            
            double s = 0;
            ret[0] = s0;

            for (Size i = 1; i < timeSteps; i++) {
                s += path[i];
                ret[i] = s0 * QL_EXP(s);
            }
            MLPutRealList(stdlink, ret, timeSteps);
        }

        delete[] ret;
    } QLML_HANDLE_EXCEPTIONS
}

void nqMonteCarloPathSVJD(double s0, double riskFreeRate,
                  double dividendYield, double volatility,
                  double volatilityOfVolatility,
                  double steadyStateVolatility,
                  double meanReversionRate,
                  double correlationUnderlyingVolatility,
                  double jumpIntensity, double jumpMean,
                  double jumpStandardDeviation,
                  double length, int timeSteps, int anti)
{
    try {
        using namespace QuantLib::MonteCarlo;

        // generate the risk neutral path
        static long seed = 0;

        if (seed == 0)
            seed = QL_TIME(0);
        else
            seed++;

        RandomNumbers::GaussianRandomSequenceGenerator gen =
            MonteCarlo::PseudoRandom::make_sequence_generator
                (timeSteps, seed);

        Handle<GaussianSVJDPathGenerator> pathGenerator =
            Handle<GaussianSVJDPathGenerator>
                  (new GaussianSVJDPathGenerator(
                  riskFreeRate, dividendYield,
                  volatility, volatilityOfVolatility,
                  steadyStateVolatility, meanReversionRate,
                  correlationUnderlyingVolatility,
                  jumpIntensity, jumpMean,
                  jumpStandardDeviation,
                  length, timeSteps, gen));

        MonteCarlo::GaussianPathGenerator::sample_type 
            pathSample = pathGenerator->next();

        MonteCarlo::Path path = pathSample.value;

        if (anti) {
            MLPutFunction(stdlink, "List", 2);
        }

        double *ret = new double[timeSteps];

        for (int i=0; i<=anti; i++) { // if antithetic run twice
            if (i==1) { 
                pathSample = pathGenerator->antithetic();
                path = pathSample.value;
            }
            
            double s = 0;
            ret[0] = s0;

            for (Size i = 1; i < timeSteps; i++) {
                s += path[i];
                ret[i] = s0 * QL_EXP(s);
            }
            MLPutRealList(stdlink, ret, timeSteps);
        }

        delete[] ret;
    } QLML_HANDLE_EXCEPTIONS
}