1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2007 Allen Kuo
5 Copyright (C) 2010 Alessandro Roveda
6 Copyright (C) 2015 Andres Hernandez
7
8 This file is part of QuantLib, a free-software/open-source library
9 for financial quantitative analysts and developers - http://quantlib.org/
10
11 QuantLib is free software: you can redistribute it and/or modify it
12 under the terms of the QuantLib license. You should have received a
13 copy of the license along with this program; if not, please email
14 <quantlib-dev@lists.sf.net>. The license is also available online at
15 <http://quantlib.org/license.shtml>.
16
17 This program is distributed in the hope that it will be useful, but WITHOUT
18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 FOR A PARTICULAR PURPOSE. See the license for more details.
20*/
21
22#include <ql/math/bernsteinpolynomial.hpp>
23#include <ql/termstructures/yield/nonlinearfittingmethods.hpp>
24#include <utility>
25
26namespace QuantLib {
27
28 ExponentialSplinesFitting::ExponentialSplinesFitting(
29 bool constrainAtZero,
30 const Array& weights,
31 const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
32 const Array& l2,
33 const Real minCutoffTime,
34 const Real maxCutoffTime,
35 const Size numCoeffs,
36 const Real fixedKappa)
37 : FittedBondDiscountCurve::FittingMethod(
38 constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
39 numCoeffs_(numCoeffs), fixedKappa_(fixedKappa)
40 {
41 QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
42 }
43
44 ExponentialSplinesFitting::ExponentialSplinesFitting(bool constrainAtZero,
45 const Array& weights,
46 const Array& l2, const Real minCutoffTime, const Real maxCutoffTime,
47 const Size numCoeffs, const Real fixedKappa)
48 : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), l2,
49 minCutoffTime, maxCutoffTime),
50 numCoeffs_(numCoeffs),fixedKappa_(fixedKappa)
51 {
52 QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
53 }
54
55 ExponentialSplinesFitting::ExponentialSplinesFitting(bool constrainAtZero,
56 const Size numCoeffs,
57 const Real fixedKappa,
58 const Array& weights )
59 : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), Array(),0.0,QL_MAX_REAL),
60 numCoeffs_(numCoeffs), fixedKappa_(fixedKappa)
61 {
62 QL_REQUIRE(ExponentialSplinesFitting::size() > 0, "At least 1 unconstrained coefficient required");
63 }
64
65 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
66 ExponentialSplinesFitting::clone() const {
67 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
68 new ExponentialSplinesFitting(*this));
69 }
70
71 Size ExponentialSplinesFitting::size() const {
72 Size N = constrainAtZero_ ? numCoeffs_ : numCoeffs_ + 1;
73
74 return (fixedKappa_ != Null<Real>()) ? N-1 : N; //One fewer optimization parameters if kappa is fixed
75 }
76
77 DiscountFactor ExponentialSplinesFitting::discountFunction(const Array& x,
78 Time t) const {
79 DiscountFactor d = 0.0;
80 Size N = size();
81 //Use the interal fixedKappa_ member if non-zero, otherwise take kappa from the passed x[] array
82 Real kappa = (fixedKappa_ != Null<Real>()) ? fixedKappa_: x[N-1];
83 Real coeff = 0;
84
85 if (!constrainAtZero_) {
86 for (Size i = 0; i < N - 1; ++i) {
87 d += x[i] * std::exp(x: -kappa * (i + 1) * t);
88 }
89 } else {
90 // notation:
91 // d(t) = coeff* exp(-kappa*1*t) + x[0]* exp(-kappa*2*t) +
92 // x[1]* exp(-kappa*3*t) + ..+ x[7]* exp(-kappa*9*t)
93 for (Size i = 0; i < N - 1; i++) {
94 d += x[i] * std::exp(x: -kappa * (i + 2) * t);
95 coeff += x[i];
96 }
97 coeff = 1.0 - coeff;
98 d += coeff * std::exp(x: -kappa * t);
99 }
100
101 return d;
102 }
103
104
105 NelsonSiegelFitting::NelsonSiegelFitting(
106 const Array& weights,
107 const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
108 const Array& l2,
109 const Real minCutoffTime,
110 const Real maxCutoffTime)
111 : FittedBondDiscountCurve::FittingMethod(
112 true, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime) {}
113
114 NelsonSiegelFitting::NelsonSiegelFitting(const Array& weights,
115 const Array& l2,
116 const Real minCutoffTime, const Real maxCutoffTime)
117 : FittedBondDiscountCurve::FittingMethod(true, weights, ext::shared_ptr<OptimizationMethod>(), l2,
118 minCutoffTime, maxCutoffTime) {}
119
120 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
121 NelsonSiegelFitting::clone() const {
122 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
123 new NelsonSiegelFitting(*this));
124 }
125
126 Size NelsonSiegelFitting::size() const {
127 return 4;
128 }
129
130 DiscountFactor NelsonSiegelFitting::discountFunction(const Array& x,
131 Time t) const {
132 Real kappa = x[size()-1];
133 Real zeroRate = x[0] + (x[1] + x[2])*
134 (1.0 - std::exp(x: -kappa*t))/
135 ((kappa+QL_EPSILON)*(t+QL_EPSILON)) -
136 (x[2])*std::exp(x: -kappa*t);
137 DiscountFactor d = std::exp(x: -zeroRate * t) ;
138 return d;
139 }
140
141
142 SvenssonFitting::SvenssonFitting(const Array& weights,
143 const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
144 const Array& l2,
145 const Real minCutoffTime,
146 const Real maxCutoffTime)
147 : FittedBondDiscountCurve::FittingMethod(
148 true, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime) {}
149
150 SvenssonFitting::SvenssonFitting(const Array& weights,
151 const Array& l2, const Real minCutoffTime, const Real maxCutoffTime)
152 : FittedBondDiscountCurve::FittingMethod(true, weights, ext::shared_ptr<OptimizationMethod>(), l2,
153 minCutoffTime, maxCutoffTime) {}
154
155 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
156 SvenssonFitting::clone() const {
157 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
158 new SvenssonFitting(*this));
159 }
160
161 Size SvenssonFitting::size() const {
162 return 6;
163 }
164
165 DiscountFactor SvenssonFitting::discountFunction(const Array& x,
166 Time t) const {
167 Real kappa = x[size()-2];
168 Real kappa_1 = x[size()-1];
169
170 Real zeroRate = x[0] + (x[1] + x[2])*
171 (1.0 - std::exp(x: -kappa*t))/
172 ((kappa+QL_EPSILON)*(t+QL_EPSILON)) -
173 (x[2])*std::exp(x: -kappa*t) +
174 x[3]* (((1.0 - std::exp(x: -kappa_1*t))/((kappa_1+QL_EPSILON)*(t+QL_EPSILON)))- std::exp(x: -kappa_1*t));
175 DiscountFactor d = std::exp(x: -zeroRate * t) ;
176 return d;
177 }
178
179
180 CubicBSplinesFitting::CubicBSplinesFitting(
181 const std::vector<Time>& knots,
182 bool constrainAtZero,
183 const Array& weights,
184 const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
185 const Array& l2,
186 const Real minCutoffTime,
187 const Real maxCutoffTime)
188 : FittedBondDiscountCurve::FittingMethod(
189 constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
190 splines_(3, knots.size() - 5, knots) {
191
192 QL_REQUIRE(knots.size() >= 8,
193 "At least 8 knots are required" );
194 Size basisFunctions = knots.size() - 4;
195
196 if (constrainAtZero) {
197 size_ = basisFunctions-1;
198
199 // Note: A small but nonzero N_th basis function at t=0 may
200 // lead to an ill conditioned problem
201 N_ = 1;
202
203 QL_REQUIRE(std::abs(splines_(N_, 0.0)) > QL_EPSILON,
204 "N_th cubic B-spline must be nonzero at t=0");
205 } else {
206 size_ = basisFunctions;
207 N_ = 0;
208 }
209 }
210
211 CubicBSplinesFitting::CubicBSplinesFitting(const std::vector<Time>& knots,
212 bool constrainAtZero,
213 const Array& weights,
214 const Array& l2,
215 const Real minCutoffTime, const Real maxCutoffTime)
216 : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights, ext::shared_ptr<OptimizationMethod>(), l2,
217 minCutoffTime, maxCutoffTime),
218 splines_(3, knots.size() - 5, knots) {
219
220 QL_REQUIRE(knots.size() >= 8,
221 "At least 8 knots are required");
222 Size basisFunctions = knots.size() - 4;
223
224 if (constrainAtZero) {
225 size_ = basisFunctions - 1;
226
227 // Note: A small but nonzero N_th basis function at t=0 may
228 // lead to an ill conditioned problem
229 N_ = 1;
230
231 QL_REQUIRE(std::abs(splines_(N_, 0.0)) > QL_EPSILON,
232 "N_th cubic B-spline must be nonzero at t=0");
233 }
234 else {
235 size_ = basisFunctions;
236 N_ = 0;
237 }
238 }
239
240 Real CubicBSplinesFitting::basisFunction(Integer i, Time t) const {
241 return splines_(i,t);
242 }
243
244 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
245 CubicBSplinesFitting::clone() const {
246 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
247 new CubicBSplinesFitting(*this));
248 }
249
250 Size CubicBSplinesFitting::size() const {
251 return size_;
252 }
253
254 DiscountFactor CubicBSplinesFitting::discountFunction(const Array& x,
255 Time t) const {
256 DiscountFactor d = 0.0;
257
258 if (!constrainAtZero_) {
259 for (Size i=0; i<size_; ++i) {
260 d += x[i] * splines_(i,t);
261 }
262 } else {
263 const Real T = 0.0;
264 Real sum = 0.0;
265 for (Size i=0; i<size_; ++i) {
266 if (i < N_) {
267 d += x[i] * splines_(i,t);
268 sum += x[i] * splines_(i,T);
269 } else {
270 d += x[i] * splines_(i+1,t);
271 sum += x[i] * splines_(i+1,T);
272 }
273 }
274 Real coeff = 1.0 - sum;
275 coeff /= splines_(N_,T);
276 d += coeff * splines_(N_,t);
277 }
278
279 return d;
280 }
281
282
283 SimplePolynomialFitting::SimplePolynomialFitting(
284 Natural degree,
285 bool constrainAtZero,
286 const Array& weights,
287 const ext::shared_ptr<OptimizationMethod>& optimizationMethod,
288 const Array& l2,
289 const Real minCutoffTime,
290 const Real maxCutoffTime)
291 : FittedBondDiscountCurve::FittingMethod(
292 constrainAtZero, weights, optimizationMethod, l2, minCutoffTime, maxCutoffTime),
293 size_(constrainAtZero ? degree : degree + 1) {}
294
295 SimplePolynomialFitting::SimplePolynomialFitting(Natural degree, bool constrainAtZero,
296 const Array& weights, const Array& l2,
297 const Real minCutoffTime, const Real maxCutoffTime)
298 : FittedBondDiscountCurve::FittingMethod(constrainAtZero, weights,
299 ext::shared_ptr<OptimizationMethod>(), l2, minCutoffTime, maxCutoffTime),
300 size_(constrainAtZero ? degree : degree + 1) {}
301
302 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
303 SimplePolynomialFitting::clone() const {
304 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
305 new SimplePolynomialFitting(*this));
306 }
307
308 Size SimplePolynomialFitting::size() const {
309 return size_;
310 }
311
312 DiscountFactor SimplePolynomialFitting::discountFunction(const Array& x,
313 Time t) const {
314 DiscountFactor d = 0.0;
315
316 if (!constrainAtZero_) {
317 for (Size i=0; i<size_; ++i)
318 d += x[i] * BernsteinPolynomial::get(i,n: i,x: t);
319 } else {
320 d = 1.0;
321 for (Size i=0; i<size_; ++i)
322 d += x[i] * BernsteinPolynomial::get(i: i+1,n: i+1,x: t);
323 }
324 return d;
325 }
326
327 SpreadFittingMethod::SpreadFittingMethod(const ext::shared_ptr<FittingMethod>& method,
328 Handle<YieldTermStructure> discountCurve,
329 const Real minCutoffTime,
330 const Real maxCutoffTime)
331 : FittedBondDiscountCurve::FittingMethod(
332 method != nullptr ? method->constrainAtZero() : true,
333 method != nullptr ? method->weights() : Array(),
334 method != nullptr ? method->optimizationMethod() : ext::shared_ptr<OptimizationMethod>(),
335 method != nullptr ? method->l2() : Array(),
336 minCutoffTime,
337 maxCutoffTime),
338 method_(method), discountingCurve_(std::move(discountCurve)) {
339 QL_REQUIRE(method, "Fitting method is empty");
340 QL_REQUIRE(!discountingCurve_.empty(), "Discounting curve cannot be empty");
341 }
342
343 std::unique_ptr<FittedBondDiscountCurve::FittingMethod>
344 SpreadFittingMethod::clone() const {
345 return std::unique_ptr<FittedBondDiscountCurve::FittingMethod>(
346 new SpreadFittingMethod(*this));
347 }
348
349 Size SpreadFittingMethod::size() const {
350 return method_->size();
351 }
352
353 DiscountFactor SpreadFittingMethod::discountFunction(const Array& x, Time t) const{
354 return method_->discount(x, t)*discountingCurve_->discount(t, extrapolate: true)/rebase_;
355 }
356
357 void SpreadFittingMethod::init(){
358 //In case discount curve has a different reference date,
359 //discount to this curve's reference date
360 if (curve_->referenceDate() != discountingCurve_->referenceDate()){
361 rebase_ = discountingCurve_->discount(d: curve_->referenceDate());
362 }
363 else{
364 rebase_ = 1.0;
365 }
366 //Call regular init
367 FittedBondDiscountCurve::FittingMethod::init();
368 }
369}
370
371

source code of quantlib/ql/termstructures/yield/nonlinearfittingmethods.cpp