/* -*- 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
 <http://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/compoundoption/analyticcompoundoptionengine.hpp>

namespace QuantLib {

    AnalyticCompoundOptionEngine::AnalyticCompoundOptionEngine(
            const boost::shared_ptr<GeneralizedBlackScholesProcess>& process)
    : process_(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);


        boost::shared_ptr<ImpliedSpotHelper> f(
                new ImpliedSpotHelper(dividendDiscount, riskFreeDiscount,
                                      vol, payoffDaughter(), strikeMother()));

        Brent solver;
        solver.setMaxEvaluations(1000);
        Real accuracy = 1.0e-6;

        Real X=0.0;
        Real sSolved=0.0;
        sSolved=solver.solve(*f, accuracy, strikeDaughter(), 1.0e-6, strikeDaughter()*1000.0);
        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 tempRes=0.0;
        Real tempDelta=0.0;
        Real tempGamma=0.0;
        Real tempVega=0.0;
        Real tempTheta=0.0;
        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);

        tempRes=phi*w*S*ddD*N2XmSM-phi*w*strD*rdD*N2X-w*strM*rdM*NX;
        tempDelta=phi*w*ddD*N2XmSM;
        tempGamma=(ddD/(vD*S))*(invMTime*nXm*NT12+w*invDTime*ndP*NeX);
        tempVega=ddD*S*((1/invMTime)*nXm*NT12+w*(1/invDTime)*ndP*NeX);

        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_.exercise->lastDate();
    }

    Date AnalyticCompoundOptionEngine::maturityMother() const{
        return (arguments_.motherOption->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());
    }


    boost::shared_ptr<PlainVanillaPayoff>
    AnalyticCompoundOptionEngine::payoffDaughter() const{
        boost::shared_ptr<PlainVanillaPayoff> dPayoff =
            boost::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);

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

        return dPayoff;
    }

    boost::shared_ptr<PlainVanillaPayoff>
    AnalyticCompoundOptionEngine::payoffMother() const{

        boost::shared_ptr<PlainVanillaPayoff> mPayoff =
            boost::dynamic_pointer_cast<PlainVanillaPayoff>(
                                         (arguments_.motherOption)->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);
    }


    ImpliedSpotHelper::ImpliedSpotHelper(
                                 DiscountFactor dividendDiscount,
                                 DiscountFactor riskFreeDiscount,
                                 Real standardDeviation,
                                 boost::shared_ptr<PlainVanillaPayoff> payoff,
                                 Real strike)
    : dividendDiscount_(dividendDiscount), riskFreeDiscount_(riskFreeDiscount),
      standardDeviation_(standardDeviation),
      strike_(strike),payoff_(payoff){}

    Real ImpliedSpotHelper::operator ()(Real spot)const{

        Real forwardPrice = spot * dividendDiscount_ / riskFreeDiscount_;

        // Should be handled more efficient! Each time the optimizer calls the operator
        // a new object is created. Need a Calc function which returns value for a given
        // spot.
        // Any better solution with current QuantLib architecture?

        boost::shared_ptr<BlackCalculator> blackCalc(
                   new BlackCalculator(payoff_, forwardPrice,
                                       standardDeviation_,riskFreeDiscount_));
        return blackCalc->value()-strike_;
    }

}