Basis functions for Longstaff-Schwartz method

[Kakhkhor Abdijalilov]

Dear All,
Class LsmBasisSystem (lsmbasissystem.hpp) doesn't look good. Its
implementation uses very inefficient and probability broken method to
generate multidimensional basis functions from single factor basis
functions. It generates many duplicates and then uses "remove-o-zap"
method to remove them. The size of basis is not the same for different
polynomial types.  You can see the test results below. I reworked the
implementation, so that basis sizes are the same for all polynomial
types and no duplicates are created during the process. I would like
to submit my code to QuantLib, but not sure how to do it.

Regards,
Kakhkhor Abdijalilov.


------------------------------------------------------------------------------------------------------------------------
The number of functions returned by multiPathBasisSystem method.

In (n-m), n - size after "remove-o-zap",
m - number of functions removed by "remove-o-zap"

dimension = 1
order; Mnm; Lgr; Hrm; Hpr; Lgn; Chb; Chb2
=========================================
0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0
1; 2-0; 1-1; 2-0; 2-0; 2-0; 2-0; 2-0
2; 3-1; 3-1; 4-0; 3-1; 4-0; 4-0; 4-0
3; 4-4; 5-3; 6-2; 4-4; 7-1; 6-2; 7-1
4; 5-11; 10-6; 10-6; 7-9; 12-4; 10-6; 12-4

dimension = 2
order; Mnm; Lgr; Hrm; Hpr; Lgn; Chb; Chb2
=========================================
0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0
1; 3-0; 1-2; 3-0; 3-0; 3-0; 3-0; 3-0
2; 6-3; 6-3; 8-1; 6-3; 8-1; 8-1; 8-1
3; 10-17; 12-15; 16-11; 10-17; 18-9; 17-10; 18-9
4; 15-66; 32-49; 32-49; 20-61; 38-43; 35-46; 38-43

dimension = 3
order; Mnm; Lgr; Hrm; Hpr; Lgn; Chb; Chb2
=========================================
0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0
1; 4-0; 1-3; 4-0; 4-0; 4-0; 4-0; 4-0
2; 10-6; 10-6; 13-3; 10-6; 13-3; 13-3; 13-3
3; 20-44; 22-42; 32-32; 20-44; 35-29; 34-30; 35-29
4; 35-221; 73-183; 74-182; 44-212; 86-170; 82-174; 86-170

dimension = 4
order; Mnm; Lgr; Hrm; Hpr; Lgn; Chb; Chb2
=========================================
0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0; 1-0
1; 5-0; 1-4; 5-0; 5-0; 5-0; 5-0; 5-0
2; 15-10; 15-10; 19-6; 15-10; 19-6; 19-6; 19-6
3; 35-90; 35-90; 55-70; 35-90; 59-66; 57-68; 59-66
4; 70-555; 140-485; 144-481; 84-541; 164-461; 154-471; 164-461
------------------------------------------------------------------------------------------------------------------------

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Klaus Spanderen-2]

Hi

> ..The size of basis is not the same for different
> polynomial types.  

I've fixed one issue in the SVN.

> .. I reworked the
> implementation, so that basis sizes are the same for all polynomial
> types and no duplicates are created during the process. I would like
> to submit my code to QuantLib, but not sure how to do it.

Please submit the code to this mailing list or via email. I'm going to merge
it into the SVN.

thanks in advance
 Klaus

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Kakhkhor Abdijalilov]

It seems that increasing the threshold in std::find_if algorithm
doesn't solve the problem. I tried it with very high threshold and
there are still many more basis functions than needed. The basis set
size is the same as it is reported in my original message. I guess
"remove-o-zap" works correctly, but something else goes wrong. The
correct number of basis functions should be (Dim+Order-1)! / [
(Dim-1)! * Order!] and it shouldn't depend on polynomial type.

To fix the problem I implemented different basis generation algorithm
which is described below.

The idea is to build multi-factor basis set by multiplying together
single factor basis functions.
Let's say we have the single factor basis:

F_0(.), F_1(.), F_2(.),...

