QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
generallinearleastsquares.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 Dirk Eddelbuettel
5 Copyright (C) 2006, 2009, 2010 Klaus Spanderen
6 Copyright (C) 2010 Kakhkhor Abdijalilov
7 Copyright (C) 2010 Slava Mazur
8
9 This file is part of QuantLib, a free-software/open-source library
10 for financial quantitative analysts and developers - http://quantlib.org/
11
12 QuantLib is free software: you can redistribute it and/or modify it
13 under the terms of the QuantLib license. You should have received a
14 copy of the license along with this program; if not, please email
15 <quantlib-dev@lists.sf.net>. The license is also available online at
16 <http://quantlib.org/license.shtml>.
17
18 This program is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
20 FOR A PARTICULAR PURPOSE. See the license for more details.
21*/
22
23/*! \file generallinearleastsquares.hpp
24 \brief general linear least square regression
25*/
26
27#ifndef quantlib_general_linear_least_squares_hpp
28#define quantlib_general_linear_least_squares_hpp
29
30#include <ql/qldefines.hpp>
32#include <ql/math/array.hpp>
33#include <vector>
34
35namespace QuantLib {
36
37 //! general linear least squares regression
38 /*! References:
39 "Numerical Recipes in C", 2nd edition,
40 Press, Teukolsky, Vetterling, Flannery,
41
42 \test the correctness of the returned values is tested by
43 checking their properties.
44 */
46 public:
47 template <class xContainer, class yContainer, class vContainer>
48 GeneralLinearLeastSquares(const xContainer & x,
49 const yContainer & y, const vContainer & v);
50
51 template<class xIterator, class yIterator, class vIterator>
52 GeneralLinearLeastSquares(xIterator xBegin, xIterator xEnd,
53 yIterator yBegin, yIterator yEnd,
54 vIterator vBegin, vIterator vEnd);
55
56 const Array& coefficients() const { return a_; }
57 const Array& residuals() const { return residuals_; }
58
59 //! standard parameter errors as given by Excel, R etc.
60 const Array& standardErrors() const { return standardErrors_; }
61 //! modeling uncertainty as definied in Numerical Recipes
62 const Array& error() const { return err_;}
63
64 Size size() const { return residuals_.size(); }
65
66 Size dim() const { return a_.size(); }
67
68 protected:
70
71 template <class xIterator, class yIterator, class vIterator>
72 void calculate(
73 xIterator xBegin, xIterator xEnd,
74 yIterator yBegin, yIterator yEnd,
75 vIterator vBegin);
76 };
77
78 template <class xContainer, class yContainer, class vContainer> inline
80 const yContainer & y,
81 const vContainer & v)
82 : a_(v.size(), 0.0),
83 err_(v.size(), 0.0),
84 residuals_(y.size()),
85 standardErrors_(v.size()) {
86 calculate(x.begin(), x.end(), y.begin(), y.end(), v.begin());
87 }
88
89 template<class xIterator, class yIterator, class vIterator> inline
91 xIterator xBegin, xIterator xEnd,
92 yIterator yBegin, yIterator yEnd,
93 vIterator vBegin, vIterator vEnd)
94 : a_(std::distance(vBegin, vEnd), 0.0),
95 err_(a_.size(), 0.0),
96 residuals_(std::distance(yBegin, yEnd)),
97 standardErrors_(a_.size()) {
98 calculate(xBegin, xEnd, yBegin, yEnd, vBegin);
99 }
100
101
102 template <class xIterator, class yIterator, class vIterator>
103 void GeneralLinearLeastSquares::calculate(xIterator xBegin, xIterator xEnd,
104 yIterator yBegin, yIterator yEnd,
105 vIterator vBegin) {
106
107 const Size n = residuals_.size();
108 const Size m = err_.size();
109
110 QL_REQUIRE( n == Size(std::distance(yBegin, yEnd)),
111 "sample set need to be of the same size");
112 QL_REQUIRE(n >= m, "sample set is too small");
113
114 Size i;
115
116 Matrix A(n, m);
117 for (i=0; i<m; ++i)
118 std::transform(xBegin, xEnd, A.column_begin(i), *vBegin++);
119
120 const SVD svd(A);
121 const Matrix& V = svd.V();
122 const Matrix& U = svd.U();
123 const Array& w = svd.singularValues();
124 const Real threshold = n * QL_EPSILON * svd.singularValues()[0];
125
126 for (i=0; i<m; ++i) {
127 if (w[i] > threshold) {
128 const Real u = std::inner_product(U.column_begin(i),
129 U.column_end(i),
130 yBegin, Real(0.0))/w[i];
131
132 for (Size j=0; j<m; ++j) {
133 a_[j] +=u*V[j][i];
134 err_[j]+=V[j][i]*V[j][i]/(w[i]*w[i]);
135 }
136 }
137 }
138 err_ = Sqrt(err_);
139 Array tmp = A*a_;
140 std::transform(tmp.begin(), tmp.end(), yBegin, residuals_.begin(), std::minus<>());
141
142 const Real chiSq
143 = std::inner_product(residuals_.begin(), residuals_.end(), residuals_.begin(), Real(0.0));
144 const Real multiplier = std::sqrt(chiSq/(n-2));
145 std::transform(err_.begin(), err_.end(), standardErrors_.begin(),
146 [=](Real x) -> Real { return x * multiplier; });
147 }
148
149}
150
151#endif
1-D array used in linear algebra.
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
general linear least squares regression
const Array & error() const
modeling uncertainty as definied in Numerical Recipes
const Array & standardErrors() const
standard parameter errors as given by Excel, R etc.
void calculate(xIterator xBegin, xIterator xEnd, yIterator yBegin, yIterator yEnd, vIterator vBegin)
GeneralLinearLeastSquares(const xContainer &x, const yContainer &y, const vContainer &v)
Matrix used in linear algebra.
Definition: matrix.hpp:41
const_column_iterator column_begin(Size i) const
Definition: matrix.hpp:415
const_column_iterator column_end(Size i) const
Definition: matrix.hpp:434
Singular value decomposition.
Definition: svd.hpp:54
const Matrix & V() const
Definition: svd.cpp:489
const Array & singularValues() const
Definition: svd.cpp:493
const Matrix & U() const
Definition: svd.cpp:485
Real a_
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_EPSILON
Definition: qldefines.hpp:178
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Definition: any.hpp:35
Array Sqrt(const Array &v)
Definition: array.hpp:847
STL namespace.
ext::shared_ptr< BlackVolTermStructure > v
Global definitions and compiler switches.
singular value decomposition