In (2.9) the estimate of ${?}_{t,x}[H_t(X_{\cdot }^x)]$ was obtained pointwise for a given location $(t,x)$ through a plain Monte Carlo approach. Alternatively, given paths $x_{t:T}^n$, $n=1,\mathrm{\dots },N$, emanating from different initial conditions $x_t^n$, one may borrow information cross-sectionally by employing a statistical regression framework. A regression model specifies the link,

where $\epsilon (t,x)$ is the mean-zero noise term with variance ${\sigma }^2(t,x)$. This noise arises from the random component of the pathwise payoff due to the internal fluctuations in $X_{t:T}$ and, hence, $\widehat{\tau }(t,x)$. Letting $y_t^{1:N}\equiv H_t(x_{\cdot }^{1:N})$ be a vector of obtained pathwise payoffs, one regresses $y_t^{1:N}$ against $x_t^{1:N}$ to fit an approximation $\widehat{C}(t,\cdot )$. Evaluating $\widehat{C}(t,\cdot )$ at any $x\in ?$ then yields the predicted continuation value from the exercise strategy ${\widehat{?}}_{t+1:T}$ conditional on $X_t=x$.

The advantage of regression is that it replaces pointwise estimation (which requires rerunning additional scenarios for each location $x$) with a single step that fuses all the information contained in ${(x_t,y_t)}^{1:N}$ to fit $\widehat{C}(t,\cdot )$. Armed with the estimated $\widehat{C}(t,\cdot )$, we take the natural estimator of the stopping region at $t$,

which yields an inductive procedure for approximating the full exercise strategy via (2.7). Note that (2.11) is a double approximation of ${?}_t$. This is due to the discrepancy between $\widehat{C}(t,\cdot )$ and the true expectation ${?}_{t,x}[H_t(X_{\cdot }^x;{\widehat{?}}_{t+1:T})]$ from (2.8). An additional approximation is caused by partially propagating any previous errors in ${\widehat{?}}_{t+1:T}$, relative to the true ${?}_{t+1:T}$.

In the machine learning and simulation optimization literatures, the problem of obtaining $\widehat{C}(t,\cdot )$ based on (2.10) is known as metamodeling or emulation (Kleijnen 2015; Ankenman et al 2010; Santner et al 2013). It consists of two fundamental steps: first, a design (ie, a grid) $?:=x^{1:N}$ is constructed and the corresponding pathwise payoffs $y^{1:N}$ are realized via simulation. Second, a regression model (2.10) is trained to link $y$ with $x$ via an estimated $\widehat{C}(t,\cdot )$. In the context of (2.3), emulation can be traced to the seminal works of Tsitsiklis and Van Roy (2001) and Longstaff and Schwartz (2001) (although the idea of applying regression originated earlier: see, for example, Carrière (1996)). The above references pioneered classical least squares regression (also known by the acronyms LSM (least squares method), LSMC and, more generally, RMC and regression-based schemes) for (2.10), that is, the use of projection onto a finite family of basis functions ${\{B_r(x)\}}_{r=1}^R$,

The projection was carried out by minimizing the $L^2$ distance between $\widehat{C}(t,\cdot )$ and the manifold ${\mathscr{H}}_R:=span(B_r)$ spanned by the basis functions (akin to the definition of conditional expectation):

where the loss function $L_2$ is a weighted $L^2$ norm based on the distribution of $X_t$,

Motivated by the weights in (2.13), the accompanying design $x^{1:N}$ is commonly generated by independent and identically distributed (iid) sampling from $p(t,\cdot \mid 0,X_0)$. As suggested in the original papers of Longstaff and Schwartz (2001) and Tsitsiklis and Van Roy (2001), this can be efficiently implemented by a single global simulation of $X$-paths $x_{0:T}^{1:N}$ at the cost of introducing further $t$-dependencies among $\widehat{C}(t,\cdot )$.

Returning to the more abstract view, least squares regression (2.12) is an instance of a metamodeling technique. Moreover, it obscures the twin issue of generating $?=x^{1:N}$. Indeed, in contrast to standard statistical setups, in metamodeling there is no a priori data per se; instead, the solver is responsible both for generating the data and for training the model. These two tasks are intrinsically intertwined. The subject of design of experiments (DoE) has developed around this issue, but it has been largely absent in RMC approaches. In the next section, we reexamine existing RMC methods and propose novel RMC algorithms that build on the DoE/metamodeling paradigm by targeting efficient designs, $?$.

Figure 1 shows the distribution of $H_t(X_{\cdot }^x)$ in a one-dimensional geometric Brownian motion Bermudan put option problem (see Section 4.4), an archetypal financial application. Figure 1(a) shows a scatter plot of $(x,H_t(X_{\cdot }^x)-h(t,x))$ as $x\in {ℝ}_{+}$ varies, while 1(b) gives a histogram of $H_t(X_{\cdot }^x)$ for a fixed initial value $X_t=x$. Two features become apparent from these plots. First, the noise variance ${\sigma }^2(t,x)$ is extremely large, swamping the actual shape of $x\mapsto T(t,x)$. Second, the noise distribution $\epsilon (t,x)$ is highly nonstandard. It involves a potent brew of (i) heavy tails, (ii) significant skewness and (iii) multimodality. In fact, $\epsilon (t,x)$ does not even have a smooth density, as the nonnegativity of the payoff $h(t,\cdot )\ge 0$ implies that $H_t(X_{\cdot })$ has a point mass at zero.

Moreover, as can be observed in Figure 1(a), the distribution of $\epsilon (t,\cdot )$ is heteroscedastic, with a state-dependent conditional skew and state-dependent point mass at zero. These phenomenons violate the assumptions of ordinary least squares (OLS) regression, making the sampling distribution of the fitted ${\widehat{\beta }}_r$ ill behaved, that is, far from Gaussian. More robust versions of (2.12), including regularized (eg, Lasso) or localized (eg, Loess or piecewise linear models) least squares frameworks, have therefore been proposed (see, for instance, Kohler 2010; Kohler and Krzyżak 2012). Further structured regression proposals may be found in Kan and Reesor (2012) and Létourneau and Stentoft (2014). Alternatively, a range of variance reduction methods, such as control variates, have been suggested to ameliorate the estimation of $\widehat{\beta }$ in (2.12) (see, for example, Agarwal et al 2016; Belomestny et al 2015; Hepperger 2013; Jain and Oosterlee 2015; Juneja and Kalra 2009).

The parametric constraints placed by the LSMC approach (2.12) on the shape of the continuation value $\widehat{C}\in {\mathscr{H}}_R$, make its performance highly sensitive to the distance between the manifold ${\mathscr{H}}_R$ and the true $C(t,\cdot )$ (Stentoft 2004). Moreover, when using (2.12), a delicate balance must be observed between the number of scenarios $N$ and the number of basis functions $R$. These concerns highlight the inflexible nature of parametric regression and become extremely challenging in higher dimensional settings with unknown geometry of $C(t,\cdot )$. One work-around is to employ hierarchical representations, for example, by partitioning $?$ into disjoint bins and employing independent linear hyperplane fits in each bin (see Bouchard and Warin 2011). Nevertheless, ongoing challenges remain in the adaptive construction of ${\mathscr{H}}_R$.

Another limitation of the standard LSMC approach is its objective function. As discussed, least squares regression can be interpreted as minimizing a global weighted mean-squared error loss function. However, as shown in (2.5), the principal aim in optimal stopping is not to learn $C(t,\cdot )$ globally in the sense of its mean-squared error but to rank it vis-à-vis $h(t,\cdot )$. Indeed, recalling the definition of the timing function, the correct loss function $L_{\mathrm{RMC}}$ (see Gramacy and Ludkovski 2015) is

Consequently, the loss criterion is localized around the contour $\partial {?}_t=\{\widehat{C}(t,x)=h(t,x)\}=\{\widehat{T}(t,x)=0\}$. Indeed, a good intuition is that optimal stopping is effectively a classification problem; in some regions of $?$, the sign of $T(t,x)$ is easy to detect and the stopping decision is “obvious”, while in other regions $T(t,x)$ is close to zero and resolving the optimal decision requires more statistical discrimination. For example, in Bermudan options, it is a priori clear that out-of-the-money (OTM) $C(t,x)>h(t,x)$, so it is optimal to continue at such $(t,x)$. This observation was made by Longstaff and Schwartz (2001), who proposed regressing in-the-money (ITM) trajectories only. However, it also belies the apparent inefficiency in the corresponding design; that is, the widely used LSMC implementation first “blindly” generates a design that covers the full domain of $X_t$, only to discard all OTM scenarios. This dramatically shrinks the effective design size, up to 80% in the case of an OTM option. The mismatch between (2.14) and (2.13) makes the terminology of LSMC, with its exclusive emphasis on the $L^2$ criterion, unfortunate.

While one could directly use (2.14) during the cross-sectional regression (to estimate the coefficients $\widehat{\beta }$ in (2.12), for example, though this would require implementing an entirely new optimization problem), a much more effective approach is to take this criterion into account during the experimental design. Intuitively, the local accuracy of any approximation of $C(t,x)$ is limited by the amount of observations in the neighborhood, that is, the density of $?$ around $x$. The form of (2.14), therefore, suggests that the ideal approach would be to place $?$ along the stopping boundary $\partial {?}_t$ in analogy to importance sampling. Of course, the obvious challenge is that knowing ${?}_t$ is equivalent to solving the stopping problem. Nevertheless, this perspective suggests multiple improvements to the baseline approach, which is entirely oblivious to (2.14) when picking $x^{1:N}$.

Motivated by the above discussion, we seek efficient metamodels for obtaining $C(t,\cdot )$ that can exploit the local nature of (2.14). This leads us to consider nonparametric, kernel-based frameworks, specifically Gaussian process regression, also known as kriging (Ankenman et al 2010; Kleijnen 2015; Williams and Rasmussen 2006). We then explore a range of DoE approaches for constructing $?$. Proposed design families include (i) uniform gridding, (ii) random and deterministic space-filling designs, and (iii) adaptive sequential designs based on expected improvement criteria. While only the last choice is truly tailored to (2.14), we explain that reasonably well-chosen versions of (i) and (ii) can already be a significant improvement on the baseline.

To formalize construction of $?$, we rephrase (2.10) as an RSM (or metamodeling) problem. RSM is concerned with the task of estimating an unknown function $x\mapsto f(x)$ that is noisily observed through a stochastic sampler,

where we remind the reader to substitute in their mind $Y(x)\equiv H_t(X_{\cdot }^x)$ and $f\equiv C(t,\cdot )$; $t$ is treated as a parameter. Two basic requirements in RSM are a flexible nonparametric approximation architecture $\mathscr{H}$ that is used to search for the outputted fit $\widehat{f}$, and global consistency: specifically, convergence to the ground truth as the number of simulations $N$ increases without bound. The experimental design problem then aims to maximize the information obtained from $N$ samples ${(x,y)}^{1:N}$ toward the end of minimizing the specified loss function $L(\widehat{f})$.

The pathwise realizations of $H_t(X_{\cdot }^x)$ involve simulation noise that must be smoothed out and, as we saw, is often highly non-Gaussian. To alleviate this issue, we mimic the plain Monte Carlo strategy in (2.9) that uses the classical law of large numbers to obtain a local estimate of $C(t,x)$ by constructing batched or replicated designs. Batching generates multiple independent samples $y^{(i)}\sim Y(x)$, $i=1,\mathrm{\dots },M$, at the same $x$; these samples are then used to compute the empirical average $\overline{y}=(1/M){\sum }_i{y}^{(i)}$. Clearly, $\overline{Y}$ follows the same statistical model (2.15) but with a signal-to-noise ratio improved by a factor of $M$, $?[{\overline{\epsilon }}^2(x)]={\sigma }^2(x)/M$. This also reduces the post-averaged design size $|?|$, which speeds up and improves the training of the metamodel.

To recap, kriging gives a flexible nonparametric representation for $C(t,\cdot )$. The fitting method automatically detects the spatial structure of the continuation value, obviating the need to specify the approximation function class. While the user must pick the kernel family, extensive experiments suggest that these have a second-order effect. Also, the probabilistic structure of kriging offers convenient goodness-of-fit diagnostics after $\widehat{C}$ is obtained, yielding a tractable quantification of posterior uncertainty. This can be exploited for DoE (see Section 5 for details of how).

To fit a kriging metamodel requires specifying ${\mathscr{H}}_{?}$, that is, fixing the kernel family and then estimating the respective kernel hyperparameters (such as $s,\overrightarrow{\theta }$ in (3.4)) that are plugged into (3.5)–(3.6). To do this, one would typically use a maximum likelihood approach; this would entail solving a nonlinear optimization problem to find the maximum likelihood estimators (MLEs) $s^{*},{\overrightarrow{\theta }}^{*}$ defining $?$. Alternatively, cross-validation techniques or expectation–maximization (EM) algorithms are also available.

While the choice of kernel family is akin to the choice of orthonormal basis in LSMC, in our experience, they play only a marginal role in the performance of the overall pricing method. As a rule of thumb, we suggest the use of either the Matérn 5/2 kernel or the squared exponential kernel. In both cases, the respective function spaces are smooth (twice differentiable for the Matérn 5/2 family (3.4) and $C^{\mathrm{\infty }}$ for the squared exponential kernel) and lead to similar fits. The performance is more sensitive to the kernel hyperparameters, which are typically estimated in a frequentist framework but could also be inferred from a Bayesian prior, or even guided by heuristics.

Inference on the kriging hyperparameters requires knowledge of the sampling noise ${\sigma }^2(x)$. Indeed, while it is possible to train $?$ and to infer a constant simulation variance ${\sigma }^2$ simultaneously (the latter is known as the “nugget” in the machine learning community (Santner et al 2013)), with state-dependent noise, $?$ is unidentifiable. Batching allows us to circumvent this challenge by providing empirical estimates of ${\sigma }^2(x)$. For each site $x$, we sample $M$ independent replicates $y^{(1)}(x),\mathrm{\dots },y^{(M)}(x)$ and estimate the conditional variance as

is the batch mean. We then use ${\tilde{\sigma }}^2(x)$ as a proxy for the unknown ${\sigma }^2(x)$ to estimate $?$. One could further improve the estimation of ${\sigma }^2(\cdot )$ by fitting an auxiliary kriging metamodel from ${\tilde{\sigma }}^2$ (see Kleijnen (2015), Chapter 5.6 and the references therein). As shown in Picheny and Ginsbourger (2013), Section 4.4.2, after fixing $?$ we can treat the $M$ samples at $x$ as a single design entry $(x,\overline{y}(x))$ with noise variance ${\tilde{\sigma }}^2(x)/M$. The end result is that the batched data set has just $N^{\prime }:=N/M$ rows ${(x,\overline{y})}^{1:N^{\prime }}$ that are fed into the kriging metamodel.

For our examples below, we have used the R package DiceKriging (Roustant et al 2012). The software takes the input locations $x^{1:N^{\prime }}$, the corresponding simulation outputs $y^{1:N^{\prime }}$ and the noise levels ${\tilde{\sigma }}^2(x^{1:N^{\prime }})$, as well as the kernel family (Matérn 5/2 (3.4) by default), and returns a trained kriging model. One then can utilize a predict call to evaluate (3.5)–(3.6) at given values of $x$. The package can also implement universal kriging. Finally, in cases where the design is not batched or the batch size $M$ is very small (below twenty or so), we resort to assuming a homoscedastic noise level ${\sigma }^2$ that is estimated as part of training the kriging model (an option available in DiceKriging). The advantage of such plug-and-play functionality is the access it offers to an independently tested, speed-optimized, debugged and community-vetted code that minimizes implementation overheads and future maintenance (while perhaps introducing some speed inefficiencies).

To quantify the accuracy of the obtained kriging estimator $\widehat{C}(t,\cdot )$, we derive an empirical estimate of the loss function $L_{\mathrm{RMC}}$. To do so, we integrate the posterior distributions ${ℳ}_x(\cdot )\sim ?(m(x),v^2(x))$ vis-à-vis the surrogate means that are proxies for $C(t,\cdot )$ using (2.14) as follows:

where $\mathrm{\Phi },\varphi$ are the standard Gaussian cumulative/probability density functions. The quantity

gives precisely the Bayesian posterior probability of making the wrong exercise decision at $(t,x)$ based on the information from $?$; thus, a lower ${\mathrm{\ell }}_{\mathrm{RMC}}(x;?)$ indicates a better performance of the design $?$ in terms of stopping optimally. Integrating ${\mathrm{\ell }}_{\mathrm{RMC}}(\cdot )$ over $?$ then gives an estimate for $L_{\mathrm{RMC}}(\widehat{C})$:

Practically, we replace the latter integral with a weighted sum over $?$,

From an approximation theory perspective, kriging is an example of a regularized kernel regression. The general formulation looks for the minimizer $\widehat{f}\in \mathscr{H}$ of the following penalized residual sum of squares (RSS) problem:

where $\lambda \ge 0$ is a chosen smoothing parameter and $\mathscr{H}$ is an RKHS. The summation in (3.10) is a measure of the closeness of $f$ to the data, while the rightmost term penalizes the fluctuations of $f$ to avoid overfitting. The representer theorem implies that the minimizer of (3.10) has an expansion in terms of the eigenfunctions

relating the prediction at $x$ to the kernel functions based at the design sites $x^{1:N}$; compare with (3.5). In kriging, the function space is ${\mathscr{H}}_{?}$, $\lambda =1/2$, and the corresponding norm $\parallel \cdot {\parallel }_{\mathscr{H}}$ has a spectral decomposition in terms of differential operators (Williams and Rasmussen 2006, Chapter 6.2).

Another popular version of (3.10) is based on the smoothing splines (also known as thin-plate splines (TPSs)) that take ${?}_{\mathrm{TPS}}(x,x^{\prime })={\parallel x-x^{\prime }\parallel }^2\log \parallel x-x^{\prime }\parallel$, where $\parallel \cdot \parallel$ denotes the Euclidean norm in ${ℝ}^d$. In this case, the RKHS ${\mathscr{H}}_{\mathrm{TPS}}$ is the set of all twice continuously differentiable functions, as per Hastie et al (2009), Chapter 5, and

This generalizes the one-dimensional penalty function $\parallel f\parallel ={\int }_{ℝ}{\{f^{″}(x)\}}^2\mathrm{d}x$. Note that under $\lambda =\mathrm{\infty }$ and (3.12) the optimization in (3.10) reduces to the traditional least squares linear fit $\widehat{f}(x)={\beta }_0+{\beta }_{1}x$, since it introduces the constraint $f^{″}(x)=0$. A common parameterization for the smoothing parameter $\lambda$ is through the effective degrees of freedom statistic, ${\mathrm{df}}_{\lambda }$; one may also select $\lambda$ adaptively via cross-validation or MLE (Hastie et al 2009, Chapter 5). The resulting optimization of (3.10), based on (3.12), gives a smooth ${?}^2$ response surface, that is, a TPS, and has the explicit form

See Kohler (2008) for an implementation of RMC via splines.

A further RKHS approach, applied to RMC in Tompaidis and Yang (2014), is furnished by the so-called (Gaussian) radial basis functions. These are based on the already mentioned Gaussian kernel ${?}_{\mathrm{RBF}}(x,x^{\prime })=\exp (-\theta {\parallel x-x^{\prime }\parallel }^2)$, but they substitute the penalization in (3.10) with ridge regression, reducing the sum in (3.11) to $N^{\prime }\ll N$ terms by identifying/optimizing the “prototype” sites $x_{n^{\prime }}$.

A standard strategy in experimental design is to spread out the locations $x^{1:N}$ to extract as much information as possible from the samples $y^{1:N}$. A space-filling design attempts to distribute $x^n$ evenly in the input space $\widetilde{?}$, here distinguished from the theoretical domain $?$ of $(X_t)$. While this approach does not directly target the objective (2.14), assuming $\widetilde{?}$ is selected reasonably, it can still be much more efficient than traditional LSMC designs. Examples of space-filling approaches include maximum entropy designs, maximin designs and maximum expected information designs (Kleijnen 2015). Alternatively, one might choose a quasi-Monte Carlo (QMC) method; approaches of this kind use a deterministic space-filling sequence, such as the Sobol, for generating $x^{1:N}$. Gridded designs represent the simplest choice, picking $x^{1:N}$ on a prespecified lattice (see Figure 2). Quasi-random/low-discrepancy sequences have already been used in American option pricing but for different purposes. Boyle et al (2013), for instance, used them for approximating one-step values in a stochastic mesh approach; while Chaudhary (2005) used them for generating the trajectories of $(X_t)$. Here, we propose using quasi-random sequences directly for constructing the simulation design (see also Yang and Tompaidis (2013), who employ a similar strategy in a stochastic control application).

One may also generate random space-filling designs. This proves advantageous in our context of running multiple response surface models (indexed by $t$), since it avoids any grid-ghosting effects. One flexible procedure is Latin hypercube sampling (LHS) (McKay et al 1979). In contrast with deterministic approaches, which are generally only feasible in hyperrectangular $\widetilde{?}$, LHS can easily be combined with an acceptance–rejection step to generate space-filling designs on any finite subspace.

Typically, the theoretical domain $?$ of $(X_t)$ is unbounded (eg, ${ℝ}_{+}^d$), so the success of a space-filling design is driven by the judicious choice of an appropriate $\widetilde{?}$. This issue is akin to specifying the domain for a numerical partial differential equation (PDE) solver and requires structural knowledge. For example, for a Bermudan put with strike $K$, we might take $\widetilde{?}=[0.6K,K]$, covering a wide swath of the ITM region, while eschewing OTM and very deep ITM scenarios. Ultimately, the problem setting determines the importance of this choice (compare the more straightforward setting of Section 6.3 below with that of Section 6.2, where finding a good $\widetilde{?}$ is much harder).

The original solution proposed by Longstaff and Schwartz (2001) was to construct the design $?$ using the conditional density of $X_t\mid X_0$, that is, to generate $N$ independent draws $x^n\sim p(t,\cdot \mid 0,X_0)$, $n=1,\mathrm{\dots },N$. This strategy has two advantages. First, it permits implementing the backward recursion for estimating ${?}_{T-1},{?}_{T-2},\mathrm{\dots }$ on a fixed global set of scenarios that requires just a single (albeit very large) simulation; hence, it saves on the simulation budget. Second, a probabilistic design automatically adapts to the dynamics of $(X_t)$, picking out the regions that $X_t$ is likely to visit and removing the aforementioned challenge of specifying $\widetilde{?}$.

Moreover, the empirical sampling of $x^n$ reflects the density of $X_t$, which helps to lower the weighted loss (2.14). Assuming that the boundary $\partial {?}_t$ is not in a region where $p(t,\cdot \mid 0,X_0)$ is very low, this generates a well-adapted design, which partly explains why the experimental design of LSMC is not entirely without its merits. Compared with a space-filling design that is sensitive to $\widetilde{?}$, a probabilistic design is primarily driven by the initial condition $X_0$. This can be a boon (as discussed) or a limitation; with OTM put contracts, for example, the vast majority of empirically sampled sites $x^n$ would be OTM, dramatically lowering the quality of an empirical $?$. A good compromise is to take $x^n\sim X_t\mid \{X_0,h(t,X_t)>0\}$, that is, to effectively condition on being ITM.

As discussed, the design sites $x^{1:N}$ need not be distinct, and the macro design can be replicated. Such batched designs offer several benefits in the context of metamodeling. First, the averaging of the simulation outputs dramatically improves the statistical properties of the averaged $\overline{y}$ compared with the raw $y$. In most cases, once $M\gg 10$, the Gaussian assumption regarding the respective noise $\overline{\epsilon }$ becomes excellent. Second, batching raises the signal-to-noise ratio and hence simplifies response surface modeling. Training the kriging kernel $?$ with very noisy samples is difficult, as the likelihood function to be optimized possesses multiple local maximums. Third, as mentioned in Section 3.1, batching allows the empirical estimation of heteroscedastic sampling noise ${\sigma }^2(x)$. Algorithm 1 summarizes the steps for building a kriging-based RMC solver with a replicated design.

Batching also connects with the general bias–variance trade-off in statistical estimation. Plain (pointwise) Monte Carlo offers a zero bias/high-variance estimator of $f(x)$, which is therefore computationally expensive. A low-complexity model, such as parametric least squares (2.12), offers a low-variance/high-bias estimator (since it requires constraining $\widehat{f}\in {\mathscr{H}}_R$). Coupled with kriging, batched designs interpolate these two extremes by reducing bias through flexible approximation classes and reducing variance through the batching and modeling of spatial response correlation.

ALGORITHM 1  Regression Monte Carlo for optimal stopping using kriging.

As the number of replications $M$ becomes very large, batching effectively eliminates all sampling noise. Consequently, one can substitute the empirical average $\overline{y}(x)$ from (3.7) for the true $f(x)$, an approach first introduced by Van Beers and Kleijnen (2003). It now remains only to interpolate these responses, reducing the metamodeling problem to the analysis of deterministic experiments. There are a number of highly efficient interpolating metamodels, including classical deterministic kriging (which takes $\mathrm{\Sigma }\equiv \mathrm{?}$ in (3.5)–(3.6)) and natural cubic splines. Assuming a smooth underlying response, deterministic metamodels often enjoy very fast convergence. Application of interpolators offers a new perspective (which is also more pleasing visually) on the RMC metamodeling problem. To sum up, batching offers a tunable spectrum of methods that blend the smoothing and interpolation of $f(\cdot )$.

Figure 2 illustrates the above idea for a one-dimensional Bermudan put in a GBM model (see Section 4.4). We construct a design $?$ of macro size $9600=1600\times 6$, which uses just six distinct sites ${?}^{\prime }=x^{1:N^{\prime }}$, with $M=1600$ replications at each site. We then interpolate the obtained values $\overline{y}(x^{1:6})$ using (i) a deterministic kriging model and (ii) a natural cubic spline. We observe that the resulting fit for $\widehat{T}(t,\cdot )$, shown in Figure 2, looks excellent and yields an accurate approximation of the exercise boundary.

To illustrate different simulation designs, we revisit the classic example of a one-dimensional Bermudan put option in a geometric Brownian motion model. More extensive experiments are given in Section 6. The lognormal $X$-dynamics are

and the put payoff is $h(t,x)={\mathrm{e}}^{-rt}{(K-x)}_{+}$. The option matures at $T=1$, and we assume twenty-five exercise opportunities spaced out evenly with $\mathrm{\Delta }t=0.04$. This means that, with a slight abuse of notation, the RMC algorithm is executed over $t=T,T-\mathrm{\Delta }t,\mathrm{\dots },\mathrm{\Delta }t,0$. The rest of the parameter values are interest rate $r=0.06$, dividend yield $\delta =0$, volatility $\sigma =0.2$ and at-the-money (ATM) strike $K=X_0=40$. In this case, $V(0,X_0)=2.314$ (see Longstaff and Schwartz 2001). In one dimension, the stopping region for the Bermudan put is an interval $[0,\underset{¯}{s}(t)]$; thus, there is a unique exercise boundary $\underset{¯}{s}(t)=\partial {?}_t$.

We implement Algorithm 1 using a batched LHS design $?$ with $N=3000$, $M=100$, so that at each step there are $N^{\prime }=N/M=30$ distinct simulation sites $x^{1:N^{\prime }}$. The domain for the experimental designs is the ITM region $\widetilde{?}=[25,40]$. To focus on the effect of the various instances of $?$, we freeze the kriging kernel $?$ across all time steps $t$ as a Matérn 5/2-type (3.4) with hyperparameters $s=1,\theta =4$, which were close to the MLE estimates obtained. This choice is maintained for the rest of this section as well as for Figures 4 and 5. Our experiments (see, in particular, Table 3) imply that the algorithm is insensitive to both the kernel family and the specific software/approach used for training the kriging kernels.

Figure 3 illustrates the effect of various experimental designs on the kriging metamodels. We compare three different batched designs with a fixed overall size $N=|?|=3000$. These designs are (1) an LHS design with small batches $M=20$, $N^{\prime }=150$, as illustrated in panels (a), (d) and (g); (2) an LHS design with a large $M=100$, $N^{\prime }=30$, as illustrated in panels (b), (e) and (h); and (3) a probabilistic density-based design also with $M=100$, $N^{\prime }=30$, as illustrated in panels (c), (f) and (i). Recall that the latter generates $N/M=30$ paths of $(X_t^x)$ and then creates $M$ subsimulations initiating at each $x_t^k$. The resulting ${(x_t,y_t)}^{1:N^{\prime }}$ at $t=0.6$ are shown in panels (a)–(c) of Figure 3, along with the obtained estimates $\widehat{T}(t,\cdot )$. Panels (d)–(f) of Figure 3 show the resulting surrogate standard deviations $v(\cdot )$. The shape of $x\mapsto v(x)$ is driven by the local density of the respective $?$ as well as the simulation noise ${\sigma }^2(x)$. Here, ${\sigma }^2(x)$ is highly heteroscedastic, being much larger near ATM $x\simeq 40$ than deep ITM. This is because ITM one is likely to stop sooner, and so the variance of $H_t(X_{\cdot }^x)$ is lower. For LHS designs, the roughly uniform density of $?$ leads to the posterior variance being proportional to the simulation noise, $v^2(x)\propto {\sigma }^2(x)$; in contrast, the empirical design reflects the higher density of $X_t$ closer to $X_0=40$, so that $v(x)$ is smallest there. Moreover, because there are very few design sites for (just five in Figure 3), the corresponding surrogate variance has a distinctive hill-and-valley shape. In all cases, the surrogate variance explodes along the edges of $?$, reflecting the lack of information about the response there.

Panels (g)–(i) of Figure 3 show the local loss metric ${\mathrm{\ell }}_{\mathrm{RMC}}(x;?)$ for the corresponding designs. We stress that ${\mathrm{\ell }}_{\mathrm{RMC}}(\cdot )$ is a function of the sampled $y^{1:N}$ and hence will vary across algorithm runs. Intuitively, ${\mathrm{\ell }}_{\mathrm{RMC}}(\cdot )$ is driven by the respective surrogate variance $v^2(x)$, weighted in terms of the distance $|m(x)-h(t,x)|$ to the stopping boundary. Consequently, large $v^2(x)$ is fine for regions far from the exercise boundary. Overall, we can make the following observations.

1. (i)

   Replicating simulations does not materially impact on surrogate variance at any given $x^{\prime }\notin ?$ and hence has little effect on ${\mathrm{\ell }}_{\mathrm{RMC}}(x^{\prime })$. Consequently, there is minimal loss of fidelity relative to using a lot of distinct design sites in $?$.

2. (ii)

   The modeled response surfaces tend to be smooth. Consequently, the kernel $?$ should enforce long-range spatial dependence (via a sufficiently large length scale) to ensure that $\widehat{T}(t,\cdot )$ is likewise smooth. Undesirable fluctuations in $\widehat{T}(t,\cdot )$ can generate spurious stopping/continuation regions (see, for example, panel (a) in Figure 3 around $x=28$).

3. (iii)

   Probabilistic designs tend to increase ${\widehat{L}}_{\mathrm{RMC}}(\widehat{C})$ because they fail to place enough design sites deep ITM, leading to extreme posterior uncertainty $v^2(x)$ there (see panel (i) in Figure 3).

4. (iv)

   Due to the very low signal-to-noise ratio that is common in financial applications, replication of the simulation aids in seeing the “shape” of $C(t,\cdot )$ and hence simplifies the task of metamodeling (see point (ii)).

To give a sense of the aggregate quality of the metamodels, Figure 4 shows the estimated $\widehat{T}(t,\cdot )$ at three different steps $t$ as the RMC implements backward recursion. The overall shape of $\widehat{T}(t,\cdot )$ remains essentially the same, being strongly positive close to the strike $K=40$, crossing zero at $\underset{¯}{s}$ and remaining slightly negative deep ITM. As $t$ decreases, the exercise boundary ${\underset{¯}{s}}_t$ also moves lower. Figure 4 shows how the cross-sectional borrowing of information (modeled by the Gaussian process correlation kernel $?$) reduces the surrogate variance $v^2(x)$ compared with the raw ${\tilde{\sigma }}^2(x)$ (see the corresponding 95% credibility intervals (CIs)). We also observe some undesirable “wiggling” of $\widehat{T}(t,\cdot )$, which can generate spurious exercise boundaries such as in the extreme left of Figure 4(c). However, with an LHS design, this effect is largely mitigated compared with the performance of the empirical design in Figure 3. This confirms the importance of spacing out the design $?$. Indeed, for the LHS designs, the phenomenon of spurious exercising only arises for $t\simeq 0$ and very small $x$, whereby the event that would lead to the wrong exercise strategy has a negligible probability of occurring.

Fixing a time step $t$, the sequential design loop is initialized with ${?}^{(N_0)}$, generated, for example, using LHS. New design sites $x^{N_0+1},x^{N_0+2},\mathrm{\dots }$ are then iteratively added by greedily optimizing an acquisition function $E(x)$ that quantifies information gains via expected improvement (EI) scores. The aim of EI scores is to identify the locations $x$ that are most promising in terms of lowering the global loss function $L_{\mathrm{RMC}}(\widehat{C})$ from (2.14). In our context, the EI scores $E_k(x)$ are based on the posterior distributions ${ℳ}_x^{(k)}(\cdot )$, which summarize information learned so far about $C(t,\cdot )$. The seminal works of Cohn (1996) and MacKay (1992) showed that a good heuristic for maximizing learning rates is to sample at sites that have high surrogate variance or high expected surrogate variance reduction, respectively. Because we target (2.14), a second objective is to preferentially explore regions close to the contour $\{\widehat{C}(t,x)=h(t,x)\}$. This is achieved by blending the distance to the contour with the above variance metrics, in analogy to the efficient global optimization approach (Jones et al 1998) in simulation optimization.

A greedy strategy augments ${?}^{(k)}$ with the location that maximizes the acquisition function

Two major concerns for implementing (5.1) are (i) the computational cost of optimizing $x^{k+1}$ over $?$ and (ii) the danger of myopic strategies. For the first aspect, we note that (5.1) introduces a whole new optimization subproblem just for augmenting the design. This can generate substantial overhead that deteriorates the running time of the entire RMC. Consequently, approximate optimality is a possible pragmatic compromise to maximize performance. For the second aspect, the myopic nature of (5.1) might lead to an undesirable concentration of the design $?$ interfering with the convergence of ${\widehat{C}}^{(k)}(t,\cdot )\to C(t,\cdot )$ as $k$ grows. It is well known that many greedy sequential schemes can get trapped sampling in some subregion of $?$, generating poor estimates elsewhere. For example, if there are multiple exercise boundaries, the EI metric might over-prefer established boundaries at the risk of missing other stopping regions that have not yet been found. A common way to circumvent this concern is to randomize (5.1), thus ensuring that $?$ grows dense uniformly on $\widetilde{?}$. In the examples below we follow Gramacy and Ludkovski (2015) and replace $\arg {sup}_{x\in ?}$ in (5.1) with $\underset{x\in ?}{\arg \max }$, where $?$ is a finite candidate set, generated using LHS over some domain $\widetilde{?}$. LHS candidates ensure that potential new $x^{k+1}$ locations are representative and well spaced out over $?$.

Kriging metamodels are especially well suited to sequential design thanks to the availability of simple updating formulas that allow us to efficiently assimilate new data points into an existing fit. If a new sample $(x^{k+1},y^{k+1})$ is added to an existing design $x^{1:k}\equiv {?}^{(k)}$ while keeping the kernel $?$ fixed, the surrogate mean and variance at location $x$ are updated via