and we want to build a 3-dimensional basis set. In order to build the
order N basis set, we need all terms of the type
F_i(x)*F_j*(y)*F_k(z)  such that i+j+k <= N. We first compute all
tuples (i, j, k) such that i+j+k <= N and then generate the basis set
from those tuples. This would avoid duplicates if the tuple set
contains only unique elements. The required tuple set can be generated
by the following induction algorithm:

Assume that we have all order N (i.e. i+j+k = N) tuples. In order to
build order N+1 tuples do the following:
1) Start with the set:  (N+1, 0, 0), (0, N+1, 0), (0, 0, N+1).
2) For each tuple (i, j, k) from the order N set add the tuples (i+1,
j, k), (i, j+1, k+2) to the order N+1 set.
3) Remove the duplicates from the order N+1 set.

The resulting order N+1 set will contain all order N+1 tuples. The
function next_order_tuples implements the induction algorithm.
Duplicates can be avoided if  we use STL container std::set to store
the tuples.

The code is below. Please feel free to ask any question you might have.

==================================================

// lsmbasissystem.hpp

#ifndef quantlib_lsm_basis_system_hpp
#define quantlib_lsm_basis_system_hpp

#include <ql/qldefines.hpp>
#include <ql/math/array.hpp>
#include <boost/function.hpp>
#include <vector>

namespace QuantLib {

    class LsmBasisSystem {
    public:
        enum PolynomType { Monomial, Laguerre, Hermite, Hyperbolic,
                           Legendre, Chebyshev, Chebyshev2nd };

        static std::vector<boost::function1<Real, Real> >
        pathBasisSystem(Size order, PolynomType polyType);

        static std::vector<boost::function1<Real, Array> >
        multiPathBasisSystem(Size dim, Size order, PolynomType polyType);
    };


} // namespace QuantLib

#endif

==================================================

// lsmbasissystem.hpp

#include <ql/methods/montecarlo/lsmbasissystem.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/integrals/gaussianquadratures.hpp>
#include <boost/bind.hpp>
#include <set>
#include <numeric>

namespace QuantLib {
namespace {

    // makes typing a little easier
    typedef std::vector<boost::function1<Real, Real> > VF_R;
    typedef std::vector<boost::function1<Real, Array> > VF_A;
    typedef std::vector<std::vector<Size> > VV;
    Real (GaussianOrthogonalPolynomial::*ptr_w)(Size, Real) const =
        &GaussianOrthogonalPolynomial::weightedValue;

    // pow(x, order)
    class MonomialFct : public std::unary_function<Real, Real> {
    public:
        MonomialFct(Size order): order_(order) {}
        inline Real operator()(const Real x) const {
            Real ret = 1.0;
            for(Size i=0; i<order_; ++i)
                ret *= x;
            return ret;
        }
    private:
        const Size order_;
    };

    /* multiplies [Real -> Real] functors
       to create [Array -> Real] functor */
    class MultiDimFct : public std::unary_function<Real, Array> {
    public:
        MultiDimFct(const VF_R b): b_(b) {
            QL_REQUIRE(b_.size()>0, "zero size basis");
        }
        inline Real operator()(const Array& a) const {
            #if defined(QL_EXTRA_SAFETY_CHECKS)
            QL_REQUIRE(b_.size()==a.size(), "wrong argument size");
            #endif
            Real ret = b_[0].operator()(a[0]);
            for(Size i=1; i<b_.size(); ++i)
                ret *= b_[i].operator()(a[i]);
            return ret;
        }
    private:
        const VF_R b_;
    };

    // constant functor [Real -> Real] (returns 1.0)
    class MultiDimConstFct : public std::unary_function<Real, Array> {
    public:
        MultiDimConstFct(Size dim): dim_(dim) {
            QL_REQUIRE(dim>0, "zero dimension");
        }
        inline Real operator()(const Array& a) const {
            #if defined(QL_EXTRA_SAFETY_CHECKS)
            QL_REQUIRE(dim_==a.size(), "wrong argument size");
            #endif
            return 1.0;
        }
    private:
        const Size dim_;
    };

