Importance Sampling Examples

from qmcpy import *
from numpy import *
import pandas as pd
pd.options.display.float_format = '{:.2e}'.format
from matplotlib import pyplot as plt
import matplotlib
%matplotlib inline

Game Example

Consider a game where $X_1, X_2 \overset{\textrm{IID}}{\sim} \mathcal{U}(0,1)$ are drawn with a payoff of

$Y = \text{payoff}(X_1,X_2) = \begin{cases} \$10, & 1.7 \le X_1 + X_2 \le 2, \\ 0, & 0 \le X_1 + X_2 < 1.7, \end{cases}$

What is the expected payoff of this game?

payoff = lambda x: 10*(x.sum(1)>1.7)
abs_tol = 1e-3

Vanilla Monte Carlo

With ordinary Monte Carlo we do the following:

$\mu = \mathbb{E}(Y) = \int_{[0,1]^2} \text{payoff}(x_1,x_2) \, \mathrm{d} x_1 \mathrm{d}x_2$

d = 2
integral = CustomFun(
    true_measure = Uniform(Lattice(d)),
    g = payoff)
solution1,data1 = CubQMCLatticeG(integral, abs_tol).integrate()
data1

LDTransformData (AccumulateData Object)
    solution        0.450
    comb_bound_low  0.449
    comb_bound_high 0.450
    comb_flags      1
    n_total         2^(16)
    n               2^(16)
    time_integrate  0.043
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
CustomFun (Integrand Object)
Uniform (TrueMeasure Object)
    lower_bound     0
    upper_bound     1
