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

/*

 Copyright (C) 2013, 2015 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.

*/

#include <ql/pricingengines/swaption/basketgeneratingengine.hpp>

#include <ql/rebatedexercise.hpp>

#include <ql/math/optimization/levenbergmarquardt.hpp>

#include <ql/math/optimization/simplex.hpp>

#include <ql/models/shortrate/calibrationhelpers/swaptionhelper.hpp>

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

#include <ql/quotes/simplequote.hpp>

#include <cmath>

using std::exp;

using std::fabs;

namespace QuantLib {

    std::vector<ext::shared_ptr<BlackCalibrationHelper>>

    BasketGeneratingEngine::calibrationBasket(

        const ext::shared_ptr<Exercise>& exercise,

        const ext::shared_ptr<SwapIndex>& standardSwapBase,

        const ext::shared_ptr<SwaptionVolatilityStructure>& swaptionVolatility,

        const CalibrationBasketType basketType) const {

        QL_REQUIRE(

            !standardSwapBase->forwardingTermStructure().empty(),

            "standard swap base forwarding term structure must not be empty.");

        QL_REQUIRE(

            !standardSwapBase->exogenousDiscount() ||

                !standardSwapBase->discountingTermStructure().empty(),

            "standard swap base discounting term structure must not be empty.");

        std::vector<ext::shared_ptr<BlackCalibrationHelper> > result;

        Date today = Settings::instance().evaluationDate();

        Size minIdxAlive = static_cast<Size>(

            std::upper_bound(exercise->dates().begin(), exercise->dates().end(),

                             today) -

            exercise->dates().begin());

        ext::shared_ptr<RebatedExercise> rebEx =

            ext::dynamic_pointer_cast<RebatedExercise>(exercise);

        for (Size i = minIdxAlive; i < exercise->dates().size(); i++) {

            Date expiry = exercise->date(i);

            Real rebate = 0.0;

            Date rebateDate = expiry;

            if (rebEx != nullptr) {

                rebate = rebEx->rebate(i);

                rebateDate = rebEx->rebatePaymentDate(i);

            }

            ext::shared_ptr<SwaptionHelper> helper;

            switch (basketType) {

            case Naive: {

                Real swapLength = swaptionVolatility->dayCounter().yearFraction(

                    standardSwapBase->valueDate(expiry), underlyingLastDate());

                ext::shared_ptr<SmileSection> sec =

                    swaptionVolatility->smileSection(

                        expiry,

                        static_cast<Size>(std::lround(swapLength * 12.0)) * Months,

                        true);

                Real atmStrike = sec->atmLevel();

                Real atmVol;

                if (atmStrike == Null<Real>())

                    atmVol = sec->volatility(0.03);

                else

                    atmVol = sec->volatility(atmStrike);

                Real shift = sec->shift();

                helper = ext::make_shared<SwaptionHelper>(

                    expiry, underlyingLastDate(),

                    Handle<Quote>(ext::make_shared<SimpleQuote>(atmVol)),

                    standardSwapBase->iborIndex(),

                    standardSwapBase->fixedLegTenor(),

                    standardSwapBase->dayCounter(),

                    standardSwapBase->iborIndex()->dayCounter(),

                    standardSwapBase->exogenousDiscount()

                        ? standardSwapBase->discountingTermStructure()

                        : standardSwapBase->forwardingTermStructure(),

                    BlackCalibrationHelper::RelativePriceError, Null<Real>(), 1.0,

                    swaptionVolatility->volatilityType() ,shift);

                break;

            }

            case MaturityStrikeByDeltaGamma: {

                // determine the npv, first and second order derivatives at

                // $y=0$ of the underlying swap

                const Real h = 0.0001; // finite difference step in $y$, make

                                       // this a parameter of the engine ?

                Real zSpreadDsc =

                    oas_.empty() ? Real(1.0)

                                 : exp(-oas_->value() *

                                       onefactormodel_->termStructure()

                                           ->dayCounter()

                                           .yearFraction(expiry, rebateDate));

                Real npvm = underlyingNpv(expiry, -h) +

                            rebate *

                                onefactormodel_->zerobond(rebateDate, expiry,

                                                          -h, discountCurve_) *

                                zSpreadDsc;

                Real npv = underlyingNpv(expiry, 0.0) +

                           rebate * onefactormodel_->zerobond(

                                        rebateDate, expiry, 0, discountCurve_) *

                               zSpreadDsc;

                Real npvp = underlyingNpv(expiry, h) +

                            rebate *

                                onefactormodel_->zerobond(rebateDate, expiry, h,

                                                          discountCurve_) *

                                zSpreadDsc;

                Real delta = (npvp - npvm) / (2.0 * h);

                Real gamma = (npvp - 2.0 * npv + npvm) / (h * h);

                QL_REQUIRE(npv * npv + delta * delta + gamma * gamma > 0.0,

                           "(npv,delta,gamma) must have a positive norm");

                // debug output

                // std::cout << "EXOTIC npv " << npv << " delta " << delta

                //           << " gamma " << gamma << std::endl;

                // Real xtmp = -5.0;

                // std::cout

                //     << "********************************************EXERCISE "

                //     << expiry << " ******************" << std::endl;

                // std::cout << "globalExoticNpv;";

                // while (xtmp <= 5.0 + QL_EPSILON) {

                //     std::cout << underlyingNpv(expiry, xtmp) << ";";

                //     xtmp += 0.1;

                // }

                // std::cout << std::endl;

                // end debug output

                // play safe, we restrict the maximum maturity so to easily fit

                // in the date class restriction

                Real maxMaturity =

                    swaptionVolatility->dayCounter().yearFraction(

                        expiry, Date::maxDate() - 365);

                ext::shared_ptr<MatchHelper> matchHelper_;

                matchHelper_ = ext::make_shared<MatchHelper>(

                    underlyingType(), npv, delta, gamma, *onefactormodel_,

                    standardSwapBase, expiry, maxMaturity, h);

                // Optimize

                Array initial = initialGuess(expiry);

                QL_REQUIRE(initial.size() == 3,

                           "initial guess must have size 3 (but is "

                               << initial.size() << ")");

                EndCriteria ec(1000, 200, 1E-8, 1E-8, 1E-8); // make these

                                                             // criteria and the

                                                             // optimizer itself

                                                             // parameters of

                                                             // the method ?

                Constraint constraint = NoConstraint();

                Problem p(*matchHelper_, constraint, initial);

                LevenbergMarquardt lm;

                EndCriteria::Type ret = lm.minimize(p, ec);

                QL_REQUIRE(ret != EndCriteria::None &&

                               ret != EndCriteria::Unknown &&

                               ret != EndCriteria::MaxIterations,

                           "optimizer returns error (" << ret << ")");

                Array solution = p.currentValue();

                Real maturity = fabs(solution[1]);

                Size years = (Size)std::floor(maturity);

                maturity -= (Real)years;

                maturity *= 12.0;

                Size months = (Size)std::floor(maturity + 0.5);

                if (years == 0 && months == 0)

                    months = 1; // ensure a maturity of at least one months

                // maturity -= (Real)months; maturity *= 365.25;

                // Size days = (Size)std::floor(maturity);

                Period matPeriod =

                    years * Years + months * Months; //+days*Days;

                ext::shared_ptr<SmileSection> sec =

                    swaptionVolatility->smileSection(expiry, matPeriod, true);

                Real shift = sec->shift();

                // we have to floor the strike of the calibration instrument,

                // see warning in the header

                solution[2] = std::max(

                    solution[2], 0.00001 - shift); // floor at 0.1bp - shift

                // also the calibrated nominal may be zero, so we floor it, too

                solution[0] =

                    std::max(solution[0], 0.000001); // float at 0.01bp

                Real vol = sec->volatility(solution[2]);

                helper = ext::make_shared<SwaptionHelper>(

                    expiry, matPeriod,

                    Handle<Quote>(ext::make_shared<SimpleQuote>(

                                      vol)),

                    standardSwapBase->iborIndex(),

                    standardSwapBase->fixedLegTenor(),

                    standardSwapBase->dayCounter(),

                    standardSwapBase->iborIndex()->dayCounter(),

                    standardSwapBase->exogenousDiscount()

                        ? standardSwapBase->discountingTermStructure()

                        : standardSwapBase->forwardingTermStructure(),

                    BlackCalibrationHelper::RelativePriceError, solution[2],

                    fabs(solution[0]), swaptionVolatility->volatilityType(), shift);

                break;

            }

            default:

                QL_FAIL("Calibration basket type not known (" << basketType

                                                              << ")");

            }

            result.push_back(helper);

        }

        return result;

    }

}