diff --git a/QuantLib.vcxproj b/QuantLib.vcxproj index eb9d2e40dfb..3ac7544e02c 100644 --- a/QuantLib.vcxproj +++ b/QuantLib.vcxproj @@ -1107,6 +1107,7 @@ + @@ -2308,6 +2309,7 @@ + diff --git a/QuantLib.vcxproj.filters b/QuantLib.vcxproj.filters index db249c711ff..cba8a6501da 100644 --- a/QuantLib.vcxproj.filters +++ b/QuantLib.vcxproj.filters @@ -1236,6 +1236,9 @@ math\randomnumbers + + math\randomnumbers + math\randomnumbers @@ -4958,6 +4961,9 @@ math\matrixutilities + + math\randomnumbers + math\randomnumbers diff --git a/ql/CMakeLists.txt b/ql/CMakeLists.txt index f124cb23c33..da79a0b8790 100644 --- a/ql/CMakeLists.txt +++ b/ql/CMakeLists.txt @@ -413,6 +413,7 @@ set(QL_SOURCES math/polynomialmathfunction.cpp math/primenumbers.cpp math/quadratic.cpp + math/randomnumbers/burley2020sobolrsg.cpp math/randomnumbers/faurersg.cpp math/randomnumbers/haltonrsg.cpp math/randomnumbers/knuthuniformrng.cpp @@ -1515,6 +1516,7 @@ set(QL_HEADERS math/primenumbers.hpp math/quadratic.hpp math/randomnumbers/boxmullergaussianrng.hpp + math/randomnumbers/burley2020sobolrsg.hpp math/randomnumbers/centrallimitgaussianrng.hpp math/randomnumbers/faurersg.hpp math/randomnumbers/haltonrsg.hpp diff --git a/ql/math/randomnumbers/Makefile.am b/ql/math/randomnumbers/Makefile.am index 21faedb2ad0..ac088630f85 100644 --- a/ql/math/randomnumbers/Makefile.am +++ b/ql/math/randomnumbers/Makefile.am @@ -5,6 +5,7 @@ this_includedir=${includedir}/${subdir} this_include_HEADERS = \ all.hpp \ boxmullergaussianrng.hpp \ + burley2020sobolrsg.hpp \ centrallimitgaussianrng.hpp \ faurersg.hpp \ haltonrsg.hpp \ @@ -27,8 +28,9 @@ this_include_HEADERS = \ xoshiro256starstaruniformrng.hpp cpp_files = \ - faurersg.cpp \ - haltonrsg.cpp \ + faurersg.cpp \ + haltonrsg.cpp \ + burley2020sobolrsg.cpp \ knuthuniformrng.cpp \ latticersg.cpp \ latticerules.cpp \ diff --git a/ql/math/randomnumbers/burley2020sobolrsg.cpp b/ql/math/randomnumbers/burley2020sobolrsg.cpp new file mode 100644 index 00000000000..ec441ee21b7 --- /dev/null +++ b/ql/math/randomnumbers/burley2020sobolrsg.cpp @@ -0,0 +1,153 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2023 Peter Caspers + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +#include +#include +#include + +namespace QuantLib { + + Burley2020SobolRsg::Burley2020SobolRsg(Size dimensionality, + unsigned long seed, + SobolRsg::DirectionIntegers directionIntegers, + unsigned long scrambleSeed) + : dimensionality_(dimensionality), seed_(seed), directionIntegers_(directionIntegers), + integerSequence_(dimensionality), sequence_(std::vector(dimensionality), 1.0) { + reset(); + group4Seeds_.resize((dimensionality_ - 1) / 4 + 1); + MersenneTwisterUniformRng mt(scrambleSeed); + for (auto& s : group4Seeds_) { + s = static_cast(mt.nextInt32()); + } + } + + void Burley2020SobolRsg::reset() const { + sobolRsg_ = ext::make_shared(dimensionality_, seed_, directionIntegers_, false); + nextSequenceCounter_ = 0; + } + + const std::vector& Burley2020SobolRsg::skipTo(std::uint32_t n) const { + reset(); + for (Size k = 0; k < n + 1; ++k) { + nextInt32Sequence(); + } + return integerSequence_; + } + + namespace { + + // for reverseBits() see http://graphics.stanford.edu/~seander/bithacks.html#BitReverseTable + + static const std::uint8_t bitReverseTable[] = { + 0u, 128u, 64u, 192u, 32u, 160u, 96u, 224u, 16u, 144u, 80u, 208u, 48u, 176u, + 112u, 240u, 8u, 136u, 72u, 200u, 40u, 168u, 104u, 232u, 24u, 152u, 88u, 216u, + 56u, 184u, 120u, 248u, 4u, 132u, 68u, 196u, 36u, 164u, 100u, 228u, 20u, 148u, + 84u, 212u, 52u, 180u, 116u, 244u, 12u, 140u, 76u, 204u, 44u, 172u, 108u, 236u, + 28u, 156u, 92u, 220u, 60u, 188u, 124u, 252u, 2u, 130u, 66u, 194u, 34u, 162u, + 98u, 226u, 18u, 146u, 82u, 210u, 50u, 178u, 114u, 242u, 10u, 138u, 74u, 202u, + 42u, 170u, 106u, 234u, 26u, 154u, 90u, 218u, 58u, 186u, 122u, 250u, 6u, 134u, + 70u, 198u, 38u, 166u, 102u, 230u, 22u, 150u, 86u, 214u, 54u, 182u, 118u, 246u, + 14u, 142u, 78u, 206u, 46u, 174u, 110u, 238u, 30u, 158u, 94u, 222u, 62u, 190u, + 126u, 254u, 1u, 129u, 65u, 193u, 33u, 161u, 97u, 225u, 17u, 145u, 81u, 209u, + 49u, 177u, 113u, 241u, 9u, 137u, 73u, 201u, 41u, 169u, 105u, 233u, 25u, 153u, + 89u, 217u, 57u, 185u, 121u, 249u, 5u, 133u, 69u, 197u, 37u, 165u, 101u, 229u, + 21u, 149u, 85u, 213u, 53u, 181u, 117u, 245u, 13u, 141u, 77u, 205u, 45u, 173u, + 109u, 237u, 29u, 157u, 93u, 221u, 61u, 189u, 125u, 253u, 3u, 131u, 67u, 195u, + 35u, 163u, 99u, 227u, 19u, 147u, 83u, 211u, 51u, 179u, 115u, 243u, 11u, 139u, + 75u, 203u, 43u, 171u, 107u, 235u, 27u, 155u, 91u, 219u, 59u, 187u, 123u, 251u, + 7u, 135u, 71u, 199u, 39u, 167u, 103u, 231u, 23u, 151u, 87u, 215u, 55u, 183u, + 119u, 247u, 15u, 143u, 79u, 207u, 47u, 175u, 111u, 239u, 31u, 159u, 95u, 223u, + 63u, 191u, 127u, 255u}; + + inline std::uint32_t reverseBits(std::uint32_t x) { + return (bitReverseTable[x & 0xff] << 24) | (bitReverseTable[(x >> 8) & 0xff] << 16) | + (bitReverseTable[(x >> 16) & 0xff] << 8) | (bitReverseTable[(x >> 24) & 0xff]); + } + + inline std::uint32_t laine_karras_permutation(std::uint32_t x, std::uint32_t seed) { + x += seed; + x ^= x * 0x6c50b47cu; + x ^= x * 0xb82f1e52u; + x ^= x * 0xc7afe638u; + x ^= x * 0x8d22f6e6u; + return x; + } + + inline std::uint32_t nested_uniform_scramble(std::uint32_t x, std::uint32_t seed) { + x = reverseBits(x); + x = laine_karras_permutation(x, seed); + x = reverseBits(x); + return x; + } + + // the results depend a lot on the details of the hash_combine() function that is used + // we use hash_combine() calling hash(), hash_mix() as implemented here: + // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L560 + // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/hash.hpp#L115 + // https://github.com/boostorg/container_hash/blob/boost-1.83.0/include/boost/container_hash/detail/hash_mix.hpp#L67 + + inline std::uint64_t local_hash_mix(std::uint64_t x) { + const std::uint64_t m = 0xe9846af9b1a615d; + x ^= x >> 32; + x *= m; + x ^= x >> 32; + x *= m; + x ^= x >> 28; + return x; + } + + inline std::uint64_t local_hash(const std::uint64_t v) { + std::uint64_t seed = 0; + seed = (v >> 32) + local_hash_mix(seed); + seed = (v & 0xFFFFFFFF) + local_hash_mix(seed); + return seed; + } + + inline std::uint64_t local_hash_combine(std::uint64_t x, const uint64_t v) { + return local_hash_mix(x + 0x9e3779b9 + local_hash(v)); + } + } + + const std::vector& Burley2020SobolRsg::nextInt32Sequence() const { + auto n = nested_uniform_scramble(nextSequenceCounter_, group4Seeds_[0]); + const auto& seq = sobolRsg_->skipTo(n); + std::copy(seq.begin(), seq.end(), integerSequence_.begin()); + Size i = 0, group = 0; + do { + std::uint64_t seed = group4Seeds_[group++]; + for (Size g = 0; g < 4 && i < dimensionality_; ++g, ++i) { + seed = local_hash_combine(seed, g); + integerSequence_[i] = + nested_uniform_scramble(integerSequence_[i], static_cast(seed)); + } + } while (i < dimensionality_); + QL_REQUIRE(++nextSequenceCounter_ != 0, + "Burley2020SobolRsg::nextIn32Sequence(): period exceeded"); + return integerSequence_; + } + + const SobolRsg::sample_type& Burley2020SobolRsg::nextSequence() const { + const std::vector& v = nextInt32Sequence(); + // normalize to get a double in (0,1) + for (Size k = 0; k < dimensionality_; ++k) { + sequence_.value[k] = static_cast(v[k]) / 4294967296.0; + } + return sequence_; + } +} diff --git a/ql/math/randomnumbers/burley2020sobolrsg.hpp b/ql/math/randomnumbers/burley2020sobolrsg.hpp new file mode 100644 index 00000000000..445fc3a7679 --- /dev/null +++ b/ql/math/randomnumbers/burley2020sobolrsg.hpp @@ -0,0 +1,64 @@ +/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ + +/* + Copyright (C) 2023 Peter Caspers + + This file is part of QuantLib, a free-software/open-source library + for financial quantitative analysts and developers - http://quantlib.org/ + + QuantLib is free software: you can redistribute it and/or modify it + under the terms of the QuantLib license. You should have received a + copy of the license along with this program; if not, please email + . The license is also available online at + . + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the license for more details. +*/ + +/*! \file burley2020sobolrsg.hpp + \brief scrambled Sobol sequence following Burley, 2020 +*/ + +#ifndef quantlib_burley2020_scrambled_sobolrsg_hpp +#define quantlib_burley2020_scrambled_sobolrsg_hpp + +#include +#include + +namespace QuantLib { + + //! Scrambled sobol sequence according to Burley, 2020 + /*! Reference: Brent Burley: Practical Hash-based Owen Scrambling, + Journal of Computer Graphics Techniques, Vol. 9, No. 4, 2020 */ + class Burley2020SobolRsg { + public: + typedef Sample> sample_type; + explicit Burley2020SobolRsg( + Size dimensionality, + unsigned long seed = 42, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel, + unsigned long scrambleSeed = 43); + const std::vector& skipTo(std::uint32_t n) const; + const std::vector& nextInt32Sequence() const; + const SobolRsg::sample_type& nextSequence() const; + const sample_type& lastSequence() const { return sequence_; } + Size dimension() const { return dimensionality_; } + + private: + void reset() const; + Size dimensionality_; + unsigned long seed_; + SobolRsg::DirectionIntegers directionIntegers_; + mutable ext::shared_ptr sobolRsg_; + mutable std::vector integerSequence_; + mutable sample_type sequence_; + mutable std::uint32_t nextSequenceCounter_; + mutable std::vector group4Seeds_; + }; + +} + + +#endif diff --git a/ql/math/randomnumbers/sobolbrownianbridgersg.cpp b/ql/math/randomnumbers/sobolbrownianbridgersg.cpp index 19198fc3687..8845d16f15d 100644 --- a/ql/math/randomnumbers/sobolbrownianbridgersg.cpp +++ b/ql/math/randomnumbers/sobolbrownianbridgersg.cpp @@ -25,27 +25,31 @@ #include namespace QuantLib { - SobolBrownianBridgeRsg::SobolBrownianBridgeRsg( - Size factors, Size steps, - SobolBrownianGenerator::Ordering ordering, - unsigned long seed, - SobolRsg::DirectionIntegers directionIntegers) - : factors_(factors), steps_(steps), dim_(factors*steps), - seq_(sample_type::value_type(factors*steps), 1.0), - gen_(factors, steps, ordering, seed, directionIntegers) { + + namespace { + void setNextSequence(SobolBrownianGeneratorBase& gen, std::vector& seq) { + gen.nextPath(); + std::vector output(gen.numberOfFactors()); + for (Size i = 0; i < gen.numberOfSteps(); ++i) { + gen.nextStep(output); + std::copy(output.begin(), output.end(), seq.begin() + i * gen.numberOfFactors()); + } + } } + SobolBrownianBridgeRsg::SobolBrownianBridgeRsg(Size factors, + Size steps, + SobolBrownianGenerator::Ordering ordering, + unsigned long seed, + SobolRsg::DirectionIntegers directionIntegers) + : seq_(sample_type::value_type(factors * steps), 1.0), + gen_(factors, steps, ordering, seed, directionIntegers) {} + const SobolBrownianBridgeRsg::sample_type& SobolBrownianBridgeRsg::nextSequence() const { - gen_.nextPath(); - std::vector output(factors_); - for (Size i=0; i < steps_; ++i) { - gen_.nextStep(output); - std::copy(output.begin(), output.end(), - seq_.value.begin()+i*factors_); - } - + setNextSequence(gen_, seq_.value); return seq_; + } const SobolBrownianBridgeRsg::sample_type& @@ -54,6 +58,32 @@ namespace QuantLib { } Size SobolBrownianBridgeRsg::dimension() const { - return dim_; + return gen_.numberOfFactors() * gen_.numberOfSteps(); + } + + Burley2020SobolBrownianBridgeRsg::Burley2020SobolBrownianBridgeRsg( + Size factors, + Size steps, + SobolBrownianGenerator::Ordering ordering, + unsigned long seed, + SobolRsg::DirectionIntegers directionIntegers, + unsigned long scrambleSeed) + : seq_(sample_type::value_type(factors * steps), 1.0), + gen_(factors, steps, ordering, seed, directionIntegers, scrambleSeed) {} + + const Burley2020SobolBrownianBridgeRsg::sample_type& + Burley2020SobolBrownianBridgeRsg::nextSequence() const { + setNextSequence(gen_, seq_.value); + return seq_; + } + + const Burley2020SobolBrownianBridgeRsg::sample_type& + Burley2020SobolBrownianBridgeRsg::lastSequence() const { + return seq_; } + + Size Burley2020SobolBrownianBridgeRsg::dimension() const { + return gen_.numberOfFactors() * gen_.numberOfSteps(); + } + } diff --git a/ql/math/randomnumbers/sobolbrownianbridgersg.hpp b/ql/math/randomnumbers/sobolbrownianbridgersg.hpp index ee699be9690..e877fe34516 100644 --- a/ql/math/randomnumbers/sobolbrownianbridgersg.hpp +++ b/ql/math/randomnumbers/sobolbrownianbridgersg.hpp @@ -45,10 +45,30 @@ namespace QuantLib { Size dimension() const; private: - const Size factors_, steps_, dim_; mutable sample_type seq_; mutable SobolBrownianGenerator gen_; }; + + class Burley2020SobolBrownianBridgeRsg { + public: + typedef Sample > sample_type; + + Burley2020SobolBrownianBridgeRsg( + Size factors, + Size steps, + SobolBrownianGenerator::Ordering ordering = SobolBrownianGenerator::Diagonal, + unsigned long seed = 42, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::JoeKuoD7, + unsigned long scrambleSeed = 43); + + const sample_type& nextSequence() const; + const sample_type& lastSequence() const; + Size dimension() const; + + private: + mutable sample_type seq_; + mutable Burley2020SobolBrownianGenerator gen_; + }; } #endif diff --git a/ql/math/randomnumbers/sobolrsg.cpp b/ql/math/randomnumbers/sobolrsg.cpp index 30264c10f7d..1e0810e843c 100644 --- a/ql/math/randomnumbers/sobolrsg.cpp +++ b/ql/math/randomnumbers/sobolrsg.cpp @@ -78479,10 +78479,14 @@ namespace QuantLib { const double SobolRsg::normalizationFactor_ = 0.5/(1UL<<(SobolRsg::bits_-1)); - SobolRsg::SobolRsg(Size dimensionality, unsigned long seed, DirectionIntegers directionIntegers) + SobolRsg::SobolRsg(Size dimensionality, + unsigned long seed, + DirectionIntegers directionIntegers, + bool useGrayCode) : dimensionality_(dimensionality), sequence_(std::vector(dimensionality), 1.0), integerSequence_(dimensionality, 0), - directionIntegers_(dimensionality, std::vector(bits_)) { + directionIntegers_(dimensionality, std::vector(bits_)), + useGrayCode_(useGrayCode) { QL_REQUIRE(dimensionality>0, "dimensionality must be greater than 0"); @@ -78765,33 +78769,60 @@ namespace QuantLib { */ // initialize the Sobol integer/double vectors - // first draw - for (k=0; k& SobolRsg::skipTo(std::uint_least32_t skip) const { + std::uint_least32_t N = skip + 1; - // Convert to Gray code - std::uint_least32_t G = N ^ (N>>1); - for (Size k=0; k> index & 1) != 0U) - integerSequence_[k] ^= directionIntegers_[k][index]; + if (useGrayCode_) { + unsigned int ops = (unsigned int)(std::log((double)N) / M_LN2) + 1; + + // Convert to Gray code + std::uint_least32_t G = N ^ (N >> 1); + for (Size k = 0; k < dimensionality_; k++) { + integerSequence_[k] = 0; + for (Size index = 0; index < ops; index++) { + if ((G >> index & 1) != 0U) + integerSequence_[k] ^= directionIntegers_[k][index]; + } + } + } else { + std::fill(integerSequence_.begin(), integerSequence_.end(), 0u); + std::uint_least32_t mask = 1; + for (Size index = 0; index < bits_; index++) { + if ((N & mask) != 0U) { + for (Size k = 0; k < dimensionality_; k++) { + integerSequence_[k] ^= directionIntegers_[k][index]; + } + } + mask = mask << 1; } } sequenceCounter_ = skip; + return integerSequence_; } - - const std::vector& SobolRsg::nextInt32Sequence() const { + + if (!useGrayCode_) { + skipTo(sequenceCounter_); + if (firstDraw_) { + firstDraw_ = false; + } else { + ++sequenceCounter_; + QL_REQUIRE(sequenceCounter_ != 0, "period exceeded"); + } + return integerSequence_; + } + if (firstDraw_) { // it was precomputed in the constructor firstDraw_ = false; diff --git a/ql/math/randomnumbers/sobolrsg.hpp b/ql/math/randomnumbers/sobolrsg.hpp index d94581f390b..ee52a5084fa 100644 --- a/ql/math/randomnumbers/sobolrsg.hpp +++ b/ql/math/randomnumbers/sobolrsg.hpp @@ -114,12 +114,21 @@ namespace QuantLib { Unit, Jaeckel, SobolLevitan, SobolLevitanLemieux, JoeKuoD5, JoeKuoD6, JoeKuoD7, Kuo, Kuo2, Kuo3 }; - /*! \pre dimensionality must be <= PPMT_MAX_DIM */ + /*! The so called generating integer is chosen to be \f$\gamma(n) = n\f$ if useGrayCode is set to false and + \f$\gamma(n) = G(n)\f$ where \f$G(n)\f$ is the Gray code of \f$n\f$ otherwise. The Sobol integers are then + constructed using formula 8.20 resp. 8.23, see "Monte Carlo Methods in Finance" by Peter Jäckel. The default + is to use the Gray code since this allows a faster sequence generation. The Burley2020SobolRsg relies on an + underlying SobolRsg not using the Gray code on the other hand due to its specific way of constructing the + integer sequence. + + \pre dimensionality must be <= PPMT_MAX_DIM + */ explicit SobolRsg(Size dimensionality, unsigned long seed = 0, - DirectionIntegers directionIntegers = Jaeckel); + DirectionIntegers directionIntegers = Jaeckel, + bool useGrayCode = true); /*! skip to the n-th sample in the low-discrepancy sequence */ - void skipTo(std::uint_least32_t n); + const std::vector& skipTo(std::uint_least32_t n) const; const std::vector& nextInt32Sequence() const; const SobolRsg::sample_type& nextSequence() const { @@ -140,6 +149,7 @@ namespace QuantLib { mutable sample_type sequence_; mutable std::vector integerSequence_; std::vector> directionIntegers_; + bool useGrayCode_; }; } diff --git a/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.cpp b/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.cpp index 20d7f6a2411..038fa5af280 100644 --- a/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.cpp +++ b/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.cpp @@ -107,13 +107,10 @@ namespace QuantLib { } - SobolBrownianGenerator::SobolBrownianGenerator(Size factors, + SobolBrownianGeneratorBase::SobolBrownianGeneratorBase(Size factors, Size steps, - Ordering ordering, - unsigned long seed, - SobolRsg::DirectionIntegers integers) + Ordering ordering) : factors_(factors), steps_(steps), ordering_(ordering), - generator_(SobolRsg(factors * steps, seed, integers), InverseCumulativeNormal()), bridge_(steps), orderedIndices_(factors, std::vector(steps)), bridgedVariates_(factors, std::vector(steps)) { @@ -133,12 +130,8 @@ namespace QuantLib { } - Real SobolBrownianGenerator::nextPath() { - typedef InverseCumulativeRsg::sample_type - sample_type; - - const sample_type& sample = generator_.nextSequence(); + Real SobolBrownianGeneratorBase::nextPath() { + const auto& sample = nextSequence(); // Brownian-bridge the variates according to the ordered indices for (Size i=0; i >& - SobolBrownianGenerator::orderedIndices() const { + SobolBrownianGeneratorBase::orderedIndices() const { return orderedIndices_; } - std::vector > SobolBrownianGenerator::transform( + std::vector > SobolBrownianGeneratorBase::transform( const std::vector >& variates) { QL_REQUIRE( (variates.size() == factors_*steps_), @@ -190,7 +183,7 @@ namespace QuantLib { return retVal; } - Real SobolBrownianGenerator::nextStep(std::vector& output) { + Real SobolBrownianGeneratorBase::nextStep(std::vector& output) { #if defined(QL_EXTRA_SAFETY_CHECKS) QL_REQUIRE(output.size() == factors_, "size mismatch"); QL_REQUIRE(lastStep_ + Burley2020SobolBrownianGeneratorFactory::create(Size factors, Size steps) const { + return ext::shared_ptr(new Burley2020SobolBrownianGenerator( + factors, steps, ordering_, seed_, integers_, scrambleSeed_)); + } } diff --git a/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp b/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp index f769bf2605d..eacb9abbae3 100644 --- a/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp +++ b/ql/models/marketmodels/browniangenerators/sobolbrowniangenerator.hpp @@ -24,6 +24,7 @@ #include #include #include +#include #include #include #include @@ -34,7 +35,7 @@ namespace QuantLib { /*! Incremental Brownian generator using a Sobol generator, inverse-cumulative Gaussian method, and Brownian bridging. */ - class SobolBrownianGenerator : public BrownianGenerator { + class SobolBrownianGeneratorBase : public BrownianGenerator { public: enum Ordering { Factors, /*!< The variates with the best quality will be @@ -46,13 +47,10 @@ namespace QuantLib { most important factors and the largest steps. */ }; - SobolBrownianGenerator( + SobolBrownianGeneratorBase( Size factors, Size steps, - Ordering ordering, - unsigned long seed = 0, - SobolRsg::DirectionIntegers directionIntegers - = SobolRsg::Jaeckel); + Ordering ordering); Real nextPath() override; Real nextStep(std::vector&) override; @@ -65,10 +63,12 @@ namespace QuantLib { std::vector > transform( const std::vector >& variates); + protected: + virtual const SobolRsg::sample_type& nextSequence() = 0; + private: Size factors_, steps_; Ordering ordering_; - InverseCumulativeRsg generator_; BrownianBridge bridge_; // work variables Size lastStep_ = 0; @@ -76,19 +76,62 @@ namespace QuantLib { std::vector > bridgedVariates_; }; + class SobolBrownianGenerator : public SobolBrownianGeneratorBase { + public: + SobolBrownianGenerator(Size factors, + Size steps, + Ordering ordering, + unsigned long seed = 0, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel); + + private: + const SobolRsg::sample_type& nextSequence() override; + InverseCumulativeRsg generator_; + }; + class SobolBrownianGeneratorFactory : public BrownianGeneratorFactory { public: - SobolBrownianGeneratorFactory( - SobolBrownianGenerator::Ordering ordering, - unsigned long seed = 0, - SobolRsg::DirectionIntegers directionIntegers - = SobolRsg::Jaeckel); + explicit SobolBrownianGeneratorFactory( + SobolBrownianGenerator::Ordering ordering, + unsigned long seed = 0, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel); + ext::shared_ptr create(Size factors, Size steps) const override; + + private: + SobolBrownianGenerator::Ordering ordering_; + unsigned long seed_; + SobolRsg::DirectionIntegers integers_; + }; + + class Burley2020SobolBrownianGenerator : public SobolBrownianGeneratorBase { + public: + Burley2020SobolBrownianGenerator( + Size factors, + Size steps, + Ordering ordering, + unsigned long seed = 42, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel, + unsigned long scrambleSeed = 43); + + private: + const Burley2020SobolRsg::sample_type& nextSequence() override; + InverseCumulativeRsg generator_; + }; + + class Burley2020SobolBrownianGeneratorFactory : public BrownianGeneratorFactory { + public: + explicit Burley2020SobolBrownianGeneratorFactory( + SobolBrownianGenerator::Ordering ordering, + unsigned long seed = 42, + SobolRsg::DirectionIntegers directionIntegers = SobolRsg::Jaeckel, + unsigned long scrambleSeed = 43); ext::shared_ptr create(Size factors, Size steps) const override; private: SobolBrownianGenerator::Ordering ordering_; unsigned long seed_; SobolRsg::DirectionIntegers integers_; + unsigned long scrambleSeed_; }; } diff --git a/test-suite/lowdiscrepancysequences.cpp b/test-suite/lowdiscrepancysequences.cpp index 5d161b03c93..fe68c19ca49 100644 --- a/test-suite/lowdiscrepancysequences.cpp +++ b/test-suite/lowdiscrepancysequences.cpp @@ -18,10 +18,12 @@ FOR A PARTICULAR PURPOSE. See the license for more details. */ +#include "speedlevel.hpp" #include "toplevelfixture.hpp" #include "utilities.hpp" #include #include +#include #include #include #include @@ -1069,6 +1071,73 @@ BOOST_AUTO_TEST_CASE(testSobolSkipping) { } } +BOOST_AUTO_TEST_CASE(testHighDimensionalIntegrals, *precondition(if_speed(Slow))) { + BOOST_TEST_MESSAGE("Testing High Dimensional Integrals"); + + /* We are running "Integration test 1, results for high dimensions" (Figure 9) from: + + Sobol, Asotsky, Kreinin, Kucherenko: Construction and Comparison of High-Dimensional Sobol’ + Generators, available at https://www.broda.co.uk/doc/HD_SobolGenerator.pdf + + We check the error of Kuo1 (using Gray code and sequential numbers) roughly against what + their graph suggests. In addition we check the error of the Burley2020-scrambled version of + Kuo1 against what we experimentally find - the error turns out to be more than one order + better than the unscrambled version. */ + + auto integrand = [](const std::vector& c, const std::vector& x) { + Real p = 1.0; + for (Size i = 0; i < c.size(); ++i) { + p *= 1.0 + c[i] * (x[i] - 0.5); + } + return p; + }; + + Size N = 30031; + + BOOST_TEST_MESSAGE("dimension,Sobol(Gray),Sobol(Seq),Burley2020"); + + std::vector dimension = {1000, 2000, 5000}; + std::vector> expectedOrderOfError = { + {-3.0, -3.0, -4.5}, {-2.5, -2.5, -4.0}, {-2.0, -2.0, -4.0}}; + + for (Size d = 0; d < dimension.size(); ++d) { + + std::vector c1(dimension[d], 0.01); + + SobolRsg s1(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, true); + SobolRsg s2(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, false); + Burley2020SobolRsg s3(dimension[d], 42, SobolRsg::DirectionIntegers::Kuo, 43); + + Real I1 = 0.0, I2 = 0.0, I3 = 0.0; + for (Size i = 0; i < N; ++i) { + I1 += integrand(c1, s1.nextSequence().value) / static_cast(N); + I2 += integrand(c1, s2.nextSequence().value) / static_cast(N); + I3 += integrand(c1, s3.nextSequence().value) / static_cast(N); + } + + Real errOrder1 = std::log10(std::abs(I1 - 1.0)); + Real errOrder2 = std::log10(std::abs(I2 - 1.0)); + Real errOrder3 = std::log10(std::abs(I3 - 1.0)); + + BOOST_TEST_MESSAGE(dimension[d] << "," << errOrder1 << "," << errOrder2 << "," + << errOrder3); + + BOOST_CHECK_MESSAGE(errOrder1 < expectedOrderOfError[d][0], + "order of error for dimension " + std::to_string(dimension[d]) + " is" + + std::to_string(errOrder1) + " expected " + + std::to_string(expectedOrderOfError[d][0])); + BOOST_CHECK_MESSAGE(errOrder2 < expectedOrderOfError[d][1], + "order of error for dimension " + std::to_string(dimension[d]) + " is" + + std::to_string(errOrder2) + " expected " + + std::to_string(expectedOrderOfError[d][1])); + BOOST_CHECK_MESSAGE(errOrder3 < expectedOrderOfError[d][2], + "order of error for dimension " + std::to_string(dimension[d]) + " is" + + std::to_string(errOrder3) + " expected " + + std::to_string(expectedOrderOfError[d][2])); + } +} + + BOOST_AUTO_TEST_SUITE_END() -BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file +BOOST_AUTO_TEST_SUITE_END()