Sturt

A nonparametric algorithm for optimal stopping based on robust optimization

A nonparametric algorithm for optimal stopping based on robust optimization

Bradley Sturt Department of Information and Decision Sciences
University of Illinois at Chicago<EMAIL_ADDRESS>
Optimal stopping is a fundamental class of stochastic dynamic optimization
problems with numerous applications in finance and operations management. We
introduce a new approach for solving computationally-demanding stochastic
optimal stopping problems with known probability distributions. The approach
uses simulation to construct a robust optimization problem that approximates
the stochastic optimal stopping problem to any arbitrary accuracy; we then
solve the robust optimization problem to obtain near-optimal Markovian
stopping rules for the stochastic optimal stopping problem.

In this paper, we focus on designing algorithms for solving the robust
optimization problems that approximate the stochastic optimal stopping
problems. These robust optimization problems are challenging to solve because
they require optimizing over the infinite-dimensional space of all Markovian
stopping rules. We overcome this challenge by characterizing the structure of
optimal Markovian stopping rules for the robust optimization problems. In
particular, we show that optimal Markovian stopping rules for the robust
optimization problems have a structure that is surprisingly simple and finite-
dimensional. We leverage this structure to develop an exact reformulation of
the robust optimization problem as a zero-one bilinear program over totally
unimodular constraints. We show that the bilinear program can be solved in
polynomial time in special cases, establish computational complexity results
for general cases, and develop polynomial-time heuristics by relating the
bilinear program to the maximal closure problem from graph theory. Numerical
experiments demonstrate that our algorithms for solving the robust
optimization problems are practical and can outperform state-of-the-art
simulation-based algorithms in the context of widely-studied stochastic
optimal stopping problems from high-dimensional option pricing.

robust optimization; optimal stopping; options pricing.

First version: March 4, 2021. Revisions submitted on June 17, 2022 and
December 21, 2022.

##  1 Introduction

Consider the following class of stochastic dynamic optimization problems: A
sequence of random states are incrementally revealed to a decision maker.
After observing the state in each period, the decision maker chooses whether
to continue to the next period or stop and receive a reward that depends on
the current state. The problem is to find a control policy, called a stopping
rule, for selecting when to stop the process to maximize the expected reward.

Such optimal stopping problems are widely studied and arise in a variety of
domains like finance, promotion planning (Feng and Gallego 1995), and organ
transplantation (David and Yechiali 1985). In particular, optimal stopping has
considerable importance to industry for the pricing of financial derivatives.
With a record trading volume that exceeded seven billion contracts in 2020,
equity options are among the most widely-traded type of financial derivative
(Reuters 2021), and financial firms depend on solving optimal stopping
problems to determine accurate prices for American-style options, the most
common type of equity option.

In this paper, we study a general class of optimal stopping problems in which
the sequence of random states is driven by a non-Markovian probability
distribution. We recall that a sequence of random states is non-Markovian if
the state in the next time period (e.g., a stock’s price tomorrow) has a
probability distribution which depends both on the state in the current time
period (e.g., the stock’s price today) as well as the states in the past
periods (e.g., the stock’s price yesterday). This class of optimal stopping
problems has witnessed a surge of interest as financial firms increasingly use
non-Markovian probability distributions to accurately model the volatility
patterns of stocks (Gatheral et al. 2018, Leão et al. 2019, Becker et al.
2019, Bezerra et al. 2020, Goudenège et al. 2020, Bayer et al. 2020). Optimal
stopping problems with non-Markovian probability distributions also occur when
using popular dimensionality-reduction techniques for pricing high-dimensional
basket options (Bayer et al. 2019, p. 372) and pricing options when the
probability distribution of underlying assets is accessed via a black-box
simulator constructed from historical data (Ciocan and Mišić 2020, §5.5).

Despite their importance in practice, non-Markovian optimal stopping (NMOS)
problems are “not easy to solve” (Leão et al. 2019, p. 982). The difficulty of
these problems arises because the optimal decision in each period may depend
on the entire history of the state process. In principle, an NMOS problem can
be transformed into an equivalent Markovian optimal stopping problem by
converting the original sequence of random states
$x_{1},\ldots,x_{T}\in\mathcal{X}$ into a Markovian stochastic process
$X_{1},\ldots,X_{T}\in\mathcal{X}^{T},$ where the new state
$X_{t}\coloneqq(x_{1},\ldots,x_{t},0,\ldots,0)\in\mathcal{X}^{T}$ in each
period $t$ includes the entire state history of the original process.
Unfortunately, the enlarged state space $\mathcal{X}^{T}$ will be high-
dimensional when the optimal stopping problem has many periods, and the
difficulty of solving a Markovian optimal stopping problem explodes in the
dimensionality of the state space.

To contend with the curse-of-dimensionality that arises in NMOS problems, a
natural approximation technique is to search only for stopping rules that are
Markovian. Rather than depending on the entire history of the original
sequence, a Markovian stopping rule makes a decision in each period $t$ based
only on the current original state $x_{t}$ and knowledge that the sequence of
random states was not stopped in any of the previous time periods. In general,
the best Markovian stopping rule for an NMOS problem is not guaranteed to be
an optimal stopping rule for the NMOS problem. However, recent numerical
evidence demonstrates that Markovian stopping rules can lead to highly
accurate approximations of optimal stopping rules in NMOS problems from
options pricing; see Goudenège et al. (2020, §9.2), Ciocan and Mišić (2020,
§5.5), Bayer et al. (2020, §3.5).

As far as we are aware, only two papers until now have suggested methods that
are theoretically capable of finding the best Markovian stopping rules to NMOS
problems. Belomestny (2011a) analyzes simulation-based methods which optimize
directly over spaces of Markovian stopping rules and suggests a nonparametric
space of Markovian stopping rules based on $\epsilon$-nets; however, he does
not propose any concrete algorithms for optimizing over this nonparametric
space of Markovian stopping rules. Ciocan and Mišić (2020) propose optimizing
over Markovian stopping rules which are restricted to decision trees with
fixed depth. Due to the computational intractability of optimizing over all
decision trees of fixed depth, the authors develop greedy heuristics which are
shown to limit the range of attainable decision trees, and overcoming this
limitation of the heuristics “is not obvious, especially in light of the
structure of the optimal stopping problem that is leveraged to efficiently
optimize split points in our construction” (Ciocan and Mišić 2020, p. 22).

We take a different approach to the aforementioned literature, and in doing so
make our contributions to optimal stopping, by drawing on the traditionally
unrelated field of _robust optimization_. Over the past two decades, robust
optimization has emerged as a leading tool in operations research for dynamic
decision-making when uncertainty is driven by unknown or ambiguous probability
distributions (Ben-Tal et al. 2009, Delage and Iancu 2015). In this paper, we
show that robust optimization can be combined with simulation to develop
algorithms for finding Markovian stopping rules to NMOS problems with known
probability distributions. Compared to Belomestny (2011a) and Ciocan and Mišić
(2020), our approach for optimal stopping does not restrict the space of
Markovian stopping rules to any parametric class, and we develop concrete
algorithms that are guaranteed to yield $\epsilon$-optimal Markovian stopping
rules for general classes of NMOS problems.

In greater detail, this paper introduces a new approach for computing
Markovian stopping rules for NMOS problems with known probability
distributions that is based on a combination of simulation and robust
optimization. At a high level, our approach is comprised of the following
steps:

Step 1:

We use Monte-Carlo simulation to generate sample paths of the sequence of
random states.

Step 2:

From those sample paths, we construct a robust optimization problem that
approximates the NMOS problem.

Step 3:

We solve the robust optimization problem to obtain stopping rules for the NMOS
problem.

The robust optimization problem constructed in Step 2 can be interpreted as a
_proxy_ or _surrogate_ for the NMOS problem. Indeed, if the number of
simulated sample paths in Step 1 is sufficiently large, then every Markovian
stopping rule that is optimal for the robust optimization problem constructed
in Step 2 is guaranteed with high probability to be an $\epsilon$-optimal
Markovian stopping rule for the NMOS problem; see §2.4. Thus, the approach
comprised of the above steps enables the task of computing near-optimal
Markovian stopping rules for an NMOS problem to be reduced to the task of
solving a robust optimization problem.

With respect to the robust optimization literature, Step 2 of our approach
follows a recent paper by Bertsimas, Shtern, and Sturt (2022, henceforth
referred to as BSS22). In that paper, the authors showed that a general class
of stochastic dynamic optimization problems with unknown probability
distributions can in principle be approximated to arbitrary accuracy by a
robust optimization problem constructed from historical data. In particular,
our work draws on BSS22 in two specific ways. First, our robust optimization
problem constructed in Step 2 is a variant of a robust optimization
formulation that is proposed in BSS22. Second, by focusing on the application
of optimal stopping, we strengthen theoretical developments from BSS22 to
prove under relatively mild and verifiable assumptions that the optimal
objective value and optimal Markovian stopping rules of the robust
optimization problem from Step 2 converge almost surely to those of the NMOS
problem as the number of simulated sample paths in Step 1 grows to infinity
(Theorems 2.2-2.4 in §2.4). For the interested reader, a discussion of our
improvements to the theoretical convergence guarantees from BSS22 is provided
at the beginning of Appendix 9.

In contrast to BSS22, the key novelty of the present paper lies _not_ in
showing that robust optimization can be used to construct approximations of
stochastic dynamic optimization problems; rather, the main contributions of
the present paper are _algorithmic_. Indeed, in order for a combination of
simulation and robust optimization to yield near-optimal algorithms for a
class of stochastic dynamic optimization problems with known probability
distributions, it is not sufficient to construct a robust optimization problem
that approximates the stochastic problem to arbitrary accuracy: one must also
have exact or provably near-optimal algorithms for solving the robust
optimization problem. However, up to this point, there have been no exact
algorithms in the literature for solving the type of robust optimization
problems proposed by BSS22 in any class of stochastic dynamic optimization
problems where uncertainty unfolds over two or more periods. In particular,
because solving the type of robust optimization problems from BSS22 requires
optimizing over infinite-dimensional spaces of control policies, it has been
unknown in any application whether optimal control policies for these robust
optimization problems even _exist_ , let alone whether they can ever be
tractably computed. These hurdles have, practically speaking, prevented the
robust optimization techniques from BSS22 from being combined with simulation
to develop algorithms for computing near-optimal control policies for any
class of stochastic dynamic optimization problems with known probability
distributions until now.

In this paper, we resolve the aforementioned gaps in the robust optimization
literature by establishing the first characterization of optimal control
policies for the type of robust optimization problems proposed by BSS22 in an
application where uncertainty unfolds over two or more periods (Theorem 3.1 in
§3). Specifically, we consider the task of solving the robust optimization
problems constructed in Step 2 over the infinite-dimensional space of all
Markovian stopping rules. Using a novel pruning technique, we prove in Theorem
3.1 that optimal Markovian stopping rules for these robust optimization
problems not only exist, but also have a structure that is simple and finite-
dimensional. In fact, our characterization reveals that optimal Markovian
stopping rules for these robust optimization problems can be compactly
parameterized by integer variables
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$, where $N$ is the number of
simulated sample paths chosen in Step 1 and $T$ is the number of periods in
the NMOS problem.

Leveraging the structure of optimal Markovian stopping rules for the robust
optimization problems constructed in Step 2, we develop exact and heuristic
algorithms for solving the robust optimization problem in Step 3.
Specifically, we make the following algorithmic contributions:

  1. (a)

We use our characterization of optimal Markovian stopping rules for the robust
optimization problems constructed in Step 2 to transform the robust
optimization problem into a finite-dimensional optimization problem over
integers $\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$ (Theorem 3.7 in
§3.3).

  2. (b)

For NMOS problems with two periods, we show that optimal Markovian stopping
rules for the robust optimization problem can be computed in
$\mathcal{O}(N^{3})$ time (Theorem 4.1 in §4.1).

  3. (c)

For NMOS problems with three or more periods, we prove that finding optimal
Markovian stopping rules for the robust optimization problem is NP-hard
(Theorem 4.2 in §4.1).

  4. (d)

We develop a general-purpose exact algorithm for computing optimal Markovian
stopping rules for the robust optimization problem that consists of solving a
zero-one bilinear program over totally unimodular constraints (Theorem 4.3 in
§4.1).

  5. (e)

We design and analyze a polynomial-time heuristic algorithm for approximately
solving the robust optimization problem by relating the bilinear program to
the maximal closure problem from graph theory (Propositions 4.6-4.13 in §4.2).

In summary, our exact and heuristic algorithms allow us to solve the robust
optimization problem constructed in Step 2. Since the robust optimization
problem constructed via Steps 1 and 2 serves as an approximation of the NMOS
problem to any arbitrary accuracy, our algorithms for solving the robust
optimization problem in turn allow us to compute near-optimal Markovian
stopping rules for the NMOS problem.

We conclude with numerical experiments that demonstrate the value of our
robust optimization-based algorithms in several settings. First, we consider a
simple one-dimensional non-Markovian optimal stopping problem with fifty
periods, and we compare the robust optimization algorithm to existing methods
based on approximate dynamic programming (Longstaff and Schwartz 2001) and
parametric stopping rules (Ciocan and Mišić 2020). The experiments show that
our method can find stopping rules that significantly outperform those found
by the other techniques, while maintaining a comparable computational cost. In
particular, the experiments reveal that our method can strictly outperform
alternative algorithms for finding Markovian stopping rules to NMOS problems
that are based on backwards recursion. Second, we consider a widely-studied
and important problem of pricing high-dimensional Bermudan barrier options
with over fifty periods. Across several variants of this problem, we
demonstrate that our combination of robust optimization and simulation can
find stopping rules that match, and in some cases significantly outperform,
those from state-of-the-art algorithms by Longstaff and Schwartz (2001),
Ciocan and Mišić (2020) as well as the duality-based pathwise optimization
method of Desai et al. (2012).

The rest of our paper has the following organization. §1.1 provides a review
of methods for solving optimal stopping problems with Markovian and non-
Markovian probability distributions. §2 formalizes the problem setting and
introduces our robust optimization-based method. §3 characterizes the
structure of optimal policies for the robust optimization problem. §4 develops
tractable algorithms and computational complexity results for the robust
optimization problem. §5 illustrates the performance of our algorithms in
numerical experiments. Unless stated otherwise, all technical proofs can be
found in the appendices.

###  1.1 Other Related Literature

Many methods based on approximate dynamic programming (ADP) have been
developed for optimal stopping problems with high-dimensional Markovian
stochastic processes. The most popular ADP methods for these optimal stopping
problems are based on Monte-Carlo simulation and regression, which originate
with Carriere (1996), Longstaff and Schwartz (2001) and Tsitsiklis and Van Roy
(2001). Given sample paths of the entire stochastic process, these methods use
backwards recursion and regression to obtain approximations of the value
function, and exercise policies are then obtained by proceeding greedily with
respect to the approximate value functions. The efficacy of regression-based
methods hinges on selecting a parametrization of basis functions for the value
function that strikes a balance between approximation quality and sample
complexity. Nonparametric choices for the basis functions, e.g., Lagurerre
polynomials, are discussed in the aforementioned works and subsequently
analyzed in works such as Clément et al. (2002), Glasserman and Yu (2004),
Egloff (2005), Belomestny (2011b) and Zanger (2020). In §5, we provide
numerical comparisons of our proposed algorithms to ADP techniques in the
context of NMOS problems.

A variety of other nonparametric methods have been developed for solving
optimal stopping problems with Markovian probability distributions, such as
quantization-based approximations of value functions (Bally and Pages 2003)
and scenario tree discretizations of the sequence of random states (Broadie
and Glasserman 1997). Recent works have also considered using deep learning to
learn the continuation function, including Becker et al. (2019) and Goudenège
et al. (2020). The efficacy of deep learning algorithms for optimal stopping
hinges on carefully selecting the topology of the neural network and choosing
the right tuning algorithm, and performing these tasks effectively in the
context of optimal stopping is an ongoing area of research (Fathan and Delage
2021).  Methods to compute upper bounds on optimal stopping problems grew in
interest due to the independent works of Haugh and Kogan (2004) and Rogers
(2002), and duality-based algorithms to obtain upper bounds which combine
simulation and suboptimal stopping rules were first proposed by Andersen and
Broadie (2004). Other works that harness dual representations to solve optimal
stopping and other stochastic dynamic optimization problems include Brown et
al. (2010), Desai et al. (2012), Belomestny (2013), and Goldberg and Chen
(2018), among many others.  In §5, we provide numerical comparisons of our
proposed algorithms to the duality-based pathwise optimization method of Desai
et al. (2012).

In the context of non-Markovian optimal stopping, methods have been developed
which address settings that are different from ours in non-trivial ways. Leão
et al. (2019) and Bezerra et al. (2020) develop discretization schemes for
NMOS problems over continuous time and restrict the class of probability
distributions to those based on the Brownian motion. In contrast to these
works, our paper develops algorithms for finding Markovian stopping rules for
discrete-time optimal stopping problems, and we do not require any parametric
assumptions on the probability distributions of the underlying stochastic
processes. NMOS problems can also be addressed by the scenario tree method of
Broadie and Glasserman (1997) and the recursive-dual algorithm of Goldberg and
Chen (2018), provided that one can perform Monte-Carlo simulation on the
conditional probability distribution of the stochastic process in each time
period. In contrast to these methods, the algorithms in this paper require
only the ability to simulate sample paths of the entire stochastic process and
are shown in numerical experiments to be practically tractable in low-
dimensional NMOS problems with dozens of time periods. Within the optimal
stopping literature, our method is most closely related to a class of
simulation-based methods which optimize directly over spaces of deterministic
stopping rules, as explored by Garcıa (2003), Andersen (1999), Belomestny
(2011a), Gemmrich (2012), Ciocan and Mišić (2020), and Glasserman (2013,
§8.2). We discuss connections between our approach and this stream of
literature in §2.3, and a discussion of the challenges of using dynamic
programming to find the best Markovian stopping rules for NMOS problems can be
found in Appendix 7.

Finally, we note that prior research in operations research and economics have
studied robust optimal stopping problems in which the goal is to find stopping
rules that perform well under worst-case probability distributions (Bayraktar
and Yao 2014, Riedel 2009) or under worst-case state trajectories (Iancu et
al. 2021). That stream of research differs significantly from ours, as that
stream of research does not consider robust optimization problems that are
approximations of stochastic optimal stopping problems with known probability
distributions.

##  2 Robust Optimization for Stochastic Optimal Stopping

###  2.1 Problem Setting

We consider stochastic optimal stopping problems defined by the following
components:

##### States:

Let $x\equiv(x_{1},\ldots,x_{T})$ denote a sequence of random states, where
the state
$x_{t}\in\mathcal{X}\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\equiv}\;\mathbb{R}^{d}$
in each period $t$ is a random vector of dimension $d$. For example, the state
in each period may represent the prices of multiple assets at that point in
time. The joint probability distribution of this stochastic process is assumed
to be known and accessible through a simulator which generates independent
sample paths of the entire stochastic process.

##### Policies:

Let $\mu\equiv(\mu_{1},\ldots,\mu_{T})$ represent a collection of exercise
policies, where the exercise policy in each period $t$ is a measurable
function of the form
$\mu_{t}:\mathcal{X}\to\\{\textsc{Stop},\textsc{Continue}\\}$. Speaking
intuitively, each exercise policy is a partitioning of the state space into
regions for stopping and continuing. From the exercise policies, the
corresponding Markovian stopping rule
$\tau_{\mu}:\mathcal{X}^{T}\to\\{1,\ldots,T\\}\cup\\{\infty\\}$ is a function
that maps a realization of the stochastic process to a stopping period:

$\displaystyle\tau_{\mu}(x)\triangleq\min\\{t\in\\{1,\ldots,T\\}:\;\mu_{t}(x_{t})=\textsc{Stop}\\}.$

Throughout this paper, a minimization problem with no feasible solutions is
defined equal to $\infty$.111We remark that the Markovian stopping rule
$\tau_{\mu}$ is a non-anticipative control policy, meaning that the event
$\\{\tau_{\mu}(x)=t\\}$ does not depend on the future states
$x_{t+1},\ldots,x_{T}$ for each period $t\in\\{1,\ldots,T\\}$. To see why this
is the case, consider any two realizations of the sequence of random states
$x\equiv(x_{1},\ldots,x_{T})$ and
$x^{\prime}\equiv(x^{\prime}_{1},\ldots,x^{\prime}_{T})$, and suppose for a
given period $t\in\\{1,\ldots,T\\}$ that the two realizations satisfy
$x_{s}=x_{s}^{\prime}$ for all $s\in\\{1,\ldots,t\\}$. Then we readily observe
from our definition of Markovian stopping rules and from algebra that the
Markovian stopping rule $\tau_{\mu}$ satisfies $\tau_{\mu}(x)=t$ if and only
if $\tau_{\mu}(x^{\prime})=t$ is satisfied.

##### Rewards:

Let $g:\\{1,\ldots,T\\}\cup\\{\infty\\}\times\mathcal{X}^{T}\to\mathbb{R}_{+}$
be a known and deterministic function that maps a stopping period and a
realization of the entire stochastic process to a reward. The assumption that
the reward function is nonnegative is common in many applications of optimal
stopping, and we assume throughout the paper that a stochastic process that is
never stopped yields a reward of zero: $g(\infty,x)\equiv 0$. It follows from
the definition of the reward function that the reward from stopping on any
period $t$ may in general depend on the states of the stochastic process in
previous or future time periods.

##### Problem:

With the above notation and inputs, the goal of this paper is to solve
stochastic optimal stopping problems of the form

$\displaystyle\begin{aligned}
&\sup_{\mu}&&\mathbb{E}\left[g(\tau_{\mu}(x),x)\right],\end{aligned}$ (OPT)

where the optimization is taken over the space of all Markovian stopping
rules. In the following sections, we introduce and analyze a new simulation-
based method for solving this class of stochastic dynamic optimization
problems.

###  2.2 The Robust Optimization Approach

Our proposed approach for solving stochastic optimal stopping problems of the
form (OPT) consists of the following steps. We first simulate sample paths of
the stochastic process $x\equiv(x_{1},\ldots,x_{T})$. Let $N$ denote the
number of sample paths, and let the values of the sample paths be denoted by

$x^{i}\equiv(x^{i}_{1},\ldots,x^{i}_{T})\text{ for each }i=1,\ldots,N.$

We assume that the sample paths are independent and identically distributed
realizations of the entire (possibly non-Markovian) stochastic process. We
next choose the following robustness parameter:

$\epsilon\geq 0.$

The purpose of the robustness parameter will become clear momentarily, and a
discussion on how to choose the number of sample paths and the robustness
parameter is deferred until §2.5. With these parameters, let the uncertainty
set around sample path $i$ on period $t$ be defined as

$\displaystyle\mathcal{U}^{i}_{t}\triangleq\left\\{y_{t}\in\mathcal{X}:\;\|y_{t}-x_{t}^{i}\|_{\infty}\leq\epsilon\right\\}.$

For notational convenience, denote the uncertainty set around sample path $i$
across all periods by

$\displaystyle\mathcal{U}^{i}$
$\displaystyle\triangleq\mathcal{U}^{i}_{1}\times\cdots\times\mathcal{U}^{i}_{T}.$

Hence, we observe that the role of the robustness parameter is to control the
size of these sets. Given the sample paths and choice of the robustness
parameter, our approach obtains an approximation of (OPT) by solving the
following robust optimization problem:

$\displaystyle\begin{aligned}
&\sup_{\mu}&&\frac{1}{N}\sum_{i=1}^{N}\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),x^{i}).\end{aligned}$
(RO)

By solving the above robust optimization problem, we obtain exercise policies
$\hat{\mu}\equiv(\hat{\mu}_{1},\ldots,\hat{\mu}_{T})$. These exercise policies
constitute our approximate solution to the stochastic optimal stopping problem
(OPT).

###### Remark 2.1

In the above formulation of the robust optimization problem (RO), the sample
paths in the objective function are allowed to be perturbed by an adversary in
the evaluation of the Markovian stopping rule, $\tau_{\mu}(y)$, but not in the
evaluation of the reward, $g(\cdot,x^{i})$. This is a deviation from the
robust optimization formulation from BSS22, presented below as (RO’), in which
the worst-case reward over each uncertainty set has the form
$\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),y)$:

$\displaystyle\begin{aligned}
&\sup_{\mu}&&\frac{1}{N}\sum_{i=1}^{N}\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),y).\end{aligned}$
(RO’)

The difference between the objective functions of (RO) and (RO’) turns out to
be inconsequential from the perspective of establishing convergence guarantees
(see §2.4) or characterizing the structure of optimal Markovian stopping rules
(see §3). However, (RO) is significantly simpler from the perspective of
algorithm design. For the interested reader, an extended discussion on the
similarities and differences between (RO) and the alternative robust
optimization formulation (RO’) can be found in Appendix 8.

###  2.3 Background and Motivation

In contrast to traditional robust optimization or distributionally robust
optimization, our motivation behind adding adversarial noise to the sample
paths in (RO) is not to find stopping rules which have worst-case performance
guarantees, are attractive in risk-averse settings, or perform well the
presence of an ambiguous probability distribution. Rather, this paper proposes
using robust optimization purely as an algorithmic tool for solving stochastic
optimal stopping problems of the form (OPT) when the joint probability
distributions are known. The present section elaborates on this motivation and
positions our use of robust optimization within the optimal stopping
literature.

For the sake of developing intuition, let us suppose for the moment that the
robustness parameter of the uncertainty sets in (RO) was set equal to zero. In
this case, for any fixed exercise policies $\mu=(\mu_{1},\ldots,\mu_{T})$, the
expected reward of those exercise policies,

$J^{*}(\mu)\triangleq\mathbb{E}\left[g(\tau_{\mu}(x),x)\right],$

would be approximated in (RO) by the sample average approximation:

$\displaystyle\widehat{J}_{N,0}(\mu)\triangleq\frac{1}{N}\sum_{i=1}^{N}\inf_{y\in\mathcal{X}^{T}:\|y-x^{i}\|_{\infty}\leq
0}g(\tau_{\mu}(y),x^{i})=\frac{1}{N}\sum_{i=1}^{N}g(\tau_{\mu}(x^{i}),x^{i}).$

For these fixed exercise policies, we observe that the sample average
approximation is a _consistent_ estimator of the expected reward. In other
words, for the fixed exercise policies $\mu$, it follows from the strong law
of large numbers under relatively mild assumptions222For example, Assumptions
2.4 and 2.4 in §2.4. that $\widehat{J}_{N,0}(\mu)$ will converge almost surely
to $J^{*}(\mu)$ as the number of simulated sample paths is taken to infinity.

However, it is well known in the optimal stopping literature that these
desirable asymptotic properties of $\widehat{J}_{N,0}(\mu)$ are generally not
retained when considering the problem of optimizing over the space of all
exercise policies. For instance, Ciocan and Mišić (2020, EC.1.2) provide
simple examples in which the following statements hold almost surely:

$\displaystyle\lim_{N\to\infty}\sup_{\mu}\widehat{J}_{N,0}(\mu)$
$\displaystyle\gg\sup_{\mu}{J^{*}}(\mu);\quad\lim_{N\to\infty}J^{*}(\hat{\mu}_{N,0})\ll\sup_{\mu}{J^{*}}(\mu).$
(1)

The asymptotic suboptimality of the optimal objective value and optimal
policies for the problem $\sup_{\mu}\widehat{J}_{N,0}(\mu)$ can be intuitively
understood as a type of overfitting. To see why line (1) occurs, we recall for
any fixed choice of exercise policies $\mu$ that the sample average
approximation $\widehat{J}_{N,0}(\mu)$ is an unbiased estimate of the expected
reward $J^{*}(\mu)$. However, when simultaneously considering the space of all
exercise policies, there exists for each $N\in\mathbb{N}$ with high
probability a collection of exercise policies that satisfies
$\widehat{J}_{N,0}(\mu)\gg J^{*}(\mu)$. The problem
$\sup_{\mu}\widehat{J}_{N,0}(\mu)$ will thus be biased towards choosing those
exercise policies, which in general will be suboptimal for the problem
$\sup_{\mu}J^{*}(\mu)$. Because the set of all $\mu$ is an infinite-
dimensional space, the gap between the objective values
$\widehat{J}_{N,0}(\mu)$ and $J^{*}(\mu)$ does not converge to zero uniformly
over the set of all $\mu$ as the number of sample paths tends to infinity.

To circumvent this overfitting in the context of optimal stopping in line (1),
a vast literature has focused on restricting the functional form of exercise
policies to a finite-dimensional space, such as Garcıa (2003), Andersen
(1999), Belomestny (2011a), Gemmrich (2012), Ciocan and Mišić (2020). In this
approach, the choice of the parameterization for the space of exercise
policies must be made very carefully. On one hand, the effective dimension of
the restricted space of exercise policies must be small relative to the number
of simulated sample paths to ensure that the sample average approximation
problem finds the parametric exercise policies that are ‘best-in-class’ with
respect to the stochastic optimal stopping problem (Belomestny 2011a, §3). On
the other hand, the parameterization must be chosen appropriately in order for
the sample average approximation problem to obtain a good approximation of
(OPT). Choosing such an appropriate parameterization “may be counterfactual in
some cases”, as explained by Garcıa (2003, p. 1859), “since we may not have a
good understanding of what the early exercise rule should depend on.”

Our approach, in view of the above discussion, provides an alternative means
to circumvent overfitting. The proposed robust optimization problem allows the
space of exercise policies to remain general, and thus relieves the decision
maker from the need to select and impose a parametric structure on the
exercise policies. Moreover, we show in the following section that our use of
robust optimization provably overcomes the asymptotic overfitting described in
line (1).

###  2.4 Optimality Guarantees

In this subsection, we establish theoretical justification for our combination
of robust optimization and simulation that is presented in §2.2. Specifically,
we strengthen convergence guarantees from BSS22 to the specific problem of
optimal stopping to prove that the optimal objective value and optimal
exercise policies of (RO) converge almost surely to those of the stochastic
optimal stopping problem (OPT) under mild and verifiable conditions.
Establishing these convergence guarantees in the context of optimal stopping
is necessary to ensure that (RO) will provide a high-quality approximation of
the stochastic optimal stopping problem (OPT) when the robustness parameter is
sufficiently small and the number of simulated sample paths is sufficiently
large.

To establish our theoretical results, we make four relatively mild assumptions
on the stochastic optimal stopping problem (OPT). Our first assumption,
denoted below by Assumption 2.4, concerns the structure of the reward
functions in the optimal stopping problem, and can be roughly interpreted as a
requirement that the reward function changes continuously as a function of the
states: $\lim\limits_{\epsilon\to 0}\Delta_{\epsilon}(x)=0$ almost surely,
where

