[go: up one dir, main page]
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2019 SoftSolutions! S.r.l.
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/*! \file globalbootstrap.hpp
21 \brief global bootstrap, with additional restrictions
22 \ingroup termstructures
23*/
24
25#ifndef quantlib_global_bootstrap_hpp
26#define quantlib_global_bootstrap_hpp
27
28#include <ql/functional.hpp>
29#include <ql/math/interpolations/linearinterpolation.hpp>
30#include <ql/math/optimization/levenbergmarquardt.hpp>
31#include <ql/termstructures/bootstraperror.hpp>
32#include <ql/termstructures/bootstraphelper.hpp>
33#include <ql/utilities/dataformatters.hpp>
34#include <utility>
35
36namespace QuantLib {
37
38//! Global boostrapper, with additional restrictions
39template <class Curve> class GlobalBootstrap {
40 typedef typename Curve::traits_type Traits; // ZeroYield, Discount, ForwardRate
41 typedef typename Curve::interpolator_type Interpolator; // Linear, LogLinear, ...
42
43 public:
44 GlobalBootstrap(Real accuracy = Null<Real>());
45 /*! The set of (alive) additional dates is added to the interpolation grid. The set of additional dates must only
46 depend on the current global evaluation date. The additionalErrors functor must yield at least as many values
47 such that
48
49 number of (usual, alive) rate helpers + number of (alive) additional values >= number of data points - 1
50
51 (note that the data points contain t=0). These values are treated as additional error terms in the optimization,
52 the usual rate helpers return marketQuote - impliedQuote here. All error terms are equally weighted in the
53 optimisation.
54
55 The additional helpers are treated like the usual rate helpers, but no standard pillar dates are added for them.
56
57 WARNING: This class is known to work with Traits Discount, ZeroYield, Forward (i.e. the usual traits for IR curves
58 in QL), it might fail for other traits - check the usage of Traits::updateGuess(), Traits::guess(),
59 Traits::minValueAfter(), Traits::maxValueAfter() in this class against them.
60 */
61 GlobalBootstrap(std::vector<ext::shared_ptr<typename Traits::helper> > additionalHelpers,
62 ext::function<std::vector<Date>()> additionalDates,
63 ext::function<Array()> additionalErrors,
64 Real accuracy = Null<Real>());
65 void setup(Curve *ts);
66 void calculate() const;
67
68 private:
69 void initialize() const;
70 Curve *ts_;
71 Real accuracy_;
72 mutable std::vector<ext::shared_ptr<typename Traits::helper> > additionalHelpers_;
73 ext::function<std::vector<Date>()> additionalDates_;
74 ext::function<Array()> additionalErrors_;
75 mutable bool initialized_ = false, validCurve_ = false;
76 mutable Size firstHelper_, numberHelpers_;
77 mutable Size firstAdditionalHelper_, numberAdditionalHelpers_;
78 mutable Size firstAdditionalDate_, numberAdditionalDates_;
79 mutable std::vector<Real> lowerBounds_, upperBounds_;
80};
81
82// template definitions
83
84template <class Curve>
85GlobalBootstrap<Curve>::GlobalBootstrap(Real accuracy) : ts_(0), accuracy_(accuracy) {}
86
87template <class Curve>
88GlobalBootstrap<Curve>::GlobalBootstrap(
89 std::vector<ext::shared_ptr<typename Traits::helper> > additionalHelpers,
90 ext::function<std::vector<Date>()> additionalDates,
91 ext::function<Array()> additionalErrors,
92 Real accuracy)
93: ts_(nullptr), accuracy_(accuracy), additionalHelpers_(std::move(additionalHelpers)),
94 additionalDates_(std::move(additionalDates)), additionalErrors_(std::move(additionalErrors)) {}
95
96template <class Curve> void GlobalBootstrap<Curve>::setup(Curve *ts) {
97 ts_ = ts;
98 for (Size j = 0; j < ts_->instruments_.size(); ++j)
99 ts_->registerWithObservables(ts_->instruments_[j]);
100 for (Size j = 0; j < additionalHelpers_.size(); ++j)
101 ts_->registerWithObservables(additionalHelpers_[j]);
102
103 // do not initialize yet: instruments could be invalid here
104 // but valid later when bootstrapping is actually required
105}
106
107template <class Curve> void GlobalBootstrap<Curve>::initialize() const {
108
109 // ensure helpers are sorted
110 std::sort(ts_->instruments_.begin(), ts_->instruments_.end(), detail::BootstrapHelperSorter());
111 std::sort(additionalHelpers_.begin(), additionalHelpers_.end(), detail::BootstrapHelperSorter());
112
113 // skip expired helpers
114 Date firstDate = Traits::initialDate(ts_);
115
116 firstHelper_ = 0;
117 if (!ts_->instruments_.empty()) {
118 while (firstHelper_ < ts_->instruments_.size() && ts_->instruments_[firstHelper_]->pillarDate() <= firstDate)
119 ++firstHelper_;
120 }
121 numberHelpers_ = ts_->instruments_.size() - firstHelper_;
122
123 // skip expired additional helpers
124 firstAdditionalHelper_ = 0;
125 if (!additionalHelpers_.empty()) {
126 while (firstAdditionalHelper_ < additionalHelpers_.size() &&
127 additionalHelpers_[firstAdditionalHelper_]->pillarDate() <= firstDate)
128 ++firstAdditionalHelper_;
129 }
130 numberAdditionalHelpers_ = additionalHelpers_.size() - firstAdditionalHelper_;
131
132 // skip expired additional dates
133 std::vector<Date> additionalDates;
134 if (additionalDates_)
135 additionalDates = additionalDates_();
136 firstAdditionalDate_ = 0;
137 if (!additionalDates.empty()) {
138 while (firstAdditionalDate_ < additionalDates.size() && additionalDates[firstAdditionalDate_] <= firstDate)
139 ++firstAdditionalDate_;
140 }
141 numberAdditionalDates_ = additionalDates.size() - firstAdditionalDate_;
142
143 QL_REQUIRE(numberHelpers_ + numberAdditionalDates_ >= Interpolator::requiredPoints - 1,
144 "not enough alive instruments (" << numberHelpers_ << ") + additional dates (" << numberAdditionalDates_
145 << ") = " << numberHelpers_ + numberAdditionalDates_ << " provided, "
146 << Interpolator::requiredPoints - 1 << " required");
147
148 // calculate dates and times
149 std::vector<Date> &dates = ts_->dates_;
150 std::vector<Time> &times = ts_->times_;
151
152 // first populate the dates vector and make sure they are sorted and there are no duplicates
153 dates.clear();
154 dates.push_back(x: firstDate);
155 for (Size j = 0; j < numberHelpers_; ++j)
156 dates.push_back(ts_->instruments_[firstHelper_ + j]->pillarDate());
157 for (Size j = firstAdditionalDate_; j < numberAdditionalDates_; ++j)
158 dates.push_back(x: additionalDates[firstAdditionalDate_ + j]);
159 std::sort(first: dates.begin(), last: dates.end());
160 auto it = std::unique(first: dates.begin(), last: dates.end());
161 QL_REQUIRE(it == dates.end(), "duplicate dates among alive instruments and additional dates");
162
163 // build times vector
164 times.clear();
165 for (auto& date : dates)
166 times.push_back(ts_->timeFromReference(date));
167
168 // determine maxDate
169 Date maxDate = firstDate;
170 for (Size j = 0; j < numberHelpers_; ++j) {
171 maxDate = std::max(ts_->instruments_[firstHelper_ + j]->latestRelevantDate(), maxDate);
172 }
173 for (Size j = 0; j < numberAdditionalDates_; ++j) {
174 maxDate = std::max(a: additionalDates[firstAdditionalDate_ + j], b: maxDate);
175 }
176 ts_->maxDate_ = maxDate;
177
178 // set initial guess only if the current curve cannot be used as guess
179 if (!validCurve_ || ts_->data_.size() != dates.size()) {
180 // ts_->data_[0] is the only relevant item,
181 // but reasonable numbers might be needed for the whole data vector
182 // because, e.g., of interpolation's early checks
183 ts_->data_ = std::vector<Real>(dates.size(), Traits::initialValue(ts_));
184 }
185 initialized_ = true;
186}
187
188template <class Curve> void GlobalBootstrap<Curve>::calculate() const {
189
190 // we might have to call initialize even if the curve is initialized
191 // and not moving, just because helpers might be date relative and change
192 // with evaluation date change.
193 // anyway it makes little sense to use date relative helpers with a
194 // non-moving curve if the evaluation date changes
195 if (!initialized_ || ts_->moving_)
196 initialize();
197
198 // setup helpers
199 for (Size j = 0; j < numberHelpers_; ++j) {
200 const ext::shared_ptr<typename Traits::helper> &helper = ts_->instruments_[firstHelper_ + j];
201 // check for valid quote
202 QL_REQUIRE(helper->quote()->isValid(), io::ordinal(j + 1)
203 << " instrument (maturity: " << helper->maturityDate()
204 << ", pillar: " << helper->pillarDate() << ") has an invalid quote");
205 // don't try this at home!
206 // This call creates helpers, and removes "const".
207 // There is a significant interaction with observability.
208 helper->setTermStructure(const_cast<Curve *>(ts_));
209 }
210
211 // setup additional helpers
212 for (Size j = 0; j < numberAdditionalHelpers_; ++j) {
213 const ext::shared_ptr<typename Traits::helper> &helper = additionalHelpers_[firstAdditionalHelper_ + j];
214 QL_REQUIRE(helper->quote()->isValid(), io::ordinal(j + 1)
215 << " additional instrument (maturity: " << helper->maturityDate()
216 << ") has an invalid quote");
217 helper->setTermStructure(const_cast<Curve *>(ts_));
218 }
219
220 Real accuracy = accuracy_ != Null<Real>() ? accuracy_ : ts_->accuracy_;
221
222 // setup optimizer and EndCriteria
223 Real optEps = accuracy;
224 LevenbergMarquardt optimizer(optEps, optEps, optEps); // FIXME hardcoded tolerances
225 EndCriteria ec(1000, 10, optEps, optEps, optEps); // FIXME hardcoded values here as well
226
227 // setup interpolation
228 if (!validCurve_) {
229 ts_->interpolation_ =
230 ts_->interpolator_.interpolate(ts_->times_.begin(), ts_->times_.end(), ts_->data_.begin());
231 }
232
233 // determine bounds, we use an unconstrained optimisation transforming the free variables to [lowerBound,upperBound]
234 std::vector<Real> lowerBounds(numberHelpers_ + numberAdditionalDates_),
235 upperBounds(numberHelpers_ + numberAdditionalDates_);
236 for (Size i = 0; i < numberHelpers_ + numberAdditionalDates_; ++i) {
237 // just pass zero as the first alive helper, it's not used in the standard QL traits anyway
238 lowerBounds[i] = Traits::minValueAfter(i + 1, ts_, validCurve_, 0);
239 upperBounds[i] = Traits::maxValueAfter(i + 1, ts_, validCurve_, 0);
240 }
241
242 // setup cost function
243 class TargetFunction : public CostFunction {
244 public:
245 TargetFunction(const Size firstHelper,
246 const Size numberHelpers,
247 ext::function<Array()> additionalErrors,
248 Curve* ts,
249 std::vector<Real> lowerBounds,
250 std::vector<Real> upperBounds)
251 : firstHelper_(firstHelper), numberHelpers_(numberHelpers),
252 additionalErrors_(std::move(additionalErrors)), ts_(ts),
253 lowerBounds_(std::move(lowerBounds)), upperBounds_(std::move(upperBounds)) {}
254
255 Real transformDirect(const Real x, const Size i) const {
256 return (std::atan(x: x) + M_PI_2) / M_PI * (upperBounds_[i] - lowerBounds_[i]) + lowerBounds_[i];
257 }
258
259 Real transformInverse(const Real y, const Size i) const {
260 return std::tan(x: (y - lowerBounds_[i]) * M_PI / (upperBounds_[i] - lowerBounds_[i]) - M_PI_2);
261 }
262
263 Real value(const Array& x) const override {
264 Array v = values(x);
265 std::transform(v.begin(), v.end(), v.begin(), [](Real x) -> Real { return x*x; });
266 return std::sqrt(x: std::accumulate(first: v.begin(), last: v.end(), init: Real(0.0)) / static_cast<Real>(v.size()));
267 }
268
269 Array values(const Array& x) const override {
270 for (Size i = 0; i < x.size(); ++i) {
271 Traits::updateGuess(ts_->data_, transformDirect(x: x[i], i), i + 1);
272 }
273 ts_->interpolation_.update();
274 std::vector<Real> result(numberHelpers_);
275 for (Size i = 0; i < numberHelpers_; ++i) {
276 result[i] = ts_->instruments_[firstHelper_ + i]->quote()->value() -
277 ts_->instruments_[firstHelper_ + i]->impliedQuote();
278 }
279 if (additionalErrors_) {
280 Array tmp = additionalErrors_();
281 result.resize(new_size: numberHelpers_ + tmp.size());
282 for (Size i = 0; i < tmp.size(); ++i) {
283 result[numberHelpers_ + i] = tmp[i];
284 }
285 }
286 return Array(result.begin(), result.end());
287 }
288
289 private:
290 Size firstHelper_, numberHelpers_;
291 ext::function<Array()> additionalErrors_;
292 Curve *ts_;
293 const std::vector<Real> lowerBounds_, upperBounds_;
294 };
295 TargetFunction cost(firstHelper_, numberHelpers_, additionalErrors_, ts_, lowerBounds, upperBounds);
296
297 // setup guess
298 Array guess(numberHelpers_ + numberAdditionalDates_);
299 for (Size i = 0; i < guess.size(); ++i) {
300 // just pass zero as the first alive helper, it's not used in the standard QL traits anyway
301 guess[i] = cost.transformInverse(Traits::guess(i + 1, ts_, validCurve_, 0), i);
302 }
303
304 // setup problem
305 NoConstraint noConstraint;
306 Problem problem(cost, noConstraint, guess);
307
308 // run optimization
309 optimizer.minimize(P&: problem, endCriteria: ec);
310
311 // evaluate target function on best value found to ensure that data_ contains the optimal value
312 Real finalTargetError = cost.value(problem.currentValue());
313
314 // check final error
315 QL_REQUIRE(finalTargetError <= accuracy,
316 "global bootstrap failed, error is " << finalTargetError << ", accuracy is " << accuracy);
317
318 // set valid flag
319 validCurve_ = true;
320}
321
322} // namespace QuantLib
323
324#endif
325

source code of quantlib/ql/termstructures/globalbootstrap.hpp