Control Variates in QMCPy¶

This notebook demonstrates QMCPy's current support for control variates.

Setup¶

In [1]:

from qmcpy import *
from numpy import *
from qmcpy import *
from numpy import *

from qmcpy import * from numpy import * from qmcpy import * from numpy import *

In [2]:

from matplotlib import pyplot
%matplotlib inline
size = 20
pyplot.rc('font', size=size)          # controls default text sizes
pyplot.rc('axes', titlesize=size)     # fontsize of the axes title
pyplot.rc('axes', labelsize=size)     # fontsize of the x and y labels
pyplot.rc('xtick', labelsize=size)    # fontsize of the tick labels
pyplot.rc('ytick', labelsize=size)    # fontsize of the tick labels
pyplot.rc('legend', fontsize=size)    # legend fontsize
pyplot.rc('figure', titlesize=size)   # fontsize of the figure title

from matplotlib import pyplot %matplotlib inline size = 20 pyplot.rc('font', size=size) # controls default text sizes pyplot.rc('axes', titlesize=size) # fontsize of the axes title pyplot.rc('axes', labelsize=size) # fontsize of the x and y labels pyplot.rc('xtick', labelsize=size) # fontsize of the tick labels pyplot.rc('ytick', labelsize=size) # fontsize of the tick labels pyplot.rc('legend', fontsize=size) # legend fontsize pyplot.rc('figure', titlesize=size) # fontsize of the figure title

In [3]:

def compare(problem,discrete_distrib,stopping_crit,abs_tol):
  g1,cvs,cvmus = problem(discrete_distrib)
  sc1 = stopping_crit(g1,abs_tol=abs_tol)
  name = type(sc1).__name__
  print('Stopping Criterion: %-15s absolute tolerance: %-5.1e'%(name,abs_tol))
  sol,data = sc1.integrate()
  print('\tW CV:  Solution %-10.2f time %-10.2f samples %.1e'%(sol,data.time_integrate,data.n_total))
  sc1 = stopping_crit(g1,abs_tol=abs_tol,control_variates=cvs,control_variate_means=cvmus)
  solcv,datacv = sc1.integrate()
  print('\tWO CV: Solution %-10.2f time %-10.2f samples %.1e'%(solcv,datacv.time_integrate,datacv.n_total))
  print('\tControl variates took %.1f%% the time and %.1f%% the samples\n'%\
        (100*datacv.time_integrate/data.time_integrate,100*datacv.n_total/data.n_total))

def compare(problem,discrete_distrib,stopping_crit,abs_tol): g1,cvs,cvmus = problem(discrete_distrib) sc1 = stopping_crit(g1,abs_tol=abs_tol) name = type(sc1).__name__ print('Stopping Criterion: %-15s absolute tolerance: %-5.1e'%(name,abs_tol)) sol,data = sc1.integrate() print('\tW CV: Solution %-10.2f time %-10.2f samples %.1e'%(sol,data.time_integrate,data.n_total)) sc1 = stopping_crit(g1,abs_tol=abs_tol,control_variates=cvs,control_variate_means=cvmus) solcv,datacv = sc1.integrate() print('\tWO CV: Solution %-10.2f time %-10.2f samples %.1e'%(solcv,datacv.time_integrate,datacv.n_total)) print('\tControl variates took %.1f%% the time and %.1f%% the samples\n'%\ (100*datacv.time_integrate/data.time_integrate,100*datacv.n_total/data.n_total))

Problem 1: Polynomial Function¶

We will integrate $$g(t) = 10t_1-5t_2^2+2t_3^3$$ with true measure $\mathcal{U}[0,2]^3$ and control variates $$\hat{g}_1(t) = t_1$$ and $$\hat{g}_2(t) = t_2^2$$ using the same true measure.

In [4]:

# parameters
def poly_problem(discrete_distrib):
  g1 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: 10*t[...,0]-5*t[...,1]**2+t[...,2]**3)
  cv1 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: t[...,0])
  cv2 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: t[...,1]**2)
  return g1,[cv1,cv2],[1,4/3]
compare(poly_problem,IIDStdUniform(3,seed=7),CubMCCLT,abs_tol=1e-2)
compare(poly_problem,IIDStdUniform(3,seed=7),CubMCCLT,abs_tol=1e-2)
compare(poly_problem,Sobol(3,seed=7),CubQMCSobolG,abs_tol=1e-8)

# parameters def poly_problem(discrete_distrib): g1 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: 10*t[...,0]-5*t[...,1]**2+t[...,2]**3) cv1 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: t[...,0]) cv2 = CustomFun(Uniform(discrete_distrib,0,2),lambda t: t[...,1]**2) return g1,[cv1,cv2],[1,4/3] compare(poly_problem,IIDStdUniform(3,seed=7),CubMCCLT,abs_tol=1e-2) compare(poly_problem,IIDStdUniform(3,seed=7),CubMCCLT,abs_tol=1e-2) compare(poly_problem,Sobol(3,seed=7),CubQMCSobolG,abs_tol=1e-8)

Stopping Criterion: CubMCCLT        absolute tolerance: 1.0e-02
	W CV:  Solution 5.33       time 1.11       samples 6.7e+06
	WO CV: Solution 5.34       time 0.12       samples 4.8e+05
	Control variates took 10.5% the time and 7.1% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 1.0e-02
	W CV:  Solution 5.33       time 0.98       samples 6.7e+06
	WO CV: Solution 5.34       time 0.09       samples 4.8e+05
	Control variates took 9.1% the time and 7.1% the samples

Stopping Criterion: CubQMCNetG      absolute tolerance: 1.0e-08
	W CV:  Solution 5.33       time 0.12       samples 2.6e+05
	WO CV: Solution 5.33       time 0.04       samples 6.6e+04
	Control variates took 33.4% the time and 25.0% the samples