    // check size and order of tuples
    void check_tuples(const VV& v, Size dim, Size order) {
        for(Size i=0; i<v.size(); ++i) {
            QL_REQUIRE(dim==v[i].size(), "wrong tuple size");
            QL_REQUIRE(order==std::accumulate(v[i].begin(), v[i].end(), 0),
                "wrong tuple order");
        }
    }

    // build order N+1 tuples from order N tuples
    VV next_order_tuples(const VV& v) {
        const Size order = std::accumulate(v[0].begin(), v[0].end(), 0);
        const Size dim = v[0].size();

        check_tuples(v, dim, order);

        // the set of unique tuples
        std::set<std::vector<Size> > tuples;

        for(Size i=0; i<dim; ++i) {
            // order+1 for i-th dimension, the rest is 0
            std::vector<Size> x(dim);
            x[i] += order+1;
            tuples.insert(x);
            // increase i-th value in every tuple by 1
            for(Size j=0; j<v.size(); ++j) {
                x = v[j];
                x[i] += 1;
                tuples.insert(x);
            }
        }

        VV ret(tuples.begin(), tuples.end());
        return ret;
    }

} // unnamed namespace

    // LsmBasisSystem static methods

    VF_R LsmBasisSystem::pathBasisSystem(Size order, PolynomType polyType) {
        VF_R ret(order+1);
        for (Size i=0; i<=order; ++i) {
            switch (polyType) {
              case Monomial:
                ret[i] = MonomialFct(i);
                break;
              case Laguerre:
                ret[i] = boost::bind(ptr_w, GaussLaguerrePolynomial(), i, _1);
                break;
              case Hermite:
                ret[i] = boost::bind(ptr_w, GaussHermitePolynomial(), i, _1);
                break;
              case Hyperbolic:
                ret[i] = boost::bind(ptr_w, GaussHyperbolicPolynomial(), i, _1);
                break;
              case Legendre:
                ret[i] = boost::bind(ptr_w, GaussLegendrePolynomial(), i, _1);
                break;
              case Chebyshev:
                ret[i] = boost::bind(ptr_w, GaussChebyshevPolynomial(), i, _1);
                break;
              case Chebyshev2nd:
                ret[i] = boost::bind(ptr_w,
GaussChebyshev2ndPolynomial(), i, _1);
                break;
              default:
                QL_FAIL("unknown regression type");
            }
        }
        return ret;
    }

    VF_A LsmBasisSystem::multiPathBasisSystem(Size dim, Size order,
PolynomType polyType) {
        QL_REQUIRE(dim>0, "zero dimension");
        // constant function
        VF_A ret(1, MultiDimConstFct(dim));
        // get single factor basis
        VF_R pathBasis = pathBasisSystem(order, polyType);
        // start with all 0 tuple
        VV tuples(1, std::vector<Size>(dim));
        // add multi-factor terms
        for(Size i=1; i<=order; ++i) {
            tuples = next_order_tuples(tuples);
            // now we have all tuples of order i
            // for each tuple add the corresponding term
            for(Size j=0; j<tuples.size(); ++j) {
                VF_R term(dim);
                for(Size k=0; k<dim; ++k)
                    term[k] = pathBasis[tuples[j][k]];
                ret.push_back(MultiDimFct(term));
            }
        }
        return ret;
    }

} // namespace QuantLib

=============================================================





On Thu, May 13, 2010 at 6:22 PM, Klaus Spanderen <[hidden email]> wrote:

> Hi
>
>> ..The size of basis is not the same for different
>> polynomial types.
>
> I've fixed one issue in the SVN.
>
>> .. I reworked the
>> implementation, so that basis sizes are the same for all polynomial
>> types and no duplicates are created during the process. I would like
>> to submit my code to QuantLib, but not sure how to do it.
>
> Please submit the code to this mailing list or via email. I'm going to merge
> it into the SVN.
>
> thanks in advance
>  Klaus
>

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Kakhkhor Abdijalilov]

