1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2006 Klaus Spanderen 5 Copyright (C) 2010 Kakhkhor Abdijalilov 6 7 This file is part of QuantLib, a free-software/open-source library 8 for financial quantitative analysts and developers - http://quantlib.org/ 9 10 QuantLib is free software: you can redistribute it and/or modify it 11 under the terms of the QuantLib license. You should have received a 12 copy of the license along with this program; if not, please email 13 <quantlib-dev@lists.sf.net>. The license is also available online at 14 <http://quantlib.org/license.shtml>. 15 16 This program is distributed in the hope that it will be useful, but WITHOUT 17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 18 FOR A PARTICULAR PURPOSE. See the license for more details. 19 */ 20 21 /*! \file lsmbasissystem.cpp 22 \brief utility classes for longstaff schwartz early exercise Monte Carlo 23 */ 24 // lsmbasissystem.hpp 25 26 #include <ql/math/integrals/gaussianquadratures.hpp> 27 #include <ql/methods/montecarlo/lsmbasissystem.hpp> 28 #include <ql/functional.hpp> 29 #include <set> 30 #include <numeric> 31 32 namespace QuantLib { 33 namespace { 34 35 // makes typing a little easier 36 typedef std::vector<ext::function<Real(Real)> > VF_R; 37 typedef std::vector<ext::function<Real(Array)> > VF_A; 38 typedef std::vector<std::vector<Size> > VV; 39 Real (GaussianOrthogonalPolynomial::*ptr_w)(Size, Real) const = 40 &GaussianOrthogonalPolynomial::weightedValue; 41 42 // pow(x, order) 43 class MonomialFct { 44 public: MonomialFct(Size order)45 explicit MonomialFct(Size order): order_(order) {} operator ()(const Real x) const46 inline Real operator()(const Real x) const { 47 Real ret = 1.0; 48 for(Size i=0; i<order_; ++i) 49 ret *= x; 50 return ret; 51 } 52 private: 53 const Size order_; 54 }; 55 56 /* multiplies [Real -> Real] functors 57 to create [Array -> Real] functor */ 58 class MultiDimFct { 59 public: MultiDimFct(const VF_R & b)60 explicit MultiDimFct(const VF_R& b): b_(b) { 61 QL_REQUIRE(!b_.empty(), "zero size basis"); 62 } operator ()(const Array & a) const63 inline Real operator()(const Array& a) const { 64 #if defined(QL_EXTRA_SAFETY_CHECKS) 65 QL_REQUIRE(b_.size()==a.size(), "wrong argument size"); 66 #endif 67 Real ret = b_[0](a[0]); 68 for(Size i=1; i<b_.size(); ++i) 69 ret *= b_[i](a[i]); 70 return ret; 71 } 72 private: 73 const VF_R b_; 74 }; 75 76 // check size and order of tuples check_tuples(const VV & v,Size dim,Size order)77 void check_tuples(const VV& v, Size dim, Size order) { 78 for(Size i=0; i<v.size(); ++i) { 79 QL_REQUIRE(dim==v[i].size(), "wrong tuple size"); 80 QL_REQUIRE(order == std::accumulate(v[i].begin(), v[i].end(), 0UL), 81 "wrong tuple order"); 82 } 83 } 84 85 // build order N+1 tuples from order N tuples next_order_tuples(const VV & v)86 VV next_order_tuples(const VV& v) { 87 const Size order = std::accumulate(v[0].begin(), v[0].end(), 0UL); 88 const Size dim = v[0].size(); 89 90 check_tuples(v, dim, order); 91 92 // the set of unique tuples 93 std::set<std::vector<Size> > tuples; 94 std::vector<Size> x; 95 for(Size i=0; i<dim; ++i) { 96 // increase i-th value in every tuple by 1 97 for(Size j=0; j<v.size(); ++j) { 98 x = v[j]; 99 x[i] += 1; 100 tuples.insert(x); 101 } 102 } 103 104 VV ret(tuples.begin(), tuples.end()); 105 return ret; 106 } 107 } 108 109 // LsmBasisSystem static methods 110 pathBasisSystem(Size order,PolynomType polyType)111 VF_R LsmBasisSystem::pathBasisSystem(Size order, PolynomType polyType) { 112 using namespace ext::placeholders; 113 VF_R ret(order+1); 114 for (Size i=0; i<=order; ++i) { 115 switch (polyType) { 116 case Monomial: 117 ret[i] = MonomialFct(i); 118 break; 119 case Laguerre: 120 ret[i] = ext::bind(ptr_w, GaussLaguerrePolynomial(), i, _1); 121 break; 122 case Hermite: 123 ret[i] = ext::bind(ptr_w, GaussHermitePolynomial(), i, _1); 124 break; 125 case Hyperbolic: 126 ret[i] = ext::bind(ptr_w, GaussHyperbolicPolynomial(), i, _1); 127 break; 128 case Legendre: 129 ret[i] = ext::bind(ptr_w, GaussLegendrePolynomial(), i, _1); 130 break; 131 case Chebyshev: 132 ret[i] = ext::bind(ptr_w, GaussChebyshevPolynomial(), i, _1); 133 break; 134 case Chebyshev2nd: 135 ret[i] = ext::bind(ptr_w,GaussChebyshev2ndPolynomial(),i, _1); 136 break; 137 default: 138 QL_FAIL("unknown regression type"); 139 } 140 } 141 return ret; 142 } 143 multiPathBasisSystem(Size dim,Size order,PolynomType polyType)144 VF_A LsmBasisSystem::multiPathBasisSystem(Size dim, Size order, 145 PolynomType polyType) { 146 QL_REQUIRE(dim>0, "zero dimension"); 147 // get single factor basis 148 VF_R pathBasis = pathBasisSystem(order, polyType); 149 VF_A ret; 150 // 0-th order term 151 VF_R term(dim, pathBasis[0]); 152 ret.push_back(MultiDimFct(term)); 153 // start with all 0 tuple 154 VV tuples(1, std::vector<Size>(dim)); 155 // add multi-factor terms 156 for(Size i=1; i<=order; ++i) { 157 tuples = next_order_tuples(tuples); 158 // now we have all tuples of order i 159 // for each tuple add the corresponding term 160 for(Size j=0; j<tuples.size(); ++j) { 161 for(Size k=0; k<dim; ++k) 162 term[k] = pathBasis[tuples[j][k]]; 163 ret.push_back(MultiDimFct(term)); 164 } 165 } 166 return ret; 167 } 168 } 169