Problem 2: Keister Function¶

This problem will integrate the Keister function while using control variates $$g_1(x) = \sin(\pi x)$$ and $$g_2(x) = -3(x-1/2)^2+1.$$ The following code does this problem in one-dimension for visualization purposes, but control variates are compatible with any dimension.

In [5]:

def keister_problem(discrete_distrib):
  k = Keister(discrete_distrib)
  cv1 = CustomFun(Uniform(discrete_distrib),lambda x: sin(pi*x).sum(-1))
  cv2 = CustomFun(Uniform(discrete_distrib),lambda x: (-3*(x-.5)**2+1).sum(-1))
  return k,[cv1,cv2],[2/pi,3/4]
compare(keister_problem,IIDStdUniform(1,seed=7),CubMCCLT,abs_tol=5e-4)
compare(keister_problem,IIDStdUniform(1,seed=7),CubMCCLT,abs_tol=4e-4)
compare(keister_problem,Sobol(1,seed=7),CubQMCSobolG,abs_tol=1e-7)

def keister_problem(discrete_distrib): k = Keister(discrete_distrib) cv1 = CustomFun(Uniform(discrete_distrib),lambda x: sin(pi*x).sum(-1)) cv2 = CustomFun(Uniform(discrete_distrib),lambda x: (-3*(x-.5)**2+1).sum(-1)) return k,[cv1,cv2],[2/pi,3/4] compare(keister_problem,IIDStdUniform(1,seed=7),CubMCCLT,abs_tol=5e-4) compare(keister_problem,IIDStdUniform(1,seed=7),CubMCCLT,abs_tol=4e-4) compare(keister_problem,Sobol(1,seed=7),CubQMCSobolG,abs_tol=1e-7)

Stopping Criterion: CubMCCLT        absolute tolerance: 5.0e-04
	W CV:  Solution 1.38       time 1.11       samples 9.5e+06
	WO CV: Solution 1.38       time 0.10       samples 3.4e+05
	Control variates took 9.1% the time and 3.5% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 4.0e-04
	W CV:  Solution 1.38       time 1.71       samples 1.5e+07
	WO CV: Solution 1.38       time 0.17       samples 6.8e+05
	Control variates took 10.2% the time and 4.6% the samples

Stopping Criterion: CubQMCNetG      absolute tolerance: 1.0e-07
	W CV:  Solution 1.38       time 0.31       samples 1.0e+06
	WO CV: Solution 1.38       time 0.41       samples 1.0e+06
	Control variates took 133.0% the time and 100.0% the samples

Problem 3: Option Pricing¶

We will use a European Call Option as a control variate for pricing the Asian Call Option using various stopping criterion, as done for problem 1

In [7]:

call_put = 'call'
start_price = 100
strike_price = 125
volatility = .75
interest_rate = .01 # 1% interest
t_final = 1 # 1 year
def option_problem(discrete_distrib):
  eurocv = FinancialOption(
    discrete_distrib,
    option="EUROPEAN",
    volatility=volatility,
    start_price=start_price,
    strike_price=strike_price,
    interest_rate=interest_rate,
    t_final=t_final,
    call_put=call_put)
  aco = FinancialOption(
    discrete_distrib,
    option="ASIAN",
    volatility=volatility,
    start_price=start_price,
    strike_price=strike_price,
    interest_rate=interest_rate,
    t_final=t_final,
    call_put=call_put)
  mu_eurocv = eurocv.get_exact_value()
  return aco,[eurocv],[mu_eurocv]
compare(option_problem,IIDStdUniform(4,seed=7),CubMCCLT,abs_tol=5e-2)
compare(option_problem,IIDStdUniform(4,seed=7),CubMCCLT,abs_tol=5e-2)
compare(option_problem,Sobol(4,seed=7),CubQMCSobolG,abs_tol=1e-3)

call_put = 'call' start_price = 100 strike_price = 125 volatility = .75 interest_rate = .01 # 1% interest t_final = 1 # 1 year def option_problem(discrete_distrib): eurocv = FinancialOption( discrete_distrib, option="EUROPEAN", volatility=volatility, start_price=start_price, strike_price=strike_price, interest_rate=interest_rate, t_final=t_final, call_put=call_put) aco = FinancialOption( discrete_distrib, option="ASIAN", volatility=volatility, start_price=start_price, strike_price=strike_price, interest_rate=interest_rate, t_final=t_final, call_put=call_put) mu_eurocv = eurocv.get_exact_value() return aco,[eurocv],[mu_eurocv] compare(option_problem,IIDStdUniform(4,seed=7),CubMCCLT,abs_tol=5e-2) compare(option_problem,IIDStdUniform(4,seed=7),CubMCCLT,abs_tol=5e-2) compare(option_problem,Sobol(4,seed=7),CubQMCSobolG,abs_tol=1e-3)

Stopping Criterion: CubMCCLT        absolute tolerance: 5.0e-02
	W CV:  Solution 9.55       time 1.66       samples 3.1e+06
	WO CV: Solution 9.54       time 0.81       samples 7.9e+05
	Control variates took 49.0% the time and 25.6% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 5.0e-02
	W CV:  Solution 9.55       time 1.44       samples 3.1e+06
	WO CV: Solution 9.54       time 0.66       samples 7.9e+05
	Control variates took 45.7% the time and 25.6% the samples

Stopping Criterion: CubQMCNetG      absolute tolerance: 1.0e-03
	W CV:  Solution 9.55       time 0.40       samples 5.2e+05
	WO CV: Solution 9.55       time 0.53       samples 5.2e+05
	Control variates took 132.0% the time and 100.0% the samples

In [ ]: