text:
\begin{document}

\title{Robust Multiple Stopping --- A Pathwise Duality Approach

}
\author{Roger J. A. Laeven\thanks{Corresponding author.}\\
{\footnotesize Dept. of Quantitative Economics}\\
{\footnotesize University of Amsterdam, EURANDOM}\\
{\footnotesize and CentER}\\
{\footnotesize {\tt R.J.A.Laeven@uva.nl}}\\[1mm]
\and John G. M. Schoenmakers\\
{\footnotesize Stochastic Algorithms and Nonparametric Statistics}\\
{\footnotesize Weierstrass Institute Berlin}\\
{\footnotesize {\tt schoenma@wias-berlin.de}}\\
\and Nikolaus Schweizer\\
{\footnotesize Dept. of Econometrics and Operations Research}\\
{\footnotesize Tilburg University}\\

{\footnotesize {\tt N.F.F.Schweizer@tilburguniversity.edu}}\\
\and Mitja Stadje\\
{\footnotesize Insurance Science and Financial Mathematics}\\
{\footnotesize Faculty of Mathematics and Economics}\\
{\footnotesize Ulm University}\\
{\footnotesize {\tt Mitja.Stadje@uni-ulm.de}}}
\date{This Version: \today}
\maketitle

\thispagestyle{empty}

\begin{abstract}
We develop a method to solve, theoretically and numerically, general optimal stopping problems.
Our general setting allows for multiple exercise rights, i.e., optimal multiple stopping,
for a robust evaluation that accounts for model uncertainty,
and for general reward processes driven by multi-dimensional jump-diffusions.
Our approach relies on first establishing robust martingale dual representation results
for the multiple stopping problem
that satisfy appealing pathwise optimality (i.e., almost sure) properties.
Next, we exploit these theoretical results to develop upper and lower bounds that,
as we formally show, not only converge to the true solution asymptotically,
but also constitute genuine pre-limiting upper and lower bounds.

We illustrate the applicability of our approach in a few examples
and analyze the impact of model uncertainty on optimal multiple stopping strategies.
\end{abstract}

\noindent\textbf{Keywords:} Optimal stopping; Multiple stopping;
Robustness; Model uncertainty; Ambiguity;
Pathwise duality; $g$-expectations; BSDEs; Regression.\\[3mm]
\noindent\textbf{AMS 2010 Subject Classification:} \textit{Primary}: 49L20; 60G40; 62L15;
\textit{Secondary}: 91B06; 91B16.\\[3mm]
\noindent\textbf{OR/MS Classification:} Dynamic programming/optimal control: Models;
Dynamic programming/optimal control: Applications;
Decision analysis: Risk.

\newpage

\clearpage\pagenumbering{arabic}

\setcounter{equation}{0}

\section{Introduction}

In this paper we analyze general optimal stopping problems of the following form:
\begin{equation}
Y^{\ast,L}_{t}:=\sup_{\substack{t\leq\tau_{1}<\cdots<\tau_{L} \\
(\tau_{1},\ldots,\tau_{L})\in\mathcal{T}_{t}(L)}}
\sup_{\mathbb{Q}\in
\mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^{L}H_{\tau_{l}}\Big|\mathcal{F}_{t}\right],\qquad
0\leq t\leq T, \label{eq:problem}
\end{equation}
where $\mathcal{T}_{t}(L)$ is a family of stopping time vectors,
$L$ is a number of exercise rights,
$T<\infty$ is a fixed time horizon,

$\mathcal{Q}$ is a family of probabilistic models,
and $H$ is a general $\mathcal{F}_{\centerdot}$-adapted reward process.
(The operator $\sup$ is to be understood as $\esssup$ if it applies to an uncountable family of random variables.)
The optimal stopping problem \eqref{eq:problem} features generality along three dimensions:
(i) it allows for optimal multiple stopping (when $L>1$),
(ii) it allows for a robust evaluation that explicitly takes probabilistic model uncertainty (ambiguity) into account
(when $\mathcal{Q}$ is not a singleton), and
(iii) it allows for general reward processes that will be driven by
multi-dimensional jump-diffusion processes.
The process $Y^{\ast,L}_{t}$ is referred to as the upper Snell envelope of $H$ due to $L$ exercise rights after the seminal work of Snell \cite{S52}.
Problems of this type, or special cases thereof, occur naturally in a wide variety of applications
in probability, operations research, economics and finance.

Our aim is to develop upper and lower bounds on $Y^{\ast,L}_{t}$ that satisfy several desirable properties.

We achieve this by first establishing
suitable martingale dual representations for problem \eqref{eq:problem}
that can be viewed as significant generalizations of the
classical additive dual representations for \textit{standard}
(i.e., $\mathcal{Q}=\{\mathbb{Q}\}$ as opposed to \textit{robust}) optimal stopping problems,
developed independently by Rogers \cite{R02} and Haugh and Kogan \cite{HK04}
(see also the early Davis and Karatzas \cite{DK94}) and their extension to
standard multiple stopping problems in Schoenmakers \cite{S12}.

(A multiplicative dual representation for standard optimal stopping problems
was proposed by Jamshidian \cite{J07}.)
Our
dual representations take the
form of an infimum over (robust) martingales,

with no appearance of stopping times.

An appealing feature---both theoretically and for numerical stability---of the
dual representations we establish is their pathwise optimality,
i.e., their almost sure property.
Already when $L=1$ these results are new and of independent interest
for robust optimal single stopping.
They are developed here in the general setting of robust optimal multiple stopping
\eqref{eq:problem}.
The almost sure nature of the dual representations suggests that finding
a `good'
martingale that is `close' to a `surely optimal' martingale will yield tight and nearly constant upper bounds.

The target can be the unique (robust) Doob martingale,
to be constructed from an approximation to the upper Snell envelope or,
more generally, a martingale for which the dual representation's infimum is attained
and the almost sure property is satisfied.
While this phenomenon of tightness and constancy is known in the case of standard, non-robust single stopping problems
(i.e., when $\mathcal{Q}=\{\mathbb{Q}\}$ and $L=1$, see
Rogers \cite{Rog10} and
Schoenmakers, Zhang and Huang \cite{SZH13}),
we will analyze
it
in our general setting of robust
stopping.

We will show in particular that a low (vanishing in probability) robust variance implies a
tight (converging in $L^{1}$) approximation.
The mathematical details of these results are delicate.

These new theoretical results justify and enable us to next develop a numerically
implementable method to obtain upper and lower bounds to $Y^{\ast,L}_{t}$ with
desirable properties.
Our lower bound, derived from the proposed exercise strategy, will, as we formally show,
not only converge to the
optimal solution asymptotically but also be `biased low' at the pre-limiting level in a Brownian-Poisson filtration.
This is not the case for the initially proposed upper bound: it converges to
the true solution but is not in general `biased high'.
We therefore also develop a second upper bound
that as we prove both converges to the true solution asymptotically
and is biased high in a Brownian-Poisson filtration.
It is based on a Lipschitzian $L^{2}$-approximation,
with a Lipschitz constant that we explicitly derive,
and a suitable (reversed) application of Jensen's inequality.
We will refer to this second upper bound as our \textit{genuine} upper bound.
The computational complexity of both upper bounds is only linear in the number of exercise rights,
and our method does not require nested simulation.

We provide extensive numerical examples, including
single and multiple stopping problems, univariate and multivariate stochastic drivers,
increasing and decreasing reward functions,
and pure diffusion and jump-diffusion models,
to illustrate the applicability and generality of our approach.
They demonstrate that our
approach yields
tight upper and lower bounds that, due to almost sure properties,
moreover have
low standard errors.
They also analyze the impact of multiple vs. single stopping rights
and reveal that employing a robust evaluation that takes ambiguity into account
is highly relevant for optimal stopping, especially in the presence of multiple exercise rights.

Embedded in a Brownian-Poisson filtration, the problems we analyze are naturally represented as stopping problems with respect to
$g$-expectations (Peng \cite{P97,P04}),
leading to backward stochastic differential equations (BSDEs).
Hence, we explicitly construct novel genuine upper and lower bounds to BSDE solutions
with positively homogeneous convex drivers
in a Brownian-Poisson filtration,
as a development of independent interest.
Bender, Schweizer and Zhuo \cite{BSZ17},
when analyzing solutions to discrete-time (reflected) BS$\Delta$Es
rather than the continuous-time BSDEs we consider,
develop upper and lower bounds applying techniques different from the techniques we employ.
Bender, G\"artner and Schweizer \cite{BGS18} construct Monte Carlo upper and lower bounds
for a class of discrete-time stochastic dynamic programs
which includes
discretizations of multiple stopping problems.
Our genuine upper and lower bounds apply directly to our original
continuous-time problem.
Our genuine lower bound takes advantage of an almost sure property of a `second kind'
that we formally establish in order to reduce its variance---`second kind' to distinguish it from the additive dual representation's almost sure property.
This almost sure property entails that the difference between the BSDEs terminal condition and the associated (robust) martingale
is constant almost surely.
Our genuine upper bound for the continuous-time problem is based on forward simulation of an approximate BSDE solution.
The construction is somewhat related to the a posteriori criterion
for error evaluation
introduced in Bender and Steiner \cite{BS13}

in a Brownian filtration,
and developed here to obtain explicit genuine upper bounds for BSDEs
in a Brownian-Poisson filtration.

The development of numerically implementable methods to obtain approximations to problems
of the type \eqref{eq:problem} but with $\mathcal{Q}$ a singleton (no ambiguity),
$L=1$ (single stopping), and with $H$ multi-dimensional but satisfying strong conditions,
started with the regression-based Monte Carlo methods of Carriere \cite{C96} and Longstaff and Schwartz \cite{LS01};
see also Tsitsiklis and Van Roy \cite{TR00} and Cl\'ement, Lamberton and Protter \cite{CLP02}.
These methods yield lower bounds to $Y^{\ast,1}_{t}$ by approximating the optimal stopping time using regression
and are commonly referred to as ``primal'' approaches.
An important example of a non-regression based primal approach
is the stochastic mesh method of Broadie and Glasserman \cite{BG04}
(see, for further details, Glasserman \cite{G04} and also Belomestny, Kaledin and Schoenmakers \cite{BKS20}).

``Dual'' algorithms that exploit additive dual representations
to numerically compute upper bounds were first proposed by Andersen and Broadie \cite{AB04}
in the standard single stopping problem and were further developed by e.g., Belomestny, Bender and Schoenmakers \cite{BBS09}
to allow for non-nested simulation.
While primal methods rely in a sense on constructing an
appropriate stopping time, dual methods rely on constructing an appropriate martingale.
Brown, Smith and Sun \cite{BSS10} in an innovative paper enlarge
the information on which an exercise decision may depend in dual optimization,
yielding tight upper bounds.

Model uncertainty, and the distinction between risk

and ambiguity,

has received much
attention in recent years.
Under the Bayesian paradigm, as adopted in Savage's
\cite{S54} subjective expected utility model, this distinction is, in a sense,
nullified, through subjective probabilities resulting from a subjective prior
probability over probabilistic models that quantifies model uncertainty.
A popular approach beyond the Bayesian paradigm is provided by the multiple
priors model of Gilboa and Schmeidler \cite{GS89}, which is a
decision-theoretic formalization of the classical Waldian maxmin decision rule
(Wald \cite{W50}; see also Huber \cite{H81}) and experimentally motivated by
the Ellsberg \cite{E61} paradox.
These models are intimately related to
coherent, convex and entropy convex measures of risk in financial risk
measurement (F\"ollmer and Schied \cite{FS02,FS11}, Frittelli and Rosazza
Gianin \cite{FG02}, Ruszczy\'nski and Shapiro \cite{RS06a,RS06b}, and Laeven
and Stadje \cite{LS13}).
They explicitly recognize that probabilistic models
may be misspecified and are often referred to as \textit{robust} approaches
(Hansen and Sargent \cite{HS07}).
The literature on robust single
stopping theory is rapidly growing; it includes Riedel \cite{R09},
Kr\"atschmer and Schoenmakers \cite{KS10}, Bayraktar, Karatzas and Yao
\cite{BKY10}, Bayraktar and Yao \cite{BY11}, Cheng and Riedel \cite{CR13},
{\O }ksendal, Sulem and Zhang \cite{OSZ14}, Belomestny and Kr\"atschmer
\cite{BK16a,BK16b},

