# Integration Examples using QMCPy package¶

In this demo, we show how to use qmcpy for performing numerical multiple integration of two built-in integrands, namely, the Keister function and the Asian call option payoff. To start, we import the qmcpy module and the function arrange() from numpy for generating evenly spaced discrete vectors in the examples.

from qmcpy import *
from numpy import arange

---------------------------------------------------------------------------

ModuleNotFoundError                       Traceback (most recent call last)

<ipython-input-1-5978b102680d> in <module>
----> 1 from qmcpy import *
2 from numpy import arange

ModuleNotFoundError: No module named 'qmcpy'


## Keister Example¶

We recall briefly the mathematical definitions of the Keister function, the Gaussian measure, and the Sobol distribution:

• Keister integrand:

• Gaussian true measure:

• Sobol discrete distribution:

integrand = Keister(Sobol(dimension=3,seed=7))
solution,data = CubQMCSobolG(integrand,abs_tol=.05).integrate()
print(data)

Solution: 2.1661
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
d               3
randomize       1
graycode        0
seed            7
mimics          StdUniform
dim0            0
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        2.166
error_bound     0.011
time_integrate  0.002


## Arithmetic-Mean Asian Put Option: Single Level¶

In this example, we want to estimate the payoff of an European Asian put option that matures at time . The key mathematical entities are defined as follows:

• Stock price at time for is a function of its initial price , interest rate , and volatility :

• Discounted put option payoff is defined as the difference of a fixed strike price and the arithmetic average of the underlying stock prices at discrete time intervals in :

• Brownian motion true measure: for

• Lattice discrete distribution:

integrand = AsianOption(
sampler = IIDStdUniform(dimension=16, seed=7),
volatility = 0.5,
start_price = 30,
strike_price = 25,
interest_rate = 0.01,
mean_type = 'arithmetic')
solution,data = CubMCCLT(integrand, abs_tol=.025).integrate()
print(data)

Solution: 6.2662
AsianOption (Integrand Object)
volatility      2^(-1)
call_put        call
start_price     30
strike_price    25
interest_rate   0.010
mean_type       arithmetic
dimensions      2^(4)
dim_fracs       0
IIDStdUniform (DiscreteDistribution Object)
d               2^(4)
seed            7
mimics          StdUniform
BrownianMotion (TrueMeasure Object)
time_vec        [0.062 0.125 0.188 ... 0.875 0.938 1.   ]
drift           0
mean            [0. 0. 0. ... 0. 0. 0.]
covariance      [[0.062 0.062 0.062 ... 0.062 0.062 0.062]
[0.062 0.125 0.125 ... 0.125 0.125 0.125]
[0.062 0.125 0.188 ... 0.188 0.188 0.188]
...
[0.062 0.125 0.188 ... 0.875 0.875 0.875]
[0.062 0.125 0.188 ... 0.875 0.938 0.938]
[0.062 0.125 0.188 ... 0.875 0.938 1.   ]]
decomp_type     pca
CubMCCLT (StoppingCriterion Object)
inflate         1.200
alpha           0.010
abs_tol         0.025
rel_tol         0
n_init          2^(10)
n_max           10000000000
MeanVarData (AccumulateData Object)
levels          1
solution        6.266
n               874562
n_total         875586
error_bound     0.026
confid_int      [6.241 6.292]
time_integrate  1.934


## Arithmetic-Mean Asian Put Option: Multi-Level¶

This example is similar to the last one except that we use a multi-level method for estimation of the option price. The main idea can be summarized as follows:

integrand = AsianOption(
sampler = IIDStdUniform(seed=7),
volatility = 0.5,
start_price = 30,
strike_price = 25,
interest_rate = 0.01,
mean_type = 'arithmetic',
multi_level_dimensions = [4,8,16])
solution,data = CubMCCLT(integrand, abs_tol=.025).integrate()
print(data)

Solution: 6.2758
AsianOption (Integrand Object)
volatility      2^(-1)
call_put        call
start_price     30
strike_price    25
interest_rate   0.010
mean_type       arithmetic
dimensions      [ 4  8 16]
dim_fracs       [0. 2. 2.]
IIDStdUniform (DiscreteDistribution Object)
d               2^(4)
seed            7
mimics          StdUniform
BrownianMotion (TrueMeasure Object)
time_vec        [0.062 0.125 0.188 ... 0.875 0.938 1.   ]
drift           0
mean            [0. 0. 0. ... 0. 0. 0.]
covariance      [[0.062 0.062 0.062 ... 0.062 0.062 0.062]
[0.062 0.125 0.125 ... 0.125 0.125 0.125]
[0.062 0.125 0.188 ... 0.188 0.188 0.188]
...
[0.062 0.125 0.188 ... 0.875 0.875 0.875]
[0.062 0.125 0.188 ... 0.875 0.938 0.938]
[0.062 0.125 0.188 ... 0.875 0.938 1.   ]]
decomp_type     pca
CubMCCLT (StoppingCriterion Object)
inflate         1.200
alpha           0.010
abs_tol         0.025
rel_tol         0
n_init          2^(10)
n_max           10000000000
MeanVarData (AccumulateData Object)
levels          3
solution        6.276
n               [1642813.   34269.   34594.]
n_total         1714748
error_bound     0.025
confid_int      [6.25  6.301]
time_integrate  1.175


## Keister Example using Bayesian Cubature¶

This examples repeats the Keister using cubBayesLatticeG and cubBayesNetG stopping criterion:

• Keister integrand:

• Gaussian true measure:

• Lattice discrete distribution:

dimension=3
abs_tol=.001
integrand = Keister(Lattice(dimension=dimension, order='linear'))
solution,data = CubBayesLatticeG(integrand,abs_tol=abs_tol).integrate()
print(data)

Solution: 2.1683
Keister (Integrand Object)
Lattice (DiscreteDistribution Object)
d               3
randomize       1
order           linear
seed            240450
mimics          StdUniform
Gaussian (TrueMeasure Object)
mean            0
covariance      2^(-1)
decomp_type     pca
CubBayesLatticeG (StoppingCriterion Object)
abs_tol         0.001
rel_tol         0
n_init          2^(8)
n_max           2^(22)
order           2^(1)
LDTransformBayesData (AccumulateData Object)
n_total         8192
solution        2.168
error_bound     4.91e-04
time_integrate  1.968

dimension=3
abs_tol=.001
integrand = Keister(Sobol(dimension=dimension, graycode=False))
solution,data = CubBayesNetG(integrand,abs_tol=abs_tol).integrate()
print(data)

Solution: 2.1682
Keister (Integrand Object)
Sobol (DiscreteDistribution Object)
d               3
randomize       1
graycode        0
seed            676033
mimics          StdUniform
dim0            0
Gaussian (TrueMeasure Object)
mean            0
covariance      2^(-1)
decomp_type     pca
CubBayesNetG (StoppingCriterion Object)
abs_tol         0.001
rel_tol         0
n_init          2^(8)
n_max           2^(22)
LDTransformBayesData (AccumulateData Object)
n_total         16384
solution        2.168
error_bound     4.99e-04
time_integrate  0.200