markovfunctional.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2013, 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
 <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/math/integrals/gaussianquadratures.hpp>
#include <ql/math/interpolations/cubicinterpolation.hpp>
#include <ql/math/solvers1d/brent.hpp>
#include <ql/models/shortrate/onefactormodels/markovfunctional.hpp>
#include <ql/termstructures/volatility/atmadjustedsmilesection.hpp>
#include <ql/termstructures/volatility/atmsmilesection.hpp>
#include <ql/termstructures/volatility/kahalesmilesection.hpp>
#include <ql/termstructures/volatility/sabrinterpolatedsmilesection.hpp>
#include <ql/termstructures/volatility/smilesection.hpp>
#include <ql/termstructures/volatility/smilesectionutils.hpp>
#include <utility>
 
namespace QuantLib {
 
    MarkovFunctional::MarkovFunctional(const Handle<YieldTermStructure>& termStructure,
                                       const Real reversion,
                                       std::vector<Date> volstepdates,
                                       std::vector<Real> volatilities,
                                       const Handle<SwaptionVolatilityStructure>& swaptionVol,
                                       const std::vector<Date>& swaptionExpiries,
                                       const std::vector<Period>& swaptionTenors,
                                       const ext::shared_ptr<SwapIndex>& swapIndexBase,
                                       MarkovFunctional::ModelSettings modelSettings)
    : Gaussian1dModel(termStructure), CalibratedModel(1), modelSettings_(std::move(modelSettings)),
      capletCalibrated_(false), reversion_(ConstantParameter(reversion, NoConstraint())),
      sigma_(arguments_[0]), volstepdates_(std::move(volstepdates)),
      volatilities_(std::move(volatilities)), swaptionVol_(swaptionVol),
      swaptionExpiries_(swaptionExpiries), swaptionTenors_(swaptionTenors),
      swapIndexBase_(swapIndexBase), iborIndex_(swapIndexBase->iborIndex()) {
 
        QL_REQUIRE(swaptionExpiries.size() == swaptionTenors.size(),
                   "number of swaption expiries ("
                       << swaptionExpiries.size()
                       << ") is differnt from number of swaption tenors ("
                       << swaptionTenors.size() << ")");
        QL_REQUIRE(!swaptionExpiries.empty(),
                   "need at least one swaption expiry to calibrate numeraire");
        QL_REQUIRE(!termStructure.empty(),
                   "yield term structure handle is empty");
        QL_REQUIRE(!swaptionVol.empty(),
                   "swaption volatility structure is empty");
        modelSettings_.validate();
        initialize();
    }
 
    MarkovFunctional::MarkovFunctional(const Handle<YieldTermStructure>& termStructure,
                                       const Real reversion,
                                       std::vector<Date> volstepdates,
                                       std::vector<Real> volatilities,
                                       const Handle<OptionletVolatilityStructure>& capletVol,
                                       const std::vector<Date>& capletExpiries,
                                       ext::shared_ptr<IborIndex> iborIndex,
                                       MarkovFunctional::ModelSettings modelSettings)
    : Gaussian1dModel(termStructure), CalibratedModel(1), modelSettings_(std::move(modelSettings)),
      capletCalibrated_(true), reversion_(ConstantParameter(reversion, NoConstraint())),
      sigma_(arguments_[0]), volstepdates_(std::move(volstepdates)),
      volatilities_(std::move(volatilities)), capletVol_(capletVol),
      capletExpiries_(capletExpiries), iborIndex_(std::move(iborIndex)) {
 
        QL_REQUIRE(!capletExpiries.empty(),
                   "need at least one caplet expiry to calibrate numeraire");
        QL_REQUIRE(!termStructure.empty(),
                   "yield term structure handle is empty");
        QL_REQUIRE(!capletVol.empty(), "caplet volatility structure is empty");
        modelSettings_.validate();
        initialize();
    }
 
    void MarkovFunctional::updateTimes() const {
        QL_MFMESSAGE(modelOutputs_,"updating times");
        updateTimes1();
        updateTimes2();
    }
 
    void MarkovFunctional::updateTimes1() const {
        volsteptimes_.clear();
        int j = 0;
        for (auto i = volstepdates_.begin(); i != volstepdates_.end(); ++i, ++j) {
            volsteptimes_.push_back(termStructure()->timeFromReference(*i));
            volsteptimesArray_[j] = volsteptimes_[j];
            if (j == 0)
                QL_REQUIRE(volsteptimes_[0] > 0.0,
                           "volsteptimes must be positive (" << volsteptimes_[0]
                                                             << ")");
            else
                QL_REQUIRE(volsteptimes_[j] > volsteptimes_[j - 1],
                           "volsteptimes must be strictly increasing ("
                               << volsteptimes_[j - 1] << "@" << (j - 1) << ", "
                               << volsteptimes_[j] << "@" << j << ")");
        }
    }
 