$\displaystyle\Delta_{\epsilon}(x)\triangleq\min_{t\in\\{1,\ldots,T\\}}\left\\{\inf\limits_{y\in\mathcal{X}^{T}:\|y-x\|_{\infty}\leq\epsilon}g(t,y)-g(t,x)\right\\}.$

From a practical standpoint, it is easy to see that Assumption 2.4 holds
whenever the functions $g(1,\cdot),\ldots,g(T,\cdot)$ are continuous, and it
can also hold in important stochastic optimal stopping problems with
discontinuous reward functions.333To illustrate, consider Robbin’s problem
(Bruss 2005), in which the reward functions
$g(t,x)=\text{rank}(x_{t};x_{1},\ldots,x_{T})$ are discontinuous and the
probability distribution is
$x_{1},\ldots,x_{T}\overset{\text{iid}}{\sim}\text{Uniform}[0,1]$. To show
that Assumption 2.4 is satisfied, we observe that the random variable
$\bar{\epsilon}\triangleq\min_{s<t}|x_{s}-x_{t}|$ is strictly positive with
probability one, which implies that $\Delta_{\epsilon}(x)=0$ for all
$\epsilon<\bar{\epsilon}$.

Our second and third assumptions concern the structure of the probability
distribution in (OPT). Specifically, Assumption 2.4 enforces that the
stochastic process has a light tail, and Assumption 2.4 says the stochastic
process is drawn from a continuous probability distribution. The stochastic
process satisfies $\mathbb{E}[\textnormal{exp}(\|x\|_{\infty}^{a})]<\infty$
for some $a>1$. The stochastic process $x\equiv(x_{1},\ldots,x_{T})$ has a
probability density function. Let us reflect on the practical restrictiveness
of these two assumptions. The second assumption, Assumption 2.4, is a standard
light-tail assumption on the stochastic process which is satisfied, for
example, if the stochastic process is bounded or has a multivariate normal
distribution. This assumption greatly simplifies our analysis, as it allows us
to invoke a convergence result by BSS22 in our proofs (see Appendix 9). We
impose the third assumption, Assumption 2.4, to ensure that there exist
arbitrarily near-optimal Markovian stopping rules for (OPT) that satisfy a
certain technical continuity structure that we can exploit in our proof. As
far as we can tell, these assumptions on the probability distribution are
relatively mild and routinely satisfied in applications of optimal stopping in
the context of the options pricing literature. Nonetheless, we do not preclude
the possibility that these assumptions on the probability distribution can be
weakened while still establishing convergence guarantees.

Our fourth and final assumption imposes boundedness on the reward function.
This assumption, presented below as Assumption 2.4, leads to a considerably
simpler proof and statement of the results, but can generally be relaxed to
reward functions bounded above by an integrable, Lipschitz-continuous
function. The reward function satisfies $0\leq g(t,y)\leq U$ for all
$t\in\\{1,\ldots,T\\}$ and $y\in\mathcal{X}^{T}$.

We emphasize that each of the aforementioned four assumptions on the
stochastic optimal stopping problem (OPT) can be verified _a priori_. In
particular, they do not require any knowledge of the structure of optimal
exercise policies for the stochastic optimal stopping problem (OPT). As a
result, each of these assumptions can be verified using the information
typically available in practice. We note that these assumptions are
considerably weaker than those in BSS22, which require knowledge of the
structure of optimal control policies to establish convergence results.

Under the above conditions, the following theorems provide justification for
using the robust optimization problem as a proxy for the stochastic optimal
stopping problem. In a nutshell, the following Theorems 2.2-2.4 show that (RO)
will, for all sufficiently small choices of the robustness parameter and all
sufficiently large choices of the number of simulated sample paths, yield a
near-optimal approximation of (OPT). Stated another way, the following
theorems show that our use of robust optimization provably overcomes the
asymptotic overfitting described in line (1) of §2.3. While the following
theorems do not specify how to choose the robustness parameter and number of
simulated sample paths for any particular optimal stopping problem, we provide
guidance (§2.5) and numerical evidence (§5) which suggest that these
parameters can be found effectively in practice. A discussion of the technical
innovations as well the proofs of the following theorems in this subsection
can be found in Appendix 9.

Our first theorem shows that the optimal objective value of the robust
optimization problem (RO) will converge almost surely to that of the
stochastic problem (OPT) as the robustness parameter tends to zero and the
number of sample paths tends to infinity. In the following result, we use the
notation $\widehat{J}_{N,\epsilon}(\mu)$ to denote the objective value of the
robust optimization problem (RO) corresponding to exercise policies $\mu$.

######  Theorem 2.2 (Consistency of optimal objective value)

Under Assumptions 2.4, 2.4, 2.4, and 2.4,

$\lim\limits_{\epsilon\to
0}\lim\limits_{N\to\infty}\sup\limits_{\mu}\widehat{J}_{N,\epsilon}(\mu)=\sup\limits_{\mu}J^{*}(\mu)\quad\textnormal{almost
surely}.$

Our second theorem shows that the expected reward of the optimal exercise
policies for the robust optimization problem (RO) will converge almost surely
to the optimal objective value of the stochastic problem (OPT). We let
$\hat{\mu}_{N,\epsilon}$ denote optimal exercise policies for (RO), and we
remark that the existence of optimal exercise policies for the robust
optimization problem will be established in §3.

######  Theorem 2.3 (Consistency of optimal policies)

Under Assumptions 2.4, 2.4, 2.4, and 2.4,

$\displaystyle\lim_{\epsilon\to
0}\liminf_{N\to\infty}J^{*}(\hat{\mu}_{N,\epsilon})=\lim_{\epsilon\to
0}\limsup_{N\to\infty}J^{*}(\hat{\mu}_{N,\epsilon})$
$\displaystyle=\sup_{\mu}J^{*}(\mu)\quad\textnormal{almost surely}.$

Because we will develop algorithms that solve the robust optimization problem
approximately as well as exactly, it is imperative for us to have theoretical
guarantees that hold for any Markovian stopping rule that can be found by the
robust optimization problem. To this end, our third and final theorem of this
section shows that the (in-sample) robust objective value will asymptotically
provide a low-bias estimate of the expected reward, and this bound holds
uniformly over all exercise policies. The result yields theoretical assurance,
provided that the robustness parameter is sufficiently small and the number of
sample paths is sufficiently large, that searching for exercise policies with
high robust objective values $\widehat{J}_{N,\epsilon}(\mu)$ will typically
result in exercise policies with high expected rewards $J^{*}(\mu)$.

######  Theorem 2.4 (Asymptotic low-bias)

Under Assumptions 2.4, 2.4, and 2.4,

$\lim\limits_{\epsilon\to
0}\liminf\limits_{N\to\infty}\inf\limits_{\mu}\left\\{J^{*}(\mu)-\widehat{J}_{N,\epsilon}(\mu)\right\\}\geq
0\quad\textnormal{almost surely}.$

###  2.5 Implementation Details

In anticipation of algorithmic techniques for solving the robust optimization
problem (RO) in the remainder of the paper, it remains to be specified how the
parameters of the robust optimization problem (the number of simulated sample
paths $N\in\mathbb{N}$ and the robustness parameter $\epsilon\geq 0$) should
be selected in practice. For the sake of concreteness, we conclude §2 by
briefly providing guidance for choosing these parameters and applying the
robust optimization approach in practice. The procedures described below are
formalized in Algorithm LABEL:fig:heuristic and implemented in our numerical
experiments in §5.

As described previously, this paper addresses stochastic optimal stopping
problems in which the probability distributions are known. Consequently, the
decision-maker is granted flexibility in choosing the number of sample paths
$N$ to simulate. On one hand, we have established in the previous section that
larger choices of the number of simulated sample paths will generally lead to
tighter approximations of the stochastic optimal stopping problem. On the
other hand, larger choices of $N$ require a greater computation cost in
performing the Monte-Carlo simulation and creates a robust optimization
problem of a larger size. To balance these tradeoffs in particular
applications, we recommend using a straightforward procedure of starting out
with a small choice of $N$ and iteratively increasing the number of simulated
sample paths until the total computational cost meets the allocated
computational budget.

Given a fixed number of sample paths, the choice of the robustness parameter
$\epsilon\geq 0$ can have a significant impact on the policies produced by the
robust optimization problem. To this end, we recommend solving (RO) over a
grid of possible choices for the robustness parameter. Because the probability
distribution is known, we can generate a second set of ‘validation’ sample
paths to select the best choice of the robustness parameter. Specifically, for
each choice of the robustness parameter, one solves the robust optimization
problem to obtain exercise policies. The expected reward of the exercise
policies is then estimated using the validation set of sample paths. Finally,
we select the value of the robustness parameter (and the corresponding
exercise policies) which maximizes the average reward with respect to the
validation set.

In summary, we have described straightforward and easy-to-implement heuristics
for choosing the parameters of the robust optimization problem. Applying the
heuristics and solving the robust optimization problem yields exercise
policies for the stochastic optimal stopping problem, and an unbiased estimate
of the expected reward of these exercise policies can similarly be obtained by
simulating a set of ‘testing’ sample paths (see Algorithm
LABEL:fig:heuristic). Because the exercise policies obtained from the robust
optimization problem are feasible for the stochastic optimal stopping problem,
the expected reward of these exercise policies is thus a lower bound on the
optimal objective value of the stochastic optimal stopping problem. Finally,
we remark that under a stronger assumption in which one has the ability to
perform conditional Monte Carlo simulation, the exercise policies obtained
from solving the robust optimization problem can be combined with the method
of Andersen and Broadie (2004) to obtain an upper bound on the optimal
objective value of the stochastic optimal stopping problem.

##  3 Characterization of Optimal Markovian Stopping Rules

In §2, we showed that our combination of robust optimization and simulation
(§2.2) can yield an arbitrarily close approximation of the stochastic optimal
stopping problem (§2.4-§2.5). In this section, we develop the key technical
result of this paper, Theorem 3.1, which will enable us to design exact and
heuristic algorithms for solving the robust optimization problem.
Specifically, Theorem 3.1 establishes the existence and characterizes the
structure of optimal Markovian stopping rules for the robust optimization
problem (RO). In §3.1 and §3.2, we present the statement of Theorem 3.1 and
provide a sketch of its proof. In §3.3, we use the characterization of optimal
Markovian stopping rules to transform (RO) from an optimization problem over
an infinite-dimensional space of exercise policies into a finite-dimensional
optimization problem over integer decision variables
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$ (Theorem 3.7).

###  3.1 Statement of Theorem 3.1

We begin by introducing the notation that will be used in our characterization
of the structure of optimal Markovian stopping rules for the robust
optimization problem (RO). Consider any instance of the robust optimization
problem (RO). For any choice of integers
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$, we define
$\mu^{\sigma^{1}\cdots\sigma^{N}}\equiv(\mu^{\sigma^{1}\cdots\sigma^{N}}_{1},\ldots,\mu^{\sigma^{1}\cdots\sigma^{N}}_{T})$
as the exercise policy that satisfies the following equality for each period
$t\in\\{1,\ldots,T\\}$ and each state $y_{t}\in\mathcal{X}$:

$\displaystyle\mu_{t}^{\sigma^{1}\cdots\sigma^{N}}(y_{t})\triangleq\begin{cases}\textsc{Stop},&\text{
if }y_{t}\in\bigcup\limits_{i:\;\sigma^{i}=t}\mathcal{U}^{i}_{t},\\\
\textsc{Continue},&\text{ if
}y_{t}\notin\bigcup\limits_{i:\;\sigma^{i}=t}\mathcal{U}^{i}_{t}.\end{cases}$
(2)

To develop intuition of the above exercise policy, we remark for each sample
path $i\in\\{1,\ldots,N\\}$ that the integer $\sigma^{i}\in\\{1,\ldots,T\\}$
can be interpreted as a selection of one of the uncertainty sets
$\mathcal{U}^{i}_{1},\ldots,\mathcal{U}^{i}_{T}$. Specifically, if
$\sigma^{i}=t$, then we observe from the above definition that the exercise
policy will satisfy
$\mu_{t}^{\sigma^{1}\cdots\sigma^{N}}(y_{t})=\textsc{Stop}$ for all
$y_{t}\in\mathcal{U}^{i}_{t}$. More generally, we observe that the exercise
policy will satisfy
$\mu_{t}^{\sigma^{1}\cdots\sigma^{N}}(y_{t})=\textsc{Stop}$ for period
$t\in\\{1,\ldots,T\\}$ and state $y_{t}\in\mathcal{X}$ if and only if there
exists an uncertainty set $\mathcal{U}^{i}_{t}$ such that the state is
contained in the uncertainty set, $y_{t}\in\mathcal{U}^{i}_{t}$, and the
integer corresponding to the $i$th sample path is equal to the current period,
$\sigma^{i}=t$. A visualization of the exercise policy
$\mu^{\sigma^{1}\cdots\sigma^{N}}$ is found in Figure 1.

(a) Sample paths and their corresponding uncertainty sets

(b) Stopping regions of
$\mu^{23}\equiv(\mu^{23}_{1},\mu^{23}_{2},\mu^{23}_{3})$

Figure 1: Visualization of a robust optimization problem with $N=2$ sample
paths and $T=3$ periods. The values of the sample paths in the robust
optimization problem are given by
$x^{1}=(x^{1}_{1},x^{1}_{2},x^{1}_{3})=(8,7,6)$ and
$x^{2}=(x^{2}_{1},x^{2}_{2},x^{2}_{3})=(3,4,3)$, the state space is the real
numbers ($\mathcal{X}=\mathbb{R}^{1}$) and the robustness parameter is
$\epsilon=2.0$. The green intervals correspond to the uncertainty sets
$\mathcal{U}^{1}_{1}=[6,10]$, $\mathcal{U}^{1}_{2}=[5,9]$, and
$\mathcal{U}^{1}_{3}=[4,8]$ around the first sample path, and the blue
intervals correspond to the uncertainty sets $\mathcal{U}^{2}_{1}=[1,5]$,
$\mathcal{U}^{2}_{2}=[2,6]$, and $\mathcal{U}^{2}_{3}=[1,5]$ around the second
sample path. In Figure 1a, we visualize the two sample paths using dashed
lines. In Figure 1b, we visualize the stopping regions induced by the exercise
policy
$\mu^{\sigma^{1}\sigma^{2}}\equiv(\mu^{\sigma^{1}\sigma^{2}}_{1},\mu^{\sigma^{1}\sigma^{2}}_{2},\mu^{\sigma^{1}\sigma^{2}}_{3})$
in the case where the integers $\sigma^{1},\sigma^{2}\in\\{1,2,3\\}$ satisfy
$\sigma^{1}=2$ and $\sigma^{2}=3$. Specifically, the black vertical lines in
Figure 1b show the stopping regions
$\\{y_{1}\in\mathbb{R}:\mu^{23}_{1}(y_{1})=\textsc{Stop}\\}$,
$\\{y_{2}\in\mathbb{R}:\mu^{23}_{2}(y_{2})=\textsc{Stop}\\}$, and
$\\{y_{3}\in\mathbb{R}:\mu^{23}_{3}(y_{3})=\textsc{Stop}\\}$ induced by
$\mu^{23}\equiv(\mu^{23}_{1},\mu^{23}_{2},\mu^{23}_{3})$.

Let the set of all exercise polices generated by integers
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$ be denoted by

$\displaystyle\mathcal{M}\triangleq\left\\{\mu\equiv(\mu_{1},\ldots,\mu_{T}):\textnormal{
there exist }\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}\textnormal{ such
that }\mu=\mu^{\sigma^{1}\cdots\sigma^{N}}\right\\}.$

It is clear from the above definition that each of the exercise policies
$\mu\in\mathcal{M}$ is parameterized by integers
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$; thus, we readily observe
that the cardinality of $\mathcal{M}$ is always finite and upper bounded by
$T^{N}$. Moreover, we observe that the definition of the set of exercise
policies is sample-path dependent, in the sense that $\mathcal{M}$ depends on
the number and realizations of the simulated sample paths
$x^{1}\equiv(x^{1}_{1},\ldots,x^{1}_{T}),\ldots,x^{N}\equiv(x^{N}_{1},\ldots,x^{N}_{T})$
and on the choice of the robustness parameter $\epsilon\geq 0$. In Figure
LABEL:fig:exhaustive, we present a visualization of the set of exercise
policies $\mathcal{M}$.

In view of the above notation, we now present our main result:

###### Theorem 3.1

There exists $\mu\in\mathcal{M}$ that is optimal for (RO).

The above theorem shows that there always exists an exercise policy in the set
$\mathcal{M}$ that is optimal for the robust optimization problem. The result
is significant because it will allow us to transform (RO) from an optimization
problem over an infinite-dimensional space of exercise policies into a finite-
dimensional combinatorial optimization problem over the set $\mathcal{M}$.
Moreover, Theorem 3.1 is important because it establishes the existence of
optimal Markovian stopping rules for the robust optimization problem (RO).

###  3.2 Proof Sketch of Theorem 3.1

Our proof of Theorem 3.1 follows an exchange argument that is rooted in a
technique that we refer to as _pruning_. The technique consists of starting
with an initial exercise policy $\mu\equiv(\mu_{1},\ldots,\mu_{T})$ and then
modifying the exercise policy to reduce the size of the stopping regions
$\\{y_{t}\in\mathcal{X}:\mu_{t}(y_{t})=\textsc{Stop}\\}$. By carefully pruning
(_i.e._ , reducing the size of) the stopping regions of the initial exercise
policy $\mu\equiv(\mu_{1},\ldots,\mu_{T})$, we will show that any initial
exercise policy can be transformed into a new exercise policy
$\mu^{\prime}\in\mathcal{M}$ with the same or better objective value in the
robust optimization problem (RO).

To discuss the proof in greater detail, we begin by stating two preliminary
lemmas:

###### Lemma 3.2

The optimal objective value of (RO) is equal to the optimal objective value of

$\displaystyle\begin{aligned}
&\sup_{\mu}&&\frac{1}{N}\sum_{i=1}^{N}\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),x^{i})\\\
&\textnormal{subject to}&&\textnormal{for each
}i\in\\{1,\ldots,N\\},\;\textnormal{ there exists }t\in\\{1,\ldots,T\\}\\\
&&&\textnormal{such that }\mu_{t}(y_{t})=\textsc{Stop}\textnormal{ for all
}y_{t}\in\mathcal{U}^{i}_{t}.\end{aligned}$ (ROT)

###### Lemma 3.3

Consider any $\mu\equiv(\mu_{1},\ldots,\mu_{T})$ that satisfies the
constraints of (ROT), and define

$\displaystyle\sigma^{i}\triangleq\min\left\\{t\in\\{1,\ldots,T\\}:\;\mu_{t}(y_{t})=\textsc{Stop}\textnormal{
for all }y_{t}\in\mathcal{U}^{i}_{t}\right\\}\quad\forall
i\in\\{1,\ldots,N\\}.$

Then the following equality holds for each $i\in\\{1,\ldots,N\\}$:

$\displaystyle\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),x^{i})=\min\limits_{t\in\\{1,\ldots,\sigma^{i}\\}}\left\\{g(t,x^{i}):\textnormal{there
exists }y_{t}\in\mathcal{U}^{i}_{t}\textnormal{ such that
}\mu_{t}(y_{t})=\textsc{Stop}\right\\}.$

The first preliminary lemma shows that a structural constraint can be imposed
onto the exercise policies in the robust optimization problem (RO) without any
loss of optimality. In particular, Lemma 3.2 says that we can restrict to
exercise policies in which the resulting Markovian stopping rule satisfies
$\tau_{\mu}(y)<\infty$ for each trajectory $y$ in each of the uncertainty sets
$\mathcal{U}^{1},\ldots,\mathcal{U}^{N}$.444Recall that the definition of a
Markovian stopping rule is
$\tau_{\mu}(y)\triangleq\min\\{t\in\\{1,\ldots,T\\}:\mu_{t}(y_{t})=\textsc{Stop}\\}$.
Therefore, we readily observe for each sample path $i$ that an exercise policy
$\mu\equiv(\mu_{1},\ldots,\mu_{T})$ satisfies $\tau_{\mu}(y)<\infty$ for all
$y\in\mathcal{U}^{i}\equiv\mathcal{U}^{i}_{1}\times\cdots\times\mathcal{U}^{i}_{T}$
if and only if there exists a period $t\in\\{1,\ldots,T\\}$ that satisfies
$\mu_{t}(y_{t})=\textsc{Stop}$ for all $y_{t}\in\mathcal{U}^{i}_{t}$.  The
second preliminary lemma develops a convenient representation of the objective
function of (ROT). Specifically, Lemma 3.3 shows for each sample path $i$ that
the quantity $\inf_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),x^{i})$ in the
objective function of (ROT) is equal to the the minimum $g(t,x^{i})$ among all
periods $t\in\\{1,\ldots,\max_{y\in\mathcal{U}^{i}}\tau_{\mu}(y)\\}$ for which
the uncertainty set $\mathcal{U}^{i}_{t}$ has a nonempty intersection with the
stopping region $\\{y_{t}\in\mathcal{X}_{t}:\mu_{t}(y_{t})=\textsc{Stop}\\}$.

Equipped with the above preliminary lemmas, we now formally define the pruning
technique that underpins our proof of Theorem 3.1.

###### Definition 3.4

_Let $\mu$ be an exercise policy that is feasible for (ROT). We say that
$\mu^{\prime}$ is a _pruned version_ of $\mu$ if the stopping regions of
$\mu^{\prime}$ are a subset of the stopping regions of $\mu$, _i.e._ ,_

$\displaystyle\left\\{y_{t}\in\mathcal{X}:\mu_{t}^{\prime}(y_{t})=\textsc{Stop}\right\\}\subseteq\left\\{y_{t}\in\mathcal{X}:\mu_{t}(y_{t})=\textsc{Stop}\right\\}\quad\forall
t\in\\{1,\ldots,T\\},$

_and if the following equality holds for each sample path
$i\in\\{1,\ldots,N\\}$:_

$\displaystyle\min\left\\{t\in\\{1,\ldots,T\\}:\;\mu_{t}^{\prime}(y_{t})=\textsc{Stop}\;\forall
y_{t}\in\mathcal{U}^{i}_{t}\right\\}=\min\left\\{t\in\\{1,\ldots,T\\}:\;\mu_{t}(y_{t})=\textsc{Stop}\;\forall
y_{t}\in\mathcal{U}^{i}_{t}\right\\}.$

Speaking intuitively, an exercise policy $\mu^{\prime}$ is a pruned version of
an exercise policy $\mu$ if the stopping regions of $\mu^{\prime}$ are a
subset of the stopping regions of $\mu$ and if
$\max_{y\in\mathcal{U}^{i}}\tau_{\mu}(y)=\max_{y\in\mathcal{U}^{i}}\tau_{\mu^{\prime}}(y)$
for each sample path $i$. The significance of Definition 3.4 in combination
with Lemma 3.3 is established by the following Lemma 3.5; specifically, the
following Lemma 3.5 shows that if $\mu^{\prime}$ is a pruned version of $\mu$,
then the robust objective value associated with $\mu^{\prime}$ is greater than
or equal to the robust objective value associated with $\mu$.

###### Lemma 3.5

If $\mu^{\prime}$ is a pruned version of $\mu$, then
$\inf\limits_{y\in\mathcal{U}^{i}}g(\tau_{\mu^{\prime}}(y),x^{i})\geq\inf\limits_{y\in\mathcal{U}^{i}}g(\tau_{\mu}(y),x^{i})$
$\forall i\in\\{1,\ldots,N\\}$.

In summary, we have shown in Lemma 3.2 that the robust optimization problem
(RO) is equivalent to the robust optimization problem (ROT). Moreover, for
every initial exercise policy $\mu$ that is feasible for (ROT), Lemma 3.5
shows that the objective value associated with $\mu$ is less than or equal to
the objective value associated with every exercise policy $\mu^{\prime}$ that
is a pruned version of $\mu$. The final step of our proof of Theorem 3.1 is
given by the following lemma:

###### Lemma 3.6

Consider any $\mu\equiv(\mu_{1},\ldots,\mu_{T})$ that satisfies the
constraints of (ROT), and define

$\displaystyle\sigma^{i}\triangleq\min\left\\{t\in\\{1,\ldots,T\\}:\;\mu_{t}(y_{t})=\textsc{Stop}\textnormal{
for all }y_{t}\in\mathcal{U}^{i}_{t}\right\\}\quad\forall
i\in\\{1,\ldots,N\\}.$

Then $\mu^{\sigma^{1}\cdots\sigma^{N}}$ is a pruned version of $\mu$.

The above lemma implies that if $\mu$ is a feasible exercise policy for (ROT),
then there always exists an exercise policy $\mu^{\prime}\in\mathcal{M}$ that
is a pruned version of $\mu$. Hence, we conclude from Lemmas 3.2, 3.5, and 3.6
that there always exists an exercise policy $\mu\in\mathcal{M}$ that is
optimal for (RO). The omitted details of the proofs of Theorem 3.1 and Lemmas
3.2, 3.3, 3.5, and 3.6 are found in Appendix 10.

###  3.3 Reformulation of (RO) as a Finite-Dimensional Optimization Problem

We conclude §3 by using our characterization of optimal Markovian stopping
rules (Theorem 3.1) to develop a reformulation of (RO) as a finite-dimensional
optimization problem over integer decision variables
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$. The reformulation from this
subsection is important because it establishes a natural combinatorial
interpretation of the robust optimization problem (RO), which will provide the
foundation for our algorithmic developments in §4.1 and §4.2.

Our finite-dimensional reformulation is presented as the following
optimization problem (IP):

$\displaystyle\underset{\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}}{\textnormal{maximize}}\;\frac{1}{N}\sum_{i=1}^{N}\min_{t\in\\{1,\ldots,\sigma^{i}\\}}\left\\{g(t,x^{i}):\textnormal{
there exists }j\textnormal{ such that
}\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\textnormal{ and
}\sigma^{j}=t\right\\}.$ (IP)

The decision variables in the above optimization problem are the integers
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$. We further observe that the
inner minimization problems in (IP) involves constraints that depend on
whether the intersection $\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}$ is
nonempty for each pair of sample paths $i,j$ and each period $t$. An important
insight is that these intersections can be precomputed; that is, given the
construction of the uncertainty sets from §2.2, we can efficiently precompute
the set of all pairs of sample paths and periods that satisfy
$\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset$.555We recall from
§2.2 that each uncertainty set $\mathcal{U}^{i}_{t}\subseteq\mathbb{R}^{d}$ is
a hypercube. Thus, determining whether
$\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}$ is nonempty for each pair of
sample paths $i,j\in\\{1,\ldots,N\\}$ and each period $t\in\\{1,\ldots,T\\}$
can be precomputed in a total of $\mathcal{O}(N^{2}Td)$ computation time.

The equivalence of (IP) and (RO) is formalized by the following theorem.

###### Theorem 3.7

The optimal objective value of (IP) is equal to the optimal objective value of
(RO). Furthermore, for any choice of integers
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$,

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\inf_{y\in\mathcal{U}^{i}}g\left(\tau_{\mu^{\sigma^{1}\cdots\sigma^{N}}}(y),x^{i}\right)$
$\displaystyle\geq\frac{1}{N}\sum_{i=1}^{N}\min_{t\in\\{1,\ldots,\sigma^{i}\\}}\left\\{g(t,x^{i}):\textnormal{
there exists }j\textnormal{ such that
}\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\textnormal{ and
}\sigma^{j}=t\right\\}.$

The proof of Theorem 3.7, which is found in Appendix 10, follows readily from
Lemmas 3.2, 3.3, 3.5, and 3.6. Stated in words, the above theorem shows, for
any given feasible solution $\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$
for the optimization problem (IP), that the objective value associated with
the exercise policy $\mu^{\sigma^{1}\cdots\sigma^{N}}$ in the optimization
problem (RO) is greater than or equal to the objective value associated with
$\sigma^{1},\ldots,\sigma^{N}$ in the optimization problem (IP). Because the
above theorem also shows that optimal objective value of (RO) is equal to the
optimal objective value of (IP), it follows immediately from Theorem 3.7 that
any optimal solution for (IP) can be transformed into an optimal solution for
(RO).

##  4 Computation of Optimal Markovian Stopping Rules

In this section, we use our characterization of optimal Markovian stopping
rules from §3 to develop exact and heuristic algorithms for solving the robust
optimization problem (RO). In §4.1, we use this finite-dimensional
reformulation of (RO) to design exact algorithms and hardness results for
solving the robust optimization problem. In §4.2, we propose and analyze an
efficient heuristic algorithm for approximately solving the robust
optimization problem. We emphasize that the algorithms and analysis in this
section are general and do not require any of the assumptions that were made
in §2.4.

###  4.1 Exact Algorithm

In §3.3, we developed a finite-dimensional optimization problem (IP) with
integer decision variables $\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$
that is equivalent to the infinite-dimensional robust optimization problem
(RO). Specifically, we showed that the optimal objective values of the
optimization problems (IP) and (RO) are always equal. Moreover, we readily
observe from Theorem 3.7 that any optimal solution for (IP) can further be
transformed into an optimal exercise policy for (RO) through the
transformation described by line (2) in §3. Harnessing this optimization
problem (IP), we now design exact algorithms for solving the robust
optimization problem (RO).

We begin our discussion by considering the case of (IP) when there are two
periods. Optimal stopping problems with two periods has been studied in the
literature as a testbed for understanding the complexity of solving optimal
stopping problems, e.g., Glasserman and Yu (2004). We defer the proof of the
following result to §4.2.4.

###### Theorem 4.1

If $T=2$, then (IP) can be solved in $\mathcal{O}(N^{3})$ time.

Continuing with our discussion on the theoretical tractability of (IP), we
next consider problems with three or more periods. The following negative
result shows that the computational tractability of (IP) in the case of two
periods established in Theorem 4.1 does not generally extend to optimal
stopping problems with three or more periods. The proof of the following
result, which is found in Appendix 11, consists of a reduction from MIN-2-SAT,
which is shown to be strongly NP-hard by Kohli et al. (1994).

###### Theorem 4.2

(IP) is strongly NP-hard for any fixed $T\geq 3$.

Motivated by the above hardness result, we proceed to develop an exact
algorithm for solving (IP) by reformulating it as a zero-one bilinear program
over totally unimodular constraints. This reformulation, presented below as
(BP), is valuable for three primary reasons. First, it is well known that
bilinear programs can generally be transformed into mixed-integer linear
optimization problems using linearization techniques; see Appendix 12 for
details. Thus, for robust optimization problems with small values of $T$ and
$N$, the exact reformulation from this section can directly solved by off-the-
shelf commercial optimization software. Second, the exact reformulation (BP)
will provide the foundation for our heuristic algorithm for the robust
optimization problem in the following §4.2. Third, the exact reformulation
(BP) is relatively compact, requiring only $\mathcal{O}(NT^{2})$ decision
variables and $\mathcal{O}(NT(T+N))$ constraints. In particular, the mild
dependence of the size of the exact reformulation on the number of simulated
sample paths is attractive in practical settings where the number of simulated
sample paths is much larger than the number of time periods.

