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

g(t) = 10t_1-5t_2^2+2t_3^3

with true measure :math:`\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.

# 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

g_1(x) = \sin(\pi x)

and

g_2(x) = -3(x-1/2)^2+1.

The following code does this problem in 1 dimension for visualization

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