    void MarkovFunctional::updateTimes2() const {
        numeraireTime_ = termStructure()->timeFromReference(numeraireDate_);
        times_.clear();
        times_.push_back(0.0);
        modelOutputs_.expiries_.clear();
        modelOutputs_.tenors_.clear();
        for (auto& calibrationPoint : calibrationPoints_) {
            times_.push_back(termStructure()->timeFromReference(calibrationPoint.first));
            modelOutputs_.expiries_.push_back(calibrationPoint.first);
            modelOutputs_.tenors_.push_back(calibrationPoint.second.tenor_);
        }
        times_.push_back(numeraireTime_);
        QL_REQUIRE(volatilities_.size() == volsteptimes_.size() + 1,
                   "there must be n+1 volatilities ("
                       << volatilities_.size()
                       << ") for n volatility step times ("
                       << volsteptimes_.size() << ")");
    }
 
 
    void MarkovFunctional::initialize() {
 
        QL_MFMESSAGE(modelOutputs_, "initializing");
        modelOutputs_.dirty_ = true;
 
        modelOutputs_.settings_ = modelSettings_;
 
        GaussHermiteIntegration gaussHermite(
            modelSettings_.gaussHermitePoints_);
        normalIntegralX_ = gaussHermite.x();
        normalIntegralW_ = gaussHermite.weights();
        for (Size i = 0; i < normalIntegralX_.size(); i++) {
            normalIntegralW_[i] *=
                std::exp(-normalIntegralX_[i] * normalIntegralX_[i]) * M_1_SQRTPI;
            normalIntegralX_[i] *= M_SQRT2;
        }
 
        volsteptimesArray_ = Array(volstepdates_.size());
 
        updateTimes1();
 
        if (capletCalibrated_) {
            for (auto capletExpirie : capletExpiries_) {
                makeCapletCalibrationPoint(capletExpirie);
            }
        } else {
            std::vector<Date>::const_iterator i;
            std::vector<Period>::const_iterator j;
            for (i = swaptionExpiries_.begin(), j = swaptionTenors_.begin();
                 i != swaptionExpiries_.end(); ++i, ++j) {
                makeSwaptionCalibrationPoint(*i, *j);
            }
        }
 
        bool done;
        numeraireDate_ = Date::minDate();
        do {
            Date numeraireKnown = numeraireDate_;
            done = true;
            for (auto i = calibrationPoints_.rbegin(); i != calibrationPoints_.rend() && done;
                 ++i) {
                if (i->second.paymentDates_.back() > numeraireDate_) {
                    numeraireDate_ = i->second.paymentDates_.back();
                    numeraireKnown = i->second.paymentDates_.back();
                    done = i == calibrationPoints_.rbegin();
                }
                for (auto j = i->second.paymentDates_.rbegin();
                     j != i->second.paymentDates_.rend() && done; ++j) {
                    if (*j < numeraireKnown) {
                        if (capletCalibrated_) {
                            makeCapletCalibrationPoint(*j);
                            done = false;
                            break;
                        } else {
                            Size months = std::max(
                                1, static_cast<Integer>(((numeraireKnown - *j) / 365.25) * 12.0));
                            while (underlyingSwap(swapIndexBase_, *j, months * Months)
                                       ->maturityDate() < numeraireKnown)
                                ++months;
                            makeSwaptionCalibrationPoint(*j, months * Months);
                            done = false;
                            break;
                        }
                    }
                }
                if (done) {
                    numeraireKnown = i->first;
                }
            }
        } while (!done);
 
        updateTimes2();
 
        sigma_ =
            PiecewiseConstantParameter(volsteptimes_, PositiveConstraint());
        for (Size i = 0; i < sigma_.size(); i++) {
            sigma_.setParam(i, volatilities_[i]);
        }
 
        stateProcess_ = ext::make_shared<MfStateProcess>(
            reversion_(0.0), volsteptimesArray_, sigma_.params());
 
        y_ = yGrid(modelSettings_.yStdDevs_, modelSettings_.yGridPoints_);
 
        discreteNumeraire_ = ext::make_shared<Matrix>(
            times_.size(), 2 * modelSettings_.yGridPoints_ + 1, 1.0);
        for (Size i = 0; i < times_.size(); i++) {
            ext::shared_ptr<Interpolation> numInt(new CubicInterpolation(
                y_.begin(), y_.end(), discreteNumeraire_->row_begin(i),
                CubicInterpolation::Spline, true, CubicInterpolation::Lagrange,
                0.0, CubicInterpolation::Lagrange, 0.0));
            numInt->enableExtrapolation();
            numeraire_.push_back(numInt);
        }
 
        registerWith(termStructure());
        if (!swaptionVol_.empty())
            registerWith(swaptionVol_);
        if (!capletVol_.empty())
            registerWith(capletVol_);
    }
 
    void MarkovFunctional::makeSwaptionCalibrationPoint(const Date &expiry,
                                                        const Period &tenor) {
 
        QL_REQUIRE(calibrationPoints_.count(expiry) == 0,
                   "swaption expiry ("
                       << expiry
                       << ") occurs more than once in calibration set");
 
        CalibrationPoint p;
        p.isCaplet_ = false;
        p.tenor_ = tenor;
 
        ext::shared_ptr<VanillaSwap> underlying = underlyingSwap(swapIndexBase_, expiry, tenor);
 
        Schedule sched = underlying->fixedSchedule();
        const Calendar& cal = sched.calendar();
        BusinessDayConvention bdc = underlying->paymentConvention();
 
        for (unsigned int k = 1; k < sched.size(); k++) {
            p.yearFractions_.push_back(
                swapIndexBase_->dayCounter().yearFraction(
                    k == 1 ? expiry : sched.date(k - 1), sched.date(k)));
            p.paymentDates_.push_back(cal.adjust(sched.date(k), bdc));
        }
        calibrationPoints_[expiry] = p;
    }
 
    void MarkovFunctional::makeCapletCalibrationPoint(const Date &expiry) {
 
        QL_REQUIRE(
            calibrationPoints_.count(expiry) == 0,
            "caplet expiry (" << expiry
                              << ") occurs more than once in calibration set");
 
        CalibrationPoint p;
        p.isCaplet_ = true;
        // p.expiry_ = expiry;
        p.tenor_ = iborIndex_->tenor();
        Date valueDate = iborIndex_->valueDate(expiry);
        Date endDate = iborIndex_->fixingCalendar().advance(
            valueDate, iborIndex_->tenor(), iborIndex_->businessDayConvention(),
            iborIndex_->endOfMonth());
        // FIXME Here we should use a calculation date calendar ?
        p.paymentDates_.push_back(endDate);
        p.yearFractions_.push_back(
            iborIndex_->dayCounter().yearFraction(expiry, endDate));
        // adjust the first period to start on expiry
        calibrationPoints_[expiry] = p;
    }
 
