/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*

 Copyright (C) 2020 Jack Gillett

 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

 <quantlib-dev@lists.sf.net>. The license is also available online at

 <https://www.quantlib.org/license.shtml>.

 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 <ql/experimental/forward/analytichestonforwardeuropeanengine.hpp>

#include <complex>

#include <utility>

namespace QuantLib {

    class P12Integrand {

      private:

        ext::shared_ptr<AnalyticHestonEngine>& engine_;

        Real logK_, phiRightLimit_;

        Time tenor_;

        std::complex<Real> i_, adj_;

      public:

        P12Integrand(ext::shared_ptr<AnalyticHestonEngine>& engine,

                     Real logK,

                     Time tenor,

                     bool P1, // true for P1, false for P2

                     Real phiRightLimit = 100) : engine_(engine), logK_(logK),

            phiRightLimit_(phiRightLimit), tenor_(tenor), i_(std::complex<Real>(0.0, 1.0)) {

            // Only difference between P1 and P2 integral is the additional term in the chF evaluation

            if (P1) {

                adj_ = std::complex<Real>(0.0, -1.0);

            } else {

                adj_ = std::complex<Real>(0.0, 0.0);

            }

        }

        // QL Gaussian Quadrature - map phi from [-1 to 1] to {0, phiRightLimit] 

        Real operator()(Real phi) const {

            Real phiDash = (0.5+1e-8+0.5*phi) * phiRightLimit_; // Map phi to full range

            return 0.5*phiRightLimit_*std::real((std::exp(-phiDash*logK_*i_) / (phiDash*i_)) * engine_->chF(phiDash+adj_, tenor_));

        }

    };

    class P12HatIntegrand {

      private:

        Time tenor_, resetTime_;

        Handle<Quote>& s0_;

        bool P1_;

        Real logK_, phiRightLimit_, nuRightLimit_;

        const AnalyticHestonForwardEuropeanEngine* const parent_;

        GaussLegendreIntegration innerIntegrator_;

      public:

        P12HatIntegrand(Time tenor,

                        Time resetTime,

                        Handle<Quote>& s0,

                        Real logK,

                        bool P1, // true for P1, false for P2

                        const AnalyticHestonForwardEuropeanEngine* const parent,

                        Real phiRightLimit,

                        Real nuRightLimit) : tenor_(tenor), resetTime_(resetTime),

            s0_(s0), P1_(P1), logK_(logK), phiRightLimit_(phiRightLimit),

            nuRightLimit_(nuRightLimit), parent_(parent), innerIntegrator_(128) {}

        Real operator()(Real nu) const {

            // Rescale nu to [-1, 1]

            Real nuDash = nuRightLimit_ * (0.5 * nu + 0.5 + 1e-8);

            // Calculate the chF from var(t) to expiry

            ext::shared_ptr<AnalyticHestonEngine> engine = parent_->forwardChF(s0_, nuDash);

            P12Integrand pIntegrand(engine, logK_, tenor_, P1_, phiRightLimit_);

            Real p1Integral = innerIntegrator_(pIntegrand);

            // Calculate the value of the propagator to nu

            Real propagator = parent_->propagator(resetTime_, nuDash);

            // Take the product, and scale integral's value back up to [0, right_lim]

            return propagator * (0.5 + p1Integral/M_PI);

        }

    };

    AnalyticHestonForwardEuropeanEngine::AnalyticHestonForwardEuropeanEngine(

        ext::shared_ptr<HestonProcess> process, Size integrationOrder)

    : process_(std::move(process)), integrationOrder_(integrationOrder), outerIntegrator_(128) {

        v0_ = process_->v0();

        rho_ = process_->rho();

        kappa_ = process_->kappa();

        theta_ = process_->theta();

        sigma_ = process_->sigma();

        s0_ = process_->s0();

        QL_REQUIRE(sigma_ > 0.1,

                   "Very low values (<~10%) for Heston Vol-of-Vol cause numerical issues" \

                   "in this implementation of the propagator function, try using" \

                   "MCForwardEuropeanHestonEngine Monte-Carlo engine instead");

        riskFreeRate_ = process_->riskFreeRate();

        dividendYield_ = process_->dividendYield();

        // Some of the required constant intermediate variables can be calculated now

        kappaHat_ = kappa_ - rho_ * sigma_;

        thetaHat_ = kappa_ * theta_ / kappaHat_;

        R_ = 4 * kappaHat_ * thetaHat_ / (sigma_ * sigma_);

    }

    void AnalyticHestonForwardEuropeanEngine::calculate() const {

        // This is a european option pricer

        QL_REQUIRE(this->arguments_.exercise->type() == Exercise::European,

                   "not an European option");

        // We only price plain vanillas

        ext::shared_ptr<PlainVanillaPayoff> payoff =

            ext::dynamic_pointer_cast<PlainVanillaPayoff>(this->arguments_.payoff);

        QL_REQUIRE(payoff, "non plain vanilla payoff given");

        Time resetTime = this->process_->time(this->arguments_.resetDate);

        Time expiryTime = this->process_->time(this->arguments_.exercise->lastDate());

        Time tenor = expiryTime - resetTime;

        Real moneyness = this->arguments_.moneyness;

        // K needs to be scaled to forward AT RESET TIME, not spot...

        Real expiryDcf = riskFreeRate_->discount(expiryTime);

        Real resetDcf = riskFreeRate_->discount(resetTime);

        Real expiryDividendDiscount = dividendYield_->discount(expiryTime);

        Real resetDividendDiscount = dividendYield_->discount(resetTime);

        Real expiryRatio = expiryDcf / expiryDividendDiscount;

        Real resetRatio = resetDcf / resetDividendDiscount;

        QL_REQUIRE(resetTime >= 0.0, "Reset Date cannot be in the past");

        QL_REQUIRE(expiryTime >= 0.0, "Expiry Date cannot be in the past");

        // Use some heuristics to decide upon phiRightLimit and nuRightLimit

        Real phiRightLimit = 100.0;

        Real nuRightLimit = std::max(2.0, 10.0 * (1+std::max(0.0, rho_)) * sigma_ * std::sqrt(resetTime * std::max(v0_, theta_)));

        // do the 2D integral calculation. For very short times, we just fall back on the standard

        // calculation, both for accuracy and because tStar==0 causes some numerical issues...

        std::pair<Real, Real> P1HatP2Hat;

        if (resetTime <= 1e-3) {

            Handle<Quote> tempQuote(ext::shared_ptr<Quote>(new SimpleQuote(s0_->value())));

            P1HatP2Hat = calculateP1P2(tenor, tempQuote, moneyness * s0_->value(), expiryRatio, phiRightLimit);

        } else {

            P1HatP2Hat = calculateP1P2Hat(tenor, resetTime, moneyness, expiryRatio/resetRatio, phiRightLimit, nuRightLimit);

        }

        // Apply the payoff functions

        Real value = 0.0;

        Real F = s0_->value() / expiryRatio;

        switch (payoff->optionType()){

            case Option::Call:

                value = expiryDcf * (F*P1HatP2Hat.first - moneyness*s0_->value()*P1HatP2Hat.second/resetRatio);

                break;

            case Option::Put:

                value = expiryDcf * (moneyness*s0_->value()*(1-P1HatP2Hat.second)/resetRatio - F*(1-P1HatP2Hat.first));

                break;

            default:

                QL_FAIL("unknown option type");

            }

        results_.value = value;

        results_.additionalResults["dcf"] = expiryDcf;

        results_.additionalResults["qf"] = expiryDividendDiscount;

        results_.additionalResults["expiryRatio"] = expiryRatio;

        results_.additionalResults["resetRatio"] = resetRatio;

        results_.additionalResults["moneyness"] = moneyness;

        results_.additionalResults["s0"] = s0_->value();

        results_.additionalResults["fwd"] = F;

        results_.additionalResults["resetTime"] = resetTime;

        results_.additionalResults["expiryTime"] = expiryTime;

        results_.additionalResults["P1Hat"] = P1HatP2Hat.first;

        results_.additionalResults["P2Hat"] = P1HatP2Hat.second;

        results_.additionalResults["phiRightLimit"] = phiRightLimit;

        results_.additionalResults["nuRightLimit"] = nuRightLimit;

    }

