lognormalcmsspreadpricer.cpp

Go to the documentation of this file.

/*
 Copyright (C) 2018 Quaternion Risk Management Ltd
 All rights reserved.
 
 This file is part of ORE, a free-software/open-source library
 for transparent pricing and risk analysis - http://opensourcerisk.org
 
 ORE is free software: you can redistribute it and/or modify it
 under the terms of the Modified BSD License.  You should have received a
 copy of the license along with this program.
 The license is also available online at <http://opensourcerisk.org>
 
 This program is distributed on the basis that it will form a useful
 contribution to risk analytics and model standardisation, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
*/
 
/*
  Copyright (C) 2014, 2015, 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 <qle/cashflows/lognormalcmsspreadpricer.hpp>
 
#include <ql/experimental/coupons/cmsspreadcoupon.hpp>
#include <ql/math/integrals/kronrodintegral.hpp>
#include <ql/pricingengines/blackformula.hpp>
#include <ql/termstructures/volatility/swaption/swaptionvolcube.hpp>
 
#include <boost/make_shared.hpp>
 
using std::sqrt;
 
using namespace QuantLib;
 
namespace QuantExt {
 
class LognormalCmsSpreadPricer::integrand_f {
    const LognormalCmsSpreadPricer* pricer;
 
public:
    explicit integrand_f(const LognormalCmsSpreadPricer* pricer) : pricer(pricer) {}
    Real operator()(Real x) const { return pricer->integrand(x); }
};
 
LognormalCmsSpreadPricer::LognormalCmsSpreadPricer(const QuantLib::ext::shared_ptr<CmsCouponPricer> cmsPricer,
                                                   const Handle<QuantExt::CorrelationTermStructure>& correlation,
                                                   const Handle<YieldTermStructure>& couponDiscountCurve,
                                                   const Size integrationPoints,
                                                   const boost::optional<VolatilityType> volatilityType,
                                                   const Real shift1, const Real shift2)
    : CmsSpreadCouponPricer2(correlation), cmsPricer_(cmsPricer), couponDiscountCurve_(couponDiscountCurve) {
 
    if (!couponDiscountCurve_.empty())
        registerWith(couponDiscountCurve_);
    registerWith(cmsPricer_);
 
    QL_REQUIRE(integrationPoints >= 4, "at least 4 integration points should be used (" << integrationPoints << ")");
    integrator_ = QuantLib::ext::make_shared<GaussHermiteIntegration>(integrationPoints);
 
    cnd_ = QuantLib::ext::make_shared<CumulativeNormalDistribution>(0.0, 1.0);
 
    if (volatilityType == boost::none) {
        QL_REQUIRE(shift1 == Null<Real>() && shift2 == Null<Real>(),
                   "if volatility type is inherited, no shifts should be "
                   "specified");
        inheritedVolatilityType_ = true;
        volType_ = cmsPricer->swaptionVolatility()->volatilityType();
    } else {
        shift1_ = shift1 == Null<Real>() ? 0.0 : shift1;
        shift2_ = shift2 == Null<Real>() ? 0.0 : shift2;
        inheritedVolatilityType_ = false;
        volType_ = *volatilityType;
    }
}
 
Real LognormalCmsSpreadPricer::integrand(const Real x) const {
 
    // this is Brigo, 13.16.2 with x = v / sqrt(2)
 
    Real v = M_SQRT2 * x;
    Real h = k_ - b_ * s2_ * std::exp((m2_ - 0.5 * v2_ * v2_) * fixingTime_ + v2_ * std::sqrt(fixingTime_) * v);
    Real phi1, phi2;
    phi1 = (*cnd_)(phi_ *
                   (std::log(a_ * s1_ / h) + (m1_ + (0.5 - rho() * rho()) * v1_ * v1_) * fixingTime_ +
                    rho() * v1_ * std::sqrt(fixingTime_) * v) /
                   (v1_ * std::sqrt(fixingTime_ * (1.0 - rho() * rho()))));
    phi2 = (*cnd_)(
        phi_ *
        (std::log(a_ * s1_ / h) + (m1_ - 0.5 * v1_ * v1_) * fixingTime_ + rho() * v1_ * std::sqrt(fixingTime_) * v) /
        (v1_ * std::sqrt(fixingTime_ * (1.0 - rho() * rho()))));
    Real f = a_ * phi_ * s1_ *
                 std::exp(m1_ * fixingTime_ - 0.5 * rho() * rho() * v1_ * v1_ * fixingTime_ +
                          rho() * v1_ * std::sqrt(fixingTime_) * v) *
                 phi1 -
             phi_ * h * phi2;
    return std::exp(-x * x) * f;
}
 
Real LognormalCmsSpreadPricer::integrand_normal(const Real x) const {
 
    // this is http://ssrn.com/abstract=2686998, 3.20 with x = s / sqrt(2)
 
    Real s = M_SQRT2 * x;
 
    Real beta = phi_ * (gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_ - k_ +
                        std::sqrt(fixingTime_) * (rho() * gearing1_ * vol1_ + gearing2_ * vol2_) * s);
    Real f = close_enough(alpha_, 0.0)
                 ? std::max(beta, 0.0)
                 : psi_ * alpha_ / (M_SQRTPI * M_SQRT2) * std::exp(-beta * beta / (2.0 * alpha_ * alpha_)) +
                       beta * (1.0 - (*cnd_)(-psi_ * beta / alpha_));
    return std::exp(-x * x) * f;
}
 
void LognormalCmsSpreadPricer::initialize(const FloatingRateCoupon& coupon) {
 
    coupon_ = dynamic_cast<const CmsSpreadCoupon*>(&coupon);
    QL_REQUIRE(coupon_, "CMS spread coupon needed");
    index_ = coupon_->swapSpreadIndex();
    gearing_ = coupon_->gearing();
    spread_ = coupon_->spread();
 
    fixingDate_ = coupon_->fixingDate();
    paymentDate_ = coupon_->date();
 
    // if no coupon discount curve is given just use the discounting curve
    // from the _first_ swap index.
    // for rate calculation this curve cancels out in the computation, so
    // e.g. the discounting
    // swap engine will produce correct results, even if the
    // couponDiscountCurve is not set here.
    // only the price member function in this class will be dependent on the
    // coupon discount curve.
 
    today_ = QuantLib::Settings::instance().evaluationDate();
 
    if (couponDiscountCurve_.empty())
        couponDiscountCurve_ = index_->swapIndex1()->exogenousDiscount()
                                   ? index_->swapIndex1()->discountingTermStructure()
                                   : index_->swapIndex1()->forwardingTermStructure();
 
    discount_ =
        paymentDate_ > couponDiscountCurve_->referenceDate() ? couponDiscountCurve_->discount(paymentDate_) : 1.0;
 
    spreadLegValue_ = spread_ * coupon_->accrualPeriod() * discount_;
 
    gearing1_ = index_->gearing1();
    gearing2_ = index_->gearing2();
 
    QL_REQUIRE(gearing1_ > 0.0 && gearing2_ < 0.0, "gearing1 (" << gearing1_ << ") should be positive while gearing2 ("
                                                                << gearing2_ << ") should be negative");
 
    c1_ = QuantLib::ext::shared_ptr<CmsCoupon>(
        new CmsCoupon(coupon_->date(), coupon_->nominal(), coupon_->accrualStartDate(), coupon_->accrualEndDate(),
                      coupon_->fixingDays(), index_->swapIndex1(), 1.0, 0.0, coupon_->referencePeriodStart(),
                      coupon_->referencePeriodEnd(), coupon_->dayCounter(), coupon_->isInArrears()));
 
    c2_ = QuantLib::ext::shared_ptr<CmsCoupon>(
        new CmsCoupon(coupon_->date(), coupon_->nominal(), coupon_->accrualStartDate(), coupon_->accrualEndDate(),
                      coupon_->fixingDays(), index_->swapIndex2(), 1.0, 0.0, coupon_->referencePeriodStart(),
                      coupon_->referencePeriodEnd(), coupon_->dayCounter(), coupon_->isInArrears()));
 
    c1_->setPricer(cmsPricer_);
    c2_->setPricer(cmsPricer_);
 
    if (fixingDate_ > today_) {
 
        fixingTime_ = cmsPricer_->swaptionVolatility()->timeFromReference(fixingDate_);
 
        swapRate1_ = c1_->indexFixing();
        swapRate2_ = c2_->indexFixing();
 
        adjustedRate1_ = c1_->adjustedFixing();
        adjustedRate2_ = c2_->adjustedFixing();
 
        QuantLib::ext::shared_ptr<SwaptionVolatilityStructure> swvol = *cmsPricer_->swaptionVolatility();
        QuantLib::ext::shared_ptr<SwaptionVolatilityCube> swcub = QuantLib::ext::dynamic_pointer_cast<SwaptionVolatilityCube>(swvol);
 
        if (inheritedVolatilityType_ && volType_ == ShiftedLognormal) {
            shift1_ = swvol->shift(fixingDate_, index_->swapIndex1()->tenor());
            shift2_ = swvol->shift(fixingDate_, index_->swapIndex2()->tenor());
        }
 
        if (swcub == NULL) {
            // not a cube, just an atm surface given, so we can
            // not easily convert volatilities and just forbid it
            QL_REQUIRE(inheritedVolatilityType_, "if only an atm surface is given, the volatility "
                                                 "type must be inherited");
            vol1_ = swvol->volatility(fixingDate_, index_->swapIndex1()->tenor(), swapRate1_);
            vol2_ = swvol->volatility(fixingDate_, index_->swapIndex2()->tenor(), swapRate2_);
        } else {
            vol1_ = swcub->smileSection(fixingDate_, index_->swapIndex1()->tenor())
                        ->volatility(swapRate1_, volType_, shift1_);
            vol2_ = swcub->smileSection(fixingDate_, index_->swapIndex2()->tenor())
                        ->volatility(swapRate2_, volType_, shift2_);
        }
 
        if (volType_ == ShiftedLognormal) {
            mu1_ = 1.0 / fixingTime_ * std::log((adjustedRate1_ + shift1_) / (swapRate1_ + shift1_));
            mu2_ = 1.0 / fixingTime_ * std::log((adjustedRate2_ + shift2_) / (swapRate2_ + shift2_));
        }
        // for the normal volatility case we do not need the drifts
        // but rather use adjusted rates directly in the integrand
 
    } else {
        // fixing is in the past or today
        adjustedRate1_ = c1_->indexFixing();
        adjustedRate2_ = c2_->indexFixing();
    }
}
 
Real LognormalCmsSpreadPricer::optionletPrice(Option::Type optionType, Real strike) const {
    // this method is only called for future fixings
    optionType_ = optionType;
    phi_ = optionType == Option::Call ? 1.0 : -1.0;
    Real res = 0.0;
    if (volType_ == ShiftedLognormal) {
        // (shifted) lognormal volatility
        if (strike >= 0.0) {
            a_ = gearing1_;
            b_ = gearing2_;
            s1_ = swapRate1_ + shift1_;
            s2_ = swapRate2_ + shift2_;
            m1_ = mu1_;
            m2_ = mu2_;
            v1_ = vol1_;
            v2_ = vol2_;
            k_ = strike + gearing1_ * shift1_ + gearing2_ * shift2_;
        } else {
            a_ = -gearing2_;
            b_ = -gearing1_;
            s1_ = swapRate2_ + shift1_;
            s2_ = swapRate1_ + shift2_;
            m1_ = mu2_;
            m2_ = mu1_;
            v1_ = vol2_;
            v2_ = vol1_;
            k_ = -strike - gearing1_ * shift1_ - gearing2_ * shift2_;
            res += phi_ * (gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_ - strike);
        }
        res += 1.0 / M_SQRTPI * (*integrator_)(integrand_f(this));
    } else {
        // normal volatility
        Real forward = gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_;
        Real stddev =
            std::sqrt(fixingTime_ * (gearing1_ * gearing1_ * vol1_ * vol1_ + gearing2_ * gearing2_ * vol2_ * vol2_ +
                                     2.0 * gearing1_ * gearing2_ * rho() * vol1_ * vol2_));
        res = bachelierBlackFormula(optionType_, strike, forward, stddev, 1.0);
    }
    return res * discount_ * coupon_->accrualPeriod();
}
 
Rate LognormalCmsSpreadPricer::swapletRate() const { return swapletPrice() / (coupon_->accrualPeriod() * discount_); }
 
Real LognormalCmsSpreadPricer::capletPrice(Rate effectiveCap) const {
    // caplet is equivalent to call option on fixing
    if (fixingDate_ <= today_) {
        // the fixing is determined
        const Rate Rs = std::max(coupon_->index()->fixing(fixingDate_) - effectiveCap, 0.);
        Rate price = gearing_ * Rs * coupon_->accrualPeriod() * discount_;
        return price;
    } else {
        Real capletPrice = optionletPrice(Option::Call, effectiveCap);
        return gearing_ * capletPrice;
    }
}
 
Rate LognormalCmsSpreadPricer::capletRate(Rate effectiveCap) const {
    return capletPrice(effectiveCap) / (coupon_->accrualPeriod() * discount_);
}
 
Real LognormalCmsSpreadPricer::floorletPrice(Rate effectiveFloor) const {
    // floorlet is equivalent to put option on fixing
    if (fixingDate_ <= today_) {
        // the fixing is determined
        const Rate Rs = std::max(effectiveFloor - coupon_->index()->fixing(fixingDate_), 0.);
        Rate price = gearing_ * Rs * coupon_->accrualPeriod() * discount_;
        return price;
    } else {
        Real floorletPrice = optionletPrice(Option::Put, effectiveFloor);
        return gearing_ * floorletPrice;
    }
}
 
Rate LognormalCmsSpreadPricer::floorletRate(Rate effectiveFloor) const {
    return floorletPrice(effectiveFloor) / (coupon_->accrualPeriod() * discount_);
}
 
Real LognormalCmsSpreadPricer::swapletPrice() const {
    return gearing_ * coupon_->accrualPeriod() * discount_ * (gearing1_ * adjustedRate1_ + gearing2_ * adjustedRate2_) +
           spreadLegValue_;
}
} // namespace QuantExt