QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
conundrumpricer.cpp
Go to the documentation of this file.
1/*
2 Copyright (C) 2006 Giorgio Facchinetti
3 Copyright (C) 2006 Mario Pucci
4 Copyright (C) 2023 Andre Miemiec
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
16 This program is distributed in the hope that it will be useful, but
17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 or FITNESS FOR A PARTICULAR PURPOSE. See the license for more details.
19 */
20
21/*! \file conundrumpricer.hpp
22 \brief
23*/
24
27#include <ql/functional.hpp>
38#include <ql/time/schedule.hpp>
39#include <utility>
40
41namespace QuantLib {
42
43//===========================================================================//
44// Market Quoted Options Pricer //
45//===========================================================================//
46
48 Rate forwardValue,
49 Date expiryDate,
50 const Period& swapTenor,
51 const ext::shared_ptr<SwaptionVolatilityStructure>& volatilityStructure) :
52 forwardValue_(forwardValue), expiryDate_(expiryDate), swapTenor_(swapTenor),
53 volatilityStructure_(volatilityStructure),
54 smile_(volatilityStructure_->smileSection(expiryDate_, swapTenor_)) {
55 QL_REQUIRE((volatilityStructure->volatilityType() == Normal) ||
56 (volatilityStructure->volatilityType() == ShiftedLognormal &&
57 close_enough(volatilityStructure->shift(expiryDate, swapTenor), 0.0)),
58 "VanillaOptionPricer: a normal or a zero-shift lognormal volatility is required");
59 }
60
62 Option::Type optionType,
63 Real deflator) const {
64 const Real variance = smile_->variance(strike);
65 if (volatilityStructure_->volatilityType() == ShiftedLognormal) {
66 return deflator *
67 blackFormula(optionType, strike, forwardValue_, std::sqrt(variance));
68 } else {
69 return deflator *
70 bachelierBlackFormula(optionType, strike, forwardValue_, std::sqrt(variance));
71 }
72 }
73
74
75//===========================================================================//
76// HaganPricer //
77//===========================================================================//
79 GFunctionFactory::YieldCurveModel modelOfYieldCurve,
80 Handle<Quote> meanReversion)
81 : CmsCouponPricer(swaptionVol), modelOfYieldCurve_(modelOfYieldCurve),
82 meanReversion_(std::move(meanReversion)) {
84 }
85
87 coupon_ = dynamic_cast<const CmsCoupon*>(&coupon);
88 QL_REQUIRE(coupon_, "CMS coupon needed");
91 Time accrualPeriod = coupon_->accrualPeriod();
92 QL_REQUIRE(accrualPeriod != 0.0, "null accrual period");
93
96 const ext::shared_ptr<SwapIndex>& swapIndex = coupon_->swapIndex();
97 rateCurve_ = swapIndex->discountingTermStructure().empty() ? *(swapIndex->forwardingTermStructure()) : *(swapIndex->discountingTermStructure());
98
100
101 if(paymentDate_ > today)
102 discount_ = rateCurve_->discount(paymentDate_);
103 else discount_= 1.;
104
105 spreadLegValue_ = spread_ * accrualPeriod * discount_;
106
107 if (fixingDate_ > today){
108 swapTenor_ = swapIndex->tenor();
109 ext::shared_ptr<VanillaSwap> swap = swapIndex->underlyingSwap(fixingDate_);
110
111 swapRateValue_ = swap->fairRate();
112
113 static const Spread bp = 1.0e-4;
114 annuity_ = std::fabs(swap->fixedLegBPS()/bp);
115
116 Size q = swapIndex->fixedLegTenor().frequency();
117 const Schedule& schedule = swap->fixedSchedule();
118 const DayCounter& dc = swapIndex->dayCounter();
119 //const DayCounter dc = coupon.dayCounter();
120 Time startTime = dc.yearFraction(rateCurve_->referenceDate(),
121 swap->startDate());
122 Time swapFirstPaymentTime =
123 dc.yearFraction(rateCurve_->referenceDate(), schedule.date(1));
124 Time paymentTime = dc.yearFraction(rateCurve_->referenceDate(),
126 Real delta = (paymentTime-startTime) / (swapFirstPaymentTime-startTime);
127
128 switch (modelOfYieldCurve_) {
131 break;
134 break;
136 Handle<Quote> nullMeanReversionQuote(ext::shared_ptr<Quote>(new SimpleQuote(0.0)));
138 }
139 break;
142 break;
143 default:
144 QL_FAIL("unknown/illegal gFunction type");
145 }
146 vanillaOptionPricer_= ext::shared_ptr<VanillaOptionPricer>(new
149 }
150 }
151
153
156 }
157
158 Real HaganPricer::capletPrice(Rate effectiveCap) const {
159 // caplet is equivalent to call option on fixing
161 if (fixingDate_ <= today) {
162 // the fixing is determined
163 const Rate Rs =
164 std::max(coupon_->swapIndex()->fixing(fixingDate_)-effectiveCap, 0.);
165 Rate price = (gearing_*Rs)*(coupon_->accrualPeriod()*discount_);
166 return price;
167 }
168 else {
169 Real capletPrice = 0.0;
170
171 if (swaptionVolatility()->volatilityType() == ShiftedLognormal)
172 {
173 Real cutoffNearZero = 1e-10;
174
175 if (effectiveCap < cutoffForCaplet_) {
176 Rate effectiveStrikeForMax = std::max(effectiveCap, cutoffNearZero);
177 capletPrice = optionletPrice(Option::Call, effectiveStrikeForMax);
178 }
179 }
180 else
181 {
182 capletPrice = optionletPrice(Option::Call, effectiveCap);
183 }
184 return gearing_ * capletPrice;
185 }
186 }
187
188 Rate HaganPricer::capletRate(Rate effectiveCap) const {
189 return capletPrice(effectiveCap)/(coupon_->accrualPeriod()*discount_);
190 }
191
192 Real HaganPricer::floorletPrice(Rate effectiveFloor) const {
193 // floorlet is equivalent to put option on fixing
195 if (fixingDate_ <= today) {
196 // the fixing is determined
197 const Rate Rs =
198 std::max(effectiveFloor-coupon_->swapIndex()->fixing(fixingDate_),0.);
199 Rate price = (gearing_*Rs)*(coupon_->accrualPeriod()*discount_);
200 return price;
201 }
202 else {
203 Real floorletPrice = 0.0;
204 if(swaptionVolatility()->volatilityType() == ShiftedLognormal)
205 {
206 Real cutoffNearZero = 1e-10;
207
208 if (effectiveFloor > cutoffForFloorlet_) {
209 Rate effectiveStrikeForMin = std::max(effectiveFloor, cutoffNearZero);
210 floorletPrice = optionletPrice(Option::Put, effectiveStrikeForMin);
211 }
212 }
213 else
214 {
215 floorletPrice = optionletPrice(Option::Put, effectiveFloor);
216 }
217
218 return gearing_ * floorletPrice;
219 }
220 }
221
222
223 Rate HaganPricer::floorletRate(Rate effectiveFloor) const {
224 return floorletPrice(effectiveFloor)/(coupon_->accrualPeriod()*discount_);
225 }
226
227//===========================================================================//
228// NumericHaganPricer //
229//===========================================================================//
230
231 namespace {
232
233 class VariableChange {
234 public:
235 VariableChange(ext::function<Real (Real)>& f,
236 Real a, Real b, Size k)
237 : a_(a), width_(b-a), f_(f), k_(k) {}
238 Real value(Real x) const {
239 Real newVar;
240 Real temp = width_;
241 for (Size i = 1; i < k_ ; ++i) {
242 temp *= x;
243 }
244 newVar = a_ + x* temp;
245 return f_(newVar) * k_* temp;
246 }
247 private:
249 ext::function<Real (Real)> f_;
251 };
252
253 class Spy {
254 public:
255 explicit Spy(ext::function<Real(Real)> f) : f_(std::move(f)) {}
256 Real value(Real x){
257 abscissas.push_back(x);
258 Real value = f_(x);
259 functionValues.push_back(value);
260 return value;
261 }
262 private:
263 ext::function<Real (Real)> f_;
264 std::vector<Real> abscissas;
265 std::vector<Real> functionValues;
266 };
267
268 }
269
271 GFunctionFactory::YieldCurveModel modelOfYieldCurve,
272 const Handle<Quote>& meanReversion,
273 Real lowerLimit,
274 Real upperLimit,
275 Real precision,
276 Real hardUpperLimit)
277 : HaganPricer(swaptionVol, modelOfYieldCurve, meanReversion),
278 lowerLimit_(lowerLimit), upperLimit_(upperLimit),
279 precision_(precision), hardUpperLimit_(hardUpperLimit) {}
280
282
283 Real result =.0;
284 //double abserr =.0;
285 //double alpha = 1.0;
286 //double epsabs = precision_;
287 //double epsrel = 1.0; // we are interested only in absolute precision
288 //size_t neval =0;
289
290 // we use the non adaptive algorithm only for semi infinite interval
291 if (a > 0) {
292
293 // we estimate the actual boundary by testing integrand values
294 Real upperBoundary = 2 * a;
295 while (integrand(upperBoundary) > precision_)
296 upperBoundary *= 2.0;
297 // sometimes b < a because of a wrong estimation of b based on stdev
298 if (b > a)
299 upperBoundary = std::min(upperBoundary, b);
300
301 ext::function<Real(Real)> f;
303 gaussKronrodNonAdaptive(precision_, 1000000, 1.0);
304 // if the integration intervall is wide enough we use the
305 // following change variable x -> a + (b-a)*(t/(a-b))^3
306 upperBoundary = std::max(a, std::min(upperBoundary, hardUpperLimit_));
307 if (upperBoundary > 2 * a) {
308 Size k = 3;
309 ext::function<Real(Real)> temp = ext::cref(integrand);
310 VariableChange variableChange(temp, a, upperBoundary, k);
311 f = [&](Real _x) { return variableChange.value(_x); };
312 result = gaussKronrodNonAdaptive(f, .0, 1.0);
313 }
314 else {
315 f = ext::cref(integrand);
316 result = gaussKronrodNonAdaptive(f, a, upperBoundary);
317 }
318
319 // if the expected precision has not been reached we use the old algorithm
320 if (!gaussKronrodNonAdaptive.integrationSuccess()) {
321 const GaussKronrodAdaptive integrator(precision_, 100000);
322 b = std::max(a, std::min(b, hardUpperLimit_));
323 result = integrator(integrand, a, b);
324 }
325 }
326 else { // if a < b we use the old algorithm
327 b = std::max(a, std::min(b, hardUpperLimit_));
328 if (swaptionVolatility()->volatilityType() == ShiftedLognormal) {
329 const GaussKronrodAdaptive integrator(precision_, 100000);
330 result = integrator(integrand, a, b);
331 }
332 else //Normal and floorlet
333 {
334 const GaussKronrodNonAdaptive integrator(precision_, 100000, 1.0);
335 //Doing an integral from -inf to strike, where strike itself is negative,
336 //the GaussKronrodAdaptive quadrature rule throws an exceptions in
337 //GaussKronrodAdaptive::integrateRecursively due to maximum number of
338 //function evaluations exceed.
339 result = integrator(integrand, a, b);
340 }
341 }
342 return result;
343 }
344
345
347 Option::Type optionType, Real strike) const {
348
349 ext::shared_ptr<ConundrumIntegrand> integrand(new
352 swapRateValue_, strike, optionType));
355 Real a, b, integralValue;
356 if (optionType==Option::Call) {
358 // while(upperLimit_ <= strike){
359 // stdDeviationsForUpperLimit_ += 1.;
360 // upperLimit_ = resetUpperLimit(stdDeviationsForUpperLimit_);
361 // }
362 integralValue = integrate(strike, upperLimit_, *integrand);
363 //refineIntegration(integralValue, *integrand);
364 } else {
366 a = std::min(strike, lowerLimit_);
367 b = strike;
368 integralValue = integrate(a, b, *integrand);
369 }
370
371 Real dFdK = integrand->firstDerivativeOfF(strike);
372 Real swaptionPrice =
373 (*vanillaOptionPricer_)(strike, optionType, annuity_);
374
375 // v. HAGAN, Conundrums..., formule 2.17a, 2.18a
377 ((1 + dFdK) * swaptionPrice + Integer(optionType) * integralValue);
378 }
379
381
383 if (fixingDate_ <= today) {
384 // the fixing is determined
385 const Rate Rs = coupon_->swapIndex()->fixing(fixingDate_);
387 return price;
388 } else {
390 Real atmFloorletPrice = optionletPrice(Option::Put, swapRateValue_);
392 + atmCapletPrice - atmFloorletPrice)
394 }
395 }
396
398 const ConundrumIntegrand& integrand) const {
399 Real percDiff = 1000.;
400 while(std::fabs(percDiff) < refiningIntegrationTolerance_){
404 Real diff = integrate(lowerLimit, upperLimit_,integrand);
405 percDiff = diff/integralValue;
406 integralValue += diff;
407 }
408 return integralValue;
409 }
410
411 Real NumericHaganPricer::resetUpperLimit(Real stdDeviationsForUpperLimit) const {
412
413 // return 1.0;
414 Real variance =
416
417 if (swaptionVolatility()->volatilityType() == ShiftedLognormal) {
418 return swapRateValue_ * std::exp(stdDeviationsForUpperLimit * std::sqrt(variance));
419 } else {
420 return swapRateValue_ + stdDeviationsForUpperLimit * std::sqrt(variance);
421 }
422 }
423
424
425 Real NumericHaganPricer::resetLowerLimit(Real stdDeviationsForUpperLimit) const {
426 // return -1.0;
427 Real variance =
429
430 if (swaptionVolatility()->volatilityType() == ShiftedLognormal) {
431 return lowerLimit_;
432 } else {
433 return swapRateValue_ - stdDeviationsForUpperLimit * std::sqrt(variance);
434 }
435 }
436
437
438//===========================================================================//
439// ConundrumIntegrand //
440//===========================================================================//
441
443 ext::shared_ptr<VanillaOptionPricer> o,
444 const ext::shared_ptr<YieldTermStructure>&,
445 ext::shared_ptr<GFunction> gFunction,
446 Date fixingDate,
447 Date paymentDate,
448 Real annuity,
449 Real forwardValue,
450 Real strike,
451 Option::Type optionType)
452 : vanillaOptionPricer_(std::move(o)), forwardValue_(forwardValue), annuity_(annuity),
453 fixingDate_(fixingDate), paymentDate_(paymentDate), strike_(strike), optionType_(optionType),
454 gFunction_(std::move(gFunction)) {}
455
457 strike_ = strike;
458 }
459
461 return strike_;
462 }
463
465 return annuity_;
466 }
467
469 return fixingDate_;
470 }
471
473 const Real Gx = (*gFunction_)(x);
474 const Real GR = (*gFunction_)(forwardValue_);
475 return (x - strike_) * (Gx/GR - 1.0);
476 }
477
479 const Real Gx = (*gFunction_)(x);
480 const Real GR = (*gFunction_)(forwardValue_) ;
481 const Real G1 = gFunction_->firstDerivative(x);
482 return (Gx/GR - 1.0) + G1/GR * (x - strike_);
483 }
484
486 const Real GR = (*gFunction_)(forwardValue_) ;
487 const Real G1 = gFunction_->firstDerivative(x);
488 const Real G2 = gFunction_->secondDerivative(x);
489 return 2.0 * G1/GR + (x - strike_) * G2/GR;
490 }
491
493 const Real option = (*vanillaOptionPricer_)(x, optionType_, annuity_);
494 return option * secondDerivativeOfF(x);
495 }
496
497
498
499//===========================================================================//
500// AnalyticHaganPricer //
501//===========================================================================//
502
504 const Handle<SwaptionVolatilityStructure>& swaptionVol,
505 GFunctionFactory::YieldCurveModel modelOfYieldCurve,
507 : HaganPricer(swaptionVol, modelOfYieldCurve, meanReversion)
508 { }
509
510 //Hagan, 3.5b, 3.5c
512 Real strike) const {
513 Real variance = swaptionVolatility()->blackVariance(fixingDate_,
516 Real firstDerivativeOfGAtForwardValue = gFunction_->firstDerivative(
518 Real price = 0.0;
519
520 Real CK = (*vanillaOptionPricer_)(strike, optionType, annuity_);
521 price += (discount_/annuity_)*CK;
522
523 if (swaptionVolatility()->volatilityType() == ShiftedLognormal) {
524 const Real sqrtSigma2T = std::sqrt(variance);
525 const Real lnRoverK = std::log(swapRateValue_ / strike);
526 const Real d32 = (lnRoverK + 1.5*variance) / sqrtSigma2T;
527 const Real d12 = (lnRoverK + .5*variance) / sqrtSigma2T;
528 const Real dminus12 = (lnRoverK - .5*variance) / sqrtSigma2T;
529
530 CumulativeNormalDistribution cumulativeOfNormal;
531 auto sign = Integer(optionType);
532 const Real N32 = cumulativeOfNormal(sign*d32);
533 const Real N12 = cumulativeOfNormal(sign*d12);
534 const Real Nminus12 = cumulativeOfNormal(sign*dminus12);
535
536 price += sign * firstDerivativeOfGAtForwardValue * annuity_ *
537 swapRateValue_ * (swapRateValue_ * std::exp(variance) * N32 -
538 (swapRateValue_ + strike) * N12 + strike * Nminus12);
539 } else {
540 const Real sqrtSigma2T = std::sqrt(variance);
541 const Real d = (swapRateValue_ - strike) / sqrtSigma2T;
542
543 CumulativeNormalDistribution cumulativeOfNormal;
544 auto sign = Integer(optionType);
545
546 const Real N = cumulativeOfNormal(sign*d);
547 price += sign * firstDerivativeOfGAtForwardValue * annuity_ * variance * N;
548 }
549
550 price *= coupon_->accrualPeriod();
551 return price;
552 }
553
554 //Hagan 3.4c
556 {
557
559 if (fixingDate_ <= today) {
560 // the fixing is determined
561 const Rate Rs = coupon_->swapIndex()->fixing(fixingDate_);
563 return price;
564 } else {
568 Real firstDerivativeOfGAtForwardValue(gFunction_->firstDerivative(
570 Real price = 0.0;
571 price += discount_*swapRateValue_;
572 if (swaptionVolatility()->volatilityType()==ShiftedLognormal) {
573 price += firstDerivativeOfGAtForwardValue * annuity_*swapRateValue_*
574 swapRateValue_*(std::exp(variance) - 1.);
575 } else {
576 price += firstDerivativeOfGAtForwardValue * annuity_*variance;
577 }
578 return (gearing_ * price +spread_*discount_)* coupon_->accrualPeriod();
579 }
580 }
581
582
583
584//===========================================================================//
585// GFunctionStandard //
586//===========================================================================//
587
589 Real n = static_cast<Real>(swapLength_) * q_;
590 return x / std::pow((1.0 + x/q_), delta_) * 1.0 /
591 (1.0 - 1.0 / std::pow((1.0 + x/q_), n));
592 }
593
595 Real n = static_cast<Real>(swapLength_) * q_;
596 Real a = 1.0 + x / q_;
597 Real AA = a - delta_/q_ * x;
598 Real B = std::pow(a,(n - delta_ - 1.0))/(std::pow(a,n) - 1.0);
599
600 Real secNum = n * x * std::pow(a,(n-1.0));
601 Real secDen = q_ * std::pow(a, delta_) * (std::pow(a, n) - 1.0) *
602 (std::pow(a, n) - 1.0);
603 Real sec = secNum / secDen;
604
605 return AA * B - sec;
606 }
607
609 Real n = static_cast<Real>(swapLength_) * q_;
610 Real a = 1.0 + x/q_;
611 Real AA = a - delta_/q_ * x;
612 Real A1 = (1.0 - delta_)/q_;
613 Real B = std::pow(a,(n - delta_ - 1.0))/(std::pow(a,n) - 1.0);
614 Real Num = (1.0 + delta_ - n) * std::pow(a, (n-delta_-2.0)) -
615 (1.0 + delta_) * std::pow(a, (2.0*n-delta_-2.0));
616 Real Den = (std::pow(a, n) - 1.0) * (std::pow(a, n) - 1.0);
617 Real B1 = 1.0 / q_ * Num / Den;
618
619 Real C = x / std::pow(a, delta_);
620 Real C1 = (std::pow(a, delta_)
621 - delta_ /q_ * x * std::pow(a, (delta_ - 1.0))) / std::pow(a, 2 * delta_);
622
623 Real D = std::pow(a, (n-1.0))/ ((std::pow(a, n) - 1.0) * (std::pow(a, n) - 1.0));
624 Real D1 = ((n - 1.0) * std::pow(a, (n-2.0)) * (std::pow(a, n) - 1.0)
625 - 2 * n * std::pow(a, (2 * (n-1.0))))
626 / (q_ * (std::pow(a, n) - 1.0)*(std::pow(a, n) - 1.0)*(std::pow(a, n) - 1.0));
627
628 return A1 * B + AA * B1 - n/q_ * (C1 * D + C * D1);
629 }
630
632 Real delta, Size swapLength) {
633 return ext::shared_ptr<GFunction>(new GFunctionStandard(q, delta, swapLength));
634 }
635
636//===========================================================================//
637// GFunctionExactYield //
638//===========================================================================//
639
641
642 const ext::shared_ptr<SwapIndex>& swapIndex = coupon.swapIndex();
643 const ext::shared_ptr<VanillaSwap>& swap =
644 swapIndex->underlyingSwap(coupon.fixingDate());
645
646 const Schedule& schedule = swap->fixedSchedule();
648 swapIndex->forwardingTermStructure();
649
650 const DayCounter& dc = swapIndex->dayCounter();
651
652 Real swapStartTime = dc.yearFraction(rateCurve->referenceDate(),
653 schedule.startDate());
654 Real swapFirstPaymentTime = dc.yearFraction(rateCurve->referenceDate(),
655 schedule.date(1));
656
657 Real paymentTime = dc.yearFraction(rateCurve->referenceDate(),
658 coupon.date());
659
660 delta_ = (paymentTime-swapStartTime) / (swapFirstPaymentTime-swapStartTime);
661
662 const Leg& fixedLeg(swap->fixedLeg());
663 Size n = fixedLeg.size();
664 accruals_.reserve(n);
665 for (Size i=0; i<n; ++i) {
666 ext::shared_ptr<Coupon> cpn =
667 ext::dynamic_pointer_cast<Coupon>(fixedLeg[i]);
668 accruals_.push_back(cpn->accrualPeriod());
669 }
670 }
671
673 Real product = 1.;
674 for (Real accrual : accruals_) {
675 product *= 1. / (1. + accrual * x);
676 }
677 return x*std::pow(1.+ accruals_[0]*x,-delta_)*(1./(1.-product));
678 }
679
681 Real c = -1.;
682 Real derC = 0.;
683 std::vector<Real> b;
684 b.reserve(accruals_.size());
685 for (Real accrual : accruals_) {
686 Real temp = 1.0 / (1.0 + accrual * x);
687 b.push_back(temp);
688 c *= temp;
689 derC += accrual * temp;
690 }
691 c += 1.;
692 c = 1./c;
693 derC *= (c-c*c);
694
695 return -delta_*accruals_[0]*std::pow(b[0],delta_+1.)*x*c+
696 std::pow(b[0],delta_)*c+ std::pow(b[0],delta_)*x*derC;
697 //Real dx = 1.0e-8;
698 //return (operator()(x+dx)-operator()(x-dx))/(2.0*dx);
699 }
700
702 Real c = -1.;
703 Real sum = 0.;
704 Real sumOfSquare = 0.;
705 std::vector<Real> b;
706 b.reserve(accruals_.size());
707 for (Real accrual : accruals_) {
708 Real temp = 1.0 / (1.0 + accrual * x);
709 b.push_back(temp);
710 c *= temp;
711 sum += accrual * temp;
712 sumOfSquare += std::pow(accrual * temp, 2.0);
713 }
714 c += 1.;
715 c = 1./c;
716 Real derC =sum*(c-c*c);
717
718 return (-delta_*accruals_[0]*std::pow(b[0],delta_+1.)*c+ std::pow(b[0],delta_)*derC)*
719 (-delta_*accruals_[0]*b[0]*x + 1. + x*(1.-c)*sum)+
720 std::pow(b[0],delta_)*c*(delta_*std::pow(accruals_[0]*b[0],2.)*x - delta_* accruals_[0]*b[0] -
721 x*derC*sum + (1.-c)*sum - x*(1.-c)*sumOfSquare);
722 //Real dx = 1.0e-8;
723 //return (firstDerivative(x+dx)-firstDerivative(x-dx))/(2.0*dx);
724 }
725
726 ext::shared_ptr<GFunction> GFunctionFactory::newGFunctionExactYield(const CmsCoupon& coupon) {
727 return ext::shared_ptr<GFunction>(new GFunctionExactYield(coupon));
728 }
729
730
731
732//===========================================================================//
733// GFunctionWithShifts //
734//===========================================================================//
735
737 Handle<Quote> meanReversion)
738 : meanReversion_(std::move(meanReversion)) {
739
740 const ext::shared_ptr<SwapIndex>& swapIndex = coupon.swapIndex();
741 const ext::shared_ptr<VanillaSwap>& swap = swapIndex->underlyingSwap(coupon.fixingDate());
742
743 swapRateValue_ = swap->fairRate();
744
745 objectiveFunction_ = ext::make_shared<ObjectiveFunction>(*this, swapRateValue_);
746
747 const Schedule& schedule = swap->fixedSchedule();
749 swapIndex->forwardingTermStructure();
750 const DayCounter& dc = swapIndex->dayCounter();
751
752 swapStartTime_ = dc.yearFraction(rateCurve->referenceDate(),
753 schedule.startDate());
754 discountAtStart_ = rateCurve->discount(schedule.startDate());
755
756 Real paymentTime = dc.yearFraction(rateCurve->referenceDate(),
757 coupon.date());
758
759 shapedPaymentTime_ = shapeOfShift(paymentTime);
760
761 const Leg& fixedLeg(swap->fixedLeg());
762 Size n = fixedLeg.size();
763 accruals_.reserve(n);
764 shapedSwapPaymentTimes_.reserve(n);
765 swapPaymentDiscounts_.reserve(n);
766 for(Size i=0; i<n; ++i) {
767 ext::shared_ptr<Coupon> cpn =
768 ext::dynamic_pointer_cast<Coupon>(fixedLeg[i]);
769 accruals_.push_back(cpn->accrualPeriod());
770 const Date paymentDate(cpn->date());
771 const Real swapPaymentTime(dc.yearFraction(rateCurve->referenceDate(), paymentDate));
772 shapedSwapPaymentTimes_.push_back(shapeOfShift(swapPaymentTime));
773 swapPaymentDiscounts_.push_back(rateCurve->discount(paymentDate));
774 }
776 }
777
779 const Real calibratedShift = calibrationOfShift(Rs);
780 return Rs* functionZ(calibratedShift);
781 }
782
784 return std::exp(-shapedPaymentTime_*x)
785 / (1.-discountRatio_*std::exp(-shapedSwapPaymentTimes_.back()*x));
786 }
787
789 Real sqrtDenominator = 0;
790 Real derSqrtDenominator = 0;
791 for(Size i=0; i<accruals_.size(); i++) {
792 sqrtDenominator += accruals_[i]*swapPaymentDiscounts_[i]
793 *std::exp(-shapedSwapPaymentTimes_[i]*x);
794 derSqrtDenominator -= shapedSwapPaymentTimes_[i]* accruals_[i]*swapPaymentDiscounts_[i]
795 *std::exp(-shapedSwapPaymentTimes_[i]*x);
796 }
797 const Real denominator = sqrtDenominator* sqrtDenominator;
798
799 Real numerator = 0;
800 numerator += shapedSwapPaymentTimes_.back()* swapPaymentDiscounts_.back()*
801 std::exp(-shapedSwapPaymentTimes_.back()*x)*sqrtDenominator;
802 numerator -= (discountAtStart_ - swapPaymentDiscounts_.back()* std::exp(-shapedSwapPaymentTimes_.back()*x))*
803 derSqrtDenominator;
804 QL_REQUIRE(denominator!=0, "GFunctionWithShifts::derRs_derX: denominator == 0");
805 return numerator/denominator;
806 }
807
809 Real denOfRfunztion = 0.;
810 Real derDenOfRfunztion = 0.;
811 Real der2DenOfRfunztion = 0.;
812 for(Size i=0; i<accruals_.size(); i++) {
813 denOfRfunztion += accruals_[i]*swapPaymentDiscounts_[i]
814 *std::exp(-shapedSwapPaymentTimes_[i]*x);
815 derDenOfRfunztion -= shapedSwapPaymentTimes_[i]* accruals_[i]*swapPaymentDiscounts_[i]
816 *std::exp(-shapedSwapPaymentTimes_[i]*x);
817 der2DenOfRfunztion+= shapedSwapPaymentTimes_[i]*shapedSwapPaymentTimes_[i]* accruals_[i]*
818 swapPaymentDiscounts_[i]*std::exp(-shapedSwapPaymentTimes_[i]*x);
819 }
820
821 const Real denominator = std::pow(denOfRfunztion, 4);
822
823 Real numOfDerR = 0;
824 numOfDerR += shapedSwapPaymentTimes_.back()* swapPaymentDiscounts_.back()*
825 std::exp(-shapedSwapPaymentTimes_.back()*x)*denOfRfunztion;
826 numOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.back()* std::exp(-shapedSwapPaymentTimes_.back()*x))*
827 derDenOfRfunztion;
828
829 const Real denOfDerR = std::pow(denOfRfunztion,2);
830
831 Real derNumOfDerR = 0.;
832 derNumOfDerR -= shapedSwapPaymentTimes_.back()*shapedSwapPaymentTimes_.back()* swapPaymentDiscounts_.back()*
833 std::exp(-shapedSwapPaymentTimes_.back()*x)*denOfRfunztion;
834 derNumOfDerR += shapedSwapPaymentTimes_.back()* swapPaymentDiscounts_.back()*
835 std::exp(-shapedSwapPaymentTimes_.back()*x)*derDenOfRfunztion;
836
837 derNumOfDerR -= (shapedSwapPaymentTimes_.back()*swapPaymentDiscounts_.back()*
838 std::exp(-shapedSwapPaymentTimes_.back()*x))* derDenOfRfunztion;
839 derNumOfDerR -= (discountAtStart_ - swapPaymentDiscounts_.back()* std::exp(-shapedSwapPaymentTimes_.back()*x))*
840 der2DenOfRfunztion;
841
842 const Real derDenOfDerR = 2*denOfRfunztion*derDenOfRfunztion;
843
844 const Real numerator = derNumOfDerR*denOfDerR -numOfDerR*derDenOfDerR;
845 QL_REQUIRE(denominator!=0, "GFunctionWithShifts::der2Rs_derX2: denominator == 0");
846 return numerator/denominator;
847 }
848
850 const Real sqrtDenominator = (1.-discountRatio_*std::exp(-shapedSwapPaymentTimes_.back()*x));
851 const Real denominator = sqrtDenominator* sqrtDenominator;
852 QL_REQUIRE(denominator!=0, "GFunctionWithShifts::derZ_derX: denominator == 0");
853
854 Real numerator = 0;
855 numerator -= shapedPaymentTime_* std::exp(-shapedPaymentTime_*x)* sqrtDenominator;
856 numerator -= shapedSwapPaymentTimes_.back()* std::exp(-shapedPaymentTime_*x)* (1.-sqrtDenominator);
857
858 return numerator/denominator;
859 }
860
862 const Real denOfZfunction = (1.-discountRatio_*std::exp(-shapedSwapPaymentTimes_.back()*x));
863 const Real derDenOfZfunction = shapedSwapPaymentTimes_.back()*discountRatio_*std::exp(-shapedSwapPaymentTimes_.back()*x);
864 const Real denominator = std::pow(denOfZfunction, 4);
865 QL_REQUIRE(denominator!=0, "GFunctionWithShifts::der2Z_derX2: denominator == 0");
866
867 Real numOfDerZ = 0;
868 numOfDerZ -= shapedPaymentTime_* std::exp(-shapedPaymentTime_*x)* denOfZfunction;
869 numOfDerZ -= shapedSwapPaymentTimes_.back()* std::exp(-shapedPaymentTime_*x)* (1.-denOfZfunction);
870
871 const Real denOfDerZ = std::pow(denOfZfunction,2);
872 const Real derNumOfDerZ = (-shapedPaymentTime_* std::exp(-shapedPaymentTime_*x)*
873 (-shapedPaymentTime_+(shapedPaymentTime_*discountRatio_-
874 shapedSwapPaymentTimes_.back()*discountRatio_)* std::exp(-shapedSwapPaymentTimes_.back()*x))
875 -shapedSwapPaymentTimes_.back()*std::exp(-shapedPaymentTime_*x)*
876 (shapedPaymentTime_*discountRatio_- shapedSwapPaymentTimes_.back()*discountRatio_)*
877 std::exp(-shapedSwapPaymentTimes_.back()*x));
878
879 const Real derDenOfDerZ = 2*denOfZfunction*derDenOfZfunction;
880 const Real numerator = derNumOfDerZ*denOfDerZ -numOfDerZ*derDenOfDerZ;
881
882 return numerator/denominator;
883 }
884
886 //Real dRs = 1.0e-8;
887 //return (operator()(Rs+dRs)-operator()(Rs-dRs))/(2.0*dRs);
888 const Real calibratedShift = calibrationOfShift(Rs);
889 return functionZ(calibratedShift) + Rs * derZ_derX(calibratedShift)/derRs_derX(calibratedShift);
890 }
891
893 //Real dRs = 1.0e-8;
894 //return (firstDerivative(Rs+dRs)-firstDerivative(Rs-dRs))/(2.0*dRs);
895 const Real calibratedShift = calibrationOfShift(Rs);
896 return 2.*derZ_derX(calibratedShift)/derRs_derX(calibratedShift) +
897 Rs * der2Z_derX2(calibratedShift)/std::pow(derRs_derX(calibratedShift),2.)-
898 Rs * derZ_derX(calibratedShift)*der2Rs_derX2(calibratedShift)/
899 std::pow(derRs_derX(calibratedShift),3.);
900 }
901
903 Real result = 0;
904 derivative_ = 0;
905 for(Size i=0; i<o_.accruals_.size(); i++) {
906 Real temp = o_.accruals_[i]*o_.swapPaymentDiscounts_[i]
907 *std::exp(-o_.shapedSwapPaymentTimes_[i]*x);
908 result += temp;
909 derivative_ -= o_.shapedSwapPaymentTimes_[i] * temp;
910 }
911 result *= Rs_;
912 derivative_ *= Rs_;
913 Real temp = o_.swapPaymentDiscounts_.back()
914 * std::exp(-o_.shapedSwapPaymentTimes_.back()*x);
915
916 result += temp-o_.discountAtStart_;
917 derivative_ -= o_.shapedSwapPaymentTimes_.back()*temp;
918 return result;
919 }
920
922 return derivative_;
923 }
924
926 Rs_ = x;
927 }
928
930 const Real x(s-swapStartTime_);
931 Real meanReversion = meanReversion_->value();
932 if(meanReversion>0) {
933 return (1.-std::exp(-meanReversion*x))/meanReversion;
934 }
935 else {
936 return x;
937 }
938 }
939
941
942 if(Rs!=tmpRs_){
943 Real initialGuess, N=0, D=0;
944 for(Size i=0; i<accruals_.size(); i++) {
945 N+=accruals_[i]*swapPaymentDiscounts_[i];
946 D+=accruals_[i]*swapPaymentDiscounts_[i]*shapedSwapPaymentTimes_[i];
947 }
948 N *= Rs;
949 D *= Rs;
950 N += accruals_.back() * swapPaymentDiscounts_.back()
951 - objectiveFunction_->gFunctionWithShifts().discountAtStart_;
952 D += accruals_.back() * swapPaymentDiscounts_.back()*
953 shapedSwapPaymentTimes_.back();
954 initialGuess = N/D;
955
956 objectiveFunction_->setSwapRateValue(Rs);
957 Newton solver;
958 solver.setMaxEvaluations(1000);
959
960 // these boundaries migth not be big enough if the volatility
961 // of big swap rate values is too high . In this case the G function
962 // is not even integrable, so better to fix the vol than increasing
963 // these values
964 const Real lower = -20, upper = 20.;
965
966 try {
967 calibratedShift_ = solver.solve(*objectiveFunction_, accuracy_,
968 std::max( std::min(initialGuess, upper*.99), lower*.99),
969 lower, upper);
970 } catch (std::exception& e) {
971 QL_FAIL("meanReversion: " << meanReversion_->value() <<
972 ", swapRateValue: " << swapRateValue_ <<
973 ", swapStartTime: " << swapStartTime_ <<
974 ", shapedPaymentTime: " << shapedPaymentTime_ <<
975 "\n error message: " << e.what());
976 }
977 tmpRs_=Rs;
978 }
979 return calibratedShift_;
980 }
981
982 ext::shared_ptr<GFunction> GFunctionFactory::newGFunctionWithShifts(const CmsCoupon& coupon,
983 const Handle<Quote>& meanReversion) {
984 return ext::shared_ptr<GFunction>(new GFunctionWithShifts(coupon, meanReversion));
985 }
986
987}
Black formula.
AnalyticHaganPricer(const Handle< SwaptionVolatilityStructure > &swaptionVol, GFunctionFactory::YieldCurveModel modelOfYieldCurve, const Handle< Quote > &meanReversion)
Real optionletPrice(Option::Type optionType, Real strike) const override
Real swapletPrice() const override
CMS coupon class.
Definition: cmscoupon.hpp:39
const ext::shared_ptr< SwapIndex > & swapIndex() const
Definition: cmscoupon.hpp:56
base pricer for vanilla CMS coupons
Handle< SwaptionVolatilityStructure > swaptionVolatility() const
Date date() const override
Definition: coupon.hpp:53
Time accrualPeriod() const
accrual period as fraction of year
Definition: coupon.cpp:44
Cumulative normal distribution function.
Concrete date class.
Definition: date.hpp:125
day counter class
Definition: daycounter.hpp:44
Time yearFraction(const Date &, const Date &, const Date &refPeriodStart=Date(), const Date &refPeriodEnd=Date()) const
Returns the period between two dates as a fraction of year.
Definition: daycounter.hpp:128
base floating-rate coupon class
virtual Date fixingDate() const
fixing date
Real gearing() const
index gearing, i.e. multiplicative coefficient for the index
Spread spread() const
spread paid over the fixing of the underlying index
Two-additive-factor gaussian model class.
Definition: g2.hpp:56
GFunctionWithShifts(const CmsCoupon &coupon, Handle< Quote > meanReversion)
ext::shared_ptr< ObjectiveFunction > objectiveFunction_
static ext::shared_ptr< GFunction > newGFunctionExactYield(const CmsCoupon &coupon)
static ext::shared_ptr< GFunction > newGFunctionStandard(Size q, Real delta, Size swapLength)
static ext::shared_ptr< GFunction > newGFunctionWithShifts(const CmsCoupon &coupon, const Handle< Quote > &meanReversion)
Integral of a 1-dimensional function using the Gauss-Kronrod methods.
Integral of a 1-dimensional function using the Gauss-Kronrod methods.
CMS-coupon pricer.
ext::shared_ptr< GFunction > gFunction_
const CmsCoupon * coupon_
Real capletPrice(Rate effectiveCap) const override
Rate floorletRate(Rate effectiveFloor) const override
Handle< Quote > meanReversion_
void initialize(const FloatingRateCoupon &coupon) override
virtual Real optionletPrice(Option::Type optionType, Real strike) const =0
ext::shared_ptr< VanillaOptionPricer > vanillaOptionPricer_
GFunctionFactory::YieldCurveModel modelOfYieldCurve_
ext::shared_ptr< YieldTermStructure > rateCurve_
Real swapletPrice() const override=0
HaganPricer(const Handle< SwaptionVolatilityStructure > &swaptionVol, GFunctionFactory::YieldCurveModel modelOfYieldCurve, Handle< Quote > meanReversion)
Rate swapletRate() const override
Real floorletPrice(Rate effectiveFloor) const override
Real meanReversion() const override
Rate capletRate(Rate effectiveCap) const override
Shared handle to an observable.
Definition: handle.hpp:41
virtual bool integrationSuccess() const
Definition: integral.cpp:71
ext::shared_ptr< SwaptionVolatilityStructure > volatilityStructure_
ext::shared_ptr< SmileSection > smile_
MarketQuotedOptionPricer(Rate forwardValue, Date expiryDate, const Period &swapTenor, const ext::shared_ptr< SwaptionVolatilityStructure > &volatilityStructure)
Real operator()(Real strike, Option::Type optionType, Real deflator) const override
Newton 1-D solver
Definition: newton.hpp:40
ConundrumIntegrand(ext::shared_ptr< VanillaOptionPricer > o, const ext::shared_ptr< YieldTermStructure > &rateCurve, ext::shared_ptr< GFunction > gFunction, Date fixingDate, Date paymentDate, Real annuity, Real forwardValue, Real strike, Option::Type optionType)
Real optionletPrice(Option::Type optionType, Rate strike) const override
Real resetLowerLimit(Real stdDeviationsForLowerLimit) const
Real resetUpperLimit(Real stdDeviationsForUpperLimit) const
NumericHaganPricer(const Handle< SwaptionVolatilityStructure > &swaptionVol, GFunctionFactory::YieldCurveModel modelOfYieldCurve, const Handle< Quote > &meanReversion, Rate lowerLimit=0.0, Rate upperLimit=1.0, Real precision=1.0e-6, Real hardUpperLimit=QL_MAX_REAL)
Real integrate(Real a, Real b, const ConundrumIntegrand &Integrand) const
Real swapletPrice() const override
Real refineIntegration(Real integralValue, const ConundrumIntegrand &integrand) const
std::pair< iterator, bool > registerWith(const ext::shared_ptr< Observable > &)
Definition: observable.hpp:228
Integer length() const
Definition: period.hpp:50
Payment schedule.
Definition: schedule.hpp:40
const Date & startDate() const
Definition: schedule.hpp:180
const Date & date(Size i) const
Definition: schedule.hpp:160
DateProxy & evaluationDate()
the date at which pricing is to be performed.
Definition: settings.hpp:147
market element returning a stored value
Definition: simplequote.hpp:33
static Settings & instance()
access to the unique instance
Definition: singleton.hpp:104
void setMaxEvaluations(Size evaluations)
Definition: solver1d.hpp:238
Real solve(const F &f, Real accuracy, Real guess, Real step) const
Definition: solver1d.hpp:84
CMS coupon.
Real a_
Real width_
ext::function< Real(Real)> f_
std::vector< Real > abscissas
std::vector< Real > functionValues
Size k_
CMS-coupon pricer.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
ext::function< Real(Real)> b
LinearInterpolation variance
Maps function, bind and cref to either the boost or std implementation.
Real Time
continuous quantity with 1-year units
Definition: types.hpp:62
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
Real Spread
spreads on interest rates
Definition: types.hpp:74
Real Rate
interest rates
Definition: types.hpp:70
std::size_t Size
size of a container
Definition: types.hpp:58
base class for interest rate indexes
Integral of a 1-dimensional function using the Gauss-Kronrod method.
Definition: any.hpp:35
Real bachelierBlackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount)
void swap(Array &v, Array &w) noexcept
Definition: array.hpp:903
Real blackFormula(Option::Type optionType, Real strike, Real forward, Real stdDev, Real discount, Real displacement)
bool close_enough(const Quantity &m1, const Quantity &m2, Size n)
Definition: quantity.cpp:182
std::vector< ext::shared_ptr< CashFlow > > Leg
Sequence of cash-flows.
Definition: cashflow.hpp:78
STL namespace.
Newton 1-D solver.
normal, cumulative and inverse cumulative distributions
ext::shared_ptr< YieldTermStructure > q
date schedule
simple quote class
Smile section base class.
swap-rate indexes
Simple fixed-rate vs Libor swap.
Interest-rate term structure.