QMCPy Documentation¶

Discrete Distribution Class¶

Abstract Discrete Distribution Class¶

class qmcpy.discrete_distribution._discrete_distribution.DiscreteDistribution¶

Discrete Distribution abstract class. DO NOT INSTANTIATE.

gen_samples(*args)¶

ABSTRACT METHOD to generate samples from this discrete distribution.

Parameters: args (tuple) – tuple of positional argument. See implementations for details
Returns: n x d array of samples
Return type: ndarray

pdf(x)¶: ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

plot(dim_x=0, dim_y=1, n=128, point_size=5, color='c', show=True, out=None)¶

Make a scatter plot from samples. Requires dimension >= 2.

Parameters

dim_x (int) – index of first dimension to be plotted on horizontal axis.
dim_y (int) – index of the second dimension to be plotted on vertical axis.
n (int) – number of samples to draw as self.gen_samples(n)
point_size (int) – ax.scatter(…,s=point_size)
color (str) – ax.scatter(…,color=color)
show (bool) – show plot or not?
out (str) – file name to output image. If None, the image is not output

Returns

fig,ax (tuple) from fig,ax = matplotlib.pyplot.subplots(…)

set_seed(seed)¶

ABSTRACT METHOD to reset the seed of the problem.

Parameters: seed (int) – new seed for generator

Lattice¶

class qmcpy.discrete_distribution.lattice.lattice.Lattice(dimension=1, randomize=True, order='natural', seed=None, z_path=None)¶

Quasi-Random Lattice nets in base 2.

>>> l = Lattice(2,seed=7)
>>> l.gen_samples(4)
array([[0.076, 0.78 ],
       [0.576, 0.28 ],
       [0.326, 0.53 ],
       [0.826, 0.03 ]])
>>> l.gen_samples(1)
array([[0.076, 0.78 ]])
>>> l
Lattice (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    order           natural
    seed            7
    mimics          StdUniform
>>> Lattice(dimension=2,randomize=False,order='natural').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])
>>> Lattice(dimension=2,randomize=False,order='linear').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.25, 0.75],
       [0.5 , 0.5 ],
       [0.75, 0.25]])
>>> Lattice(dimension=2,randomize=False,order='mps').gen_samples(4, warn=False)
array([[0.  , 0.  ],
       [0.5 , 0.5 ],
       [0.25, 0.75],
       [0.75, 0.25]])

References

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

[2] F.Y. Kuo & D. Nuyens. Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation, Foundations of Computational Mathematics, 16(6):1631-1696, 2016. springer link: https://link.springer.com/article/10.1007/s10208-016-9329-5 arxiv link: https://arxiv.org/abs/1606.06613

[3] D. Nuyens, The Magic Point Shop of QMC point generators and generating vectors. MATLAB and Python software, 2018. Available from https://people.cs.kuleuven.be/~dirk.nuyens/

[4] Constructing embedded lattice rules for multivariate integration R Cools, FY Kuo, D Nuyens - SIAM J. Sci. Comput., 28(6), 2162-2188.

[5] L’Ecuyer, Pierre & Munger, David. (2015). LatticeBuilder: A General Software Tool for Constructing Rank-1 Lattice Rules. ACM Transactions on Mathematical Software. 42. 10.1145/2754929.

__init__(dimension=1, randomize=True, order='natural', seed=None, z_path=None)¶

Parameters

dimension (int) – dimension of samples
randomize (bool) – If True, apply shift to generated samples. Note: Non-randomized lattice sequence includes the origin.
order (str) – ‘linear’, ‘natural’, or ‘mps’ ordering.
seed (int) – seed the random number generator for reproducibility
z_path (str) – path to generating vector. z_path should be formatted like ‘lattice_vec.3600.20.npy’ where ‘name.d_max.m_max.npy’ and d_max is the maximum dimenion and 2^m_max is the max number samples supported

gen_samples(n=None, n_min=0, n_max=8, warn=True, return_unrandomized=False)¶

Generate lattice samples

Parameters

n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments
n_min (int) – Starting index of sequence.
n_max (int) – Final index of sequence.
return_unrandomized (bool) – return samples without randomization as 2nd return value. Will not be returned if randomize=False.

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

Note

Lattice generates in blocks from 2**m to 2**(m+1) so generating n_min=3 to n_max=9 requires necessarily produces samples from n_min=2 to n_max=16 and automatically subsets. May be inefficient for non-powers-of-2 samples sizes.

pdf(x)¶: pdf of a standard uniform

set_seed(seed)¶: See abstract method.

Sobol’¶

qmcpy.discrete_distribution.sobol.sobol.DigitalNet¶: alias of qmcpy.discrete_distribution.sobol.sobol.Sobol

class qmcpy.discrete_distribution.sobol.sobol.Sobol(dimension=1, randomize='LMS', graycode=False, seed=None, z_path=None, dim0=0)¶

Quasi-Random Sobol nets in base 2.

>>> s = Sobol(2,seed=7)
>>> s.gen_samples(4)
array([[0.673, 0.063],
       [0.398, 0.788],
       [0.913, 0.564],
       [0.125, 0.288]])
