Coverage Report

Created: 2026-01-25 06:59

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/quantlib/ql/pricingengines/swaption/gaussian1dnonstandardswaptionengine.cpp
Line
Count
Source
1
/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3
/*
4
 Copyright (C) 2013 Peter Caspers
5
6
 This file is part of QuantLib, a free-software/open-source library
7
 for financial quantitative analysts and developers - http://quantlib.org/
8
9
 QuantLib is free software: you can redistribute it and/or modify it
10
 under the terms of the QuantLib license.  You should have received a
11
 copy of the license along with this program; if not, please email
12
 <quantlib-dev@lists.sf.net>. The license is also available online at
13
 <https://www.quantlib.org/license.shtml>.
14
15
 This program is distributed in the hope that it will be useful, but WITHOUT
16
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17
 FOR A PARTICULAR PURPOSE.  See the license for more details.
18
*/
19
20
#include <ql/pricingengines/swaption/gaussian1dnonstandardswaptionengine.hpp>
21
#include <ql/rebatedexercise.hpp>
22
#include <ql/time/daycounters/actualactual.hpp>
23
#include <ql/quotes/simplequote.hpp>
24
#include <ql/math/interpolations/cubicinterpolation.hpp>
25
#include <ql/payoff.hpp>
26
27
using std::exp;
28
29
namespace QuantLib {
30
31
    Real
32
    Gaussian1dNonstandardSwaptionEngine::underlyingNpv(const Date &expiry,
33
0
                                                       const Real y) const {
34
35
        // determine the indices on both legs representing the cashflows that
36
        // are part of the exercise into right
37
38
0
        Size fixedIdx =
39
0
            std::upper_bound(arguments_.fixedResetDates.begin(),
40
0
                             arguments_.fixedResetDates.end(), expiry - 1) -
41
0
            arguments_.fixedResetDates.begin();
42
0
        Size floatingIdx =
43
0
            std::upper_bound(arguments_.floatingResetDates.begin(),
44
0
                             arguments_.floatingResetDates.end(), expiry - 1) -
45
0
            arguments_.floatingResetDates.begin();
46
47
        // calculate the npv of these cashflows conditional on y at expiry
48
49
0
        Real type = (Real)arguments_.type;
50
51
0
        Real npv = 0.0;
52
0
        for (Size i = fixedIdx; i < arguments_.fixedResetDates.size(); i++) {
53
0
            npv -=
54
0
                arguments_.fixedCoupons[i] *
55
0
                model_->zerobond(arguments_.fixedPayDates[i], expiry, y,
56
0
                                 discountCurve_) *
57
0
                (oas_.empty()
58
0
                     ? Real(1.0)
59
0
                     : exp(-oas_->value() *
60
0
                           model_->termStructure()->dayCounter().yearFraction(
61
0
                               expiry, arguments_.fixedPayDates[i])));
62
0
        }
63
64
0
        for (Size i = floatingIdx; i < arguments_.floatingResetDates.size();
65
0
             i++) {
66
0
            Real amount;
67
0
            if (!arguments_.floatingIsRedemptionFlow[i])
68
0
                amount = (arguments_.floatingGearings[i] *
69
0
                              model_->forwardRate(
70
0
                                  arguments_.floatingFixingDates[i], expiry, y,
71
0
                                  arguments_.swap->iborIndex()) +
72
0
                          arguments_.floatingSpreads[i]) *
73
0
                         arguments_.floatingAccrualTimes[i] *
74
0
                         arguments_.floatingNominal[i];
75
0
            else
76
0
                amount = arguments_.floatingCoupons[i];
77
0
            npv +=
78
0
                amount * model_->zerobond(arguments_.floatingPayDates[i],
79
0
                                          expiry, y, discountCurve_) *
80
0
                (oas_.empty()
81
0
                     ? Real(1.0)
82
0
                     : exp(-oas_->value() *
83
0
                           model_->termStructure()->dayCounter().yearFraction(
84
0
                               expiry, arguments_.floatingPayDates[i])));
85
0
        }
86
87
0
        return type * npv;
88
0
    }
89
90
0
    Swap::Type Gaussian1dNonstandardSwaptionEngine::underlyingType() const {
91
0
        return arguments_.swap->type();
92
0
    }
93
94
    // NOLINTNEXTLINE(readability-const-return-type)
95
0
    const Date Gaussian1dNonstandardSwaptionEngine::underlyingLastDate() const {
96
0
        return arguments_.fixedPayDates.back();
97
0
    }
98
99
    // NOLINTNEXTLINE(readability-const-return-type)
100
0
    const Array Gaussian1dNonstandardSwaptionEngine::initialGuess(const Date &expiry) const {
101
102
0
        Size fixedIdx =
103
0
            std::upper_bound(arguments_.fixedResetDates.begin(),
104
0
                             arguments_.fixedResetDates.end(), expiry - 1) -
105
0
            arguments_.fixedResetDates.begin();
106
107
0
        Array initial(3);
108
0
        Real nominalSum = 0.0, weightedRate = 0.0, ind = 0.0;
109
0
        for (Size i = fixedIdx; i < arguments_.fixedResetDates.size(); i++) {
110
0
            nominalSum += arguments_.fixedNominal[i];
111
0
            Real rate = arguments_.fixedRate[i];
112
0
            if (close(rate, 0.0))
113
0
                rate = 0.03; // this value is at least better than zero
114
0
            weightedRate += arguments_.fixedNominal[i] * rate;
115
0
            if (arguments_.fixedNominal[i] > 1E-8) // exclude zero nominal periods
116
0
                ind += 1.0;
117
0
        }
118
0
        Real nominalAvg = nominalSum / ind;
119
120
0
        QL_REQUIRE(nominalSum > 0.0,
121
0
                   "sum of nominals on fixed leg must be positive ("
122
0
                       << nominalSum << ")");
123
124
0
        weightedRate /= nominalSum;
125
0
        initial[0] = nominalAvg;
126
0
        initial[1] =
127
0
            model_->termStructure()->timeFromReference(underlyingLastDate()) -
128
0
            model_->termStructure()->timeFromReference(expiry);
129
0
        initial[2] = weightedRate;
130
131
0
        return initial;
132
0
    }
133
134
0
    void Gaussian1dNonstandardSwaptionEngine::calculate() const {
135
136
0
        QL_REQUIRE(arguments_.settlementMethod != Settlement::ParYieldCurve,
137
0
                   "cash settled (ParYieldCurve) swaptions not priced with "
138
0
                   "Gaussian1dNonstandardSwaptionEngine");
139
140
0
        Date settlement = model_->termStructure()->referenceDate();
141
142
0
        if (arguments_.exercise->dates().back() <=
143
0
            settlement) { // swaption is expired, possibly generated swap is not
144
                          // valued
145
0
            results_.value = 0.0;
146
0
            return;
147
0
        }
148
149
0
        ext::shared_ptr<RebatedExercise> rebatedExercise =
150
0
            ext::dynamic_pointer_cast<RebatedExercise>(arguments_.exercise);
151
152
0
        int idx = arguments_.exercise->dates().size() - 1;
153
0
        int minIdxAlive = static_cast<int>(
154
0
            std::upper_bound(arguments_.exercise->dates().begin(),
155
0
                             arguments_.exercise->dates().end(), settlement) -
156
0
            arguments_.exercise->dates().begin());
157
158
0
        NonstandardSwap swap = *arguments_.swap;
159
0
        Option::Type type =
160
0
            arguments_.type == Swap::Payer ? Option::Call : Option::Put;
161
162
0
        Array npv0(2 * integrationPoints_ + 1, 0.0),
163
0
            npv1(2 * integrationPoints_ + 1, 0.0);
164
0
        Array z = model_->yGrid(stddevs_, integrationPoints_);
165
0
        Array p(z.size(), 0.0);
166
167
        // for probability computation
168
0
        std::vector<Array> npvp0, npvp1;
169
0
        if (probabilities_ != None) {
170
0
            for (int i = 0; i < idx - minIdxAlive + 2; ++i) {
171
0
                Array npvTmp0(2 * integrationPoints_ + 1, 0.0);
172
0
                Array npvTmp1(2 * integrationPoints_ + 1, 0.0);
173
0
                npvp0.push_back(npvTmp0);
174
0
                npvp1.push_back(npvTmp1);
175
0
            }
176
0
        }
177
        // end probability computation
178
179
0
        Date expiry1 = Date(), expiry0;
180
0
        Time expiry1Time = Null<Real>(), expiry0Time;
181
182
0
        do {
183
184
0
            if (idx == minIdxAlive - 1)
185
0
                expiry0 = settlement;
186
0
            else
187
0
                expiry0 = arguments_.exercise->dates()[idx];
188
189
0
            expiry0Time = std::max(
190
0
                model_->termStructure()->timeFromReference(expiry0), 0.0);
191
192
0
            Size j1 =
193
0
                std::upper_bound(arguments_.fixedResetDates.begin(),
194
0
                                 arguments_.fixedResetDates.end(), expiry0 - 1) -
195
0
                arguments_.fixedResetDates.begin();
196
0
            Size k1 =
197
0
                std::upper_bound(arguments_.floatingResetDates.begin(),
198
0
                                 arguments_.floatingResetDates.end(), expiry0 - 1) -
199
0
                arguments_.floatingResetDates.begin();
200
201
            // todo add openmp support later on (as in gaussian1dswaptionengine)
202
203
0
            for (Size k = 0; k < (expiry0 > settlement ? npv0.size() : 1);
204
0
                 k++) {
205
206
0
                Real price = 0.0;
207
0
                if (expiry1Time != Null<Real>()) {
208
0
                    Real zSpreadDf =
209
0
                        oas_.empty() ? Real(1.0)
210
0
                                     : std::exp(-oas_->value() *
211
0
                                                (expiry1Time - expiry0Time));
212
0
                    Array yg = model_->yGrid(stddevs_, integrationPoints_,
213
0
                                             expiry1Time, expiry0Time,
214
0
                                             expiry0 > settlement ? z[k] : 0.0);
215
0
                    CubicInterpolation payoff0(
216
0
                        z.begin(), z.end(), npv1.begin(),
217
0
                        CubicInterpolation::Spline, true,
218
0
                        CubicInterpolation::Lagrange, 0.0,
219
0
                        CubicInterpolation::Lagrange, 0.0);
220
0
                    for (Size i = 0; i < yg.size(); i++) {
221
0
                        p[i] = payoff0(yg[i], true);
222
0
                    }
223
0
                    CubicInterpolation payoff1(
224
0
                        z.begin(), z.end(), p.begin(),
225
0
                        CubicInterpolation::Spline, true,
226
0
                        CubicInterpolation::Lagrange, 0.0,
227
0
                        CubicInterpolation::Lagrange, 0.0);
228
0
                    for (Size i = 0; i < z.size() - 1; i++) {
229
0
                        price += Gaussian1dModel::gaussianShiftedPolynomialIntegral(
230
0
                                     0.0, payoff1.cCoefficients()[i],
231
0
                                     payoff1.bCoefficients()[i],
232
0
                                     payoff1.aCoefficients()[i], p[i], z[i],
233
0
                                     z[i], z[i + 1]) *
234
0
                                 zSpreadDf;
235
0
                    }
236
0
                    if (extrapolatePayoff_) {
237
0
                        if (flatPayoffExtrapolation_) {
238
0
                            price +=
239
0
                                Gaussian1dModel::gaussianShiftedPolynomialIntegral(
240
0
                                    0.0, 0.0, 0.0, 0.0, p[z.size() - 2],
241
0
                                    z[z.size() - 2], z[z.size() - 1], 100.0) *
242
0
                                zSpreadDf;
243
0
                            price += Gaussian1dModel::gaussianShiftedPolynomialIntegral(
244
0
                                         0.0, 0.0, 0.0, 0.0, p[0], z[0], -100.0,
245
0
                                         z[0]) *
246
0
                                     zSpreadDf;
247
0
                        } else {
248
0
                            if (type == Option::Call)
249
0
                                price +=
250
0
                                    Gaussian1dModel::gaussianShiftedPolynomialIntegral(
251
0
                                        0.0,
252
0
                                        payoff1.cCoefficients()[z.size() - 2],
253
0
                                        payoff1.bCoefficients()[z.size() - 2],
254
0
                                        payoff1.aCoefficients()[z.size() - 2],
255
0
                                        p[z.size() - 2], z[z.size() - 2],
256
0
                                        z[z.size() - 1], 100.0) *
257
0
                                    zSpreadDf;
258
0
                            if (type == Option::Put)
259
0
                                price +=
260
0
                                    Gaussian1dModel::gaussianShiftedPolynomialIntegral(
261
0
                                        0.0, payoff1.cCoefficients()[0],
262
0
                                        payoff1.bCoefficients()[0],
263
0
                                        payoff1.aCoefficients()[0], p[0], z[0],
264
0
                                        -100.0, z[0]) *
265
0
                                    zSpreadDf;
266
0
                        }
267
0
                    }
268
0
                }
269
270
0
                npv0[k] = price;
271
272
                // for probability computation
273
0
                if (probabilities_ != None) {
274
0
                    for (Size m = 0; m < npvp0.size(); m++) {
275
0
                        Real price = 0.0;
276
0
                        if (expiry1Time != Null<Real>()) {
277
0
                            Real zSpreadDf =
278
0
                                oas_.empty()
279
0
                                    ? Real(1.0)
280
0
                                    : std::exp(-oas_->value() *
281
0
                                               (expiry1Time - expiry0Time));
282
0
                            Array yg = model_->yGrid(
283
0
                                stddevs_, integrationPoints_, expiry1Time,
284
0
                                expiry0Time, expiry0 > settlement ? z[k] : 0.0);
285
0
                            CubicInterpolation payoff0(
286
0
                                z.begin(), z.end(), npvp1[m].begin(),
287
0
                                CubicInterpolation::Spline, true,
288
0
                                CubicInterpolation::Lagrange, 0.0,
289
0
                                CubicInterpolation::Lagrange, 0.0);
290
0
                            for (Size i = 0; i < yg.size(); i++) {
291
0
                                p[i] = payoff0(yg[i], true);
292
0
                            }
293
0
                            CubicInterpolation payoff1(
294
0
                                z.begin(), z.end(), p.begin(),
295
0
                                CubicInterpolation::Spline, true,
296
0
                                CubicInterpolation::Lagrange, 0.0,
297
0
                                CubicInterpolation::Lagrange, 0.0);
298
0
                            for (Size i = 0; i < z.size() - 1; i++) {
299
0
                                price +=
300
0
                                    Gaussian1dModel::gaussianShiftedPolynomialIntegral(
301
0
                                        0.0, payoff1.cCoefficients()[i],
302
0
                                        payoff1.bCoefficients()[i],
303
0
                                        payoff1.aCoefficients()[i], p[i], z[i],
304
0
                                        z[i], z[i + 1]) *
305
0
                                    zSpreadDf;
306
0
                            }
307
0
                            if (extrapolatePayoff_) {
308
0
                                if (flatPayoffExtrapolation_) {
309
0
                                    price +=
310
0
                                        Gaussian1dModel::gaussianShiftedPolynomialIntegral(
311
0
                                                  0.0, 0.0, 0.0, 0.0,
312
0
                                                  p[z.size() - 2],
313
0
                                                  z[z.size() - 2],
314
0
                                                  z[z.size() - 1], 100.0) *
315
0
                                        zSpreadDf;
316
0
                                    price +=
317
0
                                        Gaussian1dModel::gaussianShiftedPolynomialIntegral(
318
0
                                                  0.0, 0.0, 0.0, 0.0, p[0],
319
0
                                                  z[0], -100.0, z[0]) *
320
0
                                        zSpreadDf;
321
0
                                } else {
322
0
                                    if (type == Option::Call)
323
0
                                        price +=
324
0
                                            Gaussian1dModel::gaussianShiftedPolynomialIntegral(
325
0
                                                      0.0,
326
0
                                                      payoff1.cCoefficients()
327
0
                                                          [z.size() - 2],
328
0
                                                      payoff1.bCoefficients()
329
0
                                                          [z.size() - 2],
330
0
                                                      payoff1.aCoefficients()
331
0
                                                          [z.size() - 2],
332
0
                                                      p[z.size() - 2],
333
0
                                                      z[z.size() - 2],
334
0
                                                      z[z.size() - 1], 100.0) *
335
0
                                            zSpreadDf;
336
0
                                    if (type == Option::Put)
337
0
                                        price +=
338
0
                                            Gaussian1dModel::gaussianShiftedPolynomialIntegral(
339
0
                                                      0.0,
340
0
                                                      payoff1
341
0
                                                          .cCoefficients()[0],
342
0
                                                      payoff1
343
0
                                                          .bCoefficients()[0],
344
0
                                                      payoff1
345
0
                                                          .aCoefficients()[0],
346
0
                                                      p[0], z[0], -100.0,
347
0
                                                      z[0]) *
348
0
                                            zSpreadDf;
349
0
                                }
350
0
                            }
351
0
                        }
352
353
0
                        npvp0[m][k] = price;
354
0
                    }
355
0
                }
356
                // end probability computation
357
358
0
                if (expiry0 > settlement) {
359
0
                    Real floatingLegNpv = 0.0;
360
0
                    for (Size l = k1; l < arguments_.floatingCoupons.size();
361
0
                         l++) {
362
0
                        Real zSpreadDf =
363
0
                            oas_.empty()
364
0
                                ? Real(1.0)
365
0
                                : std::exp(
366
0
                                      -oas_->value() *
367
0
                                      (model_->termStructure()
368
0
                                           ->dayCounter()
369
0
                                           .yearFraction(
370
0
                                                expiry0,
371
0
                                                arguments_
372
0
                                                    .floatingPayDates[l])));
373
0
                        Real amount;
374
0
                        if (arguments_.floatingIsRedemptionFlow[l])
375
0
                            amount = arguments_.floatingCoupons[l];
376
0
                        else
377
0
                            amount = arguments_.floatingNominal[l] *
378
0
                                     arguments_.floatingAccrualTimes[l] *
379
0
                                     (arguments_.floatingGearings[l] *
380
0
                                          model_->forwardRate(
381
0
                                              arguments_.floatingFixingDates[l],
382
0
                                              expiry0, z[k],
383
0
                                              arguments_.swap->iborIndex()) +
384
0
                                      arguments_.floatingSpreads[l]);
385
0
                        floatingLegNpv +=
386
0
                            amount *
387
0
                            model_->zerobond(arguments_.floatingPayDates[l],
388
0
                                             expiry0, z[k], discountCurve_) *
389
0
                            zSpreadDf;
390
0
                    }
391
0
                    Real fixedLegNpv = 0.0;
392
0
                    for (Size l = j1; l < arguments_.fixedCoupons.size(); l++) {
393
0
                        Real zSpreadDf =
394
0
                            oas_.empty()
395
0
                                ? Real(1.0)
396
0
                                : std::exp(
397
0
                                      -oas_->value() *
398
0
                                      (model_->termStructure()
399
0
                                           ->dayCounter()
400
0
                                           .yearFraction(
401
0
                                                expiry0,
402
0
                                                arguments_.fixedPayDates[l])));
403
0
                        fixedLegNpv +=
404
0
                            arguments_.fixedCoupons[l] *
405
0
                            model_->zerobond(arguments_.fixedPayDates[l],
406
0
                                             expiry0, z[k], discountCurve_) *
407
0
                            zSpreadDf;
408
0
                    }
409
0
                    Real rebate = 0.0;
410
0
                    Real zSpreadDf = 1.0;
411
0
                    Date rebateDate = expiry0;
412
0
                    if (rebatedExercise != nullptr) {
413
0
                        rebate = rebatedExercise->rebate(idx);
414
0
                        rebateDate = rebatedExercise->rebatePaymentDate(idx);
415
0
                        zSpreadDf =
416
0
                            oas_.empty()
417
0
                                ? Real(1.0)
418
0
                                : std::exp(
419
0
                                      -oas_->value() *
420
0
                                      (model_->termStructure()
421
0
                                           ->dayCounter()
422
0
                                           .yearFraction(expiry0, rebateDate)));
423
0
                    }
424
0
                    Real exerciseValue =
425
0
                        ((type == Option::Call ? 1.0 : -1.0) *
426
0
                             (floatingLegNpv - fixedLegNpv) +
427
0
                         rebate * model_->zerobond(rebateDate, expiry0, z[k],
428
0
                                                   discountCurve_) *
429
0
                             zSpreadDf) /
430
0
                        model_->numeraire(expiry0Time, z[k], discountCurve_);
431
432
                    // for probability computation
433
0
                    if (probabilities_ != None) {
434
0
                        if (idx == static_cast<int>(
435
0
                                       arguments_.exercise->dates().size()) -
436
0
                                       1) // if true we are at the latest date,
437
                                          // so we init
438
                                          // the no call probability
439
0
                            npvp0.back()[k] =
440
0
                                probabilities_ == Naive
441
0
                                    ? Real(1.0)
442
0
                                    : 1.0 / (model_->zerobond(expiry0Time, 0.0,
443
0
                                                              0.0,
444
0
                                                              discountCurve_) *
445
0
                                             model_->numeraire(expiry0, z[k],
446
0
                                                               discountCurve_));
447
0
                        if (exerciseValue >= npv0[k]) {
448
0
                            npvp0[idx - minIdxAlive][k] =
449
0
                                probabilities_ == Naive
450
0
                                    ? Real(1.0)
451
0
                                    : 1.0 /
452
0
                                          (model_->zerobond(expiry0Time, 0.0,
453
0
                                                            0.0,
454
0
                                                            discountCurve_) *
455
0
                                           model_->numeraire(expiry0Time, z[k],
456
0
                                                             discountCurve_));
457
0
                            for (Size ii = idx - minIdxAlive + 1;
458
0
                                 ii < npvp0.size(); ii++)
459
0
                                npvp0[ii][k] = 0.0;
460
0
                        }
461
0
                    }
462
                    // end probability computation
463
464
0
                    npv0[k] = std::max(npv0[k], exerciseValue);
465
0
                }
466
0
            }
467
468
0
            npv1.swap(npv0);
469
470
            // for probability computation
471
0
            if (probabilities_ != None) {
472
0
                for (Size i = 0; i < npvp0.size(); i++)
473
0
                    npvp1[i].swap(npvp0[i]);
474
0
            }
475
            // end probability computation
476
477
0
            expiry1 = expiry0;
478
0
            expiry1Time = expiry0Time;
479
480
0
        } while (--idx >= minIdxAlive - 1);
481
482
0
        results_.value = npv1[0] * model_->numeraire(0.0, 0.0, discountCurve_);
483
484
        // for probability computation
485
0
        if (probabilities_ != None) {
486
0
            std::vector<Real> prob(npvp0.size());
487
0
            for (Size i = 0; i < npvp0.size(); i++) {
488
0
                prob[i] = npvp1[i][0] *
489
0
                          (probabilities_ == Naive
490
0
                               ? 1.0
491
0
                               : model_->numeraire(0.0, 0.0, discountCurve_));
492
0
            }
493
0
            results_.additionalResults["probabilities"] = prob;
494
0
        }
495
        // end probability computation
496
0
    }
497
}