QuantLib: a free/open-source library for quantitative finance
Fully annotated sources - version 1.32
Loading...
Searching...
No Matches
gaussian1dfloatfloatswaptionengine.cpp
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2013, 2015 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 <http://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/gaussian1dfloatfloatswaptionengine.hpp>
21#include <ql/experimental/coupons/swapspreadindex.hpp> // internal
22#include <ql/math/interpolations/cubicinterpolation.hpp>
23#include <ql/payoff.hpp>
24
25namespace QuantLib {
26
28
29 QL_REQUIRE(arguments_.settlementMethod != Settlement::ParYieldCurve,
30 "cash settled (ParYieldCurve) swaptions not priced with "
31 "Gaussian1dFloatFloatSwaptionEngine");
32
33 Date settlement = model_->termStructure()->referenceDate();
34
35 if (arguments_.exercise->dates().back() <=
36 settlement) { // swaption is expired,
37 // possibly generated swap
38 // is not
39 // valued
40 results_.value = 0.0;
41 return;
42 }
43
45 ext::dynamic_pointer_cast<RebatedExercise>(arguments_.exercise);
46
47 std::pair<Real, Real> result =
48 npvs(settlement, 0.0, includeTodaysExercise_, true);
49
50 results_.value = result.first;
51 results_.additionalResults["underlyingValue"] = result.second;
52 }
53
54 Real
56 const Real y) const {
57 return npvs(expiry, y, true).second;
58 }
59
61 return arguments_.swap->type();
62 }
63
64 // NOLINTNEXTLINE(readability-const-return-type)
66 Date l1 = arguments_.leg1PayDates.back();
67 Date l2 = arguments_.leg2PayDates.back();
68 return l2 >= l1 ? l2 : l1;
69 }
70
71 // NOLINTNEXTLINE(readability-const-return-type)
73
74 Size idx1 =
75 std::upper_bound(arguments_.leg1ResetDates.begin(),
76 arguments_.leg1ResetDates.end(), expiry - 1) -
77 arguments_.leg1ResetDates.begin();
78
79 // very simple initial guess
80 // check guess for nominal and weighted maturity !
81
82 Real nominalSum1 = 0.0;
83 for (Size i = idx1; i < arguments_.leg1ResetDates.size(); i++) {
84 nominalSum1 += arguments_.nominal1[i];
85 }
86 Real nominalAvg1 = nominalSum1 /
87 (arguments_.leg1ResetDates.size() - idx1);
88 Real weightedMaturity1 = 0.0;
89 for (Size i = idx1; i < arguments_.leg1ResetDates.size(); i++) {
90 weightedMaturity1 +=
91 arguments_.leg1AccrualTimes[i] * arguments_.nominal1[i];
92 }
93 weightedMaturity1 /= nominalAvg1;
94
95 return {
96 nominalAvg1,
97 weightedMaturity1,
98 0.03 // ???
99 };
100 }
101
102 // calculate npv and underlying npv as of expiry date
103 std::pair<Real, Real>
105 const Real y,
106 const bool includeExerciseOnExpiry,
107 const bool considerProbabilities) const {
108
109 // pricing
110
111 // event dates are coupon fixing dates and exercise dates
112 // we explicitly estimate cms and also libor coupons (although
113 // the latter could be calculated analytically) to make the code
114 // simpler
115
116 std::vector<Date> events;
117 events.insert(events.end(), arguments_.exercise->dates().begin(),
118 arguments_.exercise->dates().end());
119 events.insert(events.end(), arguments_.leg1FixingDates.begin(),
120 arguments_.leg1FixingDates.end());
121 events.insert(events.end(), arguments_.leg2FixingDates.begin(),
122 arguments_.leg2FixingDates.end());
123 std::sort(events.begin(), events.end());
124
125 auto it = std::unique(events.begin(), events.end());
126 events.resize(std::distance(events.begin(), it));
127
128 // only events on or after expiry are of interest by definition of the
129 // deal part that is exericsed into.
130
131 auto filit = std::upper_bound(events.begin(), events.end(),
132 expiry - (includeExerciseOnExpiry ? 1 : 0));
133 events.erase(events.begin(), filit);
134
135 int idx = events.size() - 1;
136
138 Option::Type type =
140
141 Array npv0(2 * integrationPoints_ + 1, 0.0),
142 npv1(2 * integrationPoints_ + 1, 0.0); // arrays for npvs of the
143 // option
144 Array npv0a(2 * integrationPoints_ + 1, 0.0),
145 npv1a(2 * integrationPoints_ + 1, 0.0); // arrays for npvs of the
146 // underlying
148 Array p(z.size(), 0.0), pa(z.size(), 0.0);
149
150 // for probability computation
151 std::vector<Array> npvp0, npvp1;
152 // how many active exercise dates are there ?
153 Size noEx = arguments_.exercise->dates().size() -
154 (std::upper_bound(arguments_.exercise->dates().begin(),
155 arguments_.exercise->dates().end(),
156 expiry - (includeExerciseOnExpiry ? 1 : 0)) -
157 arguments_.exercise->dates().begin());
158 Size exIdx = noEx; // current exercise index
159 if (considerProbabilities && probabilities_ != None) {
160 for (Size i = 0; i < noEx+1 ; ++i) {
161 Array npvTmp0(2 * integrationPoints_ + 1, 0.0);
162 Array npvTmp1(2 * integrationPoints_ + 1, 0.0);
163 npvp0.push_back(npvTmp0);
164 npvp1.push_back(npvTmp1);
165 }
166 }
167 // end probability computation
168
169 Date event1 = Null<Date>(), event0;
170 Time event1Time = Null<Real>(), event0Time;
171
172 ext::shared_ptr<IborIndex> ibor1 =
173 ext::dynamic_pointer_cast<IborIndex>(arguments_.index1);
174 ext::shared_ptr<SwapIndex> cms1 =
175 ext::dynamic_pointer_cast<SwapIndex>(arguments_.index1);
176 ext::shared_ptr<SwapSpreadIndex> cmsspread1 =
177 ext::dynamic_pointer_cast<SwapSpreadIndex>(arguments_.index1);
178 ext::shared_ptr<IborIndex> ibor2 =
179 ext::dynamic_pointer_cast<IborIndex>(arguments_.index2);
180 ext::shared_ptr<SwapIndex> cms2 =
181 ext::dynamic_pointer_cast<SwapIndex>(arguments_.index2);
182 ext::shared_ptr<SwapSpreadIndex> cmsspread2 =
183 ext::dynamic_pointer_cast<SwapSpreadIndex>(arguments_.index2);
184
185 QL_REQUIRE(ibor1 != nullptr || cms1 != nullptr || cmsspread1 != nullptr,
186 "index1 must be ibor or swap or swap spread index");
187 QL_REQUIRE(ibor2 != nullptr || cms2 != nullptr || cmsspread2 != nullptr,
188 "index2 must be ibor or swap or swap spread index");
189
190 do {
191
192 // we are at event0 date, which can be a structured coupon fixing
193 // date or an exercise date or both.
194
195 bool isEventDate = true;
196 if (idx == -1) {
197 event0 = expiry;
198 isEventDate = false;
199 } else {
200 event0 = events[idx];
201 if (event0 == expiry)
202 idx = -1; // avoid double roll back if expiry equal to
203 // earliest event date
204 }
205
206 bool isExercise =
207 std::find(arguments_.exercise->dates().begin(), arguments_.exercise->dates().end(),
208 event0) != arguments_.exercise->dates().end();
209
210 bool isLeg1Fixing =
211 std::find(arguments_.leg1FixingDates.begin(), arguments_.leg1FixingDates.end(),
212 event0) != arguments_.leg1FixingDates.end();
213
214 bool isLeg2Fixing =
215 std::find(arguments_.leg2FixingDates.begin(), arguments_.leg2FixingDates.end(),
216 event0) != arguments_.leg2FixingDates.end();
217
218 event0Time = std::max(
219 model_->termStructure()->timeFromReference(event0), 0.0);
220
221 // todo add openmp support later on (as in gaussian1dswaptionengine)
222
223 for (Size k = 0; k < (event0 > expiry ? npv0.size() : 1); k++) {
224
225 // roll back
226
227 Real price = 0.0, pricea = 0.0;
228 if (event1Time != Null<Real>()) {
229 Real zSpreadDf = oas_.empty()
230 ? Real(1.0)
231 : std::exp(-oas_->value() *
232 (event1Time - event0Time));
233 Array yg =
234 model_->yGrid(stddevs_, integrationPoints_, event1Time,
235 event0Time, event0 > expiry ? z[k] : y);
236 CubicInterpolation payoff0(
237 z.begin(), z.end(), npv1.begin(),
241 CubicInterpolation payoff0a(
242 z.begin(), z.end(), npv1a.begin(),
246 for (Size i = 0; i < yg.size(); i++) {
247 p[i] = payoff0(yg[i], true);
248 pa[i] = payoff0a(yg[i], true);
249 }
250 CubicInterpolation payoff1(
251 z.begin(), z.end(), p.begin(),
255 CubicInterpolation payoff1a(
256 z.begin(), z.end(), pa.begin(),
260 for (Size i = 0; i < z.size() - 1; i++) {
262 0.0, payoff1.cCoefficients()[i],
263 payoff1.bCoefficients()[i],
264 payoff1.aCoefficients()[i], p[i], z[i],
265 z[i], z[i + 1]) *
266 zSpreadDf;
268 0.0, payoff1a.cCoefficients()[i],
269 payoff1a.bCoefficients()[i],
270 payoff1a.aCoefficients()[i], pa[i], z[i],
271 z[i], z[i + 1]) *
272 zSpreadDf;
273 }
274 if (extrapolatePayoff_) {
276 price +=
278 0.0, 0.0, 0.0, 0.0, p[z.size() - 2],
279 z[z.size() - 2], z[z.size() - 1], 100.0) *
280 zSpreadDf;
282 0.0, 0.0, 0.0, 0.0, p[0], z[0], -100.0,
283 z[0]) *
284 zSpreadDf;
285 pricea +=
287 0.0, 0.0, 0.0, 0.0, pa[z.size() - 2],
288 z[z.size() - 2], z[z.size() - 1], 100.0) *
289 zSpreadDf;
291 0.0, 0.0, 0.0, 0.0, pa[0], z[0],
292 -100.0, z[0]) *
293 zSpreadDf;
294 } else {
295 if (type == Option::Call)
296 price +=
298 0.0,
299 payoff1.cCoefficients()[z.size() - 2],
300 payoff1.bCoefficients()[z.size() - 2],
301 payoff1.aCoefficients()[z.size() - 2],
302 p[z.size() - 2], z[z.size() - 2],
303 z[z.size() - 1], 100.0) *
304 zSpreadDf;
305 if (type == Option::Put)
306 price +=
308 0.0, payoff1.cCoefficients()[0],
309 payoff1.bCoefficients()[0],
310 payoff1.aCoefficients()[0], p[0], z[0],
311 -100.0, z[0]) *
312 zSpreadDf;
313 if (type == Option::Call)
314 pricea +=
316 0.0,
317 payoff1a.cCoefficients()[z.size() - 2],
318 payoff1a.bCoefficients()[z.size() - 2],
319 payoff1a.aCoefficients()[z.size() - 2],
320 pa[z.size() - 2], z[z.size() - 2],
321 z[z.size() - 1], 100.0) *
322 zSpreadDf;
323 if (type == Option::Put)
324 pricea +=
326 0.0, payoff1a.cCoefficients()[0],
327 payoff1a.bCoefficients()[0],
328 payoff1a.aCoefficients()[0], pa[0],
329 z[0], -100.0, z[0]) *
330 zSpreadDf;
331 }
332 }
333 }
334
335 npv0[k] = price;
336 npv0a[k] = pricea;
337
338 // for probability computation
339 if (considerProbabilities && probabilities_ != None) {
340 for (Size m = 0; m < npvp0.size(); m++) {
341 Real price = 0.0;
342 if (event1Time != Null<Real>()) {
343 Real zSpreadDf =
344 oas_.empty()
345 ? Real(1.0)
346 : std::exp(-oas_->value() *
347 (event1Time - event0Time));
348 Array yg = model_->yGrid(
349 stddevs_, integrationPoints_, event1Time,
350 event0Time, event0 > expiry ? z[k] : 0.0);
351 CubicInterpolation payoff0(
352 z.begin(), z.end(), npvp1[m].begin(),
356 for (Size i = 0; i < yg.size(); i++) {
357 p[i] = payoff0(yg[i], true);
358 }
359 CubicInterpolation payoff1(
360 z.begin(), z.end(), p.begin(),
364 for (Size i = 0; i < z.size() - 1; i++) {
365 price +=
367 0.0, payoff1.cCoefficients()[i],
368 payoff1.bCoefficients()[i],
369 payoff1.aCoefficients()[i], p[i], z[i],
370 z[i], z[i + 1]) *
371 zSpreadDf;
372 }
373 if (extrapolatePayoff_) {
375 price +=
377 0.0, 0.0, 0.0, 0.0,
378 p[z.size() - 2],
379 z[z.size() - 2],
380 z[z.size() - 1], 100.0) *
381 zSpreadDf;
382 price +=
384 0.0, 0.0, 0.0, 0.0, p[0],
385 z[0], -100.0, z[0]) *
386 zSpreadDf;
387 } else {
388 if (type == Option::Call)
389 price +=
391 0.0,
392 payoff1.cCoefficients()
393 [z.size() - 2],
394 payoff1.bCoefficients()
395 [z.size() - 2],
396 payoff1.aCoefficients()
397 [z.size() - 2],
398 p[z.size() - 2],
399 z[z.size() - 2],
400 z[z.size() - 1], 100.0) *
401 zSpreadDf;
402 if (type == Option::Put)
403 price +=
405 0.0,
406 payoff1
407 .cCoefficients()[0],
408 payoff1
409 .bCoefficients()[0],
410 payoff1
411 .aCoefficients()[0],
412 p[0], z[0], -100.0,
413 z[0]) *
414 zSpreadDf;
415 }
416 }
417 }
418
419 npvp0[m][k] = price;
420 }
421 }
422 // end probability computation
423
424 // event date calculations
425
426 if (isEventDate) {
427
428 Real zk = event0 > expiry ? z[k] : y;
429
430 if (isLeg1Fixing) { // if event is a fixing date and
431 // exercise date,
432 // the coupon is part of the exercise into right (by
433 // definition)
434 Size j = std::find(arguments_.leg1FixingDates.begin(),
435 arguments_.leg1FixingDates.end(),
436 event0) -
437 arguments_.leg1FixingDates.begin();
438 Real zSpreadDf =
439 oas_.empty()
440 ? Real(1.0)
441 : std::exp(
442 -oas_->value() *
443 (model_->termStructure()
444 ->dayCounter()
445 .yearFraction(
446 event0,
447 arguments_.leg1PayDates[j])));
448 bool done = false;
449 do {
450 Real amount;
451 if (arguments_.leg1IsRedemptionFlow[j]) {
452 amount = arguments_.leg1Coupons[j];
453 } else {
454 Real estFixing = 0.0;
455 if (ibor1 != nullptr) {
456 estFixing = model_->forwardRate(
457 arguments_.leg1FixingDates[j], event0,
458 zk, ibor1);
459 }
460 if (cms1 != nullptr) {
461 estFixing = model_->swapRate(
462 arguments_.leg1FixingDates[j],
463 cms1->tenor(), event0, zk, cms1);
464 }
465 if (cmsspread1 != nullptr)
466 estFixing =
467 cmsspread1->gearing1() *
468 model_->swapRate(
469 arguments_.leg1FixingDates[j],
470 cmsspread1->swapIndex1()
471 ->tenor(),
472 event0, zk,
473 cmsspread1->swapIndex1()) +
474 cmsspread1->gearing2() *
475 model_->swapRate(
476 arguments_.leg1FixingDates[j],
477 cmsspread1->swapIndex2()
478 ->tenor(),
479 event0, zk,
480 cmsspread1->swapIndex2());
481 Real rate =
482 arguments_.leg1Spreads[j] +
483 arguments_.leg1Gearings[j] * estFixing;
484 if (arguments_.leg1CappedRates[j] !=
485 Null<Real>())
486 rate = std::min(
487 arguments_.leg1CappedRates[j], rate);
488 if (arguments_.leg1FlooredRates[j] !=
489 Null<Real>())
490 rate = std::max(
491 arguments_.leg1FlooredRates[j], rate);
492 amount = rate * arguments_.nominal1[j] *
493 arguments_.leg1AccrualTimes[j];
494 }
495
496 npv0a[k] -=
497 amount *
498 model_->zerobond(arguments_.leg1PayDates[j],
499 event0, zk, discountCurve_) /
500 model_->numeraire(event0Time, zk,
502 zSpreadDf;
503
504 if (j < arguments_.leg1FixingDates.size() - 1) {
505 j++;
506 done =
507 (event0 != arguments_.leg1FixingDates[j]);
508 } else
509 done = true;
510
511 } while (!done);
512 }
513
514 if (isLeg2Fixing) { // if event is a fixing date and
515 // exercise date,
516 // the coupon is part of the exercise into right (by
517 // definition)
518 Size j = std::find(arguments_.leg2FixingDates.begin(),
519 arguments_.leg2FixingDates.end(),
520 event0) -
521 arguments_.leg2FixingDates.begin();
522 Real zSpreadDf =
523 oas_.empty()
524 ? Real(1.0)
525 : std::exp(
526 -oas_->value() *
527 (model_->termStructure()
528 ->dayCounter()
529 .yearFraction(
530 event0,
531 arguments_.leg2PayDates[j])));
532 bool done;
533 do {
534 Real amount;
535 if (arguments_.leg2IsRedemptionFlow[j]) {
536 amount = arguments_.leg2Coupons[j];
537 } else {
538 Real estFixing = 0.0;
539 if (ibor2 != nullptr)
540 estFixing = model_->forwardRate(arguments_.leg2FixingDates[j],event0,zk,ibor2);
541 if (cms2 != nullptr)
542 estFixing = model_->swapRate(arguments_.leg2FixingDates[j],cms2->tenor(),event0,zk,cms2);
543 if (cmsspread2 != nullptr)
544 estFixing =
545 cmsspread2->gearing1() *
546 model_->swapRate(
547 arguments_.leg2FixingDates[j],
548 cmsspread2->swapIndex1()
549 ->tenor(),
550 event0, zk,
551 cmsspread2->swapIndex1()) +
552 cmsspread2->gearing2() *
553 model_->swapRate(
554 arguments_.leg2FixingDates[j],
555 cmsspread2->swapIndex2()
556 ->tenor(),
557 event0, zk,
558 cmsspread2->swapIndex2());
559 Real rate =
560 arguments_.leg2Spreads[j] +
561 arguments_.leg2Gearings[j] * estFixing;
562 if (arguments_.leg2CappedRates[j] !=
563 Null<Real>())
564 rate = std::min(
565 arguments_.leg2CappedRates[j], rate);
566 if (arguments_.leg2FlooredRates[j] !=
567 Null<Real>())
568 rate = std::max(
569 arguments_.leg2FlooredRates[j], rate);
570 amount = rate * arguments_.nominal2[j] *
571 arguments_.leg2AccrualTimes[j];
572 }
573
574 npv0a[k] +=
575 amount *
576 model_->zerobond(arguments_.leg2PayDates[j],
577 event0, zk, discountCurve_) /
578 model_->numeraire(event0Time, zk,
580 zSpreadDf;
581 if (j < arguments_.leg2FixingDates.size() - 1) {
582 j++;
583 done =
584 (event0 != arguments_.leg2FixingDates[j]);
585 } else
586 done = true;
587
588 } while (!done);
589 }
590
591 if (isExercise) {
592 Size j = std::find(arguments_.exercise->dates().begin(),
593 arguments_.exercise->dates().end(),
594 event0) -
595 arguments_.exercise->dates().begin();
596 Real rebate = 0.0;
597 Real zSpreadDf = 1.0;
598 Date rebateDate = event0;
599 if (rebatedExercise_ != nullptr) {
600 rebate = rebatedExercise_->rebate(j);
601 rebateDate = rebatedExercise_->rebatePaymentDate(j);
602 zSpreadDf =
603 oas_.empty()
604 ? Real(1.0)
605 : std::exp(-oas_->value() *
606 (model_->termStructure()
607 ->dayCounter()
608 .yearFraction(event0,
609 rebateDate)));
610 }
611 Real exerciseValue =
612 (type == Option::Call ? 1.0 : -1.0) * npv0a[k] +
613 rebate * model_->zerobond(rebateDate, event0) *
614 zSpreadDf / model_->numeraire(event0Time, zk,
616
617 if (considerProbabilities && probabilities_ != None) {
618 if (exIdx == noEx) {
619 // if true we are at the latest date,
620 // so we init
621 // the no call probability
622 npvp0.back()[k] =
624 ? Real(1.0)
625 : 1.0 / (model_->zerobond(
626 event0Time, 0.0, 0.0,
628 model_->numeraire(
629 event0, z[k],
631 }
632 if (exerciseValue >= npv0[k]) {
633 npvp0[exIdx-1][k] =
635 ? Real(1.0)
636 : 1.0 / (model_->zerobond(
637 event0Time, 0.0, 0.0,
639 model_->numeraire(
640 event0Time, z[k],
642 for (Size ii = exIdx; ii < noEx+1; ++ii)
643 npvp0[ii][k] = 0.0;
644 }
645 }
646 // end probability computation
647
648 npv0[k] = std::max(npv0[k], exerciseValue);
649 }
650 }
651 }
652
653 if(isExercise)
654 --exIdx;
655
656 npv1.swap(npv0);
657 npv1a.swap(npv0a);
658
659 // for probability computation
660 if(considerProbabilities && probabilities_ != None) {
661 for(Size i=0;i<npvp0.size();++i) {
662 npvp1[i].swap(npvp0[i]);
663 }
664 }
665 // end probability computation
666
667 event1 = event0;
668 event1Time = event0Time;
669
670 } while (--idx >= -1);
671
672 std::pair<Real, Real> res(
673 npv1[0] * model_->numeraire(event1Time, y, discountCurve_),
674 npv1a[0] * model_->numeraire(event1Time, y, discountCurve_) *
675 (type == Option::Call ? 1.0 : -1.0));
676
677 // for probability computation
678 if (considerProbabilities && probabilities_ != None) {
679 std::vector<Real> prob(noEx+1);
680 for (Size i = 0; i < noEx+1; i++) {
681 prob[i] = npvp1[i][0] *
683 ? 1.0
684 : model_->numeraire(0.0, 0.0, discountCurve_));
685 }
686 results_.additionalResults["probabilities"] = prob;
687 }
688 // end probability computation
689
690 return res;
691 }
692}
1-D array used in linear algebra.
Definition: array.hpp:52
void swap(Array &) noexcept
Definition: array.hpp:546
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
Cubic interpolation between discrete points.
const std::vector< Real > & cCoefficients() const
const std::vector< Real > & bCoefficients() const
const std::vector< Real > & aCoefficients() const
Concrete date class.
Definition: date.hpp:125
std::pair< Real, Real > npvs(const Date &expiry, Real y, bool includeExerciseOnxpiry, bool considerProbabilities=false) const
const Array initialGuess(const Date &expiry) const override
Real underlyingNpv(const Date &expiry, Real y) const override
static Real gaussianShiftedPolynomialIntegral(Real a, Real b, Real c, Real d, Real e, Real h, Real x0, Real x1)
template class providing a null value for a given type.
Definition: null.hpp:76
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
void swap(Array &v, Array &w) noexcept
Definition: array.hpp:903