1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2008 Simon Ibbotson
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 localbootstrap.hpp
21 \brief localised-term-structure bootstrapper for most curve types.
22*/
23
24#ifndef quantlib_local_bootstrap_hpp
25#define quantlib_local_bootstrap_hpp
26
27#include <ql/termstructures/bootstraphelper.hpp>
28#include <ql/math/optimization/costfunction.hpp>
29#include <ql/math/optimization/constraint.hpp>
30#include <ql/math/optimization/armijo.hpp>
31#include <ql/math/optimization/levenbergmarquardt.hpp>
32#include <ql/math/optimization/problem.hpp>
33#include <ql/utilities/dataformatters.hpp>
34#include <ql/shared_ptr.hpp>
35
36namespace QuantLib {
37
38 // penalty function class for solving using a multi-dimensional solver
39 template <class Curve>
40 class PenaltyFunction : public CostFunction {
41 typedef typename Curve::traits_type Traits;
42 typedef typename Traits::helper helper;
43 typedef
44 typename std::vector< ext::shared_ptr<helper> >::const_iterator
45 helper_iterator;
46 public:
47 PenaltyFunction(Curve* curve,
48 Size initialIndex,
49 helper_iterator rateHelpersStart,
50 helper_iterator rateHelpersEnd)
51 : curve_(curve), initialIndex_(initialIndex),
52 localisation_(std::distance(rateHelpersStart, rateHelpersEnd)),
53 rateHelpersStart_(rateHelpersStart), rateHelpersEnd_(rateHelpersEnd) {}
54
55 Real value(const Array& x) const override;
56 Array values(const Array& x) const override;
57
58 private:
59 Curve* curve_;
60 Size initialIndex_;
61 Size localisation_;
62 helper_iterator rateHelpersStart_;
63 helper_iterator rateHelpersEnd_;
64 };
65
66
67 //! Localised-term-structure bootstrapper for most curve types.
68 /*! This algorithm enables a localised fitting for non-local
69 interpolation methods.
70
71 As in the similar class (IterativeBootstrap) the input term
72 structure is solved on a number of market instruments which
73 are passed as a vector of handles to BootstrapHelper
74 instances. Their maturities mark the boundaries of the
75 interpolated segments.
76
77 Unlike the IterativeBootstrap class, the solution for each
78 interpolated segment is derived using a local
79 approximation. This restricts the risk profile s.t. the risk
80 is localised. Therefore, we obtain a local IR risk profile
81 whilst using a smoother interpolation method. Particularly
82 good for the convex-monotone spline method.
83 */
84 template <class Curve>
85 class LocalBootstrap {
86 typedef typename Curve::traits_type Traits;
87 typedef typename Curve::interpolator_type Interpolator;
88 public:
89 LocalBootstrap(Size localisation = 2,
90 bool forcePositive = true,
91 Real accuracy = Null<Real>());
92 void setup(Curve* ts);
93 void calculate() const;
94
95 private:
96 mutable bool validCurve_ = false;
97 Curve* ts_;
98 Size localisation_;
99 bool forcePositive_;
100 Real accuracy_;
101 };
102
103
104
105 // template definitions
106
107 template <class Curve>
108 LocalBootstrap<Curve>::LocalBootstrap(Size localisation, bool forcePositive, Real accuracy)
109 : ts_(nullptr), localisation_(localisation), forcePositive_(forcePositive),
110 accuracy_(accuracy) {}
111
112 template <class Curve>
113 void LocalBootstrap<Curve>::setup(Curve* ts) {
114
115 ts_ = ts;
116
117 Size n = ts_->instruments_.size();
118 QL_REQUIRE(n >= Interpolator::requiredPoints,
119 "not enough instruments: " << n << " provided, " <<
120 Interpolator::requiredPoints << " required");
121
122 QL_REQUIRE(n > localisation_,
123 "not enough instruments: " << n << " provided, " <<
124 localisation_ << " required.");
125
126 for (Size i=0; i<n; ++i){
127 ts_->registerWithObservables(ts_->instruments_[i]);
128 }
129 }
130
131 template <class Curve>
132 void LocalBootstrap<Curve>::calculate() const {
133
134 validCurve_ = false;
135 Size nInsts = ts_->instruments_.size();
136
137 // ensure rate helpers are sorted
138 std::sort(ts_->instruments_.begin(), ts_->instruments_.end(),
139 detail::BootstrapHelperSorter());
140
141 // check that there is no instruments with the same maturity
142 for (Size i=1; i<nInsts; ++i) {
143 Date m1 = ts_->instruments_[i-1]->pillarDate(),
144 m2 = ts_->instruments_[i]->pillarDate();
145 QL_REQUIRE(m1 != m2,
146 "two instruments have the same pillar date ("<<m1<<")");
147 }
148
149 // check that there is no instruments with invalid quote
150 for (Size i=0; i<nInsts; ++i)
151 QL_REQUIRE(ts_->instruments_[i]->quote()->isValid(),
152 io::ordinal(i+1) << " instrument (maturity: " <<
153 ts_->instruments_[i]->maturityDate() << ", pillar: " <<
154 ts_->instruments_[i]->pillarDate() <<
155 ") has an invalid quote");
156
157 // setup instruments
158 for (Size i=0; i<nInsts; ++i) {
159 // don't try this at home!
160 // This call creates instruments, and removes "const".
161 // There is a significant interaction with observability.
162 ts_->instruments_[i]->setTermStructure(const_cast<Curve*>(ts_));
163 }
164 // set initial guess only if the current curve cannot be used as guess
165 if (validCurve_)
166 QL_ENSURE(ts_->data_.size() == nInsts+1,
167 "dimension mismatch: expected " << nInsts+1 <<
168 ", actual " << ts_->data_.size());
169 else {
170 ts_->data_ = std::vector<Rate>(nInsts+1);
171 ts_->data_[0] = Traits::initialValue(ts_);
172 }
173
174 // calculate dates and times
175 ts_->dates_ = std::vector<Date>(nInsts+1);
176 ts_->times_ = std::vector<Time>(nInsts+1);
177 ts_->dates_[0] = Traits::initialDate(ts_);
178 ts_->times_[0] = ts_->timeFromReference(ts_->dates_[0]);
179 for (Size i=0; i<nInsts; ++i) {
180 ts_->dates_[i+1] = ts_->instruments_[i]->pillarDate();
181 ts_->times_[i+1] = ts_->timeFromReference(ts_->dates_[i+1]);
182 if (!validCurve_)
183 ts_->data_[i+1] = ts_->data_[i];
184 }
185
186 Real accuracy = accuracy_ != Null<Real>() ? accuracy_ : ts_->accuracy_;
187
188 LevenbergMarquardt solver(accuracy,
189 accuracy,
190 accuracy);
191 EndCriteria endCriteria(100, 10, 0.00, accuracy, 0.00);
192 PositiveConstraint posConstraint;
193 NoConstraint noConstraint;
194 Constraint& solverConstraint = forcePositive_ ?
195 static_cast<Constraint&>(posConstraint) :
196 static_cast<Constraint&>(noConstraint);
197
198 // now start the bootstrapping.
199 Size iInst = localisation_-1;
200
201 Size dataAdjust = Curve::interpolator_type::dataSizeAdjustment;
202
203 do {
204 Size initialDataPt = iInst+1-localisation_+dataAdjust;
205 Array startArray(localisation_+1-dataAdjust);
206 for (Size j = 0; j < startArray.size()-1; ++j)
207 startArray[j] = ts_->data_[initialDataPt+j];
208
209 // here we are extending the interpolation a point at a
210 // time... but the local interpolator can make an
211 // approximation for the final localisation period.
212 // e.g. if the localisation is 2, then the first section
213 // of the curve will be solved using the first 2
214 // instruments... with the local interpolator making
215 // suitable boundary conditions.
216 ts_->interpolation_ =
217 ts_->interpolator_.localInterpolate(
218 ts_->times_.begin(),
219 ts_->times_.begin()+(iInst + 2),
220 ts_->data_.begin(),
221 localisation_,
222 ts_->interpolation_,
223 nInsts+1);
224
225 if (iInst >= localisation_) {
226 startArray[localisation_-dataAdjust] =
227 Traits::guess(iInst, ts_, false, 0); // ?
228 } else {
229 startArray[localisation_-dataAdjust] = ts_->data_[0];
230 }
231
232 PenaltyFunction<Curve> currentCost(
233 ts_,
234 initialDataPt,
235 ts_->instruments_.begin() + ((iInst+1) - localisation_),
236 ts_->instruments_.begin() + (iInst+1));
237
238 Problem toSolve(currentCost, solverConstraint, startArray);
239
240 EndCriteria::Type endType = solver.minimize(P&: toSolve, endCriteria);
241
242 // check the end criteria
243 QL_REQUIRE(endType == EndCriteria::StationaryFunctionAccuracy ||
244 endType == EndCriteria::StationaryFunctionValue,
245 "Unable to strip yieldcurve to required accuracy " );
246 ++iInst;
247 } while ( iInst < nInsts );
248 validCurve_ = true;
249 }
250
251
252 template <class Curve>
253 Real PenaltyFunction<Curve>::value(const Array& x) const {
254 Size i = initialIndex_;
255 Array::const_iterator guessIt = x.begin();
256 while (guessIt != x.end()) {
257 Traits::updateGuess(curve_->data_, *guessIt, i);
258 ++guessIt;
259 ++i;
260 }
261
262 curve_->interpolation_.update();
263
264 Real penalty = 0.0;
265 helper_iterator instIt = rateHelpersStart_;
266 while (instIt != rateHelpersEnd_) {
267 Real quoteError = (*instIt)->quoteError();
268 penalty += std::fabs(x: quoteError);
269 ++instIt;
270 }
271 return penalty;
272 }
273
274 template <class Curve>
275 Array PenaltyFunction<Curve>::values(const Array& x) const {
276 Array::const_iterator guessIt = x.begin();
277 Size i = initialIndex_;
278 while (guessIt != x.end()) {
279 Traits::updateGuess(curve_->data_, *guessIt, i);
280 ++guessIt;
281 ++i;
282 }
283
284 curve_->interpolation_.update();
285
286 Array penalties(localisation_);
287 helper_iterator instIt = rateHelpersStart_;
288 Array::iterator penIt = penalties.begin();
289 while (instIt != rateHelpersEnd_) {
290 Real quoteError = (*instIt)->quoteError();
291 *penIt = std::fabs(x: quoteError);
292 ++instIt;
293 ++penIt;
294 }
295 return penalties;
296 }
297
298}
299
300#endif
301

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