# Control Variates in QMCPy¶

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

## Setup¶

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

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

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

# 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 0.75       samples 7.1e+06
WO CV: Solution 5.33       time 0.09       samples 4.8e+05
Control variates took 11.8% the time and 6.7% the samples

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

Stopping Criterion: CubQMCSobolG    absolute tolerance: 1.0e-08
W CV:  Solution 5.33       time 0.08       samples 2.6e+05
WO CV: Solution 5.33       time 0.05       samples 1.3e+05
Control variates took 64.0% the time and 50.0% the samples


## Problem 2: Keister Function¶

This problem will integrate the Keister function while using control variates

purposes, but control variates are compatible with any dimension.

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.35       samples 9.2e+06
WO CV: Solution 1.38       time 0.13       samples 3.6e+05
Control variates took 9.7% the time and 3.9% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 4.0e-04
W CV:  Solution 1.38       time 2.05       samples 1.4e+07
WO CV: Solution 1.38       time 0.22       samples 7.1e+05
Control variates took 10.8% the time and 4.9% the samples

Stopping Criterion: CubQMCSobolG    absolute tolerance: 1.0e-07
W CV:  Solution 1.38       time 0.39       samples 1.0e+06
WO CV: Solution 1.38       time 0.43       samples 1.0e+06
Control variates took 110.5% 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

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 = EuropeanOption(discrete_distrib,volatility,start_price,strike_price,interest_rate,t_final,call_put)
aco = AsianOption(discrete_distrib,volatility,start_price,strike_price,interest_rate,t_final,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.58       time 2.11       samples 3.4e+06
WO CV: Solution 9.57       time 0.92       samples 8.2e+05
Control variates took 43.8% the time and 24.3% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 5.0e-02
W CV:  Solution 9.58       time 2.11       samples 3.4e+06
WO CV: Solution 9.57       time 0.94       samples 8.2e+05
Control variates took 44.6% the time and 24.3% the samples

Stopping Criterion: CubQMCSobolG    absolute tolerance: 1.0e-03
W CV:  Solution 9.55       time 0.45       samples 5.2e+05
WO CV: Solution 9.55       time 0.58       samples 5.2e+05
Control variates took 127.4% the time and 100.0% the samples