Our reformulation of (IP) as a zero-one bilinear program over totally
unimodular constraints (BP) requires the following additional notation. For
each sample path $i\in\\{1,\ldots,N\\}$, we define the following set:

$\displaystyle{\mathcal{K}}^{i}$ $\displaystyle\triangleq\left\\{\kappa:\text{
there exists }t\in\\{1,\ldots,T\\}\text{ such that
}g(t,x^{i})=\kappa\right\\}\cup\\{0\\}.$

The above set can be interpreted as the set of distinct values among
$g(1,x^{i}),\ldots,g(T,x^{i})$ and $0$. For notational convenience, let the
elements of $\mathcal{K}^{i}$ be indexed in ascending order,
$\kappa^{i}_{1}<\cdots<\kappa^{i}_{|\mathcal{K}^{i}|}$, and let
$L^{i}_{t}\in\\{1,\ldots,|\mathcal{K}^{i}|\\}$ be defined for each period
$t\in\\{1,\ldots,T\\}$ as the unique index that satisfies the equality
$g(t,x^{i})=\kappa^{i}_{L^{i}_{t}}$. We readily observe that the quantities
$L^{i}_{t}$ and $\kappa^{i}_{1},\ldots,\kappa^{i}_{L^{i}_{t}}$ can be
efficiently precomputed for each sample path $i$ and period $t$. With the
above notation, our exact reformulation of (IP) is stated as follows:

$\displaystyle\underset{b,w}{\textnormal{maximize}}$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{\ell=1}^{L^{i}_{t}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})b^{i}_{t}(1-w^{i}_{t\ell})$
(BP) subject to $\displaystyle\begin{aligned} &w^{i}_{t\ell}\leq
w^{i}_{t+1,\ell}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T-1\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}\\\
&w^{i}_{t\ell}\leq w^{i}_{t,\ell+1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|-1\\}\\\
&b^{i}_{t}\leq w^{i}_{t+1,1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T-1\\}\\\ &b^{j}_{t}\leq
w^{i}_{t\ell}&&\textnormal{for all }i,j\in\\{1,\ldots,N\\}\textnormal{ and
}t\in\\{1,\dots,T\\}\textnormal{ such that }g(t,x^{i})=\kappa^{i}_{\ell}\\\
&&&\textnormal{and }\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\\\
&b^{i}_{t}\in\\{0,1\\}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\dots,T\\}\\\
&w^{i}_{t\ell}\in\mathbb{R}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}.\end{aligned}$

Let us reflect on the structure of the above optimization problem. First, we
observe that (BP) is a bilinear program because its objective function is the
sum of products of decision variables $b^{i}_{t}$ and $w^{i}_{t\ell}$. Second,
we observe that the bilinear program (BP) is comprised of binary decision
variables $b^{i}_{t}\in\\{0,1\\}$ as well as continuous decision variables
$w^{i}_{t\ell}\in\mathbb{R}$. In particular, we readily observe from the
structure of the constraints of (BP) and from the fact that each quantity
$\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell}$ is strictly positive that there always
exists an optimal solution for (BP) in which each decision variable
$w^{i}_{t\ell}$ is equal to zero or one. For this reason, (BP) will be
henceforth referred to as a zero-one bilinear program without any ambiguity in
terminology. Finally, we remark that the constraints of (BP) are totally
unimodular, which implies that every extreme point of the polyhedron defined
by the constraints of (BP) is integral (Conforti et al. 2014, §4.2). We will
utilize the total unimodularity of (BP) in §4.2 when designing our heuristic
algorithm for the robust optimization problem.

We next show that (BP) is indeed equivalent to (IP). In the following Theorem
4.3, we establish this equivalence and show that any feasible solution for
(BP) can be transformed into a feasible solution for (IP) with the same or
greater objective value. The proof of Theorem 4.3, as well as the proofs of
the subsequent lemmas in §4.1, are found in Appendix 13.

###### Theorem 4.3

The optimal objective value of (IP) is equal to the optimal objective value of
(BP). Furthermore, if $b,w$ is a feasible solution for (BP) and if
$\sigma^{i}\triangleq\min\\{\min\\{t\in\\{1,\ldots,T\\}:b^{i}_{t}=1\\},T\\}$
for each $i\in\\{1,\ldots,N\\}$, then

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\min_{t\in\\{1,\ldots,\sigma^{i}\\}}\left\\{g(t,x^{i}):\textnormal{
there exists }j\textnormal{ such that
}\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\textnormal{ and
}\sigma^{j}=t\right\\}$
$\displaystyle\geq\frac{1}{N}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{\ell=1}^{L^{i}_{t}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})b^{i}_{t}(1-w^{i}_{t\ell}).$

The above theorem shows for each sample path $i$ that the decision variables
$b^{i}_{1},\ldots,b^{i}_{T}\in\\{0,1\\}$ in (BP) can be interpreted as an
encoding of an integer $\sigma^{i}\in\\{1,\ldots,T\\}$ for (IP). In other
words, given any feasible solution $b,w$ for (BP), Theorem 4.3 shows how to
construct a feasible solution
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$ for (IP) such that the
objective value associated with $\sigma^{1},\ldots,\sigma^{N}$ in the
optimization problem (IP) is greater than or equal to the objective value
associated with $b,w$ in the optimization problem (BP). Because the above
theorem also shows that optimal objective value of (IP) is equal to the
optimal objective value of (BP), it follows immediately from Theorem 4.3 that
any optimal solution for (BP) can be transformed into an optimal solution for
(IP).

We conclude §4.1 by providing intuition for the role of the decision variables
$w^{i}_{t\ell}$ in the optimization problem (BP). To do this, we outline the
two key steps of our derivation of (BP). Indeed, recall that the optimization
problem (IP) is comprised of integer decision variables
$\sigma^{1},\ldots,\sigma^{N}\in\\{1,\ldots,T\\}$. The first step in our
derivation of (BP) is to introduce binary decision variables
$b^{i}_{1},\ldots,b^{i}_{T}\in\\{0,1\\}$ for each sample path $i$ to represent
each integer decision variable $\sigma^{i}$. Specifically, consider the
following nonlinear binary optimization problem:

$\displaystyle\underset{b}{\textnormal{maximize}}$
$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\psi^{i}\left(b\right)$ (BP-1) subject
to $\displaystyle b^{i}_{t}\in\\{0,1\\}\quad\textnormal{for all
}i\in\\{1,\ldots,N\\}\text{ and }t\in\\{1,\ldots,T\\},$

where the function $\psi^{i}\left(b\right)$ is defined for each sample path
$i\in\\{1,\ldots,N\\}$ as

$\displaystyle\psi^{i}\left(b\right)$
$\displaystyle\triangleq\sum_{t=1}^{T}b^{i}_{t}\left(\prod_{s=1}^{t-1}(1-b^{i}_{s})\right)\min_{s\in\\{1,\ldots,t\\}}\left\\{g(s,x^{i}):\textnormal{there
exists }j\textnormal{ such that
}\mathcal{U}^{i}_{s}\cap\mathcal{U}^{j}_{s}\neq\emptyset\textnormal{ and
}b^{j}_{s}=1\right\\}.$

To make sense of the function $\psi^{i}(b)$, we remark that each quantity
$b^{i}_{t}\prod_{s=1}^{t-1}(1-b^{i}_{s})$ will be equal to one if and only if
$t$ is the earliest period for sample path $i$ that satisfies the equality
$b^{i}_{t}=1$. With this observation, the first step in our derivation of (BP)
is comprised of the following lemma:

###### Lemma 4.4

The optimal objective value of (IP) is equal to the optimal objective value of
(BP-1). Furthermore, if $b$ is a feasible solution for (BP-1) and if
$\sigma^{i}\triangleq\min\\{\min\\{t\in\\{1,\ldots,T\\}:b^{i}_{t}=1\\},T\\}$
for each $i\in\\{1,\ldots,N\\}$, then

$\displaystyle\frac{1}{N}\sum_{i=1}^{N}\min_{t\in\\{1,\ldots,\sigma^{i}\\}}\left\\{g(t,x^{i}):\textnormal{
there exists }j\textnormal{ such that
}\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\textnormal{ and
}\sigma^{j}=t\right\\}\geq\frac{1}{N}\sum_{i=1}^{N}\psi^{i}(b).$

The above lemma establishes the equivalence of (IP) of (BP-1), and it shows
that any feasible solution for (BP-1) can be transformed into a feasible
solution for (IP) with the same or greater objective value. In the second and
final step in our derivation of (BP), denoted below by Lemma 4.5, we show that
each function $\psi^{i}(b)$ can be represented as the optimal objective
function of a linear optimization problem.

###### Lemma 4.5

For each sample path $i\in\\{1,\ldots,N\\}$, $\psi^{i}(b)$ is equal to the
optimal objective value of the following linear optimization problem:

$\displaystyle\underset{w^{i}}{\textnormal{maximize}}$
$\displaystyle\sum_{t=1}^{T}\sum_{\ell=1}^{L^{i}_{t}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})b^{i}_{t}(1-w^{i}_{t\ell})$
(3) subject to $\displaystyle\begin{aligned} &w^{i}_{t\ell}\leq
w^{i}_{t+1,\ell}&&\textnormal{for all
}t\in\\{1,\ldots,T-1\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}\\\
&w^{i}_{t\ell}\leq w^{i}_{t,\ell+1}&&\textnormal{for all
}t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|-1\\}\\\
&b^{j}_{t}\leq w^{i}_{t\ell}&&\textnormal{for all
}j\in\\{1,\ldots,N\\}\textnormal{ and }t\in\\{1,\dots,T\\}\\\
&&&\textnormal{such that }g(t,x^{i})=\kappa^{i}_{\ell}\textnormal{ and
}\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\\\ &b^{i}_{t}\leq
w^{i}_{t+1,1}&&\textnormal{for all }\;t\in\\{1,\ldots,T-1\\}\\\
&w^{i}_{t\ell}\in\mathbb{R}&&\textnormal{for all
}\;t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}.\end{aligned}$

Lemma 4.5 thus shows that the purpose of the decision variables $w^{i}_{s}$
and the constraints of (BP) is to encode the value of $\psi^{i}(b)$ for each
sample path $i\in\\{1,\ldots,N\\}$. Finally, we readily observe that the
combination of Lemmas 4.4 and 4.5 immediately yields the proof of Theorem 4.3.

###  4.2 Heuristic Algorithm

In §4.1, we developed an exact reformulation of the robust optimization
problem (RO) as a zero-one bilinear program over totally unimodular
constraints (BP). That bilinear program can be solved by off-the-shelf
software for mixed-integer linear optimization; see Appendix 12. Thus, the
previous subsection can be viewed as a concrete and easily implementable exact
algorithm for solving the robust optimization problem.

Building upon (BP), we now turn to the task of developing efficient algorithms
that can approximately solve the robust optimization problem. Specifically,
the main contribution of §4.2 is a practical heuristic for solving the robust
optimization problem with computation time that is polynomial in both the
number of sample paths $N$ and the number of periods $T$. The proposed
heuristic is thus significantly more computational tractable than solving
(BP), which was shown in Theorem 4.2 to be NP-hard. Moreover, we will provide
theoretical and empirical evidence that the proposed heuristic can perform
surprisingly well, both with respect to approximation quality and
computational tractability. All omitted proofs of results from §4.2 can be
found in Appendix 14.

This subsection is organized as follows. In §4.2.1, we discuss the high-level
motivation and intuition behind our proposed heuristic for the robust
optimization problem. In §4.2.2 and §4.2.3, we formalize the heuristic and
offer a strongly polynomial-time algorithm for implementing it. In §4.2.4, we
provide theoretical guarantees which show that the approximation gap of the
heuristic cannot be arbitrarily bad and, in some cases, is guaranteed to be
equal to zero. We study the empirical performance of the heuristic in the
subsequent §5.

####  4.2.1 Preliminaries.

Our heuristic for approximately solving the robust optimization problem is
motivated by the structure of the exact reformulation (BP) from §4.1 for the
robust optimization problem. Recall that (BP) is a zero-one optimization
problem with a bilinear objective function and totally unimodular constraints.
Because the constraints of (BP) are totally unimodular, the computational
difficulty of solving (BP) can thus be attributed to the nonlinearity of its
objective function. The proposed heuristic, which is formalized in §4.2.2,
aims to contend with this difficulty by approximating the nonlinear objective
function of (BP) with a linear function.

To describe the motivation behind our heuristic in greater detail, consider
any function $f(b,w)$ that is linear in the vectors $b,w$ and is a lower bound
on the objective function of the optimization problem (BP) for all $b,w$ that
satisfy the constraints of (BP). That is, let $f(b,w)$ be a linear function
that satisfies the following condition:

$\displaystyle
f(b,w)\leq\frac{1}{N}\sum_{i=1}^{N}\sum_{t=1}^{T}\sum_{\ell=1}^{L^{i}_{t}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})b^{i}_{t}(1-w^{i}_{t\ell})\;\;\;\;\textnormal{for
all }b,w\textnormal{ that are feasible for \eqref{prob:bp}.}$ (4)

Then it follows immediately from the condition on line (4) that a
conservative, lower-bound approximation of the optimization problem (BP) is
given by the following optimization problem:

$\displaystyle\underset{b,w}{\textnormal{maximize}}$ $\displaystyle f(b,w)$
(H) subject to $\displaystyle\begin{aligned} &w^{i}_{t\ell}\leq
w^{i}_{t+1,\ell}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T-1\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}\\\
&w^{i}_{t\ell}\leq w^{i}_{t,\ell+1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|-1\\}\\\
&b^{i}_{t}\leq w^{i}_{t+1,1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T-1\\}\\\ &b^{j}_{t}\leq
w^{i}_{t\ell}&&\textnormal{for all }i,j\in\\{1,\ldots,N\\}\textnormal{ and
}t\in\\{1,\dots,T\\}\textnormal{ such that }g(t,x^{i})=\kappa^{i}_{\ell}\\\
&&&\textnormal{and }\mathcal{U}^{i}_{t}\cap\mathcal{U}^{j}_{t}\neq\emptyset\\\
&b^{i}_{t}\in\\{0,1\\}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\dots,T\\}\\\
&w^{i}_{t\ell}\in\mathbb{R}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\\{1,\ldots,T\\},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}.\end{aligned}$

Indeed, we observe that the constraints of (H) are identical to the
constraints of (BP). The only difference between these two optimization
problems is that the nonlinear objective function of (BP) has been replaced
with the linear function $f(b,w)$ that satisfies the condition from line (4).
Therefore, we conclude that every optimal solution for the optimization
problem (H) will be a feasible solution for the optimization problem (BP), and
the optimal objective value of (H) will always be less than or equal to the
optimal objective value of (BP).

Given any linear function $f(b,w)$ that satisfies the condition from line (4),
the above optimization problem (H) is useful because it provides a means to
computing approximate solutions for (BP). Indeed, the fact that (H) and (BP)
have the same constraints implies that any optimal solution for (H) is a
feasible solution for (BP). Moreover, it is always theoretically possible to
choose the function $f(b,w)$ such that every optimal solution for (H) is an
optimal solution for (BP), as shown by the following Proposition 4.6. We thus
conclude that solving the optimization problem (H) has the potential to yield
high-quality approximate solutions for (BP).

