1 /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ 2 3 /* 4 Copyright (C) 2003, 2004 Ferdinando Ametrano 5 Copyright (C) 2006 Richard Gould 6 Copyright (C) 2007 Mark Joshi 7 8 This file is part of QuantLib, a free-software/open-source library 9 for financial quantitative analysts and developers - http://quantlib.org/ 10 11 QuantLib is free software: you can redistribute it and/or modify it 12 under the terms of the QuantLib license. You should have received a 13 copy of the license along with this program; if not, please email 14 <quantlib-dev@lists.sf.net>. The license is also available online at 15 <http://quantlib.org/license.shtml>. 16 17 This program is distributed in the hope that it will be useful, but WITHOUT 18 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 19 FOR A PARTICULAR PURPOSE. See the license for more details. 20 */ 21 22 /*! \file sobolrsg.hpp 23 \brief Sobol low-discrepancy sequence generator 24 */ 25 26 #ifndef quantlib_sobol_ld_rsg_hpp 27 #define quantlib_sobol_ld_rsg_hpp 28 29 #include <ql/methods/montecarlo/sample.hpp> 30 #include <boost/cstdint.hpp> 31 #include <vector> 32 33 namespace QuantLib { 34 35 //! Sobol low-discrepancy sequence generator 36 /*! A Gray code counter and bitwise operations are used for very 37 fast sequence generation. 38 39 The implementation relies on primitive polynomials modulo two 40 from the book "Monte Carlo Methods in Finance" by Peter 41 Jäckel. 42 43 21 200 primitive polynomials modulo two are provided in QuantLib. 44 Jäckel has calculated 8 129 334 polynomials: if you need that many 45 dimensions you can replace the primitivepolynomials.cpp file included 46 in QuantLib with the one provided in the CD of the "Monte Carlo 47 Methods in Finance" book. 48 49 The choice of initialization numbers (also know as free direction 50 integers) is crucial for the homogeneity properties of the sequence. 51 Sobol defines two homogeneity properties: Property A and Property A'. 52 53 The unit initialization numbers suggested in "Numerical 54 Recipes in C", 2nd edition, by Press, Teukolsky, Vetterling, 55 and Flannery (section 7.7) fail the test for Property A even 56 for low dimensions. 57 58 Bratley and Fox published coefficients of the free direction 59 integers up to dimension 40, crediting unpublished work of 60 Sobol' and Levitan. See Bratley, P., Fox, B.L. (1988) 61 "Algorithm 659: Implementing Sobol's quasirandom sequence 62 generator," ACM Transactions on Mathematical Software 63 14:88-100. These values satisfy Property A for d<=20 and d = 64 23, 31, 33, 34, 37; Property A' holds for d<=6. 65 66 Jäckel provides in his book (section 8.3) initialization 67 numbers up to dimension 32. Coefficients for d<=8 are the same 68 as in Bradley-Fox, so Property A' holds for d<=6 but Property 69 A holds for d<=32. 70 71 The implementation of Lemieux, Cieslak, and Luttmer includes 72 coefficients of the free direction integers up to dimension 73 360. Coefficients for d<=40 are the same as in Bradley-Fox. 74 For dimension 40<d<=360 the coefficients have 75 been calculated as optimal values based on the "resolution" 76 criterion. See "RandQMC user's guide - A package for 77 randomized quasi-Monte Carlo methods in C," by C. Lemieux, 78 M. Cieslak, and K. Luttmer, version January 13 2004, and 79 references cited there 80 (http://www.math.ucalgary.ca/~lemieux/randqmc.html). 81 The values up to d<=360 has been provided to the QuantLib team by 82 Christiane Lemieux, private communication, September 2004. 83 84 For more info on Sobol' sequences see also "Monte Carlo 85 Methods in Financial Engineering," by P. Glasserman, 2004, 86 Springer, section 5.2.3 87 88 The Joe--Kuo numbers and the Kuo numbers are due to Stephen Joe 89 and Frances Kuo. 90 91 S. Joe and F. Y. Kuo, Constructing Sobol sequences with better 92 two-dimensional projections, preprint Nov 22 2007 93 94 See http://web.maths.unsw.edu.au/~fkuo/sobol/ for more information. 95 96 The Joe-Kuo numbers are available under a BSD-style license 97 available at the above link. 98 99 Note that the Kuo numbers were generated to work with a 100 different ordering of primitive polynomials for the first 40 101 or so dimensions which is why we have the Alternative 102 Primitive Polynomials. 103 104 \test 105 - the correctness of the returned values is tested by 106 reproducing known good values. 107 - the correctness of the returned values is tested by checking 108 their discrepancy against known good values. 109 */ 110 class SobolRsg { 111 public: 112 typedef Sample<std::vector<Real> > sample_type; 113 enum DirectionIntegers { 114 Unit, Jaeckel, SobolLevitan, SobolLevitanLemieux, 115 JoeKuoD5, JoeKuoD6, JoeKuoD7, 116 Kuo, Kuo2, Kuo3 }; 117 /*! \pre dimensionality must be <= PPMT_MAX_DIM */ 118 explicit SobolRsg(Size dimensionality, 119 unsigned long seed = 0, 120 DirectionIntegers directionIntegers = Jaeckel); 121 /*! skip to the n-th sample in the low-discrepancy sequence */ 122 void skipTo(boost::uint_least32_t n); 123 const std::vector<boost::uint_least32_t>& nextInt32Sequence() const; 124 nextSequence() const125 const SobolRsg::sample_type& nextSequence() const { 126 const std::vector<boost::uint_least32_t>& v = nextInt32Sequence(); 127 // normalize to get a double in (0,1) 128 for (Size k=0; k<dimensionality_; ++k) 129 sequence_.value[k] = v[k] * normalizationFactor_; 130 return sequence_; 131 } lastSequence() const132 const sample_type& lastSequence() const { return sequence_; } dimension() const133 Size dimension() const { return dimensionality_; } 134 private: 135 static const int bits_; 136 static const double normalizationFactor_; 137 Size dimensionality_; 138 mutable boost::uint_least32_t sequenceCounter_; 139 mutable bool firstDraw_; 140 mutable sample_type sequence_; 141 mutable std::vector<boost::uint_least32_t> integerSequence_; 142 std::vector<std::vector<boost::uint_least32_t> > directionIntegers_; 143 }; 144 145 } 146 147 #endif 148