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

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib 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.
*/

/*
 NesQuant is an extension to QuantLib
 http://www.nielses.dk/quantlib/nesquant
*/

#include <nq/PricingEngines/lsmsvjdengine.hpp>
#include <iostream>
using namespace std;

namespace NesQuant {

    void LSMSVJDEngine::calculate() const {

        double s0 = arguments_.underlying,
               v0 = arguments_.volatility * arguments_.volatility;
        double maturity = arguments_.maturity;
        PlainVanillaPayoff payoff(arguments_.type, arguments_.strike);

        GaussianRandomSequenceGenerator gen =
            PseudoRandom::make_sequence_generator(timeSteps_, seed_);

        Handle<GaussianSVJDPathGenerator> pathGenerator =
            Handle<GaussianSVJDPathGenerator>
                       (new GaussianSVJDPathGenerator(
                           arguments_.riskFreeRate,
                           arguments_.dividendYield,
                           arguments_.volatility,
                           arguments_.volatilityOfVolatility,
                           arguments_.steadyStateVolatility,
                           arguments_.meanReversionRate,
                           arguments_.correlationUnderlyingVolatility,
                           arguments_.jumpIntensity,
                           arguments_.jumpMean,
                           arguments_.jumpStandardDeviation,
                           maturity, timeSteps_,
                           gen));

        // number of time steps between values used as exercise times
        Size timeDistance = timeSteps_ / exerciseTimes_;

        // if antithetic paths are used, twice as many paths are used
        Size modSamples; // modified samples
        if (antitheticVariate_)
            modSamples = 2 * samples_;
        else
            modSamples = samples_;

        // Initialize control variate
        double controlVariateValue;
        std::vector<double> controlVariatePathValue(modSamples);
        EuropeanPathPricer_old controlPP(arguments_.type, s0, arguments_.strike,
                              QL_EXP(-arguments_.riskFreeRate * maturity), false);
        if (controlVariate_) {
            SVJDEngine controlPE;

            SVJDArguments* controlArguments =
                dynamic_cast<SVJDArguments*>(
                    controlPE.arguments());
            *controlArguments = arguments_;
            controlPE.calculate();

            const OptionValue* controlResults =
                dynamic_cast<const OptionValue*>(controlPE.results());
            controlVariateValue = controlResults->value;
        }

        // generate paths
        std::vector<std::vector<double> >
            spaths(modSamples, std::vector<double>(exerciseTimes_));
        std::vector<std::vector<double> >
            vpaths(modSamples, std::vector<double>(exerciseTimes_));

        GaussianSVJDPathGenerator::sample_type
            pathSample(pathGenerator->timeGrid(), 1.0);
        GaussianSVJDPathGenerator::sample_type
            volPathSample(pathGenerator->timeGrid(), 1.0);

        for (Size i = 0; i < modSamples; i++) {
            if (antitheticVariate_ && i%2 == 1)
                pathSample = pathGenerator->antithetic();
            else
                pathSample = pathGenerator->next();
            volPathSample = pathGenerator->volatility();

            double s = s0, v = v0, k = 0;
            for (Size j = 0; j < timeSteps_; j++) {
                s *= QL_EXP(pathSample.value[j]);
                v += volPathSample.value[j];
                if ((j + 1) % timeDistance == 0) {
                    spaths[i][k] = s;
                    vpaths[i][k] = v;
                    k++;
                }
            }

            if (controlVariate_) {
                controlVariatePathValue[i] = controlPP(pathSample.value);
            }
        }

        // set values at maturity
        std::vector<double> U(modSamples);
        exercisePoints_ = Matrix(modSamples, 3, -1.0); // (t, St, Vt)

        Size k = exerciseTimes_;
        DiscountFactor discount = QL_EXP(-arguments_.riskFreeRate * maturity);

        for (i = 0; i < modSamples; i++) {
            U[i] = discount * payoff(spaths[i][k-1]);
        }

        // move backwards
        double dt = maturity / exerciseTimes_;

        for (k = exerciseTimes_ - 1; k > 0; k--) {
            cout << "k : " << k << endl;

            // find in-the-money paths and (x, y) regression data
            std::vector<Size> itmpaths;
            std::vector<double> xdata, vdata, fdata;

            for (i = 0; i < modSamples; i++) {
                double s = spaths[i][k-1];
                if (payoff(s) > 0) {
                    itmpaths.push_back(i);
                    xdata.push_back(s);
                    vdata.push_back(vpaths[i][k-1]);
                    fdata.push_back(U[i]);
                }
            }

            if (itmpaths.size() > 0) {
                // do regression
                PolynomialFit pf(xdata, vdata, fdata, degree_);
                const Array contval = pf.value();

                // compare continuation value with exercise value
                discount = QL_EXP(-arguments_.riskFreeRate * k * dt);
                for (Size l = 0; l < itmpaths.size(); l++) {
                    i = itmpaths[l];
                    double exerval = discount * payoff(xdata[l]);

                    if (exerval > contval[l]) {
                        U[i] = exerval;
                        exercisePoints_[i][0] = k * dt;
                        exercisePoints_[i][1] = xdata[l];
                        exercisePoints_[i][2] = QL_SQRT(vdata[l]);
                    }
                }
            }
        }

        // adjust for control variate
        if (controlVariate_)
            for (i = 0; i < modSamples; i++) {
                U[i] = U[i] + 1.0 * (controlVariateValue - controlVariatePathValue[i]);
            }

        // calculate average
        if (antitheticVariate_)
            for (i = 0; i < modSamples; i+=2) {
                sampleAccumulator_.add((U[i] + U[i+1]) / 2);
            }
        else
            for (i = 0; i < samples_; i++) {
                sampleAccumulator_.add(U[i]);
            }

        double value = sampleAccumulator_.mean();

        // compare with exercise at time 0
        if (payoff(s0) > value) {
            value = payoff(s0);
            sampleAccumulator_.reset();
            for (i = 0; i < modSamples; i++) {
                if ((!antitheticVariate_) || i%2==0)
                    sampleAccumulator_.add(value);
                exercisePoints_[i][0] = 0;
                exercisePoints_[i][1] = s0;
                exercisePoints_[i][2] = arguments_.volatility;
            }
        }

        // return value
        results_.value = value;
        results_.errorEstimate = sampleAccumulator_.errorEstimate();
    }

}