http://quantlib.414.s1.nabble.com/Basis-functions-for-Longstaff-Schwartz-method-tp13169p13172.html
types is equivalent to setting to 0. But in the final version of the
> 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>
>
>