    void MarkovFunctional::updateSmiles() const {
 
        QL_MFMESSAGE(modelOutputs_, "updating smiles");
        modelOutputs_.dirty_ = true;
 
        arbitrageIndices_.clear();
 
        Size pointIndex = 0;
 
        for (auto i = calibrationPoints_.rbegin(); i != calibrationPoints_.rend(); ++i) {
 
            ext::shared_ptr<SmileSection> smileSection;
            if (i->second.isCaplet_) {
                i->second.annuity_ =
                    i->second.yearFractions_[0] *
                    termStructure()->discount(i->second.paymentDates_[0], true);
                i->second.atm_ = (termStructure()->discount(i->first, true) -
                                  termStructure()->discount(
                                      i->second.paymentDates_[0], true)) /
                                 i->second.annuity_;
                smileSection = capletVol_->smileSection(i->first, true);
            } else {
                Real annuity = 0.0;
                for (unsigned int k = 0; k < i->second.paymentDates_.size();
                     k++) {
                    annuity += i->second.yearFractions_[k] *
                               termStructure()->discount(
                                   i->second.paymentDates_[k], true);
                }
                i->second.annuity_ = annuity;
                i->second.atm_ = (termStructure()->discount(i->first, true) -
                                  termStructure()->discount(
                                      i->second.paymentDates_.back(), true)) /
                                 annuity;
                smileSection = swaptionVol_->smileSection(
                    i->first, i->second.tenor_, true);
            }
 
            i->second.rawSmileSection_ = ext::shared_ptr<SmileSection>(
                new AtmSmileSection(smileSection, i->second.atm_));
 
            int forcedLeftIndex = -1;
            int forcedRightIndex = QL_MAX_INTEGER;
            if(forcedArbitrageIndices_.size() > pointIndex) {
                forcedLeftIndex = forcedArbitrageIndices_[pointIndex].first;
                forcedRightIndex = forcedArbitrageIndices_[pointIndex].second;
            }
 
            if ((modelSettings_.adjustments_ & ModelSettings::KahaleSmile) != 0) {
 
                i->second.smileSection_ = ext::make_shared<KahaleSmileSection>(
                    
                        i->second.rawSmileSection_, i->second.atm_,
                        (modelSettings_.adjustments_ &
                         ModelSettings::KahaleInterpolation) != 0,
                        (modelSettings_.adjustments_ &
                         ModelSettings::SmileExponentialExtrapolation) != 0,
                        (modelSettings_.adjustments_ &
                         ModelSettings::SmileDeleteArbitragePoints) != 0,
                        modelSettings_.smileMoneynessCheckpoints_,
                        modelSettings_.digitalGap_,
                        forcedLeftIndex, forcedRightIndex);
 
                arbitrageIndices_.push_back(
                    ext::dynamic_pointer_cast<KahaleSmileSection>(
                        i->second.smileSection_)->coreIndices());
 
            } else {
 
                if ((modelSettings_.adjustments_ & ModelSettings::SabrSmile) != 0) {
 
                    SmileSectionUtils ssutils(
                        *i->second.rawSmileSection_,
                        modelSettings_.smileMoneynessCheckpoints_);
                    std::vector<Real> k = ssutils.strikeGrid();
                    k.erase(k.begin()); // the first strike is zero which we do
                                        // not want in the sabr calibration
                    QL_REQUIRE(i->second.rawSmileSection_->volatilityType() ==
                                   ShiftedLognormal,
                               "MarkovFunctional: SABR calibration to normal "
                               "input volatilities is not supported");
                    QL_REQUIRE(
                        k.size() >= 4,
                        "for sabr calibration at least 4 points are needed (is "
                            << k.size() << ")");
                    std::vector<Real> v;
                    v.reserve(k.size());
                    for (Real j : k) {
                        v.push_back(i->second.rawSmileSection_->volatility(j));
                    }
 
                    // TODO should we fix beta to avoid numerical instabilities
                    // during calibration ?
                    ext::shared_ptr<SabrInterpolatedSmileSection> sabrSection(
                        new SabrInterpolatedSmileSection(
                            i->first, i->second.atm_, k, false,
                            i->second.rawSmileSection_->volatility(
                                i->second.atm_),
                            v, 0.03, 0.80, 0.50, 0.00, false, false, false,
                            false, true, ext::shared_ptr<EndCriteria>(),
                            ext::shared_ptr<OptimizationMethod>(),
                            Actual365Fixed(),
                                i->second.rawSmileSection_->shift()));
 
                    // we make the sabr section arbitrage free by superimposing
                    // a kahalesection
 
                    i->second.smileSection_ = ext::make_shared<
                        KahaleSmileSection>(
                        sabrSection, i->second.atm_, false,
                        (modelSettings_.adjustments_ &
                         ModelSettings::SmileExponentialExtrapolation) != 0,
                        (modelSettings_.adjustments_ &
                         ModelSettings::SmileDeleteArbitragePoints) != 0,
                        modelSettings_.smileMoneynessCheckpoints_,
                        modelSettings_.digitalGap_,
                        forcedLeftIndex, forcedRightIndex);
 
                    arbitrageIndices_.push_back(
                        ext::dynamic_pointer_cast<KahaleSmileSection>(
                            i->second.smileSection_)->coreIndices());
 
                } else if ((modelSettings_.adjustments_ & ModelSettings::CustomSmile) != 0) {
 
                    // Custom smile section is af by assumption
                    i->second.smileSection_ =
                        modelSettings_.customSmileFactory_->smileSection(
                            i->second.rawSmileSection_, i->second.atm_);
                    arbitrageIndices_.emplace_back(Null<Size>(), Null<Size>());
                } else { // no smile pretreatment
 
                    i->second.smileSection_ = i->second.rawSmileSection_;
                }
            }
 
            // custom smile will take care of this itself
            if ((modelSettings_.adjustments_ & ModelSettings::CustomSmile) == 0) {
                i->second.minRateDigital_ =
                    i->second.smileSection_->digitalOptionPrice(
                        modelSettings_.lowerRateBound_ -
                            i->second.smileSection_->shift(),
                        Option::Call, i->second.annuity_,
                        modelSettings_.digitalGap_);
                i->second.maxRateDigital_ =
                    i->second.smileSection_->digitalOptionPrice(
                        modelSettings_.upperRateBound_ -
                            i->second.smileSection_->shift(),
                        Option::Call, i->second.annuity_,
                        modelSettings_.digitalGap_);
            }
 
            ++pointIndex;
        }
    }
 
