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

/*

 Copyright (C) 2009 Dimitri Reiswich

 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/pricingengines/exotic/analyticcompoundoptionengine.hpp>

#include <ql/math/solvers1d/brent.hpp>

#include <ql/pricingengines/blackformula.hpp>

#include <utility>

namespace QuantLib {

    namespace {

        // Helper Class needed to solve an implicit problem of finding a

        // spot to a corresponding option price.

        class ImpliedSpotHelper {

          public:

            ImpliedSpotHelper(DiscountFactor dividendDiscount,

                              DiscountFactor riskFreeDiscount,

                              Real standardDeviation,

                              ext::shared_ptr<PlainVanillaPayoff> payoff,

                              Real strike)

            : dividendDiscount_(dividendDiscount), riskFreeDiscount_(riskFreeDiscount),

              standardDeviation_(standardDeviation), strike_(strike), payoff_(std::move(payoff)) {}

            Real operator()(Real spot) const {

                Real forwardPrice = spot*dividendDiscount_/riskFreeDiscount_;

                Real value = blackFormula(payoff_, forwardPrice,

                                          standardDeviation_,riskFreeDiscount_);

                return value - strike_;

            }

          private:

            DiscountFactor dividendDiscount_;

            DiscountFactor riskFreeDiscount_;

            Real standardDeviation_;

            Real strike_;

            ext::shared_ptr<PlainVanillaPayoff> payoff_;

        };

    }

    AnalyticCompoundOptionEngine::AnalyticCompoundOptionEngine(

        ext::shared_ptr<GeneralizedBlackScholesProcess> process)

    : process_(std::move(process)) {

        registerWith(process_);

    }

    void AnalyticCompoundOptionEngine::calculate() const {

        QL_REQUIRE(strikeDaughter()>0.0,

                   "Daughter strike must be positive");

        QL_REQUIRE(strikeMother()>0.0,

                   "Mother strike must be positive");

        QL_REQUIRE(spot() > 0.0, "negative or null underlying given");

        /* Solver Setup ***************************************************/

        Date helpDate(process_->riskFreeRate()->referenceDate());

        Date helpMaturity=helpDate+(maturityDaughter()-maturityMother())*Days;

        Real vol =process_->blackVolatility()->blackVol(helpMaturity,

                                                        strikeDaughter());

        Time helpTimeToMat=process_->time(helpMaturity);

        vol=vol*std::sqrt(helpTimeToMat);

        DiscountFactor dividendDiscount =

            process_->dividendYield()->discount(helpMaturity);

        DiscountFactor riskFreeDiscount =

            process_->riskFreeRate()->discount(helpMaturity);

        auto f = ext::make_shared<ImpliedSpotHelper>(dividendDiscount, riskFreeDiscount,

                                      vol, payoffDaughter(), strikeMother());

        Brent solver;

        solver.setMaxEvaluations(1000);

        Real accuracy = 1.0e-6;

        Real sSolved=solver.solve(*f, accuracy, strikeDaughter(), 1.0e-6, strikeDaughter()*1000.0);

        Real X=transformX(sSolved); // transform stock to return as in Wystup's book

        /* Solver Setup Finished*****************************************/

        Real phi=typeDaughter(); // -1 or 1

        Real w=typeMother(); // -1 or 1

        Real rho=std::sqrt(residualTimeMother()/residualTimeDaughter());

        BivariateCumulativeNormalDistributionDr78 N2(w*rho) ;

        DiscountFactor ddD=dividendDiscountDaughter();

        DiscountFactor rdD=riskFreeDiscountDaughter();

        //DiscountFactor ddM=dividendDiscountMother();

        DiscountFactor rdM=riskFreeDiscountMother();

        Real XmSM=X-stdDeviationMother();

        Real S=spot();

        Real dP=dPlus();

        Real dPT12=dPlusTau12(sSolved);

        Real vD=volatilityDaughter();

        Real dM=dMinus();

        Real strD=strikeDaughter();

        Real strM=strikeMother();

        Real rTM=residualTimeMother();

        Real rTD=residualTimeDaughter();

        Real rD=riskFreeRateDaughter();

        Real dD=dividendRateDaughter();

        Real N2XmSM=N2(-phi*w*XmSM,phi*dP);

        Real N2X=N2(-phi*w*X,phi*dM);

        Real NeX=N_(-phi*w*e(X));

        Real NX=N_(-phi*w*X);

        Real NT12=N_(phi*dPT12);

        Real ndP=n_(dP);

        Real nXm=n_(XmSM);

        Real invMTime=1/std::sqrt(rTM);

        Real invDTime=1/std::sqrt(rTD);

        Real tempRes=phi*w*S*ddD*N2XmSM-phi*w*strD*rdD*N2X-w*strM*rdM*NX;

        Real tempDelta=phi*w*ddD*N2XmSM;

        Real tempGamma=(ddD/(vD*S))*(invMTime*nXm*NT12+w*invDTime*ndP*NeX);

        Real tempVega=ddD*S*((1/invMTime)*nXm*NT12+w*(1/invDTime)*ndP*NeX);

        Real tempTheta=phi*w*dD*S*ddD*N2XmSM-phi*w*rD*strD*rdD*N2X-w*rD*strM*rdM*NX;

        tempTheta-=0.5*vD*S*ddD*(invMTime*nXm*NT12+w*invDTime*ndP*NeX);

        results_.value=tempRes;

        results_.delta=tempDelta;

        results_.gamma=tempGamma;

        results_.vega=tempVega;

        results_.theta=tempTheta;

    }

    Real AnalyticCompoundOptionEngine::typeDaughter() const {

        // returns -1 or 1 according to put or call

        return (Real) payoffDaughter()->optionType();

    }

    Real AnalyticCompoundOptionEngine::typeMother() const {

        return (Real) payoffMother()->optionType();

    }

    Date AnalyticCompoundOptionEngine::maturityDaughter() const {

        return arguments_.daughterExercise->lastDate();

    }

    Date AnalyticCompoundOptionEngine::maturityMother() const {

        return arguments_.exercise->lastDate();

    }

    Time AnalyticCompoundOptionEngine::residualTimeDaughter() const {

        return process_->time(maturityDaughter());

    }

    Time AnalyticCompoundOptionEngine::residualTimeMother() const {

        return process_->time(maturityMother());

    }

    Time AnalyticCompoundOptionEngine::residualTimeMotherDaughter() const {

        return residualTimeDaughter()-residualTimeMother();

    }

    Real AnalyticCompoundOptionEngine::volatilityDaughter() const {

        return process_->blackVolatility()->blackVol(maturityDaughter(),

                                                     strikeDaughter());

    }

    Real AnalyticCompoundOptionEngine::volatilityMother() const {

        return process_->blackVolatility()->blackVol(maturityMother(),

                                                     strikeMother());

    }

    Real AnalyticCompoundOptionEngine::stdDeviationDaughter() const {

        return volatilityDaughter()*std::sqrt(residualTimeDaughter());

    }

    Real AnalyticCompoundOptionEngine::stdDeviationMother() const {

        return volatilityMother()*std::sqrt(residualTimeMother());

    }

    ext::shared_ptr<PlainVanillaPayoff>

    AnalyticCompoundOptionEngine::payoffDaughter() const {

        ext::shared_ptr<PlainVanillaPayoff> dPayoff =

            ext::dynamic_pointer_cast<PlainVanillaPayoff>(

                                                   arguments_.daughterPayoff);

        QL_REQUIRE(dPayoff, "non-plain payoff given");

        return dPayoff;

    }

    ext::shared_ptr<PlainVanillaPayoff>

    AnalyticCompoundOptionEngine::payoffMother() const {

        ext::shared_ptr<PlainVanillaPayoff> mPayoff =

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

        QL_REQUIRE(mPayoff, "non-plain payoff given");

        return mPayoff;

    }

    Real AnalyticCompoundOptionEngine::strikeMother() const {

        return payoffMother()->strike();

    }

    Real AnalyticCompoundOptionEngine::strikeDaughter() const {

        return payoffDaughter()->strike();

    }

    DiscountFactor AnalyticCompoundOptionEngine::riskFreeDiscountDaughter() const {

        return process_->riskFreeRate()->discount(residualTimeDaughter());

    }

    DiscountFactor AnalyticCompoundOptionEngine::riskFreeDiscountMother() const {

        return process_->riskFreeRate()->discount(residualTimeMother());

    }

    DiscountFactor AnalyticCompoundOptionEngine::riskFreeDiscountMotherDaughter() const {

        return process_->riskFreeRate()->discount(residualTimeMotherDaughter());

    }

    DiscountFactor AnalyticCompoundOptionEngine::dividendDiscountDaughter() const {

        return process_->dividendYield()->discount(residualTimeDaughter());

    }

    DiscountFactor AnalyticCompoundOptionEngine::dividendDiscountMother() const {

        return process_->dividendYield()->discount(residualTimeMother());

    }

    DiscountFactor AnalyticCompoundOptionEngine::dividendDiscountMotherDaughter() const {

        return process_->dividendYield()->discount(residualTimeMotherDaughter());

    }

    Real AnalyticCompoundOptionEngine::dPlus() const {

        Real forward = spot() * dividendDiscountDaughter() / riskFreeDiscountDaughter();

        Real sd=stdDeviationDaughter();

        return std::log(forward/strikeDaughter())/sd+0.5*sd;

    }

    Real AnalyticCompoundOptionEngine::dMinus() const {

        return dPlus()-stdDeviationDaughter();

    }

    Real AnalyticCompoundOptionEngine::dPlusTau12(Real S) const {

        Real forward = S * dividendDiscountMotherDaughter() / riskFreeDiscountMotherDaughter();

        Real sd=volatilityDaughter()*std::sqrt(residualTimeMotherDaughter());

        return std::log(forward/strikeDaughter())/sd+0.5*sd;

    }

    Real AnalyticCompoundOptionEngine::spot() const {

        return process_->x0();

    }

    Real AnalyticCompoundOptionEngine::riskFreeRateDaughter() const {

        return process_->riskFreeRate()->zeroRate(residualTimeDaughter(),

                                                  Continuous,

                                                  NoFrequency);

    }

    Real AnalyticCompoundOptionEngine::dividendRateDaughter() const {

        return process_->dividendYield()->zeroRate(residualTimeDaughter(),

                                                   Continuous,

                                                   NoFrequency);

    }

    Real AnalyticCompoundOptionEngine::transformX(Real X) const {

        Real sd=stdDeviationMother();

        Real resX=riskFreeDiscountMother()*X/(spot()*dividendDiscountMother());

        resX=resX*std::exp(0.5*sd*sd);

        resX=std::log(resX);

        return resX/sd;

    }

    Real AnalyticCompoundOptionEngine::e(Real X) const {

        Real rtM=residualTimeMother();

        Real rtD=residualTimeDaughter();

        return (X*std::sqrt(rtD)+std::sqrt(rtM)*dMinus())/std::sqrt(rtD-rtM);

    }

}