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

/*

  Copyright (C) 2014, 2015, 2018 Peter Caspers

  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.

*/

/*! \file lognormalcmsspreadpricer.cpp

*/

#include <ql/experimental/coupons/cmsspreadcoupon.hpp>

#include <ql/experimental/coupons/lognormalcmsspreadpricer.hpp>

#include <ql/math/integrals/kronrodintegral.hpp>

#include <ql/pricingengines/blackformula.hpp>

#include <ql/termstructures/volatility/swaption/swaptionvolcube.hpp>

#include <ql/optional.hpp>

#include <utility>

using std::sqrt;

namespace QuantLib {

    class LognormalCmsSpreadPricer::integrand_f {

        const LognormalCmsSpreadPricer* pricer;

      public:

        explicit integrand_f(const LognormalCmsSpreadPricer* pricer)

        : pricer(pricer) {}

        Real operator()(Real x) const {

            return pricer->integrand(x);

        }

    };

    LognormalCmsSpreadPricer::LognormalCmsSpreadPricer(

        const ext::shared_ptr<CmsCouponPricer>& cmsPricer,

        const Handle<Quote>& correlation,

        Handle<YieldTermStructure> couponDiscountCurve,

        const Size integrationPoints,

        const ext::optional<VolatilityType>& volatilityType,

        const Real shift1,

        const Real shift2)

    : CmsSpreadCouponPricer(correlation), cmsPricer_(cmsPricer),

      couponDiscountCurve_(std::move(couponDiscountCurve)) {

        registerWith(correlation);

        if (!couponDiscountCurve_.empty())

            registerWith(couponDiscountCurve_);

        registerWith(cmsPricer_);

        QL_REQUIRE(integrationPoints >= 4,

                   "at least 4 integration points should be used ("

                       << integrationPoints << ")");

        integrator_ =

            ext::make_shared<GaussHermiteIntegration>(integrationPoints);

        cnd_ = ext::make_shared<CumulativeNormalDistribution>(0.0, 1.0);

        if (!volatilityType) {

            QL_REQUIRE(shift1 == Null<Real>() && shift2 == Null<Real>(),

                       "if volatility type is inherited, no shifts should be "

                       "specified");

            inheritedVolatilityType_ = true;

            volType_ = cmsPricer->swaptionVolatility()->volatilityType();

        } else {

            shift1_ = shift1 == Null<Real>() ? 0.0 : shift1;

            shift2_ = shift2 == Null<Real>() ? 0.0 : shift2;

            inheritedVolatilityType_ = false;

            volType_ = *volatilityType;

        }

    }

Unexecuted instantiation: QuantLib::LognormalCmsSpreadPricer::LognormalCmsSpreadPricer(boost::shared_ptr<QuantLib::CmsCouponPricer> const&, QuantLib::Handle<QuantLib::Quote> const&, QuantLib::Handle<QuantLib::YieldTermStructure>, unsigned long, std::__1::optional<QuantLib::VolatilityType> const&, double, double)

Unexecuted instantiation: QuantLib::LognormalCmsSpreadPricer::LognormalCmsSpreadPricer(boost::shared_ptr<QuantLib::CmsCouponPricer> const&, QuantLib::Handle<QuantLib::Quote> const&, QuantLib::Handle<QuantLib::YieldTermStructure>, unsigned long, std::__1::optional<QuantLib::VolatilityType> const&, double, double)

    Real LognormalCmsSpreadPricer::integrand(const Real x) const {

        // this is Brigo, 13.16.2 with x = v / sqrt(2)

        Real v = M_SQRT2 * x;

        Real h =

            k_ - b_ * s2_ * std::exp((m2_ - 0.5 * v2_ * v2_) * fixingTime_ +

                                     v2_ * std::sqrt(fixingTime_) * v);

        Real phi1, phi2;

        phi1 = (*cnd_)(

            phi_ * (std::log(a_ * s1_ / h) +

                    (m1_ + (0.5 - rho_ * rho_) * v1_ * v1_) * fixingTime_ +

                    rho_ * v1_ * std::sqrt(fixingTime_) * v) /

            (v1_ * std::sqrt(fixingTime_ * (1.0 - rho_ * rho_))));

        phi2 = (*cnd_)(

            phi_ * (std::log(a_ * s1_ / h) +

                    (m1_ - 0.5 * v1_ * v1_) * fixingTime_ +

                    rho_ * v1_ * std::sqrt(fixingTime_) * v) /

            (v1_ * std::sqrt(fixingTime_ * (1.0 - rho_ * rho_))));

        Real f = a_ * phi_ * s1_ *

                     std::exp(m1_ * fixingTime_ -

                              0.5 * rho_ * rho_ * v1_ * v1_ * fixingTime_ +

                              rho_ * v1_ * std::sqrt(fixingTime_) * v) *

                     phi1 -

                 phi_ * h * phi2;

        return std::exp(-x * x) * f;

    }

    Real LognormalCmsSpreadPricer::integrand_normal(const Real x) const {

        // this is http://ssrn.com/abstract=2686998, 3.20 with x = s / sqrt(2)

        Real s = M_SQRT2 * x;

        Real beta =

            phi_ *

            (gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_ - k_ +

             std::sqrt(fixingTime_) *

                 (rho_ * gearing1_ * vol1_ + gearing2_ * vol2_) * s);

        Real f =

            close_enough(alpha_, 0.0)

                ? Real(std::max(beta, 0.0))

                : psi_ * alpha_ / (M_SQRTPI * M_SQRT2) *

                          std::exp(-beta * beta / (2.0 * alpha_ * alpha_)) +

                      beta * (1.0 - (*cnd_)(-psi_ * beta / alpha_));

        return std::exp(-x * x) * f;

    }

    void

    LognormalCmsSpreadPricer::initialize(const FloatingRateCoupon &coupon) {

        coupon_ = dynamic_cast<const CmsSpreadCoupon *>(&coupon);

        QL_REQUIRE(coupon_, "CMS spread coupon needed");

        index_ = coupon_->swapSpreadIndex();

        gearing_ = coupon_->gearing();

        spread_ = coupon_->spread();

        fixingDate_ = coupon_->fixingDate();

        paymentDate_ = coupon_->date();

        // if no coupon discount curve is given just use the discounting curve

        // from the _first_ swap index.

        // for rate calculation this curve cancels out in the computation, so

        // e.g. the discounting

        // swap engine will produce correct results, even if the

        // couponDiscountCurve is not set here.

        // only the price member function in this class will be dependent on the

        // coupon discount curve.

        today_ = QuantLib::Settings::instance().evaluationDate();

        if (couponDiscountCurve_.empty())

            couponDiscountCurve_ =

                index_->swapIndex1()->exogenousDiscount()

                    ? index_->swapIndex1()->discountingTermStructure()

                    : index_->swapIndex1()->forwardingTermStructure();

        discount_ = paymentDate_ > couponDiscountCurve_->referenceDate()

                        ? couponDiscountCurve_->discount(paymentDate_)

                        : 1.0;

        spreadLegValue_ = spread_ * coupon_->accrualPeriod() * discount_;

        gearing1_ = index_->gearing1();

        gearing2_ = index_->gearing2();

        QL_REQUIRE(gearing1_ > 0.0 && gearing2_ < 0.0,

                   "gearing1 (" << gearing1_

                                << ") should be positive while gearing2 ("

                                << gearing2_ << ") should be negative");

        c1_ = ext::make_shared<CmsCoupon>(

            coupon_->date(), coupon_->nominal(), coupon_->accrualStartDate(),

            coupon_->accrualEndDate(), coupon_->fixingDays(),

            index_->swapIndex1(), 1.0, 0.0, coupon_->referencePeriodStart(),

            coupon_->referencePeriodEnd(), coupon_->dayCounter(),

            coupon_->isInArrears());

        c2_ = ext::make_shared<CmsCoupon>(

            coupon_->date(), coupon_->nominal(), coupon_->accrualStartDate(),

            coupon_->accrualEndDate(), coupon_->fixingDays(),

            index_->swapIndex2(), 1.0, 0.0, coupon_->referencePeriodStart(),

            coupon_->referencePeriodEnd(), coupon_->dayCounter(),

            coupon_->isInArrears());

        c1_->setPricer(cmsPricer_);

        c2_->setPricer(cmsPricer_);

        if (fixingDate_ > today_) {

            fixingTime_ = cmsPricer_->swaptionVolatility()->timeFromReference(

                fixingDate_);

            swapRate1_ = c1_->indexFixing();

            swapRate2_ = c2_->indexFixing();

            adjustedRate1_ = c1_->adjustedFixing();

            adjustedRate2_ = c2_->adjustedFixing();

            ext::shared_ptr<SwaptionVolatilityStructure> swvol =

                *cmsPricer_->swaptionVolatility();

            ext::shared_ptr<SwaptionVolatilityCube> swcub =

                ext::dynamic_pointer_cast<SwaptionVolatilityCube>(swvol);

            if(inheritedVolatilityType_ && volType_ == ShiftedLognormal) {

                shift1_ =

                    swvol->shift(fixingDate_, index_->swapIndex1()->tenor());

                shift2_ =

                    swvol->shift(fixingDate_, index_->swapIndex2()->tenor());

            }

            if (swcub == nullptr) {

                // not a cube, just an atm surface given, so we can

                // not easily convert volatilities and just forbid it

                QL_REQUIRE(inheritedVolatilityType_,

                           "if only an atm surface is given, the volatility "

                           "type must be inherited");

                vol1_ = swvol->volatility(

                    fixingDate_, index_->swapIndex1()->tenor(), swapRate1_);

                vol2_ = swvol->volatility(

                    fixingDate_, index_->swapIndex2()->tenor(), swapRate2_);

            } else {

                vol1_ = swcub->smileSection(fixingDate_,

                                            index_->swapIndex1()->tenor())

                            ->volatility(swapRate1_, volType_, shift1_);

                vol2_ = swcub->smileSection(fixingDate_,

                                            index_->swapIndex2()->tenor())

                            ->volatility(swapRate2_, volType_, shift2_);

            }

            if(volType_ == ShiftedLognormal) {

                mu1_ = 1.0 / fixingTime_ * std::log((adjustedRate1_ + shift1_) /

                                                    (swapRate1_ + shift1_));

                mu2_ = 1.0 / fixingTime_ * std::log((adjustedRate2_ + shift2_) /

                                                    (swapRate2_ + shift2_));

            }

            // for the normal volatility case we do not need the drifts

            // but rather use adjusted rates directly in the integrand

            rho_ = std::max(std::min(correlation()->value(), 0.9999),

                            -0.9999); // avoid division by zero in integrand

        } else {

            // fixing is in the past or today

            adjustedRate1_ = c1_->indexFixing();

            adjustedRate2_ = c2_->indexFixing();

        }

    }

    Real LognormalCmsSpreadPricer::optionletPrice(Option::Type optionType,

                                                  Real strike) const {

        // this method is only called for future fixings

        optionType_ = optionType;

        phi_ = optionType == Option::Call ? 1.0 : -1.0;

        Real res = 0.0;

        if (volType_ == ShiftedLognormal) {

            // (shifted) lognormal volatility

            if (strike >= 0.0) {

                a_ = gearing1_;

                b_ = gearing2_;

                s1_ = swapRate1_ + shift1_;

                s2_ = swapRate2_ + shift2_;

                m1_ = mu1_;

                m2_ = mu2_;

                v1_ = vol1_;

                v2_ = vol2_;

                k_ = strike + gearing1_ * shift1_ + gearing2_ * shift2_;

            } else {

                a_ = -gearing2_;

                b_ = -gearing1_;

                s1_ = swapRate2_ + shift1_;

                s2_ = swapRate1_ + shift2_;

                m1_ = mu2_;

                m2_ = mu1_;

                v1_ = vol2_;

                v2_ = vol1_;

                k_ = -strike - gearing1_ * shift1_ - gearing2_ * shift2_;

                res += phi_ * (gearing1_ * adjustedRate1_ +

                               gearing2_ * adjustedRate2_ - strike);

            }

            res +=

                1.0 / M_SQRTPI * (*integrator_)(integrand_f(this));

        } else {

            // normal volatility

            Real forward = gearing1_ * adjustedRate1_ +

                gearing2_ * adjustedRate2_;

            Real stddev =

                std::sqrt(fixingTime_ *

                          (gearing1_ * gearing1_ * vol1_ * vol1_ +

                           gearing2_ * gearing2_ * vol2_ * vol2_ +

                           2.0 * gearing1_ * gearing2_ * rho_ * vol1_ * vol2_));

            res =

                bachelierBlackFormula(optionType_, strike, forward, stddev, 1.0);

        }

        return res * discount_ * coupon_->accrualPeriod();

    }

    Rate LognormalCmsSpreadPricer::swapletRate() const {

        return swapletPrice() / (coupon_->accrualPeriod() * discount_);

    }

    Real LognormalCmsSpreadPricer::capletPrice(Rate effectiveCap) const {

        // caplet is equivalent to call option on fixing

        if (fixingDate_ <= today_) {

            // the fixing is determined

            const Rate Rs = std::max(

                coupon_->index()->fixing(fixingDate_) - effectiveCap, 0.);

            Rate price = gearing_ * Rs * coupon_->accrualPeriod() * discount_;

            return price;

        } else {

            Real capletPrice = optionletPrice(Option::Call, effectiveCap);

            return gearing_ * capletPrice;

        }

    }

    Rate LognormalCmsSpreadPricer::capletRate(Rate effectiveCap) const {

        return capletPrice(effectiveCap) /

               (coupon_->accrualPeriod() * discount_);

    }

    Real LognormalCmsSpreadPricer::floorletPrice(Rate effectiveFloor) const {

        // floorlet is equivalent to put option on fixing

        if (fixingDate_ <= today_) {

            // the fixing is determined

            const Rate Rs = std::max(

                effectiveFloor - coupon_->index()->fixing(fixingDate_), 0.);

            Rate price = gearing_ * Rs * coupon_->accrualPeriod() * discount_;

            return price;

        } else {

            Real floorletPrice = optionletPrice(Option::Put, effectiveFloor);

            return gearing_ * floorletPrice;

        }

    }

    Rate LognormalCmsSpreadPricer::floorletRate(Rate effectiveFloor) const {

        return floorletPrice(effectiveFloor) /

               (coupon_->accrualPeriod() * discount_);

    }

    Real LognormalCmsSpreadPricer::swapletPrice() const {

        return gearing_ * coupon_->accrualPeriod() * discount_ *

                   (gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_) +

               spreadLegValue_;

    }

}