    void MarkovFunctional::updateNumeraireTabulation() const {
 
        QL_MFMESSAGE(modelOutputs_, "updating numeraire tabulation");
        modelOutputs_.dirty_ = true;
 
        modelOutputs_.adjustmentFactors_.clear();
        modelOutputs_.digitalsAdjustmentFactors_.clear();
 
        int idx = times_.size() - 2;
 
        for (auto i = calibrationPoints_.rbegin(); i != calibrationPoints_.rend(); ++i, --idx) {
 
            ext::shared_ptr<CustomSmileSection> mfSec;
            if ((modelSettings_.adjustments_ & ModelSettings::CustomSmile) != 0) {
                mfSec = ext::dynamic_pointer_cast<CustomSmileSection>(
                    i->second.smileSection_);
                QL_REQUIRE(mfSec,
                           "no CustomSmileSection given, this is unexpected...");
            }
 
            Array discreteDeflatedAnnuities(y_.size(), 0.0);
            Array deflatedFinalPayments;
 
            Real numeraire0 = termStructure()->discount(numeraireTime_, true);
            Real normalization =
                termStructure()->discount(times_[idx], true) / numeraire0;
 
            for (unsigned int k = 0; k < i->second.paymentDates_.size(); k++) {
                deflatedFinalPayments =
                    deflatedZerobondArray(termStructure()->timeFromReference(
                                              i->second.paymentDates_[k]),
                                          times_[idx], y_);
                discreteDeflatedAnnuities +=
                    deflatedFinalPayments * i->second.yearFractions_[k];
            }
 
            CubicInterpolation deflatedAnnuities(
                y_.begin(), y_.end(), discreteDeflatedAnnuities.begin(),
                CubicInterpolation::Spline, true, CubicInterpolation::Lagrange,
                0.0, CubicInterpolation::Lagrange, 0.0);
            deflatedAnnuities.enableExtrapolation();
 
            Real digitalsCorrectionFactor = 1.0;
            modelOutputs_.digitalsAdjustmentFactors_.insert(
                modelOutputs_.digitalsAdjustmentFactors_.begin(),
                digitalsCorrectionFactor);
 
            Real digital = 0.0, swapRate, swapRate0;
 
            for (int c = 0;
                 c == 0 ||
                 (c == 1 && ((modelSettings_.adjustments_ & ModelSettings::AdjustDigitals) != 0));
                 c++) {
 
                if (c == 1) {
                    digitalsCorrectionFactor = i->second.annuity_ / digital;
                    modelOutputs_.digitalsAdjustmentFactors_.front() =
                        digitalsCorrectionFactor;
                }
 
                digital = 0.0;
                swapRate0 =
                    modelSettings_.upperRateBound_ / 2.0; // initial guess
                for (int j = y_.size() - 1; j >= 0; j--) {
 
                    Real integral = 0.0;
 
                    if (j == (int)(y_.size() - 1)) {
                        if ((modelSettings_.adjustments_ &
                             ModelSettings::NoPayoffExtrapolation) == 0) {
                            if ((modelSettings_.adjustments_ &
                                 ModelSettings::ExtrapolatePayoffFlat) != 0) {
                                integral = gaussianShiftedPolynomialIntegral(
                                    0.0, 0.0, 0.0, 0.0,
                                    discreteDeflatedAnnuities[j - 1], y_[j - 1],
                                    y_[j], 100.0);
                            } else {
                                Real ca =
                                    deflatedAnnuities.aCoefficients()[j - 1];
                                Real cb =
                                    deflatedAnnuities.bCoefficients()[j - 1];
                                Real cc =
                                    deflatedAnnuities.cCoefficients()[j - 1];
                                integral = gaussianShiftedPolynomialIntegral(
                                    0.0, cc, cb, ca,
                                    discreteDeflatedAnnuities[j - 1], y_[j - 1],
                                    y_[j], 100.0);
                            }
                        }
                    } else {
                        Real ca = deflatedAnnuities.aCoefficients()[j];
                        Real cb = deflatedAnnuities.bCoefficients()[j];
                        Real cc = deflatedAnnuities.cCoefficients()[j];
                        integral = gaussianShiftedPolynomialIntegral(
                            0.0, cc, cb, ca, discreteDeflatedAnnuities[j],
                            y_[j], y_[j], y_[j + 1]);
                    }
 
                    if (integral < 0) {
                        QL_MFMESSAGE(modelOutputs_,
                                     "WARNING: integral for digitalPrice is "
                                     "negative for j="
                                         << j << " (" << integral
                                         << ") --- reset it to zero.");
                        integral = 0.0;
                    }
 
                    digital += integral * numeraire0 * digitalsCorrectionFactor;
 
                    bool check = true;
                    if ((modelSettings_.adjustments_ & ModelSettings::CustomSmile) != 0) {
                        swapRate = mfSec->inverseDigitalCall(
                            digital, i->second.annuity_);
                    } else if (digital >= i->second.minRateDigital_) {
                        swapRate = modelSettings_.lowerRateBound_ -
                                   i->second.rawSmileSection_->shift();
                        check = false;
                    } else if (digital <= i->second.maxRateDigital_) {
                        swapRate = modelSettings_.upperRateBound_;
                        check = false;
                    } else {
                        swapRate = marketSwapRate(
                            i->first, i->second, digital, swapRate0,
                            i->second.rawSmileSection_->shift());
                    }
                    if (check && j < (int)y_.size() - 1 &&
                        swapRate > swapRate0) {
                        QL_MFMESSAGE(
                            modelOutputs_,
                            "WARNING: swap rate is decreasing in y for "
                            "t=" << times_[idx]
                                 << ", j=" << j << " (y, swap rate) is ("
                                 << y_[j] << "," << swapRate << ") but for j="
                                 << j + 1 << " it is (" << y_[j + 1] << ","
                                 << swapRate0 << ") --- reset rate to "
                                 << swapRate0 << " in node j=" << j);
                        swapRate = swapRate0;
                    }
                    swapRate0 = swapRate;
                    Real numeraire =
                        1.0 / std::max(swapRate * discreteDeflatedAnnuities[j] +
                                       deflatedFinalPayments[j], 1E-6);
                    (*discreteNumeraire_)[idx][j] = numeraire * normalization;
                }
            }
 
            if ((modelSettings_.adjustments_ & ModelSettings::AdjustYts) != 0) {
                numeraire_[idx]->update();
                Real modelDeflatedZerobond = deflatedZerobond(times_[idx], 0.0);
                Real marketDeflatedZerobond =
                    termStructure()->discount(times_[idx], true) /
                    termStructure()->discount(numeraireTime_, true);
                for (int j = y_.size() - 1; j >= 0; j--) {
                    (*discreteNumeraire_)[idx][j] *=
                        modelDeflatedZerobond / marketDeflatedZerobond;
                }
                modelOutputs_.adjustmentFactors_.insert(
                    modelOutputs_.adjustmentFactors_.begin(),
                    modelDeflatedZerobond / marketDeflatedZerobond);
            } else {
                modelOutputs_.adjustmentFactors_.insert(
                    modelOutputs_.adjustmentFactors_.begin(), 1.0);
            }
 
            numeraire_[idx]->update();
        }
    }
 
    const MarkovFunctional::ModelOutputs &
    MarkovFunctional::modelOutputs() const {
 
        if (modelOutputs_.dirty_) {
 
            calculate();
 
            // yield term structure
            modelOutputs_.marketZerorate_.clear();
            modelOutputs_.modelZerorate_.clear();
            for (Size i = 1; i < times_.size() - 1; i++) {
                modelOutputs_.marketZerorate_.push_back(
                    termStructure()->zeroRate(times_[i], QuantLib::Continuous,
                                              QuantLib::Annual));
                // we need to put a small positive time here since the zerobond
                // implementation optimizes the case t=0.0 then using the
                // initial yts
                modelOutputs_.modelZerorate_.push_back(
                    -std::log(zerobond(times_[i], 1.0E-10)) / times_[i]);
            }
 
            // volatility surface
            modelOutputs_.smileStrikes_.clear();
            modelOutputs_.marketCallPremium_.clear();
            modelOutputs_.marketPutPremium_.clear();
            modelOutputs_.modelCallPremium_.clear();
            modelOutputs_.modelPutPremium_.clear();
            modelOutputs_.marketVega_.clear();
            modelOutputs_.marketRawCallPremium_.clear();
            modelOutputs_.marketRawPutPremium_.clear();
 
            for (auto& calibrationPoint : calibrationPoints_) {
                modelOutputs_.atm_.push_back(calibrationPoint.second.atm_);
                modelOutputs_.annuity_.push_back(calibrationPoint.second.annuity_);
                ext::shared_ptr<SmileSection> sec = calibrationPoint.second.smileSection_;
                ext::shared_ptr<SmileSection> rawSec = calibrationPoint.second.rawSmileSection_;
                SmileSectionUtils ssutils(*sec, modelSettings_.smileMoneynessCheckpoints_,
                                          calibrationPoint.second.atm_);
                Real shift = sec->shift();
                std::vector<Real> money = ssutils.moneyGrid();
                std::vector<Real> strikes, marketCall, marketPut, modelCall,
                    modelPut, marketVega, marketRawCall, marketRawPut;
                for (Size j = 0; j < money.size(); j++) {
                    strikes.push_back(sec->volatilityType() == Normal ?
                                          Real(calibrationPoint.second.atm_ + money[j]) :
                                          money[j] * (calibrationPoint.second.atm_ + shift) -
                                              shift);
                    try {
                        marketRawCall.push_back(rawSec->optionPrice(
                            strikes[j], Option::Call, calibrationPoint.second.annuity_));
                        marketRawPut.push_back(rawSec->optionPrice(
                            strikes[j], Option::Put, calibrationPoint.second.annuity_));
                    }
                    catch (Error&) {
                        // the smile section might not be able to output an
                        // option price because it has no atm level
                        marketRawCall.push_back(0.0);
                        marketRawPut.push_back(0.0);
                    }
                    marketCall.push_back(sec->optionPrice(strikes[j], Option::Call,
                                                          calibrationPoint.second.annuity_));
                    marketPut.push_back(sec->optionPrice(strikes[j], Option::Put,
                                                         calibrationPoint.second.annuity_));
                    modelCall.push_back(
                        calibrationPoint.second.isCaplet_ ?
                            capletPriceInternal(Option::Call, calibrationPoint.first, strikes[j],
                                                Null<Date>(), 0.0, true) :
                            swaptionPriceInternal(Option::Call, calibrationPoint.first,
                                                  calibrationPoint.second.tenor_, strikes[j],
                                                  Null<Date>(), 0.0, true));
                    modelPut.push_back(
                        calibrationPoint.second.isCaplet_ ?
                            capletPriceInternal(Option::Put, calibrationPoint.first, strikes[j],
                                                Null<Date>(), 0.0, true) :
                            swaptionPriceInternal(Option::Put, calibrationPoint.first,
                                                  calibrationPoint.second.tenor_, strikes[j],
                                                  Null<Date>(), 0.0, true));
                    marketVega.push_back(sec->vega(strikes[j], calibrationPoint.second.annuity_));
                }
                modelOutputs_.smileStrikes_.push_back(strikes);
                modelOutputs_.marketCallPremium_.push_back(marketCall);
                modelOutputs_.marketPutPremium_.push_back(marketPut);
                modelOutputs_.modelCallPremium_.push_back(modelCall);
                modelOutputs_.modelPutPremium_.push_back(modelPut);
                modelOutputs_.marketVega_.push_back(marketVega);
                modelOutputs_.marketRawCallPremium_.push_back(marketRawCall);
                modelOutputs_.marketRawPutPremium_.push_back(marketRawPut);
            }
 
            modelOutputs_.dirty_ = false;
        }
 
        return modelOutputs_;
    }
 
    Array MarkovFunctional::numeraireArray(const Time t, const Array& y) const {
 
        calculate();
        Array res(y.size(), termStructure()->discount(numeraireTime_, true));
        if (t < QL_EPSILON)
            return res;
 
        Real inverseNormalization =
            termStructure()->discount(numeraireTime_, true) /
            termStructure()->discount(t, true);
 
        Time tz = std::min(t, times_.back());
        Size i = std::min<Size>(
            std::upper_bound(times_.begin(), times_.end() - 1, t) -
                times_.begin(),
            times_.size() - 1);
 
        Real ta = times_[i - 1];
        Real tb = times_[i];
        Real dt = tb - ta;
 
        for (Size j = 0; j < y.size(); j++) {
            Real yv = y[j];
            if (yv < y_.front())
                yv = y_.front();
            // FIXME flat extrapolation should be incoperated into interpolation
            // object, see above
            if (yv > y_.back())
                yv = y_.back();
            Real na = (*numeraire_[i - 1])(yv);
            Real nb = (*numeraire_[i])(yv);
            res[j] =
                inverseNormalization / ((tz - ta) / nb + (tb - tz) / na) * dt;
            // linear in reciprocal of normalized numeraire
        }
 
        return res;
    }
 
    Array MarkovFunctional::zerobondArray(const Time T, const Time t, const Array& y) const {
 
        return deflatedZerobondArray(T, t, y) * numeraireArray(t, y);
    }
 
    Array MarkovFunctional::deflatedZerobondArray(const Time T, const Time t, const Array& y) const {
 
        calculate();
 
        Array result(y.size(), 0.0);
 
        // Gauss Hermite
 
        Real stdDev_0_t = stateProcess_->stdDeviation(0.0, 0.0, t);
        // we use that the standard deviation is independent of $x$ here
        Real stdDev_0_T = stateProcess_->stdDeviation(0.0, 0.0, T);
        Real stdDev_t_T = stateProcess_->stdDeviation(t, 0.0, T - t);
 
        for (Size j = 0; j < y.size(); j++) {
            Array ya(modelSettings_.gaussHermitePoints_);
            for (Size i = 0; i < modelSettings_.gaussHermitePoints_; i++) {
                ya[i] = (y[j] * stdDev_0_t + stdDev_t_T * normalIntegralX_[i]) /
                        stdDev_0_T;
            }
            Array res = numeraireArray(T, ya);
            for (Size i = 0; i < modelSettings_.gaussHermitePoints_; i++) {
                result[j] += normalIntegralW_[i] / res[i];
            }
        }
 
        return result;
    }
 
    Real MarkovFunctional::numeraireImpl(
        const Time t, const Real y,
        const Handle<YieldTermStructure> &yts) const {
 
        if (t == 0)
            return yts.empty()
                       ? this->termStructure()->discount(numeraireTime(), true)
                       : yts->discount(numeraireTime());
 
        Array ya(1, y);
        return numeraireArray(t, ya)[0] *
               (yts.empty() ? Real(1.0)
                            : (yts->discount(numeraireTime()) /
                               yts->discount(t) * termStructure()->discount(t) /
                               termStructure()->discount(numeraireTime())));
    }
 
    Real
    MarkovFunctional::zerobondImpl(const Time T, const Time t, const Real y,
                                   const Handle<YieldTermStructure> &yts) const {
 
        if (t == 0.0)
            return yts.empty() ? this->termStructure()->discount(T, true)
                               : yts->discount(T, true);
        Array ya(1, y);
        return zerobondArray(T, t, ya)[0] *
               (yts.empty() ? Real(1.0) : (yts->discount(T) / yts->discount(t) *
                                     termStructure()->discount(t) /
                                     termStructure()->discount(T)));
    }
 
    Real MarkovFunctional::deflatedZerobond(Time T, Time t,
                                            Real y) const {
 
        Array ya(1, y);
        return deflatedZerobondArray(T, t, ya)[0];
    }
 
    Real MarkovFunctional::marketSwapRate(const Date &expiry,
                                          const CalibrationPoint &p,
                                          const Real digitalPrice,
                                          const Real guess,
                                          const Real shift) const {
 
        ZeroHelper z(this, expiry, p, digitalPrice);
        Brent b;
        Real solution = b.solve(
            z, modelSettings_.marketRateAccuracy_,
            std::max(std::min(guess, modelSettings_.upperRateBound_ - 0.00001),
                     modelSettings_.lowerRateBound_ - shift + 0.00001),
            modelSettings_.lowerRateBound_ - shift, modelSettings_.upperRateBound_);
        return solution;
    }
 
    Real MarkovFunctional::marketDigitalPrice(const Date &expiry,
                                              const CalibrationPoint &p,
                                              const Option::Type &type,
                                              const Real strike) const {
 
        return p.smileSection_->digitalOptionPrice(strike, type, p.annuity_,
                                                   modelSettings_.digitalGap_);
    }
 
    std::ostream &operator<<(std::ostream &out,
                             const MarkovFunctional::ModelOutputs &m) {
        out << "Markov functional model trace output " << std::endl;
        out << "Model settings" << std::endl;
        out << "Grid points y        : " << m.settings_.yGridPoints_
            << std::endl;
        out << "Std devs y           : " << m.settings_.yStdDevs_ << std::endl;
        out << "Lower rate bound     : " << m.settings_.lowerRateBound_
            << std::endl;
        out << "Upper rate bound     : " << m.settings_.upperRateBound_
            << std::endl;
        out << "Gauss Hermite points : " << m.settings_.gaussHermitePoints_
            << std::endl;
        out << "Digital gap          : " << m.settings_.digitalGap_
            << std::endl;
        out << "Adjustments          : "
            << ((m.settings_.adjustments_ & MarkovFunctional::ModelSettings::AdjustDigitals) != 0 ?
                    "Digitals " :
                    "")
            << ((m.settings_.adjustments_ & MarkovFunctional::ModelSettings::AdjustYts) != 0 ?
                    "Yts " :
                    "")
            << ((m.settings_.adjustments_ &
                 MarkovFunctional::ModelSettings::ExtrapolatePayoffFlat) != 0 ?
                    "FlatPayoffExt " :
                    "")
            << ((m.settings_.adjustments_ &
                 MarkovFunctional::ModelSettings::NoPayoffExtrapolation) != 0 ?
                    "NoPayoffExt " :
                    "")
            << ((m.settings_.adjustments_ & MarkovFunctional::ModelSettings::KahaleSmile) != 0 ?
                    "Kahale " :
                    "")
            << ((m.settings_.adjustments_ &
                 MarkovFunctional::ModelSettings::SmileExponentialExtrapolation) != 0 ?
                    "SmileExp " :
                    "")
            << ((m.settings_.adjustments_ & MarkovFunctional::ModelSettings::KahaleInterpolation) !=
                        0 ?
                    "KahaleInt " :
                    "")
            << ((m.settings_.adjustments_ &
                 MarkovFunctional::ModelSettings::SmileDeleteArbitragePoints) != 0 ?
                    "SmileDelArb " :
                    "")
            << ((m.settings_.adjustments_ & MarkovFunctional::ModelSettings::SabrSmile) != 0 ?
                    "Sabr" :
                    "")
            << std::endl;
        out << "Smile moneyness checkpoints: ";
        for (Size i = 0; i < m.settings_.smileMoneynessCheckpoints_.size(); i++)
            out << m.settings_.smileMoneynessCheckpoints_[i]
                << (i < m.settings_.smileMoneynessCheckpoints_.size() - 1 ? ";"
                                                                          : "");
        out << std::endl;
 
        QL_REQUIRE(!m.dirty_, "model outputs are dirty");
 
        if (m.expiries_.empty())
            return out; // no trace information was collected so no output
        out << std::endl;
        out << "Messages:" << std::endl;
        for (const auto& message : m.messages_)
            out << message << std::endl;
        out << std::endl << std::setprecision(16);
        out << "Yield termstructure fit:" << std::endl;
        out << "expiry;tenor;atm;annuity;digitalAdj;ytsAdj;marketzerorate;"
               "modelzerorate;diff(bp)" << std::endl;
        for (Size i = 0; i < m.expiries_.size(); i++) {
            out << m.expiries_[i] << ";" << m.tenors_[i] << ";" << m.atm_[i]
                << ";" << m.annuity_[i] << ";"
                << m.digitalsAdjustmentFactors_[i] << ";"
                << m.adjustmentFactors_[i] << ";" << m.marketZerorate_[i] << ";"
                << m.modelZerorate_[i] << ";"
                << (m.marketZerorate_[i] - m.modelZerorate_[i]) * 10000.0
                << std::endl;
        }
        out << std::endl;
        out << "Volatility smile fit:" << std::endl;
        for (Size i = 0; i < m.expiries_.size(); i++) {
            std::ostringstream os;
            os << m.expiries_[i] << "/" << m.tenors_[i];
            std::string p = os.str();
            out << "strike(" << p << ");marketCallRaw(" << p << ";marketCall("
                << p << ");modelCall(" << p << ");marketPutRaw(" << p
                << ");marketPut(" << p << ");modelPut(" << p << ");marketVega("
                << p << ")" << (i < m.expiries_.size() - 1 ? ";" : "");
        }
        out << std::endl;
        for (Size j = 0; j < m.smileStrikes_[0].size(); j++) {
            for (Size i = 0; i < m.expiries_.size(); i++) {
                out << m.smileStrikes_[i][j] << ";"
                    << m.marketRawCallPremium_[i][j] << ";"
                    << m.marketCallPremium_[i][j] << ";"
                    << m.modelCallPremium_[i][j] << ";"
                    << m.marketRawPutPremium_[i][j] << ";"
                    << m.marketPutPremium_[i][j] << ";"
                    << m.modelPutPremium_[i][j] << ";" << m.marketVega_[i][j]
                    << (i < m.expiries_.size() - 1 ? ";" : "");
            }
            out << std::endl;
        }
        return out;
    }
 
    Real MarkovFunctional::forwardRateInternal(
        const Date &fixing, const Date &referenceDate, const Real y,
        const bool zeroFixingDays, ext::shared_ptr<IborIndex> iborIdx) const {
 
        calculate();
 
        if (!iborIdx)
            iborIdx = iborIndex_;
 
        Date valueDate = zeroFixingDays ? fixing : iborIdx->valueDate(fixing);
        Date endDate = iborIdx->fixingCalendar().advance(
            iborIdx->valueDate(fixing), iborIdx->tenor(),
            iborIdx->businessDayConvention(),
            iborIdx->endOfMonth()); // FIXME Here we should use the calculation
                                    // date calendar ?
        Real dcf = iborIdx->dayCounter().yearFraction(valueDate, endDate);
 
        return (zerobond(valueDate, referenceDate, y) -
                zerobond(endDate, referenceDate, y)) /
               (dcf * zerobond(endDate, referenceDate, y));
    }
 
    Real
    MarkovFunctional::swapRateInternal(const Date &fixing, const Period &tenor,
                                       const Date &referenceDate, const Real y,
                                       bool zeroFixingDays,
                                       ext::shared_ptr<SwapIndex> swapIdx) const {
 
        calculate();
 
        if (!swapIdx)
            swapIdx = swapIndexBase_;
        QL_REQUIRE(swapIdx, "No swap index given");
 
        ext::shared_ptr<VanillaSwap> underlying = underlyingSwap(swapIdx, fixing, tenor);
 
        Schedule sched = underlying->fixedSchedule();
        Real annuity = swapAnnuityInternal(fixing, tenor, referenceDate, y,
                                      zeroFixingDays, swapIdx);
        Rate atm =
            (zerobond(zeroFixingDays ? fixing : sched.dates().front(),
                      referenceDate, y) -
             zerobond(sched.calendar().adjust(sched.dates().back(),
                                              underlying->paymentConvention()),
                      referenceDate, y)) /
            annuity;
        return atm;
    }
 
    Real MarkovFunctional::swapAnnuityInternal(
        const Date &fixing, const Period &tenor, const Date &referenceDate,
        const Real y, const bool zeroFixingDays,
        ext::shared_ptr<SwapIndex> swapIdx) const {
 
        calculate();
 
        if (!swapIdx)
            swapIdx = swapIndexBase_;
        QL_REQUIRE(swapIdx, "No swap index given");
 
        ext::shared_ptr<VanillaSwap> underlying = underlyingSwap(swapIdx, fixing, tenor);
 
        Schedule sched = underlying->fixedSchedule();
 
        Real annuity = 0.0;
        for (unsigned int j = 1; j < sched.size(); j++) {
            annuity +=
                zerobond(sched.calendar().adjust(
                             sched.date(j), underlying->paymentConvention()),
                         referenceDate, y) *
                swapIdx->dayCounter().yearFraction(
                    j == 1 && zeroFixingDays ? fixing : sched.date(j - 1),
                    sched.date(j));
        }
        return annuity;
    }
 