Lattice (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       1
    order           natural
    entropy         260768472417330885254602359583380448047
    spawn_key       ()

Monte Carlo with Importance Sampling

We may add the importance sampling to increase the number of samples with positive payoffs. Let

$\boldsymbol{Z} = (X_1^{1/(p+1)}, X_2^{1/(p+1)}), \qquad \boldsymbol{X} \sim \mathcal{U}(0,1)^2$

This means that $Z_1$ and $Z_2$ are IID with common CDF $F(z) =z^{p+1}$ and common PDF $\varrho(z) = (p+1)z^{p}$ . Thus,

$\mu = \mathbb{E}(Y) = \int_{[0,1]^2} \frac{\text{payoff}(z_1,z_2)}{(p+1)^2(z_1z_2)^{p}} \, \varrho(z_1) \varrho(z_2) \, \mathrm{d} z_1 \mathrm{d}z_2 \\ = \int_{[0,1]^2} \frac{\text{payoff}(x_1^{1/(p+1)},x_2^{1/(p+1)})}{(p+1)^2(x_1x_2)^{p/(p+1)}} \, \mathrm{d} x_1 \mathrm{d}x_2$

p = 1
d = 2
integral = CustomFun(
    true_measure = Uniform(Lattice(d)),
    g = lambda x: payoff(x**(1/(p+1))) / ((p+1)**2 * (x.prod(1))**(p/(p+1))))
solution2,data2 = CubQMCLatticeG(integral, abs_tol).integrate()
data2

LDTransformData (AccumulateData Object)
    solution        0.451
    comb_bound_low  0.450
    comb_bound_high 0.451
    comb_flags      1
    n_total         2^(14)
    n               2^(14)
    time_integrate  0.019
CubQMCLatticeG (StoppingCriterion Object)
    abs_tol         0.001
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
CustomFun (Integrand Object)
Uniform (TrueMeasure Object)
    lower_bound     0
    upper_bound     1
Lattice (DiscreteDistribution Object)
    d               2^(1)
    dvec            [0 1]
    randomize       1
    order           natural
    entropy         101794984569737015308869089738144162915
    spawn_key       ()

print('Imporance Sampling takes %.3f the time and %.3f the samples'%\
     (data2.time_integrate/data1.time_integrate,data2.n_total/data1.n_total))

Imporance Sampling takes 0.451 the time and 0.250 the samples

Asian Call Option Example

The stock price must raise significantly for the payoff to be positive. So we will give a upward drift to the Brownian motion that defines the stock price path. We can think of the option price as the multidimensional integral

$\mu = \mathbb{E}[f(\boldsymbol{X})] = \int_{\mathbb{R}^d} f(\boldsymbol{x}) \frac{\exp\bigl(-\frac{1}{2} \boldsymbol{x}^T\mathsf{\Sigma}^{-1} \boldsymbol{x}\bigr)} {\sqrt{(2 \pi)^{d} \det(\mathsf{\Sigma})}} \, \mathrm{d} \boldsymbol{x}$

where

$\boldsymbol{X} \sim \mathcal{N}(\boldsymbol{0}, \mathsf{\Sigma}), \qquad \mathsf{\Sigma} = \bigl(\min(j,k)T/d \bigr)_{j,k=1}^d, \qquad d = 13, \\ f(\boldsymbol{x}) = \max\biggl(K - \frac 1d \sum_{j=1}^d S(jT/d,\boldsymbol{x}), 0 \biggr) \mathrm{e}^{-rT}, \\ S(jT/d,\boldsymbol{x}) = S(0) \exp\bigl((r - \sigma^2/2) jT/d + \sigma x_j\bigr).$

We will replace $\boldsymbol{X}$ by

$\boldsymbol{Z} \sim \mathcal{N}(\boldsymbol{a}, \mathsf{\Sigma}), \qquad \boldsymbol{a} = (aT/d)(1, \ldots, d)$

where a positive $a$ will create more positive payoffs. This corresponds to giving our Brownian motion a drift. To do this we re-write the integral as

$\mu = \mathbb{E}[f_{\mathrm{new}}(\boldsymbol{Z})] = \int_{\mathbb{R}^d} f_{\mathrm{new}}(\boldsymbol{z}) \frac{\exp\bigl(-\frac{1}{2} (\boldsymbol{z}-\boldsymbol{a})^T \mathsf{\Sigma}^{-1} (\boldsymbol{z} - \boldsymbol{a}) \bigr)} {\sqrt{(2 \pi)^{d} \det(\mathsf{\Sigma})}} \, \mathrm{d} \boldsymbol{z} , \\ f_{\mathrm{new}}(\boldsymbol{z}) = f(\boldsymbol{z}) \frac{\exp\bigl(-\frac{1}{2} \boldsymbol{z}^T \mathsf{\Sigma}^{-1} \boldsymbol{z} \bigr)} {\exp\bigl(-\frac{1}{2} (\boldsymbol{z}-\boldsymbol{a})^T \mathsf{\Sigma}^{-1} (\boldsymbol{z} - \boldsymbol{a}) \bigr)} = f(\boldsymbol{z}) \exp\bigl((\boldsymbol{a}/2 - \boldsymbol{z})^T \mathsf{\Sigma}^{-1}\boldsymbol{a} \bigr)$

Finally note that

$\mathsf{\Sigma}^{-1}\boldsymbol{a} = \begin{pmatrix} 0 \\ 0 \\ \vdots \\ 0 \\ a \end{pmatrix}, \qquad f_{\mathrm{new}}(\boldsymbol{z}) = f(\boldsymbol{z}) \exp\bigl((aT/2 - z_d)a \bigr)$

This drift in the Brownian motion may be implemented by changing the drift input to the BrownianMotion object.

abs_tol = 1e-2
dimension = 32

Vanilla Monte Carlo

integrand = AsianOption(Sobol(dimension))
solution1,data1 = CubQMCSobolG(integrand, abs_tol).integrate()
data1

LDTransformData (AccumulateData Object)
    solution        1.788
    comb_bound_low  1.782
    comb_bound_high 1.794
    comb_flags      1
    n_total         2^(12)
    n               2^(12)
    time_integrate  0.014
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dim_frac        0
BrownianMotion (TrueMeasure Object)
    time_vec        [0.031 0.062 0.094 ... 0.938 0.969 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.031 0.031 0.031 ... 0.031 0.031 0.031]
                    [0.031 0.062 0.062 ... 0.062 0.062 0.062]
                    [0.031 0.062 0.094 ... 0.094 0.094 0.094]
                    ...
                    [0.031 0.062 0.094 ... 0.938 0.938 0.938]
                    [0.031 0.062 0.094 ... 0.938 0.969 0.969]
                    [0.031 0.062 0.094 ... 0.938 0.969 1.   ]]
    decomp_type     PCA
