1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2 
3 /*
4  Copyright (C) 2015 Andres Hernandez
5 
6  This file is part of QuantLib, a free-software/open-source library
7  for financial quantitative analysts and developers - http://quantlib.org/
8 
9  QuantLib is free software: you can redistribute it and/or modify it
10  under the terms of the QuantLib license.  You should have received a
11  copy of the license along with this program; if not, please email
12  <quantlib-dev@lists.sf.net>. The license is also available online at
13  <http://quantlib.org/license.shtml>.
14 
15  This program is distributed in the hope that it will be useful, but WITHOUT
16  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17  FOR A PARTICULAR PURPOSE.  See the license for more details.
18 */
19 
20 /*! \file isotropicrandomwalk.hpp
21     \brief Isotropic random walk
22 */
23 
24 #ifndef quantlib_isotropic_random_walk_hpp
25 #define quantlib_isotropic_random_walk_hpp
26 
27 #include <ql/mathconstants.hpp>
28 #include <ql/math/randomnumbers/mt19937uniformrng.hpp>
29 #include <ql/math/array.hpp>
30 #include <boost/random/variate_generator.hpp>
31 
32 namespace QuantLib {
33 
34     //! Isotropic random walk
35     /*! A variate is used to draw from a random element of a
36         probability distribution. The draw corresponds to the
37         radius of a d-dimensional sphere. The position on the
38         surface of the d-dimensional sphere is randomly chosen
39         with all points on the surface having the same probability,
40         i.e. all directions are isotropic and the step is randomly
41         drawn from the given variate.
42     */
43     template <class Distribution, class Engine>
44     class IsotropicRandomWalk {
45       public:
46         typedef boost::variate_generator<Engine, Distribution> VariateGenerator;
IsotropicRandomWalk(const Engine & eng,Distribution dist,Size dim,const Array & weights=Array (),unsigned long seed=0)47         IsotropicRandomWalk(const Engine& eng, Distribution dist, Size dim,
48                             const Array& weights = Array(),
49                             unsigned long seed = 0) :
50             variate_(eng, dist), rng_(seed),
51             weights_(weights), dim_(dim) {
52             if (weights_.empty())
53                 weights_ = Array(dim, 1.0);
54             else
55                 QL_REQUIRE(dim_ == weights_.size(), "Invalid weights");
56         }
57         template <class InputIterator>
nextReal(InputIterator first) const58         inline void nextReal(InputIterator first) const {
59             Real radius = variate_();
60             Array::const_iterator weight = weights_.begin();
61             if (dim_ > 1) {
62                 //Isotropic random direction
63                 Real phi = M_PI*rng_.nextReal();
64                 for (Size i = 0; i < dim_ - 2; i++) {
65                     *first++ = radius*cos(phi)*(*weight++);
66                     radius *= sin(phi);
67                     phi = M_PI*rng_.nextReal();
68                 }
69                 *first++ = radius*cos(2.0*phi)*(*weight++);
70                 *first = radius*sin(2.0*phi)*(*weight);
71             }
72             else {
73                 if (rng_.nextReal() < 0.5)
74                     *first = -radius*(*weight);
75                 else
76                     *first = radius*(*weight);
77             }
78         }
setDimension(Size dim)79         inline void setDimension(Size dim) {
80             dim_ = dim;
81             weights_ = Array(dim, 1.0);
82         }
setDimension(Size dim,const Array & weights)83         inline void setDimension(Size dim, const Array& weights) {
84             QL_REQUIRE(dim == weights.size(), "Invalid weights");
85             dim_ = dim;
86             weights_ = weights;
87         }
88         /*!
89         The isotropic random walk will not adjust its draw to be within the lower and upper bounds,
90         but if the limits are provided, they are used to rescale the sphere so as to make it to an
91         ellipsoid, with different radius in different dimensions.
92         */
setDimension(Size dim,const Array & lowerBound,const Array & upperBound)93         inline void setDimension(Size dim,
94             const Array& lowerBound, const Array& upperBound) {
95             QL_REQUIRE(dim == lowerBound.size(),
96                 "Incompatible dimension and lower bound");
97             QL_REQUIRE(dim == upperBound.size(),
98                 "Incompatible dimension and upper bound");
99             //Find largest bound
100             Array bounds = upperBound - lowerBound;
101             Real maxBound = bounds[0];
102             for (Size j = 1; j < dim; j++) {
103                 if (bounds[j] > maxBound) maxBound = bounds[j];
104             }
105             //weights by dimension is the size of the bound
106             //divided by the largest bound
107             maxBound = 1.0 / maxBound;
108             bounds *= maxBound;
109             setDimension(dim, bounds);
110         }
111       protected:
112         mutable VariateGenerator variate_;
113         MersenneTwisterUniformRng rng_;
114         Array weights_;
115         Size dim_;
116     };
117 }
118 #endif
119