    Real MarkovFunctional::swaptionPriceInternal(const Option::Type& type,
                                                 const Date& expiry,
                                                 const Period& tenor,
                                                 const Rate strike,
                                                 const Date& referenceDate,
                                                 const Real y,
                                                 const bool zeroFixingDays,
                                                 const ext::shared_ptr<SwapIndex>& swapIdx) const {
 
        calculate();
 
        Time fixingTime = termStructure()->timeFromReference(expiry);
        Time referenceTime =
            referenceDate == Null<Date>()
                ? 0.0
                : termStructure()->timeFromReference(referenceDate);
 
        Array yg = yGrid(modelSettings_.yStdDevs_, modelSettings_.yGridPoints_,
                         fixingTime, referenceTime, y);
        Array z = yGrid(modelSettings_.yStdDevs_, modelSettings_.yGridPoints_);
        Array p(yg.size());
 
        for (Size i = 0; i < yg.size(); i++) {
            Real annuity = swapAnnuityInternal(expiry, tenor, expiry, yg[i],
                                          zeroFixingDays, swapIdx);
            Rate atm = swapRateInternal(expiry, tenor, expiry, yg[i], zeroFixingDays,
                                   swapIdx);
            p[i] = annuity * std::max((type == Option::Call ? 1.0 : -1.0) *
                                          (atm - strike),
                                      0.0) /
                   numeraire(fixingTime, yg[i]);
        }
 
        CubicInterpolation payoff(z.begin(), z.end(), p.begin(),
                                  CubicInterpolation::Spline, true,
                                  CubicInterpolation::Lagrange, 0.0,
                                  CubicInterpolation::Lagrange, 0.0);
 
        Real price = 0.0;
        for (Size i = 0; i < z.size() - 1; i++) {
            price += gaussianShiftedPolynomialIntegral(
                0.0, payoff.cCoefficients()[i], payoff.bCoefficients()[i],
                payoff.aCoefficients()[i], p[i], z[i], z[i], z[i + 1]);
        }
        if ((modelSettings_.adjustments_ &
             ModelSettings::NoPayoffExtrapolation) == 0) {
            if ((modelSettings_.adjustments_ &
                 ModelSettings::ExtrapolatePayoffFlat) != 0) {
                price += gaussianShiftedPolynomialIntegral(
                    0.0, 0.0, 0.0, 0.0, p[z.size() - 2], z[z.size() - 2],
                    z[z.size() - 1], 100.0);
                price += gaussianShiftedPolynomialIntegral(
                    0.0, 0.0, 0.0, 0.0, p[0], z[0], -100.0, z[0]);
            } else {
                if (type == Option::Call)
                    price += gaussianShiftedPolynomialIntegral(
                        0.0, payoff.cCoefficients()[z.size() - 2],
                        payoff.bCoefficients()[z.size() - 2],
                        payoff.aCoefficients()[z.size() - 2], p[z.size() - 2],
                        z[z.size() - 2], z[z.size() - 1], 100.0);
                if (type == Option::Put)
                    price += gaussianShiftedPolynomialIntegral(
                        0.0, payoff.cCoefficients()[0],
                        payoff.bCoefficients()[0], payoff.aCoefficients()[0],
                        p[0], z[0], -100.0, z[0]);
            }
        }
 
        return numeraire(referenceTime, y) * price;
    }
 
    Real MarkovFunctional::capletPriceInternal(
        const Option::Type &type, const Date &expiry, const Rate strike,
        const Date &referenceDate, const Real y, const bool zeroFixingDays,
        ext::shared_ptr<IborIndex> iborIdx) const {
 
        calculate();
 
        if (!iborIdx)
            iborIdx = iborIndex_;
 
        Time fixingTime = termStructure()->timeFromReference(expiry);
        Time referenceTime =
            referenceDate == Null<Date>()
                ? 0.0
                : termStructure()->timeFromReference(referenceDate);
 
        Array yg = yGrid(modelSettings_.yStdDevs_, modelSettings_.yGridPoints_,
                         fixingTime, referenceTime, y);
        Array z = yGrid(modelSettings_.yStdDevs_, modelSettings_.yGridPoints_);
        Array p(yg.size());
 
        Date valueDate = iborIdx->valueDate(expiry);
        Date endDate = iborIdx->fixingCalendar().advance(
            valueDate, iborIdx->tenor(), iborIdx->businessDayConvention(),
            iborIdx->endOfMonth()); // FIXME Here we should use the calculation
                                    // date calendar ?
        Real dcf = iborIdx->dayCounter().yearFraction(
            zeroFixingDays ? expiry : valueDate, endDate);
 
        for (Size i = 0; i < yg.size(); i++) {
            Real annuity = zerobond(endDate, expiry, yg[i]) * dcf;
            Rate atm =
                forwardRateInternal(expiry, expiry, yg[i], zeroFixingDays, iborIdx);
            p[i] = annuity * std::max((type == Option::Call ? 1.0 : -1.0) *
                                          (atm - strike),
                                      0.0) /
                   numeraire(fixingTime, yg[i]);
        }
 
        CubicInterpolation payoff(z.begin(), z.end(), p.begin(),
                                  CubicInterpolation::Spline, true,
                                  CubicInterpolation::Lagrange, 0.0,
                                  CubicInterpolation::Lagrange, 0.0);
 
        Real price = 0.0;
        for (Size i = 0; i < z.size() - 1; i++) {
            price += gaussianShiftedPolynomialIntegral(
                0.0, payoff.cCoefficients()[i], payoff.bCoefficients()[i],
                payoff.aCoefficients()[i], p[i], z[i], z[i], z[i + 1]);
        }
        if ((modelSettings_.adjustments_ &
             ModelSettings::NoPayoffExtrapolation) == 0) {
            if ((modelSettings_.adjustments_ &
                 ModelSettings::ExtrapolatePayoffFlat) != 0) {
                price += gaussianShiftedPolynomialIntegral(
                    0.0, 0.0, 0.0, 0.0, p[z.size() - 2], z[z.size() - 2],
                    z[z.size() - 1], 100.0);
                price += gaussianShiftedPolynomialIntegral(
                    0.0, 0.0, 0.0, 0.0, p[0], z[0], -100.0, z[0]);
            } else {
                if (type == Option::Call)
                    price += gaussianShiftedPolynomialIntegral(
                        0.0, payoff.cCoefficients()[z.size() - 2],
                        payoff.bCoefficients()[z.size() - 2],
                        payoff.aCoefficients()[z.size() - 2], p[z.size() - 2],
                        z[z.size() - 2], z[z.size() - 1], 100.0);
                if (type == Option::Put)
                    price += gaussianShiftedPolynomialIntegral(
                        0.0, payoff.cCoefficients()[0],
                        payoff.bCoefficients()[0], payoff.aCoefficients()[0],
                        p[0], z[0], -100.0, z[0]);
            }
        }
 
        return numeraire(referenceTime, y) * price;
    }
}