Bayraktar and Yao \cite{BY14,BY15a,BY15b}, Ekren, Touzi, and Zhang
\cite{ETZ14}, Matoussi, Piozin, and Possama\"{\i} \cite{MPP14}, Matoussi,
Possama\"{\i}, and Zhou \cite{MPZ13}, and Nutz and Zhang \cite{NZ15}.
However, numerically implementable methods to solve general optimal stopping problems
of the form \eqref{eq:problem} have not been well-developed as yet.
Kr\"atschmer \textit{et al.} \cite{KLLSS18} propose a numerically
implementable method for single stopping problems under uncertainty in drift
and jump intensity.
Their approach is dual but not path-wise, i.e., it does not rely on a dual representation with the appealing almost sure property,
and cannot handle multiple stopping problems.

The multiple stopping problem can be viewed as $L$ nested single stopping problems,
where the decision-maker first chooses between stopping at time $\tau_{1}$ on the one hand,
thus collecting the reward and entering into a new contract with $L-1$ exercise rights,
and retaining $L$ exercise rights on the other hand, and so on.
Multiple exercise rights occur naturally in many applications across various fields.
For example,
in environmental economics,
a swing option gives the investor the right to change his purchased
energy quantity a number of times per time period;
in finance, a flexible interest rate cap gives the investor
the right to exercise at each interest rate reset date a number of times over the life of the contract;
and in insurance, a partial surrender option provides a payoff to the policyholder
each time he partially surrenders his life insurance contract;
see e.g., Carmona and Dayanik \cite{CD08} and Carmona and Touzi \cite{CT08} and the references therein.

Kobylanski, Quenez and Rouy-Mironescu \cite{KQR11} analyze the standard multiple stopping
problem (without ambiguity) allowing the payoff to be
a general functional of an ordered sequence of stopping times.
Bender, Schoenmakers and Zhang \cite{BSZ15} develop a dual approach to
generalized multiple stopping problems
with respect to standard conditional expectations
that is intimately related to the information
relaxation approach of Brown, Smith and Sun \cite{BSS10}.

A primal-dual algorithm for standard multiple stopping with respect to standard conditional expectations
in the context of flexible interest rate caps has been proposed in Balder, Mahayni and Schoenmakers \cite{BMS13}.

