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

/*
 Copyright (C) 2014 Master IMAFA - Polytech'Nice Sophia - Université de Nice Sophia Antipolis

 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/exercise.hpp>
#include <ql/pricingengines/exotic/analyticholderextensibleoptionengine.hpp>
#include <ql/math/distributions/bivariatenormaldistribution.hpp>
#include <utility>

using std::pow;
using std::log;
using std::exp;
using std::sqrt;

namespace QuantLib {

    AnalyticHolderExtensibleOptionEngine::AnalyticHolderExtensibleOptionEngine(
        ext::shared_ptr<GeneralizedBlackScholesProcess> process)
    : process_(std::move(process)) {
        registerWith(process_);
    }

    void AnalyticHolderExtensibleOptionEngine::calculate() const {
        //Spot
        Real S = process_->x0();
        Real r = riskFreeRate();
        Real b = r - dividendYield();
        Real X1 = strike();
        Real X2 = arguments_.secondStrike;
        Time T2 = secondExpiryTime();
        Time t1 = firstExpiryTime();
        Real A = arguments_.premium;


        Real z1 = this->z1();

        Real z2 = this->z2();

        Real rho = sqrt(t1 / T2);


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

        //QuantLib requires sigma * sqrt(T) rather than just sigma/volatility
        Real vol = volatility();

        //calculate dividend discount factor assuming continuous compounding (e^-rt)
        DiscountFactor growth = dividendDiscount(t1);
        //calculate payoff discount factor assuming continuous compounding
        DiscountFactor discount = riskFreeDiscount(t1);
        Real result = 0;
        constexpr double minusInf = -std::numeric_limits<double>::infinity();

        Real y1 = this->y1(payoff->optionType()),
             y2 = this->y2(payoff->optionType());
        if (payoff->optionType() == Option::Call) {
            //instantiate payoff function for a call
            ext::shared_ptr<PlainVanillaPayoff> vanillaCallPayoff =
                ext::make_shared<PlainVanillaPayoff>(Option::Call, X1);
            Real BSM = BlackScholesCalculator(vanillaCallPayoff, S, growth, vol*sqrt(t1), discount).value();
            result = BSM
                + S*exp((b - r)*T2)*M2(y1, y2, minusInf, z1, rho)
                - X2*exp(-r*T2)*M2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1), minusInf, z1 - vol*sqrt(T2), rho)
                - S*exp((b - r)*t1)*N2(y1, z2) + X1*exp(-r*t1)*N2(y1 - vol*sqrt(t1), z2 - vol*sqrt(t1))
                - A*exp(-r*t1)*N2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1));
        } else {
            //instantiate payoff function for a call
            ext::shared_ptr<PlainVanillaPayoff> vanillaPutPayoff =
                ext::make_shared<PlainVanillaPayoff>(Option::Put, X1);
            result = BlackScholesCalculator(vanillaPutPayoff, S, growth, vol*sqrt(t1), discount).value()
                - S*exp((b - r)*T2)*M2(y1, y2, minusInf, -z1, rho)
                + X2*exp(-r*T2)*M2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1), minusInf, -z1 + vol*sqrt(T2), rho)
                + S*exp((b - r)*t1)*N2(z2, y2) - X1*exp(-r*t1)*N2(z2 - vol*sqrt(t1), y2 - vol*sqrt(t1))
                - A*exp(-r*t1)*N2(y1 - vol*sqrt(t1), y2 - vol*sqrt(t1));
        }
        this->results_.value = result;
    }

    Real AnalyticHolderExtensibleOptionEngine::I1Call() const {
        Real Sv = process_->x0();
        Real A = arguments_.premium;

        if(A==0)
        {
            return 0;
        }
        else
        {
            BlackScholesCalculator bs = bsCalculator(Sv, Option::Call);
            Real ci = bs.value();
            Real dc = bs.delta();

            Real yi = ci - A;
            //da/ds = 0
            Real di = dc - 0;
            Real epsilon = 0.001;

            //Newton-Raphson process
            while (std::fabs(yi) > epsilon){
                Sv = Sv - yi / di;

                bs = bsCalculator(Sv, Option::Call);
                ci = bs.value();
                dc = bs.delta();

                yi = ci - A;
                di = dc - 0;
            }
            return Sv;
        }
    }

    Real AnalyticHolderExtensibleOptionEngine::I2Call() const {
        Real Sv = process_->x0();
        Real X1 = strike();
        Real X2 = arguments_.secondStrike;
        Real A = arguments_.premium;
        Time T2 = secondExpiryTime();
        Time t1 = firstExpiryTime();
        Real r=riskFreeRate();

        Real val=X1-X2*std::exp(-r*(T2-t1));
        if(A< val){
            return std::numeric_limits<Real>::infinity();
        } else {
            BlackScholesCalculator bs = bsCalculator(Sv, Option::Call);
            Real ci = bs.value();
            Real dc = bs.delta();

            Real yi = ci - A - Sv + X1;
            //da/ds = 1
            Real di = dc - 1;
            Real epsilon = 0.001;

            //Newton-Raphson process
            while (std::fabs(yi) > epsilon){
                Sv = Sv - yi / di;

                bs = bsCalculator(Sv, Option::Call);
                ci = bs.value();
                dc = bs.delta();

                yi = ci - A - Sv + X1;
                di = dc - 1;
            }
            return Sv;
        }
    }

    Real AnalyticHolderExtensibleOptionEngine::I1Put() const {
        Real Sv = process_->x0();
        //Srtike
        Real X1 = strike();
        //Premium
        Real A = arguments_.premium;

        BlackScholesCalculator bs = bsCalculator(Sv, Option::Put);
        Real pi = bs.value();
        Real dc = bs.delta();

        Real yi = pi - A + Sv - X1;
        //da/ds = 1
        Real di = dc - 1;
        Real epsilon = 0.001;

        //Newton-Raphson prosess
        while (std::fabs(yi) > epsilon){
            Sv = Sv - yi / di;

            bs = bsCalculator(Sv, Option::Put);
            pi = bs.value();
            dc = bs.delta();

            yi = pi - A + Sv - X1;
            di = dc - 1;
        }
        return Sv;
    }

    Real AnalyticHolderExtensibleOptionEngine::I2Put() const {
        Real Sv = process_->x0();
        Real A = arguments_.premium;
        if(A==0){
            return std::numeric_limits<Real>::infinity();
        }
        else{
            BlackScholesCalculator bs = bsCalculator(Sv, Option::Put);
            Real pi = bs.value();
            Real dc = bs.delta();

            Real yi = pi - A;
            //da/ds = 0
            Real di = dc - 0;
            Real epsilon = 0.001;

            //Newton-Raphson prosess
            while (std::fabs(yi) > epsilon){
                Sv = Sv - yi / di;

                bs = bsCalculator(Sv, Option::Put);
                pi = bs.value();
                dc = bs.delta();

                yi = pi - A;
                di = dc - 0;
            }
            return Sv;
        }
    }


    BlackScholesCalculator AnalyticHolderExtensibleOptionEngine::bsCalculator(
                                    Real spot, Option::Type optionType) const {
        //Real spot = process_->x0();
        Real vol;
        DiscountFactor growth;
        DiscountFactor discount;
        Real X2 = arguments_.secondStrike;
        Time T2 = secondExpiryTime();
        Time t1 = firstExpiryTime();
        Time t = T2 - t1;

        //payoff
        ext::shared_ptr<PlainVanillaPayoff > vanillaPayoff =
            ext::make_shared<PlainVanillaPayoff>(optionType, X2);

        //QuantLib requires sigma * sqrt(T) rather than just sigma/volatility
        vol = volatility() * std::sqrt(t);
        //calculate dividend discount factor assuming continuous compounding (e^-rt)
        growth = dividendDiscount(t);
        //calculate payoff discount factor assuming continuous compounding
        discount = riskFreeDiscount(t);

        BlackScholesCalculator bs(vanillaPayoff, spot, growth, vol, discount);
        return bs;
    }

    Real AnalyticHolderExtensibleOptionEngine::M2(Real a, Real b, Real c, Real d, Real rho) const {
        BivariateCumulativeNormalDistributionDr78 CmlNormDist(rho);
        return CmlNormDist(b, d) - CmlNormDist(a, d) - CmlNormDist(b, c) + CmlNormDist(a,c);
    }

    Real AnalyticHolderExtensibleOptionEngine::N2(Real a, Real b) const {
        CumulativeNormalDistribution  NormDist;
        return NormDist(b) - NormDist(a);
    }

    Real AnalyticHolderExtensibleOptionEngine::strike() const {
        ext::shared_ptr<PlainVanillaPayoff> payoff =
            ext::dynamic_pointer_cast<PlainVanillaPayoff>(arguments_.payoff);
        QL_REQUIRE(payoff, "non-plain payoff given");
        return payoff->strike();
    }

    Time AnalyticHolderExtensibleOptionEngine::firstExpiryTime() const {
        return process_->time(arguments_.exercise->lastDate());
    }

    Time AnalyticHolderExtensibleOptionEngine::secondExpiryTime() const {
        return process_->time(arguments_.secondExpiryDate);
    }

    Volatility AnalyticHolderExtensibleOptionEngine::volatility() const {
        return process_->blackVolatility()->blackVol(firstExpiryTime(), strike());
    }
    Rate AnalyticHolderExtensibleOptionEngine::riskFreeRate() const {
        return process_->riskFreeRate()->zeroRate(firstExpiryTime(), Continuous,
            NoFrequency);
    }
    Rate AnalyticHolderExtensibleOptionEngine::dividendYield() const {
        return process_->dividendYield()->zeroRate(firstExpiryTime(),
            Continuous, NoFrequency);
    }

    DiscountFactor AnalyticHolderExtensibleOptionEngine::dividendDiscount(Time t) const {
        return process_->dividendYield()->discount(t);
    }

    DiscountFactor AnalyticHolderExtensibleOptionEngine::riskFreeDiscount(Time t) const {
        return process_->riskFreeRate()->discount(t);
    }

    Real AnalyticHolderExtensibleOptionEngine::y1(Option::Type type) const {
        Real S = process_->x0();
        Real I2 = (type == Option::Call) ? I2Call() : I2Put();

        Real b = riskFreeRate() - dividendYield();
        Real vol = volatility();
        Time t1 = firstExpiryTime();

        return (log(S / I2) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1));
    }

    Real AnalyticHolderExtensibleOptionEngine::y2(Option::Type type) const {
        Real S = process_->x0();
        Real I1 = (type == Option::Call) ? I1Call() : I1Put();

        Real b = riskFreeRate() - dividendYield();
        Real vol = volatility();
        Time t1 = firstExpiryTime();

        return (log(S / I1) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1));
    }

    Real AnalyticHolderExtensibleOptionEngine::z1() const {
        Real S = process_->x0();
        Real X2 = arguments_.secondStrike;
        Real b = riskFreeRate() - dividendYield();
        Real vol = volatility();
        Time T2 = secondExpiryTime();

        return (log(S / X2) + (b + pow(vol, 2) / 2)*T2) / (vol*sqrt(T2));
    }

    Real AnalyticHolderExtensibleOptionEngine::z2() const {
        Real S = process_->x0();
        Real X1 = strike();

        Real b = riskFreeRate() - dividendYield();
        Real vol = volatility();
        Time t1 = firstExpiryTime();

        return (log(S / X1) + (b + pow(vol, 2) / 2)*t1) / (vol*sqrt(t1));
    }

}