In mludkov/mlOSP: Machine Learning and Regression Monte Carlo Algorithms for Optimal Stopping

This is a demo file to generate Figure 6 in the paper "Adaptive Batching for Gaussian Process Surrogates with Application in Noisy Level Set Estimation". The plot shows the fitted exercise boundary with its 95\% credible interval (solid line and dashed line) obtained with Gaussian Process for two-dimensional basket put Bermudan option. Two batch heuristics, ABSUR and ADSA, are used to select the location of inputs and their replications, which are shown as the dots and their color/size.

knitr::opts_chunk$set(echo = TRUE)
library(ks)
library(fields) # for plotting purposes
library(mlOSP)
library(DiceKriging)
library(tgp)  # use lhs from there
library(randtoolbox)  # use sobol and halton QMC sequences
library(hetGP)
library(laGP)
library(ggplot2)
library(pander)
data("int_2d")

We consider adaptive batching for a Gaussian GP metamodel for a two-dimensional basket Put Bermudan option. The asset follows log-normal dynamics $$ {Z}{t+\Delta t} = {Z}{t} \cdot\exp \bigg((r-\frac{1}{2} diag{{\Xi}})\Delta t + \sqrt{\Delta t} \cdot {\Xi} (\Delta {W}{t})\bigg) $$ where $\Delta W_t$ are independent Gaussians, and the payoff is $h{Put}(t,{z}) = e^{-r t}( {\cal K} - z^1 - z^2)_+$.

Set up the model for Two-dim basket put with parameters in Table 4

The code below sets up the model for an arithmetic Basket Put with parameters in Table 4 of the article, including the total simulations $N_T$ stored in \textit{total.budget}, initial design size $k_0$ in \textit{init.size}, initial batch size $r_0$ in \textit{batch.nrep} and kernel function $K$ in \textit{kernel.family}. The parameters of the model are initialized as a list, which is later used as input in the main function \textbf{osp.seq.batch.design}.

model2d <- list(x0 = rep(40,2),K=40,sigma=rep(0.2,2),r=0.06,div=0,T=1,dt=0.04,dim=2,sim.func=sim.gbm, payoff.func=put.payoff)
model2d$pilot.nsims <- 1000
model2d$look.ahead <- 1
model2d$cand.len <- 1000   # size of candidate set m_0 for acquisition function
model2d$max.lengthscale <- c(50,50)
model2d$min.lengthscale <- c(5,5)
model2d$tmse.eps <- 0

model2d$ucb.gamma <- 1.96

model2d$seq.design.size <- 100  # budget N = r_0 * k = 2500
model2d$batch.nrep <- 25  # initial replication r_0
model2d$total.budget <- 2500

model2d$init.size <- 20   # initial design size k_0
model2d$init.grid <- int_2d

model2d$tmse.eps <- 0
model2d$kernel.family <- "Matern5_2"  # kernel function for Gaussian Process
model2d$ucb.gamma <- 1.96
model2d$update.freq <- 10   # number of sequential design steps to update the GP surrogate
model2d$r.cand <- c(20, 30,40,50,60, 80, 120, 160) # r_L

Compare results for GP with ABSUR and ADSA

To implement adaptive batching we need to specify the batch heuristic via \textit{batch.heuristic} and the sequential design framework with \textit{ei.func}. The function \textbf{osp.seq.batch.design.simplified} then fits the GP simulator. Fitting is done through the \texttt{DiceKriging} library, namely the main \texttt{km} function there. Below we apply Adaptive Design with Sequential Allocation, ADSA (which relies on the MCU acquisition function) and Adaptive Batching with Stepwise Uncertainty Reduction, ABSUR. The method parameter \texttt{trainkm} means that the GP metamodels are trained using the default \texttt{DiceKriging::km} MLE optimizer.

### GP + ADSA
set.seed(110)
model2d$batch.heuristic <- 'adsa'
model2d$ei.func <- "amcu"
oos.obj.adsa <- osp.seq.batch.design(model2d, method="hetgp")

### GP + ABSUR
set.seed(122)
model2d$batch.heuristic <- 'absur'
model2d$ei.func <- 'absur'
oos.obj.absur <- osp.seq.batch.design(model2d, method="hetgp")

### plot Figure 6
plt.2d.surf.batch(oos.obj.adsa$fit[[15]], oos.obj.absur$fit[[15]], oos.obj.adsa$batches[1:oos.obj.adsa$ndesigns[15] - 1, 15], oos.obj.absur$batches[1:oos.obj.absur$ndesigns[15] - 1, 15], "ADSA", "ABSUR", x=seq(25,50,len=201),y = seq(25, 50,len=201))

The objects \texttt{oos.obj.xxx} are lists that contain the GP metamodels for each time-step of the Bermudan option problem ($M=25$ in the example above). We visualize the fitted timing values at $t=0.6 = 15 \Delta t$. The respective zero contour is the \textit{exercise boundary}.

The function \textbf{plt.2d.surf.with.batch} plots the fitted exercise boundary together with its 95\% credible interval (solid line and dashed curves, respectively). The inputs are shown as the dots and their color indicates the replication amounts $r_i$'s. The first argument is the fitted GP emulator (including the fitted model and the input sites), and the second argument is the batch sizes corresponding to the selected input sites.

First we plot the figure for results obtained with ADSA.

oos.obj.adsa$ndesigns[15]  # number of unique designs
### plot Figure 6 right panel - ADSA
plt.2d.surf.with.batch(oos.obj.adsa$fit[[15]], 
                       oos.obj.adsa$batches[1:oos.obj.adsa$ndesigns[15] - 1, 15])

The second one is for ABSUR

oos.obj.absur$ndesigns[15]  # number of unique designs
### plot Figure 6 left panel - ABSUR
plt.2d.surf.with.batch(oos.obj.absur$fit[[15]], 
                       oos.obj.absur$batches[1:oos.obj.absur$ndesigns[15] - 1, 15])

mludkov/mlOSP documentation built on April 29, 2023, 7:56 p.m.