As an important application, our approach may be used for robust no-arbitrage
pricing (Hansen and Jagannathan \cite{HJ91},
Cochrane and Sa\'{a}-Requejo \cite{CS00})
of American-style derivatives with possibly multiple exercise rights,
via superhedging.
This entails a significant advancement of the standard
approach, where in a usually incomplete market the corresponding stopping
problem is solved with respect to an arbitrarily chosen (local equivalent
martingale) measure.
Our results can also be applied to indifference
valuation (seller's perspective; Carmona \cite{C09}, Laeven and Stadje \cite{LS14}) of general optimally stopped reward processes
under the multiple priors model.
Another application is that
of robust risk measurement (Ben-Tal and Nemirovski \cite{BN98}, Bertsimas and Brown \cite{BB09}, F\"ollmer and Schied \cite{FS11})
to determine e.g., the risk capital required to cover
optimally stopped reward processes.

The remainder of this paper is organized as follows.
In Section \ref{sec:robstop} we recall some basic notions, establish some general properties,
introduce the robust optimal multiple stopping problem,
and provide some examples.
In Section \ref{sec:pathdual}, we present our pathwise dual representations
and establish our results on surely optimal (robust) martingales.
In Section \ref{sec:MSAL}, we outline a general primal-dual algorithm and prove its convergence.
Section \ref{sec:BPcon} presents explicit upper and lower bounds in a Brownian-Poisson filtration.

Section \ref{sec:num} provides extensive numerical results.

All proofs and several auxiliary results are in the Online Appendix.

\setcounter{equation}{0}

\section{Robust Optimal Multiple Stopping}

\label{sec:robstop}

\subsection{Basic Notions and General Properties}

We start by considering a general stochastic setup.
We let $\big(\Omega,\left(\mathcal{F}_{t}\right)_{t\in\left\{0,\ldots,T\right\}},\mathbb{P}\big)$
be a filtered probability space and
let $\mathfrak{X}$ be a linear subspace of $L^{0}(\Omega,\mathcal{F},\mathbb{P})$
with $\mathcal{F}:=\mathcal{F}_{T}$.
We further assume that $\mathfrak{X}$ has a lattice structure,
i.e., $\mathfrak{X}$ is closed under the operations
$\wedge$ (min) and $\vee$ (max), and that $\mathfrak{X}$ contains all
indicator functions $1_{A},$ $A\in\mathcal{F}$.
((In)equalities
between random variables are understood in the $\mathbb{P}$-almost sure sense,
often without explicit mention.)

To represent preferences, we consider
a family of mappings $\rho:=\left(\rho_{t}\right)_{t=0,\ldots,T}$,
\begin{equation*}
\rho_{t}:\mathfrak{X}\rightarrow\mathfrak{X}\cap L^{0}(\Omega,\mathcal{F}_{t},\mathbb{P}).
\end{equation*}
It is referred to as a monotone, regular, recursive, conditional translation invariant
\emph{dynamic monetary utility functional},
henceforth \emph{DMU} for short,
if it satisfies the following conditions:

\begin{enumerate}
\item[(C1)] $\rho_{t}\left(  X\right)  \leq\rho_{t}\left(  Y\right)  $ for all
$X,Y\in\mathfrak{X}$ with $X\leq Y$
and $t\in\left\{  0,\ldots,T\right\}
$ \emph{(monotonicity)}.

\item[(C2)] $\rho_{t}\left(  1_{A}X\right)  =1_{A}\rho_{t}\left(  X\right)  $
for all $X\in\mathfrak{X}$, $A$ $\in\mathcal{F}_{t}$ and $t\in\left\{
0,\ldots,T\right\}  $ \emph{(regularity).}

\item[(C3)] $\rho_{t}=\rho_{t}\circ\rho_{t+1}$ for all $t\in\left\{
0,\ldots,T-1\right\}  $ \emph{(recursiveness).}

\item[(C4)] $\rho_{t}\left(  X+Y\right)  =\rho_{t}\left(  X\right)  +Y$ for
all $X,Y\in\mathfrak{X}$ with $Y$ $\in\mathcal{F}_{t}$ and $t\in\left\{
0,\ldots,T\right\}  $ \emph{(conditional translation invariance).}
\end{enumerate}

As additional
properties we consider:

\begin{enumerate}
\item[(P1)] $\rho_{t}(X+Y)\leq\rho_{t}(X)+\rho_{t}(Y)$ for all $X,Y\in\mathfrak{X}$
and $t\in\left\{0,\ldots,T\right\}$ (\emph{subadditivity}).

\item[(P2)]

$\left[  X\leq0\text{ \ \ and \ \ }\rho_{t}(X)\geq0\right]  \Longrightarrow
X=0,\text{ \ \ for all }X\in\mathfrak{X}\text{ and }t\in\left\{
0,\ldots,T\right\}$ (\emph{sensitivity}).

\item[(P3)]
$\rho_{t}\left(\lambda X\right) =\lambda\rho_{t}\left(X\right)$ for all $X\in\mathfrak{X},$
$\lambda\geq0$ and $t\in\left\{0,\ldots,T\right\}$

(\emph{positive homogeneity}).
\end{enumerate}
Conditions (C1)--(C4) will always be assumed.
In the sequel, we will mention explicitly which of the properties (P1)--(P3) is required.
Properties (P1) and (P2) also entail the implication
$\left[  X\geq0\text{ \ \ and \ \ }\rho_{t}(X)\leq0\right]  \Longrightarrow X=0$;
see Lemma \ref{lem:P2a} in Appendix \ref{sec:A1}.
DMUs that satisfy (P1)--(P3), in addition to (C1)--(C4), take the form of robust, or worst case, expectations
and have been widely used in applied probability, operations research, economics and finance;
see the references in the Introduction and Section \ref{sec:examples} below.

In this paper, we will frequently use the following implications of (C2) and (C4):
\begin{enumerate}
\item[(C5)] $\rho_{t}\left(0\right)=0$ for all $t\in\left\{0,\ldots,T\right\}$
(\emph{normalization}).

\item[(C6)] $\rho_{t}\left(X\right)=X$ for all $X\in\mathfrak{X}$ with $X$
$\in\mathcal{F}_{t}$ and $t\in\left\{0,\ldots,T\right\}$
(\emph{$\mathcal{F}_{t}$-invariance}).
\end{enumerate}

Let $\mathcal{H}$ be the set of adapted processes $\left(U_{t}\right)_{t\in\{0,\ldots,T\}}$
such that $U_{t}\in\mathfrak{X}\cap L^{0}(\Omega,\mathcal{F}_{t},\mathbb{P})$.
A process $M=\left(M_{t}\right)_{t\in\{0,\ldots,T\}}\in\mathcal{H}$ is said to be a $\rho$-martingale if
\begin{equation}
M_{t}=\rho_{t}(M_{t+1}),\qquad 0\leq t<T.
\label{eq:rhomartingale}
\end{equation}

We present two auxiliary lemmas.
The first lemma provides a generalization of Doob's optional sampling theorem
towards our setup:
\begin{lemma}
\label{lem:Doob} (Doob) Suppose $\rho$ satisfies (C1)--(C4).
Then, for any $\rho$-martingale $M$
and any stopping time $\tau_{i}$, $i\leq\tau_{i}\leq T$, it holds that $\rho_{i}\left(
M_{\tau_{i}}\right)  =M_{i},$ $0\leq i\leq T.$
\end{lemma}
Due to the next lemma, the properties of recursiveness (C3)
and conditional translation invariance (C4) carry over to stopping times, as we will exploit later:
\begin{lemma}
\label{L3} Let $\rho$ satisfy (C1)--(C4), and let $t\in\left\{  0,\ldots,T\right\}$ be fixed.
Consider, for any stopping time $\tau$, $t\leq\tau\leq T$, the functional
\begin{equation*}
\rho_{\tau}(X):=\sum_{j=t}^{T}1_{\tau=j}\rho_{j}(X).
\end{equation*}
Then, $\rho_{\tau}$ acts from $\mathcal{F}_{T}\rightarrow\mathcal{F}
_{\tau}\supset\mathcal{F}_{t}$, and
\begin{itemize}
\item[(i)] $\rho_{\tau}$ satisfies $\rho_{t} = \rho_{t}\circ\rho_{\tau}$;
\item[(ii)] $\rho_{\tau}(X+Y) = X+\rho_{\tau}(Y),
\text{ \ for \ }X\in\mathcal{F}_{\tau},\text{ \ \ }Y\in\mathcal{F}_{T}$.
\end{itemize}
\end{lemma}

\subsection{The Stopping Problem}

Consider a fixed adapted reward, or (discounted) cash-flow, process
$H=\left(H_{t}\right)_{t\in\{0,\ldots,T\}}\in\mathcal{H}$
and a DMU decision-maker with $L$ exercise rights
that have to be exercised at different exercise dates.
For each fixed $t$ and $L$, $0\leq t\leq T$,
let $\mathcal{T}_{t}(L)$ be the family of stopping vectors
$\left(\tau_{1},\ldots,\tau_{L}\right)$
such that $\tau_{1}\geq t$ and $\tau_{l}\geq\tau_{l-1}+1$ for all $l$, $1<l\leq L$.
The decision-maker faces the following robust optimal multiple stopping problem:

\begin{equation}
Y_{t}^{\ast, L}:=

\underset{t\leq\tau_{1}<\tau_{2}<\cdots<\tau_{L}}{\esssup }\text{ }\rho_{t}
\left(\sum_{l=1}^{L}H_{\tau_{l}}\right),\qquad t\in\{0,\ldots,T\},
\label{eq:multiplestop}
\end{equation}
for a DMU functional $\rho$ that satisfies (C1)--(C4).
(We note that problem \eqref{eq:multiplestop} is even slightly more general than problem \eqref{eq:problem},
which arises when additionally (P1)--(P3) are satisfied.)
Henceforth, we write $\sup$ (and $\inf$) instead of $\esssup$ (and $\essinf$) for convenience,
understanding that they apply to an uncountable family of random variables.
For a clean formulation of the multiple stopping problem (\ref{eq:multiplestop}),
we extend the cash-flow process by setting $H_{j}\equiv 0$ and $\mathcal{F}_{j}\equiv\mathcal{F}_{T}$, for $j=T+1,T+2,\ldots$.
That is, the subset of rights $l$, $l=2,\ldots,L$, not exercised by time $T$ become valueless.
Hence, for any $\rho$-martingale $M$, $M_{j}=M_{T}$, $j>T$.

When $L\equiv 1$, the single stopping problem
\begin{equation}
Y_{t}^{\ast}\equiv Y_{t}^{\ast,1}=\sup_{\tau\in\mathcal{T}_{t}}\rho_{t}(H_{\tau}),\qquad
t\in\{0,\ldots,T\},
\label{eq:OST}
\end{equation}
occurs as a special case,
where the family of stopping times $\mathcal{T}_{t}\equiv\mathcal{T}_{t}(1)$
takes values in the set $\left\{t,\ldots,T\right\}$.

The multiple stopping problem can be viewed as $L$ nested single stopping problems with only a single exercise right.
Indeed, setting $Y^{\ast,0}\equiv 0$, $Y^{\ast,1}\equiv Y^{\ast}$
is the upper Snell envelope of $H$ due to a single exercise right.
Then, for multiple exercise rights $L\geq 1$,
$Y^{\ast,L}$ can be viewed as the upper Snell envelope
of the process
\begin{equation*}
H_{t}+\rho_{t}\left(Y_{t+1}^{\ast,L-1}\right),\qquad t\in\{0,\ldots,T-1\},
\end{equation*}
due to only a single exercise right.

Let us denote the set of $\rho$-martingales $M$ with $M_{0}=0$ by $\mathcal{M}_{0}^{\rho}$.
There exists a unique $\rho$-martingale $M^{\ast\rho}\in\mathcal{M}_{0}^{\rho}$
and a non-decreasing predictable $A^{\ast\rho}\in\mathcal{H}$ such that
\begin{equation}
Y_{t}^{\ast}=Y_{0}^{\ast}+M_{t}^{\ast\rho}-A_{t}^{\ast\rho},\qquad
t\in\{0,\ldots,T\},
\label{eq:DoobMeyerdisc}
\end{equation}
which represents the $\rho$-Doob decomposition of
$Y^{\ast}=\left(Y_{t}^{\ast}\right)_{t\in\{0,\ldots,T\}}$.
It is easy to verify that, for $t\in\{0,\ldots,T-1\}$,
\begin{equation}
M_{t+1}^{\ast\rho}-M_{t}^{\ast\rho}=Y_{t+1}^{\ast}-\rho_{t}\left(Y_{t+1}^{\ast}\right),
\text{ \ and \ }
A_{t+1}^{\ast\rho}-A_{t}^{\ast\rho}=Y_{t}^{\ast}-\rho_{t}\left(  Y_{t+1}^{\ast}\right).
\label{eq:DD}
\end{equation}
Henceforth, the $\rho$-martingale $M^{\ast\rho}$ will often be referred to as the $\rho$-Doob martingale
and we often suppress its superscript $\rho$ to simplify notation.

In Appendix \ref{sec:singlestop}, we
establish some auxiliary results for problem \eqref{eq:OST}
that will be exploited in the proofs of the results that follow.

\subsection{Examples}

\label{sec:examples}

We provide the following examples in which specific versions of the robust optimal multiple stopping problem
of the general form \eqref{eq:multiplestop}

occur naturally:

\begin{enumerate}
\item[(A.)] No-arbitrage pricing: Let $\mathcal{Q}$ be the set
of local equivalent martingale measures.
(Only if markets are complete $\mathcal{Q}$ is a singleton,
i.e., $\mathcal{Q}=\{\mathbb{Q}\}$.)
Then, the \textit{superhedging} price $\pi^{L}$ of a contract with $L\geq 1$ exercise rights
and associated payoff $\sum_{l=1}^L H_{\tau_l}$ is given by
\begin{equation*}
\pi^{L}=\sup_{\tau_1 < \tau_2 < \cdots < \tau_L}\sup_{\mathbb{Q}\in\mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^L H_{\tau_l}\right].
\end{equation*}
Many different approaches to no-arbitrage pricing have been proposed in the
literature; see, e.g., the good-deal bounds of Cochrane and Sa\'{a}-Requejo \cite{CS00},
Hansen and Jagannathan \cite{HJ91} and Bj\"ork and Slinko \cite{BS06},
or the acceptable opportunities of Carr, Geman and Madan \cite{CGM01}.
All these approaches yield prices of the form
\begin{equation*}
\tilde{\pi}^{L}=\sup_{\tau_1 < \tau_2 < \cdots < \tau_L}\sup_{\mathbb{Q}\in\mathcal{Q}_{\mathrm{restricted}}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^L H_{\tau_l}\right],
\end{equation*}
where $\mathcal{Q}_{\mathrm{restricted}}\subset\mathcal{Q}$.

Prototypical situations leading to single and multiple stopping problems in economics and finance
are the pricing and exercising of American-style, Bermudan-style, and swing options.
American options give the holder the
right to exercise the option (once)
on any preferred trading day before expiration.
Different from American options, Bermudan options prescribe a set of trading days
on which the option can be exercised (once).
Swing options, more generally, give the holder the right to exercise the option multiple times,
at a pre-specified set of exercise dates.
With $L\geq 1$ exercise rights, exercised at $\tau_{1}<\tau_{2}<\cdots<\tau_{L}$,
the payoff equals $\sum_{l=1}^{L}H_{\tau_{l}}$ for a cash-flow process $H\in\mathcal{H}$.
Swing options are particularly popular in energy markets
to manage the risk of fluctuations in oil, gas, or electricity prices.

\item[(B.)] Indifference valuation---the seller's perspective:
Suppose that the seller of a contract has a max-min utility functional
of the form
\begin{equation*}
U(H)=\inf_{\mathbb{Q} \in \mathcal{Q}} \mathbb{E}_{\mathbb{Q}}[H],
\end{equation*}
for a family of probabilistic models (i.e., priors) $\mathcal{Q}$
and adopts a utility indifference valuation approach (Carmona \cite{C09}, Laeven and Stadje \cite{LS14}).
Then, the value $V^{L}$ of a contract with $L\geq 1$ exercise rights
and associated payoff $\sum_{l=1}^L H_{\tau_l}$
is determined from the indifference relation
\begin{equation*}
U(0)=\inf_{\tau_1 < \tau_2 < \cdots < \tau_L} U\left(-\sum_{l=1}^L H_{\tau_l}+ V^{L}\right)
=\inf_{\tau_1 < \tau_2 < \cdots < \tau_L}\inf_{\mathbb{Q} \in \mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[-\sum_{l=1}^L H_{\tau_l}+V^{L}\right].
\end{equation*}
Hence,
\begin{equation*}
V^{L}=\sup_{\tau_1 < \tau_2 < \cdots < \tau_L}\sup_{\mathbb{Q} \in \mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^L H_{\tau_l}\right].
\end{equation*}

\item[(C.)] Robust risk measurement:
Suppose that $\rho$ is a robust, or worst case, expectation, that is,
\begin{equation}
\rho(H)=\sup_{\mathbb{Q} \in \mathcal{Q}} \mathbb{E}_{\mathbb{Q}}[H],
\label{eq:coherentrm}
\end{equation}
for a family of probabilistic models $\mathcal{Q}$.
In financial risk measurement, \eqref{eq:coherentrm} is referred to as a coherent risk measure
and $\mathcal{Q}$ as a set of generalized scenarios
(Artzner \textit{et al.} \cite{ADEH99} and F\"ollmer and Schied \cite{FS11});
see also Ben-Tal and Nemirovski \cite{BN98} for the intimately connected
robust optimization paradigm.
It determines the minimal amount of risk capital required to be added to the financial position $H$
to make it `safe' from the viewpoint of the regulatory authority.
Applications of coherent risk measures and generalized scenarios to
decision and optimization include
Lesnevski, Nelson and Staum \cite{LNS07},
Bertsimas and Brown \cite{BB09},

Choi, Ruszczy\'nski and Zhao \cite{CRZ11},
Philpott, de Matos and Finardi \cite{PMF13}
and Tekaya, Shapiro, Soares and da Costa \cite{TSPP13}.
Assume now that $H_{\tau_{l}}$ is a payout obligation (i.e., liability)
at time $\tau_{l}$,
where $\tau_{l}$ is a stopping time,
due to e.g., a flexible interest rate cap in interest rate markets

or a partial surrender option in life insurance
to be paid to a policyholder who decides to partially surrender his insurance contract.
Then, the required amount of risk capital due to $L\geq 1$ stopping rights is given by
\begin{equation*}
\sup_{\tau_1 < \tau_2 < \cdots < \tau_L}\rho\left(\sum_{l=1}^L H_{\tau_l}\right)
=\sup_{\tau_1 < \tau_2 < \cdots < \tau_L}\sup_{\mathbb{Q} \in \mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^L H_{\tau_l}\right].
\end{equation*}
\end{enumerate}

\setcounter{equation}{0}

\section{Pathwise Duality}

\label{sec:pathdual}

\subsection{Pathwise Dual Representation}

\label{sec:pathdualrep}

The following theorem establishes our pathwise (i.e., almost sure)
additive dual representation for general multiple stopping problems of the form \eqref{eq:multiplestop}.
\begin{theorem}
\label{th:mulsub}
Suppose
$\rho$ satisfies (C1)--(C4) and is subadditive (P1).
Then,
for any adapted process $H=\left(H_{t}\right)_{t\in\{0,\ldots,T\}}\in\mathcal{H}$
and each fixed $t\in\left\{ 0,\ldots,T\right\}$,
\begin{itemize}
\item[(i)] we have the dual representation
\begin{equation}
Y_{t}^{\ast, L}=\underset{M^{(1)},\ldots,M^{(L)}\in\mathcal{M}_{0}^{\rho}}{\inf}\rho_{t}\left(\max_{t\leq j_{1}<j_{2}<\cdots< j_{L}}\sum_{k=1}^{L}
\left(H_{j_{k}}+M_{j_{k-1}}^{(k)}-M_{j_{k}}^{(k)}\right)\right);
\end{equation}
\item[(ii)] the dual representation's infimum is attained:
\begin{equation}
Y_{t}^{\ast, L}=\rho_{t}\left(\max_{t\leq j_{1}<j_{2}<\cdots< j_{L}}
\sum_{k=1}^{L}\left(H_{j_{k}}+M_{j_{k-1}}^{\ast, L-k+1}-M_{j_{k}}^{\ast, L-k+1}\right)\right);
\end{equation}
\item[(iii)] if in addition $\rho$ is sensitive (P2), we have the pathwise dual representation
\begin{equation}
Y_{t}^{\ast, L}=\max_{t\leq j_{1}<j_{2}<\cdots< j_{L}}
\sum_{k=1}^{L}\left(H_{j_{k}}+M_{j_{k-1}}^{\ast, L-k+1}-M_{j_{k}}^{\ast, L-k+1}\right),\qquad\mathrm{almost\ surely};
\end{equation}
\end{itemize}
where the $\rho$-martingales $M^{\ast, L-k+1}$ satisfy
\begin{equation}
M_{r+1}^{\ast, L-k+1}-M_{r}^{\ast, L-k+1}=Y_{r+1}^{\ast, L-k+1}-\rho_{r}
\left(Y_{r+1}^{\ast, L-k+1}\right),
\label{eq:DoobM}
\end{equation}
and $Y^{\ast, L-k+1}$ is the upper Snell envelope due to $L-k+1$ exercise rights,
satisfying the Bellman principle,
\begin{equation}
Y_{r}^{\ast, L-k+1}=\max\left[  H_{r}+\rho_{r}\left(Y_{r+1}^{\ast,L-k}\right),
\rho_{r}\left(Y_{r+1}^{\ast, L-k+1}\right)\right].
\label{eq:Bellman}
\end{equation}
\end{theorem}

Already when $L\equiv 1$, Theorem \ref{th:mulsub} is new and of significant independent interest.
In this case it simplifies to:
\begin{corollary}
\label{P1} Suppose
$\rho$ satisfies (C1)--(C4) and
(P1).
Then, for any adapted process $H\in\mathcal{H}$ and each
fixed $t\in\left\{ 0,\ldots,T\right\}$, we have the
dual representation
\begin{align}
Y_{t}^{\ast}  &  =\inf_{M\in\mathcal{M}_{0}^{\rho}}\rho_{t}\left(  \max_{t\leq
j\leq T}\left(  H_{j}+M_{t}-M_{j}\right)  \right) \label{as0}\\
&  =\rho_{t}\left(  \max_{t\leq j\leq T}\left(  H_{j}+M_{t}^{\ast}-M_{j}
^{\ast}\right)  \right), \label{as00}
\end{align}
where $M^{\ast}$
is the $\rho$-Doob martingale in Eqn.~\eqref{eq:DoobMeyerdisc}.
If in addition $\rho$ is sensitive (P2), we have
the almost sure property:
\begin{equation}
Y_{t}^{\ast}=\max_{t\leq j\leq T}\left(  H_{j}+M_{t}^{\ast}-M_{j}^{\ast
}\right),\qquad\mathrm{almost\ surely.} \label{as}
\end{equation}
\end{corollary}

\begin{remark}
We note that the
single stopping problem also admits
an alternative but non-pathwise additive dual representation
for functionals $\rho$ satisfying (C1)--(C4);

see Proposition~\ref{prop:ADP} in the Appendix.
In Theorem~\ref{th:mulsub} and Corollary~\ref{P1}, the subadditivity property (P1) of $\rho$ is required,
and exploited through application of Lemma~\ref{L2} in the proof of Theorem~\ref{th:mulsub}.

A dual representation theorem in the spirit of Theorem~\ref{th:mulsub} and Corollary~\ref{P1}
without assuming (P1) seems not possible to us.
\end{remark}

\subsection{Surely Optimal \texorpdfstring{$\rho$-}{rho-}Martingales}

\label{sec:suroptmart}

The $\rho$-Doob martingale in Eqn.~\eqref{eq:DoobMeyerdisc} plays a special role
in the set of $\rho$-martingales $\mathcal{M}_{0}^{\rho}$
as its appearance in Corollary~\ref{P1}, Eqn.~\eqref{as00}
(and indirect appearance in Theorem~\ref{th:mulsub}, (ii))
confirms.
In our numerically implementable method developed and applied in Sections~\ref{sec:MSAL}--\ref{sec:num}
we rely on the $\rho$-Doob martingale.
From a theoretical perspective, however,
and as a general justification of our pathwise dual, martingale-based approach,
we develop in this section several results
on so-called \textit{surely optimal $\rho$-martingales}.
To achieve this, we generalize the concept of standard surely optimal martingales
(see Schoenmakers, Zhang and Huang \cite{SZH13}
in the context of standard conditional expectations and optimal single stopping problems)
to subadditive DMU functionals.
The results in this section show formally that
if a general $\rho$-martingale---not necessarily the $\rho$-Doob martingale---induces `small' (robust) variance,
then the associated bounds obtained from the dual representation can be expected to be `tight'
and nearly constant.

Our results on surely optimal $\rho$-martingales can also serve as a diagnostic device
to assess the quality of the estimated $\rho$-Doob martingale,
derived from an (input) approximation to the upper Snell envelope.
If the (robust) variance the estimate induces fails to be small, then it must be far from the $\rho$-Doob martingale.
If, on the other hand, this variance is small, then the estimate will be close to an optimal $\rho$-martingale
(attaining the dual representation's infimum),
even though not necessarily close to the $\rho$-Doob martingale.

For ease of exposition, we focus attention first on optimal single stopping problems.
The next theorem generalizes the analogous key measurability result for standard conditional expectations and optimal single stopping problems
to DMU functionals satisfying (C1)--(C4) and (P1).
\begin{theorem}
\label{th:sure1}
Let $Y_{i}^{\ast}$ be the upper Snell envelope of the cash-flow process $H$
with respect to a subadditive DMU functional $\rho$ satisfying (C1)--(C4) and (P1) as in Corollary~\ref{P1}
and let $M$ be a $\rho$-martingale.
Then, for any $i\in{\{0,\ldots,T\}}$,
\begin{equation*}
\max_{i\leq j\leq T}(H_{j}-M_{j}+M_{i})\in\mathcal{F}_{i}\Rightarrow
\max_{i\leq j\leq T}(H_{j}-M_{j}+M_{i})=Y_{i}^{\ast}.
\end{equation*}
\end{theorem}

The following lemma will later allow for a generalization of the results in this section to multiple stopping.
\begin{lemma}
\label{th:sure1_cor}
Let $Y^{\ast}$, $H$, $M$ and $\rho$ be as in Theorem~\ref{th:sure1}.
Then, for any fixed $0\leq i<T$,
\begin{align*}
\theta_{i+} &  :=\max_{i<j\leq T}(H_{j}-M_{j}+M_{i})\in\mathcal{F}
_{i}\Rightarrow\\
(i)\text{ \ }\theta_{i+} &  =\rho_{i}\left(  Y_{i+1}^{\ast}\right)  \text{
\ \ and \ \ }(ii)\text{ \ }M_{i+1}-M_{i}=Y_{i+1}^{\ast}-\rho_{i}\left(
Y_{i+1}^{\ast}\right),\text{ }
\end{align*}
hence $M_{i+1}-M_{i}$ is a $\rho$-Doob martingale increment.
(Note the strict first inequality under the max operator.)
Thus, in particular, if
$\theta_{i+}\in\mathcal{F}_{i}$ for every $0\leq i<T,$ then $M$ is the $\rho$-Doob martingale.
\end{lemma}

Let us define the conditional $\rho$-variance as follows:
\begin{equation}
\mathrm{Var}_{\rho_{i}}\left(  X\right)  :=\rho_{i}\left(  \left(  X-\rho
_{i}\left(  X\right)  \right)  ^{2}\right). \label{cv}
\end{equation}
It admits a conditional Chebyshev inequality,
exploited in the proof of Theorem \ref{th:Th14} below,
as follows:
\begin{proposition}
\label{prop:ch}
Assume (C1)--(C4).
If $\rho$ is positively homogeneous (P3),
then
\begin{equation}
\rho_{i}\left(  1_{\left\vert X-\rho_{i}\left(  X\right)  \right\vert
\geq\epsilon}\right)  \leq\frac{\mathrm{Var}_{\rho_{i}}\left(  X\right)
}{\epsilon^{2}}. \label{ch}
\end{equation}
\end{proposition}
Next, we state the following lemma:
\begin{lemma}
\label{sensv}
Assume (C1)--(C4).
Let $\rho$ be subadditive (P1) and sensitive (P2).
Then,
\begin{equation*}
\mathrm{Var}_{\rho_{i}}\left(  X\right)  =0\Longleftrightarrow\text{ }
X\in\mathcal{F}_{i}.
\end{equation*}
\end{lemma}

By virtue of Lemma~\ref{sensv}, Theorem~\ref{th:sure1} implies that if a $\rho$-martingale
$M$ is such that, for some $i\leq j\leq T$, the conditional $\rho$-variance
\begin{equation*}
\mathrm{Var}_{\rho_{i}}\left(  \theta_{i}(M)\right)  :=\mathrm{Var}_{\rho_{i}
}\left(  \max_{i\leq j\leq T}(H_{j}-M_{j}+M_{i})\right)
\end{equation*}
is zero a.s., then $\theta_{i}(M)=Y_{i}^{\ast}$ a.s.
In that case, we say that the $\rho$-martingale $M$ is \textit{surely optimal} at $i$.
(Note that, in particular, the $\rho$-Doob martingale in (\ref{eq:DoobMeyerdisc}) is surely optimal.)

We then present a stability result for $\rho$-martingales $M$ that are,
in loose terms, `close' to be surely optimal,
in the sense that the conditional $\rho$-variance $\mathrm{Var}_{\rho_{i}}(\theta_{i}(M))$ is `small'.
In particular, for a sequence of $\rho$-martingales $(M^{(n)})_{n\geq1}$
that induces vanishing conditional $\rho$-variance,
we establish weak conditions
guaranteeing that the corresponding upper bounds converge to the upper Snell
envelope (in $L^{1}$),
even though the sequence of $\rho$-martingales $(M^{(n)})_{n\geq1}$ itself does not necessarily converge.

\begin{theorem}
\label{th:Th14}
Assume (C1)--(C4).
Let $\rho$ be subadditive (P1) and positively homogeneous (P3).
Suppose that
\begin{equation*}
\mathrm{Var}_{\rho_{i}}\left(  \theta_{i}^{(n)}\right)  \overset{\text{P}}{\rightarrow}0,\quad
\mathrm{with}\quad
\theta_{i}^{(n)}=\max_{i\leq j\leq T}\left(  H_{j}-M_{j}^{(n)}+M_{i}^{(n)}\right).
\end{equation*}
If, in addition, for every $i$ and every $\epsilon>0$ there exists
$K_{\epsilon}>0$ such that
\begin{equation}
\sup_{n\geq1}\mathbb{E}\left[\rho_{i}\left(  \left\vert M_{i}^{(n)}\right\vert
1_{\left\vert M_{i}^{(n)}\right\vert >K_{\epsilon}}\right)\right]  <\epsilon,
\label{gUI}
\end{equation}
then
\begin{equation*}
\rho_{i}\left(  \theta_{i}^{(n)}\right)  \overset{L_{1}}{\longrightarrow}Y_{i}^{\ast}.
\end{equation*}
\end{theorem}

Note that, if $\rho_{i}\equiv\mathbb{E}_{i}$, \eqref{gUI} boils down to a standard
uniform integrability condition.
More generally, we have the following:

\begin{proposition}
\label{prop:EM}
Assume (C1)--(C4).
Let $\rho$ be subadditive (P1) and positively homogeneous (P3).
If, for some $\eta>0$,
\begin{equation*}
\sup_{n\geq1}\mathbb{E}\left[\rho_{i}\left(\left\vert M_{i}^{(n)}\right\vert
^{1+\eta}\right)\right]  <\infty,
\end{equation*}
then $\left(M_{i}^{(n)}\right)_{n\geq1}$ satisfies (\ref{gUI}).
\end{proposition}

Under an additional Lipschitz continuity condition,
Theorem~\ref{th:Th14} may be readily applied as follows.
Let us assume that, for some number $p$,
\begin{equation}
\mathbb{E}\left[\left\vert \rho_{i}\left(  Z\right)  \right\vert ^{p}\right]\leq
C_{p}\mathbb{E}\left[  \left\vert Z\right\vert ^{p}\right],  \label{l1}
\end{equation}
with $C_{p}>0.$
In particular, if (\ref{l1}) holds for $p=1$, one obviously has
\begin{equation}
\mathbb{E}\left[\mathrm{Var}_{\rho_{i}}\left(X\right)\right]:=\mathbb{E}\left[
\rho_{i}\left(\left(X-\rho_{i}\left(X\right)\right)^{2}\right)
\right]  \leq C_{1}\widetilde{\mathrm{Var}}_{i}\left(X\right),
\label{vc1}
\end{equation}
where

$\widetilde{\mathrm{Var}}_{i}(X):=\mathbb{E}_{\mathcal{F}_{i}}\left[  \left(
X-\rho_{i}\left(  X\right)  \right)  ^{2}\right]$.

That is, if we achieve in an algorithm that

$\widetilde{\mathrm{Var}}_{i}\left(\theta_{i}^{(n)}\right)
\overset{\text{P}}{\rightarrow}0$,

and the $\left(  M_{i}^{(n)}\right)_{n\geq1}$ are standard uniformly
integrable, i.e.,
\begin{equation*}
\sup_{n\geq1}\mathbb{E}\left[  \left\vert M_{i}^{(n)}\right\vert 1_{\left\vert
M_{i}^{(n)}\right\vert >K_{\epsilon}}\right]  <\epsilon,
\end{equation*}
then on the one hand,
\begin{align*}
\sup_{n\geq1}\mathbb{E}\left[\rho_{i}\left(  \left\vert M_{i}^{(n)}\right\vert
1_{\left\vert M_{i}^{(n)}\right\vert >K_{\epsilon}}\right)\right]   &  \leq C_{1}
\sup_{n\geq1}\mathbb{E}\left[\mathbb{E}_{i}\left\vert M_{i}^{(n)}\right\vert 1_{\left\vert
M_{i}^{(n)}\right\vert >K_{\epsilon}}\right]\\
&  =C_{1}\sup_{n\geq1}\mathbb{E}\left[\left\vert M_{i}^{(n)}\right\vert
1_{\left\vert M_{i}^{(n)}\right\vert >K_{\epsilon}}\right]<\epsilon,
\end{align*}
i.e., the $\left(  M_{i}^{(n)}\right)_{n\geq1}$ satisfy the notion of
$\rho_{i}$-uniform integrability,
and on the other hand we have due to (\ref{vc1}) that
\begin{equation*}
\mathrm{Var}_{\rho_{i}}\left(  \theta_{i}^{(n)}\right)  \overset{\text{P}
}{\rightarrow}0.
\end{equation*}
Then, Theorem~\ref{th:Th14} implies that
\begin{equation*}
\rho_{i}\left(  \theta_{i}^{(n)}\right)  \overset{L_{1}}{\rightarrow}
Y_{i}^{\ast}.
\end{equation*}

The next theorem generalizes Theorem~\ref{th:sure1} and Lemma~\ref{th:sure1_cor} to multiple stopping:
\begin{theorem}
\label{th:sure1_cor_mult}
Assume (C1)--(C4) and (P1).
Let us define for a set of $\rho$-martingales $M^{(k)}$, $k=1,\ldots,L$,
\begin{equation*}
\Theta_{i+}^{q}:=\max_{i<j_{1}<j_{2}<\cdots<j_{q}}\sum_{k=1}
^{q}\left(H_{j_{k}}+M_{j_{k-1}}^{(q-k+1)}-M_{j_{k}}^{(q-k+1)}\right)\text{ \ \ for
\ \ }q=1,\ldots,L,
\end{equation*}
with $j_{0}:=i.$
(Note the strict first inequality under the max operator.)
Then it holds that
\begin{align*}
\Theta_{i+}^{q}  &  \in\mathcal{F}_{i}\text{ \ \ for \ \ }q=1,\ldots,L,\text{
\ \ }0\leq i<T\text{ \ \ \ }\Longrightarrow\\
&\left\{
\begin{array}
[c]{l}
(i)\text{ \ \ }\Theta_{i+}^{q}=\rho_{i}\left(  Y_{i+1}^{\ast, q}\right) \\
(ii)\text{ \ }M_{i+1}^{(q)}-M_{i}^{(q)}=Y_{i+1}^{\ast, q}-\rho_{i}\left(
Y_{i+1}^{\ast, q}\right)
\end{array}
\right.  \text{for \ \ }q    =1,\ldots,L,\text{ \ \ }0\leq i<T.
\end{align*}
\end{theorem}

\begin{remark}
Without doubt it is also possible to derive a version of Theorem~\ref{th:Th14} for the multiple stopping setting.
However, as our algorithm in Section~\ref{sec:MSAL} below aims at approximative construction
of $\rho$-Doob martingale increments associated with
the upper Snell envelopes $Y^{\ast, l}$ of the
generalized cash-flows
\begin{equation}
U_{j}^{\ast, l}:=H_{j}+\rho_{j}\left(  Y_{j+1}^{\ast, l-1}\right)  ,\text{
\ \ }l=1,\ldots,L, \label{genc0}
\end{equation}
respectively, rather than approximative construction of merely surely optimal
$\rho$-martingales, we refrain from such an analysis.
\end{remark}

\setcounter{equation}{0}

\section{A General Primal-Dual Pseudo Algorithm}

\label{sec:MSAL}

In this section, we develop a primal-dual pseudo algorithm for robust multiple stopping (henceforth called \textit{algorithm} for short).
Our treatment in this section applies to DMUs satisfying (C1)--(C4), (P1) and weak continuity conditions,
and to general reward processes in a Markovian environment;
in particular, our treatment in this section is not restricted to $g$-expectations.
The following lemma will serve as a cornerstone in our construction.

\begin{lemma}
\label{var+}
Let $\rho$ satisfy (C1)--(C4), (P1)
and be Lipschitz continuous in the sense of (\ref{l1}) for $p=2$.
Furthermore, let $\mathcal{C}^{N},\mathcal{C},\mathcal{U}\in\mathcal{F}_{j}$,
$\mathcal{Y}\in\mathcal{F}_{j+1}$,
and let $\mathfrak{m}^{N}\in\mathcal{F}_{j+1}$ be a $\rho$-martingale increment,
that is, $\rho_{j}\left(\mathfrak{m}^{N}\right)=0,$ for $j=0,\ldots,T$,
$N\in\mathbb{N}$,
such that
\begin{align}
\mathbb{E}\left[  \left(  \mathcal{Y}-\mathfrak{m}^{N}-\mathcal{C}^{N}\right)
^{2}\right]   &  \rightarrow0,\text{ \ \ and}\label{coy}\\
\mathbb{E}\left[  \left(  \mathcal{C}^{N}-\mathcal{C}\right)  ^{2}\right]   &
\rightarrow0,\text{ \ \ for }N\rightarrow\infty.\nonumber
\end{align}
Then,
\begin{equation*}
\mathcal{C}=\rho_{j}\left(  \mathcal{Y}\right)  \text{ \ \ and \ \ }
\mathfrak{m}^{N}\overset{L_{2}}{\rightarrow}\mathcal{Y}-\rho_{j}\left(
\mathcal{Y}\right).
\end{equation*}
\end{lemma}

\begin{corollary}
\label{loc}
Let $Y_{j+1}^{l}$ be an already constructed approximation to a random variable $Y_{j+1}^{\ast, l}$.
Furthermore, let $\rho$,
the $\rho$-martingale increment $\mathfrak{m}_{j+1}^{l,N}\in\mathcal{F}_{j+1}$,
and $\mathcal{C}_{j}^{l,N},\mathcal{C}_{j}^{l}\in\mathcal{F}_{j}$ be as in Lemma~\ref{var+},
such that
\begin{align*}
\mathbb{E}\left[  \left(  Y_{j+1}^{l}\,-\mathfrak{m}_{j+1}^{l,N}
-\mathcal{C}_{j}^{l,N}\right)  ^{2}\right]\rightarrow0,
\quad\text{and}\quad
\mathbb{E}\left[  \left(  \mathcal{C}_{j}^{l,N}-\mathcal{C}_{j}^{l}\right)
^{2}\right]\rightarrow0,
\quad\text{for}\ N\rightarrow\infty.
\end{align*}
Then,
\begin{equation*}
Y_{j+1}^{l}\,-\mathfrak{m}_{j+1}^{l,N}\overset{L_{2}}{\rightarrow}\rho
_{j}\left(  Y_{j+1}^{l}\right).
\end{equation*}
\end{corollary}

Guided by Lemma~\ref{var+} and Corollary~\ref{loc}, we now develop
a primal-dual algorithm in the context of a Markovian underlying process $X$
with state space $\mathbb{R}^{d}$, possibly in continuous time, that is
monitored at the exercise dates as $X_{j},$ $j=0,\ldots,T$.
As usual, we assume that $\mathcal{F}_{j}$ is the $\sigma$-field generated
(directly or, as in the next section, indirectly)
by the process
$X$ up to exercise date $j$.
Furthermore, we assume that the cash-flows are of the form
\begin{equation*}
H_{j}=f_{j}(X_{j}),\text{ \ \ for \ \ }j=0,\ldots,T,
\end{equation*}
where $f_{j}:\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq0}$, $j=0,\ldots,T$, are
given nonnegative payoff functions
such that $H\in\mathcal{H}$.
Note that, due to the Bellman principle (\ref{eq:Bellman}),
$Y^{\ast, l}$ can be seen as the upper Snell envelope corresponding
to the generalized cash-flow
\begin{equation}
U_{j}^{\ast, l}:=H_{j}+\rho_{j}\left(Y_{j+1}^{\ast, l-1}\right)
=:\ H_{j}+c_{j}^{\ast, l}(X_{j}),\text{\ \ }l=1,\ldots,L,\quad j=0,\ldots,T-1,
\label{genc}
\end{equation}
due to a single exercise right,
where the so-called continuation functions $c_{j}^{\ast, l}$ exist by
Markovianity.

We also assume to have a set of Monte Carlo simulated training
trajectories
\begin{equation*}
X^{(n)}\equiv X^{n},\qquad n=1,\ldots,N.
\end{equation*}

We proceed in the following steps:
\begin{itemize}
\item[1.] Initialize $\overline{M}^{0}$ $=$ $\overline{Y}^{0}$
$=\overline{c}^{0}$ $=$ $0$, for $l=0$.
\item[2.] Suppose that, for a particular $l$ with $0\leq l<L$:\newline
(i) we have constructed a set of (approximate) continuation functions
$\overline{c}_{j}^{l}:\mathbb{R}^{d}\rightarrow\mathbb{R}_{\geq0},$ $1\leq j\leq T$,
hence an (approximate) continuation value process (for up to $l$ exercise times)
of the form
\begin{equation*}
\overline{C}_{j}^{l}=\overline{c}_{j}^{l}(X_{j});
\end{equation*}
(ii) we have constructed a (true) $\rho$-martingale $\overline{M}_{j}^{l}$;
and \newline
(iii) we have constructed, on each trajectory $n$, a path
\begin{align}

\overline{Y}_{j}^{l,n}\,  &  :=\max\left[ \overline{U}_{j}^{l,n}, \overline{c}_{j}^{l}(X_j^n) \right]
 \mbox{ if } j \in \{0,1,\ldots,T\}, \qquad\text{\ where }  \label{bar}\\
\overline{U}_{j}^{l,n}  &  :=\left\{
\begin{array}
[c]{l}
f_{j}(X_{j}^{n})+\overline{c}_{j}^{l-1}(X_{j}^{n})\text{ \ \ if
\ \ \ }l>0\\
0\text{ \ \ \ if \ \ }l=0
\end{array}
\right.  ,\qquad 0\leq j\leq T,

\end{align}
as an approximation to $Y^{\ast, l,n}$.
\item[3.] Now construct, using these trajectories,
a subsequent (true) $\rho$-martingale $\overline{M}^{l+1}$,
a subsequent set of continuation functions $\overline{c}_{j}^{l+1}$, $j=0,\ldots,T$,
and subsequent trajectories $\overline{Y}^{l+1,n}$, $n=1,\ldots,N$
(as approximations to $M^{\ast, l+1}$ and $Y^{\ast, l+1}$, respectively)
such that (\ref{bar}) holds for $l+1$.
To this end, we carry out the
following backward procedure, or \textquotedblleft backward
subroutine\textquotedblright, at level $l+1$:
\begin{itemize}
\item As initialization, set
$\overline{Y}_{T}^{l+1}=H_{T},$ $\overline{c}_{T}^{l+1}=0$.
(We also set $\overline{c}^0_T = 0$.)
\item Suppose that, for $0<j+1\leq T$,
the values $\overline{Y}_{j+1}^{l+1,n}$, $n=1,\ldots,N$,
the set of $\rho$-martingale increments $\left(\overline{M}_{r}^{l+1}-\overline{M}_{j+1}^{l+1}\right)_{j+1<r\leq T}$
(which is empty if $j+1=T$), and the continuation function $\overline{c}_{j+1}^{l+1}$ have been
constructed.
\item Then construct, according to the regression subroutine in Section~\ref{sec:PD} below,
a continuation function $\overline{c}_{j}^{l+1}$,
a $\rho$-martingale increment $\overline{\mathfrak{m}}_{j+1}^{l+1}$ $\in\mathcal{F}_{j+1}$
with $\rho_{j}(\overline{\mathfrak{m}}_{j+1}^{l+1})=0$,
and set $\left(\overline{M}_{r}^{l+1}-\overline{M}_{j}^{l+1}\right)

=\left(\overline{M}_{r}^{l+1}-\overline{M}_{j+1}^{l+1}
+\overline{\mathfrak{m}}_{j+1}^{l+1}\right)$,

for
\begin{align}
\overline{Y}_{j}^{l+1,n}  &  =\max\left[  \overline{U}_{j}^{l+1,n},
\overline{c}_{j}^{l+1}(X_j^n) \right], \label{algorecursion} \text{ \ \ with}\\
\overline{U}_{j}^{l+1,n}  &  =f_{j}(X_{j}^{n})
+\overline{c}_{j}^{l}(X_{j}^{n}),\text{ \ \ and \ \ }n=1,\ldots,N. \nonumber
\end{align}
\end{itemize}
\end{itemize}

Proceeding this way:
\begin{itemize}
\item[(a.)] Working forward from $l=0,\ldots,L$ thus yields
a family of continuation functions $\overline{c}^{l}$ and
a family of (true) $\rho$-martingales $\overline{M}^{l}$,
respectively:
\begin{equation*}
\overline{c}_{j}^{l}(\cdot),\text{ \ \ and \ \ }
\overline{M}_{j}^{l}:=\sum_{r=1}^{j}\overline{\mathfrak{m}}_{r}^{l},\text{ \ \ \ }l=1,\ldots
,L,\text{ }j=0,\ldots,T.
\end{equation*}

\item[(b.)] An upper bound for the solution to the robust multiple stopping problem
at $t=0$ due to $L$ exercise rights
is now given by (cf. Theorem~\ref{th:mulsub}, (ii)):
\begin{equation}
Y_{0}^{\text{upp},L}:=\rho_{0}\left(  \max_{0\leq j_{1}<j_{2}<\cdots< j_{L}}
\sum_{l=1}^{L}\left(  f_{j_{l}}(X_{j_{l}})-\overline{M}_{j_{l}}^{L-l+1}
+\overline{M}_{j_{l-1}}^{L-l+1}\right)  \right),
\label{eq:upp}
\end{equation}
which needs to be estimated by a separate (Monte Carlo) procedure.

\item[(c.)] A lower bound for the solution to the robust multiple stopping problem
at $t=0$ due to $L$ exercise rights
may next be obtained from the family of stopping times
\begin{equation}
\label{taul}\tau^{l}:=\min\left\{  j:\tau^{l-1}<j\leq T,\text{ }f_{j}
(X_{j})+\overline{c}_{j}^{l-1}(X_{j})\geq\overline{c}_{j}^{l}(X_{j})\right\},
\end{equation}
via a (Monte Carlo) estimation of:
\begin{equation}
\label{eq:low}
Y_{0}^{\text{low},L}:=\rho_{0}\left(  \sum_{l=1}^{L}\text{ }f_{\tau^{l}}
(X_{\tau^{l}})\right).
\end{equation}
\end{itemize}

\subsection{Regression Subroutine}

\label{sec:PD}

Let there be given a collection of
`elementary'
$\rho$-martingale increments
$\mathcal{E}_{j}^\beta(X)$,
with $\beta=(\beta_1,\ldots,\beta_{K^{\prime}}) \in \mathbb{R}^{K^{\prime}}$ and ${K^{\prime}}\in \mathbb{N}$,
that is,
\begin{equation*}
\rho_{j}\left(  \mathcal{E}_{j}^\beta(X) \right)=0,\qquad j=0,\ldots,T-1,
\end{equation*}
and a collection of basis functions
$\psi_{1},\ldots,\psi_{K^{\prime\prime}}:\mathbb{R}^{d}\rightarrow\mathbb{R}$.
We assume that the set of $\rho$-martingale increments
\begin{equation*}
\bigg\{ \mathcal{E}^{(\beta_1,\ldots,\beta_{K^{\prime}})}_j(X)  \bigg\},
\end{equation*}
is $L^2$-dense among the $\mathcal{F}_{j+1}$-measurable square-integrable random variables $\mathcal{E}_j$
such that $\rho_{j}(\mathcal{E}_j)=0$.
We then solve, in view of Lemma~\ref{var+} and Corollary~\ref{loc},
for fixed $N$ and $K^{\prime}$ and $K^{\prime\prime}$ the least squares problem
\begin{align}
\text{MSE} &:= \sum_{n=1}^{N}\left(  Y_{j+1}^{l+1,n}-\mathcal{E}_{j}^\beta(X)-\sum_{k=1}^{K^{\prime\prime}}\gamma_{k}\psi_{k}(X_{j}^{n})\right)
^{2} = \sum_{n=1}^{N} \left(  Y_{j+1}^{l+1,n}-\mathcal{E}
_{j}^\beta(X)-\gamma\psi(X_{j}^{n})\right)
^{2}\nonumber\\
&\longrightarrow \argmin_{\beta\in\mathbb{R}^{K^{\prime}},\gamma\in\mathbb{R}^{K^{\prime\prime}}}=:\left[  \beta
^{l+1,j,K^{\prime},N},\gamma^{l+1,j,K^{\prime\prime},N}\right], \label{argmin}
\end{align}
where we used vector notation $\gamma=(\gamma_1,\ldots,\gamma_{K^{\prime\prime}})$
and $\psi=(\psi_1,\ldots,\psi_{K^{\prime\prime}})^{\transpose}$.
For the algorithm to converge, it is actually sufficient that the MSE above
converges to zero as $K^{\prime},K^{\prime\prime} \rightarrow \infty$ for our choice of $\beta$ and $\gamma$.
We will suppress the superscripts $K^{\prime},K^{\prime\prime}$ and $N$ whenever there is no ambiguity.
We set
\begin{equation*}
\overline{\mathfrak{m}}_{j+1}^{l+1}:=\overline{\mathfrak{m}}_{j+1}^{l+1}(X):=
\mathcal{E}_{j}^{\beta^{l+1,j}}(X),\text{ \ \ \ }\overline{c}_{j}^{l+1}
(\cdot):=\sum_{k=1}^{K^{\prime\prime}}\gamma_{k}^{l+1,j}\psi_{k}(\cdot).
\end{equation*}

\subsection{Convergence Theorem}

We state the following theorem:

\begin{theorem}
\label{theoConv}
Let $\rho$ be subadditive (P1)
and Lipschitz continuous in the sense of (\ref{l1}) for $p=2$.
We set $K=\min{\{K^{\prime},K^{\prime\prime}\}}$ and denote by
$\overline{M}_{j}^{l,K,N}:=\overline{M}_{j}^{l,K,N}(X_j):=\overline{M}_{j}^{l}$,
$\overline{c}_{j}^{l,K,N}:=\overline{c}_{j}^{l,K,N}(X_j):=\overline{c}_{j}^{l}$ and
$\overline{Y}_{j}^{l,K,N}:=\overline{Y}_{j}^{l,K,N}(X_j):=\overline{Y}_{j}^{l,N}$
the functions constructed in the algorithm above.
Then,
\begin{equation}
\lim_{K\rightarrow \infty} \lim_{N\rightarrow \infty} \overline{M}_j^{l,K,N} = M_j^{\ast, l} \mbox{ in } L^2, \label{Mconv}
\end{equation}
\begin{equation}
\lim_{K\rightarrow \infty} \lim_{N\rightarrow \infty} \overline{c}_j^{l,K,N} = c_j^{\ast, l} \mbox{ in } L^2, \label{Cconv}
\end{equation}
\begin{equation}
\lim_{K\rightarrow \infty} \lim_{N\rightarrow \infty} \overline{Y}_j^{l,K,N} = Y_j^{\ast, l} \mbox{ in } L^2, \label{Yconv}
\end{equation}
for all $j= T, T-1, \ldots, 0$ and $l= 1,\ldots,L$.
Furthermore,
\begin{align*}
\left\vert Y_{j}^{\ast, l} \right. \left. -\rho_{j}\left(  \max_{j\leq r\leq T}\left(
\overline{U}_{r}^{l}-\overline{M}_{r}^{l}\right)  \right)  \right\vert
&\leq
\rho_{j}\left(  \max_{j\leq r\leq T}\left\vert M_{r}^{\ast, l}-\overline{M}_{r}^{l}\right\vert \right)
+
\rho_{j}\left(  \max_{j\leq r\leq T}\,\left\vert \overline{c}_{r}^{l-1}-c_{r}^{\ast, l-1}\right\vert \right) \\
&\longrightarrow_{K\rightarrow \infty, N\rightarrow \infty} 0.
\end{align*}
\end{theorem}

\subsection{Complexity}

At first sight, the path-wise maximum in \eqref{eq:upp} would require the evaluation
of $T!/(L!(T-L)!)$ terms.
Fortunately, due to following proposition, it only requires $O(LT)$ evaluations.

\begin{proposition}
Define, for $1\leq q\leq L$ and $0\leq i\leq T,$
\begin{equation*}
\Theta_{i}^{q}:=\max_{i\leq j_{1}<j_{2}<\cdots< j_{q}}\sum_{l=1}^{q}\left(
f_{j_{l}}(X_{j_{l}})-\overline{M}_{j_{l}}^{q-l+1}+\overline{M}_{j_{l-1}
}^{q-l+1}\right)  , \text{ \ \ with \ }j_{0}=i,
\end{equation*}
and naturally $\Theta_{i}^{q}=0$, $i>T$.
Then,
\begin{equation}
\Theta_{i}^{q}=\max\left[  f_{i}(X_{i})+\overline{M}_{i}^{q-1}-\overline
{M}_{i+1}^{q-1}+\Theta_{i+1}^{q-1},\overline{M}_{i}^{q}-\overline{M}_{i+1}
^{q}+\Theta_{i+1}^{q}\right]. \label{rec}
\end{equation}
\label{prop:com}
\end{proposition}

Thus, the evaluation of \eqref{eq:upp} may be described as follows:

\paragraph{Recursive evaluation of \eqref{eq:upp}}

\begin{itemize}
\item[1.] Initialize $\Theta_{i}^{0}=0,$ for $i=0,\ldots,T;$

\item[2.] Suppose that, for $0\le q-1<L$ and for $i=0,\ldots,T$, the construction
of $\Theta_{i}^{q-1}$ has been conducted;

\item[3.] Backward subroutine: Initialize $\Theta_{T}^{q}=f_{T}(X_{T})$.
When $\Theta_{i+1}^{q}$ has been constructed for $i+1\leq T$,
compute $\Theta_{i}^{q}$ via (\ref{rec}).
\end{itemize}

\setcounter{equation}{0}

\section{Explicit
Construction
in a Brownian-Poisson Filtration}

\label{sec:BPcon}

In the sequel, we assume that we have a completed continuous-time filtration $\mathbb{F}=(\mathcal{F}_t)_{t\in[0,T]}$
on a filtered probability space $(\Omega, \mathcal{F}, \mathbb{F}, \mathbb{P})$
generated by a $d_1$-dimensional standard (i.e., zero mean and unit variance)
Brownian motion $W=(W_1,\ldots,W_{d_1})^{\transpose}$
and a $d_2$-dimensional Poisson process $N=(N_1,\ldots,N_{d_2})^{\transpose}$
with arrival intensity $\lambda_{\mathbb{P}}=(\lambda_{\mathbb{P}}^{1},\ldots,\lambda_{\mathbb{P}}^{d_{2}})^{\transpose}$.
As usual, we define the compensated counterpart of $N$ as $\tilde{N}_{t}=N_{t}-\lambda_{\mathbb{P}}t$.

The components of the processes $W$ and $N$ are assumed to be independent.
The stochastic drivers $W$ and $N$ generate
the underlying Markovian adapted reward process $(X_{t})_{t\in[0,T]}$ with state space $\mathbb{R}^{d}$
of Section~\ref{sec:MSAL}.

Furthermore, we assume that $\rho$ satisfies (the continuous-time analogs of) (C1)--(C4) and (P1)--(P3).
This means, in particular, that $\rho$ is a coherent risk measure.
By classical duality results (e.g., F\"ollmer and Schied \cite{FS11}),
the robust multiple stopping problem at time $t$ is then given by
\begin{align}
\label{eq:prob}
Y_{t}^{\ast,L}= \sup_{\substack{t\leq\tau_{1}<\cdots<\tau_{L} \\
(\tau_{1},\ldots,\tau_{L})\in\mathcal{T}_{t}(L)}}
\rho_{t}\left(  \sum_{l=1}^{L}H_{\tau_{l}}\right)
=\sup_{\substack{t\leq\tau_{1}<\cdots<\tau_{L} \\
(\tau_{1},\ldots,\tau_{L})\in\mathcal{T}_{t}(L)}} \sup_{\mathbb{Q}\in\mathcal{Q}}
\mathbb{E}_{\mathbb{Q}}\left[\sum_{l=1}^{L}H_{\tau_{l}}\Big|\mathcal{F}_{t}\right],
\quad 0\leq t\leq T,\end{align}
with $\mathcal{T}_{t}(L)$ our family of stopping vectors,
and $\mathcal{Q}$ a closed convex set of probability measures absolutely continuous with respect to $\mathbb{P}$
and satisfying a stability assumption.
In such a continuous-time setting, it is known that every recursive coherent risk measure
can be identified with a solution to a backwards stochastic differential equation (BSDE) also called a $g$-expectation,
modulo a compactness assumption;
see Section \ref{sec:tcdp} for the precise definitions and results.
Exploiting our algorithm presented in Section \ref{sec:MSAL},
this section constructs explicit upper and lower bounds to $Y^{\ast,L}$
with desirable properties.

\subsection{Bellman's Principle, the Set of Priors, and BSDE drivers}

\label{sec:tcdp}

A probability measure change from $\mathbb{P}$ to an absolutely continuous measure $\mathbb{Q}\in\mathcal{Q}$
admits an explicit representation in our Brownian-Poisson setting.
Consider the Radon-Nikodym derivative
\begin{equation*}
D_{t}:=\mathbb{E}\left[\frac{d\mathbb{Q}}{d\mathbb{P}}|\mathcal{F}_{t}\right],\qquad t\in[0,T].
\end{equation*}
As is well-known, $D_{t}$ has the Dol\'eans-Dade exponential form
\begin{equation}
\label{eq:Radon}
D_{t}=\exp\left(\int_{0}^{t}q_{s}dW_{s}
+ \int_{0}^{t}\log\left(\frac{\lambda_{s}}{\lambda_{\mathbb{P}}}\right)dN_{s}
-\int_{0}^{t}\left(\frac{|q_{s}|^{2}}{2} + \lambda_{s}-\lambda_{\mathbb{P}}\right)ds\right),
\end{equation}
where $q$ is a predictable, $\mathbb{R}^{d_{1}}$-valued, stochastic drift,
and $\lambda$ is a positive, predictable, $\mathbb{R}^{d_{2}}$-valued process
with $\frac{\lambda_{s}}{\lambda_{\mathbb{P}}}:=(\frac{\lambda^{1}_{s}}{\lambda^{1}_{\mathbb{P}}},
\ldots,\frac{\lambda^{d_{2}}_{s}}{\lambda^{d_{2}}_{\mathbb{P}}})^{\transpose}$,
which jointly uniquely characterize $\mathbb{Q}$.
From Girsanov's theorem, we know that
$W^{\mathbb{Q}}_{t}:=W_{t}-\int_{0}^{t} q_{s} ds$ is a standard Brownian motion under $\mathbb{Q}$
while the process ${N}_{t}$ has arrival intensity $\lambda_{t}$.
In particular, the reference model $\mathbb{P}$ corresponds to $q\equiv 0$ and $\lambda\equiv \lambda_{\mathbb{P}}$.
The stochastic drift $q$ may be given the interpretation of a drift in the diffusive component
that is misspecified to be absent by the reference model $\mathbb{P}$.
Similarly, $\lambda_{s}-\lambda_{\mathbb{P}}$ represents a deviation
from the misspecified arrival intensity $\lambda_{\mathbb{P}}$ under $\mathbb{P}$.

In our dynamic setting,
with DMU evaluations satisfying the continuous-time analogs of (C1)--(C4) and (P1)--(P3),
time-consistency of choice under uncertainty is satisfied
as it is equivalent to recursiveness
or Bellman's dynamic programming principle.
Time-consistency of a dynamic evaluation $(\rho_{t}(H))_{t\in[0,T]}$ requires---according to its usual definition,
also referred to as `strong' time-consistency---that $\rho_{s}(H_{1}) \ge\rho_{s}(H_{2})$
whenever $\rho_{t}(H_{1}) \ge\rho_{t}(H_{2})$, $t\geq s$.
That is, if $H_{2}$ is preferred over $H_{1}$,
in each state of nature at time $t$,
then the same preference necessarily applies prior to time $t$;
see e.g.,
Riedel \cite{R04}, Ruszczy\'nski and Shapiro \cite{RS06b},
Shapiro, Dentcheva and Ruszczy\'nski \cite{SDR09}, Chapter 6,
Ruszczy\'nski \cite{R10} and Shapiro \cite{S16}.
Indeed, requiring recursiveness or Bellman's dynamic programming principle
is equivalent to requiring time-consistency for
$\rho_{t}(H)=\sup_{\{\mathbb{Q}\sim \mathbb{P}|\mathbb{Q}=\mathbb{P}\mbox{ on }\mathcal{F}_{t}\}}
\mathbb{E}_{\mathbb{Q}}[H|\mathcal{F}_{t}]$, $t\in[0,T]$,
which is, in turn, equivalent to the set of priors $\mathcal{Q}$ being $m$-stable; see Delbaen \cite{D06}.
More formally, the following statements are equivalent
(see Lemma 11.11 of F\"ollmer and Schied \cite{FS11} for the equivalence (i)--(ii),
Delbaen \cite{D06} and Delbaen, Peng and Rosazza Gianin \cite{DPR10} for (ii)--(iii) in a Brownian setting,
and Tang and Wei \cite{TW12} and Laeven and Stadje \cite{LS14} for (ii)--(iii) in a general semi-martingale setting):

\begin{itemize}
\item[(i)] $\rho$ is recursive, i.e., $\rho$ satisfies Bellman's dynamic
programming principle $\rho_{0}(\rho_{t}(H)I_{A})=\rho_{0}(H I_{A})$
for every $t\in[0,T]$, $A\in\mathcal{F}_{t}$, and bounded $H$.

\item[(ii)] $\rho$ is time-consistent over bounded rewards.

\item[(iii)] There exists a closed, convex, set-valued predictable mapping $C$
taking values in $\mathbb{R}^{d_{1}}\times(-\lambda_{\mathbb{P}}^{1},\infty
)\times\cdots\times(-\lambda_{\mathbb{P}}^{d_{2}},\infty)$
such that
\begin{equation*}
\rho_{t}(H)=\sup_{(q,\lambda)\in C}\mathbb{E}_{\mathbb{Q}}\left[H|\mathcal{F}_{t}\right],\qquad t\in[0,T].
\end{equation*}
\end{itemize}

\vskip -0.2cm As in our continuous-time Brownian-Poisson setting recursiveness (C3) is equivalent to (iii),
we assume henceforth:
\begin{itemize}
\item[(A1)] $C=(C_{t})_{t\in[0,T]}\subset[0,T]\times\mathbb{R}^{d_{1}}
\times(-\lambda_{\mathbb{P}}^{1}+\varepsilon,\infty)\times\cdots\times(-\lambda_{\mathbb{P}}^{d_{2}}+\varepsilon,\infty)$
with $\varepsilon> 0$ is compact.

\end{itemize}
We note that (A1) also implies that $\rho$ is recursive and time-consistent over square-integrable rewards.

For $t\in[0,T]$, $z\in\mathbb{R}^{1\times d_{1}}$ and $\tilde{z}\in\mathbb{R}^{1\times d_{2}}$,
and $C$ satisfying Assumption (A1),
let us define a function $g$ via Fenchel's duality:
\begin{equation}
\label{g1}
g(t,z,\tilde{z}):=\sup_{(q,\lambda-{\lambda}_{\mathbb{P}})\in C_{t}}
\{zq+\tilde{z} (\lambda-\lambda_{\mathbb{P}})\}.
\end{equation}
One easily verifies that $g$ is convex, positively homogeneous and Lipschitz continuous.

Then, from Kr\"atschmer \textit{et al.} \cite{KLLSS18}, we have the following statement,
which is essentially (with (A1)) equivalent to (i)--(iii) above:

\begin{itemize}
\item[(iv)] For every $H\in L^{2}(\mathcal{F}_{j+1})$,
there exists a unique square-integrable predictable $(Z,\tilde{Z})$ such that
\begin{align*}
d\rho_{t}(H_{j+1})=-g(t,Z_{t},\tilde{Z}_{t})dt
+Z_{t}dW_{t} + \tilde{Z}_{t}d\tilde{N}_{t},\ \mbox{for
		}\,t\in[j,j+1],\ j\in\{0,\ldots,T-1\}.
\end{align*}
\end{itemize}
In particular, there exists a unique square-integrable predictable
$(Z^{\ast},\tilde{Z}^{\ast})$ such that
\begin{equation}
\label{eq:BSDE2}
d\rho_{t}(Y^{\ast}_{j+1})=-g(t,Z^{\ast}_{t},\tilde{Z}^{\ast}_{t})dt
+Z^{\ast}_{t}dW_{t} + \tilde{Z}^{\ast}_{t}d\tilde{N}_{t},\ \mbox{for
        }\,t\in[j,j+1],\ j\in\{0,\ldots,T-1\}.
\end{equation}
Furthermore, for $t\in[0,T]$, the $(Z^{\ast},\tilde{Z}^{\ast})$ in \eqref{eq:BSDE2}
recover---and later allow to practically compute---the $\rho$-Doob martingale as follows
(cf. Eqns.~\eqref{eq:DoobMeyerdisc}--\eqref{eq:DD}):
\begin{equation}
\label{sttar}M^{\ast}_{t} =\rho_{t}(M^{\ast}_{T})=-\int_{0}^{t} g(s,Z^{\ast}_{s},
\tilde{Z}^{\ast}_{s})ds + \int_{0}^{t} Z^{\ast}_{s}\,dW_{s}+ \int_{0}^{t}\tilde
{Z}^{\ast}_{s}\,d\tilde{N}_{s}.
\end{equation}
Here, $j=0,\ldots,L-1$ should be interpreted as exercise dates
and $t\in[j,j+1]$ as the continuous embedding.

Because $\rho_{j+1}(Y^{\ast}_{j+1})=Y^{\ast}_{j+1}$, by (iv), for $t\in[j,j+1]$,
\begin{equation}
\label{ruck}
\rho_{t}(Y^{\ast}_{j+1})=Y^{\ast}_{j+1}+\int_{t}^{j+1}g(s,Z^{\ast}_{s},\tilde{Z}^{\ast}_{s})ds
-\int_{t}^{j+1}Z^{\ast}_{s}dW_{s}
-\int_{t}^{j+1}\tilde{Z}^{\ast}_{s}d\tilde{N}_{s}.
\end{equation}

Eqn. (\ref{eq:BSDE2}) is referred to as a backward stochastic differential equation (BSDE).
Formally, given a terminal payoff $H\in L^{2}$
and a function $g:[0,T]\times\mathbb{R}^{d_{1}}\times\mathbb{R}^{d_{2}}\to\mathbb{R}$,
referred to as a driver, the solution to the corresponding BSDE is a
triple of square-integrable and suitably measurable processes $(Y,Z,\tilde{Z})$
that satisfies
\begin{equation*}
dY_{t}=-g(t,Z_{t},\tilde{Z}_{t})dt+Z_{t} dW_{t}+\tilde{Z}_{t}d\tilde{N}
_{t},\quad\mbox{and}\quad Y_{T}=H.
\end{equation*}
The solution is often referred to as a (conditional) $g$-expectation;
see, e.g., Peng \cite{P04}.

As a means of illustrating the generality of our setup given by (\ref{eq:prob}) with (A1)
we provide a few examples.

\begin{example}
\label{g-ex}
\begin{itemize}
\item[(1.)] Ball scenarios:
Consider a decision-maker endowed with a set of priors
constituting a small ball environment surrounding $\mathbb{P}$
all deemed equally plausible.
Then,
\begin{align*}
\mathcal{Q}=\Big\{\mathbb{Q}^{(q,\lambda)} \ll \mathbb{P} \Big||q_{t}|\leq\delta_{1},\quad|\lambda_{t}-{\lambda}_{\mathbb{P}}|\leq\delta_{2},
\mbox{ for Lebesgue-a.s. all }t\Big\},\qquad \delta_{1},\delta_{2}>0,
\end{align*}
and

$C_{t}= \big\{(q,\lambda) \big{|} |q|\leq\delta_{1}, |\lambda-{\lambda}_{\mathbb{P}}|\leq\delta_{2}\big\}$.

Suppose without losing generality that $|{\lambda}_{\mathbb{P}}|\geq\delta_{2}$.
Then, from \eqref{g1}, in explicit form,

$g(t,z,\tilde{z})=\delta_{1}|z| + \delta_{2}|\tilde{z}|$.

\item[(2.)] Discrete scenarios:
Imagine a decision-maker who considers,
at each time $t>0$,
finite-dimensional families $\{q_{1,t},\ldots,q_{m,t}\}$ and $\{\lambda_{1,t},\ldots,\lambda_{m,t}\}$,
$m\in\mathbb{N}$, with all elements deemed equally plausible.
Then,
\begin{align*}
\mathcal{Q}=\Big\{\mathbb{Q}^{(q,\lambda)} \ll \mathbb{P} \Big|(q_{t},\lambda_{t})\in\{(q_{i,t},\lambda_{j,t}), i,j\in\{1,\ldots,m\}\},
\mbox{ for Lebesgue-a.s. all }t\Big\},
\end{align*}
and
$C_{t}= \big\{(q,\lambda) \big{|} (q,\lambda)\in\mathrm{conv}\left(\{(q_{i,t},\lambda_{j,t}), i,j\in\{1,\ldots,m\}\}\right)\big\}$,
with $\mathrm{conv}(\cdot)$ the convex hull.
We can assume that $0\in\mathrm{conv}\left(\{(q_{i,t},\lambda_{j,t}), i,j\in\{1,\ldots,m\}\}\right)$
without losing generality, upon redefining the reference measure.
Furthermore,

$g(t,z,\tilde{z})=\max_{i=1,\ldots,m}q_{i,t}z+\max_{j=1,\ldots,m}\lambda_{j,t}\tilde{z}$.

\end{itemize}
\end{example}

To obtain genuine upper and lower bounds to the optimal solution of the stopping problem \eqref{eq:prob},
we henceforth impose the following additional assumption:

\begin{itemize}
\item[(A2)] $H_{t} = f_{t}(X_{t})$ and we can simulate i.i.d. copies of
$(X_{t})_{t\in[0,T]}$.

\end{itemize}

Once, we have constructed a `good' family of $\rho$-martingales $\overline{M}$,
we are faced with the computation of $Y_{0}^{\text{upp},L}$ in \eqref{eq:upp}.

An important advantage of this dual approach for numerical stability is that

we have, in fact, the pathwise dual representation
\begin{equation*}
Y_{0}^{\ast,L}:= \max_{0\leq j_{1}<j_{2}<\cdots< j_{L}}\sum_{l=1}^{L}\left(
f_{j_{l}}(X_{j_{l}})-M_{j_{l}}^{\ast,L-l+1}+M_{j_{l-1}}^{\ast,L-l+1}\right)
,\qquad\text{almost surely}.
\end{equation*}

Indeed, we will obtain an estimate with `low' variance
provided the $\rho$-martingale $\overline{M}$ is `good'.
Furthermore, to obtain a lower bound we employ $Y_{0}^{\text{low},L}$ in \eqref{eq:low}.

By the results discussed in this subsection, for a square-integrable payoff $U$, $\rho$ has a representation of the form
\begin{align}
\rho_{t}(U)  &  = \sup_{(q,\lambda) \in C} \mathbb{E}_{\mathbb{Q}^{(q,\lambda)}}\left[
U \big|\mathcal{F}_{t}\right] \label{bs2}\\
&  =U+\int_{t}^{T}g(s,\mathcal{Z}_{s},\mathcal{\tilde{Z}}_{s})ds-\int_{t}
^{T}\mathcal{Z}_{s}dW_{s} - \int_{t}^{T}\mathcal{\tilde{Z}}_{s}d\tilde{N}_{s}.
\label{bs}
\end{align}
Let us now first consider the question of how to explicitly obtain
a `good' family of $\rho$-martingales $\overline{M}$
in our Brownian-Poisson filtration with (C1)--(C4), (P1)--(P3) and (A1)--(A2).
Next, we develop in Sections~\ref{sec:LB} and~\ref{sec:UB} explicit genuine
lower and upper bounds.

\subsection{Parameterization of the \texorpdfstring{$\rho$-}{rho-}Martingale Increments}

\label{sec:algorithm}

As in Section \ref{sec:MSAL},
we assume that $(X_{j})_{0\leq j\leq T}$ is an $\mathcal{F}_{j}$-adapted
$d$-dimensional underlying Markovian process, now in a Brownian-Poisson filtration,
and that our cash-flow process has a structure of the form
$H_{j}=f_{j}(X_{j})$, $j=0,\ldots,T$ such that $H\in\mathcal{H}$.
We further assume that the resulting random variables $H_{j}$ are square integrable.

We are going to construct a $\rho$-martingale backwardly.
To this end, we consider, between two exercise dates $j$ and $j+1$,
the (fine) grid $\pi_{j}=\{s_{(j-1)n_{0}}=j,s_{j1},\ldots,s_{jn_{0}}=j+1\}$,
where $s_{jp}=j+p\Delta$ with $\Delta=n_{0}^{-1}$.
We also define
\begin{equation*}
\Pi:=\{s_{00}=0,s_{01},\ldots,s_{0n_{0}}=1,s_{11},\ldots,
s_{(T-2)n_{0}}=T-1,s_{(T-1)1},\ldots,s_{(T-1)n_{0}}=T\},
\end{equation*}
and sometimes use the notation $\Pi=\{t_{0}=0,t_{1},t_{2},\ldots,t_{n_{1}}=T\}$,
where the $t_{i}$ are simply the enumerated $s_{jp}$.
For our numerical schemes we always assume that Assumption (A2) is in place,
next to (A1).
In particular, we can also simulate i.i.d. copies
of the $(Z_{s_{jp}},\tilde{Z}_{s_{jp}})_{s_{jp}\in\Pi}$.

We formally initialize
$\left(\overline{M}_{j}^{l+1}- \overline{M}_{T}^{l+1}\right)_{j\geq T}=(0)_{T}$.
Suppose that for some (fixed) $j<T$, an approximation $Y_{j+1}^{l+1}$ to the upper
Snell envelope $Y_{j+1}^{\ast,l+1}$ and the set of $\rho$-martingale increments
$\left(\overline{M}_{q}^{l+1}-\overline{M}_{j+1}^{l+1}\right)_{j+1\leq q\leq T}$ have been constructed.
Then, we carry out the following loop.
For $p=n_{0}$, we initialize $0\leq U_{n_{0}}=Y_{j+1}^{l+1}\approx\rho_{j+1}\left(Y_{j+1}^{\ast,l+1}\right)$.
Now, if $0\leq U_{p+1}$, $p<n_{0}$, has been constructed,
we solve the piecewise linear minimization problem
\begin{align}
\label{defag}

\left[  \gamma^{N_1}_{s_{jp}},\beta^{N_1}_{s_{jp}},\tilde{\beta}^{N_1}_{s_{jp}}\right]
=&\underset
{\gamma,\beta,\tilde{\beta}\in\mathbb{R}^{K_1}}{\arg\min}\,\,\frac{1}{N_1}\sum_{n=1}^{N_1}  \bigg(
U^n_{p+1}- \sum_{k=1}^{K_1}\gamma_{k} \psi_{k}(s_{jp},X^n_{s_{jp}})\nonumber\\
\hspace{0.2cm}&  + g\Big(s_{jp},\sum_{k=1}^{K_1}
\beta_{k}\varphi_{k}(s_{jp},X^n_{s_{jp}}),\sum_{k=1}^{K_1}
\tilde{\beta}_{k}\tilde{\varphi}_{k}(s_{jp},X^n_{s_{jp}})\Big)(s_{j(p+1)}-s_{jp}) \nonumber\\
\hspace{0.2cm}& - \sum_{k=1}^{K_1}
\beta_{k}\varphi_{k}(s_{jp},X^n_{s_{jp}})\Delta W^n_{s_{jp}}- \sum_{k=1}^{K_1}
\tilde{\beta}_{k}\tilde{\varphi}_{k}(s_{jp},X^n_{s_{jp}})\Delta\tilde{N}^n_{s_{jp}}  \bigg)^{2},
\end{align}


TITLE: Advice on second order non-linear ordinary differential equation
QUESTION [3 upvotes]: I'm currently working on some problems concerning the calculus of variations and I have come up with the following differential equation that I now want to solve: $$1 + y'(x)^2 - y''(x)(y(x)-\lambda) = 0.$$ I'm only used to solving linear differential equations so I would appreciate some advice on where to start.

REPLY [1 votes]: As abel noted, with $u = y - \lambda$ we get $1 + (u')^2 - u'' u = 0$.
Now if you recall that $\cosh^2 - \sinh^2 = 1$ and $\cosh'' = \cosh$, you
are led to consider something of the form $u = a \cosh(b x + c)$, and we find that this works if $b = 1/a$.  Thus we get the solution
$$ y = \lambda + a \cosh(x/a + c)$$
where $a \ne 0$ and $c$ are arbitrary parameters.


TITLE: abelian groups of order $10^5$
QUESTION [0 upvotes]: I was trying to solve the following problem:
Find number of abelian groups of order $10^5$, up to isomorphism ?
Could someone point me in the right direction?
Thanks in advance for your time.

REPLY [1 votes]: If $n=\prod_{i=1}^rp_i^{k_i}$, then the number of distinct abelian groups of order $n$ is given by
$$
\prod_{i=1}^rp(k_i),
$$
where $p(k)$ denotes the number of partitions of $k$. Now $n=10^5=2^5\cdot 5^5$, and $p(5)^2=7^2=49$.

REPLY [1 votes]: Hint : The number of abelian groups of order $$p_1^{a_1}\cdots p_n^{a_n}$$
upto isomorphism is $$p(a_1)\cdots p(a_n)$$ where $p(n)$ denotes the number of partitions of $n$


TITLE: Generalizing the conjugate gradient like this works?
QUESTION [4 upvotes]: Given $A \in \mathbb{R}^{n \times n}$, a SPD matrix, and a vector $b \in \mathbb{R}^n$, it is possible to solve the problem
$$\min_x \| Ax - b\|$$
with the conjugate gradient method. Its algorithm basically is
$r_0 = b - Ax_0\\
p_0 = r_0\\
k = 0\\
\text{while } r_k \text{ is not small enough}\\
\ \ \ \ \ \ \alpha_k = \displaystyle \frac{ \|r_k\|^2 }{\langle p_k, Ap_k \rangle}\\
\ \ \ \ \ \ x_{k+1} = x_k + \alpha_k p_k\\
\ \ \ \ \ \ r_{k+1} = r_k - \alpha_k Ap_k\\
\ \ \ \ \ \ \beta_k = \frac{ \| r_{k+1} \|^2}{ \| r_k \|^2 }\\
\ \ \ \ \ \ p_{k+1} = r_{k+1} + \beta_k p_k\\
\ \ \ \ \ \ k = k+1\\
\text{end while}\\
\text{return } x_{k+1}$
I wonder what happens if I used the same algorithm for a function like $f(x) = Ax + r(x)$, where $r(x)$ usually is much smaller then $f(x)$ (in the norm sense). To be more precise, I use the same algorithm but changing the formulas for $\alpha_k$ and $r_{k+1}$ by
$\alpha_k = \displaystyle \frac{ \|r_k\|^2 }{\langle p_k, f(p_k) \rangle},\\
r_{k+1} = r_k - \alpha_k f(p_k).$
Is the conjugate gradient method flexible at the point I can do this? What condition I need on $r$ to make it work? One possibility is $r$ being a SPD matrix, but this is an obvious one.

REPLY [1 votes]: Unlikely without special assumptions placed on $f$. The CG algorithm requires $x^T A x \not = 0$ to avoid division by zero. If $A$ is symmetric positive definite and if $v$ is an unit eigenvector of $A$ corresponding to the eigenvalue $\lambda$, then $B = A - \lambda vv^T$ is a symmetric positive semi-definite matrix. In particular, $v^T B v = 0$. If $\lambda$ is the smallest eigenvalue of $A$, then many would consider $f(x) = Bx$ to be only a small perturbation of the original map, especially if the largest eigenvalue of $A$ is much larger than $\lambda$.


TITLE: Solving a characteristics system
QUESTION [2 upvotes]: Problem Sheet Question:

A PDE initial value problem has three associated characteristics equations to solve $x,y,u$. The first, $y$, is a simple ODE and can be solved by integration with a constant.
The other two have the characteristics system as follows:
\begin{align*}
\frac{dx}{dt} &= u\\
\frac{du}{dt} &= x
\end{align*}
Where $x,u$ are both functions of the parametrised variables $s$ and $t$.

Solution:
Quoting from the mark scheme:
  "The above form a system which can be reduced to one equation by eliminating one of the functions. For example, if we eliminate $u$, we will have
  \begin{align*}
\frac{d^2x}{dt^2} &= x\\
\end{align*}
  Which can be integrated using the methods relevant to ODE's. Thus, the general solution is given by:
  $$\ x = C_{1}e^{t} + C_{2}e^{-t} , u= C_{1}e^{t} - C_{2}e^{-t}$$
  Where $C_{1}$ and $C_{2}$ are constants of integration."

My Question:
what is the method in full to obtain $x$ and $u$ using the solution above?

REPLY [1 votes]: First think of it like this:
\begin{align*}
\frac{dx}{dt} &= u\\
\frac{du}{dt} &= x
\end{align*}
can be written as
\begin{align*}
x\frac{d}{dt} &= u\\
u\frac{d}{dt} &= x
\end{align*}
Substituting $u$ in both systems seperatly will give:
\begin{align*}
\frac{d^2x}{dt^2} &= u\\
\frac{d^2u}{dt^2} &= x
\end{align*}
On rearranging we have:
\begin{align*}
\frac{d^2x}{dt^2} & - u= 0\\
\frac{d^2u}{dt^2} & - x= 0
\end{align*}
Which should remind you of a general method for solving 2nd order homogeneous linear ODEs with constant coefficients. Hence we need to solve the quadratic characterisc equation:
$$m^2-1 = 0$$
Which has the roots:
$$m_{1,2}=1,-1$$
Since the roots are real and distinct we have for x:
$$x=C_1e^{M_1t}+C_2e^{M_2t}$$
Which gives you the solution your looking for $x$. However, i don't know why $u$ is negative, i could only find one solution.


TITLE: First and second order conditions - show that they are equivalent
QUESTION [2 upvotes]: Function $ f:\mathbb{R}^n\rightarrow \mathbb{R} $ , given by $ f(x)=(1/2)x^T P x + q^T x + r$, where $P\succeq 0$, Prove this first and second order conditions are equivalent.
In S.Boyd's lecture:

First-order: f is convex if $f(y)\geq f(x)+\nabla f(x)^T(y-x)$

Second-order: f is convex if Hessian $\geq 0$

My attempt:
Second order: Calculate Hessian.
$$
\nabla f(x) = \frac12(P+P^T)x+ q = Px+q
$$
$$
\nabla^2 f(x) = P
$$
If $P\succeq 0$, then convex (Okay, so far).
Now, How do I do the same with the first-order?
I know that first-order and second-order are equivalent. How I do?
My attempt:
First-order: $f(y)\geq f(x)+\nabla f(x)^T(y-x)$
Then,
$$
f(y)\geq (1/2)x^T P x + q^T x + r+ (Px+q)^T(y-x)
$$
I have difficulty simplifying first order

REPLY [0 votes]: $$
f(y) - f(x) - \nabla f(x)(y-x) \ge 0
$$
$$
\varepsilon  = \frac12 (y-x)^TP(y-x)\ge 0
$$
Now, you demonstrate the derivative and find P.