1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008, 2011, 2015 Ferdinando Ametrano
5 Copyright (C) 2007 Chris Kenyon
6 Copyright (C) 2007 StatPro Italia srl
7 Copyright (C) 2015 Paolo Mazzocchi
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23/*! \file iterativebootstrap.hpp
24 \brief universal piecewise-term-structure boostrapper.
25*/
26
27#ifndef quantlib_iterative_bootstrap_hpp
28#define quantlib_iterative_bootstrap_hpp
29
30#include <ql/termstructures/bootstraphelper.hpp>
31#include <ql/termstructures/bootstraperror.hpp>
32#include <ql/math/interpolations/linearinterpolation.hpp>
33#include <ql/math/solvers1d/finitedifferencenewtonsafe.hpp>
34#include <ql/math/solvers1d/brent.hpp>
35#include <ql/utilities/dataformatters.hpp>
36
37namespace QuantLib {
38
39namespace detail {
40
41 /*! If \c dontThrow is \c true in IterativeBootstrap and on a given pillar the bootstrap fails when
42 searching for a helper root between \c xMin and \c xMax, we use this function to return the value that
43 gives the minimum absolute helper error in the interval between \c xMin and \c xMax inclusive.
44 */
45 template <class Curve>
46 Real dontThrowFallback(const BootstrapError<Curve>& error,
47 Real xMin, Real xMax, Size steps) {
48
49 QL_REQUIRE(xMin < xMax, "Expected xMin to be less than xMax");
50
51 // Set the initial value of the result to xMin and store the absolute bootstrap error at xMin
52 Real result = xMin;
53 Real absError = std::abs(error(xMin));
54 Real minError = absError;
55
56 // Step out to xMax
57 Real stepSize = (xMax - xMin) / steps;
58 for (Size i = 0; i < steps; i++) {
59
60 // Get absolute bootstrap error at updated x value
61 xMin += stepSize;
62 absError = std::abs(error(xMin));
63
64 // If this absolute bootstrap error is less than the minimum, update result and minError
65 if (absError < minError) {
66 result = xMin;
67 minError = absError;
68 }
69 }
70
71 return result;
72 }
73
74}
75
76 //! Universal piecewise-term-structure boostrapper.
77 template <class Curve>
78 class IterativeBootstrap {
79 typedef typename Curve::traits_type Traits;
80 typedef typename Curve::interpolator_type Interpolator;
81 public:
82 /*! Constructor
83 \param accuracy Accuracy for the bootstrap stopping criterion. If it is set to
84 \c Null<Real>(), its value is taken from the termstructure's accuracy.
85 \param minValue Allow to override the initial minimum value coming from traits.
86 \param maxValue Allow to override the initial maximum value coming from traits.
87 \param maxAttempts Number of attempts on each iteration. A number greater than 1 implies retries.
88 \param maxFactor Factor for max value retry on each iteration if there is a failure.
89 \param minFactor Factor for min value retry on each iteration if there is a failure.
90 \param dontThrow If set to \c true, the bootstrap doesn't throw and returns a <em>fall back</em>
91 result.
92 \param dontThrowSteps If \p dontThrow is \c true, this gives the number of steps to use when searching
93 for a fallback curve pillar value that gives the minimum bootstrap helper error.
94 */
95 IterativeBootstrap(Real accuracy = Null<Real>(),
96 Real minValue = Null<Real>(),
97 Real maxValue = Null<Real>(),
98 Size maxAttempts = 1,
99 Real maxFactor = 2.0,
100 Real minFactor = 2.0,
101 bool dontThrow = false,
102 Size dontThrowSteps = 10,
103 Size maxEvaluations = MAX_FUNCTION_EVALUATIONS);
104 void setup(Curve* ts);
105 void calculate() const;
106 private:
107 void initialize() const;
108 Real accuracy_;
109 Real minValue_, maxValue_;
110 Size maxAttempts_;
111 Real maxFactor_;
112 Real minFactor_;
113 bool dontThrow_;
114 Size dontThrowSteps_;
115 Curve* ts_;
116 Size n_ = 0;
117 Brent firstSolver_;
118 FiniteDifferenceNewtonSafe solver_;
119 mutable bool initialized_ = false, validCurve_ = false, loopRequired_;
120 mutable Size firstAliveHelper_ = 0, alive_ = 0;
121 mutable std::vector<Real> previousData_;
122 mutable std::vector<ext::shared_ptr<BootstrapError<Curve> > > errors_;
123 };
124
125
126 // template definitions
127
128 template <class Curve>
129 IterativeBootstrap<Curve>::IterativeBootstrap(Real accuracy,
130 Real minValue,
131 Real maxValue,
132 Size maxAttempts,
133 Real maxFactor,
134 Real minFactor,
135 bool dontThrow,
136 Size dontThrowSteps,
137 Size maxEvaluations)
138 : accuracy_(accuracy), minValue_(minValue), maxValue_(maxValue), maxAttempts_(maxAttempts),
139 maxFactor_(maxFactor), minFactor_(minFactor), dontThrow_(dontThrow),
140 dontThrowSteps_(dontThrowSteps), ts_(nullptr), loopRequired_(Interpolator::global) {
141 QL_REQUIRE(maxFactor_ >= 1.0, "Expected that maxFactor would be at least 1.0 but got " << maxFactor_);
142 QL_REQUIRE(minFactor_ >= 1.0, "Expected that minFactor would be at least 1.0 but got " << minFactor_);
143 firstSolver_.setMaxEvaluations(maxEvaluations);
144 solver_.setMaxEvaluations(maxEvaluations);
145 }
146
147 template <class Curve>
148 void IterativeBootstrap<Curve>::setup(Curve* ts) {
149 ts_ = ts;
150 n_ = ts_->instruments_.size();
151 QL_REQUIRE(n_ > 0, "no bootstrap helpers given");
152 for (Size j=0; j<n_; ++j)
153 ts_->registerWithObservables(ts_->instruments_[j]);
154
155 // do not initialize yet: instruments could be invalid here
156 // but valid later when bootstrapping is actually required
157 }
158
159 template <class Curve>
160 void IterativeBootstrap<Curve>::initialize() const {
161 // ensure helpers are sorted
162 std::sort(ts_->instruments_.begin(), ts_->instruments_.end(),
163 detail::BootstrapHelperSorter());
164 // skip expired helpers
165 Date firstDate = Traits::initialDate(ts_);
166 QL_REQUIRE(ts_->instruments_[n_-1]->pillarDate()>firstDate,
167 "all instruments expired");
168 firstAliveHelper_ = 0;
169 while (ts_->instruments_[firstAliveHelper_]->pillarDate() <= firstDate)
170 ++firstAliveHelper_;
171 alive_ = n_-firstAliveHelper_;
172 Size nodes = alive_+1;
173 QL_REQUIRE(nodes >= Interpolator::requiredPoints,
174 "not enough alive instruments: " << alive_ <<
175 " provided, " << Interpolator::requiredPoints-1 <<
176 " required");
177
178 // calculate dates and times, create errors_
179 std::vector<Date>& dates = ts_->dates_;
180 std::vector<Time>& times = ts_->times_;
181 dates.resize(new_size: alive_+1);
182 times.resize(new_size: alive_+1);
183 errors_.resize(alive_+1);
184 dates[0] = firstDate;
185 times[0] = ts_->timeFromReference(dates[0]);
186
187 Date latestRelevantDate, maxDate = firstDate;
188 // pillar counter: i
189 // helper counter: j
190 for (Size i=1, j=firstAliveHelper_; j<n_; ++i, ++j) {
191 const ext::shared_ptr<typename Traits::helper>& helper =
192 ts_->instruments_[j];
193 dates[i] = helper->pillarDate();
194 times[i] = ts_->timeFromReference(dates[i]);
195 // check for duplicated pillars
196 QL_REQUIRE(dates[i-1]!=dates[i],
197 "more than one instrument with pillar " << dates[i]);
198
199 latestRelevantDate = helper->latestRelevantDate();
200 // check that the helper is really extending the curve, i.e. that
201 // pillar-sorted helpers are also sorted by latestRelevantDate
202 QL_REQUIRE(latestRelevantDate > maxDate,
203 io::ordinal(j+1) << " instrument (pillar: " <<
204 dates[i] << ") has latestRelevantDate (" <<
205 latestRelevantDate << ") before or equal to "
206 "previous instrument's latestRelevantDate (" <<
207 maxDate << ")");
208 maxDate = latestRelevantDate;
209
210 // when a pillar date is different from the last relevant date the
211 // convergence loop is required even if the Interpolator is local
212 if (dates[i] != latestRelevantDate)
213 loopRequired_ = true;
214
215 errors_[i] = ext::shared_ptr<BootstrapError<Curve> >(new
216 BootstrapError<Curve>(ts_, helper, i));
217 }
218 ts_->maxDate_ = maxDate;
219
220 // set initial guess only if the current curve cannot be used as guess
221 if (!validCurve_ || ts_->data_.size()!=alive_+1) {
222 // ts_->data_[0] is the only relevant item,
223 // but reasonable numbers might be needed for the whole data vector
224 // because, e.g., of interpolation's early checks
225 ts_->data_ = std::vector<Real>(alive_+1, Traits::initialValue(ts_));
226 previousData_.resize(new_size: alive_+1);
227 validCurve_ = false;
228 }
229 initialized_ = true;
230 }
231
232 template <class Curve>
233 void IterativeBootstrap<Curve>::calculate() const {
234
235 // we might have to call initialize even if the curve is initialized
236 // and not moving, just because helpers might be date relative and change
237 // with evaluation date change.
238 // anyway it makes little sense to use date relative helpers with a
239 // non-moving curve if the evaluation date changes
240 if (!initialized_ || ts_->moving_)
241 initialize();
242
243 // setup helpers
244 for (Size j=firstAliveHelper_; j<n_; ++j) {
245 const ext::shared_ptr<typename Traits::helper>& helper =
246 ts_->instruments_[j];
247 // check for valid quote
248 QL_REQUIRE(helper->quote()->isValid(),
249 io::ordinal(j + 1) << " instrument (maturity: " <<
250 helper->maturityDate() << ", pillar: " <<
251 helper->pillarDate() << ") has an invalid quote");
252 // don't try this at home!
253 // This call creates helpers, and removes "const".
254 // There is a significant interaction with observability.
255 helper->setTermStructure(const_cast<Curve*>(ts_));
256 }
257
258 const std::vector<Time>& times = ts_->times_;
259 const std::vector<Real>& data = ts_->data_;
260 Real accuracy = accuracy_ != Null<Real>() ? accuracy_ : ts_->accuracy_;
261
262 Size maxIterations = Traits::maxIterations()-1;
263
264 // there might be a valid curve state to use as guess
265 bool validData = validCurve_;
266
267 for (Size iteration=0; ; ++iteration) {
268 previousData_ = ts_->data_;
269
270 // Store min value and max value at each pillar so that we can expand search if necessary.
271 std::vector<Real> minValues(alive_+1, Null<Real>());
272 std::vector<Real> maxValues(alive_+1, Null<Real>());
273 std::vector<Size> attempts(alive_+1, 1);
274
275 for (Size i=1; i<=alive_; ++i) { // pillar loop
276
277 // shorter aliases for readability and to avoid duplication
278 Real& min = minValues[i];
279 Real& max = maxValues[i];
280
281 // bracket root and calculate guess
282 if (min == Null<Real>()) {
283 // First attempt; we take min and max either from
284 // explicit constructor parameter or from traits
285 min = (minValue_ != Null<Real>() ? minValue_ :
286 Traits::minValueAfter(i, ts_, validData, firstAliveHelper_));
287 max = (maxValue_ != Null<Real>() ? maxValue_ :
288 Traits::maxValueAfter(i, ts_, validData, firstAliveHelper_));
289 } else {
290 // Extending a previous attempt. A negative min
291 // is enlarged; a positive one is shrunk towards 0.
292 min = (min < 0.0 ? Real(min * minFactor_) : Real(min / minFactor_));
293 // The opposite holds for the max.
294 max = (max > 0.0 ? Real(max * maxFactor_) : Real(max / maxFactor_));
295 }
296 Real guess = Traits::guess(i, ts_, validData, firstAliveHelper_);
297
298 // adjust guess if needed
299 if (guess >= max)
300 guess = max - (max - min) / 5.0;
301 else if (guess <= min)
302 guess = min + (max - min) / 5.0;
303
304 // extend interpolation if needed
305 if (!validData) {
306 try { // extend interpolation a point at a time
307 // including the pillar to be boostrapped
308 ts_->interpolation_ = ts_->interpolator_.interpolate(
309 times.begin(), times.begin()+i+1, data.begin());
310 } catch (...) {
311 if (!Interpolator::global)
312 throw; // no chance to fix it in a later iteration
313
314 // otherwise use Linear while the target
315 // interpolation is not usable yet
316 ts_->interpolation_ = Linear().interpolate(
317 xBegin: times.begin(), xEnd: times.begin()+i+1, yBegin: data.begin());
318 }
319 ts_->interpolation_.update();
320 }
321
322 try {
323 if (validData)
324 solver_.solve(*errors_[i], accuracy, guess, min, max);
325 else
326 firstSolver_.solve(*errors_[i], accuracy, guess, min, max);
327 } catch (std::exception &e) {
328 if (validCurve_) {
329 // the previous curve state might have been a
330 // bad guess, so we retry without using it.
331 // This would be tricky to do here (we're
332 // inside multiple nested for loops, we need
333 // to re-initialize...), so we invalidate the
334 // curve, make a recursive call and then exit.
335 validCurve_ = initialized_ = false;
336 calculate();
337 return;
338 }
339
340 // If we have more attempts left on this iteration, try again. Note that the max and min
341 // bounds will be widened on the retry.
342 if (attempts[i] < maxAttempts_) {
343 attempts[i]++;
344 i--;
345 continue;
346 }
347
348 if (dontThrow_) {
349 // Use the fallback value
350 ts_->data_[i] = detail::dontThrowFallback(*errors_[i], min, max, dontThrowSteps_);
351
352 // Remember to update the interpolation. If we don't and we are on the last "i", we will still
353 // have the last attempted value in the solver being used in ts_->interpolation_.
354 ts_->interpolation_.update();
355 } else {
356 QL_FAIL(io::ordinal(iteration + 1) << " iteration: failed "
357 "at " << io::ordinal(i) << " alive instrument, "
358 "pillar " << errors_[i]->helper()->pillarDate() <<
359 ", maturity " << errors_[i]->helper()->maturityDate() <<
360 ", reference date " << ts_->dates_[0] <<
361 ": " << e.what());
362 }
363 }
364 }
365
366 if (!loopRequired_)
367 break;
368
369 // exit condition
370 Real change = std::fabs(x: data[1]-previousData_[1]);
371 for (Size i=2; i<=alive_; ++i)
372 change = std::max(a: change, b: std::fabs(x: data[i]-previousData_[i]));
373 if (change<=accuracy) // convergence reached
374 break;
375
376 // If we hit the max number of iterations and dontThrow is true, just use what we have
377 if (iteration == maxIterations) {
378 if (dontThrow_) {
379 break;
380 } else {
381 QL_FAIL("convergence not reached after " << iteration <<
382 " iterations; last improvement " << change <<
383 ", required accuracy " << accuracy);
384 }
385 }
386
387 validData = true;
388 }
389 validCurve_ = true;
390 }
391
392}
393
394#endif
395

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