QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
cubicinterpolation.hpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2000, 2001, 2002, 2003 RiskMap srl
5 Copyright (C) 2001, 2002, 2003 Nicolas Di Césaré
6 Copyright (C) 2004, 2008, 2009, 2011 Ferdinando Ametrano
7 Copyright (C) 2009 Sylvain Bertrand
8 Copyright (C) 2013 Peter Caspers
9 Copyright (C) 2016 Nicholas Bertocchi
10
11 This file is part of QuantLib, a free-software/open-source library
12 for financial quantitative analysts and developers - http://quantlib.org/
13
14 QuantLib is free software: you can redistribute it and/or modify it
15 under the terms of the QuantLib license. You should have received a
16 copy of the license along with this program; if not, please email
17 <quantlib-dev@lists.sf.net>. The license is also available online at
18 <http://quantlib.org/license.shtml>.
19
20 This program is distributed in the hope that it will be useful, but WITHOUT
21 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
22 FOR A PARTICULAR PURPOSE. See the license for more details.
23*/
24
25/*! \file cubicinterpolation.hpp
26 \brief cubic interpolation between discrete points
27*/
28
29#ifndef quantlib_cubic_interpolation_hpp
30#define quantlib_cubic_interpolation_hpp
31
32#include <ql/math/matrix.hpp>
35#include <vector>
36
37namespace QuantLib {
38
39 namespace detail {
40
42 public:
44 : n_(n), primitiveConst_(n-1), a_(n-1), b_(n-1), c_(n-1),
46 virtual ~CoefficientHolder() = default;
48 // P[i](x) = y[i] +
49 // a[i]*(x-x[i]) +
50 // b[i]*(x-x[i])^2 +
51 // c[i]*(x-x[i])^3
52 std::vector<Real> primitiveConst_, a_, b_, c_;
53 std::vector<bool> monotonicityAdjustments_;
54 };
55
56 template <class I1, class I2> class CubicInterpolationImpl;
57
58 }
59
60 //! %Cubic interpolation between discrete points.
61 /*! Cubic interpolation is fully defined when the ${f_i}$ function values
62 at points ${x_i}$ are supplemented with ${f^'_i}$ function derivative
63 values.
64
65 Different type of first derivative approximations are implemented,
66 both local and non-local. Local schemes (Fourth-order, Parabolic,
67 Modified Parabolic, Fritsch-Butland, Akima, Kruger) use only $f$ values
68 near $x_i$ to calculate each $f^'_i$. Non-local schemes (Spline with
69 different boundary conditions) use all ${f_i}$ values and obtain
70 ${f^'_i}$ by solving a linear system of equations. Local schemes
71 produce $C^1$ interpolants, while the spline schemes generate $C^2$
72 interpolants.
73
74 Hyman's monotonicity constraint filter is also implemented: it can be
75 applied to all schemes to ensure that in the regions of local
76 monotoniticity of the input (three successive increasing or decreasing
77 values) the interpolating cubic remains monotonic. If the interpolating
78 cubic is already monotonic, the Hyman filter leaves it unchanged
79 preserving all its original features.
80
81 In the case of $C^2$ interpolants the Hyman filter ensures local
82 monotonicity at the expense of the second derivative of the interpolant
83 which will no longer be continuous in the points where the filter has
84 been applied.
85
86 While some non-linear schemes (Modified Parabolic, Fritsch-Butland,
87 Kruger) are guaranteed to be locally monotonic in their original
88 approximation, all other schemes must be filtered according to the
89 Hyman criteria at the expense of their linearity.
90
91 See R. L. Dougherty, A. Edelman, and J. M. Hyman,
92 "Nonnegativity-, Monotonicity-, or Convexity-Preserving CubicSpline and
93 Quintic Hermite Interpolation"
94 Mathematics Of Computation, v. 52, n. 186, April 1989, pp. 471-494.
95
96 \todo implement missing schemes (FourthOrder and ModifiedParabolic) and
97 missing boundary conditions (Periodic and Lagrange).
98
99 \test to be adapted from old ones.
100
101 \ingroup interpolations
102 \warning See the Interpolation class for information about the
103 required lifetime of the underlying data.
104 */
106 public:
108 /*! Spline approximation (non-local, non-monotonic, linear[?]).
109 Different boundary conditions can be used on the left and right
110 boundaries: see BoundaryCondition.
111 */
113
114 //! Overshooting minimization 1st derivative
116
117 //! Overshooting minimization 2nd derivative
119
120 //! Fourth-order approximation (local, non-monotonic, linear)
122
123 //! Parabolic approximation (local, non-monotonic, linear)
125
126 //! Fritsch-Butland approximation (local, monotonic, non-linear)
128
129 //! Akima approximation (local, non-monotonic, non-linear)
131
132 //! Kruger approximation (local, monotonic, non-linear)
134
135 //! Weighted harmonic mean approximation (local, monotonic, non-linear)
137 };
139 //! Make second(-last) point an inactive knot
141
142 //! Match value of end-slope
144
145 //! Match value of second derivative at end
147
148 //! Match first and second derivative at either end
150
151 /*! Match end-slope to the slope of the cubic that matches
152 the first four data at the respective end
153 */
155 };
156 /*! \pre the \f$ x \f$ values must be sorted. */
157 template <class I1, class I2>
158 CubicInterpolation(const I1& xBegin,
159 const I1& xEnd,
160 const I2& yBegin,
162 bool monotonic,
164 Real leftConditionValue,
166 Real rightConditionValue) {
167 impl_ = ext::shared_ptr<Interpolation::Impl>(new
168 detail::CubicInterpolationImpl<I1,I2>(xBegin, xEnd, yBegin,
169 da,
170 monotonic,
171 leftCond,
172 leftConditionValue,
173 rightCond,
174 rightConditionValue));
175 impl_->update();
176 }
177 const std::vector<Real>& primitiveConstants() const {
178 return coeffs().primitiveConst_;
179 }
180 const std::vector<Real>& aCoefficients() const { return coeffs().a_; }
181 const std::vector<Real>& bCoefficients() const { return coeffs().b_; }
182 const std::vector<Real>& cCoefficients() const { return coeffs().c_; }
183 const std::vector<bool>& monotonicityAdjustments() const {
185 }
186 private:
188 return *dynamic_cast<detail::CoefficientHolder*>(impl_.get());
189 }
190 };
191
192
193 // convenience classes
194
196 public:
197 /*! \pre the \f$ x \f$ values must be sorted. */
198 template <class I1, class I2>
199 CubicNaturalSpline(const I1& xBegin,
200 const I1& xEnd,
201 const I2& yBegin)
202 : CubicInterpolation(xBegin, xEnd, yBegin,
203 Spline, false,
204 SecondDerivative, 0.0,
205 SecondDerivative, 0.0) {}
206 };
207
209 public:
210 /*! \pre the \f$ x \f$ values must be sorted. */
211 template <class I1, class I2>
213 const I1& xEnd,
214 const I2& yBegin)
215 : CubicInterpolation(xBegin, xEnd, yBegin,
216 Spline, true,
217 SecondDerivative, 0.0,
218 SecondDerivative, 0.0) {}
219 };
220
222 public:
223 /*! \pre the \f$ x \f$ values must be sorted. */
224 template <class I1, class I2>
226 const I1& xEnd,
227 const I2& yBegin)
228 : CubicInterpolation(xBegin, xEnd, yBegin,
229 SplineOM1, false,
230 SecondDerivative, 0.0,
231 SecondDerivative, 0.0) {}
232 };
233
235 public:
236 /*! \pre the \f$ x \f$ values must be sorted. */
237 template <class I1, class I2>
239 const I1& xEnd,
240 const I2& yBegin)
241 : CubicInterpolation(xBegin, xEnd, yBegin,
242 SplineOM2, false,
243 SecondDerivative, 0.0,
244 SecondDerivative, 0.0) {}
245 };
246
248 public:
249 /*! \pre the \f$ x \f$ values must be sorted. */
250 template <class I1, class I2>
251 AkimaCubicInterpolation(const I1& xBegin,
252 const I1& xEnd,
253 const I2& yBegin)
254 : CubicInterpolation(xBegin, xEnd, yBegin,
255 Akima, false,
256 SecondDerivative, 0.0,
257 SecondDerivative, 0.0) {}
258 };
259
261 public:
262 /*! \pre the \f$ x \f$ values must be sorted. */
263 template <class I1, class I2>
264 KrugerCubic(const I1& xBegin,
265 const I1& xEnd,
266 const I2& yBegin)
267 : CubicInterpolation(xBegin, xEnd, yBegin,
268 Kruger, false,
269 SecondDerivative, 0.0,
270 SecondDerivative, 0.0) {}
271 };
272
274 public:
275 /*! \pre the \f$ x \f$ values must be sorted. */
276 template <class I1, class I2>
277 HarmonicCubic(const I1& xBegin,
278 const I1& xEnd,
279 const I2& yBegin)
280 : CubicInterpolation(xBegin, xEnd, yBegin,
281 Harmonic, false,
282 SecondDerivative, 0.0,
283 SecondDerivative, 0.0) {}
284 };
285
287 public:
288 /*! \pre the \f$ x \f$ values must be sorted. */
289 template <class I1, class I2>
290 FritschButlandCubic(const I1& xBegin,
291 const I1& xEnd,
292 const I2& yBegin)
293 : CubicInterpolation(xBegin, xEnd, yBegin,
294 FritschButland, true,
295 SecondDerivative, 0.0,
296 SecondDerivative, 0.0) {}
297 };
298
300 public:
301 /*! \pre the \f$ x \f$ values must be sorted. */
302 template <class I1, class I2>
303 Parabolic(const I1& xBegin,
304 const I1& xEnd,
305 const I2& yBegin)
306 : CubicInterpolation(xBegin, xEnd, yBegin,
308 SecondDerivative, 0.0,
309 SecondDerivative, 0.0) {}
310 };
311
313 public:
314 /*! \pre the \f$ x \f$ values must be sorted. */
315 template <class I1, class I2>
316 MonotonicParabolic(const I1& xBegin,
317 const I1& xEnd,
318 const I2& yBegin)
319 : CubicInterpolation(xBegin, xEnd, yBegin,
320 Parabolic, true,
321 SecondDerivative, 0.0,
322 SecondDerivative, 0.0) {}
323 };
324
325 //! %Cubic interpolation factory and traits
326 /*! \ingroup interpolations */
327 class Cubic {
328 public:
331 bool monotonic = false,
334 Real leftConditionValue = 0.0,
337 Real rightConditionValue = 0.0)
338 : da_(da), monotonic_(monotonic),
339 leftType_(leftCondition), rightType_(rightCondition),
340 leftValue_(leftConditionValue), rightValue_(rightConditionValue) {}
341 template <class I1, class I2>
342 Interpolation interpolate(const I1& xBegin,
343 const I1& xEnd,
344 const I2& yBegin) const {
345 return CubicInterpolation(xBegin, xEnd, yBegin,
349 }
350 static const bool global = true;
351 static const Size requiredPoints = 2;
352 private:
357 };
358
359
360 namespace detail {
361
362 template <class I1, class I2>
364 public Interpolation::templateImpl<I1,I2> {
365 public:
366 CubicInterpolationImpl(const I1& xBegin,
367 const I1& xEnd,
368 const I2& yBegin,
370 bool monotonic,
372 Real leftConditionValue,
374 Real rightConditionValue)
375 : CoefficientHolder(xEnd-xBegin),
376 Interpolation::templateImpl<I1,I2>(xBegin, xEnd, yBegin,
377 Cubic::requiredPoints),
378 da_(da),
379 monotonic_(monotonic),
380 leftType_(leftCondition), rightType_(rightCondition),
381 leftValue_(leftConditionValue),
382 rightValue_(rightConditionValue),
383 tmp_(n_), dx_(n_-1), S_(n_-1), L_(n_) {
386 QL_REQUIRE((xEnd-xBegin) >= 4,
387 "Lagrange boundary condition requires at least "
388 "4 points (" << (xEnd-xBegin) << " are given)");
389 }
390 }
391
392 void update() override {
393
394 for (Size i=0; i<n_-1; ++i) {
395 dx_[i] = this->xBegin_[i+1] - this->xBegin_[i];
396 S_[i] = (this->yBegin_[i+1] - this->yBegin_[i])/dx_[i];
397 }
398
399 // first derivative approximation
401 for (Size i=1; i<n_-1; ++i) {
402 L_.setMidRow(i, dx_[i], 2.0*(dx_[i]+dx_[i-1]), dx_[i-1]);
403 tmp_[i] = 3.0*(dx_[i]*S_[i-1] + dx_[i-1]*S_[i]);
404 }
405
406 // left boundary condition
407 switch (leftType_) {
409 // ignoring end condition value
410 L_.setFirstRow(dx_[1]*(dx_[1]+dx_[0]),
411 (dx_[0]+dx_[1])*(dx_[0]+dx_[1]));
412 tmp_[0] = S_[0]*dx_[1]*(2.0*dx_[1]+3.0*dx_[0]) +
413 S_[1]*dx_[0]*dx_[0];
414 break;
416 L_.setFirstRow(1.0, 0.0);
417 tmp_[0] = leftValue_;
418 break;
420 L_.setFirstRow(2.0, 1.0);
421 tmp_[0] = 3.0*S_[0] - leftValue_*dx_[0]/2.0;
422 break;
424 QL_FAIL("this end condition is not implemented yet");
426 L_.setFirstRow(1.0, 0.0);
428 this->xBegin_[0],this->xBegin_[1],
429 this->xBegin_[2],this->xBegin_[3],
430 this->yBegin_[0],this->yBegin_[1],
431 this->yBegin_[2],this->yBegin_[3],
432 this->xBegin_[0]);
433 break;
434 default:
435 QL_FAIL("unknown end condition");
436 }
437
438 // right boundary condition
439 switch (rightType_) {
441 // ignoring end condition value
442 L_.setLastRow(-(dx_[n_-2]+dx_[n_-3])*(dx_[n_-2]+dx_[n_-3]),
443 -dx_[n_-3]*(dx_[n_-3]+dx_[n_-2]));
444 tmp_[n_-1] = -S_[n_-3]*dx_[n_-2]*dx_[n_-2] -
445 S_[n_-2]*dx_[n_-3]*(3.0*dx_[n_-2]+2.0*dx_[n_-3]);
446 break;
448 L_.setLastRow(0.0, 1.0);
449 tmp_[n_-1] = rightValue_;
450 break;
452 L_.setLastRow(1.0, 2.0);
453 tmp_[n_-1] = 3.0*S_[n_-2] + rightValue_*dx_[n_-2]/2.0;
454 break;
456 QL_FAIL("this end condition is not implemented yet");
458 L_.setLastRow(0.0,1.0);
460 this->xBegin_[n_-4],this->xBegin_[n_-3],
461 this->xBegin_[n_-2],this->xBegin_[n_-1],
462 this->yBegin_[n_-4],this->yBegin_[n_-3],
463 this->yBegin_[n_-2],this->yBegin_[n_-1],
464 this->xBegin_[n_-1]);
465 break;
466 default:
467 QL_FAIL("unknown end condition");
468 }
469
470 // solve the system
473 Matrix T_(n_-2, n_, 0.0);
474 for (Size i=0; i<n_-2; ++i) {
475 T_[i][i]=dx_[i]/6.0;
476 T_[i][i+1]=(dx_[i+1]+dx_[i])/3.0;
477 T_[i][i+2]=dx_[i+1]/6.0;
478 }
479 Matrix S_(n_-2, n_, 0.0);
480 for (Size i=0; i<n_-2; ++i) {
481 S_[i][i]=1.0/dx_[i];
482 S_[i][i+1]=-(1.0/dx_[i+1]+1.0/dx_[i]);
483 S_[i][i+2]=1.0/dx_[i+1];
484 }
485 Matrix Up_(n_, 2, 0.0);
486 Up_[0][0]=1;
487 Up_[n_-1][1]=1;
488 Matrix Us_(n_, n_-2, 0.0);
489 for (Size i=0; i<n_-2; ++i)
490 Us_[i+1][i]=1;
491 Matrix Z_ = Us_*inverse(T_*Us_);
492 Matrix I_(n_, n_, 0.0);
493 for (Size i=0; i<n_; ++i)
494 I_[i][i]=1;
495 Matrix V_ = (I_-Z_*T_)*Up_;
496 Matrix W_ = Z_*S_;
497 Matrix Q_(n_, n_, 0.0);
498 Q_[0][0]=1.0/(n_-1)*dx_[0]*dx_[0]*dx_[0];
499 Q_[0][1]=7.0/8*1.0/(n_-1)*dx_[0]*dx_[0]*dx_[0];
500 for (Size i=1; i<n_-1; ++i) {
501 Q_[i][i-1]=7.0/8*1.0/(n_-1)*dx_[i-1]*dx_[i-1]*dx_[i-1];
502 Q_[i][i]=1.0/(n_-1)*dx_[i]*dx_[i]*dx_[i]+1.0/(n_-1)*dx_[i-1]*dx_[i-1]*dx_[i-1];
503 Q_[i][i+1]=7.0/8*1.0/(n_-1)*dx_[i]*dx_[i]*dx_[i];
504 }
505 Q_[n_-1][n_-2]=7.0/8*1.0/(n_-1)*dx_[n_-2]*dx_[n_-2]*dx_[n_-2];
506 Q_[n_-1][n_-1]=1.0/(n_-1)*dx_[n_-2]*dx_[n_-2]*dx_[n_-2];
507 Matrix J_ = (I_-V_*inverse(transpose(V_)*Q_*V_)*transpose(V_)*Q_)*W_;
508 Array Y_(n_);
509 for (Size i=0; i<n_; ++i)
510 Y_[i]=this->yBegin_[i];
511 Array D_ = J_*Y_;
512 for (Size i=0; i<n_-1; ++i)
513 tmp_[i]=(Y_[i+1]-Y_[i])/dx_[i]-(2.0*D_[i]+D_[i+1])*dx_[i]/6.0;
514 tmp_[n_-1]=tmp_[n_-2]+D_[n_-2]*dx_[n_-2]+(D_[n_-1]-D_[n_-2])*dx_[n_-2]/2.0;
515
517 Matrix T_(n_-2, n_, 0.0);
518 for (Size i=0; i<n_-2; ++i) {
519 T_[i][i]=dx_[i]/6.0;
520 T_[i][i+1]=(dx_[i]+dx_[i+1])/3.0;
521 T_[i][i+2]=dx_[i+1]/6.0;
522 }
523 Matrix S_(n_-2, n_, 0.0);
524 for (Size i=0; i<n_-2; ++i) {
525 S_[i][i]=1.0/dx_[i];
526 S_[i][i+1]=-(1.0/dx_[i+1]+1.0/dx_[i]);
527 S_[i][i+2]=1.0/dx_[i+1];
528 }
529 Matrix Up_(n_, 2, 0.0);
530 Up_[0][0]=1;
531 Up_[n_-1][1]=1;
532 Matrix Us_(n_, n_-2, 0.0);
533 for (Size i=0; i<n_-2; ++i)
534 Us_[i+1][i]=1;
535 Matrix Z_ = Us_*inverse(T_*Us_);
536 Matrix I_(n_, n_, 0.0);
537 for (Size i=0; i<n_; ++i)
538 I_[i][i]=1;
539 Matrix V_ = (I_-Z_*T_)*Up_;
540 Matrix W_ = Z_*S_;
541 Matrix Q_(n_, n_, 0.0);
542 Q_[0][0]=1.0/(n_-1)*dx_[0];
543 Q_[0][1]=1.0/2*1.0/(n_-1)*dx_[0];
544 for (Size i=1; i<n_-1; ++i) {
545 Q_[i][i-1]=1.0/2*1.0/(n_-1)*dx_[i-1];
546 Q_[i][i]=1.0/(n_-1)*dx_[i]+1.0/(n_-1)*dx_[i-1];
547 Q_[i][i+1]=1.0/2*1.0/(n_-1)*dx_[i];
548 }
549 Q_[n_-1][n_-2]=1.0/2*1.0/(n_-1)*dx_[n_-2];
550 Q_[n_-1][n_-1]=1.0/(n_-1)*dx_[n_-2];
551 Matrix J_ = (I_-V_*inverse(transpose(V_)*Q_*V_)*transpose(V_)*Q_)*W_;
552 Array Y_(n_);
553 for (Size i=0; i<n_; ++i)
554 Y_[i]=this->yBegin_[i];
555 Array D_ = J_*Y_;
556 for (Size i=0; i<n_-1; ++i)
557 tmp_[i]=(Y_[i+1]-Y_[i])/dx_[i]-(2.0*D_[i]+D_[i+1])*dx_[i]/6.0;
558 tmp_[n_-1]=tmp_[n_-2]+D_[n_-2]*dx_[n_-2]+(D_[n_-1]-D_[n_-2])*dx_[n_-2]/2.0;
559 } else { // local schemes
560 if (n_==2)
561 tmp_[0] = tmp_[1] = S_[0];
562 else {
563 switch (da_) {
565 QL_FAIL("FourthOrder not implemented yet");
566 break;
568 // intermediate points
569 for (Size i=1; i<n_-1; ++i)
570 tmp_[i] = (dx_[i-1]*S_[i]+dx_[i]*S_[i-1])/(dx_[i]+dx_[i-1]);
571 // end points
572 tmp_[0] = ((2.0*dx_[ 0]+dx_[ 1])*S_[ 0] - dx_[ 0]*S_[ 1]) / (dx_[ 0]+dx_[ 1]);
573 tmp_[n_-1] = ((2.0*dx_[n_-2]+dx_[n_-3])*S_[n_-2] - dx_[n_-2]*S_[n_-3]) / (dx_[n_-2]+dx_[n_-3]);
574 break;
576 // intermediate points
577 for (Size i=1; i<n_-1; ++i) {
578 Real Smin = std::min(S_[i-1], S_[i]);
579 Real Smax = std::max(S_[i-1], S_[i]);
580 if(Smax+2.0*Smin == 0){
581 if (Smin*Smax < 0)
582 tmp_[i] = QL_MIN_REAL;
583 else if (Smin*Smax == 0)
584 tmp_[i] = 0;
585 else
586 tmp_[i] = QL_MAX_REAL;
587 }
588 else
589 tmp_[i] = 3.0*Smin*Smax/(Smax+2.0*Smin);
590 }
591 // end points
592 tmp_[0] = ((2.0*dx_[ 0]+dx_[ 1])*S_[ 0] - dx_[ 0]*S_[ 1]) / (dx_[ 0]+dx_[ 1]);
593 tmp_[n_-1] = ((2.0*dx_[n_-2]+dx_[n_-3])*S_[n_-2] - dx_[n_-2]*S_[n_-3]) / (dx_[n_-2]+dx_[n_-3]);
594 break;
596 tmp_[0] = (std::abs(S_[1]-S_[0])*2*S_[0]*S_[1]+std::abs(2*S_[0]*S_[1]-4*S_[0]*S_[0]*S_[1])*S_[0])/(std::abs(S_[1]-S_[0])+std::abs(2*S_[0]*S_[1]-4*S_[0]*S_[0]*S_[1]));
597 tmp_[1] = (std::abs(S_[2]-S_[1])*S_[0]+std::abs(S_[0]-2*S_[0]*S_[1])*S_[1])/(std::abs(S_[2]-S_[1])+std::abs(S_[0]-2*S_[0]*S_[1]));
598 for (Size i=2; i<n_-2; ++i) {
599 if ((S_[i-2]==S_[i-1]) && (S_[i]!=S_[i+1]))
600 tmp_[i] = S_[i-1];
601 else if ((S_[i-2]!=S_[i-1]) && (S_[i]==S_[i+1]))
602 tmp_[i] = S_[i];
603 else if (S_[i]==S_[i-1])
604 tmp_[i] = S_[i];
605 else if ((S_[i-2]==S_[i-1]) && (S_[i-1]!=S_[i]) && (S_[i]==S_[i+1]))
606 tmp_[i] = (S_[i-1]+S_[i])/2.0;
607 else
608 tmp_[i] = (std::abs(S_[i+1]-S_[i])*S_[i-1]+std::abs(S_[i-1]-S_[i-2])*S_[i])/(std::abs(S_[i+1]-S_[i])+std::abs(S_[i-1]-S_[i-2]));
609 }
610 tmp_[n_-2] = (std::abs(2*S_[n_-2]*S_[n_-3]-S_[n_-2])*S_[n_-3]+std::abs(S_[n_-3]-S_[n_-4])*S_[n_-2])/(std::abs(2*S_[n_-2]*S_[n_-3]-S_[n_-2])+std::abs(S_[n_-3]-S_[n_-4]));
611 tmp_[n_-1] = (std::abs(4*S_[n_-2]*S_[n_-2]*S_[n_-3]-2*S_[n_-2]*S_[n_-3])*S_[n_-2]+std::abs(S_[n_-2]-S_[n_-3])*2*S_[n_-2]*S_[n_-3])/(std::abs(4*S_[n_-2]*S_[n_-2]*S_[n_-3]-2*S_[n_-2]*S_[n_-3])+std::abs(S_[n_-2]-S_[n_-3]));
612 break;
614 // intermediate points
615 for (Size i=1; i<n_-1; ++i) {
616 if (S_[i-1]*S_[i]<0.0)
617 // slope changes sign at point
618 tmp_[i] = 0.0;
619 else
620 // slope will be between the slopes of the adjacent
621 // straight lines and should approach zero if the
622 // slope of either line approaches zero
623 tmp_[i] = 2.0/(1.0/S_[i-1]+1.0/S_[i]);
624 }
625 // end points
626 tmp_[0] = (3.0*S_[0]-tmp_[1])/2.0;
627 tmp_[n_-1] = (3.0*S_[n_-2]-tmp_[n_-2])/2.0;
628 break;
630 // intermediate points
631 for (Size i=1; i<n_-1; ++i) {
632 Real w1 = 2*dx_[i]+dx_[i-1];
633 Real w2 = dx_[i]+2*dx_[i-1];
634 if (S_[i-1]*S_[i]<=0.0)
635 // slope changes sign at point
636 tmp_[i] = 0.0;
637 else
638 // weighted harmonic mean of S_[i] and S_[i-1] if they
639 // have the same sign; otherwise 0
640 tmp_[i] = (w1+w2)/(w1/S_[i-1]+w2/S_[i]);
641 }
642 // end points [0]
643 tmp_[0] = ((2 * dx_[0] + dx_[1])*S_[0] - dx_[0] * S_[1]) / (dx_[1] + dx_[0]);
644 if (tmp_[0]*S_[0]<0.0) {
645 tmp_[0] = 0;
646 }
647 else if (S_[0]*S_[1]<0) {
648 if (std::fabs(tmp_[0])>std::fabs(3*S_[0])) {
649 tmp_[0] = 3*S_[0];
650 }
651 }
652 // end points [n-1]
653 tmp_[n_-1] = ((2*dx_[n_-2]+dx_[n_-3])*S_[n_-2]-dx_[n_-2]*S_[n_-3])/(dx_[n_-3]+dx_[n_-2]);
654 if (tmp_[n_-1]*S_[n_-2]<0.0) {
655 tmp_[n_-1] = 0;
656 }
657 else if (S_[n_-2]*S_[n_-3]<0) {
658 if (std::fabs(tmp_[n_-1])>std::fabs(3*S_[n_-2])) {
659 tmp_[n_-1] = 3*S_[n_-2];
660 }
661 }
662 break;
663 default:
664 QL_FAIL("unknown scheme");
665 }
666 }
667 }
668
669 std::fill(monotonicityAdjustments_.begin(),
670 monotonicityAdjustments_.end(), false);
671 // Hyman monotonicity constrained filter
672 if (monotonic_) {
673 Real correction;
674 Real pm, pu, pd, M;
675 for (Size i=0; i<n_; ++i) {
676 if (i==0) {
677 if (tmp_[i]*S_[0]>0.0) {
678 correction = tmp_[i]/std::fabs(tmp_[i]) *
679 std::min<Real>(std::fabs(tmp_[i]),
680 std::fabs(3.0*S_[0]));
681 } else {
682 correction = 0.0;
683 }
684 if (correction!=tmp_[i]) {
685 tmp_[i] = correction;
686 monotonicityAdjustments_[i] = true;
687 }
688 } else if (i==n_-1) {
689 if (tmp_[i]*S_[n_-2]>0.0) {
690 correction = tmp_[i]/std::fabs(tmp_[i]) *
691 std::min<Real>(std::fabs(tmp_[i]),
692 std::fabs(3.0*S_[n_-2]));
693 } else {
694 correction = 0.0;
695 }
696 if (correction!=tmp_[i]) {
697 tmp_[i] = correction;
698 monotonicityAdjustments_[i] = true;
699 }
700 } else {
701 pm=(S_[i-1]*dx_[i]+S_[i]*dx_[i-1])/
702 (dx_[i-1]+dx_[i]);
703 M = 3.0 * std::min(std::min(std::fabs(S_[i-1]),
704 std::fabs(S_[i])),
705 std::fabs(pm));
706 if (i>1) {
707 if ((S_[i-1]-S_[i-2])*(S_[i]-S_[i-1])>0.0) {
708 pd=(S_[i-1]*(2.0*dx_[i-1]+dx_[i-2])
709 -S_[i-2]*dx_[i-1])/
710 (dx_[i-2]+dx_[i-1]);
711 if (pm*pd>0.0 && pm*(S_[i-1]-S_[i-2])>0.0) {
712 M = std::max<Real>(M, 1.5*std::min(
713 std::fabs(pm),std::fabs(pd)));
714 }
715 }
716 }
717 if (i<n_-2) {
718 if ((S_[i]-S_[i-1])*(S_[i+1]-S_[i])>0.0) {
719 pu=(S_[i]*(2.0*dx_[i]+dx_[i+1])-S_[i+1]*dx_[i])/
720 (dx_[i]+dx_[i+1]);
721 if (pm*pu>0.0 && -pm*(S_[i]-S_[i-1])>0.0) {
722 M = std::max<Real>(M, 1.5*std::min(
723 std::fabs(pm),std::fabs(pu)));
724 }
725 }
726 }
727 if (tmp_[i]*pm>0.0) {
728 correction = tmp_[i]/std::fabs(tmp_[i]) *
729 std::min(std::fabs(tmp_[i]), M);
730 } else {
731 correction = 0.0;
732 }
733 if (correction!=tmp_[i]) {
734 tmp_[i] = correction;
735 monotonicityAdjustments_[i] = true;
736 }
737 }
738 }
739 }
740
741
742 // cubic coefficients
743 for (Size i=0; i<n_-1; ++i) {
744 a_[i] = tmp_[i];
745 b_[i] = (3.0*S_[i] - tmp_[i+1] - 2.0*tmp_[i])/dx_[i];
746 c_[i] = (tmp_[i+1] + tmp_[i] - 2.0*S_[i])/(dx_[i]*dx_[i]);
747 }
748
749 primitiveConst_[0] = 0.0;
750 for (Size i=1; i<n_-1; ++i) {
752 + dx_[i-1] *
753 (this->yBegin_[i-1] + dx_[i-1] *
754 (a_[i-1]/2.0 + dx_[i-1] *
755 (b_[i-1]/3.0 + dx_[i-1] * c_[i-1]/4.0)));
756 }
757 }
758 Real value(Real x) const override {
759 Size j = this->locate(x);
760 Real dx_ = x-this->xBegin_[j];
761 return this->yBegin_[j] + dx_*(a_[j] + dx_*(b_[j] + dx_*c_[j]));
762 }
763 Real primitive(Real x) const override {
764 Size j = this->locate(x);
765 Real dx_ = x-this->xBegin_[j];
766 return primitiveConst_[j]
767 + dx_*(this->yBegin_[j] + dx_*(a_[j]/2.0
768 + dx_*(b_[j]/3.0 + dx_*c_[j]/4.0)));
769 }
770 Real derivative(Real x) const override {
771 Size j = this->locate(x);
772 Real dx_ = x-this->xBegin_[j];
773 return a_[j] + (2.0*b_[j] + 3.0*c_[j]*dx_)*dx_;
774 }
775 Real secondDerivative(Real x) const override {
776 Size j = this->locate(x);
777 Real dx_ = x-this->xBegin_[j];
778 return 2.0*b_[j] + 6.0*c_[j]*dx_;
779 }
780
781 private:
786 mutable Array tmp_;
787 mutable std::vector<Real> dx_, S_;
789
791 Real a, Real b, Real c, Real d,
792 Real u, Real v, Real w, Real z, Real x) const {
793 return (-((((a-c)*(b-c)*(c-x)*z-(a-d)*(b-d)*(d-x)*w)*(a-x+b-x)
794 +((a-c)*(b-c)*z-(a-d)*(b-d)*w)*(a-x)*(b-x))*(a-b)+
795 ((a-c)*(a-d)*v-(b-c)*(b-d)*u)*(c-d)*(c-x)*(d-x)
796 +((a-c)*(a-d)*(a-x)*v-(b-c)*(b-d)*(b-x)*u)
797 *(c-x+d-x)*(c-d)))/
798 ((a-b)*(a-c)*(a-d)*(b-c)*(b-d)*(c-d));
799 }
800 };
801
802 }
803
804}
805
806#endif
AkimaCubicInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
1-D array used in linear algebra.
Definition: array.hpp:52
Cubic interpolation factory and traits
CubicInterpolation::DerivativeApprox da_
static const bool global
CubicInterpolation::BoundaryCondition rightType_
CubicInterpolation::BoundaryCondition leftType_
static const Size requiredPoints
Cubic(CubicInterpolation::DerivativeApprox da=CubicInterpolation::Kruger, bool monotonic=false, CubicInterpolation::BoundaryCondition leftCondition=CubicInterpolation::SecondDerivative, Real leftConditionValue=0.0, CubicInterpolation::BoundaryCondition rightCondition=CubicInterpolation::SecondDerivative, Real rightConditionValue=0.0)
Interpolation interpolate(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin) const
Cubic interpolation between discrete points.
const std::vector< Real > & cCoefficients() const
const std::vector< bool > & monotonicityAdjustments() const
const std::vector< Real > & bCoefficients() const
CubicInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, CubicInterpolation::DerivativeApprox da, bool monotonic, CubicInterpolation::BoundaryCondition leftCond, Real leftConditionValue, CubicInterpolation::BoundaryCondition rightCond, Real rightConditionValue)
const std::vector< Real > & aCoefficients() const
const std::vector< Real > & primitiveConstants() const
@ Kruger
Kruger approximation (local, monotonic, non-linear)
@ SplineOM2
Overshooting minimization 2nd derivative.
@ Parabolic
Parabolic approximation (local, non-monotonic, linear)
@ Harmonic
Weighted harmonic mean approximation (local, monotonic, non-linear)
@ FourthOrder
Fourth-order approximation (local, non-monotonic, linear)
@ SplineOM1
Overshooting minimization 1st derivative.
@ Akima
Akima approximation (local, non-monotonic, non-linear)
@ FritschButland
Fritsch-Butland approximation (local, monotonic, non-linear)
@ SecondDerivative
Match value of second derivative at end.
@ Periodic
Match first and second derivative at either end.
@ NotAKnot
Make second(-last) point an inactive knot.
@ FirstDerivative
Match value of end-slope.
const detail::CoefficientHolder & coeffs() const
CubicNaturalSpline(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
CubicSplineOvershootingMinimization1(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
CubicSplineOvershootingMinimization2(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
FritschButlandCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
HarmonicCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
basic template implementation
templateImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const int requiredPoints=2)
base class for 1-D interpolations.
ext::shared_ptr< Impl > impl_
KrugerCubic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Matrix used in linear algebra.
Definition: matrix.hpp:41
MonotonicCubicNaturalSpline(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
MonotonicParabolic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Parabolic(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin)
Base implementation for tridiagonal operator.
Array solveFor(const Array &rhs) const
solve linear system for a given right-hand side
void setMidRow(Size, Real, Real, Real)
CubicInterpolation::DerivativeApprox da_
Real cubicInterpolatingPolynomialDerivative(Real a, Real b, Real c, Real d, Real u, Real v, Real w, Real z, Real x) const
CubicInterpolation::BoundaryCondition rightType_
CubicInterpolation::BoundaryCondition leftType_
CubicInterpolationImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, CubicInterpolation::DerivativeApprox da, bool monotonic, CubicInterpolation::BoundaryCondition leftCondition, Real leftConditionValue, CubicInterpolation::BoundaryCondition rightCondition, Real rightConditionValue)
Real secondDerivative(Real x) const override
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_FAIL(message)
throw an error (possibly with file and line information)
Definition: errors.hpp:92
Date d
ext::function< Real(Real)> b
#define QL_MAX_REAL
Definition: qldefines.hpp:176
#define QL_MIN_REAL
Definition: qldefines.hpp:175
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
base class for 1-D interpolations
matrix used in linear algebra.
Definition: any.hpp:35
Matrix transpose(const Matrix &m)
Definition: matrix.hpp:700
Matrix inverse(const Matrix &m)
Definition: matrix.cpp:44
ext::shared_ptr< BlackVolTermStructure > v
tridiagonal operator