>>> s.gen_samples(1)
array([[0.673, 0.063]])
>>> s
Sobol (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    graycode        0
    seed            [61616 58565]
    mimics          StdUniform
    dim0            0
>>> Sobol(dimension=2,randomize=False,graycode=True).gen_samples(n_min=2,n_max=4)
array([[0.75, 0.25],
       [0.25, 0.75]])
>>> Sobol(dimension=2,randomize=False,graycode=False).gen_samples(n_min=2,n_max=4)
array([[0.25, 0.75],
       [0.75, 0.25]])

References

[1] Marius Hofert and Christiane Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. https://CRAN.R-project.org/package=qrng.

[2] Faure, Henri, and Christiane Lemieux. “Implementation of Irreducible Sobol’ Sequences in Prime Power Bases.” Mathematics and Computers in Simulation 161 (2019): 13–22. Crossref. Web.

[3] F.Y. Kuo & D. Nuyens. Application of quasi-Monte Carlo methods to elliptic PDEs with random diffusion coefficients - a survey of analysis and implementation, Foundations of Computational Mathematics, 16(6):1631-1696, 2016. springer link: https://link.springer.com/article/10.1007/s10208-016-9329-5 arxiv link: https://arxiv.org/abs/1606.06613

[4] D. Nuyens, The Magic Point Shop of QMC point generators and generating vectors. MATLAB and Python software, 2018. Available from https://people.cs.kuleuven.be/~dirk.nuyens/

[5] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., … Chintala, S. (2019). PyTorch: An Imperative Style, High-Performance Deep Learning Library. In H. Wallach, H. Larochelle, A. Beygelzimer, F. d extquotesingle Alch'e-Buc, E. Fox, & R. Garnett (Eds.), Advances in Neural Information Processing Systems 32 (pp. 8024–8035). Curran Associates, Inc. Retrieved from http://papers.neurips.cc/paper/9015-pytorch-an-imperative-style-high-performance-deep-learning-library.pdf

[6] I.M. Sobol’, V.I. Turchaninov, Yu.L. Levitan, B.V. Shukhman: “Quasi-Random Sequence Generators” Keldysh Institute of Applied Mathematics, Russian Acamdey of Sciences, Moscow (1992).

[7] Sobol, Ilya & Asotsky, Danil & Kreinin, Alexander & Kucherenko, Sergei. (2011). Construction and Comparison of High-Dimensional Sobol’ Generators. Wilmott. 2011. 10.1002/wilm.10056.

[8] Paul Bratley and Bennett L. Fox. 1988. Algorithm 659: Implementing Sobol’s quasirandom sequence generator. ACM Trans. Math. Softw. 14, 1 (March 1988), 88–100. DOI:https://doi.org/10.1145/42288.214372

__init__(dimension=1, randomize='LMS', graycode=False, seed=None, z_path=None, dim0=0)¶

Parameters

dimension (int) – dimension of samples
randomize (bool) – Apply randomization? True defaults to LMS. Can also explicitly pass in ‘LMS’: Linear matrix scramble with DS ‘DS’: Just Digital Shift
graycode (bool) – indicator to use graycode ordering (True) or natural ordering (False)
seeds (list) – int seed of list of seeds, one for each dimension.
z_path (str) – path to generating matricies. z_path sould be formatted like gen_mat.21201.32.msb.npy with name.d_max.m_max.msb_or_lsb.npy
dim0 (int) – first dimension

gen_samples(n=None, n_min=0, n_max=8, warn=True, return_jlms=False)¶

Generate samples

Parameters

n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments
n_min (int) – Starting index of sequence.
n_max (int) – Final index of sequence.
return_jlms (bool) – return the LMS matrix without digital shift. Only applies when randomize=’LMS’ (the default).

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

pdf(x)¶: pdf of a standard uniform

set_dim0(dim0)¶

Reset the first dimension

Parameters: dim0 (int) – first dimension

set_graycode(graycode)¶

Reset the graycode

Parameters: graycode (bool) – use graycode?

set_randomize(randomize)¶

Reset the randomization

Parameters: randomize (str) – randomization type. Either ‘LMS’: linear matrix scramble with digital shift ‘DS’: just the digital shift

set_seed(seeds)¶

Reset the seeds

Parameters: seeds (int/list/None) – new seeds

Halton¶

class qmcpy.discrete_distribution.halton.halton.Halton(dimension=1, generalize=True, randomize=True, seed=None)¶

Quasi-Random Halton nets.

>>> h = Halton(2,seed=7)
>>> h.gen_samples(4)
array([[0.166, 0.363],
       [0.666, 0.696],
       [0.416, 0.03 ],
       [0.916, 0.474]])
>>> h.gen_samples(1)
array([[0.166, 0.363]])
>>> h
Halton (DiscreteDistribution Object)
    d               2^(1)
    generalize      1
    randomize       1
    seed            7
    mimics          StdUniform

References

[1] Marius Hofert and Christiane Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. https://CRAN.R-project.org/package=qrng.

[2] Owen, A. B. “A randomized Halton algorithm in R,” 2017. arXiv:1706.02808 [stat.CO]

__init__(dimension=1, generalize=True, randomize=True, seed=None)¶

Parameters

dimension (int) – dimension of samples
generalize (bool) – generalize the Halton sequence?
randomize (bool/str) – If False, does not randomize Halton points. If True, will use ‘QRNG’ randomization as in [1]. Supports max dimension 360. You can also set radnomize=’QRNG’ or randomize=’Halton’ to explicitly select a randomization method. Halton OWEN supports up to 1000 dimensions.
seed (int) – seed the random number generator for reproducibility

Note

See References [1] and [2] for specific randomization methods and differences.

gen_samples(n=None, n_min=0, n_max=8, warn=True)¶

Generate samples

Parameters

n (int) – if n is supplied, generate from n_min=0 to n_max=n samples. Otherwise use the n_min and n_max explicitly supplied as the following 2 arguments
n_min (int) – Starting index of sequence.
n_max (int) – Final index of sequence.

Returns

(n_max-n_min) x d (dimension) array of samples

Return type

ndarray

pdf(x)¶: ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

set_seed(seed)¶

ABSTRACT METHOD to reset the seed of the problem.

Parameters: seed (int) – new seed for generator

Korobov¶

class qmcpy.discrete_distribution.korobov.korobov.Korobov(dimension=1, generator=[1], randomize=True, seed=None)¶

Quasi-Random Korobov nets.

>>> k = Korobov(2,generator=[1,3],seed=7)
>>> k.gen_samples(4)
array([[0.982, 0.883],
       [0.232, 0.633],
       [0.482, 0.383],
       [0.732, 0.133]])
>>> k.gen_samples(4)
array([[0.982, 0.883],
       [0.232, 0.633],
       [0.482, 0.383],
       [0.732, 0.133]])
>>> k
Korobov (DiscreteDistribution Object)
    d               2^(1)
    generator       [1 3]
    randomize       1
    seed            7
    mimics          StdUniform
>>> Korobov(2,generator=[3,1],seed=7).gen_samples(4)
array([[0.982, 0.883],
       [0.732, 0.133],
       [0.482, 0.383],
       [0.232, 0.633]])

References

[1] Marius Hofert and Christiane Lemieux (2019). qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. https://CRAN.R-project.org/package=qrng.

__init__(dimension=1, generator=[1], randomize=True, seed=None)¶

Parameters

dimension (int) – dimension of samples
generator (ndarray of ints) – generator in {1,..,n-1} either a vector of length d or a single number (which is appropriately extended)
randomize (bool) – randomize the Korobov sequence? Note: Non-randomized Korobov sequence includes origin
seed (int) – seed the random number generator for reproducibility

gen_samples(n=None, n_min=0, n_max=8, warn=True)¶

Generate samples

Parameters: n (int) – number of samples
Returns: n x d (dimension) array of samples
Return type: ndarray

pdf(x)¶: pdf of a standard uniform

set_seed(seed)¶: See abstract method.

IID Standard Uniform¶

class qmcpy.discrete_distribution.iid_std_uniform.IIDStdUniform(dimension=1, seed=None)¶

A wrapper around NumPy’s IID Standard Uniform generator numpy.random.rand.

>>> dd = IIDStdUniform(dimension=2,seed=7)
>>> dd.gen_samples(4)
array([[0.076, 0.78 ],
       [0.438, 0.723],
       [0.978, 0.538],
       [0.501, 0.072]])
>>> dd
IIDStdUniform (DiscreteDistribution Object)
    d               2^(1)
    seed            7
    mimics          StdUniform

__init__(dimension=1, seed=None)¶

Parameters

dimension (int) – dimension of samples
seed (int) – seed the random number generator for reproducibility

gen_samples(n)¶

Generate samples

Parameters: n (int) – Number of observations to generate
Returns: n x self.d array of samples
Return type: ndarray

pdf(x)¶: ABSTRACT METHOD to evaluate pdf of distribution the samples mimic at locations of x.

True Measure Class¶

Abstract Measure Class¶

class qmcpy.true_measure._true_measure.TrueMeasure¶

True Measure abstract class. DO NOT INSTANTIATE.

gen_samples(*args, **kwargs)¶

Generate samples from the discrete distribution and transform them via the transform method.

Parameters

args (tuple) – positional arguments to the discrete distributions gen_samples method
kwargs (dict) – keyword arguments to the discrete distributions gen_samples method

Returns

n x d matrix of transformed samples

Return type

ndarray

Uniform¶

class qmcpy.true_measure.uniform.Uniform(sampler, lower_bound=0.0, upper_bound=1.0)¶

>>> u = Uniform(Sobol(2,seed=7),lower_bound=[0,.5],upper_bound=[2,3])
>>> u.gen_samples(4)
array([[1.346, 0.658],
       [0.797, 2.471],
       [1.826, 1.909],
       [0.25 , 1.22 ]])
>>> u
Uniform (TrueMeasure Object)
    lower_bound     [0.  0.5]
    upper_bound     [2 3]

__init__(sampler, lower_bound=0.0, upper_bound=1.0)¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
lower_bound (float) – a for Uniform(a,b)
upper_bound (float) – b for Uniform(a,b)

Gaussian¶

class qmcpy.true_measure.gaussian.Gaussian(sampler, mean=0.0, covariance=1.0, decomp_type='PCA')¶

Normal Measure.

>>> g = Gaussian(Sobol(2,seed=7),mean=[1,2],covariance=[[9,4],[4,5]])
>>> g.gen_samples(4)
array([[-1.567,  3.268],
       [ 2.411,  1.376],
       [-2.784, -0.638],
       [ 3.846,  4.804]])
>>> g
Gaussian (TrueMeasure Object)
    mean            [1 2]
    covariance      [[9 4]
                    [4 5]]
    decomp_type     pca

__init__(sampler, mean=0.0, covariance=1.0, decomp_type='PCA')¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
mean (float) – mu for Normal(mu,sigma^2)
covariance (ndarray) – sigma^2 for Normal(mu,sigma^2). A float or d (dimension) vector input will be extended to covariance*eye(d)
decomp_type (str) – method of decomposition either “PCA” for principal component analysis or “Cholesky” for cholesky decomposition.

Brownian Motion¶

class qmcpy.true_measure.brownian_motion.BrownianMotion(sampler, t_final=1, drift=0, decomp_type='PCA')¶

Geometric Brownian Motion.

>>> bm = BrownianMotion(Sobol(4,seed=7),t_final=2,drift=2)
>>> bm.gen_samples(2)
array([[0.583, 0.514, 2.726, 3.967],
       [0.958, 2.973, 3.259, 3.951]])
>>> bm
BrownianMotion (TrueMeasure Object)
    time_vec        [0.5 1.  1.5 2. ]
    drift           2^(1)
    mean            [1. 2. 3. 4.]
    covariance      [[0.5 0.5 0.5 0.5]
                    [0.5 1.  1.  1. ]
                    [0.5 1.  1.5 1.5]
                    [0.5 1.  1.5 2. ]]
    decomp_type     pca

__init__(sampler, t_final=1, drift=0, decomp_type='PCA')¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
t_final (float) – end time for the Brownian Motion.
drift (int) – Gaussian mean is time_vec*drift
decomp_type (str) – method of decomposition either “PCA” for principal component analysis or “Cholesky” for cholesky decomposition.

Lebesgue¶

class qmcpy.true_measure.lebesgue.Lebesgue(sampler)¶

>>> Lebesgue(Gaussian(Sobol(2,seed=7)))
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      1
                       decomp_type     pca
>>> Lebesgue(Uniform(Sobol(2,seed=7)))
Lebesgue (TrueMeasure Object)
    transform       Uniform (TrueMeasure Object)
                       lower_bound     0
                       upper_bound     1

__init__(sampler)¶

Parameters: sampler (TrueMeasure) – A true measure by which to compose a transform.

Kumaraswamy¶

class qmcpy.true_measure.kumaraswamy.Kumaraswamy(sampler, a=2, b=2)¶

For $\boldsymbol{a},\boldsymbol{b}>0$ we have

PDF: $f(\boldsymbol{x}) = \prod_{j=1}^d [a_j b_j x_j^{a_j-1}(1-x_j^{a_j})^{b_j-1}]$

CDF: $F(\boldsymbol{x}) = \prod_{j=1}^d [1-(1-x_j^{a_j})^{b_j}]$

Inverse CDF: $\Psi_j(x_j) = (1-(1-x_j)^{1/b_j})^{1/a_j}$

__init__(sampler, a=2, b=2)¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
a (ndarray) – a
b (ndarray) – b

Integrand Class¶

Abstract Integrand Class¶

class qmcpy.integrand._integrand.Integrand¶

Integrand abstract class. DO NOT INSTANTIATE.

f(x, *args, **kwargs)¶

Evalute transformed integrand based on true measures and discrete distribution

Parameters

x (ndarray) – n x d array of samples from a discrete distribution
*args – other ordered args to g
**kwargs (dict) – other keyword args to g

Returns

length n vector of funciton evaluations

Return type

ndarray

f_periodized(x, ptransform='NONE', *args, **kwargs)¶

Periodized transformed integrand.

Parameters

x (ndarray) – n x d array of samples from a discrete distribution
ptransform (str) – periodization transform.
*args – other ordered args to g
**kwargs (dict) – other keyword args to g

Returns

length n vector of funciton evaluations

Return type

ndarray

g(t, *args, **kwargs)¶

ABSTRACT METHOD for original integrand to be integrated.

Parameters: t (ndarray) – n x d array of samples to be intput into orignal integrand.
Returns: n vector of function evaluations
Return type: ndarray

Keister Function¶

class qmcpy.integrand.keister.Keister(sampler)¶

$f(\boldsymbol{x}) = \pi^{d/2} \cos(\| \boldsymbol{x} \|)$ .

The standard example integrates the Keister integrand with respect to an IID Gaussian distribution with variance 1./2.

>>> k = Keister(Sobol(2,seed=7))
>>> x = k.discrete_distrib.gen_samples(2**10)
>>> y = k.f(x)
>>> y.mean()
1.807...
>>> k.true_measure
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
>>> k = Keister(Gaussian(Sobol(2,seed=7),mean=0,covariance=2))
>>> x = k.discrete_distrib.gen_samples(2**10)
>>> y = k.f(x)
>>> y.mean()
1.808...
>>> yp = k.f_periodized(x,'c2sin')
>>> yp.mean()
1.805...

References

[1] B. D. Keister, Multidimensional Quadrature Algorithms, Computers in Physics, 10, pp. 119-122, 1996.

__init__(sampler)¶

Parameters: sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

g(x)¶

ABSTRACT METHOD for original integrand to be integrated.

Parameters: t (ndarray) – n x d array of samples to be intput into orignal integrand.
Returns: n vector of function evaluations
Return type: ndarray

Custom Function¶

class qmcpy.integrand.custom_fun.CustomFun(true_measure, g)¶

Custom user-supplied function handle.

>>> cf = CustomFun(
...     true_measure = Gaussian(Sobol(2,seed=7),mean=[1,2]),
...     g = lambda x: x[:,0]**2*x[:,1])
>>> x = cf.discrete_distrib.gen_samples(2**10)
>>> y = cf.f(x)
>>> y.mean()
4.009...

__init__(true_measure, g)¶

Parameters

true_measure (TrueMeasure) – a TrueMeasure instance.
g (function) – a function handle.

European Option¶

class qmcpy.integrand.european_option.EuropeanOption(sampler, volatility=0.5, start_price=30, strike_price=35, interest_rate=0, t_final=1, call_put='call')¶

European financial option.

>>> eo = EuropeanOption(Sobol(4,seed=7),call_put='put')
>>> eo
EuropeanOption (Integrand Object)
    volatility      2^(-1)
    call_put        put
    start_price     30
    strike_price    35
    interest_rate   0
>>> x = eo.discrete_distrib.gen_samples(2**12)
>>> y = eo.f(x)
>>> y.mean()
9.210...
>>> eo = EuropeanOption(BrownianMotion(Sobol(4,seed=7),drift=1),call_put='put')
>>> x = eo.discrete_distrib.gen_samples(2**12)
>>> y = eo.f(x)
>>> y.mean()
9.220...

__init__(sampler, volatility=0.5, start_price=30, strike_price=35, interest_rate=0, t_final=1, call_put='call')¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
volatility (float) – sigma, the volatility of the asset
start_price (float) – S(0), the asset value at t=0
strike_price (float) – strike_price, the call/put offer
interest_rate (float) – r, the annual interest rate
t_final (float) – exercise time
call_put (str) – ‘call’ or ‘put’ option

g(x)¶: See abstract method.

get_exact_value()¶

Get the fair price of a European call/put option.

Returns: fair price
Return type: float

Asian Option¶

class qmcpy.integrand.asian_option.AsianOption(sampler, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, t_final=1, call_put='call', mean_type='arithmetic', multi_level_dimensions=None)¶

Asian financial option.

>>> ac = AsianOption(Sobol(4,seed=7))
>>> ac
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dimensions      2^(2)
    dim_fracs       0
>>> x = ac.discrete_distrib.gen_samples(2**10)
>>> y = ac.f(x)
>>> y.mean()
1.773...
>>> level_dims = [2,4,8]
>>> ac2 = AsianOption(Sobol(seed=7),multi_level_dimensions=level_dims)
>>> ac2
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dimensions      [2 4 8]
    dim_fracs       [0. 2. 2.]
>>> y2 = 0
>>> for l in range(len(level_dims)):
...     new_dim = ac2._dim_at_level(l)
...     ac2.true_measure._set_dimension_r(new_dim)
...     x2 = ac2.discrete_distrib.gen_samples(2**10)
...     level_est = ac2.f(x2,l=l).mean()
...     y2 += level_est
>>> y2
1.793...

__init__(sampler, volatility=0.5, start_price=30.0, strike_price=35.0, interest_rate=0.0, t_final=1, call_put='call', mean_type='arithmetic', multi_level_dimensions=None)¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
volatility (float) – sigma, the volatility of the asset
start_price (float) – S(0), the asset value at t=0
strike_price (float) – strike_price, the call/put offer
interest_rate (float) – r, the annual interest rate
t_final (float) – exercise time
mean_type (string) – ‘arithmetic’ or ‘geometric’ mean
multi_level_dimensions (list of ints) – list of dimensions at each level. Leave as None for single-level problems

g(x, l=0)¶: See abstract method.

Multilevel Call Options with Milstein Discretization¶

class qmcpy.integrand.ml_call_options.MLCallOptions(sampler, option='european', volatility=0.2, start_strike_price=100.0, interest_rate=0.05, t_final=1.0)¶

Various call options from finance using Milstein discretization with $2^l$ timesteps on level $l$ .

>>> mlco = MLCallOptions(Sobol(seed=7))
>>> mlco
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
>>> y = 0
>>> for level in range(4):
...     new_dim = mlco._dim_at_level(level)
...     mlco.true_measure._set_dimension_r(new_dim)
...     x = mlco.discrete_distrib.gen_samples(2**10)
...     y += mlco.f(x,l=level).mean()
>>> y
10.406...

References:

[1] M.B. Giles. Improved multilevel Monte Carlo convergence using the Milstein scheme. 343-358, in Monte Carlo and Quasi-Monte Carlo Methods 2006, Springer, 2008. http://people.maths.ox.ac.uk/~gilesm/files/mcqmc06.pdf.

__init__(sampler, option='european', volatility=0.2, start_strike_price=100.0, interest_rate=0.05, t_final=1.0)¶

Parameters

sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform
option_type (str) – type of option in [“European”,”Asian”]
volatility (float) – sigma, the volatility of the asset
start_strike_price (float) – S(0), the asset value at t=0, and K, the strike price. Assume start_price = strike_price
interest_rate (float) – r, the annual interest rate
t_final (float) – exercise time

g(samples, l)¶

Parameters

samples (ndarray) – Gaussian(0,1^2) samples
l (int) – level

Returns

First, an ndarray of length 6 vector of summary statistic sums. Second, a float of cost on this level.

Return type

tuple

get_exact_value()¶: Print exact analytic value, based on s0=k.

Linear Function¶

class qmcpy.integrand.linear0.Linear0(sampler)¶

>>> l = Linear0(Sobol(100,seed=7))
>>> x = l.discrete_distrib.gen_samples(2**10)
>>> y = l.f(x)
>>> y.mean()
-2.654...e-08
>>> ytf = l.f_periodized(x,'C1SIN')
>>> ytf.mean()
8.288...e-10

__init__(sampler)¶

Parameters: sampler (DiscreteDistribution/TrueMeasure) – A discrete distribution from which to transform samples or a true measure by which to compose a transform

g(x)¶

ABSTRACT METHOD for original integrand to be integrated.

Parameters: t (ndarray) – n x d array of samples to be intput into orignal integrand.
Returns: n vector of function evaluations
Return type: ndarray

Stopping Criterion Algorithms¶

Abstract Stopping Criterion Class¶

class qmcpy.stopping_criterion._stopping_criterion.StoppingCriterion(allowed_levels, allowed_distribs)¶

Stopping Criterion abstract class. DO NOT INSTANTIATE.

__init__(allowed_levels, allowed_distribs)¶

Parameters

distribution (DiscreteDistribution) – a DiscreteDistribution
allowed_levels (list) – which integrand types are supported: ‘single’, ‘fixed-multi’, ‘adaptive-multi’
allowed_distribs (list) – list of names (strings) of compatible distributions

integrate()¶

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:

solution (float): approximation to the integral
data (AccumulateData): an AccumulateData object

Return type

tuple

plot(*args, **kwargs)¶: Create a plot relevant to the stopping criterion object.

set_tolerance(*args, **kwargs)¶: ABSTRACT METHOD to reset the absolute tolerance.

Guaranteed Lattice Cubature (qMC)¶

class qmcpy.stopping_criterion.cub_qmc_lattice_g.CubQMCLatticeG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCLatticeG.<lambda>>, check_cone=False, ptransform='Baker')¶

