/*
 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.cpp
    \brief SVJD pricing engine for European options

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

#include <nq/Math/KronrodIntegral.hpp>
#include <nq/PricingEngines/svjdengine.hpp>

namespace NesQuant {

    namespace PricingEngines {

        void SVJDEngine::calculate() const {
            // translate argument names to mathematical names
            S_ = arguments_.underlying;
            K_ = arguments_.strike;
            r_ = arguments_.riskFreeRate;
            q_ = arguments_.dividendYield;
            T_ = arguments_.residualTime;
            V_ = arguments_.volatility * arguments_.volatility;
            sigmav_ = arguments_.volatilityOfVolatility;
            betastar_ = arguments_.meanReversionRate;
            alpha_ = (arguments_.steadyStateVolatility
                        * arguments_.steadyStateVolatility)
                     * betastar_;
            rho_ = arguments_.correlationUnderlyingVolatility;
            lambdastar_ = arguments_.jumpIntensity;
            kmeanstar_ = arguments_.jumpMean;
            delta_ = arguments_.jumpStandardDeviation;

            // calculate discount factors
            DiscountFactor dividendDiscount = QL_EXP(-q_ * T_);
            DiscountFactor riskFreeDiscount = QL_EXP(-r_ * T_);

            // calculate probabilities
            double MP1, MP2; // modified P1, P2

            switch (arguments_.type) {
              case Option::Call:
                MP1 = P(1);
                MP2 = P(2);
                break;
              case Option::Put:
                MP1 = P(1) - 1.0;
                MP2 = P(2) - 1.0;
                break;
              case Option::Straddle:
                MP1 = 2.0 * P(1) - 1.0;
                MP2 = 2.0 * P(2) - 1.0;
                break;
              default:
                throw IllegalArgumentError("SVJDEngine: invalid option type");
            }

            // calculate option price
            results_.value = S_ * dividendDiscount * MP1
                             - K_ * riskFreeDiscount * MP2;
        }

        complex<double> SVJDEngine::gamma(int j, complex<double> u) const {
            return QL_COMPLEX_SQRT( (rho_ * sigmav_ * u - beta(j))
                                      * (rho_ * sigmav_ * u - beta(j))
                                    - 2.0 * sigmav_ * sigmav_
                                      * (my(j) * u + (u * u) / 2.0));
        }

        complex<double> SVJDEngine::C(int j, complex<double> u) const {
            complex<double> expgammajT = QL_COMPLEX_EXP(gamma(j, u) * T_);

            return (r_ - q_ - lambdastar_ * kmeanstar_) * u * T_
                   - ((alpha_ * T_) / (sigmav_ * sigmav_))
                     * (rho_ * sigmav_ * u - beta(j) - gamma(j, u))
                   - ((2.0 * alpha_) / (sigmav_ * sigmav_))
           * QL_COMPLEX_LOG(1.0 + (rho_ * sigmav_ * u - beta(j) - gamma(j, u))
                   / 2.0
                     * ((1.0 - expgammajT) / gamma(j, u)));
        }

        complex<double> SVJDEngine::D(int j, complex<double> u) const {
            complex<double> expgammajT = QL_COMPLEX_EXP(gamma(j, u) * T_);

            return -2.0 * (my(j) * u + (u * u) / 2.0)
                      / (rho_ * sigmav_ * u - beta(j)
                         + gamma(j, u) * (1.0 + expgammajT)
                           / (1.0 - expgammajT));
        }

        complex<double> SVJDEngine::E(int j, complex<double> u) const {
            double temp1;
            complex<double> temp2;

            temp1 = pow(1.0 + kmeanstar_, my(j) + 0.5);
            temp2 = QL_COMPLEX_POW(1.0 + kmeanstar_, u);

            return lambdastar_ * T_ * temp1
                   * (temp2 * QL_COMPLEX_EXP(
                        delta_ * delta_ * (my(j) * u + (u * u) / 2.0))
                     - 1.0);
        }

        complex<double> SVJDEngine::F(int j, complex<double> u) const
        {
            return QL_COMPLEX_EXP(C(j, u) + D(j, u) * V_ + E(j, u));
        }

        double SVJDIntegrand::operator()(double u) const {
            const complex<double> I(0,1);
            double temp = QL_LOG(K_/S_);
            return QL_IMAG(engine_->F(j_, I * u)
                   * QL_COMPLEX_EXP(-I * u * temp))/u;
        }

        double SVJDEngine::P(int j) const {
            using NesQuant::Math::KronrodIntegral;
            KronrodIntegral integrator(0.00001);

            SVJDIntegrand integrand(j, S_, K_, this);

            double integral = integrator(integrand, 0, 300);

            return 0.5 + M_1_PI * integral;
        }

        double SVJDPdfIntegrand::operator()(double u) const {
            const complex<double> I(0,1);
            return QL_REAL(engine_->F(2, I * u)
                   *QL_COMPLEX_EXP(-I * u * z_));
        }

        double SVJDEngine::pdf(double z) const {
            using NesQuant::Math::KronrodIntegral;
            KronrodIntegral integrator(0.00001);

            SVJDPdfIntegrand integrand(z, this);

            double integral = integrator(integrand, 0, 300);

            return M_1_PI * integral;
        }

    }

}