1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4   Copyright (C) 2014, 2016 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 
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 lineartsrpricer.cpp
22 */
23 
24 #include <ql/cashflows/lineartsrpricer.hpp>
25 #include <ql/cashflows/fixedratecoupon.hpp>
26 #include <ql/cashflows/iborcoupon.hpp>
27 #include <ql/cashflows/cmscoupon.hpp>
28 #include <ql/termstructures/yieldtermstructure.hpp>
29 #include <ql/quotes/simplequote.hpp>
30 #include <ql/indexes/iborindex.hpp>
31 #include <ql/time/schedule.hpp>
32 #include <ql/instruments/vanillaswap.hpp>
33 #include <ql/math/solvers1d/brent.hpp>
34 #include <ql/math/integrals/kronrodintegral.hpp>
35 #include <ql/pricingengines/blackformula.hpp>
36 #include <ql/termstructures/volatility/atmsmilesection.hpp>
37 
38 namespace QuantLib {
39 
40     class LinearTsrPricer::integrand_f {
41         const LinearTsrPricer* pricer;
42       public:
integrand_f(const LinearTsrPricer * pricer)43         explicit integrand_f(const LinearTsrPricer* pricer) : pricer(pricer) {}
operator ()(Real x) const44         Real operator()(Real x) const {
45             return pricer->integrand(x);
46         }
47     };
48 
49     const Real LinearTsrPricer::defaultLowerBound = 0.0001,
50              LinearTsrPricer::defaultUpperBound = 2.0000;
51 
LinearTsrPricer(const Handle<SwaptionVolatilityStructure> & swaptionVol,const Handle<Quote> & meanReversion,const Handle<YieldTermStructure> & couponDiscountCurve,const Settings & settings,const ext::shared_ptr<Integrator> & integrator)52     LinearTsrPricer::LinearTsrPricer(
53         const Handle<SwaptionVolatilityStructure> &swaptionVol,
54         const Handle<Quote> &meanReversion,
55         const Handle<YieldTermStructure> &couponDiscountCurve,
56         const Settings &settings,
57         const ext::shared_ptr<Integrator> &integrator)
58         : CmsCouponPricer(swaptionVol), meanReversion_(meanReversion),
59           couponDiscountCurve_(couponDiscountCurve), settings_(settings),
60           volDayCounter_(swaptionVol->dayCounter()), integrator_(integrator) {
61 
62         if (!couponDiscountCurve_.empty())
63             registerWith(couponDiscountCurve_);
64 
65         if (integrator_ == NULL)
66             integrator_ =
67                 ext::make_shared<GaussKronrodNonAdaptive>(1E-10, 5000, 1E-10);
68     }
69 
GsrG(const Date & d) const70     Real LinearTsrPricer::GsrG(const Date &d) const {
71 
72         Real yf = volDayCounter_.yearFraction(fixingDate_, d);
73         if (std::fabs(meanReversion_->value()) < 1.0E-4)
74             return yf;
75         else
76             return (1.0 - std::exp(-meanReversion_->value() * yf)) /
77                    meanReversion_->value();
78     }
79 
singularTerms(const Option::Type type,const Real strike) const80     Real LinearTsrPricer::singularTerms(const Option::Type type,
81                                         const Real strike) const {
82 
83         Real omega = (type == Option::Call ? 1.0 : -1.0);
84         Real s1 = std::max(omega * (swapRateValue_ - strike), 0.0) *
85                   (a_ * swapRateValue_ + b_);
86         Real s2 = (a_ * strike + b_) *
87                   smileSection_->optionPrice(strike, strike < swapRateValue_
88                                                          ? Option::Put
89                                                          : Option::Call);
90         return s1 + s2;
91     }
92 
integrand(const Real strike) const93     Real LinearTsrPricer::integrand(const Real strike) const {
94         return 2.0 * a_ * smileSection_->optionPrice(
95                               strike, strike < swapRateValue_ ? Option::Put
96                                                               : Option::Call);
97     }
98 
initialize(const FloatingRateCoupon & coupon)99     void LinearTsrPricer::initialize(const FloatingRateCoupon &coupon) {
100 
101         coupon_ = dynamic_cast<const CmsCoupon *>(&coupon);
102         QL_REQUIRE(coupon_, "CMS coupon needed");
103         gearing_ = coupon_->gearing();
104         spread_ = coupon_->spread();
105 
106         fixingDate_ = coupon_->fixingDate();
107         paymentDate_ = coupon_->date();
108         swapIndex_ = coupon_->swapIndex();
109 
110         forwardCurve_ = swapIndex_->forwardingTermStructure();
111         if (swapIndex_->exogenousDiscount())
112             discountCurve_ = swapIndex_->discountingTermStructure();
113         else
114             discountCurve_ = forwardCurve_;
115 
116         // if no coupon discount curve is given just use the discounting curve
117         // from the swap index. for rate calculation this curve cancels out in
118         // the computation, so e.g. the discounting swap engine will produce
119         // correct results, even if the couponDiscountCurve is not set here.
120         // only the price member function in this class will be dependent on the
121         // coupon discount curve.
122 
123         today_ = QuantLib::Settings::instance().evaluationDate();
124 
125         if (paymentDate_ > today_ && !couponDiscountCurve_.empty())
126             couponDiscountRatio_ =
127                 couponDiscountCurve_->discount(paymentDate_) /
128                 discountCurve_->discount(paymentDate_);
129         else
130             couponDiscountRatio_ = 1.;
131 
132         spreadLegValue_ = spread_ * coupon_->accrualPeriod() *
133                           discountCurve_->discount(paymentDate_) *
134                           couponDiscountRatio_;
135 
136         if (fixingDate_ > today_) {
137 
138             swapTenor_ = swapIndex_->tenor();
139             swap_ = swapIndex_->underlyingSwap(fixingDate_);
140 
141             swapRateValue_ = swap_->fairRate();
142             annuity_ = 1.0E4 * std::fabs(swap_->fixedLegBPS());
143 
144             ext::shared_ptr<SmileSection> sectionTmp =
145                 swaptionVolatility()->smileSection(fixingDate_, swapTenor_);
146 
147             adjustedLowerBound_ = settings_.lowerRateBound_;
148             adjustedUpperBound_ = settings_.upperRateBound_;
149 
150             if(sectionTmp->volatilityType() == Normal) {
151                 // adjust lower bound if it was not set explicitly
152                 if(settings_.defaultBounds_)
153                     adjustedLowerBound_ = std::min(adjustedLowerBound_, -adjustedUpperBound_);
154             } else {
155                 // adjust bounds by section's shift
156                 adjustedLowerBound_ -= sectionTmp->shift();
157                 adjustedUpperBound_ -= sectionTmp->shift();
158             }
159 
160             // if the section does not provide an atm level, we enhance it to
161             // have one, no need to exit with an exception ...
162 
163             if (sectionTmp->atmLevel() == Null<Real>())
164                 smileSection_ = ext::make_shared<AtmSmileSection>(
165                     sectionTmp, swapRateValue_);
166             else
167                 smileSection_ = sectionTmp;
168 
169             // compute linear model's parameters
170 
171             Real gx = 0.0, gy = 0.0;
172             for (Size i = 0; i < swap_->fixedLeg().size(); i++) {
173                 ext::shared_ptr<Coupon> c =
174                     ext::dynamic_pointer_cast<Coupon>(swap_->fixedLeg()[i]);
175                 Real yf = c->accrualPeriod();
176                 Date d = c->date();
177                 Real pv = yf * discountCurve_->discount(d);
178                 gx += pv * GsrG(d);
179                 gy += pv;
180             }
181 
182             Real gamma = gx / gy;
183             Date lastd = swap_->fixedLeg().back()->date();
184 
185             a_ = discountCurve_->discount(paymentDate_) *
186                  (gamma - GsrG(paymentDate_)) /
187                  (discountCurve_->discount(lastd) * GsrG(lastd) +
188                   swapRateValue_ * gy * gamma);
189 
190             b_ = discountCurve_->discount(paymentDate_) / gy -
191                  a_ * swapRateValue_;
192         }
193     }
194 
strikeFromVegaRatio(Real ratio,Option::Type optionType,Real referenceStrike) const195     Real LinearTsrPricer::strikeFromVegaRatio(Real ratio,
196                                               Option::Type optionType,
197                                               Real referenceStrike) const {
198 
199         Real a, b, min, max, k;
200         if (optionType == Option::Call) {
201             a = swapRateValue_;
202             min = referenceStrike;
203             b = max = k =
204                 std::min(smileSection_->maxStrike(), adjustedUpperBound_);
205         } else {
206             a = min = k =
207                 std::max(smileSection_->minStrike(), adjustedLowerBound_);
208             b = swapRateValue_;
209             max = referenceStrike;
210         }
211 
212         VegaRatioHelper h(&*smileSection_,
213                           smileSection_->vega(swapRateValue_) * ratio);
214         Brent solver;
215 
216         try {
217             k = solver.solve(h, 1.0E-5, (a + b) / 2.0, a, b);
218         }
219         catch (...) {
220             // use default value set above
221         }
222 
223         return std::min(std::max(k, min), max);
224     }
225 
strikeFromPrice(Real price,Option::Type optionType,Real referenceStrike) const226     Real LinearTsrPricer::strikeFromPrice(Real price, Option::Type optionType,
227                                           Real referenceStrike) const {
228 
229         Real a, b, min, max, k;
230         if (optionType == Option::Call) {
231             a = swapRateValue_;
232             min = referenceStrike;
233             b = max = k =
234                 std::min(smileSection_->maxStrike(), adjustedUpperBound_);
235         } else {
236             a = min = k =
237                 std::max(smileSection_->minStrike(), adjustedLowerBound_);
238             b = swapRateValue_;
239             max = referenceStrike;
240         }
241 
242         PriceHelper h(&*smileSection_, optionType, price);
243         Brent solver;
244 
245         try {
246             k = solver.solve(h, 1.0E-5, swapRateValue_, a, b);
247         }
248         catch (...) {
249             // use default value set above
250         }
251 
252         return std::min(std::max(k, min), max);
253     }
254 
optionletPrice(Option::Type optionType,Real strike) const255     Real LinearTsrPricer::optionletPrice(Option::Type optionType,
256                                          Real strike) const {
257 
258         if (optionType == Option::Call && strike >= adjustedUpperBound_)
259             return 0.0;
260         if (optionType == Option::Put && strike <= adjustedLowerBound_)
261             return 0.0;
262 
263         // determine lower or upper integration bound (depending on option type)
264 
265         Real lower = strike, upper = strike;
266 
267         switch (settings_.strategy_) {
268 
269         case Settings::RateBound: {
270             if (optionType == Option::Call)
271                 upper = adjustedUpperBound_;
272             else
273                 lower = adjustedLowerBound_;
274             break;
275         }
276 
277         case Settings::VegaRatio: {
278             // strikeFromVegaRatio ensures that returned strike is on the
279             // expected side of strike
280             Real bound =
281                 strikeFromVegaRatio(settings_.vegaRatio_, optionType, strike);
282             if (optionType == Option::Call)
283                 upper = std::min(bound, adjustedUpperBound_);
284             else
285                 lower = std::max(bound, adjustedLowerBound_);
286             break;
287         }
288 
289         case Settings::PriceThreshold: {
290             // strikeFromPrice ensures that returned strike is on the expected
291             // side of strike
292             Real bound =
293                 strikeFromPrice(settings_.vegaRatio_, optionType, strike);
294             if (optionType == Option::Call)
295                 upper = std::min(bound, adjustedUpperBound_);
296             else
297                 lower = std::max(bound, adjustedLowerBound_);
298             break;
299         }
300 
301         case Settings::BSStdDevs : {
302             Real atm = smileSection_->atmLevel();
303             Real atmVol = smileSection_->volatility(atm);
304             Real shift = smileSection_->shift();
305             Real lowerTmp, upperTmp;
306             if (smileSection_->volatilityType() == ShiftedLognormal) {
307                 upperTmp = (atm + shift) *
308                                std::exp(settings_.stdDevs_ * atmVol -
309                                         0.5 * atmVol * atmVol *
310                                             smileSection_->exerciseTime()) -
311                            shift;
312                 lowerTmp = (atm + shift) *
313                                std::exp(-settings_.stdDevs_ * atmVol -
314                                         0.5 * atmVol * atmVol *
315                                             smileSection_->exerciseTime()) -
316                            shift;
317             } else {
318                 Real tmp = settings_.stdDevs_ * atmVol *
319                            std::sqrt(smileSection_->exerciseTime());
320                 upperTmp = atm + tmp;
321                 lowerTmp = atm - tmp;
322             }
323             upper = std::min(upperTmp - shift, adjustedUpperBound_);
324             lower = std::max(lowerTmp - shift, adjustedLowerBound_);
325             break;
326         }
327 
328         default:
329             QL_FAIL("Unknown strategy (" << settings_.strategy_ << ")");
330         }
331 
332         // compute the relevant integral
333 
334         Real result = 0.0;
335         Real tmpBound;
336         if (upper > lower) {
337             tmpBound = std::min(upper, swapRateValue_);
338             if (tmpBound > lower) {
339                 result += (*integrator_)(integrand_f(this),
340                                          lower, tmpBound);
341             }
342             tmpBound = std::max(lower, swapRateValue_);
343             if (upper > tmpBound) {
344                 result += (*integrator_)(integrand_f(this),
345                                          tmpBound, upper);
346             }
347             result *= (optionType == Option::Call ? 1.0 : -1.0);
348         }
349 
350         result += singularTerms(optionType, strike);
351 
352         return annuity_ * result * couponDiscountRatio_ *
353                coupon_->accrualPeriod();
354     }
355 
meanReversion() const356     Real LinearTsrPricer::meanReversion() const { return meanReversion_->value(); }
357 
swapletRate() const358     Rate LinearTsrPricer::swapletRate() const {
359         return swapletPrice() /
360                (coupon_->accrualPeriod() *
361                 discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
362     }
363 
capletPrice(Rate effectiveCap) const364     Real LinearTsrPricer::capletPrice(Rate effectiveCap) const {
365         // caplet is equivalent to call option on fixing
366         if (fixingDate_ <= today_) {
367             // the fixing is determined
368             const Rate Rs = std::max(
369                 coupon_->swapIndex()->fixing(fixingDate_) - effectiveCap, 0.);
370             Rate price =
371                 (gearing_ * Rs) *
372                 (coupon_->accrualPeriod() *
373                  discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
374             return price;
375         } else {
376             Real capletPrice = optionletPrice(Option::Call, effectiveCap);
377             return gearing_ * capletPrice;
378         }
379     }
380 
capletRate(Rate effectiveCap) const381     Rate LinearTsrPricer::capletRate(Rate effectiveCap) const {
382         return capletPrice(effectiveCap) /
383                (coupon_->accrualPeriod() *
384                 discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
385     }
386 
floorletPrice(Rate effectiveFloor) const387     Real LinearTsrPricer::floorletPrice(Rate effectiveFloor) const {
388         // floorlet is equivalent to put option on fixing
389         if (fixingDate_ <= today_) {
390             // the fixing is determined
391             const Rate Rs = std::max(
392                 effectiveFloor - coupon_->swapIndex()->fixing(fixingDate_), 0.);
393             Rate price =
394                 (gearing_ * Rs) *
395                 (coupon_->accrualPeriod() *
396                  discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
397             return price;
398         } else {
399             Real floorletPrice = optionletPrice(Option::Put, effectiveFloor);
400             return gearing_ * floorletPrice;
401         }
402     }
403 
floorletRate(Rate effectiveFloor) const404     Rate LinearTsrPricer::floorletRate(Rate effectiveFloor) const {
405         return floorletPrice(effectiveFloor) /
406                (coupon_->accrualPeriod() *
407                 discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
408     }
409 
swapletPrice() const410     Real LinearTsrPricer::swapletPrice() const {
411         if (fixingDate_ <= today_) {
412             // the fixing is determined
413             const Rate Rs = coupon_->swapIndex()->fixing(fixingDate_);
414             Rate price =
415                 (gearing_ * Rs + spread_) *
416                 (coupon_->accrualPeriod() *
417                  discountCurve_->discount(paymentDate_) * couponDiscountRatio_);
418             return price;
419         } else {
420             Real atmCapletPrice = optionletPrice(Option::Call, swapRateValue_);
421             Real atmFloorletPrice = optionletPrice(Option::Put, swapRateValue_);
422             return gearing_ * (coupon_->accrualPeriod() *
423                                    discountCurve_->discount(paymentDate_) *
424                                    swapRateValue_ * couponDiscountRatio_ +
425                                atmCapletPrice - atmFloorletPrice) +
426                    spreadLegValue_;
427         }
428     }
429 }
430