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 $\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.79       samples 6.7e+06
    WO CV: Solution 5.34       time 0.09       samples 4.2e+05
    Control variates took 11.5% the time and 6.2% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 1.0e-02
    W CV:  Solution 5.33       time 0.68       samples 6.7e+06
    WO CV: Solution 5.34       time 0.07       samples 4.2e+05
    Control variates took 10.7% the time and 6.2% the samples

Stopping Criterion: CubQMCSobolG    absolute tolerance: 1.0e-08
    W CV:  Solution 5.33       time 0.12       samples 2.6e+05
    WO CV: Solution 5.33       time 0.08       samples 1.3e+05
    Control variates took 65.5% 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.21       samples 9.5e+06
    WO CV: Solution 1.38       time 0.14       samples 4.5e+05
    Control variates took 11.8% the time and 4.8% the samples

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

Stopping Criterion: CubQMCSobolG    absolute tolerance: 1.0e-07
    W CV:  Solution 1.38       time 0.44       samples 1.0e+06
    WO CV: Solution 1.38       time 0.47       samples 1.0e+06
    Control variates took 108.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

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.54       time 0.89       samples 2.1e+06
    WO CV: Solution 9.55       time 0.82       samples 1.1e+06
    Control variates took 91.7% the time and 51.9% the samples

Stopping Criterion: CubMCCLT        absolute tolerance: 5.0e-02
    W CV:  Solution 9.54       time 0.84       samples 2.1e+06
    WO CV: Solution 9.55       time 0.82       samples 1.1e+06
    Control variates took 97.0% the time and 51.9% the samples

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