Sorry, in my previous message  the formula for size of basis was wrong.
(Dim+Order-1)! / [(Dim-1)! * Order!] is the number of terms of a given
order, not the basis size.

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Klaus Spanderen-2]

Hi

Using your version I'm getting  the following sizes of the basis system.

dim = 1:
order  Mnm      Lgr     Hrm     Hpr     Lgn     Chb     Chb2
0        1           1         1         1         1         1        1
1        2           2         2         2         2         2        2
2        3           3         3         3         3         3        3
3        4           4         4         4         4         4        4
4        5           5         5         5         5         5        5

dim = 2:
order  Mnm      Lgr     Hrm     Hpr     Lgn     Chb     Chb2
0        1            1         1          1         1         1         1
1        3            3         3          3         3         3         3
2        6            6         6          6         6         6         6
3        10          10       10        10       10       10       10
4        15          15       15       15       15        15       15

I guess this is what you intended. The smaller basis sets should increase the
run time efficiency of the algorithm.

best regards
 Klaus

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev

[Kakhkhor Abdijalilov]

Both versions are equivalent,  because default constractor for numeric
types is equivalent to setting to 0. But in the final version of the
code it doesn't matter, since x is overwritten through assignment
inside the nested loop.


On Sat, May 15, 2010 at 5:10 PM, Klaus Spanderen <[hidden email]> wrote:

> Hi
>
> short question, shouldn't it be
>             std::vector<Size> x(dim, 0u);
> instead of
>             std::vector<Size> x(dim);
>
> I guess C++ vector aren't initialized by zero by default(?)
>
> regards
>  Klaus
>
> On Friday 14 May 2010 03:30:19 Kakhkhor Abdijalilov wrote:
>> It seems that increasing the threshold in std::find_if algorithm
>> doesn't solve the problem. I tried it with very high threshold and
>> there are still many more basis functions than needed. The basis set
>> size is the same as it is reported in my original message. I guess
>> "remove-o-zap" works correctly, but something else goes wrong. The
>> correct number of basis functions should be (Dim+Order-1)! / [
>> (Dim-1)! * Order!] and it shouldn't depend on polynomial type.
>>
>> To fix the problem I implemented different basis generation algorithm
>> which is described below.
>>
>> The idea is to build multi-factor basis set by multiplying together
>> single factor basis functions.
>> Let's say we have the single factor basis:
>>
>> F_0(.), F_1(.), F_2(.),...
>>
>> and we want to build a 3-dimensional basis set. In order to build the
>> order N basis set, we need all terms of the type
>> F_i(x)*F_j*(y)*F_k(z)  such that i+j+k <= N. We first compute all
>> tuples (i, j, k) such that i+j+k <= N and then generate the basis set
>> from those tuples. This would avoid duplicates if the tuple set
>> contains only unique elements. The required tuple set can be generated
>> by the following induction algorithm:
>>
>> Assume that we have all order N (i.e. i+j+k = N) tuples. In order to
>> build order N+1 tuples do the following:
>> 1) Start with the set:  (N+1, 0, 0), (0, N+1, 0), (0, 0, N+1).
>> 2) For each tuple (i, j, k) from the order N set add the tuples (i+1,
>> j, k), (i, j+1, k+2) to the order N+1 set.
>> 3) Remove the duplicates from the order N+1 set.
>>
>> The resulting order N+1 set will contain all order N+1 tuples. The
>> function next_order_tuples implements the induction algorithm.
>> Duplicates can be avoided if  we use STL container std::set to store
>> the tuples.
>>
>> The code is below. Please feel free to ask any question you might have.
>>
>> ==================================================
>>
>> // lsmbasissystem.hpp
>>
>> #ifndef quantlib_lsm_basis_system_hpp
>> #define quantlib_lsm_basis_system_hpp
>>
>> #include <ql/qldefines.hpp>
>> #include <ql/math/array.hpp>
>> #include <boost/function.hpp>
>> #include <vector>
>>
>> namespace QuantLib {
>>
>>     class LsmBasisSystem {
>>     public:
>>         enum PolynomType { Monomial, Laguerre, Hermite, Hyperbolic,
>>                            Legendre, Chebyshev, Chebyshev2nd };
>>
>>         static std::vector<boost::function1<Real, Real> >
>>         pathBasisSystem(Size order, PolynomType polyType);
>>
>>         static std::vector<boost::function1<Real, Array> >
>>         multiPathBasisSystem(Size dim, Size order, PolynomType polyType);
>>     };
>>
>>
>> } // namespace QuantLib
>>
>> #endif
>>
>> ==================================================
>>
>> // lsmbasissystem.hpp
>>
>> #include <ql/methods/montecarlo/lsmbasissystem.hpp>
>> #include <ql/math/functional.hpp>
>> #include <ql/math/integrals/gaussianquadratures.hpp>
>> #include <boost/bind.hpp>
>> #include <set>
>> #include <numeric>
>>
>> namespace QuantLib {
>> namespace {
>>
>>     // makes typing a little easier
>>     typedef std::vector<boost::function1<Real, Real> > VF_R;
>>     typedef std::vector<boost::function1<Real, Array> > VF_A;
>>     typedef std::vector<std::vector<Size> > VV;
>>     Real (GaussianOrthogonalPolynomial::*ptr_w)(Size, Real) const =
>>         &GaussianOrthogonalPolynomial::weightedValue;
>>
>>     // pow(x, order)
>>     class MonomialFct : public std::unary_function<Real, Real> {
>>     public:
>>         MonomialFct(Size order): order_(order) {}
>>         inline Real operator()(const Real x) const {
>>             Real ret = 1.0;
>>             for(Size i=0; i<order_; ++i)
>>                 ret *= x;
>>             return ret;
>>         }
>>     private:
>>         const Size order_;
>>     };
>>
>>     /* multiplies [Real -> Real] functors
>>        to create [Array -> Real] functor */
>>     class MultiDimFct : public std::unary_function<Real, Array> {
>>     public:
>>         MultiDimFct(const VF_R b): b_(b) {
>>             QL_REQUIRE(b_.size()>0, "zero size basis");
>>         }
>>         inline Real operator()(const Array& a) const {
>>             #if defined(QL_EXTRA_SAFETY_CHECKS)
>>             QL_REQUIRE(b_.size()==a.size(), "wrong argument size");
>>             #endif
>>             Real ret = b_[0].operator()(a[0]);
>>             for(Size i=1; i<b_.size(); ++i)
>>                 ret *= b_[i].operator()(a[i]);
>>             return ret;
>>         }
>>     private:
>>         const VF_R b_;
>>     };
>>
>>     // constant functor [Real -> Real] (returns 1.0)
>>     class MultiDimConstFct : public std::unary_function<Real, Array> {
>>     public:
>>         MultiDimConstFct(Size dim): dim_(dim) {
>>             QL_REQUIRE(dim>0, "zero dimension");
>>         }
>>         inline Real operator()(const Array& a) const {
>>             #if defined(QL_EXTRA_SAFETY_CHECKS)
>>             QL_REQUIRE(dim_==a.size(), "wrong argument size");
>>             #endif
>>             return 1.0;
>>         }
>>     private:
>>         const Size dim_;
>>     };
>>
>>     // check size and order of tuples
>>     void check_tuples(const VV& v, Size dim, Size order) {
>>         for(Size i=0; i<v.size(); ++i) {
>>             QL_REQUIRE(dim==v[i].size(), "wrong tuple size");
>>             QL_REQUIRE(order==std::accumulate(v[i].begin(), v[i].end(), 0),
>>                 "wrong tuple order");
>>         }
>>     }
>>
>>     // build order N+1 tuples from order N tuples
>>     VV next_order_tuples(const VV& v) {
>>         const Size order = std::accumulate(v[0].begin(), v[0].end(), 0);
>>         const Size dim = v[0].size();
>>
>>         check_tuples(v, dim, order);
>>
>>         // the set of unique tuples
>>         std::set<std::vector<Size> > tuples;
>>
>>         for(Size i=0; i<dim; ++i) {
>>             // order+1 for i-th dimension, the rest is 0
>>             std::vector<Size> x(dim);
>>             x[i] += order+1;
>>             tuples.insert(x);
>>             // increase i-th value in every tuple by 1
>>             for(Size j=0; j<v.size(); ++j) {
>>                 x = v[j];
>>                 x[i] += 1;
>>                 tuples.insert(x);
>>             }
>>         }
>>
>>         VV ret(tuples.begin(), tuples.end());
>>         return ret;
>>     }
>>
>> } // unnamed namespace
>>
>>     // LsmBasisSystem static methods
>>
>>     VF_R LsmBasisSystem::pathBasisSystem(Size order, PolynomType polyType)
>> { VF_R ret(order+1);
>>         for (Size i=0; i<=order; ++i) {
>>             switch (polyType) {
>>               case Monomial:
>>                 ret[i] = MonomialFct(i);
>>                 break;
>>               case Laguerre:
>>                 ret[i] = boost::bind(ptr_w, GaussLaguerrePolynomial(), i,
>> _1); break;
>>               case Hermite:
>>                 ret[i] = boost::bind(ptr_w, GaussHermitePolynomial(), i,
>> _1); break;
>>               case Hyperbolic:
>>                 ret[i] = boost::bind(ptr_w, GaussHyperbolicPolynomial(), i,
>> _1); break;
>>               case Legendre:
>>                 ret[i] = boost::bind(ptr_w, GaussLegendrePolynomial(), i,
>> _1); break;
>>               case Chebyshev:
>>                 ret[i] = boost::bind(ptr_w, GaussChebyshevPolynomial(), i,
>> _1); break;
>>               case Chebyshev2nd:
>>                 ret[i] = boost::bind(ptr_w,
>> GaussChebyshev2ndPolynomial(), i, _1);
>>                 break;
>>               default:
>>                 QL_FAIL("unknown regression type");
>>             }
>>         }
>>         return ret;
>>     }
>>
>>     VF_A LsmBasisSystem::multiPathBasisSystem(Size dim, Size order,
>> PolynomType polyType) {
>>         QL_REQUIRE(dim>0, "zero dimension");
>>         // constant function
>>         VF_A ret(1, MultiDimConstFct(dim));
>>         // get single factor basis
>>         VF_R pathBasis = pathBasisSystem(order, polyType);
>>         // start with all 0 tuple
>>         VV tuples(1, std::vector<Size>(dim));
>>         // add multi-factor terms
>>         for(Size i=1; i<=order; ++i) {
>>             tuples = next_order_tuples(tuples);
>>             // now we have all tuples of order i
>>             // for each tuple add the corresponding term
>>             for(Size j=0; j<tuples.size(); ++j) {
>>                 VF_R term(dim);
>>                 for(Size k=0; k<dim; ++k)
>>                     term[k] = pathBasis[tuples[j][k]];
>>                 ret.push_back(MultiDimFct(term));
>>             }
>>         }
>>         return ret;
>>     }
>>
>> } // namespace QuantLib
>>
>> =============================================================
>>
>> On Thu, May 13, 2010 at 6:22 PM, Klaus Spanderen <[hidden email]> wrote:
>> > Hi
>> >
>> >> ..The size of basis is not the same for different
>> >> polynomial types.
>> >
>> > I've fixed one issue in the SVN.
>> >
>> >> .. I reworked the
>> >> implementation, so that basis sizes are the same for all polynomial
>> >> types and no duplicates are created during the process. I would like
>> >> to submit my code to QuantLib, but not sure how to do it.
>> >
>> > Please submit the code to this mailing list or via email. I'm going to
>> > merge it into the SVN.
>> >
>> > thanks in advance
>> >  Klaus
>>
>> ---------------------------------------------------------------------------
>>---
>>
>> _______________________________________________
>> QuantLib-dev mailing list
>> [hidden email]
>> https://lists.sourceforge.net/lists/listinfo/quantlib-dev
>
>
>

------------------------------------------------------------------------------

_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev