generallinearleastsquares.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
 Copyright (C) 2009 Dirk Eddelbuettel
 Copyright (C) 2006, 2009, 2010 Klaus Spanderen
 Copyright (C) 2010 Kakhkhor Abdijalilov
 Copyright (C) 2010 Slava Mazur
 
 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/
 
 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <quantlib-dev@lists.sf.net>. The license is also available online at
 <http://quantlib.org/license.shtml>.
 
 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/
 
#ifndef quantlib_general_linear_least_squares_hpp
#define quantlib_general_linear_least_squares_hpp
 
#include <ql/qldefines.hpp>
#include <ql/math/matrixutilities/svd.hpp>
#include <ql/math/array.hpp>
#include <vector>
 
namespace QuantLib {
 
 
    class GeneralLinearLeastSquares {
    public:
        template <class xContainer, class yContainer, class vContainer>
        GeneralLinearLeastSquares(const xContainer & x,
                                  const yContainer & y, const vContainer & v);
 
        template<class xIterator, class yIterator, class vIterator>
        GeneralLinearLeastSquares(xIterator xBegin, xIterator xEnd,
                                  yIterator yBegin, yIterator yEnd,
                                  vIterator vBegin, vIterator vEnd);
 
        const Array& coefficients()   const { return a_; }
        const Array& residuals()      const { return residuals_; }
 
        const Array& standardErrors() const { return standardErrors_; }
        const Array& error()          const { return err_;}
 
        Size size() const { return residuals_.size(); }
 
        Size dim() const { return a_.size(); }
 
    protected:
        Array a_, err_, residuals_, standardErrors_;
 
        template <class xIterator, class yIterator, class vIterator>
        void calculate(
            xIterator xBegin, xIterator xEnd,
            yIterator yBegin, yIterator yEnd,
            vIterator vBegin);
    };
 
    template <class xContainer, class yContainer, class vContainer> inline
    GeneralLinearLeastSquares::GeneralLinearLeastSquares(const xContainer & x,
                                                         const yContainer & y,
                                                         const vContainer & v)
    : a_(v.size(), 0.0),
      err_(v.size(), 0.0),
      residuals_(y.size()),
      standardErrors_(v.size()) {
        calculate(x.begin(), x.end(), y.begin(), y.end(), v.begin());
    }
 
    template<class xIterator, class yIterator, class vIterator> inline
    GeneralLinearLeastSquares::GeneralLinearLeastSquares(
                                            xIterator xBegin, xIterator xEnd,
                                            yIterator yBegin, yIterator yEnd,
                                            vIterator vBegin, vIterator vEnd)
    : a_(std::distance(vBegin, vEnd), 0.0),
      err_(a_.size(), 0.0),
      residuals_(std::distance(yBegin, yEnd)),
      standardErrors_(a_.size()) {
        calculate(xBegin, xEnd, yBegin, yEnd, vBegin);
    }
 
 
    template <class xIterator, class yIterator, class vIterator>
    void GeneralLinearLeastSquares::calculate(xIterator xBegin, xIterator xEnd,
                                              yIterator yBegin, yIterator yEnd,
                                              vIterator vBegin) {
 
        const Size n = residuals_.size();
        const Size m = err_.size();
 
        QL_REQUIRE( n == Size(std::distance(yBegin, yEnd)),
            "sample set need to be of the same size");
        QL_REQUIRE(n >= m, "sample set is too small");
 
        Size i;
 
        Matrix A(n, m);
        for (i=0; i<m; ++i)
            std::transform(xBegin, xEnd, A.column_begin(i), *vBegin++);
 
        const SVD svd(A);
        const Matrix& V = svd.V();
        const Matrix& U = svd.U();
        const Array& w = svd.singularValues();
        const Real threshold = n * QL_EPSILON * svd.singularValues()[0];
 
        for (i=0; i<m; ++i) {
            if (w[i] > threshold) {
                const Real u = std::inner_product(U.column_begin(i),
                    U.column_end(i),
                    yBegin, Real(0.0))/w[i];
 
                for (Size j=0; j<m; ++j) {
                    a_[j]  +=u*V[j][i];
                    err_[j]+=V[j][i]*V[j][i]/(w[i]*w[i]);
                }
            }
        }
        err_      = Sqrt(err_);
        Array tmp = A*a_;
        std::transform(tmp.begin(), tmp.end(), yBegin, residuals_.begin(), std::minus<>());
 
        const Real chiSq
            = std::inner_product(residuals_.begin(), residuals_.end(), residuals_.begin(), Real(0.0));
        const Real multiplier = std::sqrt(chiSq/(n-2));
        std::transform(err_.begin(), err_.end(), standardErrors_.begin(),
                       [=](Real x) -> Real { return x * multiplier; });
    }
 
}
 
#endif