###### Proposition 4.6

For any instance of (BP), there exists a linear function $f(b,w)$ satisfying
the condition from line (4) such that the optimal objective value of (H) is
equal to the optimal objective value of (BP), and every optimal solution for
(H) is also an optimal solution for (BP).

The optimization problem (H) is ultimately attractive compared to (BP) from
the perspective of computational tractability. Indeed, we observe that (H) is
a zero-one linear optimization problem over totally unimodular constraints.
This implies that the integrality constraints on the decision variables can be
relaxed without loss of generality; that is, (H) can be solved as a linear
optimization problem in which each constraints $b^{i}_{t}\in\\{0,1\\}$ are
replaced with $0\leq b^{i}_{t}\leq 1$ (Conforti et al. 2014, §4.2). The
optimization problem (H) is thus particularly convenient from an
implementation standpoint: as a linear optimization problem, (H) can be easily
formulated and solved directly by commercial linear optimization solvers such
as CPLEX or Gurobi.

In fact, the optimization problem (H) can be solved _very_ efficiently by
exploiting the relationship between (H) and the maximal closure problem. The
maximal closure problem is a problem from combinatorial optimization that has
been widely studied in the operations research literature dating back to Rhys
(1970) and Picard (1976), with applications ranging from project scheduling to
open-pit mining.666For further background on the maximal closure problem, we
refer the interested reader to Hochbaum (2004). The goal in the maximal
closure problem is to find a maximum-weight closure in a directed graph, where
a closure is defined as a subset of vertices without edges that leave the
subset. In particular, it follows immediately from Picard (1976, §3) that the
linear optimization problem (H) with $\mathcal{O}(NT^{2})$ decision variables
and $\mathcal{O}(NT^{2})$ constraints is equivalent to computing the maximal
closure of a directed graph with $\mathcal{O}(NT^{2})$ nodes and
$\mathcal{O}(NT(T+N))$ edges. An important algorithmic property is that
maximal closure problems can be solved by computing the maximum flow in the
directed graph (Picard 1976, §4). Thus, as we will formalize in §4.2.3, the
optimization problem (H) can be solved with running time that is strongly
polynomial with respect to both the number of simulated sample paths $N$ and
the number of time periods $T$.777We note that Proposition 4.6 does not
contradict our hardness result for the robust optimization problem (RO) from
Theorem 4.2, as the right choice of the linear function $f(b,w)$ in
Proposition 4.6 may be difficult to identify.

####  4.2.2 Description of Heuristic.

In §4.2.1, we discussed the high-level motivation and intuition behind our
proposed heuristic. In view of that motivation, we now formally describe our
proposed heuristic for approximately solving the optimization problem (BP).
Specifically, our proposed heuristic consists of solving the optimization
problem (H) with a particular linear objective function, denoted below by
$\bar{f}(b,w)$.

To formally describe our proposed heuristic, we require the following
notation. Recall that $T$ is the number of periods in the optimal stopping
problem. For each sample path $i$, let
$T^{i}\triangleq\argmax_{t\in\\{1,\ldots,T\\}}\;{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g(t,x^{i})}$
be defined as the period in which the sample path $i$ achieves its maximum
reward, and if there are multiple optimal solutions, we choose the optimal
solution that is smallest. To simplify our notation, let us define
$\mathcal{T}^{i}\triangleq\\{T^{i}\\}\cup\\{T\\}$ as the set that contains the
period in which sample path $i$ achieves its maximum reward as well as the
last period of the optimal stopping period.

With this additional notation, the proposed heuristic obtains approximate
solutions $b,w$ for the optimization problem (BP) by solving the optimization
problem (H) with the objective function

$\displaystyle\bar{f}(b,w)\triangleq\frac{1}{N}\sum_{i=1}^{N}\sum_{t\in\mathcal{T}^{i}}\sum_{\ell=1}^{L^{i}_{t}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})\left(b^{i}_{t}-w^{i}_{t\ell}\right).$

In other words, our heuristic for solving the optimization problem (BP)
consists of solving the optimization problem (H) with objective function given
by $f(b,w)=\bar{f}(b,w)$. The optimal solution for (H) with objective function
$\bar{f}(b,w)$ constitutes a feasible solution for (BP), since the constraints
of (H) and (BP) are identical. This solution can thus be transformed into a
Markovian stopping rule using the transformations described in Theorem 3.7 and
4.3.

###### Remark 4.7

We readily observe that the function $\bar{f}(b,w)$ is linear in $b$ and $w$.
Moreover, it follows from algebra that $\bar{f}(b,w)$ is less than or equal to
the objective function of (BP) for all feasible solutions of (BP); that is,
the linear function $\bar{f}(b,w)$ satisfies the condition from line (4).888To
see why $\bar{f}(b,w)$ is a lower bound on the objective function of (BP),
consider any arbitrary vectors $b,w$ that satisfy the constraints of (BP).
Since feasibility for the optimization problem (BP) implies that $b$ is a
binary vector, we observe that the equality
$b^{i}_{t}(1-w^{i}_{t\ell})=\max\left\\{b^{i}_{t}-w^{i}_{t\ell},0\right\\}$
holds for each $i\in\\{1,\ldots,N\\}$, $t\in\\{1,\ldots,T\\}$, and
$\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}$. Therefore, a lower-bound
approximation of the objective function of the optimization problem (BP) can
be obtained by replacing each term $b^{i}_{t}(1-w^{i}_{t\ell})$ with
$b^{i}_{t}-w^{i}_{t\ell}$ if $t\in\mathcal{T}^{i}$ and with $0$ if
$t\notin\mathcal{T}^{i}$. Thus, we conclude that solving the optimization
problem (H) with objective function $f(b,w)=\bar{f}(b,w)$ will provide a
lower-bound approximation of the optimization problem (BP), and any optimal
solution for (H) will be a feasible solution for (BP).

Our motivations for using the heuristic outlined above are two-fold. First, we
find that the optimization problem (H) with the objective function
$\bar{f}(b,w)$ is highly tractable from both a theoretical and empirical
standpoint on realistic sizes of robust optimization problems. Second, we
provide theoretical and empirical evidence that the above heuristic can find
high-quality and in some cases optimal solutions for (BP). We elaborate on the
theoretical aspects of these two motivations in §4.2.3 and §4.2.4, and we
explore the empirical performance of the heuristic in the subsequent §5.

####  4.2.3 Algorithms for Implementing Heuristic.

In §4.2.1, we argued generally that the optimization problem (H) with any
linear objective function $f(b,w)$ can be solved in strongly polynomial time
by reducing (H) to a maximal closure problem. We now formalize those earlier
arguments to derive an explicit, strongly polynomial-time algorithm for
solving the optimization problem (H) in the particular case where
$f(b,w)=\bar{f}(b,w)$. Specifically, the main contribution of §4.2.3 is the
development of an algorithm that achieves the running time that is specified
in the following proposition:

###### Proposition 4.8

If $f(b,w)=\bar{f}(b,w)$, then (H) can be solved in $\mathcal{O}(N^{2}T(N+T))$
time.

The above proposition establishes the computational tractability of our
heuristic by demonstrating that an optimal solution for the optimization
problem (H) with objective function $f(b,w)=\bar{f}(b,w)$ can be computed with
running time that is polynomial in both the number of simulated sample paths
as well as the number of time periods. In particular, if the number of periods
$T$ is held constant, then the running time in Proposition 4.8 scales
cubically in the number of simulated sample paths $N$. Such a tractability
guarantee is ultimately important from a practical perspective: indeed, in the
following §5, we will present numerical experiments which show that our
heuristic can run in seconds on realistic problem sizes with over fifty
periods and thousands of sample paths. An outline of the proof of Proposition
4.8 is found throughout the rest of §4.2.3.

To develop an algorithm with the computation time that is specified in
Proposition 4.8, we begin by deriving a compact reformulation of (H) when the
objective function satisfies $f(b,w)=\bar{f}(b,w)$. The compact reformulation
is denoted below as the optimization problem ($\bar{\textnormal{H}}$). To
derive this compact reformulation, we first show that many of the decision
variables in the optimization problem (H) can, without loss of generality, be
removed from (H) when the objective function satisfies $f(b,w)=\bar{f}(b,w)$.
In particular, our compact reformulation will utilize the following lemma:

###### Lemma 4.9

If $f(b,w)=\bar{f}(b,w)$, then there exists an optimal solution $b,w$ for (H)
that satisfies $b^{1}_{T}=\cdots=b^{N}_{T}=1$ and satisfies $b^{i}_{t}=0$ for
each sample path $i\in\\{1,\ldots,N\\}$ and period
$t\in\\{1,\ldots,T\\}\setminus\mathcal{T}^{i}$.

The above lemma demonstrates that the optimal values for many of the decision
variables of (H) can be known in advance when the objective function satisfies
$f(b,w)=\bar{f}(b,w)$. Hence, Lemma 4.9 implies that each decision variable
$b^{i}_{t}$ can be fixed to its optimal value and removed from (H) when
$f(b,w)=\bar{f}(b,w)$, except for the decision variables $b^{i}_{t}$ in which
the sample path $i\in\\{1,\ldots,N\\}$ and period $t\in\\{1,\ldots,T\\}$
satisfy $t=T^{i}$ and $T^{i}<T$. In view of Lemma 4.9, we now consider the
following optimization problem:

$\displaystyle\underset{b,w}{\textnormal{maximize}}$
$\displaystyle\frac{1}{N}\sum_{i\in\\{1,\ldots,N\\}:T^{i}<T}\sum_{\ell=1}^{L^{i}_{T^{i}}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})\left(b^{i}_{T^{i}}-w^{i}_{T^{i}\ell}\right)+\frac{1}{N}\sum_{i=1}^{N}\sum_{\ell=1}^{L^{i}_{T}-1}(\kappa^{i}_{\ell+1}-\kappa^{i}_{\ell})\left(1-w^{i}_{T\ell}\right)$
($\bar{\textnormal{H}}$) subject to $\displaystyle\begin{aligned}
&w^{i}_{t\ell}\leq w^{i}_{t,\ell+1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\mathcal{T}^{i},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|-1\\}\\\
&b^{i}_{T^{i}}\leq w^{i}_{T1}&&\textnormal{for all
}i\in\\{1,\ldots,N\\}\textnormal{ such that }T^{i}<T\\\ &b^{j}_{T^{j}}\leq
w^{i}_{t\ell}&&\textnormal{for all }i,j\in\\{1,\ldots,N\\}\textnormal{ and
}t\in\mathcal{T}^{i}\textnormal{ such that
}g(T^{j},x^{i})=\kappa^{i}_{\ell},\\\ &&&T^{j}\leq
t,\;T^{j}<T,\;\textnormal{and
}\mathcal{U}^{i}_{T^{j}}\cap\mathcal{U}^{j}_{T^{j}}\neq\emptyset\\\
&b^{i}_{T^{i}}\in\\{0,1\\}&&\textnormal{for all
}i\in\\{1,\ldots,N\\}\textnormal{ such that }T^{i}<T\\\
&w^{i}_{t\ell}\in\mathbb{R}&&\textnormal{for all
}i\in\\{1,\ldots,N\\},\;t\in\mathcal{T}^{i},\;\ell\in\\{1,\ldots,|\mathcal{K}^{i}|\\}.\end{aligned}$

It can be shown using straightforward algebra that the above optimization
problem ($\bar{\textnormal{H}}$) is equivalent to the optimization problem (H)
when the objective function of (H) satisfies $f(b,w)=\bar{f}(b,w)$ and when
the decision variables in (H) are restricted to satisfy
$b^{1}_{T}=\cdots=b^{N}_{T}=1$ and $b^{i}_{t}=0$ for all
$i\in\\{1,\ldots,N\\}$ and $t\in\\{1,\ldots,T\\}\setminus\mathcal{T}^{i}$; see
proof of the following Lemma 4.10 for details. Hence, it follows from Lemma
4.9 that any optimal solution for ($\bar{\textnormal{H}}$) can be transformed
into an optimal solution for (H) when the objective function satisfies
$f(b,w)=\bar{f}(b,w)$. We formalize the equivalence of (H) and
($\bar{\textnormal{H}}$) in the following Lemma 4.10:

###### Lemma 4.10

Let $f(b,w)=\bar{f}(b,w)$. Then, the optimal objective value of
($\bar{\textnormal{H}}$) is equal to the optimal objective value of (H).
Moreover, let $\bar{b},\bar{w}$ denote an optimal solution for
($\bar{\textnormal{H}}$). Then there exists an optimal solution $b,w$ for (H)
that satisfies the following equality for each sample path
$i\in\\{1,\ldots,N\\}$ and period $t\in\\{1,\ldots,T\\}$:

$\displaystyle{b}^{i}_{t}$
$\displaystyle=\begin{cases}\bar{b}^{i}_{t},&\textnormal{if
}t=T^{i}\textnormal{ and }T^{i}<T,\\\ 1,&\textnormal{if }t=T,\\\
0,&\textnormal{otherwise}.\end{cases}$

To summarize, we have shown in Lemmas 4.9 and 4.10 that solving the
optimization problem (H) can be reduced to solving the optimization problem
($\bar{\textnormal{H}}$) when the objective function of (H) satisfies
$f(b,w)=\bar{f}(b,w)$. More precisely, we have shown that any optimal solution
for ($\bar{\textnormal{H}}$) can be transformed into an optimal choice for the
decision variables ${b}$ in (H) when the objective function of (H) satisfies
$f(b,w)=\bar{f}(b,w)$. By applying the transformations described in Theorems
3.7 and 4.3, we can construct exercise policies
$\mu\equiv(\mu_{1},\ldots,\mu_{T})$ from the decision variables $b$ whose
robust objective value in the robust optimization problem
$\widehat{J}_{N,\epsilon}(\mu)$ is greater than or equal to the optimal
objective value of ($\bar{\textnormal{H}}$).

From a computational tractability standpoint, the optimization problem
($\bar{\textnormal{H}}$) is ultimately attractive because it has significantly
fewer decision variables and fewer constraints than the optimization problem
(H). Indeed, we observe from inspection that ($\bar{\textnormal{H}}$) has
$\mathcal{O}(NT)$ decision variables and $\mathcal{O}(NT+N^{2})$ constraints,
whereas (H) has $\mathcal{O}(NT^{2})$ decision variables and
$\mathcal{O}(NT^{2}+N^{2}T)$ constraints. Because of the same reasoning as
given in §4.2.1, the integrality constraints on the decision variables can
also be relaxed without loss of generality; that is, (H) can be solved as a
linear optimization problem in which each constraints $b^{i}_{t}\in\\{0,1\\}$
are replaced with $0\leq b^{i}_{t}\leq 1$. Thus, ($\bar{\textnormal{H}}$) can
be easily formulated and solved directly by commercial linear optimization
solvers.

We conclude §4.2.3 by using the above optimization problem
($\bar{\textnormal{H}}$) to outline the proof of Proposition 4.8.
Specifically, our proof of Proposition 4.8 consists of reformulating
($\bar{\textnormal{H}}$) as a maximal closure problem in a directed graph. We
then use the efficient maximum-flow algorithm of Orlin (2013) to solve the
maximal closure problem. The remaining details for the proof of Proposition
4.8 are found in Appendix 14.

####  4.2.4 Approximation Guarantees.

We conclude §4.2 by performing a theoretical analysis of the approximation
quality of the proposed heuristic ($\bar{\textnormal{H}}$) for the
optimization problem (BP). Our motivation here is to understand whether the
proposed heuristic is ever guaranteed to find optimal solutions for (BP) and,
conversely, whether it is ever possible for the gap between the optimal
objective values of ($\bar{\textnormal{H}}$) and (BP) to be arbitrarily large.
Our answers to these theoretical questions are presented below in Propositions
4.11, 4.12, and 4.13.

We begin by focusing on optimal stopping problems in which the number of
periods is equal to two. We recall from Theorem 4.1 in §4.1 that any instance
of the robust optimization problem with two periods can be solved exactly in
strongly polynomial-time. We will now prove Theorem 4.1 by establishing that
the approximation gap between ($\bar{\textnormal{H}}$) and (BP) is always
equal to zero when the number of periods is equal to two:

###### Proposition 4.11

If $T=2$, then the optimal objective values of ($\bar{\textnormal{H}}$) and
(BP) are equal.

Our takeaways from the above proposition are two-fold. First, Proposition 4.11
shows that there indeed exist settings in which the proposed heuristic is
guaranteed to find optimal solutions for (BP). Second, Proposition 4.11 in
combination with Theorems 3.7 and 4.3 implies that solving the robust
optimization problem (RO) can reduced to solving the optimization problem
($\bar{\textnormal{H}}$) when $T=2$. Therefore, we observe that the proof of
Theorem 4.1 follows immediately from Propositions 4.8 and 4.11.

We next develop a general bound on the gap between the optimal objective
values of ($\bar{\textnormal{H}}$) and (BP) that holds for any fixed number of
periods $T$. For notational convenience, we let
$J^{\textnormal{\ref{prob:bp}}}$ denote the optimal objective value of (BP)
and $J^{\textnormal{\ref{prob:h_bar}}}$ denote the optimal objective value of
($\bar{\textnormal{H}}$).

###### Proposition 4.12

$J^{\textnormal{\ref{prob:h_bar}}}\geq\frac{1}{T}J^{\textnormal{\ref{prob:bp}}}$.

The above proposition can be viewed as valuable from a theoretical
perspective, as it shows that the gap between the optimal objective values of
($\bar{\textnormal{H}}$) and (BP) can never be arbitrarily large. Indeed, we
recall from Remark 4.7 and Lemma 4.10 that the optimal objective value of
($\bar{\textnormal{H}}$) is always less than or equal to the optimal objective
value of (BP). Thus, Proposition 4.12 establishes that there exist bounds on
the gap between the optimal objective values of ($\bar{\textnormal{H}}$) and
(BP) that are independent of the choice of the number of simulated sample
paths $N$, the choice of the robustness parameter $\epsilon$ that controls the
size of the uncertainty sets, the reward function, etc.

Our third theoretical guarantee in §4.2.4, stated below as Proposition 4.13,
provides a bound on the gap between the optimal objective values of
($\bar{\textnormal{H}}$) and (BP) that does not depend on the number of
periods $T$. Our theoretical bound in Proposition 4.13 will require that the
reward function of the NMOS problem satisfies the following regularity
condition:

There exists a constant $L>0$ such that, for each period
$t\in\\{1,\ldots,T\\}$, the reward function satisfies
$|g(t,y)-g(t,y^{\prime})|\leq L\|y_{t}-y_{t}^{\prime}\|_{\infty}$ for all
$y\equiv(y_{1},\ldots,y_{T}),y^{\prime}\equiv(y^{\prime}_{1},\ldots,y^{\prime}_{T})\in\mathcal{X}^{T}$.
The above assumption stipulates for each period $t\in\\{1,\ldots,T\\}$ that
the reward function $y\mapsto g(t,y)$ is Lipschitz-continuous with respect to
the state $y_{t}\in\mathcal{X}$. The above assumption is relatively mild in
applications such as options pricing when the reward function depends only on
the current state in each period.999For example, we observe that Assumption
4.2.4 is satisfied by the reward function in the multi-dimensional barrier
option pricing problem from §5.2 when the state space in the optimal stopping
problem is augmented as
$(x_{t},q_{t})\equiv(\max_{a\in\\{1,\ldots,d\\}}\xi_{t,a},1-\mathbb{I}\\{\max_{a\in\\{1,\ldots,d\\}}\xi_{s,a}\leq
B(s)\text{ for all }s\in\\{1,\ldots,t\\}\\})$.  We will impose Assumption
4.2.4 in the following Proposition 4.13 to eliminate pathological situations
in which slight perturbations of sample paths lead to drastic changes in
reward. Once again, for notational convenience, we let
$J^{\textnormal{\ref{prob:h_bar}}}$ and $J^{\textnormal{\ref{prob:bp}}}$
denote the optimal objective values of the optimization problems
($\bar{\textnormal{H}}$) and (BP), respectively.

###### Proposition 4.13

If Assumption 4.2.4 holds, then
$J^{\textnormal{\ref{prob:h_bar}}}\geq\frac{1}{\log
N+1}J^{\textnormal{\ref{prob:bp}}}-2\epsilon L$.

Compared to the bound from Proposition 4.12, we observe that the bound from
Proposition 4.13 is attractive when the robustness parameter $\epsilon$ is
relatively small and when the number of simulated sample paths $N$ is
relatively small compared to the number of periods $T$. For example, we
observe that $\frac{1}{\log N+1}\ll\frac{1}{T}$ will hold whenever the number
of sample paths $N$ is subexponential in the number of periods $T$.

In summary, Propositions 4.11, 4.12, and 4.13 establish under mild and
verifiable conditions that the approximation gap between the proposed
heuristic ($\bar{\textnormal{H}}$) and the robust optimization problem (RO)
cannot be arbitrarily large and, in some cases, is guaranteed to be equal to
zero. Ultimately, the practical value of the proposed heuristic
($\bar{\textnormal{H}}$) lies in its performance in the context of optimal
stopping applications. In the following §5, we provide numerical evidence that
($\bar{\textnormal{H}}$) can indeed find high-quality stopping rules for
realistic stochastic optimal stopping problems in practical computation times.

##  5 Numerical Experiments

In this section, we perform numerical experiments to compare our robust
optimization approach and three state-of-the-art benchmarks from the
literature (Longstaff and Schwartz 2001, Ciocan and Mišić 2020, Desai et al.
2012). The first two benchmarks serve as representatives of two classes of
approximation methods for stochastic optimal stopping problems (approximate
dynamic programming and parametric exercise policies) which, similarly as the
robust optimization approach, only require the ability to simulate sample
paths of the entire sequence of random states. The third benchmark is
representative of state-of-the-art duality-based methods for stochastic
optimal stopping problems. All experiments were conducted on a 2.6 GHz 6-Core
Intel Core i7 processor with 16 GB of memory. All methods are implemented in
the Julia programming language and solved using the JuMP library and Gurobi
optimization software.

###  5.1 A Simple Non-Markovian Problem

To demonstrate the value of our robust optimization approach, we begin by
investigating a simple, one-dimensional stochastic optimal stopping problem
with a non-Markovian probability distribution. The optimal stopping problem of
consideration involves a state space which is equal to the real numbers
($\mathcal{X}=\mathbb{R}^{1}$) and a reward function of $g(t,x)=x_{t}$ in each
period $t\in\\{1,\ldots,T\\}$. For any fixed duration $\Delta\in\mathbb{N}$,
the joint probability distribution of the stochastic process is given by

$\displaystyle x_{t}$
$\displaystyle\sim\text{Uniform}[0,1]+\frac{2\theta}{T}\mathbb{I}\left\\{\theta\leq
t\leq\theta+\Delta\right\\}\quad\text{for all }t=1,\ldots,T,$

where the random parameter $\theta\sim\text{Uniform}\\{1,2,\ldots,T-\Delta\\}$
is selected once per sample path and is unobserved. Simulated sample paths of
this non-Markovian stochastic process $x=(x_{1},\ldots,x_{T})$ are visualized
in Figure LABEL:fig:nonmarkov_trajectories. We perform numerical experiments
on the following methods:

  * •

Robust Optimization (RO): The robust optimization approach is used here to
approximate (OPT) over the non-Markovian stochastic process
$x\equiv(x_{1},\ldots,x_{T})$, and thus aims to find the best Markovian
stopping rules for the stochastic optimal stopping problem. The method was run
with robustness parameters $\epsilon\in\\{0,0.01,\ldots,0.2\\}$ for $N<1000$
and $\epsilon\in\\{0,0.01,\ldots,0.1\\}$ otherwise, and the robust
optimization problem is solved approximately using our proposed heuristic
($\bar{\textnormal{H}}$) from §4.2.3.

  * •

Least-Squares Regression (LS): We implement the method of Longstaff and
Schwartz, which employs least-squares regression to approximate the
continuation value function (i.e., the expected reward from not stopping) in
each period using backwards recursion. To apply this method, we first
transform each non-Markovian sample path
$x^{i}\equiv(x^{i}_{1},\ldots,x^{i}_{T})$ into an augmented Markovian sample
path of the form $X^{i}\equiv(X^{i}_{1},\ldots,X^{i}_{T})$ by adding the full
history into the state in each period:
$X^{i}_{t}\triangleq(x^{i}_{1},\ldots,x^{i}_{t},0,\ldots,0)$. The regression
step requires a specification of basis functions, and we consider the
following categories of basis functions.

##### Full-History:

This category of basis functions uses the entire vector $X^{i}_{t}$ of states
observed up to that point. The basis functions that we consider in this
category are One (the constant function 1), Prices (the states observed up to
that point, $X^{i}_{t}\in\mathbb{R}^{T}$), and Prices2 (the product of each
pair of states observed up to that point,
$X^{i}_{t}(X^{i}_{t})^{\intercal}\in\mathbb{R}^{T\times T}$).

##### Markovian:

This category of basis functions uses only the current state $x^{i}_{t}$ in
each period. The purpose of considering these basis functions is to analyze
the performance when, like the robust optimization approach, the method of
Longstaff-Schwartz is restricted to stopping rules which depend only on the
current state in each period. We consider basis functions based on the
Laguerre polynomials, where Laguerre-k is the polynomial
$\sum_{\ell=0}^{k}\binom{k}{\ell}\frac{(-1)^{\ell}}{\ell!}(x^{i}_{t})^{\ell}$.

  * •

Tree Method (Tree): The method of Ciocan and Mišić approximates the stochastic
optimal stopping problem (OPT) by restricting the space of exercise policies
to decision trees. Like the robust optimization approach, the tree approach is
used to find Markovian stopping rules for the non-Markovian optimal stopping
problem (OPT). We apply the method with the same information as the robust
optimization approach at each time period (the current state and the time
period) and with a splitting parameter of 0.005.

In our computational experiments, we consider a time horizon of $T=50$
stopping periods with a duration parameter of $\Delta=5$. All methods were run
using simulated training datasets of sizes
$N\in\\{10^{2},10^{2.1},\ldots,10^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}3.9}},10^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}4}}\\}$.
The robustness parameters in the “RO” method were selected using a validation
set of size $\bar{N}=10^{3}$; see §2.5 for more details. All methods were
evaluated on a common and independent testing dataset of $\tilde{N}=10^{5}$
sample paths, and experiments were repeated over 10 replications.

Figures LABEL:fig:nonmarkov_trajectories and LABEL:fig:nonmarkov_results
visualize the Markovian stopping rules found by the robust optimization
approach, the expected rewards of the stopping rules obtained by the various
methods, and the computation times of the various methods. The results of this
experiment show that the robust optimization approach outperforms the other
approaches, producing stopping rules with an expected reward of approximately
1.62 from training datasets of $N=10^{3}$ simulated sample paths and with a
computation time of approximately 0.713 seconds. We reflect below on the main
differences between our robust optimization approach and the benchmark
methods.

For the least-squares regression method with full-history, selecting a good
choice of basis functions is found to be a first-order challenge. For example,
the basis functions (One, Prices) turn out to be insufficiently rich to
provide an accurate approximation of the continuation function with the full
state history. When the number of sample paths is sufficiently large, we
expect that the basis functions (One, Prices, Prices2) should provide a better
approximation of the stochastic optimal stopping problem. However, the
computational cost resulting from this rich class of basis functions precluded
its practicality on sufficiently large training datasets. In contrast, the
robust optimization approach only searches for Markovian stopping rules, and
in doing so, has a reduction in sample complexity. In addition to achieving a
significantly better expected reward in this example, the Markovian stopping
rules found by the robust optimization approach are considerably more
interpretable than those which have full-history dependance, as illustrated in
Figure LABEL:fig:nonmarkov_trajectories.

Compared to the tree method, the complex structure of the Markovian stopping
rules found by the robust optimization approach demonstrates the value of
algorithms that do not impose parametric restrictions on the exercise
policies. In theory, decision trees with sufficient depth are capable of
approximating the best Markovian stopping rules to the NMOS problem to
arbitrary accuracy (Ciocan and Mišić 2020, Theorem 2). However, the greedy
heuristic that is proposed by Ciocan and Mišić (2020) to efficiently optimize
over decision trees is not able to find a good approximation of the optimal
Markovian stopping rule in this example, even as the number of sample paths is
large. We note that this issue is not unique to the decision trees _per se_ ;
the computational intractability of optimizing over parametric spaces of
Markovian stopping rules has resulted in heuristics for other settings as well
(Glasserman 2013, §8.2). This example provides evidence that our robust
optimization approach, in conjunction with the proposed approximation from
§4.2, can yield high-quality Markovian stopping rules to complex (low-
dimensional) non-Markovian optimal stopping problems.

Like the robust optimization and tree methods, “LS, Markovian” aims to find
Markovian stopping rules for the non-Markovian optimal stopping problem. In
particular, we observe that “LS, Markovian” was implemented with basis
functions of Laguerre polynomials with a rather large degree.101010We were
unable to run the “LS, Markovian (Laguerre 0-K)” for degrees $K>15$ due to the
numerical precision required to encode the coefficients in the basis
functions. The expected reward of “LS, Markovian (Laguerre 0-15)” may thus be
interpreted as an estimate for the best Markovian stopping rule that one could
achieve using dynamic programming. However, the optimality of backwards
induction no longer holds for exercise policies that depend on the current
state when the stochastic process is non-Markovian. This explains the inferior
performance of this method in this example and motivates the use of methods
for addressing non-Markovian optimal stopping problems that optimize directly
over the exercise policies in all periods simultaneously. A further discussion
on the limitations of dynamic programming in finding Markovian stopping rules
for NMOS problems can be found in Appendix 7.

In Figure LABEL:fig:nonmarkov_best_epsilon, we present the relationship
between the best choice of the robustness parameter for our heuristic
($\bar{\textnormal{H}}$) and the number of simulated sample paths $N$. The
results in this figure show that the best choice of the robustness parameter
$\epsilon$ decreases as the number of simulated sample paths $N$ increases.
Although a formal analysis of the relationship between these two quantities is
outside the scope of this paper, Figure LABEL:fig:nonmarkov_best_epsilon
provides empirical evidence that only smaller choices of the robustness
parameter need to be considered when using larger numbers of simulated sample
paths.

In Figure LABEL:fig:nonmarkov_epsilon, we compare the (in-sample) robust
objective values and the expected rewards of stopping rules obtained by our
heuristic ($\bar{\textnormal{H}}$) and our exact reformulation (BP) of the
robust optimization problem. The comparison in Figure
LABEL:fig:nonmarkov_epsilon is performed on a smaller instance where $T=20$,
$\Delta=2$, and the robust optimization problem is constructed from a training
set of $N=100$ simulated sample paths. We observe the gap in objective values
between the two algorithms is relatively small across choices of the
robustness parameter, and the gap is equal to zero when the robustness
parameter is set equal to zero. Moreover, the expected rewards of stopping
rules obtained from these two algorithms were similarly relatively close
across choices of the robustness parameter. We view these results as promising
as they suggest, at least for the present example, that the significantly more
tractable optimization problem ($\bar{\textnormal{H}}$) can provide a close
approximation of the robust optimization problem (RO).

In summary, the experiments from this subsection reveal our first setting in
which the robust optimization approach is attractive. In non-Markovian optimal
stopping problems, there can exist high-quality stopping rules which are
Markovian. However, finding such stopping rules is a non-trivial challenge, as
the best Markovian stopping rules will not necessarily exhibit a simple
structure, and cannot in general be found using dynamic programming. The
robust optimization approach thus provides a practical nonparametric approach
for finding them. In view of these observations, the next subsection considers
a benchmark optimal stopping problem from the options pricing literature.

###  5.2 Pricing Multi-Dimensional Barrier Options

Building upon the previous section, we next consider the well-studied problem
of pricing discretely-monitored barrier call options over multiple assets.
These are stochastic optimal stopping problems in which the random state at
each period, $\xi_{t}\in\Xi\equiv\mathbb{R}^{d}$, is comprised of $d$ non-
dividend paying assets. Provided that the option has not been ‘knocked-out’,
the reward from stopping on period $t$ is essentially an increasing function
of the maximum value of the current assets,
$\max_{a\in\\{1,\ldots,d\\}}\xi_{t,a}$. However, if the maximum value of the
assets exceeds a prespecified barrier, the option becomes ‘knocked-out’ and
the reward becomes equal to zero for the remainder of the time horizon. In
this section, we compare various methods on several instances of these
stochastic optimal stopping problems in which the barrier threshold changes
over time (Kunitomo and Ikeda 1992).

To formalize the problem setting, let the components of the state vector
$\xi_{t}\in\Xi\equiv\mathbb{R}^{d}$ be referenced through the notation
$\xi_{t}\equiv(\xi_{t,1},\ldots,\xi_{t,d}),$ where $\xi_{t,a}\in\mathbb{R}$
represents the value of asset $a$ at exercise opportunity
$t\in\\{1,\ldots,T\\}$. The exercise opportunities are evenly spaced over a
calendar of $Y$ years and thus, defining $\lambda\triangleq Y/T$, it follows
that exercise opportunity $t$ occurs at the calendar time $\lambda t$. Let
$r\in[0,1)$ be the annualized discount rate, $K\geq 0$ be the strike price,
$B(t)\triangleq B_{0}e^{\delta\lambda t}$ be the barrier threshold on exercise
opportunity $t$, and $\xi\equiv(\xi_{1},\ldots,\xi_{T})$ be the sequence of
random states. With this notation, the reward function of this stochastic
optimal stopping problem is defined as

$\displaystyle g(t,\xi)$ $\displaystyle=\begin{cases}e^{-r\lambda
t}\max\left\\{0,\max\limits_{a\in\\{1,\ldots,d\\}}\xi_{t,a}-K\right\\},&\text{if
}\max\limits_{a\in\\{1,\ldots,d\\}}\xi_{s,a}\leq B(s)\text{ for all
}s\in\\{1,\ldots,t\\},\\\\[10.0pt] 0,&\text{otherwise}.\end{cases}$

Under the standard Black-Scholes setup, the sequence of random states
$\xi\equiv(\xi_{1},\ldots,\xi_{T})$ obeys a multidimensional geometric
Brownian motion where each asset $a\in\\{1,\ldots,d\\}$ has an initial value
of $\bar{x}$, drift equal to the risk-free rate $r$, and annualized volatility
equal to $\sigma_{a}$. The value of asset $a$ at exercise opportunity $t$ is
thus given by
$\xi_{t,a}\triangleq\bar{x}e^{\left(r-\sigma_{a}^{2}/2\right)\lambda
t+\sigma_{a}\lambda W_{t,a}}$, where each $W_{t,a}$ is a standard Brownian
motion process and the instantaneous correlation of $W_{\cdot,a}$ and
$W_{\cdot,a^{\prime}}$ is equal to $\rho_{a,a^{\prime}}$. We perform numerical
experiments on the following methods:

  * •

Robust Optimization (RO): The robust optimization approach is used here to
approximate (OPT) over the projected stochastic process
$x\equiv(x_{1},\ldots,x_{T})$, where
$x_{t}\triangleq\max_{a\in\\{1,\ldots,d\\}}\xi_{t,a}$ denotes the maximum
value of the assets in each exercise opportunity $t$. Note that this
stochastic process does not include information on whether the option has been
knocked-out. The robust optimization approach is run with robustness
parameters
$\epsilon\in\\{0\\}\cup\\{0.01,\ldots,0.09\\}\cup\\{0.1,\ldots,0.9\\}\cup\\{1,\ldots,10\\}$
and solved approximately using our heuristic ($\bar{\textnormal{H}}$) from
§4.2.3.

  * •

Least-Squares Regression (LS): The method of Longstaff and Schwartz is applied
to stochastic optimal stopping problem (OPT) over the Markovian stochastic
process $X\equiv(X_{1},\ldots,X_{T})$, where the state at each exercise
opportunity $X_{t}\triangleq(\xi_{t,1},\ldots,\xi_{t,d},q_{t})$ consists of
the values of the assets at the current exercise opportunity as well as an
indicator variable which equals one if the option has been knocked out:
$q_{t}\triangleq 1-\mathbb{I}\\{\max_{a\in\\{1,\ldots,d\\}}\xi_{s,a}\leq
B(s)\text{ for all }s\in\\{1,\ldots,t\\}\\}$. We consider the following basis
functions:

    * –

One: The constant function, 1.

    * –

KOind: An indicator variable for whether the option was knocked out, $q_{t}$.

    * –

Prices: The values of the assets, $\xi_{t,1},\ldots,\xi_{t,d}$.

    * –

PricesKO: The assets multiplied by the indicator variable,
$\xi_{t,1}(1-q_{t}),\ldots,\xi_{t,d}(1-q_{t})$.

    * –

MaxPrice: The maximum asset value,
$x_{t}\triangleq\max_{a\in\\{1,\ldots,d\\}}\xi_{t,a}$.

    * –

Payoff: The reward of exercising the option, $(1-q_{t})e^{-r\lambda
t}\max\left\\{0,x_{t}-K\right\\}$.

  * •

Pathwise Optimization (PO): The duality-based method of Desai et al. is
applied to the stochastic optimal stopping problem (OPT) over the Markovian
stochastic process $X\equiv(X_{1},\ldots,X_{T})$, where the state at each
exercise opportunity $X_{t}\triangleq(\xi_{t,1},\ldots,\xi_{t,d},q_{t})$
consists of the values of the assets at the current exercise opportunity as
well as an indicator variable which equals one if the option has been knocked
out: $q_{t}\triangleq 1-\mathbb{I}\\{\max_{a\in\\{1,\ldots,d\\}}\xi_{s,a}\leq
B(s)\text{ for all }s\in\\{1,\ldots,t\\}\\}$. We consider the same possible
basis functions for PO as were considered for LS.

  * •

Tree Method (Tree): The method of Ciocan and Mišić approximates the stochastic
optimal stopping problem (OPT) by restricting the space of exercise policies
to decision trees. Like the robust optimization approach, the tree approach is
used to find Markovian stopping rules for the optimal stopping problem over
the projected stochastic process $x\equiv(x_{1},\ldots,x_{T})$. We apply the
method with the same information as the robust optimization approach at each
time period (the current state $x_{t}$ and the time period $t$) and with a
splitting parameter of 0.005.

Our experiments and parameter settings closely parallel those of Ciocan and
Mišić (2020, §5), albeit with two primary differences. First, our experiments
consider a barrier threshold that changes as a function of time. Second, we
perform experiments in which the underlying assets have symmetrical as well as
asymmetrical annualized volatilities (i.e.,
$\sigma_{a}\neq\sigma_{a^{\prime}}$ for $a\neq a^{\prime}$). Our motivations
behind these differences are to compare the various methods in settings in
which the best Markovian stopping rules may not have a simple structure and in
which the low-dimensional projections $x_{1},\ldots,x_{T}$ are not sufficient
statistics for the optimal stopping rules (Broadie and Detemple 1997, §7). In
other words, our experiments aim to evaluate the performance of the proposed
robust optimization approach in settings where the best Markovian stopping
rules are not guaranteed to be optimal stopping rules for the optimal stopping
problem. Similarly as Ciocan and Mišić (2020, §5), we perform experiments on
various numbers of assets ($d\in\\{8,16,32\\}$) and various initial prices for
the assets ($\bar{x}\in\\{90,100,110\\}$).

In our implementation, we experimented with a variety of combinations of basis
functions, and, for the sake of brevity, we report results only for a subset
of combinations with the best performance. The expected rewards of the
stopping rules obtained by the various methods were estimated using a common
and independent testing dataset of $\tilde{N}=10^{5}$ sample paths, and all
experiments were repeated over 10 replications. The RO method is run in each
experiment with $N=10^{3}$ training sample paths and $\bar{N}=10^{3}$
validation sample paths (see §2.5 for additional details), the LS and Tree
methods are run in each experiment with $10^{5}$ training sample paths, and
the PO method is run in each experiment with $10^{3}$ outer sample paths and
$500$ inner sample paths. We use smaller numbers of training and outer sample
paths for RO and PO because the computation times of these methods scales
relatively quickly with the number of sample paths.

Figure 2:  Plot shows exercise policies obtained from solving a robust
optimization problem constructed from a training dataset of size $N=10^{3}$
and with the robustness parameter selected using a validation set of size
$\tilde{N}=10^{3}$. The problem parameters are $\bar{x}=100$ and $d=16$. The
remaining parameters are the same as those shown in Table
LABEL:tbl:barrier_symmetric. The thick black rectangles are the stopping
regions and the thin black line shows the barrier threshold.

In Table LABEL:tbl:barrier_symmetric, we present the expected rewards of
stopping rules obtained from the various methods in experiments with
symmetrical annualized volatilities.111111Expected rewards and computation
times for experiments with asymmetrical annualized volatility can be found in
Tables LABEL:tbl:barrier_asymmetric and LABEL:tbl:barrier_time_asymmetric in
Appendix 15. Across the parameter settings, the RO method yields stopping
rules with either the best or close-to-best expected reward relative to the
alternative state-of-the-art methods. The improvements of the robust
optimization approach over the state-of-the-art benchmarks are viewed as
particularly encouraging given the practical significance of this class of
options pricing problems. In Figure 2, we present visualizations of the
Markovian stopping rules produced by the RO method. These visualizations
provide interpretability to the stopping rules found by the RO method, which
is not possible for the LS method due to the high-dimensional state space.
Further visualizations of the stopping rules obtained by the RO method under
additional problem parameters are provided in Appendix 15.

(a)

(b)

Figure 3:  Each plot shows performance metrics of the robust optimization
problems constructed from a training dataset of size $N=10^{3}$ as a function
of the robustness parameter $\epsilon$. The left plot shows the robust
objective value and expected reward of policies obtained by solving the robust
optimization problem. The right plot shows the computation time for solving
the robust optimization problem. In both plots, the robust optimization
problem is solved approximately using our heuristic algorithm
($\bar{\textnormal{H}}$) from §4.2.3. The problem parameters are $\bar{x}=100$
and $d=16$, and the remaining parameters are the same as those shown in Table
LABEL:tbl:barrier_symmetric.

In Table LABEL:tbl:barrier_time_symmetric, we present the computation times of
the various methods. We observe from Table LABEL:tbl:barrier_time_symmetric
that the computation times of the RO method are highly competitive compared
with the computation times of the alternative methods when the initial price
is $\bar{x}=90$ or $\bar{x}=100$. Interestingly, the computation time of the
RO method increases significantly when the initial price is $\bar{x}=110$. To
understand why this is the case, we recall that the computation time for the
RO method is equal to the sum of the computation times for solving the
heuristic ($\bar{\textnormal{H}}$) over all considered choices of the
robustness parameter. Moreover, in Figure 3, we show that the computation time
of solving ($\bar{\textnormal{H}}$) is highly dependent on the robustness
parameter; indeed, most of the computation time for the RO method comes from
solving the robust optimization problem with unnecessarily large choices of
the robustness parameter. Finally, we show in Table 5.2 that the best choice
of the robustness parameter is smallest in the numerical experiments with
$\bar{x}=110$, even though the set of robustness parameters remains constant
across all of the numerical experiments. From these results, we conclude that
a dynamic search over the space of robustness parameters has the potential to
significantly reduce the computation times of the RO method, particularly the
case where $\bar{x}=110$. Additional numerical experiments which show the
relationship between the robustness parameter and the robust optimization
approach are provided in Appendix 15.

Barrier Option (Symmetric) - Best Choice of Robustness Parameter.  Initial
Price $d$ $\bar{x}=90$ $\bar{x}=100$ $\bar{x}=110$ 8 8.0 (2.45) 6.7 (1.25) 3.7
(1.06) 16 8.3 (1.89) 5.5 (1.27) 2.8 (1.03) 32 6.8 (1.32) 3.4 (0.84) 2.6 (0.70)
Best choice of robustness parameter $\epsilon$ found using the validation
method from §2.5 in the robust optimization problems constructed from a
training dataset of size $N=10^{3}$ and validation set of size
$\bar{N}=10^{3}$. The remaining parameters are the same as those shown in
Table LABEL:tbl:barrier_symmetric.

##  6 Conclusion

Over the past two decades, dynamic robust optimization has experienced a surge
of algorithmic advances. Until now, these advances from robust optimization
have not been harnessed to develop algorithms for stochastic dynamic
optimization problems with known probability distributions. In this paper, we
showed the value of bridging these traditionally separate fields of research,
in application to the classical and widely-studied problem of optimal
stopping. In this context, we devised new and theoretically-justified
algorithms for non-Markovian optimal stopping problems and highlighted the
performance of these new algorithms on stylized and well-studied problems from
options pricing. Along the way, we also developed novel theoretical and
computational results for solving dynamic robust optimization problems which
average over multiple uncertainty sets. We believe this work takes a
meaningful step towards broadening the impact of robust optimization to
address stochastic dynamic optimization problems of importance to industry.

## References

  * Adams and Sherali (1986) Adams WP, Sherali HD (1986) A tight linearization and an algorithm for zero-one quadratic programming problems. _Management Science_ 32(10):1274–1290.
  * Andersen (1999) Andersen L (1999) A simple approach to the pricing of bermudan swaptions in the multi-factor libor market model. _Journal of Computational Finance_ 3(2):5–32.
  * Andersen and Broadie (2004) Andersen L, Broadie M (2004) Primal-dual simulation algorithm for pricing multidimensional american options. _Management Science_ 50(9):1222–1234.
  * Bally and Pages (2003) Bally V, Pages G (2003) A quantization algorithm for solving multidimensional discrete-time optimal stopping problems. _Bernoulli_ 9(6):1003–1049.
  * Bayer et al. (2019) Bayer C, Häppölä J, Tempone R (2019) Implied stopping rules for american basket options from markovian projection. _Quantitative Finance_ 19(3):371–390.
  * Bayer et al. (2020) Bayer C, Tempone R, Wolfers S (2020) Pricing american options by exercise rate optimization. _Quantitative Finance_ 1–12.
  * Bayraktar and Yao (2014) Bayraktar E, Yao S (2014) On the robust optimal stopping problem. _SIAM Journal on Control and Optimization_ 52(5):3135–3175.