Sobol (DiscreteDistribution Object)
    d               2^(5)
    dvec            [ 0  1  2 ... 29 30 31]
    randomize       LMS_DS
    graycode        0
    entropy         296048682239793520955192894339316320391
    spawn_key       ()

Monte Carlo with Importance Sampling

drift = 1
integrand = AsianOption(BrownianMotion(Sobol(dimension),drift=drift))
solution2,data2 = CubQMCSobolG(integrand, abs_tol).integrate()
data2

LDTransformData (AccumulateData Object)
    solution        1.783
    comb_bound_low  1.776
    comb_bound_high 1.790
    comb_flags      1
    n_total         2^(11)
    n               2^(11)
    time_integrate  0.014
CubQMCSobolG (StoppingCriterion Object)
    abs_tol         0.010
    rel_tol         0
    n_init          2^(10)
    n_max           2^(35)
AsianOption (Integrand Object)
    volatility      2^(-1)
    call_put        call
    start_price     30
    strike_price    35
    interest_rate   0
    mean_type       arithmetic
    dim_frac        0
BrownianMotion (TrueMeasure Object)
    time_vec        [0.031 0.062 0.094 ... 0.938 0.969 1.   ]
    drift           0
    mean            [0. 0. 0. ... 0. 0. 0.]
    covariance      [[0.031 0.031 0.031 ... 0.031 0.031 0.031]
                    [0.031 0.062 0.062 ... 0.062 0.062 0.062]
                    [0.031 0.062 0.094 ... 0.094 0.094 0.094]
                    ...
                    [0.031 0.062 0.094 ... 0.938 0.938 0.938]
                    [0.031 0.062 0.094 ... 0.938 0.969 0.969]
                    [0.031 0.062 0.094 ... 0.938 0.969 1.   ]]
    decomp_type     PCA
    transform       BrownianMotion (TrueMeasure Object)
                       time_vec        [0.031 0.062 0.094 ... 0.938 0.969 1.   ]
                       drift           1
                       mean            [0.031 0.062 0.094 ... 0.938 0.969 1.   ]
                       covariance      [[0.031 0.031 0.031 ... 0.031 0.031 0.031]
                                       [0.031 0.062 0.062 ... 0.062 0.062 0.062]
                                       [0.031 0.062 0.094 ... 0.094 0.094 0.094]
                                       ...
                                       [0.031 0.062 0.094 ... 0.938 0.938 0.938]
                                       [0.031 0.062 0.094 ... 0.938 0.969 0.969]
                                       [0.031 0.062 0.094 ... 0.938 0.969 1.   ]]
                       decomp_type     PCA
Sobol (DiscreteDistribution Object)
    d               2^(5)
    dvec            [ 0  1  2 ... 29 30 31]
    randomize       LMS_DS
    graycode        0
    entropy         144971093458566807495711477416245436939
    spawn_key       ()

print('Imporance Sampling takes %.3f the time and %.3f the samples'%\
     (data2.time_integrate/data1.time_integrate,data2.n_total/data1.n_total))

Imporance Sampling takes 1.009 the time and 0.500 the samples

Importance Sampling MC vs QMC

Test Parameters

dimension = 16
abs_tol = .025
trials = 3

df = pd.read_csv('../workouts/mc_vs_qmc/out/importance_sampling.csv')
df['Problem'] = df['Stopping Criterion'] + ' ' + df['Distribution'] + ' (' + df['MC/QMC'] + ')'
df = df.drop(['Stopping Criterion','Distribution','MC/QMC'],axis=1)
problems = ['CubMCCLT IIDStdUniform (MC)',
            'CubMCG IIDStdUniform (MC)',
            'CubQMCCLT Sobol (QMC)',
            'CubQMCLatticeG Lattice (QMC)',
            'CubQMCSobolG Sobol (QMC)']