Stopping Criterion quasi-Monte Carlo method using rank-1 Lattices cubature over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Fourier coefficients cone decay assumptions.

>>> k = Keister(Lattice(2,seed=7))
>>> sc = CubQMCLatticeG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.806...
>>> data
Solution: 1.8068
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    order           natural
    seed            7
    mimics          StdUniform
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(10)
    solution        1.807
    error_bound     0.005
    time_integrate  ...

Original Implementation:

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubLattice_g.m

References

[1] Lluis Antoni Jimenez Rugama and Fred J. Hickernell, “Adaptive multidimensional integration based on rank-1 lattices,” Monte Carlo and Quasi-Monte Carlo Methods: MCQMC, Leuven, Belgium, April 2014 (R. Cools and D. Nuyens, eds.), Springer Proceedings in Mathematics and Statistics, vol. 163, Springer-Verlag, Berlin, 2016, arXiv:1411.1966, pp. 407-422.

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee: This algorithm computes the integral of real valued functions in $[0,1]^d$ with a prescribed generalized error tolerance. The Fourier coefficients of the integrand are assumed to be absolutely convergent. If the algorithm terminates without warning messages, the output is given with guarantees under the assumption that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the Fourier coefficients. For integration over domains other than $[0,1]^d$ , this cone condition applies to $f \circ \psi$ (the composition of the functions) where $\psi$ is the transformation function for $[0,1]^d$ to the desired region. For more details on how the cone is defined, please refer to the references below.

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCLatticeG.<lambda>>, check_cone=False, ptransform='Baker')¶

Parameters

integrand (Integrand) – an instance of Integrand
abs_tol (float) – absolute error tolerance
rel_tol (float) – relative error tolerance
n_init (int) – initial number of samples
n_max (int) – maximum number of samples
fudge (function) – positive function multiplying the finite sum of Fast Fourier coefficients specified in the cone of functions
check_cone (boolean) – check if the function falls in the cone

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Guaranteed Sobol’ Cubature (qMC)¶

class qmcpy.stopping_criterion.cub_qmc_sobol_g.CubQMCSobolG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCSobolG.<lambda>>, check_cone=False)¶

Quasi-Monte Carlo method using Sobol’ cubature over the d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Walsh-Fourier coefficients cone decay assumptions.

>>> k = Keister(Sobol(2,seed=7))
>>> sc = CubQMCSobolG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.807...
>>> data
Solution: 1.8074
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    graycode        0
    seed            [61616 58565]
    mimics          StdUniform
    dim0            0
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
LDTransformData (AccumulateData Object)
    n_total         2^(10)
    solution        1.807
    error_bound     0.005
    time_integrate  ...

Original Implementation:

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubSobol_g.m

References

[1] Fred J. Hickernell and Lluis Antoni Jimenez Rugama, Reliable adaptive cubature using digital sequences, 2014. Submitted for publication: arXiv:1410.8615.

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee:: This algorithm computes the integral of real valued functions in $[0,1]^d$ with a prescribed generalized error tolerance. The Fourier coefficients of the integrand are assumed to be absolutely convergent. If the algorithm terminates without warning messages, the output is given with guarantees under the assumption that the integrand lies inside a cone of functions. The guarantee is based on the decay rate of the Fourier coefficients. For integration over domains other than $[0,1]^d$ , this cone condition applies to $f \circ \psi$ (the composition of the functions) where $\psi$ is the transformation function for $[0,1]^d$ to the desired region. For more details on how the cone is defined, please refer to the references below.

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=34359738368.0, fudge=<function CubQMCSobolG.<lambda>>, check_cone=False)¶

