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