df = df[df['Problem'].isin(problems)]
mean_shifts = df.mean_shift.unique()
df_samples = df.groupby(['Problem'])['n_samples'].apply(list).reset_index(name='n')
df_times = df.groupby(['Problem'])['time'].apply(list).reset_index(name='time')
df.loc[(df.mean_shift==0) | (df.mean_shift==1)].set_index('Problem')
# Note: mean_shift==0 --> NOT using importance sampling

	mean_shift	solution	n_samples	time
Problem
CubMCCLT IIDStdUniform (MC)	0.00e+00	1.80e+00	3.20e+05	4.90e-01
CubMCCLT IIDStdUniform (MC)	1.00e+00	1.77e+00	8.15e+04	1.30e-01
CubMCG IIDStdUniform (MC)	0.00e+00	1.80e+00	4.18e+05	6.43e-01
CubMCG IIDStdUniform (MC)	1.00e+00	1.77e+00	1.20e+05	1.91e-01
CubQMCCLT Sobol (QMC)	0.00e+00	1.78e+00	4.10e+03	7.69e-03
CubQMCCLT Sobol (QMC)	1.00e+00	1.78e+00	4.10e+03	7.79e-03
CubQMCLatticeG Lattice (QMC)	0.00e+00	1.78e+00	2.05e+03	8.56e-03
CubQMCLatticeG Lattice (QMC)	1.00e+00	1.78e+00	1.02e+03	4.98e-03
CubQMCSobolG Sobol (QMC)	0.00e+00	1.78e+00	1.02e+03	2.06e-03
CubQMCSobolG Sobol (QMC)	1.00e+00	1.78e+00	1.02e+03	1.99e-03

fig,ax = plt.subplots(nrows=1, ncols=2, figsize=(20, 6))
idx = arange(len(problems))
width = .35
ax[0].barh(idx+width,log(df.loc[df.mean_shift==0]['n_samples'].values),width)
ax[0].barh(idx,log(df.loc[df.mean_shift==1]['n_samples'].values),width)
ax[1].barh(idx+width,df.loc[df.mean_shift==0]['time'].values,width)
ax[1].barh(idx,df.loc[df.mean_shift==1]['time'].values,width)
fig.suptitle('Importance Sampling Comparison by Stopping Criterion on Asian Call Option')
xlabs = ['log(Samples)','Time']
for i in range(len(ax)):
    ax[i].set_xlabel(xlabs[i])
    ax[i].spines['top'].set_visible(False)
    ax[i].spines['bottom'].set_visible(False)
    ax[i].spines['right'].set_visible(False)
    ax[i].spines['left'].set_visible(False)
    ax[1].legend(['Vanilla Monte Carlo','Importance Sampling\nMean Shift=1'],loc='upper right',frameon=False)
ax[1].get_yaxis().set_ticks([])
ax[0].set_yticks(idx)
ax[0].set_yticklabels(problems)
plt.tight_layout();

fig,ax = plt.subplots(nrows=1, ncols=2, figsize=(22, 8))
df_samples.apply(lambda row: ax[0].plot(mean_shifts,row.n,label=row['Problem']),axis=1)
df_times.apply(lambda row: ax[1].plot(mean_shifts,row.time,label=row['Problem']),axis=1)
ax[1].legend(frameon=False, loc=(-.85,.7),ncol=1)
ax[0].set_ylabel('Total Samples')
ax[0].set_yscale('log')
ax[1].set_yscale('log')
ax[1].set_ylabel('Runtime')
for i in range(len(ax)):
    ax[i].set_xlabel('mean shift')
    ax[i].spines['top'].set_visible(False)
    ax[i].spines['right'].set_visible(False)
fig.suptitle('Comparing Mean Shift Across Problems');