Parameters

integrand (Integrand) – an instance of Integrand
abs_tol (float) – absolute error tolerance
rel_tol (float) – relative error tolerance
n_init (int) – initial number of samples
n_max (int) – maximum number of samples
fudge (function) – positive function multiplying the finite sum of Fast Fourier coefficients specified in the cone of functions
check_cone (boolean) – check if the function falls in the cone

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Bayesian Lattice Cubature (qMC)¶

class qmcpy.stopping_criterion.cub_qmc_bayes_lattice_g.CubBayesLatticeG(integrand, abs_tol=0.01, rel_tol=0, n_init=256, n_max=4194304, alpha=0.01, ptransform='C1sin')¶

Stopping criterion for Bayesian Cubature using rank-1 Lattice sequence with guaranteed accuracy over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Bayesian assumptions.

>>> k = Keister(Lattice(2, order='linear', seed=123456789))
>>> sc = CubBayesLatticeG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.808...
>>> data
Solution: 1.8082
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    order           linear
    seed            123456789
    mimics          StdUniform
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubBayesLatticeG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(22)
LDTransformBayesData (AccumulateData Object)
    n_total         256
    solution        1.808
    error_bound     7.36e-04
    time_integrate  ...

Adapted from

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesLattice_g.m

Reference

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee

This algorithm attempts to calculate the integral of function f over the hyperbox [0,1]^d to a prescribed error tolerance tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level, e.g., 99% when alpha=0.5%. If the algorithm terminates without showing any warning messages and provides an answer Q, then the following inequality would be satisfied:

Pr(| Q - I | <= tolfun) = 99%

This Bayesian cubature algorithm guarantees for integrands that are considered to be an instance of a gaussian process that fall in the middle of samples space spanned. Where The sample space is spanned by the covariance kernel parametrized by the scale and shape parameter inferred from the sampled values of the integrand. For more details on how the covariance kernels are defined and the parameters are obtained, please refer to the references below.

integrate()¶

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:

solution (float): approximation to the integral
data (AccumulateData): an AccumulateData object

Return type

tuple

Bayesian Digital Net Cubature (qMC)¶

class qmcpy.stopping_criterion.cub_qmc_bayes_net_g.CubBayesNetG(integrand, abs_tol=0.01, rel_tol=0, n_init=256, n_max=4194304, alpha=0.01)¶

Stopping criterion for Bayesian Cubature using digital net (Sobol) sequence with guaranteed accuracy over a d-dimensional region to integrate within a specified generalized error tolerance with guarantees under Bayesian assumptions.

>>> k = Keister(Sobol(2, seed=123456789))
>>> sc = CubBayesNetG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.811...
>>> data
Solution: 1.8110
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
    d               2^(1)
    randomize       1
    graycode        0
    seed            [77469 77107]
    mimics          StdUniform
    dim0            0
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubBayesNetG (StoppingCriterion Object)
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(22)
LDTransformBayesData (AccumulateData Object)
    n_total         256
    solution        1.811
    error_bound     0.015
    time_integrate  ...

Adapted from

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/cubBayesNet_g.m

Reference

[1] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee

This algorithm attempts to calculate the integral of function f over the hyperbox [0,1]^d to a prescribed error tolerance tolfun:= max(abstol,reltol*| I |) with guaranteed confidence level, e.g., 99% when alpha=0.5%. If the algorithm terminates without showing any warning messages and provides an answer Q, then the following inequality would be satisfied:

Pr(| Q - I | <= tolfun) = 99%

This Bayesian cubature algorithm guarantees for integrands that are considered to be an instance of a gaussian process that fall in the middle of samples space spanned. Where The sample space is spanned by the covariance kernel parametrized by the scale and shape parameter inferred from the sampled values of the integrand. For more details on how the covariance kernels are defined and the parameters are obtained, please refer to the references below.

integrate()¶

ABSTRACT METHOD to determine the number of samples needed to satisfy the tolerance.

Returns

tuple containing:

solution (float): approximation to the integral
data (AccumulateData): an AccumulateData object

Return type

tuple

Multilevel qMC Cubature¶

class qmcpy.stopping_criterion.cub_qmc_ml.CubQMCML(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0)¶

Stopping criterion based on multi-level quasi-Monte Carlo.

>>> mlco = MLCallOptions(Lattice(seed=7))
>>> sc = CubQMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
10.442...
>>> data
Solution: 10.4424
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
Lattice (DiscreteDistribution Object)
    d               2^(6)
    randomize       1
    order           natural
    seed            985802
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     pca
CubQMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    n_max           10000000000
    replications    2^(5)
MLQMCData (AccumulateData Object)
    levels          7
    dimensions      [ 1.  2.  4.  8. 16. 32. 64.]
    n_level         [4096.  512.  256.  256.  256.  256.  256.]
    mean_level      [1.005e+01 1.807e-01 1.033e-01 5.482e-02 2.823e-02 1.397e-02 7.290e-03]
    var_level       [8.376e-05 2.660e-05 1.911e-05 1.594e-05 3.660e-06 1.478e-06 3.424e-07]
    bias_estimate   0.007
    n_total         188416
    time_integrate  ...

References

[1] M.B. Giles and B.J. Waterhouse. ‘Multilevel quasi-Monte Carlo path simulation’. pp.165-181 in Advanced Financial Modelling, in Radon Series on Computational and Applied Mathematics, de Gruyter, 2009. http://people.maths.ox.ac.uk/~gilesm/files/radon.pdf

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, replications=32.0)¶

Parameters

integrand (Integrand) – integrand with multi-level g method
abs_tol (float) – absolute tolerance
alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
rmse_tol (float) – root mean squared error If supplied (not None), then absolute tolerance and alpha are ignored in favor of the rmse tolerance
n_max (int) – maximum number of samples
replications (int) – number of replications on each level

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)¶

See abstract method.

Parameters

integrand (Integrand) – integrand with multi-level g method
abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

CLT qMC Cubature (with Replications)¶

class qmcpy.stopping_criterion.cub_qmc_clt.CubQMCCLT(integrand, abs_tol=0.01, rel_tol=0.0, n_init=256.0, n_max=1073741824, inflate=1.2, alpha=0.01, replications=16.0)¶

Stopping criterion based on Central Limit Theorem for multiple replications.

>>> k = Keister(Lattice(seed=7))
>>> sc = CubQMCCLT(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.380...
>>> data
Solution: 1.3807
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
    d               1
    randomize       1
    order           natural
    seed            1093
    mimics          StdUniform
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubQMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(8)
    n_max           2^(30)
MeanVarDataRep (AccumulateData Object)
    replications    2^(4)
    solution        1.381
    sighat          4.61e-04
    n_total         2^(12)
    error_bound     3.56e-04
    confid_int      [1.38  1.381]
    time_integrate  ...

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=256.0, n_max=1073741824, inflate=1.2, alpha=0.01, replications=16.0)¶

Parameters

integrand (Integrand) – an instance of Integrand
inflate (float) – inflation factor when estimating variance
alpha (float) – significance level for confidence interval
abs_tol (float) – absolute error tolerance
rel_tol (float) – relative error tolerance
n_max (int) – maximum number of samples
replications (int) – number of replications

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Multilevel MC Cubature¶

class qmcpy.stopping_criterion.cub_mc_ml.CubMCML(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2.0, levels_max=10.0, alpha0=- 1.0, beta0=- 1.0, gamma0=- 1.0)¶

Stopping criterion based on multi-level monte carlo.

>>> mlco = MLCallOptions(IIDStdUniform(seed=7))
>>> sc = CubMCML(mlco,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
10.440...
>>> data
Solution: 10.4406
MLCallOptions (Integrand Object)
    option          european
    sigma           0.200
    k               100
    r               0.050
    t               1
    b               85
IIDStdUniform (DiscreteDistribution Object)
    d               2^(6)
    seed            7
    mimics          StdUniform
Gaussian (TrueMeasure Object)
    mean            0
    covariance      1
    decomp_type     pca
CubMCML (StoppingCriterion Object)
    rmse_tol        0.019
    n_init          2^(8)
    levels_min      2^(1)
    levels_max      10
    theta           2^(-2)
MLMCData (AccumulateData Object)
    levels          7
    dimensions      [ 1.  2.  4.  8. 16. 32. 64.]
    n_level         [7.795e+05 1.496e+04 5.916e+03 2.244e+03 7.560e+02 2.770e+02 1.060e+02]
    mean_level      [1.005e+01 1.777e-01 1.071e-01 5.340e-02 2.990e-02 1.167e-02 7.812e-03]
    var_level       [1.959e+02 1.441e-01 4.433e-02 1.093e-02 2.940e-03 7.304e-04 2.261e-04]
    cost_per_sample [ 1.  2.  4.  8. 16. 32. 64.]
    n_total         804118
    alpha           0.942
    beta            1.893
    gamma           1.000
    time_integrate  ...

Original Implementation:

http://people.maths.ox.ac.uk/~gilesm/mlmc/#MATLAB

References

[1] M.B. Giles. ‘Multi-level Monte Carlo path simulation’. Operations Research, 56(3):607-617, 2008. http://people.maths.ox.ac.uk/~gilesm/files/OPRE_2008.pdf.

__init__(integrand, abs_tol=0.05, alpha=0.01, rmse_tol=None, n_init=256.0, n_max=10000000000.0, levels_min=2.0, levels_max=10.0, alpha0=- 1.0, beta0=- 1.0, gamma0=- 1.0)¶

Parameters

integrand (Integrand) – integrand with multi-level g method
abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.
n_init (int) – initial number of samples
n_max (int) – maximum number of samples
levels_min (int) – minimum level of refinement >= 2
levels_max (int) – maximum level of refinement >= Lmin
alpha0 (float) – weak error is O(2^{-alpha0*level})
beta0 (float) – variance is O(2^{-bet0a*level})
gamma0 (float) – sample cost is O(2^{gamma0*level})

Note

if alpha, beta, gamma are not positive, then they will be estimated

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, alpha=0.01, rmse_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
alpha (float) – uncertaintly level. If rmse_tol not supplied, then rmse_tol = abs_tol/norm.ppf(1-alpha/2)
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not. Takes priority over aboluste tolerance and alpha if supplied.

Guaranteed MC Cubature¶

class qmcpy.stopping_criterion.cub_mc_g.CubMCG(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01)¶

Stopping criterion with guaranteed accuracy.

>>> k = Keister(IIDStdUniform(2,seed=7))
>>> sc = CubMCG(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.803...
>>> data
Solution: 1.8039
Keister (Integrand Object)
IIDStdUniform (DiscreteDistribution Object)
    d               2^(1)
    seed            7
    mimics          StdUniform
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubMCG (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
MeanVarData (AccumulateData Object)
    levels          1
    solution        1.804
    n               13473
    n_total         14497
    error_bound     0.050
    confid_int      [1.754 1.854]
    time_integrate  ...

Original Implementation:

https://github.com/GailGithub/GAIL_Dev/blob/master/Algorithms/IntegrationExpectation/meanMC_g.m

References

[1] Fred J. Hickernell, Lan Jiang, Yuewei Liu, and Art B. Owen, “Guaranteed conservative fixed width confidence intervals via Monte Carlo sampling,” Monte Carlo and Quasi-Monte Carlo Methods 2012 (J. Dick, F. Y. Kuo, G. W. Peters, and I. H. Sloan, eds.), pp. 105-128, Springer-Verlag, Berlin, 2014. DOI: 10.1007/978-3-642-41095-6_5

[2] Sou-Cheng T. Choi, Yuhan Ding, Fred J. Hickernell, Lan Jiang, Lluis Antoni Jimenez Rugama, Da Li, Jagadeeswaran Rathinavel, Xin Tong, Kan Zhang, Yizhi Zhang, and Xuan Zhou, GAIL: Guaranteed Automatic Integration Library (Version 2.3) [MATLAB Software], 2019. Available from http://gailgithub.github.io/GAIL_Dev/

Guarantee:: This algorithm attempts to calculate the mean, mu, of a random variable to a prescribed error tolerance, _tol_fun:= max(abstol,reltol*|mu|), with guaranteed confidence level 1-alpha. If the algorithm terminates without showing any warning messages and provides an answer tmu, then the follow inequality would be satisfied: $\P(| mu - tmu | \le tol\_fun) \ge 1 - \alpha$ .

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01)¶

Parameters

integrand (Integrand) – an instance of Integrand
inflate – inflation factor when estimating variance
alpha – significance level for confidence interval
abs_tol – absolute error tolerance
rel_tol – relative error tolerance
n_init – initial number of samples
n_max – maximum number of samples

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

CLT MC Cubature¶

class qmcpy.stopping_criterion.cub_mc_clt.CubMCCLT(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01)¶

Stopping criterion based on the Central Limit Theorem.

>>> k = Keister(IIDStdUniform(2,seed=7))
>>> sc = CubMCCLT(k,abs_tol=.05)
>>> solution,data = sc.integrate()
>>> solution
1.801...
>>> data
Solution: 1.8010
Keister (Integrand Object)
IIDStdUniform (DiscreteDistribution Object)
    d               2^(1)
    seed            7
    mimics          StdUniform
Lebesgue (TrueMeasure Object)
    transform       Gaussian (TrueMeasure Object)
                       mean            0
                       covariance      2^(-1)
                       decomp_type     pca
CubMCCLT (StoppingCriterion Object)
    inflate         1.200
    alpha           0.010
    abs_tol         0.050
    rel_tol         0
    n_init          2^(10)
    n_max           10000000000
MeanVarData (AccumulateData Object)
    levels          1
    solution        1.801
    n               5741
    n_total         6765
    error_bound     0.051
    confid_int      [1.75  1.852]
    time_integrate  ...
>>> ac = AsianOption(IIDStdUniform(),
...     multi_level_dimensions = [2,4,8])
>>> sc = CubMCCLT(ac,abs_tol=.05)
>>> solution,data = sc.integrate()

__init__(integrand, abs_tol=0.01, rel_tol=0.0, n_init=1024.0, n_max=10000000000.0, inflate=1.2, alpha=0.01)¶

Parameters

integrand (Integrand) – an instance of Integrand
inflate (float) – inflation factor when estimating variance
alpha (float) – significance level for confidence interval
abs_tol (float) – absolute error tolerance
rel_tol (float) – relative error tolerance
n_max (int) – maximum number of samples

integrate()¶: See abstract method.

set_tolerance(abs_tol=None, rel_tol=None)¶

See abstract method.

Parameters

abs_tol (float) – absolute tolerance. Reset if supplied, ignored if not.
rel_tol (float) – relative tolerance. Reset if supplied, ignored if not.

Accumulate Data Class¶

Abstract Accumulate Data Class¶

class qmcpy.accumulate_data._accumulate_data.AccumulateData¶

Accumulated Data abstract class. DO NOT INSTANTIATE.

__init__()¶: Initialize data instance

update_data()¶: ABSTRACT METHOD to update the accumulated data.

LD Sequence Transform Data (qMC)¶

class qmcpy.accumulate_data.ld_transform_data.LDTransformData(stopping_crit, integrand, true_measure, discrete_distrib, basis_transform, m_min, m_max, fudge, check_cone, ptransform)¶

Update and store transformation data based on low-discrepancy sequences. See the stopping criterion that utilize this object for references.

__init__(stopping_crit, integrand, true_measure, discrete_distrib, basis_transform, m_min, m_max, fudge, check_cone, ptransform)¶

Parameters

stopping_crit (StoppingCriterion) – a StoppingCriterion instance
integrand (Integrand) – an Integrand instance
true_measure (TrueMeasure) – A TrueMeasure instance
discrete_distrib (DiscreteDistribution) – a DiscreteDistribution instance
basis_transform (method) – Transform ynext, combine with y, and then transform all points. For cub_lattice this is Fast Fourier Transform (FFT). For cub_sobol this is Fast Walsh Transform (FWT)
m_min (int) – initial n == 2^m_min
m_max (int) – max n == 2^m_max
fudge (function) – positive function multiplying the finite sum of basis coefficients specified in the cone of functions
check_cone (boolean) – check if the function falls in the cone

update_data()¶: See abstract method.

Mean Variance qMC Data (for Replications)¶

class qmcpy.accumulate_data.mean_var_data_rep.MeanVarDataRep(stopping_crit, integrand, true_measure, discrete_distrib, n_init, replications)¶

Update and store mean and variance estimates with repliations. See the stopping criterion that utilize this object for references.

__init__(stopping_crit, integrand, true_measure, discrete_distrib, n_init, replications)¶

Parameters

stopping_crit (StoppingCriterion) – a StoppingCriterion instance
integrand (Integrand) – an Integrand instance
true_measure (TrueMeasure) – A TrueMeasure instance
discrete_distrib (DiscreteDistribution) – a DiscreteDistribution instance
n_init (int) – initial number of samples
replications (int) – number of replications

update_data()¶: See abstract method.

Multilevel qMC Data¶

class qmcpy.accumulate_data.mlqmc_data.MLQMCData(stopping_crit, integrand, true_measure, discrete_distrib, n_init, replications)¶

Accumulated data for quasi-random sequence calculations, and store multi-level, multi-replications mean, variance, and cost values. See the stopping criterion that utilize this object for references.

__init__(stopping_crit, integrand, true_measure, discrete_distrib, n_init, replications)¶

Initialize data instance

Parameters

stopping_crit (StoppingCriterion) – a StoppingCriterion instance
integrand (Integrand) – an Integrand instance
true_measure (TrueMeasure) – A TrueMeasure instance
discrete_distrib (DiscreteDistribution) – a DiscreteDistribution instance
replications (int) – number of replications on each level

update_data()¶: See abstract method.

Multilevel MC Data¶

class qmcpy.accumulate_data.mlmc_data.MLMCData(stopping_crit, integrand, true_measure, discrete_distrib, levels_init, n_init, alpha0, beta0, gamma0)¶

Accumulated data for IIDDistribution calculations, and store multi-level mean, variance, and cost values. See the stopping criterion that utilize this object for references.

__init__(stopping_crit, integrand, true_measure, discrete_distrib, levels_init, n_init, alpha0, beta0, gamma0)¶

Initialize data instance

Parameters

stopping_crit (StoppingCriterion) – a StoppingCriterion instance
integrand (Integrand) – an Integrand instance
true_measure (TrueMeasure) – A TrueMeasure instance
discrete_distrib (DiscreteDistribution) – a DiscreteDistribution instance
levels_init (int) – initial number of levels
n_init (int) – initial number of samples per level
alpha0 (float) – weak error is O(2^{-alpha0*level})
beta0 (float) – variance is O(2^{-beta0*level})
gamma0 (float) – sample cost is O(2^{gamma0*level})

update_data()¶: See abstract method.

Mean Variance MC Data¶

class qmcpy.accumulate_data.mean_var_data.MeanVarData(stopping_crit, integrand, true_measure, discrete_distrib, n_init)¶

Update and store mean and variance estimates.

__init__(stopping_crit, integrand, true_measure, discrete_distrib, n_init)¶

Parameters

stopping_crit (StoppingCriterion) – a StoppingCriterion instance
integrand (Integrand) – an Integrand instance
true_measure (TrueMeasure) – A TrueMeasure instance
discrete_distrib (DiscreteDistribution) – a DiscreteDistribution instance
n_init (int) – initial number of samples

update_data()¶: See abstract method.

Utilities¶

qmcpy.util.latnetbuilder_linker.latnetbuilder_linker(lnb_dir='./', out_dir='./', fout_prefix='lnb4qmcpy')¶

Parameters

lnb_dir (str) – relative path to directory where outputMachine.txt is stored e.g. ‘my_lnb/poly_lat/’
out_dir (str) – relative path to directory where output should be stored e.g. ‘my_lnb/poly_lat_qmcpy/’
fout_prefix (str) – start of output file name. e.g. ‘my_poly_lat_vec’

Returns

path to file which can be passed into QMCPy’s Lattice or Sobol’ in order to use: the linked latnetbuilder generating vector/matrix e.g. ‘my_poly_lat_vec.10.16.npy’

Return type

str

Adapted from latnetbuilder parser:: https://github.com/umontreal-simul/latnetbuilder/blob/master/python-wrapper/latnetbuilder/parse_output.py#L74