Coverage Report

Created: 2025-08-05 06:45

/src/quantlib/ql/pricingengines/credit/isdacdsengine.cpp
Line
Count
Source (jump to first uncovered line)
1
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3
/*
4
 Copyright (C) 2014 Jose Aparicio
5
 Copyright (C) 2014 Peter Caspers
6
7
 This file is part of QuantLib, a free-software/open-source library
8
 for financial quantitative analysts and developers - http://quantlib.org/
9
10
 QuantLib is free software: you can redistribute it and/or modify it
11
 under the terms of the QuantLib license.  You should have received a
12
 copy of the license along with this program; if not, please email
13
 <quantlib-dev@lists.sf.net>. The license is also available online at
14
 <http://quantlib.org/license.shtml>.
15
16
 This program is distributed in the hope that it will be useful, but WITHOUT
17
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18
 FOR A PARTICULAR PURPOSE.  See the license for more details.
19
*/
20
21
#include <ql/cashflows/fixedratecoupon.hpp>
22
#include <ql/instruments/claim.hpp>
23
#include <ql/math/interpolations/forwardflatinterpolation.hpp>
24
#include <ql/pricingengines/credit/isdacdsengine.hpp>
25
#include <ql/termstructures/credit/flathazardrate.hpp>
26
#include <ql/termstructures/credit/piecewisedefaultcurve.hpp>
27
#include <ql/termstructures/yield/flatforward.hpp>
28
#include <ql/termstructures/yield/piecewiseyieldcurve.hpp>
29
#include <ql/time/calendars/weekendsonly.hpp>
30
#include <ql/time/daycounters/actual360.hpp>
31
#include <ql/optional.hpp>
32
#include <utility>
33
34
namespace QuantLib {
35
36
    IsdaCdsEngine::IsdaCdsEngine(Handle<DefaultProbabilityTermStructure> probability,
37
                                 Real recoveryRate,
38
                                 Handle<YieldTermStructure> discountCurve,
39
                                 const ext::optional<bool>& includeSettlementDateFlows,
40
                                 const NumericalFix numericalFix,
41
                                 const AccrualBias accrualBias,
42
                                 const ForwardsInCouponPeriod forwardsInCouponPeriod)
43
0
    : probability_(std::move(probability)), recoveryRate_(recoveryRate),
44
0
      discountCurve_(std::move(discountCurve)),
45
0
      includeSettlementDateFlows_(includeSettlementDateFlows), numericalFix_(numericalFix),
46
0
      accrualBias_(accrualBias), forwardsInCouponPeriod_(forwardsInCouponPeriod) {
47
48
0
        registerWith(probability_);
49
0
        registerWith(discountCurve_);
50
0
    }
51
52
0
    void IsdaCdsEngine::calculate() const {
53
54
0
        QL_REQUIRE(numericalFix_ == None || numericalFix_ == Taylor,
55
0
                   "numerical fix must be None or Taylor");
56
0
        QL_REQUIRE(accrualBias_ == HalfDayBias || accrualBias_ == NoBias,
57
0
                   "accrual bias must be HalfDayBias or NoBias");
58
0
        QL_REQUIRE(forwardsInCouponPeriod_ == Flat ||
59
0
                       forwardsInCouponPeriod_ == Piecewise,
60
0
                   "forwards in coupon period must be Flat or Piecewise");
61
62
        // it would be possible to handle the cases which are excluded below,
63
        // but the ISDA engine is not explicitly specified to handle them,
64
        // so we just forbid them too
65
66
0
        Actual365Fixed dc;
67
0
        Actual360 dc1;
68
0
        Actual360 dc2(true);
69
70
0
        Date evalDate = Settings::instance().evaluationDate();
71
72
        // check if given curves are ISDA compatible
73
        // (the interpolation is checked below)
74
75
0
        QL_REQUIRE(!discountCurve_.empty(), "no discount term structure set");
76
0
        QL_REQUIRE(!probability_.empty(), "no probability term structure set");
77
0
        QL_REQUIRE(discountCurve_->dayCounter() == dc,
78
0
                   "yield term structure day counter ("
79
0
                       << discountCurve_->dayCounter()
80
0
                       << ") should be Act/365(Fixed)");
81
0
        QL_REQUIRE(probability_->dayCounter() == dc,
82
0
                   "probability term structure day counter ("
83
0
                       << probability_->dayCounter() << ") should be "
84
0
                       << "Act/365(Fixed)");
85
0
        QL_REQUIRE(discountCurve_->referenceDate() == evalDate,
86
0
                   "yield term structure reference date ("
87
0
                       << discountCurve_->referenceDate()
88
0
                       << " should be evaluation date (" << evalDate << ")");
89
0
        QL_REQUIRE(probability_->referenceDate() == evalDate,
90
0
                   "probability term structure reference date ("
91
0
                       << probability_->referenceDate()
92
0
                       << " should be evaluation date (" << evalDate << ")");
93
0
        QL_REQUIRE(arguments_.settlesAccrual,
94
0
                   "ISDA engine not compatible with non accrual paying CDS");
95
0
        QL_REQUIRE(arguments_.paysAtDefaultTime,
96
0
                   "ISDA engine not compatible with end period payment");
97
0
        QL_REQUIRE(ext::dynamic_pointer_cast<FaceValueClaim>(arguments_.claim) != nullptr,
98
0
                   "ISDA engine not compatible with non face value claim");
99
100
0
        Date maturity = arguments_.maturity;
101
0
        Date effectiveProtectionStart =
102
0
            std::max<Date>(arguments_.protectionStart, evalDate + 1);
103
104
        // collect nodes from both curves and sort them
105
0
        std::vector<Date> yDates, cDates;
106
107
        // the calls to dates() below might not trigger bootstrap (because
108
        // they will call the InterpolatedCurve methods, not the ones from
109
        // PiecewiseYieldCurve or PiecewiseDefaultCurve) so we force it here
110
0
        discountCurve_->discount(0.0);
111
0
        probability_->defaultProbability(0.0);
112
113
0
        if(ext::shared_ptr<InterpolatedDiscountCurve<LogLinear> > castY1 =
114
0
            ext::dynamic_pointer_cast<
115
0
                InterpolatedDiscountCurve<LogLinear> >(*discountCurve_)) {
116
0
            yDates = castY1->dates();
117
0
        } else if(ext::shared_ptr<InterpolatedForwardCurve<BackwardFlat> >
118
0
        castY2 = ext::dynamic_pointer_cast<
119
0
            InterpolatedForwardCurve<BackwardFlat> >(*discountCurve_)) {
120
0
            yDates = castY2->dates();
121
0
        } else if(ext::shared_ptr<InterpolatedForwardCurve<ForwardFlat> >
122
0
        castY3 = ext::dynamic_pointer_cast<
123
0
            InterpolatedForwardCurve<ForwardFlat> >(*discountCurve_)) {
124
0
            yDates = castY3->dates();
125
0
        } else if(ext::shared_ptr<FlatForward> castY4 =
126
0
            ext::dynamic_pointer_cast<FlatForward>(*discountCurve_)) {
127
            // no dates to extract
128
0
        } else {
129
0
            QL_FAIL("Yield curve must be flat forward interpolated");
130
0
        }
131
132
0
        if(ext::shared_ptr<InterpolatedSurvivalProbabilityCurve<LogLinear> >
133
0
        castC1 = ext::dynamic_pointer_cast<
134
0
            InterpolatedSurvivalProbabilityCurve<LogLinear> >(
135
0
            *probability_)) {
136
0
            cDates = castC1->dates();
137
0
        } else if(
138
0
        ext::shared_ptr<InterpolatedHazardRateCurve<BackwardFlat> > castC2 =
139
0
            ext::dynamic_pointer_cast<
140
0
            InterpolatedHazardRateCurve<BackwardFlat> >(*probability_)) {
141
0
            cDates = castC2->dates();
142
0
        } else if(
143
0
        ext::shared_ptr<FlatHazardRate> castC3 =
144
0
            ext::dynamic_pointer_cast<FlatHazardRate>(*probability_)) {
145
            // no dates to extract
146
0
        } else{
147
0
            QL_FAIL("Credit curve must be flat forward interpolated");
148
0
        }
149
150
0
        std::vector<Date> nodes;
151
0
        std::set_union(yDates.begin(), yDates.end(), cDates.begin(), cDates.end(), std::back_inserter(nodes));
152
153
154
0
        if(nodes.empty()){
155
0
            nodes.push_back(maturity);
156
0
        }
157
0
        const Real nFix = (numericalFix_ == None ? 1E-50 : 0.0);
158
159
        // protection leg pricing (npv is always negative at this stage)
160
0
        Real protectionNpv = 0.0;
161
162
0
        Date d0 = effectiveProtectionStart-1;
163
0
        Real P0 = discountCurve_->discount(d0);
164
0
        Real Q0 = probability_->survivalProbability(d0);
165
0
        Date d1;
166
0
        auto it =
167
0
            std::upper_bound(nodes.begin(), nodes.end(), effectiveProtectionStart);
168
169
0
        for(;it != nodes.end(); ++it) {
170
0
            if(*it > maturity) {
171
0
                d1 = maturity;
172
0
                it = nodes.end() - 1; //early exit
173
0
            } else {
174
0
                d1 = *it;
175
0
            }
176
0
            Real P1 = discountCurve_->discount(d1);
177
0
            Real Q1 = probability_->survivalProbability(d1);
178
179
0
            Real fhat = std::log(P0) - std::log(P1);
180
0
            Real hhat = std::log(Q0) - std::log(Q1);
181
0
            Real fhphh = fhat + hhat;
182
183
0
            if (fhphh < 1E-4 && numericalFix_ == Taylor) {
184
0
                Real fhphhq = fhphh * fhphh;
185
0
                protectionNpv +=
186
0
                    P0 * Q0 * hhat * (1.0 - 0.5 * fhphh + 1.0 / 6.0 * fhphhq -
187
0
                                      1.0 / 24.0 * fhphhq * fhphh +
188
0
                                      1.0 / 120 * fhphhq * fhphhq);
189
0
            } else {
190
0
                protectionNpv += hhat / (fhphh + nFix) * (P0 * Q0 - P1 * Q1);
191
0
            }
192
0
            d0 = d1;
193
0
            P0 = P1;
194
0
            Q0 = Q1;
195
0
        }
196
0
        protectionNpv *= arguments_.claim->amount(
197
0
            Date(), arguments_.notional, recoveryRate_);
198
199
0
        results_.defaultLegNPV = protectionNpv;
200
201
        // premium leg pricing (npv is always positive at this stage)
202
203
0
        Real premiumNpv = 0.0, defaultAccrualNpv = 0.0;
204
0
        for (auto& i : arguments_.leg) {
205
0
            ext::shared_ptr<FixedRateCoupon> coupon = ext::dynamic_pointer_cast<FixedRateCoupon>(i);
206
207
0
            QL_REQUIRE(coupon->dayCounter() == dc ||
208
0
                           coupon->dayCounter() == dc1 ||
209
0
                           coupon->dayCounter() == dc2,
210
0
                       "ISDA engine requires a coupon day counter Act/365Fixed "
211
0
                           << "or Act/360 (" << coupon->dayCounter() << ")");
212
213
            // premium coupons
214
0
            if (!i->hasOccurred(effectiveProtectionStart, includeSettlementDateFlows_)) {
215
0
                premiumNpv +=
216
0
                    coupon->amount() *
217
0
                    discountCurve_->discount(coupon->date()) *
218
0
                    probability_->survivalProbability(coupon->date()-1);
219
0
            }
220
221
            // default accruals
222
223
0
            if (!detail::simple_event(coupon->accrualEndDate())
224
0
                     .hasOccurred(effectiveProtectionStart, false)) {
225
0
                Date start = std::max<Date>(coupon->accrualStartDate(),
226
0
                                            effectiveProtectionStart)-1;
227
0
                Date end = coupon->date()-1;
228
0
                Real tstart =
229
0
                    discountCurve_->timeFromReference(coupon->accrualStartDate()-1) -
230
0
                    (accrualBias_ == HalfDayBias ? 1.0 / 730.0 : 0.0);
231
0
                std::vector<Date> localNodes;
232
0
                localNodes.push_back(start);
233
                //add intermediary nodes, if any
234
0
                if (forwardsInCouponPeriod_ == Piecewise) {
235
0
                    auto it0 =
236
0
                        std::upper_bound(nodes.begin(), nodes.end(), start);
237
0
                    auto it1 =
238
0
                        std::lower_bound(nodes.begin(), nodes.end(), end);
239
0
                    localNodes.insert(localNodes.end(), it0, it1);
240
0
                }
241
0
                localNodes.push_back(end);
242
243
0
                Real defaultAccrThisNode = 0.;
244
0
                auto node = localNodes.begin();
245
0
                Real t0 = discountCurve_->timeFromReference(*node);
246
0
                Real P0 = discountCurve_->discount(*node);
247
0
                Real Q0 = probability_->survivalProbability(*node);
248
249
0
                for (++node; node != localNodes.end(); ++node) {
250
0
                    Real t1 = discountCurve_->timeFromReference(*node);
251
0
                    Real P1 = discountCurve_->discount(*node);
252
0
                    Real Q1 = probability_->survivalProbability(*node);
253
0
                    Real fhat = std::log(P0) - std::log(P1);
254
0
                    Real hhat = std::log(Q0) - std::log(Q1);
255
0
                    Real fhphh = fhat + hhat;
256
0
                    if (fhphh < 1E-4 && numericalFix_ == Taylor) {
257
                        // see above, terms up to (f+h)^3 seem more than enough,
258
                        // what exactly is implemented in the standard isda C
259
                        // code ?
260
0
                        Real fhphhq = fhphh * fhphh;
261
0
                        defaultAccrThisNode +=
262
0
                            hhat * P0 * Q0 *
263
0
                            ((t0 - tstart) *
264
0
                                 (1.0 - 0.5 * fhphh + 1.0 / 6.0 * fhphhq -
265
0
                                  1.0 / 24.0 * fhphhq * fhphh) +
266
0
                             (t1 - t0) *
267
0
                                 (0.5 - 1.0 / 3.0 * fhphh + 1.0 / 8.0 * fhphhq -
268
0
                                  1.0 / 30.0 * fhphhq * fhphh));
269
0
                    } else {
270
0
                        defaultAccrThisNode +=
271
0
                            (hhat / (fhphh + nFix)) *
272
0
                            ((t1 - t0) * ((P0 * Q0 - P1 * Q1) / (fhphh + nFix) -
273
0
                                          P1 * Q1) +
274
0
                             (t0 - tstart) * (P0 * Q0 - P1 * Q1));
275
0
                    }
276
277
0
                    t0 = t1;
278
0
                    P0 = P1;
279
0
                    Q0 = Q1;
280
0
                }
281
0
                defaultAccrualNpv += defaultAccrThisNode * arguments_.notional *
282
0
                    coupon->rate() * 365. / 360.;
283
0
      }
284
0
        }
285
286
287
0
        results_.couponLegNPV = premiumNpv + defaultAccrualNpv;
288
289
        // upfront flow npv
290
291
0
        Real upfPVO1 = 0.0;
292
0
        results_.upfrontNPV = 0.0;
293
0
        if (!arguments_.upfrontPayment->hasOccurred(
294
0
                evalDate, includeSettlementDateFlows_)) {
295
0
            upfPVO1 =
296
0
                discountCurve_->discount(arguments_.upfrontPayment->date());
297
0
            if(arguments_.upfrontPayment->amount() != 0.) {
298
0
                results_.upfrontNPV = upfPVO1 * arguments_.upfrontPayment->amount();
299
0
            }
300
0
        }
301
302
0
        results_.accrualRebateNPV = 0.;
303
        // NOLINTNEXTLINE(readability-implicit-bool-conversion)
304
0
        if (arguments_.accrualRebate && arguments_.accrualRebate->amount() != 0. &&
305
0
            !arguments_.accrualRebate->hasOccurred(evalDate, includeSettlementDateFlows_)) {
306
0
            results_.accrualRebateNPV =
307
0
                discountCurve_->discount(arguments_.accrualRebate->date()) *
308
0
                arguments_.accrualRebate->amount();
309
0
        }
310
311
0
        Real upfrontSign = 1.0;
312
0
        switch (arguments_.side) {
313
0
          case Protection::Seller:
314
0
            results_.defaultLegNPV *= -1.0;
315
0
            results_.accrualRebateNPV *= -1.0;
316
0
            break;
317
0
          case Protection::Buyer:
318
0
            results_.couponLegNPV *= -1.0;
319
0
            results_.upfrontNPV   *= -1.0;
320
0
            upfrontSign = -1.0;
321
0
            break;
322
0
          default:
323
0
            QL_FAIL("unknown protection side");
324
0
        }
325
326
0
        results_.value = results_.defaultLegNPV + results_.couponLegNPV +
327
0
                         results_.upfrontNPV + results_.accrualRebateNPV;
328
329
0
        results_.errorEstimate = Null<Real>();
330
331
0
        if (results_.couponLegNPV != 0.0) {
332
0
            results_.fairSpread =
333
0
                -results_.defaultLegNPV * arguments_.spread /
334
0
                (results_.couponLegNPV + results_.accrualRebateNPV);
335
0
        } else {
336
0
            results_.fairSpread = Null<Rate>();
337
0
        }
338
339
0
        Real upfrontSensitivity = upfPVO1 * arguments_.notional;
340
0
        if (upfrontSensitivity != 0.0) {
341
0
            results_.fairUpfront =
342
0
                -upfrontSign * (results_.defaultLegNPV + results_.couponLegNPV +
343
0
                                results_.accrualRebateNPV) /
344
0
                upfrontSensitivity;
345
0
        } else {
346
0
            results_.fairUpfront = Null<Rate>();
347
0
        }
348
349
0
        static const Rate basisPoint = 1.0e-4;
350
351
0
        if (arguments_.spread != 0.0) {
352
0
            results_.couponLegBPS =
353
0
                results_.couponLegNPV * basisPoint / arguments_.spread;
354
0
        } else {
355
0
            results_.couponLegBPS = Null<Rate>();
356
0
        }
357
358
        // NOLINTNEXTLINE(readability-implicit-bool-conversion)
359
0
        if (arguments_.upfront && *arguments_.upfront != 0.0) {
360
0
            results_.upfrontBPS =
361
0
                results_.upfrontNPV * basisPoint / (*arguments_.upfront);
362
0
        } else {
363
0
            results_.upfrontBPS = Null<Rate>();
364
0
        }
365
0
    }
366
}