JOSS 2026¶
Setup¶
In [1]:
Copied!
try:
import qmcpy
except ImportError:
!pip install -q qmcpy
try:
import seaborn as sns
except ImportError:
!pip install -q seaborn
try:
import qmcpy
except ImportError:
!pip install -q qmcpy
try:
import seaborn as sns
except ImportError:
!pip install -q seaborn
In [2]:
Copied!
import os
import matplotlib
from matplotlib import pyplot
from tueplots import bundles
import seaborn as sns
matplotlib.rcParams['figure.dpi'] = 256
pyplot.rcParams.update(bundles.probnum2025())
COLORS = ['xkcd:purple', 'xkcd:green', 'xkcd:blue', 'xkcd:orange']
MARKERS = ["*","^","o","s"]
MARKERSIZE = 5
OUTDIR = "JOSS2026.outputs"
FIGWIDTH = 500/72
assert os.path.isdir(OUTDIR)
import os
import matplotlib
from matplotlib import pyplot
from tueplots import bundles
import seaborn as sns
matplotlib.rcParams['figure.dpi'] = 256
pyplot.rcParams.update(bundles.probnum2025())
COLORS = ['xkcd:purple', 'xkcd:green', 'xkcd:blue', 'xkcd:orange']
MARKERS = ["*","^","o","s"]
MARKERSIZE = 5
OUTDIR = "JOSS2026.outputs"
FIGWIDTH = 500/72
assert os.path.isdir(OUTDIR)
Listings¶
In [3]:
Copied!
import qmcpy as qp
qp.__version__
import qmcpy as qp
qp.__version__
Out[3]:
'2.1'
In [4]:
Copied!
dnb2 = qp.DigitalNetB2(dimension = 2, randomize = "LMS_DS")
x = dnb2(2**7)
print(x.shape)
dnb2 = qp.DigitalNetB2(dimension = 2, randomize = "LMS_DS")
x = dnb2(2**7)
print(x.shape)
(128, 2)
In [5]:
Copied!
dnb2 = qp.DigitalNetB2(dimension = 2, replications=5, randomize = "LMS_DS")
x = dnb2(2**7)
print(x.shape)
dnb2 = qp.DigitalNetB2(dimension = 2, replications=5, randomize = "LMS_DS")
x = dnb2(2**7)
print(x.shape)
(5, 128, 2)
In [6]:
Copied!
dnb2 = qp.DigitalNetB2( # special case of a digital net
dimension = 12, # monthly monitoring
seed = 7, # for reproducibility
)
asian_option = qp.FinancialOption(
sampler = dnb2 ,
option = "ASIAN",
call_put = "CALL",
asian_mean = "GEOMETRIC",
asian_mean_quadrature_rule = "RIGHT",
volatility = 0.5,
start_price = 30.,
strike_price = 35.,
interest_rate = 0.01, # 1% interest rate
t_final = 1, # 1 year
)
qmc_algorithm = qp.CubQMCNetG(
integrand = asian_option,
abs_tol = 1e-3,
rel_tol = 0,
n_init = 2**8,
)
approx_value,data = qmc_algorithm.integrate()
print(approx_value)
print(data)
dnb2 = qp.DigitalNetB2( # special case of a digital net
dimension = 12, # monthly monitoring
seed = 7, # for reproducibility
)
asian_option = qp.FinancialOption(
sampler = dnb2 ,
option = "ASIAN",
call_put = "CALL",
asian_mean = "GEOMETRIC",
asian_mean_quadrature_rule = "RIGHT",
volatility = 0.5,
start_price = 30.,
strike_price = 35.,
interest_rate = 0.01, # 1% interest rate
t_final = 1, # 1 year
)
qmc_algorithm = qp.CubQMCNetG(
integrand = asian_option,
abs_tol = 1e-3,
rel_tol = 0,
n_init = 2**8,
)
approx_value,data = qmc_algorithm.integrate()
print(approx_value)
print(data)
1.7665677027586677
Data (Data)
solution 1.767
comb_bound_low 1.766
comb_bound_high 1.767
comb_bound_diff 0.002
comb_flags 1
n_total 2^(15)
n 2^(15)
time_integrate 0.015
CubQMCNetG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 0
n_init 2^(8)
n_limit 2^(35)
FinancialOption (AbstractIntegrand)
option ASIAN
call_put CALL
volatility 2^(-1)
start_price 30
strike_price 35
interest_rate 0.010
t_final 1
asian_mean GEOMETRIC
GeometricBrownianMotion (AbstractTrueMeasure)
time_vec [0.083 0.167 0.25 ... 0.833 0.917 1. ]
drift 0.010
diffusion 2^(-2)
mean_gbm [30.025 30.05 30.075 ... 30.251 30.276 30.302]
covariance_gbm [[ 18.978 18.994 19.01 ... 19.121 19.137 19.153]
[ 18.994 38.42 38.452 ... 38.677 38.709 38.742]
[ 19.01 38.452 58.336 ... 58.677 58.726 58.775]
...
[ 19.121 38.677 58.677 ... 211.965 212.141 212.318]
[ 19.137 38.709 58.726 ... 212.141 236.085 236.282]
[ 19.153 38.742 58.775 ... 212.318 236.282 260.787]]
decomp_type PCA
DigitalNetB2 (AbstractLDDiscreteDistribution)
d 12
replications 1
randomize LMS DS
gen_mats_source joe_kuo.6.21201.txt
order RADICAL INVERSE
t 63
alpha 1
n_limit 2^(32)
entropy 7
In [7]:
Copied!
import numpy as np
import scipy.stats
def borehole(t): # t.shape == (..., n, d) for n samples in d dimensions
""" https://www.sfu.ca/~ssurjano/borehole.html """
rw, r, Tu, Hu, Tl, Hl, L, Kw = np.moveaxis(t, -1, 0)
numer = 2 * np.pi * Tu * (Hu-Hl) / np.log(r / rw)
denom = 1 + 2 * L * Tu / (np.log(r / rw) * rw**2 * Kw) + Tu / Tl
y = numer / denom
return y
indep_distribs = [
scipy.stats.norm(loc=0.10, scale=0.0161812),
scipy.stats.lognorm(scale=np.exp(7.71), s=1.0056),
scipy.stats.uniform(loc=63070, scale=115600-63070),
scipy.stats.uniform(loc=990, scale=1110-990),
scipy.stats.uniform(loc=63.1, scale=116-63.1),
scipy.stats.uniform(loc=700, scale=820-700),
scipy.stats.uniform(loc=1120, scale=1680-1120),
scipy.stats.uniform(loc=9855, scale=12045-9855),
]
discrete_distrib = qp.Lattice(dimension = 8, seed = 11)
true_measure = qp.SciPyWrapper(discrete_distrib, indep_distribs)
integrand = qp.CustomFun(true_measure, borehole)
qmc_algorithm = qp.CubQMCBayesLatticeG(
integrand = integrand,
ptransform = "BAKER", # periodization transform
abs_tol = 1e-3,
rel_tol = 1e-5,
error_fun = "BOTH" # abs and rel tol satisfied, default is "EITHER",
)
approx_value,data = qmc_algorithm.integrate()
print(approx_value)
print(data)
import numpy as np
import scipy.stats
def borehole(t): # t.shape == (..., n, d) for n samples in d dimensions
""" https://www.sfu.ca/~ssurjano/borehole.html """
rw, r, Tu, Hu, Tl, Hl, L, Kw = np.moveaxis(t, -1, 0)
numer = 2 * np.pi * Tu * (Hu-Hl) / np.log(r / rw)
denom = 1 + 2 * L * Tu / (np.log(r / rw) * rw**2 * Kw) + Tu / Tl
y = numer / denom
return y
indep_distribs = [
scipy.stats.norm(loc=0.10, scale=0.0161812),
scipy.stats.lognorm(scale=np.exp(7.71), s=1.0056),
scipy.stats.uniform(loc=63070, scale=115600-63070),
scipy.stats.uniform(loc=990, scale=1110-990),
scipy.stats.uniform(loc=63.1, scale=116-63.1),
scipy.stats.uniform(loc=700, scale=820-700),
scipy.stats.uniform(loc=1120, scale=1680-1120),
scipy.stats.uniform(loc=9855, scale=12045-9855),
]
discrete_distrib = qp.Lattice(dimension = 8, seed = 11)
true_measure = qp.SciPyWrapper(discrete_distrib, indep_distribs)
integrand = qp.CustomFun(true_measure, borehole)
qmc_algorithm = qp.CubQMCBayesLatticeG(
integrand = integrand,
ptransform = "BAKER", # periodization transform
abs_tol = 1e-3,
rel_tol = 1e-5,
error_fun = "BOTH" # abs and rel tol satisfied, default is "EITHER",
)
approx_value,data = qmc_algorithm.integrate()
print(approx_value)
print(data)
73.73882925995098
Data (Data)
solution 73.739
comb_bound_low 73.739
comb_bound_high 73.739
comb_bound_diff 2.76e-04
comb_flags 1
n_total 2^(16)
n 2^(16)
time_integrate 0.275
CubQMCBayesLatticeG (AbstractStoppingCriterion)
abs_tol 0.001
rel_tol 1.00e-05
n_init 2^(8)
n_limit 2^(22)
order 2^(1)
CustomFun (AbstractIntegrand)
SciPyWrapper (AbstractTrueMeasure)
Lattice (AbstractLDDiscreteDistribution)
d 2^(3)
replications 1
randomize SHIFT
gen_vec_source kuo.lattice-33002-1024-1048576.9125.txt
order RADICAL INVERSE
n_limit 2^(20)
entropy 11
Points¶
In [8]:
Copied!
ns = np.array([0,2**5,2**6,2**7])
assert len(COLORS)>=(len(ns)-1) and len(MARKERS)>=(len(ns)-1)
nmax = ns.max()
data = [
("IID",qp.IIDStdUniform(2,seed=7)(nmax)),
("LD Lattice Seq.",qp.Lattice(2,seed=7)(nmax)),
("LD Digital Seq.",qp.DigitalNetB2(2,seed=7,randomize="NUS")(nmax)),
("LD Halton Seq.",qp.Halton(2,seed=7,randomize="LMS_DS")(nmax)),
]
nrows = 1
ncols = len(data)
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,sharey=True,sharex=True,figsize=(FIGWIDTH,FIGWIDTH/ncols))
for i,(name,x) in enumerate(data):
for j in range(len(ns)-1):
nmin = ns[j]
nmax = ns[j+1]
ax[i].scatter(x[nmin:nmax,0],x[nmin:nmax,1],marker=MARKERS[j],color=COLORS[j],s=MARKERSIZE);
ax[i].set_xlim([0,1]);
ax[i].set_ylim([0,1]);
ax[i].set_xticks([0,1/4,1/2,3/4,1]); ax[i].set_xticklabels([r"$0$",r"$1/4$",r"$1/2$",r"$3/4$",r"$1$"]);
ax[i].set_yticks([0,1/4,1/2,3/4,1]); ax[i].set_yticklabels([r"$0$",r"$1/4$",r"$1/2$",r"$3/4$",r"$1$"]);
ax[i].grid();
ax[i].set_aspect(1);
ax[i].set_title(name);
fig.savefig(OUTDIR+"/points.pdf");
fig.savefig(OUTDIR+"/points.svg",transparent=True);
fig.savefig(OUTDIR+"/points.png",transparent=True);
ns = np.array([0,2**5,2**6,2**7])
assert len(COLORS)>=(len(ns)-1) and len(MARKERS)>=(len(ns)-1)
nmax = ns.max()
data = [
("IID",qp.IIDStdUniform(2,seed=7)(nmax)),
("LD Lattice Seq.",qp.Lattice(2,seed=7)(nmax)),
("LD Digital Seq.",qp.DigitalNetB2(2,seed=7,randomize="NUS")(nmax)),
("LD Halton Seq.",qp.Halton(2,seed=7,randomize="LMS_DS")(nmax)),
]
nrows = 1
ncols = len(data)
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,sharey=True,sharex=True,figsize=(FIGWIDTH,FIGWIDTH/ncols))
for i,(name,x) in enumerate(data):
for j in range(len(ns)-1):
nmin = ns[j]
nmax = ns[j+1]
ax[i].scatter(x[nmin:nmax,0],x[nmin:nmax,1],marker=MARKERS[j],color=COLORS[j],s=MARKERSIZE);
ax[i].set_xlim([0,1]);
ax[i].set_ylim([0,1]);
ax[i].set_xticks([0,1/4,1/2,3/4,1]); ax[i].set_xticklabels([r"$0$",r"$1/4$",r"$1/2$",r"$3/4$",r"$1$"]);
ax[i].set_yticks([0,1/4,1/2,3/4,1]); ax[i].set_yticklabels([r"$0$",r"$1/4$",r"$1/2$",r"$3/4$",r"$1$"]);
ax[i].grid();
ax[i].set_aspect(1);
ax[i].set_title(name);
fig.savefig(OUTDIR+"/points.pdf");
fig.savefig(OUTDIR+"/points.svg",transparent=True);
fig.savefig(OUTDIR+"/points.png",transparent=True);
Convergence¶
In [9]:
Copied!
tag = "FULL"
force = False
datapath = "%s/%s.npy"%(OUTDIR,tag)
if (not os.path.isfile(datapath)) or force:
d = 12
n_init = 2**8
eps = 10**np.linspace(np.log10(5e-1),np.log10(1e-4),10)
trials = 100
with np.printoptions(formatter={"float":lambda x: "%.1e"%x}): print("eps = %s"%str(eps))
problem = lambda discrete_distrib: \
qp.FinancialOption(discrete_distrib,
option = "ASIAN",
call_put = "CALL",
asian_mean = "GEOMETRIC",
asian_mean_quadrature_rule = "RIGHT",
volatility = 0.5,
start_price = 30.,
strike_price = 35.,
interest_rate = 0.01,
t_final = 1)
exact_value = problem(qp.IIDStdUniform(d)).get_exact_value()
print("exact value = %.3f\n"%exact_value)
algorithms = {
"IID MC": [1e-2, lambda atol: qp.CubMCG(problem(qp.IIDStdUniform(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
"QMC Digital Net'": [0, lambda atol: qp.CubQMCNetG(problem(qp.DigitalNetB2(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
"QMC Lattice": [0, lambda atol: qp.CubQMCLatticeG(problem(qp.Lattice(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
}
keys = list(algorithms.keys())
n_totals = {}
times = {}
approxs = {}
for key,(eps_threshold,sc_constsruct) in algorithms.items():
print(key)
n_totals[key] = np.nan*np.ones((len(eps),trials))
times[key] = np.nan*np.ones((len(eps),trials))
approxs[key] = np.nan*np.ones((len(eps),trials))
for j in range(len(eps)):
if eps[j]<eps_threshold: break
print("\teps[j] = %.1e"%eps[j])
for t in range(trials):
sc = sc_constsruct(eps[j])
sol,data = sc.integrate()
approxs[key][j,t] = sol
n_totals[key][j,t] = data.n_total
times[key][j,t] = data.time_integrate
_data = {
"eps": eps,
"n_totals": n_totals,
"times": times,
"approxs": approxs,
"keys": keys,
"exact_value": exact_value,
}
np.save(datapath,_data)
else:
_data = np.load(datapath,allow_pickle=True).item()
eps = _data["eps"]
n_totals = _data["n_totals"]
times = _data["times"]
approxs = _data["approxs"]
keys = _data["keys"]
exact_value = _data["exact_value"]
tag = "FULL"
force = False
datapath = "%s/%s.npy"%(OUTDIR,tag)
if (not os.path.isfile(datapath)) or force:
d = 12
n_init = 2**8
eps = 10**np.linspace(np.log10(5e-1),np.log10(1e-4),10)
trials = 100
with np.printoptions(formatter={"float":lambda x: "%.1e"%x}): print("eps = %s"%str(eps))
problem = lambda discrete_distrib: \
qp.FinancialOption(discrete_distrib,
option = "ASIAN",
call_put = "CALL",
asian_mean = "GEOMETRIC",
asian_mean_quadrature_rule = "RIGHT",
volatility = 0.5,
start_price = 30.,
strike_price = 35.,
interest_rate = 0.01,
t_final = 1)
exact_value = problem(qp.IIDStdUniform(d)).get_exact_value()
print("exact value = %.3f\n"%exact_value)
algorithms = {
"IID MC": [1e-2, lambda atol: qp.CubMCG(problem(qp.IIDStdUniform(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
"QMC Digital Net'": [0, lambda atol: qp.CubQMCNetG(problem(qp.DigitalNetB2(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
"QMC Lattice": [0, lambda atol: qp.CubQMCLatticeG(problem(qp.Lattice(d)),abs_tol=atol,rel_tol=0,n_init=n_init)],
}
keys = list(algorithms.keys())
n_totals = {}
times = {}
approxs = {}
for key,(eps_threshold,sc_constsruct) in algorithms.items():
print(key)
n_totals[key] = np.nan*np.ones((len(eps),trials))
times[key] = np.nan*np.ones((len(eps),trials))
approxs[key] = np.nan*np.ones((len(eps),trials))
for j in range(len(eps)):
if eps[j]
eps = [5.0e-01 1.9e-01 7.5e-02 2.9e-02 1.1e-02 4.4e-03 1.7e-03 6.6e-04 2.6e-04 1.0e-04] exact value = 1.767 IID MC eps[j] = 5.0e-01 eps[j] = 1.9e-01 eps[j] = 7.5e-02 eps[j] = 2.9e-02 eps[j] = 1.1e-02 QMC Digital Net' eps[j] = 5.0e-01 eps[j] = 1.9e-01 eps[j] = 7.5e-02 eps[j] = 2.9e-02 eps[j] = 1.1e-02 eps[j] = 4.4e-03 eps[j] = 1.7e-03 eps[j] = 6.6e-04 eps[j] = 2.6e-04 eps[j] = 1.0e-04 QMC Lattice eps[j] = 5.0e-01 eps[j] = 1.9e-01 eps[j] = 7.5e-02 eps[j] = 2.9e-02 eps[j] = 1.1e-02 eps[j] = 4.4e-03 eps[j] = 1.7e-03 eps[j] = 6.6e-04 eps[j] = 2.6e-04 eps[j] = 1.0e-04
Convergence Visualization¶
This visualization compares three methods—MC IID, QMC Digital Net, and QMC Lattice—across 100 trials for each error tolerance. The figure includes three subplots:
Left Panel - Sample Complexity: Illustrates the number of samples needed for each error tolerance. Solid lines indicate the median (50th percentile) across trials, while shaded areas represent the inter-decile range (10th to 90th percentiles), showing trial variability. Both axes use a logarithmic scale, with the $x$-axis reversed to highlight lower tolerances on the right. QMC methods require significantly fewer samples than Monte Carlo for equivalent tolerances.
Middle Panel - Computational Time: Shows the wall-clock time to achieve each error tolerance, using median values and 10th-90th percentile bands. Both axes are log-scaled with an inverted $x$-axis. QMC methods demonstrate faster convergence than MC IID.
Right Panel - Error Distribution (Violin Plots): The right panel displays error distributions for a specific tolerance level (eps[jstar], the 5th tolerance). Each method's distribution is represented by a violin plot. The violin shape indicates error probability density across 100 trials; wider sections show more trials within a given error range. Inner box plots display quartiles (25th, 50th, 75th percentiles). The logarithmic $x$-axis accommodates a broad range of error magnitudes.
These plots illustrate each method's consistency and reliability. Narrower distributions suggest more predictable performance, while the plot's position indicates the typical error size. Quasi-Monte Carlo (QMC) methods generally exhibit tighter, more consistent error distributions than Monte Carlo (MC) methods. This demonstrates superior average performance and more reliable convergence.
In [10]:
Copied!
nrows = 1
ncols = 3
alpha = .25
alpha2 = .75
qlow = .1
qhigh = .9
jstar = 4
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,figsize=(7/8*FIGWIDTH,7/8*FIGWIDTH/ncols),sharey=False,sharex=False)
for i,key in enumerate(keys):
n_total = n_totals[key]
time = times[key]
approx = approxs[key]
color = COLORS[i]
ax[0].plot(eps,np.quantile(n_total,.5,axis=1),label=key,color=color,marker=MARKERS[i],alpha=alpha2);
ax[0].fill_between(eps,np.quantile(n_total,qlow,axis=1),np.quantile(n_total,qhigh,axis=1),color=color,alpha=alpha);
ax[1].plot(eps,np.quantile(time[:,:],.5,axis=1),color=color,marker=MARKERS[i],alpha=alpha2);
ax[1].fill_between(eps,np.quantile(time,qlow,axis=1),np.quantile(time,qhigh,axis=1),color=color,alpha=alpha);
errors = np.stack([np.abs(exact_value-approxs[key][jstar]) for key in keys],axis=1)
sns.violinplot(data=errors,ax=ax[2],log_scale=True,orient="h",palette=COLORS[:len(keys)],inner="box");
for i in range(2):
ax[i].set_xscale("log",base=10);
ax[i].set_yscale("log",base=10);
ax[i].set_xlabel("error tolerance");
ax[i].xaxis.set_inverted(True);
for i in range(3):
for spine in ["top","right"]: ax[i].spines[spine].set_visible(False);
ax[2].axvline(eps[jstar],color="xkcd:gray",linestyle=(0,(1,1)),label="absolute error tolerance = %.2e"%eps[jstar]);
ax[0].set_ylabel("samples");
ax[1].set_ylabel("time (sec)");
ax[2].set_xlabel("true error");
ax[2].get_yaxis().set_visible(False);
ax[2].spines["left"].set_visible(False);
fig.legend(frameon=False,loc="upper center",ncol=4,bbox_to_anchor=(.5,1.1));
fig.savefig(OUTDIR+"/stopping_crit.pdf");
fig.savefig(OUTDIR+"/stopping_crit.svg",transparent=True);
fig.savefig(OUTDIR+"/stopping_crit.png",transparent=True);
nrows = 1
ncols = 3
alpha = .25
alpha2 = .75
qlow = .1
qhigh = .9
jstar = 4
fig,ax = pyplot.subplots(nrows=nrows,ncols=ncols,figsize=(7/8*FIGWIDTH,7/8*FIGWIDTH/ncols),sharey=False,sharex=False)
for i,key in enumerate(keys):
n_total = n_totals[key]
time = times[key]
approx = approxs[key]
color = COLORS[i]
ax[0].plot(eps,np.quantile(n_total,.5,axis=1),label=key,color=color,marker=MARKERS[i],alpha=alpha2);
ax[0].fill_between(eps,np.quantile(n_total,qlow,axis=1),np.quantile(n_total,qhigh,axis=1),color=color,alpha=alpha);
ax[1].plot(eps,np.quantile(time[:,:],.5,axis=1),color=color,marker=MARKERS[i],alpha=alpha2);
ax[1].fill_between(eps,np.quantile(time,qlow,axis=1),np.quantile(time,qhigh,axis=1),color=color,alpha=alpha);
errors = np.stack([np.abs(exact_value-approxs[key][jstar]) for key in keys],axis=1)
sns.violinplot(data=errors,ax=ax[2],log_scale=True,orient="h",palette=COLORS[:len(keys)],inner="box");
for i in range(2):
ax[i].set_xscale("log",base=10);
ax[i].set_yscale("log",base=10);
ax[i].set_xlabel("error tolerance");
ax[i].xaxis.set_inverted(True);
for i in range(3):
for spine in ["top","right"]: ax[i].spines[spine].set_visible(False);
ax[2].axvline(eps[jstar],color="xkcd:gray",linestyle=(0,(1,1)),label="absolute error tolerance = %.2e"%eps[jstar]);
ax[0].set_ylabel("samples");
ax[1].set_ylabel("time (sec)");
ax[2].set_xlabel("true error");
ax[2].get_yaxis().set_visible(False);
ax[2].spines["left"].set_visible(False);
fig.legend(frameon=False,loc="upper center",ncol=4,bbox_to_anchor=(.5,1.1));
fig.savefig(OUTDIR+"/stopping_crit.pdf");
fig.savefig(OUTDIR+"/stopping_crit.svg",transparent=True);
fig.savefig(OUTDIR+"/stopping_crit.png",transparent=True);
