Posted by
Kakhkhor Abdijalilov on
URL: http://quantlib.414.s1.nabble.com/Basis-functions-for-Longstaff-Schwartz-method-tp13169p13171.html
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