/*
 Copyright (C) 2003 Niels Elken Sønderby
 Copyright (C) 2002, 2003 Ferdinando Ametrano
 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl

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

/*! \file svjdpathgenerator.hpp
    \brief stochastic volatility jump diffusion (SVJD) path generator
*/

#ifndef nesquant_montecarlo_svjdpath_generator_h
#define nesquant_montecarlo_svjdpath_generator_h

#include <ql/diffusionprocess.hpp>
#include <ql/MonteCarlo/path.hpp>
#include <ql/RandomNumbers/randomarraygenerator.hpp>

using namespace QuantLib::MonteCarlo;

namespace QuantLib {

    namespace MonteCarlo {

        template <class SG>
        class SVJDPathGenerator : public PathGenerator<SG> {
          public:
            SVJDPathGenerator(double riskFreeRate,
                  double dividendYield,
                  double volatility, double volatilityOfVolatility,
                  double steadyStateVolatility,
                  double meanReversionRate,
                  double correlationUnderlyingVolatility,
                  double jumpIntensity, double jumpMean,
                  double jumpStandardDeviation,
                  Time length, Size timeSteps,
                  const SG& generator);
            const sample_type& next() const;
            const sample_type& volatility() const;
            const sample_type& antithetic() const;
          private:
            double r_, q_;
            double v0_;
            double sigmav_;
            double alphabeta_; // alpha / betastar
            double betastar_;
            double rho_;
            double lambda_;
            double kmeanstar_;
            double delta_;

            mutable sample_type volatility_;

            RandomNumbers::MersenneTwisterUniformRng urng_;
            RandomNumbers::ICGaussianRng<
                               RandomNumbers::MersenneTwisterUniformRng,
                               Math::InverseCumulativeNormal> grng_;

            mutable std::vector<double> random1_;
            mutable std::vector<double> random2_;
            mutable std::vector<bool> isJump_;
            mutable std::vector<double> randomj_;
       };

        template <class SG>
        SVJDPathGenerator<SG>::SVJDPathGenerator(
                double riskFreeRate, double dividendYield,
                double volatility, double volatilityOfVolatility,
                double steadyStateVolatility, double meanReversionRate,
                double correlationUnderlyingVolatility,
                double jumpIntensity, double jumpMean,
                double jumpStandardDeviation,
                Time length, Size timeSteps,
                const SG& generator)
        : PathGenerator<SG>
            (Handle<DiffusionProcess>
                 (new SquareRootProcess(1, 1, 1)), // dummy
             length, timeSteps, generator),
            volatility_(Path(timeGrid_),1.0),
            r_(riskFreeRate), q_(dividendYield),
            v0_(volatility * volatility),
            sigmav_(volatilityOfVolatility),
            alphabeta_(steadyStateVolatility * steadyStateVolatility),
            betastar_(meanReversionRate),
            rho_(correlationUnderlyingVolatility),
            lambda_(jumpIntensity),
            kmeanstar_(jumpMean), delta_(jumpStandardDeviation),
            random1_(next_.value.size()), random2_(next_.value.size()),
            isJump_(next_.value.size()), randomj_(next_.value.size())
            { }

        template <class SG>
        inline const typename SVJDPathGenerator<SG>::sample_type&
        SVJDPathGenerator<SG>::next() const {
            typedef typename SG::sample_type sequence_type;

            const sequence_type& sequence1 = generator_.nextSequence();
            std::copy(sequence1.value.begin(), sequence1.value.end(),
                      random1_.begin());
            const sequence_type& sequence2 = generator_.nextSequence();
            std::copy(sequence2.value.begin(), sequence2.value.end(),
                      random2_.begin());

            double V = v0_;
            double dW, dWv, dJ = 0;
            double Zdrift, Zdiff;
            double Vdrift, Vdiff;

            double dt;
            Time t;

            for (Size i = 0; i < next_.value.size(); i++) {
                t = timeGrid_[i+1];
                dt = timeGrid_.dt(i);

                dW = random1_[i];
                dWv = rho_ * dW + QL_SQRT(1 - rho_ * rho_) * random2_[i];

                Vdrift = (betastar_ * (alphabeta_ - V)) * dt;
                Vdiff = sigmav_ * QL_SQRT(V * dt) * dWv;
                V += Vdrift + Vdiff;
                if (V < 0.00001) {
                    Vdiff += -V + 0.00001;
                    V = 0.00001;
                }

                if (lambda_ > 0) {
                    isJump_[i] = (urng_.next().value < lambda_ * dt);
                    if (isJump_[i]) {
                        randomj_[i] = grng_.next().value;
                        dJ = QL_LOG(1 + kmeanstar_)
                             - delta_ * delta_ / 2
                             + delta_ * randomj_[i];
                    } else {
                        dJ = 0;
                    }
                }

                Zdrift = (r_ - q_ - lambda_ * kmeanstar_ - 0.5 * V) * dt;
                Zdiff = QL_SQRT(V * dt) * dW + dJ;

                next_.value.drift()[i] = Zdrift;
                next_.value.diffusion()[i] = Zdiff;
                volatility_.value.drift()[i] = Vdrift;
                volatility_.value.diffusion()[i] = Vdiff;
            }

            return next_;
        }

        template <class SG>
        inline const typename SVJDPathGenerator<SG>::sample_type&
        SVJDPathGenerator<SG>::antithetic() const {
            double V = v0_;
            double dW, dWv, dJ = 0;
            double Zdrift, Zdiff;
            double Vdrift, Vdiff;

            double dt;
            Time t;

            for (Size i = 0; i < next_.value.size(); i++) {
                t = timeGrid_[i+1];
                dt = timeGrid_.dt(i);

                dW = - random1_[i];
                dWv = - (rho_ * dW + QL_SQRT(1 - rho_ * rho_) * random2_[i]);

                Vdrift = (betastar_ * (alphabeta_ - V)) * dt;
                Vdiff = sigmav_ * QL_SQRT(V * dt) * dWv;
                V += Vdrift + Vdiff;
                if (V < 0.00001) {
                    Vdiff += -V + 0.00001;
                    V = 0.00001;
                }

                if (lambda_ > 0) {
                    if (isJump_[i]) {
                        dJ = QL_LOG(1 + kmeanstar_)
                             - delta_ * delta_ / 2
                             + delta_ * (- randomj_[i]);
                    } else {
                        dJ = 0;
                    }
                }

                Zdrift = (r_ - q_ - lambda_ * kmeanstar_ - 0.5 * V) * dt;
                Zdiff = QL_SQRT(V * dt) * dW + dJ;

                next_.value.drift()[i] = Zdrift;
                next_.value.diffusion()[i] = Zdiff;
                volatility_.value.drift()[i] = Vdrift;
                volatility_.value.diffusion()[i] = Vdiff;
            }

            return next_;
        }

        template <class SG>
        inline const typename SVJDPathGenerator<SG>::sample_type&
        SVJDPathGenerator<SG>::volatility() const {
            return volatility_;
        }

        typedef SVJDPathGenerator
                <RandomNumbers::GaussianRandomSequenceGenerator>
                     GaussianSVJDPathGenerator;

    }

}

#endif