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

/*! \file svjdengine.hpp
    \brief SVJD pricing engine for European options

    \fullpath
    nq/PricingEngines/%svjdengine.hpp
*/


#ifndef nesquant_svjdengine_h
#define nesquant_svjdengine_h

#include <ql/option.hpp>
#include <ql/PricingEngines/GenericEngine.hpp>

#include <complex>
using std::complex;

// really belongs in quantlib.h
#define QL_IMAG std::imag
#define QL_REAL std::real

// this won't work!
// #define QL_EXP std::exp
// #define QL_LOG std::log
// so we treat complex<double> seperately
#define QL_COMPLEX_EXP std::exp
#define QL_COMPLEX_LOG std::log

// this works, but issues a warning-message (macro redefinition)
// #define QL_SQRT std::sqrt
// #define QL_POW std::pow
// so we make a complex version as well
#define QL_COMPLEX_SQRT std::sqrt
#define QL_COMPLEX_POW std::pow


using namespace QuantLib;
using namespace QuantLib::PricingEngines;


namespace NesQuant {

    //! parameters for plain option calculation
    class PlainOptionParameters : public virtual Arguments {
      public:
        PlainOptionParameters() : type(Option::Type(-1)),
                                  underlying(Null<double>()),
                                  strike(Null<double>()),
                                  dividendYield(Null<double>()),
                                  riskFreeRate(Null<double>()),
                                  maturity(Null<double>()),
                                  volatility(Null<double>()) {}

        Option::Type type;
        double underlying, strike;
        Spread dividendYield;
        Rate riskFreeRate;
        Time maturity;
        double volatility;
    };

    //! parameters for stochastic volatility model (SV) (Heston 1993)
    class StochasticVolatilityArguments : public virtual Arguments {
      public:
        StochasticVolatilityArguments() :
                 volatilityOfVolatility(Null<double>()),
                 steadyStateVolatility(Null<double>()),
                 meanReversionRate(Null<double>()),
                 correlationUnderlyingVolatility(Null<double>()) {}
        double volatilityOfVolatility, steadyStateVolatility,
               meanReversionRate, correlationUnderlyingVolatility;
    };

    //! parameters for jump diffusion model (JD) (Bates 1996)
    class JumpArguments : public virtual Arguments {
      public:
        JumpArguments() : jumpIntensity(Null<double>()),
                          jumpMean(Null<double>()),
                          jumpStandardDeviation(Null<double>()) {}
        double jumpIntensity, jumpMean, jumpStandardDeviation;
    };

    //! parameters for stochastic volatility model and jumps (SVJD) (Bates 1996)
    class SVJDArguments : public PlainOptionParameters,
                          public StochasticVolatilityArguments,
                          public JumpArguments {
        void validate() const {
            //! \TODO Validation of arguments
        };
    };

    //! \brief SVJD option pricing engine for European options
    /*! Pricing engine for european options in the stochastic volatility and
        jumps model (SVJD) introduced by Bates (1996).
        Contains the Heston (1993) SV square root model as a special case. \n

        Reference: \n
        Bates, David S. (1996): \n
        Jumps and stochastic volatility:
        exchange rate processes implicit in deutsche mark options \n
        Review of Financial Studies 9: p. 69-107 \n
        http://rfs.oupjournals.org/cgi/content/abstract/9/1/69
    */

    class SVJDEngine;

    class SVJDIntegrand
    {
        public:
            SVJDIntegrand(int j, double S, double K, const SVJDEngine* engine)
                : j_(j), S_(S), K_(K), engine_(engine) {}

            double operator()(double u) const;
        private:
            int j_;
            double S_, K_;
            const SVJDEngine* engine_;
    };

    class SVJDPdfIntegrand
    {
        public:
            SVJDPdfIntegrand(double z, const SVJDEngine* engine)
                : z_(z), engine_(engine) {}

            double operator()(double u) const;
        private:
            double z_;
            const SVJDEngine* engine_;
    };

    class SVJDEngine : public GenericEngine<SVJDArguments, OptionValue> {

      public:
        void calculate() const;

        //! gives the value of the probability density function for the stock return z = log(ST / S0)
        double pdf(double z) const;

      // these two need access to F
      friend double SVJDIntegrand::operator()(double u) const;
      friend double SVJDPdfIntegrand::operator()(double u) const;

      protected:
        // helper function
        void translateArguments() const;

        // mathematical symbols to make code more readable
        mutable double S_, K_;
        mutable Rate r_;
        mutable Spread q_;
        mutable Time T_;
        mutable double V_, sigmav_, betastar_, alpha_, rho_;
        mutable double lambdastar_, kmeanstar_, delta_;

        // auxiliary variables / functions
        inline double my(int j) const
            { return (3.0 - 2.0 * j) / 2.0; }

        inline double beta(int j) const
            { return betastar_ + rho_ * sigmav_ * (j - 2); }

        complex<double> gamma(int j, complex<double> u) const;

        complex<double> C(int j, complex<double> u) const;

        complex<double> D(int j, complex<double> u) const;

        complex<double> E(int j, complex<double> u) const;

        // moment generating functions (F1 and F2)
        complex<double> F(int j, complex<double> u) const;

        // probabilities (P1 and P2)
        double P(int j) const;
    };

}

#endif