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/math/array.hpp>
28#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
29#include <ql/mathconstants.hpp>
30#include <utility>
31
32namespace 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 IsotropicRandomWalk(Engine eng,
47 Distribution dist,
48 Size dim,
49 Array weights = Array(),
50 unsigned long seed = 0)
51 : engine_(std::move(eng)), distribution_(std::move(dist)), rng_(seed), weights_(std::move(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>
58 inline void nextReal(InputIterator first) {
59 Real radius = distribution_(engine_);
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(x: phi)*(*weight++);
66 radius *= sin(x: phi);
67 phi = M_PI*rng_.nextReal();
68 }
69 *first++ = radius*cos(x: 2.0*phi)*(*weight++);
70 *first = radius*sin(x: 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 }
79 inline void setDimension(Size dim) {
80 dim_ = dim;
81 weights_ = Array(dim, 1.0);
82 }
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 */
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 Engine engine_;
113 Distribution distribution_;
114 MersenneTwisterUniformRng rng_;
115 Array weights_;
116 Size dim_;
117 };
118}
119#endif
120

source code of quantlib/ql/experimental/math/isotropicrandomwalk.hpp