gaussian1dmodel.cpp

/* -*- 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
 <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/models/shortrate/onefactormodels/gaussian1dmodel.hpp>
#include <ql/math/interpolations/cubicinterpolation.hpp>
#include <ql/payoff.hpp>
#include <cmath>
 
namespace QuantLib {
 
    Real Gaussian1dModel::forwardRate(const Date& fixing,
                                      const Date& referenceDate,
                                      const Real y,
                                      const ext::shared_ptr<IborIndex>& iborIdx) const {
 
        QL_REQUIRE(iborIdx != nullptr, "no ibor index given");
 
        calculate();
 
        if (fixing <= (evaluationDate_ + (enforcesTodaysHistoricFixings_ ? 0 : -1)))
            return iborIdx->fixing(fixing);
 
        Handle<YieldTermStructure> yts = iborIdx->forwardingTermStructure(); // might be empty, then
                                                                             // use model curve
 
        Date valueDate = iborIdx->valueDate(fixing);
        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(valueDate, endDate);
 
        return (zerobond(valueDate, referenceDate, y, yts) -
                zerobond(endDate, referenceDate, y, yts)) /
               (dcf * zerobond(endDate, referenceDate, y, yts));
}
 
Real Gaussian1dModel::swapRate(const Date& fixing,
                               const Period& tenor,
                               const Date& referenceDate,
                               const Real y,
                               const ext::shared_ptr<SwapIndex>& swapIdx) const {
 
    QL_REQUIRE(swapIdx != nullptr, "no swap index given");
 
    calculate();
 
    if (fixing <= (evaluationDate_ + (enforcesTodaysHistoricFixings_ ? 0 : -1)))
        return swapIdx->fixing(fixing);
 
    Handle<YieldTermStructure> ytsf =
        swapIdx->iborIndex()->forwardingTermStructure();
    Handle<YieldTermStructure> ytsd =
        swapIdx->discountingTermStructure(); // either might be empty, then
                                             // use model curve
 
    Schedule sched, floatSched;
 
    ext::shared_ptr<VanillaSwap> underlying =
        underlyingSwap(swapIdx, fixing, tenor);
 
    sched = underlying->fixedSchedule();
 
    ext::shared_ptr<OvernightIndexedSwapIndex> oisIdx =
        ext::dynamic_pointer_cast<OvernightIndexedSwapIndex>(swapIdx);
    if (oisIdx != nullptr) {
        floatSched = sched;
    } else {
        floatSched = underlying->floatingSchedule();
    }
 
    Real annuity = swapAnnuity(fixing, tenor, referenceDate, y,
                               swapIdx);  // should be fine for
                                          // overnightindexed swap indices as
                                          // well
    Rate floatleg = 0.0;
    if (ytsf.empty() && ytsd.empty()) { // simple 100-formula can be used
                                        // only in one curve setup
        floatleg =
            (zerobond(sched.dates().front(), referenceDate, y) -
             zerobond(sched.calendar().adjust(sched.dates().back(),
                                              underlying->paymentConvention()),
                      referenceDate, y));
    } else {
        for (Size i = 1; i < floatSched.size(); i++) {
            floatleg +=
                (zerobond(floatSched[i - 1], referenceDate, y, ytsf) /
                     zerobond(floatSched[i], referenceDate, y, ytsf) -
                 1.0) *
                zerobond(floatSched.calendar().adjust(
                             floatSched[i], underlying->paymentConvention()),
                         referenceDate, y, ytsd);
        }
    }
    return floatleg / annuity;
}
 
Real Gaussian1dModel::swapAnnuity(const Date& fixing,
                                  const Period& tenor,
                                  const Date& referenceDate,
                                  const Real y,
                                  const ext::shared_ptr<SwapIndex>& swapIdx) const {
 
    QL_REQUIRE(swapIdx != nullptr, "no swap index given");
 
    calculate();
 
    Handle<YieldTermStructure> ytsd =
        swapIdx->discountingTermStructure(); // might be empty, then use
                                             // model curve
 
    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, ytsd) *
                   swapIdx->dayCounter().yearFraction(sched.date(j - 1),
                                                      sched.date(j));
    }
    return annuity;
}
 
Real Gaussian1dModel::zerobondOption(
    const Option::Type &type, const Date &expiry, const Date &valueDate,
    const Date &maturity, const Rate strike, const Date &referenceDate,
    const Real y, const Handle<YieldTermStructure> &yts, const Real yStdDevs,
    const Size yGridPoints, const bool extrapolatePayoff,
    const bool flatPayoffExtrapolation) const {
 
    calculate();
 
    Time fixingTime = termStructure()->timeFromReference(expiry);
    Time referenceTime =
        referenceDate == Null<Date>()
            ? 0.0
            : termStructure()->timeFromReference(referenceDate);
 
    Array yg = yGrid(yStdDevs, yGridPoints, fixingTime, referenceTime, y);
    Array z = yGrid(yStdDevs, yGridPoints);
 
    Array p(yg.size());
 
    for (Size i = 0; i < yg.size(); i++) {
        Real expValDsc = zerobond(valueDate, expiry, yg[i], yts);
        Real discount =
            zerobond(maturity, expiry, yg[i], yts) / expValDsc;
        p[i] =
            std::max((type == Option::Call ? 1.0 : -1.0) * (discount - strike),
                     0.0) /
            numeraire(fixingTime, yg[i], yts) * expValDsc;
    }
 
    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 (extrapolatePayoff) {
        if (flatPayoffExtrapolation) {
            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, yts) * price;
}
 
Real Gaussian1dModel::gaussianPolynomialIntegral(
    const Real a, const Real b, const Real c, const Real d, const Real e,
    const Real y0, const Real y1) {
 
#ifdef GAUSS1D_ENABLE_NTL
    const boost::math::ntl::RR aa = 4.0 * a, ba = 2.0 * M_SQRT2 * b,
                               ca = 2.0 * c, da = M_SQRT2 * d;
    const boost::math::ntl::RR x0 = y0 * M_SQRT1_2, x1 = y1 * M_SQRT1_2;
    const boost::math::ntl::RR res =
        (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x1) -
         1.0 / (4.0 * M_SQRTPI) * std::exp(-x1 * x1) *
             (2.0 * aa * x1 * x1 * x1 + 3.0 * aa * x1 +
              2.0 * ba * (x1 * x1 + 1.0) + 2.0 * ca * x1 + 2.0 * da)) -
        (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x0) -
         1.0 / (4.0 * M_SQRTPI) * std::exp(-x0 * x0) *
             (2.0 * aa * x0 * x0 * x0 + 3.0 * aa * x0 +
              2.0 * ba * (x0 * x0 + 1.0) + 2.0 * ca * x0 + 2.0 * da));
    return NTL::to_double(res.value());
#else
    const Real aa = 4.0 * a, ba = 2.0 * M_SQRT2 * b, ca = 2.0 * c,
               da = M_SQRT2 * d;
    const Real x0 = y0 * M_SQRT1_2, x1 = y1 * M_SQRT1_2;
    return (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x1) -
            1.0 / (4.0 * M_SQRTPI) * std::exp(-x1 * x1) *
                (2.0 * aa * x1 * x1 * x1 + 3.0 * aa * x1 +
                 2.0 * ba * (x1 * x1 + 1.0) + 2.0 * ca * x1 + 2.0 * da)) -
           (0.125 * (3.0 * aa + 2.0 * ca + 4.0 * e) * std::erf(x0) -
            1.0 / (4.0 * M_SQRTPI) * std::exp(-x0 * x0) *
                (2.0 * aa * x0 * x0 * x0 + 3.0 * aa * x0 +
                 2.0 * ba * (x0 * x0 + 1.0) + 2.0 * ca * x0 + 2.0 * da));
#endif
}
 
Real Gaussian1dModel::gaussianShiftedPolynomialIntegral(
    const Real a, const Real b, const Real c, const Real d, const Real e,
    const Real h, const Real x0, const Real x1) {
    return gaussianPolynomialIntegral(
        a, -4.0 * a * h + b, 6.0 * a * h * h - 3.0 * b * h + c,
        -4 * a * h * h * h + 3.0 * b * h * h - 2.0 * c * h + d,
        a * h * h * h * h - b * h * h * h + c * h * h - d * h + e, x0, x1);
}
 
Array Gaussian1dModel::yGrid(
    const Real stdDevs, const int gridPoints, const Real T, const Real t, const Real y) const {
 
    // we use that the standard deviation is independent of $x$ here !
 
    QL_REQUIRE(stateProcess_ != nullptr, "state process not set");
 
    Array result(2 * gridPoints + 1, 0.0);
 
    Real x_t, e_0_t, e_t_T, stdDev_0_t, stdDev_t_T;
    Real stdDev_0_T = stateProcess_->stdDeviation(0.0, 0.0, T);
    Real e_0_T = stateProcess_->expectation(0.0, 0.0, T);
 
    if (t < QL_EPSILON) {
        // stdDev_0_t = 0.0;
        stdDev_t_T = stdDev_0_T;
        // e_0_t = 0.0;
        // x_t = 0.0;
        e_t_T = e_0_T;
    } else {
        stdDev_0_t = stateProcess_->stdDeviation(0.0, 0.0, t);
        stdDev_t_T = stateProcess_->stdDeviation(t, 0.0, T - t);
        e_0_t = stateProcess_->expectation(0.0, 0.0, t);
        x_t = y * stdDev_0_t + e_0_t;
        e_t_T = stateProcess_->expectation(t, x_t, T - t);
    }
 
    Real h = stdDevs / ((Real)gridPoints);
 
    for (int j = -gridPoints; j <= gridPoints; j++) {
        result[j + gridPoints] =
            (e_t_T + stdDev_t_T * ((Real)j) * h - e_0_T) / stdDev_0_T;
    }
 
    return result;
}
}