QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
kernelinterpolation2d.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/*! \file kernelinterpolation2d.hpp
21 \brief 2D Kernel interpolation
22*/
23
24#ifndef quantlib_kernel_interpolation2D_hpp
25#define quantlib_kernel_interpolation2D_hpp
26
29#include <utility>
30
31/*
32 Grid Explanation:
33
34 Grid=[ (x1,y1) (x1,y2) (x1,y3)... (x1,yM);
35 (x2,y1) (x2,y2) (x2,y3)... (x2,yM);
36 .
37 .
38 .
39 (xN,y1) (xN,y2) (xN,y3)... (xN,yM);
40 ]
41
42 The Passed variables are:
43 - x which is N dimensional
44 - y which is M dimensional
45 - zData which is NxM dimensional and has the z values
46 corresponding to the grid above.
47 - kernel is a template which needs a Real operator()(Real x) implementation
48*/
49
50
51namespace QuantLib {
52
53 namespace detail {
54
55 template <class I1, class I2, class M, class Kernel>
57 : public Interpolation2D::templateImpl<I1,I2,M> {
58
59 public:
60 KernelInterpolation2DImpl(const I1& xBegin,
61 const I1& xEnd,
62 const I2& yBegin,
63 const I2& yEnd,
64 const M& zData,
65 Kernel kernel)
66 : Interpolation2D::templateImpl<I1, I2, M>(xBegin, xEnd, yBegin, yEnd, zData),
67 xSize_(Size(xEnd - xBegin)), ySize_(Size(yEnd - yBegin)), xySize_(xSize_ * ySize_),
68 alphaVec_(xySize_), yVec_(xySize_), M_(xySize_, xySize_), kernel_(std::move(kernel)) {
69
71 "Z value matrix has wrong number of rows");
73 "Z value matrix has wrong number of columns");
74 }
75
76 void calculate() override { updateAlphaVec(); }
77
78 Real value(Real x1, Real x2) const override {
79
80 Real res=0.0;
81
82 Array X(2),Xn(2);
83 X[0]=x1;X[1]=x2;
84
85 Size cnt=0; // counter
86
87 for( Size j=0; j< ySize_;++j){
88 for( Size i=0; i< xSize_;++i){
89 Xn[0]=this->xBegin_[i];
90 Xn[1]=this->yBegin_[j];
91 res+=alphaVec_[cnt]*kernelAbs(X,Xn);
92 cnt++;
93 }
94 }
95 return res/gammaFunc(X);
96 }
97
98 // the calculation will solve y=M*a for a. Due to
99 // singularity or rounding errors the recalculation
100 // M*a may not give y. Here, a failure will be thrown if
101 // |M*a-y|>=invPrec_
103 invPrec_=invPrec;
104 }
105
106 private:
107
108 // returns K(||X-Y||) where X,Y are vectors
109 Real kernelAbs(const Array& X, const Array& Y)const{
110 return kernel_(Norm2(X-Y));
111 }
112
113 Real gammaFunc(const Array& X)const{
114
115 Real res=0.0;
116 Array Xn(X.size());
117
118 for(Size j=0; j< ySize_;++j){
119 for(Size i=0; i< xSize_;++i){
120 Xn[0]=this->xBegin_[i];
121 Xn[1]=this->yBegin_[j];
122 res+=kernelAbs(X,Xn);
123 }
124 }
125
126 return res;
127 }
128
130 // Function calculates the alpha vector with given
131 // fixed pillars+values
132
133 Array Xk(2),Xn(2);
134
135 Size rowCnt=0,colCnt=0;
136 Real tmpVar=0.0;
137
138 // write y-vector and M-Matrix
139 for(Size j=0; j< ySize_;++j){
140 for(Size i=0; i< xSize_;++i){
141
142 yVec_[rowCnt]=this->zData_[i][j];
143 // calculate X_k
144 Xk[0]=this->xBegin_[i];
145 Xk[1]=this->yBegin_[j];
146
147 tmpVar=1/gammaFunc(Xk);
148 colCnt=0;
149
150 for(Size jM=0; jM< ySize_;++jM){
151 for(Size iM=0; iM< xSize_;++iM){
152 Xn[0]=this->xBegin_[iM];
153 Xn[1]=this->yBegin_[jM];
154 M_[rowCnt][colCnt]=kernelAbs(Xk,Xn)*tmpVar;
155 colCnt++; // increase column counter
156 }// end iM
157 }// end jM
158 rowCnt++; // increase row counter
159 } // end i
160 }// end j
161
163
164 // check if inversion worked up to a reasonable precision.
165 // I've chosen not to check determinant(M_)!=0 before solving
166
167 Array diffVec=Abs(M_*alphaVec_ - yVec_);
168 for (Real i : diffVec) {
169 QL_REQUIRE(i < invPrec_, "inversion failed in 2d kernel interpolation");
170 }
171 }
172
173
175 Real invPrec_ = 1.0e-10;
178 Kernel kernel_;
179 };
180
181 } // end namespace detail
182
183
184 /*! Implementation of the 2D kernel interpolation approach, which
185 can be found in "Foreign Exchange Risk" by Hakala, Wystup page
186 256.
187
188 The kernel in the implementation is kept general, although a
189 Gaussian is considered in the cited text.
190
191 \ingroup interpolations
192 \warning See the Interpolation class for information about the
193 required lifetime of the underlying data.
194 */
196 public:
197 /*! \pre the \f$ x \f$ values must be sorted.
198 \pre kernel needs a Real operator()(Real x) implementation
199 */
200 template <class I1, class I2, class M, class Kernel>
201 KernelInterpolation2D(const I1& xBegin, const I1& xEnd,
202 const I2& yBegin, const I2& yEnd,
203 const M& zData,
204 const Kernel& kernel) {
205
206 impl_ = ext::shared_ptr<Interpolation2D::Impl>(new
208 yBegin, yEnd,
209 zData, kernel));
210 this->update();
211 }
212 };
213}
214
215#endif
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
basic template implementation
const Matrix & zData() const override
templateImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const I2 &yEnd, const M &zData)
base class for 2-D interpolations.
const Matrix & zData() const
ext::shared_ptr< Impl > impl_
KernelInterpolation2D(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const I2 &yEnd, const M &zData, const Kernel &kernel)
Matrix used in linear algebra.
Definition: matrix.hpp:41
Size rows() const
Definition: matrix.hpp:504
Size columns() const
Definition: matrix.hpp:508
Real kernelAbs(const Array &X, const Array &Y) const
Real value(Real x1, Real x2) const override
KernelInterpolation2DImpl(const I1 &xBegin, const I1 &xEnd, const I2 &yBegin, const I2 &yEnd, const M &zData, Kernel kernel)
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
abstract base classes for 2-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
Real Norm2(const Array &v)
Definition: array.hpp:560
STL namespace.
QR decomposition.