integralntdengine.cpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2008 Roland Lichters
 
 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/experimental/credit/integralntdengine.hpp>
#include <ql/cashflows/fixedratecoupon.hpp>
#include <ql/termstructures/yieldtermstructure.hpp>
#include <ql/experimental/credit/basket.hpp>
#include <numeric>
 
namespace QuantLib {
 
    void IntegralNtdEngine::calculate() const {
        Date today = Settings::instance().evaluationDate();
 
        results_.errorEstimate = Null<Real>();
        results_.value = 0.0;
        results_.premiumValue = 0.0;
        results_.upfrontPremiumValue = 0.;
        Real accrualValue = 0.0;
        Real claimValue = 0.0;
        Date d0;
        /* Given the expense of probsBeingNthEvent both in integrable and 
        monte carlo algorithms this engine tests who to call.
        Warning: This is not entirely a basket property but of the model too.
        The basket has to have all notionals equal but it is the model which
        determines the recovery; having all the market recoveries equal is not
        enough since we might be using a loss model which is stochastic in the
        recovery rates.
        */
        bool basketIsHomogeneous = true;// hardcoded by now
 
        for (auto& i : arguments_.premiumLeg) {
            ext::shared_ptr<FixedRateCoupon> coupon = ext::dynamic_pointer_cast<FixedRateCoupon>(i);
            Date d = i->date();
            if (d > discountCurve_->referenceDate()) {
                /*
                std::vector<Probability> probsTriggering =
                    arguments_.basket->probsBeingNthEvent(arguments_.ntdOrder, 
                        d);
                Probability defaultProb = 
                    std::accumulate(probsTriggering.begin(), 
                    probsTriggering.end(), Real(0.));
                // OVERKILL???? 1-probAtLeastNEvents is enough
 
*/
                // prob of contract not having been triggered by date of payment
                Probability probNonTriggered = 
                    1. - arguments_.basket->probAtLeastNEvents(
                        arguments_.ntdOrder, d);
 
                results_.premiumValue +=
                    i->amount() * discountCurve_->discount(d) * probNonTriggered;
 
                if (coupon->accrualStartDate() >= 
                    discountCurve_->referenceDate())
                    d = coupon->accrualStartDate();
                else
                    d = discountCurve_->referenceDate();
 
                // do steps of specified size
                d0 = d;
                Period stepSize = integrationStepSize_;
/*
                probsTriggering =
                    arguments_.basket->probsBeingNthEvent(arguments_.ntdOrder, 
                        d0);
                Probability defProb0 = std::accumulate(probsTriggering.begin(), 
                    probsTriggering.end(), Real(0.));
*/
                Probability defProb0 = arguments_.basket->probAtLeastNEvents(
                        arguments_.ntdOrder, d0);
                std::vector<Probability> probsTriggering, probsTriggering1;
                do {
                    DiscountFactor disc = discountCurve_->discount(d);
 
                    Probability defProb1;
                    if(basketIsHomogeneous) {//take test out of the while loop
                        defProb1 = arguments_.basket->probAtLeastNEvents(
                            arguments_.ntdOrder, d);
                        claimValue -= (defProb1-defProb0)
                            * arguments_.basket->claim()->amount(d, 
                                arguments_.notional, 
                                arguments_.basket->recoveryRate(d, 0))
                            * disc;
 
                    }else{
                        probsTriggering1 =
                            arguments_.basket->probsBeingNthEvent(
                                arguments_.ntdOrder, d);
                        defProb1 = std::accumulate(probsTriggering1.begin(), 
                            probsTriggering1.end(), Real(0.));
                        /*Recoveries might differ along names, depending on 
                        which name is triggering the contract the loss will be 
                        different  
                        There is an issue here; MC engines can still be used 
                        since the prob of triggering the contract can be 
                        extracted from the simulation from the 
                        probsBeingNthEvent statistic. Yet, when the RR is 
                        stochastic the realized value of the RR is the expected 
                        one subject/conditional to the contract being triggered;
                        not simply the expected value. For this reason the MC 
                        can not be used through the statistic but has to consume
                        the simulations directly.
                        */
                        for(Size iName=0; 
                            iName<arguments_.basket->remainingSize(); 
                            iName++) 
                        {
                            claimValue -= (probsTriggering1[iName]-
                                probsTriggering[iName])
                                * arguments_.basket->claim()->amount(d, 
                                    arguments_.notional,// [iName]! 
                                    arguments_.basket->recoveryRate(d, iName))
                                * disc;
                        }
                        probsTriggering = probsTriggering1;
                    }
 
                    Probability dcfdd = defProb1 - defProb0;
                    defProb0 = defProb1;
 
                    if (arguments_.settlePremiumAccrual)
                        accrualValue += coupon->accruedAmount(d)*disc*dcfdd;
 
                    d0 = d;
                    d = d0 + stepSize;
                    // reduce step size ?
                    if (stepSize != 1*Days && d > coupon->accrualEndDate()) {
                        stepSize = 1*Days;
                        d = d0 + stepSize;
                    }
                }
                while (d <= coupon->accrualEndDate());
            }
        }
 
        // The upfront might be due before the curve ref date...
        if (!arguments_.premiumLeg[0]->hasOccurred(today))
            results_.upfrontPremiumValue =
                arguments_.basket->remainingNotional() 
                    * arguments_.upfrontRate
                    * discountCurve_->discount(
                        ext::dynamic_pointer_cast<FixedRateCoupon>(
                            arguments_.premiumLeg[0])->accrualStartDate());
        if (arguments_.side == Protection::Buyer) {
            results_.premiumValue *= -1;
            accrualValue *= -1;
            claimValue *= -1;
            results_.upfrontPremiumValue *= -1;
        }
 
        results_.value = results_.premiumValue + accrualValue + claimValue + 
            results_.upfrontPremiumValue;
 
        results_.fairPremium = -arguments_.premiumRate * claimValue 
            / (results_.premiumValue + accrualValue);
        // alternatively use results buffers and omit locals.
        results_.protectionValue = claimValue;
 
        results_.additionalResults["fairPremium"] = results_.fairPremium;
        results_.additionalResults["premiumLegNPV"] = 
            Real(results_.premiumValue + results_.upfrontPremiumValue);
        results_.additionalResults["protectionLegNPV"] = 
            results_.protectionValue;
    }
 
}