QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
kernelinterpolation.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) 2009 Dimitri Reiswich
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
21#ifndef quantlib_kernel_interpolation_hpp
22#define quantlib_kernel_interpolation_hpp
23
26#include <utility>
27
28/*! \file kernelinterpolation.hpp
29 \brief Kernel interpolation
30*/
31
32namespace QuantLib {
33
34 namespace detail {
35
36 template <class I1, class I2, class Kernel>
38 : public Interpolation::templateImpl<I1,I2> {
39 public:
40 KernelInterpolationImpl(const I1& xBegin,
41 const I1& xEnd,
42 const I2& yBegin,
43 Kernel kernel,
44 const Real epsilon)
45 : Interpolation::templateImpl<I1, I2>(xBegin, xEnd, yBegin),
46 xSize_(Size(xEnd - xBegin)), invPrec_(epsilon), M_(xSize_, xSize_), alphaVec_(xSize_),
47 yVec_(xSize_), kernel_(std::move(kernel)) {}
48
49 void update() override { updateAlphaVec(); }
50
51 Real value(Real x) const override {
52
53 Real res=0.0;
54
55 for( Size i=0; i< xSize_;++i){
56 res+=alphaVec_[i]*kernelAbs(x,this->xBegin_[i]);
57 }
58
59 return res/gammaFunc(x);
60 }
61
62 Real primitive(Real) const override {
63 QL_FAIL("Primitive calculation not implemented "
64 "for kernel interpolation");
65 }
66
67 Real derivative(Real) const override {
68 QL_FAIL("First derivative calculation not implemented "
69 "for kernel interpolation");
70 }
71
72 Real secondDerivative(Real) const override {
73 QL_FAIL("Second derivative calculation not implemented "
74 "for kernel interpolation");
75 }
76
77 private:
78
79 Real kernelAbs(Real x1, Real x2)const{
80 return kernel_(std::fabs(x1-x2));
81 }
82
84
85 Real res=0.0;
86
87 for(Size i=0; i< xSize_;++i){
88 res+=kernelAbs(x,this->xBegin_[i]);
89 }
90 return res;
91 }
92
94 // Function calculates the alpha vector with given
95 // fixed pillars+values
96
97 // Write Matrix M
98 Real tmp=0.0;
99
100 for(Size rowIt=0; rowIt<xSize_;++rowIt){
101
102 yVec_[rowIt]=this->yBegin_[rowIt];
103 tmp=1.0/gammaFunc(this->xBegin_[rowIt]);
104
105 for(Size colIt=0; colIt<xSize_;++colIt){
106 M_[rowIt][colIt]=kernelAbs(this->xBegin_[rowIt],
107 this->xBegin_[colIt])*tmp;
108 }
109 }
110
111 // Solve y=M*\alpha for \alpha
113
114 // check if inversion worked up to a reasonable precision.
115 // I've chosen not to check determinant(M_)!=0 before solving
116
117 Array diffVec=Abs(M_*alphaVec_ - yVec_);
118
119 for (Real i : diffVec) {
120 QL_REQUIRE(i < invPrec_, "Inversion failed in 1d kernel interpolation");
121 }
122 }
123
128 Kernel kernel_;
129 };
130
131 } // end namespace detail
132
133
134 //! Kernel interpolation between discrete points
135 /*! Implementation of the kernel interpolation approach, which can
136 be found in "Foreign Exchange Risk" by Hakala, Wystup page
137 256.
138
139 The kernel in the implementation is kept general, although a Gaussian
140 is considered in the cited text.
141
142 \ingroup interpolations
143 \warning See the Interpolation class for information about the
144 required lifetime of the underlying data.
145 */
147 public:
148 /*! \pre the \f$ x \f$ values must be sorted.
149 \pre kernel needs a Real operator()(Real x) implementation
150
151 The calculation will solve \f$ y = Ma \f$ for \f$a\f$.
152 Due to singularity or rounding errors the recalculation
153 \f$ Ma \f$ may not give \f$ y\f$. Here, a failure will
154 be thrown if
155 \f[
156 \left\| Ma-y \right\|_\infty \geq \epsilon
157 \f] */
158 template <class I1, class I2, class Kernel>
159 KernelInterpolation(const I1& xBegin, const I1& xEnd,
160 const I2& yBegin,
161 const Kernel& kernel,
162 const double epsilon = 1.0E-7) {
163 impl_ = ext::shared_ptr<Interpolation::Impl>(new
165 yBegin, kernel,
166 epsilon));
167 impl_->update();
168 }
169
170 };
171}
172
173#endif
1-D array used in linear algebra.
Definition: array.hpp:52
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_
Kernel interpolation between discrete points.
KernelInterpolation(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const Kernel &kernel, const double epsilon=1.0E-7)
Matrix used in linear algebra.
Definition: matrix.hpp:41
KernelInterpolationImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, Kernel kernel, const Real epsilon)
#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
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
Definition: any.hpp:35
Array qrSolve(const Matrix &a, const Array &b, bool pivot, const Array &d)
QR Solve.
Array Abs(const Array &v)
Definition: array.hpp:833
STL namespace.
QR decomposition.