    std::pair<Real, Real> AnalyticHestonForwardEuropeanEngine::calculateP1P2Hat(Time tenor,

                                                                                Time resetTime,

                                                                                Real moneyness,

                                                                                Real ratio,

                                                                                Real phiRightLimit,

                                                                                Real nuRightLimit) const {

        Handle<Quote> unitQuote(ext::shared_ptr<Quote>(new SimpleQuote(1.0)));

        // Re-expressing moneyness in terms of the forward here (strike fixes to spot, but in

        // our pricing calculation we need to compare it to the future at expiry)

        Real logMoneyness = std::log(moneyness*ratio);

        P12HatIntegrand p1HatIntegrand(tenor, resetTime, unitQuote, logMoneyness, true, this, phiRightLimit, nuRightLimit);

        P12HatIntegrand p2HatIntegrand(tenor, resetTime, unitQuote, logMoneyness, false, this, phiRightLimit, nuRightLimit);

        Real p1HatIntegral = 0.5 * nuRightLimit * outerIntegrator_(p1HatIntegrand);

        Real p2HatIntegral = 0.5 * nuRightLimit * outerIntegrator_(p2HatIntegrand);

        std::pair<Real, Real> P1HatP2Hat(p1HatIntegral, p2HatIntegral);

        return P1HatP2Hat;

    }

    Real AnalyticHestonForwardEuropeanEngine::propagator(Time resetTime,

                                                         Real varReset) const {

        Real B, Lambda, term1, term2, term3;

        B = 4 * kappaHat_ / (sigma_ * sigma_ * (1 - std::exp(-kappaHat_ * resetTime)));

        Lambda = B * std::exp(-kappaHat_ * resetTime) * v0_;

        // Now construct equation (18) from the paper term-by-term

        term1 = std::exp(-0.5*(B * varReset + Lambda)) * B / 2;

        term2 = std::pow(B * varReset / Lambda, 0.5*(R_/2 - 1));

        term3 = modifiedBesselFunction_i(Real(R_/2 - 1),Real(std::sqrt(Lambda * B * varReset)));

        return term1 * term2 * term3;

    }

    ext::shared_ptr<AnalyticHestonEngine> AnalyticHestonForwardEuropeanEngine::forwardChF(

                                      Handle<Quote>& spotReset,

                                      Real varReset) const {

        // Probably a wasteful implementation here, could be improved by importing

        // only the CF-generating parts of the AnalyticHestonEngine (currently private)

        ext::shared_ptr<HestonProcess> hestonProcess(new

            HestonProcess(riskFreeRate_, dividendYield_, spotReset,

                varReset, kappa_, theta_, sigma_, rho_));

        ext::shared_ptr<HestonModel> hestonModel(new HestonModel(hestonProcess));

        ext::shared_ptr<AnalyticHestonEngine> analyticHestonEngine(

            new AnalyticHestonEngine(hestonModel, integrationOrder_));

        // Not sure how to pass only the chF, so just pass the whole thing for now!

        return analyticHestonEngine;

    }

    std::pair<Real, Real> AnalyticHestonForwardEuropeanEngine::calculateP1P2(Time tenor,

                                                                             Handle<Quote>& St,

                                                                             Real K,

                                                                             Real ratio,

                                                                             Real phiRightLimit) const {

        ext::shared_ptr<AnalyticHestonEngine> engine = forwardChF(St, v0_);

        Real logK = std::log(K*ratio/St->value());

        // Integrate the CF and the complex integrand over positive phi

        GaussLegendreIntegration integrator = GaussLegendreIntegration(128);

        P12Integrand p1Integrand(engine, logK, tenor, true, phiRightLimit);

        P12Integrand p2Integrand(engine, logK, tenor, false, phiRightLimit);

        Real p1Integral = integrator(p1Integrand);

        Real p2Integral = integrator(p2Integrand);

        std::pair<Real, Real> P1P2(0.5 + p1Integral/M_PI, 0.5 + p2Integral/M_PI);

        return P1P2;

    }

}