\section{Introduction}
\label{intro}

The classical decision problem of optimal
stopping has found amazing range of applications
from  finance, statistics, marketing, phase transitions to engineering.
While over the past decades many efficient methods have been
developed for its numerical resolution, until recently
essentially all computational approaches in high-dimensions
first approximate the maximal value.
The quantity that is equally
important  is the associated stopping rule, and
the state is divided into two regions depending
on whether at a point it is optimal to stop or to continue.
Although
approximate optimal actions  can either be
calculated directly, or are always
available from the calculated value function
in a greedy policy,
a characterization and computation
of these regions through
interpretable criterion would
provide an immediate
link between the policy and the state.

In many important
applications,
continuation and stopping regions are separated
by graphs of functions that can be used to construct interpretable
optimal stopping decisions, and
our goal is to directly  compute these
 \emph{stopping boundaries}  whose performances
are close to the maximal ones.
In lower dimensions, they
 can be   accurately approximated by methods
based on  nonlinear partial differential equations or dynamic programming.
However, in models with many degrees of freedom,
this numerical problem
is still challenging.

Our approach is based on a
method proposed
by E, Jentzen \& Han \cite{HanE,HanJentzenE},
which we call
\emph{deep empirical risk minimization} (deep ERM).
In this framework,  the controller is replaced
by a deep artificial neural  network
and the random system dynamics is simulated
via Monte-Carlo.
Recent advances in optimization tools
and computational capabilities enable us to
optimize this system efficiently,
thus providing an accurate tool
potentially also applicable to other
problems with stopping boundaries, including models with regime change
or other singularities
 in financial economics,  phase transitions in continuum mechanics,
and Stefan type problems in solidification.

Deep ERM applies essentially to all problems that
can be formulated as dynamic stochastic
optimization, maybe after problem specific modifications.
This and similar algorithms
\cite{BHLP,Becker1,Becker2,Becker4,BGTW,BGTWM,GMS,HPBL,RW,War}
have been applied to many
classical problems from quantitative finance,
financial economics, and  to nonlinear partial
differential equations that  have stochastic
control representations.
Also, the application to optimal stopping
is carried out by Becker, Jentzen \& Cheridito \cite{Becker1,Becker2}
to compute an approximation
of the optimal stopping rule and
hence, the stopping region as well as the maximal value.
They
accurately price many American options
of practical importance in very
high dimensions with confidence
intervals, showcasing
the power and the flexibility of deep ERM
in optimal stopping as well.
We refer the readers to our accompanying papers \cite{RS, RST}
for an
introduction of deep ERM, and to
the  excellent surveys \cite{FMW,HPW,RWs}
and the references therein for more information.

The stopping regions are exactly the
sets of all points at which the value and the pay-off
functions agree, and they can be derived from
the calculated value or
computed directly.  However,
 qualitatively these approximations  may not
 always have
accurate  geometric properties,
and lack interpretability.
Thus a direct computation
of these regions
with their known structures
is desirable,
and it is the main goal
of our approach.
Indeed, in many
financial applications they are characterized
as the epi or hypo-graph of certain functions  in a natural
but problem specific
coordinate system,
described
in Assumption \ref{asm:star}.
We first outline
the algorithm and its properties under this assumption,
and  then verify it in Theorem \ref{t.fb}
for all examples that we consider.
In contrast to methods mapping each state point to stopping decisions, our representation of stopping boundaries as graphs provides topological guarantees, even when approximated.

In our method, the graph separating the
continuation and stopping sets is
approximated by a deep artificial neural network
which then provides a stopping rule
used to calculate an empirical reward function.
However,  it is not possible to directly use any gradient
method for  optimization
over these ``stop-or-go'' decisions.
Therefore,  although at the high-level our approach
is deep ERM, it requires one important modification
to overcome this difficulty.  We replace these sharp  rules
based on the hitting times, by stopping probabilities.
In this relaxed formulation,
when the state is not
close to the boundary,
we still stop or continue the process
as in the non-relaxed case.
But when it  is close to the boundary,
in a region called the \emph{fuzzy} (\emph{stopping}) \emph{boundary}, we stop
with a probability proportional to the distance to the
boundary.
This is analogous to mushy regions in solidification
or phase fields models in continuum mechanics
 \cite{BSS,B-mushy,Mort,osher,Soner,V-mushy}.
Similarly  \cite{Becker1,Becker2}
uses stopping probabilities as relaxed control variables
but without any connection to the geometric structures.

We numerically study several
American and Bermudan options
with this methodology.
The numerical results for the American
put options in Black \& Scholes, and in Heston models
verify the effectiveness of the method.
The Bermudan  max-call option, studied in high dimensions,
shows once again the power of  deep ERM.  The numerical
results for the look-back options
is an example of the flexibility of
this approach in handling path-dependent options.
When possible all our computations
is benchmarked to previous studies,
and shown to be comparable.  Also, confidence intervals and
upper bounds
 computed in \cite{Becker1,Becker3,Becker2}, provide
computable guarantees.
Additionally, the stopping regions,
Figures \ref{t.2}, \ref{t.3}, \ref{fig:asymCall},
for the two-dimensional
max-call options
are qualitatively similar
to those obtained in \cite{BD}.

Direct computation of the stopping boundary
has been the object of several other studies as well.
In low dimensions, techniques from
nonlinear partial differential equations
modeling obstacle problems and phase transitions
allow for accurate calculations.
Also, Garcia \cite{Garcia} proposes a method
similar to ours based on a parsimonious parametrization
of the stopping region.
A recent study by
Ciocan \& Mi\u{s}i\'{c} \cite{CM} proposes
a tree-based method to compute this partition.
Alternatively, deep Galerkin method for the differential
equations is used by
\cite{SS} to compute the value
function of basket options
at every point and the stopping boundary
is derived from this function.
\cite{WP} studies the classical
Stephan problem of melting ice in low dimensions,
by the same approach.  We also refer the readers
to \cite{BTW} for a review
of numerical methods based on boundary parametrization
as well as other approaches to optimal stopping.

For American option pricing with many degrees of freedom,
the main alternative computational tool to deep ERM
are the Monte-Carlo based regression methods described in
Glasserman \cite{Gbook},
Longstaff \& Schwartz \cite{LS} and Tsitsiklis \& van Roy
\cite{TvR}.
The classical book of Detemple \cite{Dbook}
and the recent article of Ludkovski \cite{Lu}
provide extensive information on simulation based computational techniques
for optimal stopping.
The key difference between these approaches, and the
value computation through deep ERM
is the choice of the set of basis functions that grows rapidly
with the underlying dimension.
Deep artificial neural networks
do not require an  \emph{a priori} specification of the basis.
Therefore, they are likely to be more effective
for problems with many states, as supported
by the experiments. Indeed, \cite{Becker3,Kohler}
exploit this property and replace the
preset hypothesis class
by an artificial neural network in their
regression.

The power of these new approaches
are more apparent in high dimensional settings,
and an interesting example
is the American options  with rough volatility models.
These are infinite-dimensional models,
and their numerical analysis
is given in  Bayer \emph{et al.}~\cite{BTW,BQY} and
in Chevalier \emph{at.~al.}~\cite{CPZ}
by alternative methods.
In a different class of problems with many states,
\cite{GMPW} exploits the symmetry of the
underlying problem to build
an appropriate architecture of the neural networks.
This approach is particularly relevant for mean-field
games which models many agents that are identical.
Additionally, \cite{Becker1,Becker3,Becker2} consider
several models with many variables
including the fractional Brownian
motion.

The paper is organized as follows.  The problem and
the value function is defined in the next section.
 The algorithm is introduced in Section \ref{s.algo},
 and the financial structure is outlined in Section \ref{s.f}.
 Section \ref{s.architecture} details the network
 architecture and parameters used in our experiments.
 One dimensional put options in the Black \& Scholes model
 is studied in Section \ref{s.putBS},
 in the  Heston model
in Section \ref{s.putH}, Bermudan
 max-call options are the topic
 of Section \ref{s.maxcall},
and the look-back options are the foci of Section \ref{s.look-back}.
Appendix \ref{a.existence} proves Theorem \ref{t.fb},
Appendix \ref{a.re} proves the convergence
of the reward function $\cR_\epsilon$
as the width of the mushy region tends to zero.
Appendix \ref{a.algo} provides the details of the algorithm.

\section{Optimal Stopping}
\label{s.stopping}

Let $T>0$ be the finite time horizon, and
$(\Omega, \Q)$ be the probability space
with the filtration $ \F=(\cF_t)_{t \in [0,T]}$.
The state process $X = (X_t)_{t\in [0,T]}$
is an $\F$-adapted, continuous, Markov process
taking values in a Euclidean space $\cX$.
For  a fixed subset $\cT \subset [0,T]$,
$\vartheta = \vartheta(\cT)$ is the set of
$\F$-stopping times taking values in $\calT$.
Then, the \emph{optimal stopping} problem,
corresponding to a given \emph{reward function}
$\varphi$,
is
\begin{equation}
\label{eq:OS}
\text{to maximize}
\ \
v(\tau):=  \E[ \varphi(\tau,X_{\tau})],
\qquad
\text{over all}\  \tau \, \in \, \vartheta,
\end{equation}
where for some $a>0$, $\varphi(t,\cdot) \in \cL_a$ defined as
\begin{equation}
\label{e.lca}
\cL_{a}:= \left\{ \phi :\cX \to \R\ :\ \text{continuous and}\
|\phi(x) | \le C[1+|x|^a]\ \
\text{for some}\ C>0 \right\}.
\end{equation}

In financial examples, $a$ is almost
always equal to one.  Also, $\cT=[0,T]$
corresponds to an American option,
while $\cT$ is a finite set for Bermudan
ones.  However, since for numerical calculations
one has to discretize the time variable, from the onset
we assume that
$\cT=\{t_0, t_1, \ldots,t_{n-1}, T\}$,
with $n$ being large for the American ones.

We next define  the value function
which is a central tool in Markov optimal control  \cite{FS}.
For $t \in \cT$ and $x \in \cX$,
let $\vartheta_t$ be the set
of all $[t,T]\cap \cT$-valued stopping times,
and set
$$
v(t,x):= \sup_{\tau \in \vartheta_t}\ v(t,x,\tau)
\quad \text{where}
\quad
v(t,x,\tau):=
\E\left[  \varphi(\tau,X_\tau) \mid X_t=x\right].
$$

Finally, we introduce our notation  $\R_+=(0,\infty)$,
and
$$
\cT^\circ:= \cT \setminus\{T\}
=\{t_0, t_1, \ldots,t_{n-1}\}.
$$

\subsection{Stopping region}
\label{ss.cS}

The following subset of the state space is the
\emph{stopping region} at time $t\in \cT$:
$$
\cS_t := \left\{ x \in \cX \ :\ \varphi(t,x) = v(t,x) \right\}.
$$
As  $v(T,x)=\varphi(T,x)$ for every $x\in \cX$, we have $\cS_T=\cX$.

To obtain some regularity of the value function
and the stopping regions,
we make the following mild assumption
satisfied in all applications. Recall
that  $\cL_a$ is defined in \eqref{e.lca}.
\begin{asm}{\rm{\bf{(Growth and continuity)}}}
\label{a.growth}
There exists $a>0$,
such that for every $t \in \cT$,
$\varphi(t,\cdot) \in \cL_a$.
Moreover, for every $\phi \in \cL_a$, $ t < \bar{t} \in \cT$,
$v_\phi (t,\cdot,\bar{t}) \in \cL_a$, where
$$
v_\phi (t,x,\bar{t}):=
\E[\phi(X_{\bar{t}}) \mid X_{t}=x],
\qquad x \in \cX.
$$
\end{asm}
\noindent
The above assumed continuity essentially follows form
the continuous dependence of the distribution
of the state
on its initial value, and it is always satisfied
in all our applications.

The following are direct consequences.

\begin{proposition}
\label{p.vlsc} Under the  growth and continuity
Assumption \ref{a.growth},
for each $ t\in \cT$,
the value function $v(t,\cdot) \in \cL_a$,
and the stopping region $\cS_t$ is a relatively closed subset of $\cX$.
\end{proposition}
\begin{proof}
To prove the first statement,
we use backward induction in time
$t_0<t_1<\ldots<t_n=T$.
First observe that by Assumption \ref{a.growth},
$v(T,\cdot)=\varphi(T,\cdot) \in \cL_a$.
For the induction hypothesis,
suppose that $v(t_{k+1},\cdot) \in \cL_a$
for some $t_{k+1}\in \cT$.
By dynamic
programming,
$$
v(t_k,x)= \max\left\{ \varphi(t_k,x),\ c(t_k,x) \right\},
\qquad\text{where}
\qquad
c(t_k,x):= \E[v(t_{k+1},X_{t_{k+1}}) \mid X_{t}=x].
$$
As $v(t_{k+1},\cdot)$ is assumed to be in $\cL_{a}$,
 Assumption \ref{a.growth}
implies that
$c(t_k,x)$
is also in $\cL_a$.
Since the space $\cL_a$ is closed
under maximization, the above equation
implies that $v(t_k,\cdot) \in  \cL_a$.
Hence, $v(t,\cdot) \in \cL_a$ for every $t\in \cT$, and in particular,
it is continuous.

The stopping region $\cS_t$ is the zero level set
of  $v(t,\cdot)-\varphi(t,\cdot)$ which is shown to be continuous.
Hence, it is  relatively closed in $\cX$.
\end{proof}

Because the stopping regions are closed (hence Borel) and
$\cS_T=\cX$, the hitting time
\begin{equation}
\label{e.taustar}
\tau^* := \min \left\{ u \in \cT\     :\
X_u \in \cS_u \right\}
\end{equation}
is a well-defined stopping time and belongs to $\vartheta$.
Moreover,   $\tau^*$
is optimal, or equivalently,  with $v(\tau)$  as in \eqref{eq:OS},
$$
v(\tau^*)=\sup_{\tau \in \vartheta} v(\tau).
$$

\subsection{Stopping boundary}
\label{ss.free}

In all problems that we consider, the
boundary of stopping regions are
given as the graphs of certain functions,
but possibly in a different coordinate system.
To outline the algorithm in a general manner,
we first assume the existence of
this geometry.
We then verify it for the
financial examples, in Theorem \ref{t.fb} below,
under mild natural structural conditions.

Although the description of the following assumption is technical,
its interpretation is quite natural.
In words, we assume that there exists a stopping boundary
that can be characterized as the graph of some function
$f^*$ in the appropriate coordinates $(\alpha, \Xi)$.
In all applications, the choice of these functions
is immediate as discussed in
Example \ref{ex.maxcall} below,
and in Sections \ref{s.putBS}, \ref{s.putH},  \ref{s.maxcall},
\ref{s.look-back}. In particular,
the positive real-valued $\alpha$
is related to the payoff of the option, and
$f^*$ in \eqref{e.st} describes the stopping boundary.
The function $\Xi$ is more abstract and is a transformation
into coordinates that are meaningful to the specific problem.

\begin{asm}{\rm{\bf{(Existence of a stopping boundary)}}}
\label{asm:star}
There exist measurable functions
$$
\alpha: \cX \to \R_+,
\quad
\Xi: \cX \to \Xi(\cX) \subset \cX,
\quad
 f^* : \cT^\circ \times \ \Xi(\cX) \to [0,\infty],
 $$
 and $\eta \in \{-1,+1\}$,
 so that the map $x \in \cX \mapsto (\alpha(x), \Xi(x))
 \in \R_+ \times \ \Xi(\cX)$ is a homeomorphism
 {\rm{(}}i.e., one-to-one, onto,  continuous with a continuous inverse{\rm{)}}
 and for every $t \in \cT^\circ$,
\begin{equation}
\label{e.st}
\cS_t = \{  x \in  \cX \ :\
\eta \left(f^*(t,\Xi(x))-\alpha(x)\right) \le 0  \}.
\end{equation}
\end{asm}

The pair $(\alpha(x),\Xi(x))$ is
a representation  of $x$
that is adapted to the problem.
The space
$\Xi(\calX)$ is a \textit{latent space}
with  (Hausdorff) dimension strictly
less than that of $\calX$.
The parameter $\eta$ determines
whether the stopping region is an
epigraph or a hypograph of $f^*$.
In the financial applications,
typically $\eta=1$ corresponds to a call
type option and the stopping region is an epigraph.
While $\eta=-1$ is related to put options.

The surjectivity of $\alpha$ ensures
that this characterization is non-trivial,
and it additionally implies that for every $c>0$, the point corresponding
to the pair $(f^*(t,\Xi(x))+\eta c, \Xi(x))$ is in
the stopping region $\cS_t$, and respectively, to
the pair $(f^*(t,\Xi(x))-\eta c, \Xi(x))$ is not in $\cS_t$.  Therefore,
the graph of $f^*$ in the coordinates
$(\alpha(x),\Xi(x))$ given by $$\{x\in \cX : \alpha(x)=f^*(t,\Xi(x))\},$$
is included in the boundary of $\cS_t$.
In fact, in all examples
it separates the stopping and continuation regions.
Thus, we call $f^*$ the optimal \emph{stopping boundary}.

\begin{example}
\label{ex.maxcall}
In the financial problem of an American
max-call option with two stocks
and a given strike $K$, we have
$\varphi(t,x)= e^{-rt}(\max\{ x_1,x_2\} - K)^+$
 for $x=(x_1,x_2) \in \cX$.
With Markovian dynamics, the state space is $\cX =\R_+^2$.
In Theorem \ref{t.fb} below, we show that
 there exists a stopping boundary $f^*$
such that the triplet
$(\alpha, \Xi, f^*)$
with
$$
\alpha(x)= \max\{ x_1,x_2\},
\qquad \Xi(x)=\frac{x}{\alpha(x)},
\qquad \eta=1,
$$
satisfies the above assumption.
Figures \ref{t.2}, \ref{t.3}, \ref{fig:asymCall} below provide
examples of these regions.
In our numerical approach,
the latent variable
$\Xi(x)$ is the input
to the network approximating the stopping boundary.
As $\Xi(x) = \Xi(\gamma x)$
for all $x\in \R_+^2,$  $\gamma > 0$,
every ray emanating from the origin is
mapped to a single point by the map $\Xi$,
and the threshold $f^*(t,\Xi(x))$
separates the stopping and continuation
regions on this ray.

In the simpler example of an American put
option with one Markovian stock, one simply takes
$\cX=\R_+$, $\alpha(x)=x$, $\Xi\equiv 0$,
and $\eta=-1$.  So that  the
stopping region is the hypograph of some
stopping boundary $f^*$.  Several other examples are discussed
in the sections below.

\qed
\end{example}

As an immediate consequence of the above assumption,
we can restrict the maximization in \eqref{eq:OS} to
stopping times given by stopping boundaries.  For future
reference, we record this fact.  Let $\frakF$ be the set of
all measurable functions $f :\cT^\circ \times \ \Xi(\cX) \to [0,\infty]$.
For $f\in \frakF$, the corresponding
stopping time is given by,
$$
\tau_f:=  \min \left\{ u \in \cT\ :\
\eta\left( f(u,\Xi(X_u))-\alpha(X_u)\right) \le 0 \ \ \text{or}\ \ u=T\right\}
\in \vartheta.
$$

\begin{lemma}
\label{l.restrict} Under the growth and continuity Assumption \ref{a.growth},
and the existence of a stopping boundary Assumption
\ref{asm:star},
$$
\sup_{\tau \in \vartheta} v(\tau) = \sup_{f \in \frakF} v(\tau_f).
$$
\end{lemma}
\begin{proof}
As $\tau_f \in \vartheta$,
we have $\sup_{\tau \in \vartheta} v(\tau) \ge
v(\tau_f)$  for all  $f \in \frakF$.
Let $f^*$ be as in Assumption \ref{asm:star}.
Then, for $u \in \cT$, we have $X_u \in \cS_u$ if and only if,
either $\eta\left(f^*(u,\Xi(X_u))-\alpha(X_u) \right) \le 0$, or $u=T$.
Therefore, $\tau_{f^*}$ is equal to the
optimal hitting time $\tau^*$ defined in \eqref{e.taustar}.
These imply that
$\sup_{\tau \in \vartheta} v(\tau) \ge \sup_{f \in \frakF} v(\tau_f)
\ge v(\tau_{f^*}) = v(\tau^*)
= \sup_{\tau \in \vartheta} v(\tau)$.
\end{proof}

\section{The Algorithm}
\label{s.algo}
Our goal is to construct an efficient
algorithm for the calculation
of the interpretable
stopping decisions based on boundaries
approximated by artificial neural networks.
We train the networks by the \emph{deep
empirical risk minimization} algorithm proposed
by Weinan E, Jiequn Han, and Arnulf Jentzen
\cite{HanE,HanJentzenE}.

This method
approximates the expected value of a reward
function by Monte-Carlo simulations and
optimizes it by stochastic gradient ascent.
For the problem of optimal stopping,
the natural choice for the reward function  is
$\varphi(\tau_{\theta},X_{\tau_{\theta}})$. However, the map $\theta \mapsto \tau_\theta$
is piece-wise constant, making
optimization by gradient ascent impossible.
Therefore, we replace the hitting times
by stopping probabilities based on
fuzzy boundaries that are also used in  some
problems of solidification  \cite{B-mushy,Mort,Soner,V-mushy}.

\subsection{Fuzzy Stopping Boundary}
\label{ss.fuzzy}
Let $\eta\in \{-1,1\}$ be as in Assumption \ref{asm:star}.
For  $f \in \frakF$, set
$$
d(t,x;f):= \eta \left(f(t,\Xi(x))-\alpha(x)\right), \qquad
t \in \cT^\circ, \ x \in \cX.
$$
In the original problem, once a boundary $f$ is chosen,
we stop
the process only when $d(t,X_t;f)\le0$.
As training is not possible
with this sharp ``stop-or-continue'' rule,
we introduce the fuzzy boundaries.  Namely, we fix a
\emph{tuning parameter} $\epsilon>0$, and
at time $t \in \cT^\circ$,
define
$$
\frakB_t^\epsilon(f) :=\{ x\in \cX \ :\ |d(t,x;f)|<\epsilon\}
$$
be the
\emph{fuzzy region}.
\begin{figure}[H]
\begin{center}
\includegraphics[height=1.8in,width=3.2in]{Figures/fuzzyBdry1DCall.pdf}
\caption{
Fuzzy region of an American call (in purple)
connecting the stopping region (in red)
to the continuation region (in blue). }

\label{fig:fuzzy}
\end{center}
\end{figure}
\vspace{-0.5em}
\noindent
An illustration of the fuzzy region is given in \cref{fig:fuzzy} for an
American call option.  If the state is not in
$\frakB_t^\epsilon(f)$,
we continue or stop, as before, with probability one.
But if we are in the fuzzy region,
we stop with  probability $p_t(X_t,f)$, where
$$
p_t(x,f):= \left(\frac{\epsilon - d(t,x,f)}{2\epsilon}\right)^+\wedge 1,
\qquad   t \in \cT^\circ,\  f \in \frakF, \ x \in \cX.
$$
At maturity, we stop at all states, or equivalently,
we set $p_T(\cdot,f)\equiv 1$.
Precisely, $p_t(X_t,f)$
is the probability of stopping  at time $t \in \cT$,
conditioned on the event that the process is not stopped
prior to $t$.  The relaxed control problem is then defined
using the reward function given by
\begin{equation}
\label{e.ell}
\cR_\epsilon(X,f):=
\sum_{t\in \cT}\ p_t(X_t,f)\ b_t(X,f)\ \varphi(t,X_t),
\end{equation}
where  $b_t(X,f)$ is the probability
of not stopping strictly before time $t$
and it  is
obtained as the solution of the following difference equations,
$$
b_{t+1}(X,f)
=b_{t}(X,f) (1-p_{t}(X_t,f)), \qquad t \in \cT^\circ,
$$
with $b_0(X,f)=1$.
In other words,
$b_t(X,f)$ is the unused
\emph{stopping-budget} remaining.
The process $b_t(X,f)$ is non-increasing in
time and takes values in $[0,1]$.
Connection
between $\E[\cR_\epsilon(X,f)]$ and
$v(\tau_f)$ is
discussed in Appendix \ref{a.re}.

We refer to \cite{RST} for a different exposition of this relaxation.

\subsection{Deep ERM}
\label{ss.derm}

In this method, we
approximate the stopping boundary
by a set of deep neural networks
that we abstractly parameterize
by $\theta$ taking values in a finite-dimensional space $\Theta$
as follows:
$$
\frakN=\left\{\ g(\cdot; \theta) :\cT^\circ \times \ \Xi(\cX) \to \R_+\ :
\ \theta \in \Theta\ \right\} \subset \frakF.
$$
Let $\tau_{\theta}=\tau_{g(\cdot;\theta)}$
be the hitting time corresponding to
the stopping boundary $g(\cdot;\theta)$.
The following is the  pseudocode of the deep
empirical risk minimization based on
the  fuzzy stopping boundary.
\vspace{10pt}

\noindent
\begin{enumerate}
\item \emph{Initialize}  $\theta \in \Theta$.
\item \emph{Simulate}
independent state trajectories,  $\{X^{(1)}, \ldots, X^{(B)}\}$,
where $B$ is the \emph{batch size}.
\item Compute the \emph{empirical reward}  function:
$$
R_\epsilon(\theta):=\frac1B\ \sum_{i=1}^{B}\ \cR_\epsilon(X^{(i)},g(\cdot;\theta)),
$$
where $\cR_\epsilon(X,f)$ is as in \eqref{e.ell}. A
modification $\cR_{\epsilon,\lambda}$, given in \eqref{e.lambda},
allows for \emph{important sampling} discussed in Section
\ref{ss.important} below.
\item \emph{Optimize and update by stochastic gradient ascent}
with the  learning rate process $\zeta$:
$$
\theta \, \leftarrow \,
 \theta + \zeta  \nabla_\theta R_\epsilon(\theta).
 $$

\item \emph{Stop} after  $M$ \emph{number of iterations}.
\item \emph{Compute the initial price}
given by the training network, using $J$ many
Monte-Carlo simulations with a sharp boundary.
\end{enumerate}
\vspace{10pt}

A formal discussion of the algorithm is given in the next subsection.
In  Appendix \ref{a.algo}, we also provide more details
of the above algorithm for the benefit
of readers interested in the implementation.
\cite{Garcia} proposes a similar
approach by assuming that the stopping region
has a finite dimensional representation
without a formal justification.
Also,
Chapter 8.2 of \cite{Gbook} discusses a general
approximation for the value calculation by parametrizing
the stopping times. In our approach
we approximate the stopping boundary directly.

\subsection{Discussion}
\label{ss.discuss}

The general theory of stochastic gradient ascent  \cite{SGD} implies that
the above algorithm constructs an approximation of
$$
\theta_\epsilon^* \in \arg\max_{\theta \in \Theta}\
\E\left[\cR_\epsilon(X,g(\cdot;\theta))\right].
$$
In the context of optimal
stopping, \cite{bel}
also provides rates of convergence.
Moreover, in Lemma \ref{l.conv} below, we
show that for each $\theta \in \Theta$,
$$
\lim_{\epsilon \downarrow 0}\  \E[\cR_\epsilon(X,g(\cdot;\theta))]=
\E[\varphi(\tau_\theta,X_{\tau_\theta})]= v(\tau_\theta).
$$
Therefore,
$$
\E\left[\cR_\epsilon(X,g(\cdot;\theta_\epsilon^*))\right]
= \sup_{\theta \in \Theta}\
\E\left[\cR_\epsilon(X,g(\cdot;\theta))\right]
\approx
 \sup_{\theta \in \Theta}\ v(\tau_\theta).
 $$
Also, by the celebrated universal approximation theorem \cite{Cybenko,Hornik},
continuous stopping boundaries are well approximated
by the hypothesis class $\frakN$, and in all our applications
the optimal stopping boundaries are continuous.
Thus, in view of this
approximation capability of $\frakN$ and Lemma~\ref{l.restrict},
we formally have
$$
\E\left[\cR_\epsilon(X,g(\cdot;\theta_\epsilon^*))\right]
= \sup_{\theta \in \Theta}\
\E\left[\cR_\epsilon(X,g(\cdot;\theta))\right]
\approx \sup_{\theta \in \Theta}\ v(\tau_\theta) \approx
\sup_{f \in \frakF}\ v(\tau_f)
=\sup_{\tau \in \vartheta}\ v(\tau).
$$

Hence, the above algorithm
constructs an artificial neural network
with an  asymptotic performance formally close to the
optimal value, justifying
the algorithm.
However, there are several
technical issues in front a rigorous statement
and we postpone this convergence analysis
to a future study.

\section{Financial Examples}
\label{s.f}
We numerically study several
American or Bermudan options to illustrate
and also assess the proposed algorithm.
In these applications, $\Q$ is a risk-neutral measure and the Markov state process
$X=(S,Y,Z)$ has three components: $S$ is
the prices of the underlying stocks
possibly including the past values to make it Markov,
a factor process $Y$ used in models like Heston,
and a functional $Z$ of the
stock process, making
the pay-off of a path-dependent option a function of
the current value of $X$.  Typical examples
of $Z$ are the running maximum or minimum in look-back options,
and  the average stock price for the Asian options.
Depending on the model and the option, $Y$, $Z$, or neither may
be present.

In our examples,
the pay-off function
always has the following form:
\begin{equation}
\label{eq:payoff}
\varphi(t,x)
= e^{-rt} \left( \eta  (\alpha(s,z) -\beta(s,z)-K)\right)^+,
\qquad t \in \cT, \ x=(s,y,z) \in \cX,
\end{equation}
where $\alpha>0, \beta \ge 0$
are  \emph{positively homogenous of degree one},
and the parameter $\eta\in\{-1,+1\}$
determines the option type:
$\eta=1$ is a call and $\eta=-1$  is a put.
The choice of $\alpha$ and $\beta$ may not be
unique, and we would choose $\beta\equiv0$
when possible. It is clear that $\varphi$
grows linearly.  Additionally, in all our examples,
the stock price distribution depends
smoothly on the initial data.  Hence, the growth
and continuity Assumption \ref{a.growth} holds with $a=1$.

The following stopping boundary result is proven in the Appendix.
For max-call options a similar result is proved
in Proposition A.5 in \cite{BD}, and for look-back options in \cite{DaiKwok}.
\begin{theorem}
\label{t.fb} Suppose that in \eqref{eq:payoff},
 $\alpha : \calX \to \R_+$,  $\beta: \calX \to [0,\infty)$
are positively homogenous of degree one, and that
 $X$ scales
linearly in the initial values $S_t=s, Z_t=z$, i.e., for
$x=(s,y,z)$, $\gamma>0$,
and a bounded, continuous, integrable function $\phi$,
the following holds:
$$
\E\left[ \phi(u,X_u) \mid  X_t =(\gamma s, y, \gamma z)\right]
= \E\left[ \phi(u,\gamma X_u) \mid  X_t =(s, y,z)\right].
$$
Further assume that the growth and continuity
Assumption \ref{a.growth} holds, and
the price of the European option is strictly positive
for all $t<T$.  Then, the existence of a stopping boundary Assumption
\ref{asm:star} is satisfied with  $\alpha $
as in \eqref{eq:payoff},
\begin{equation}
\label{e.xi}
\Xi(s,y,z)= \left(\frac{s}{\alpha(s,z)}\,
,\, y\, ,\, \frac{z}{\alpha(s,z)}\right),
\end{equation}
and
\begin{equation}
\label{e.fstar}
f^*(t,\xi):=
\left\{
\begin{array}{ll}
\sup\left\{ \alpha(s,z)\ :\ x\in \cS_t\ \ \text{s.t.}\ \ \Xi(x)=\xi \right\}, \quad
 & \text{if}\ \ \eta=- 1,\\
\, \inf\left\{ \alpha(s,z)\ :\ x\in \cS_t\ \ \text{s.t.}\ \ \Xi(x)=\xi \right\}, \quad
 & \text{if}\ \ \eta=1,
 \end{array}
 \right.
\end{equation}
where we use the convention that the sup over the empty set is zero and
inf is infinity.
\end{theorem}

The positivity of the European price
ensures that on stopping regions $\varphi>0$.

\subsection{Importance sampling}
\label{ss.important}

To effectively train the network  $g(\cdot;\theta)$
by deep ERM,
it is  crucial  that the simulated paths
$X_t$ hit the fuzzy  boundary $\frakB_t^\epsilon(g(\cdot;\theta))$
 frequently
enough for each exercise date.
We achieve this by importance sampling
that we now outline.

In all our examples,
$S$ evolves according to a stochastic differential equation,
\begin{align*}
 dS_t = S_t\left(\mu_t dt + \sigma_t dW_t\right),
\quad s_0 \in \R^d,
\end{align*}
where the processes $\mu \in \R^d$,
$\sigma \in \R^{d\times d}$
may depend on $t,S_t$ and other exogenous factors.
For $\lambda=(\lambda_1,\dots,\lambda_d)  \in \R^d$, let
$\Q_\lambda$ be the measure obtained by
 the Girsanov transformation so that
 $W_t^{\lambda} :=  W_t +\lambda t $ is a
 $\Q_{\lambda}$ Brownian motion.
 Hence, under $\Q_\lambda$,
$$
dS_t =  S_t\left([\mu_t -   \sigma_t \lambda ] dt
+ \sigma_t dW^{\lambda}_t\right).
$$
We then adjust the hitting
frequency of the simulated trajectories,
by choosing $\lambda$ appropriately,
thus allowing for more efficient training.
Although in all our experiments we have employed constant $\lambda$,
one could also use time dependent
ones as well.

However, as the probability measure is modified,
we need to account for it in the reward function.
This is achieved by the Radon-Nikodym derivative
given by,
$$
\cZ^{\lambda}_t (S):=
\frac{d\Q}{d\Q_{\lambda}} \big |_{\calF_t}
= \exp( \lambda \cdot W_t^{\lambda}  -\frac12
|  \lambda|^2 t).
$$
Since $W^\lambda$ can be expressed in terms of $S$,
the process $\cZ^\lambda$
can be viewed as a function of the state process $S$.
We then modify the reward function given in \eqref{e.ell}
by,
\begin{equation}
\label{e.lambda}
\cR_{\epsilon,\lambda}(X,f):=
\sum_{t\in \cT}\ \cZ^{\lambda}_t(S)\ p_t(X_t,f)\ b_t(X,f)\ \varphi(t,X_t).
\end{equation}

\section{Network Architecture and Parameters}
\label{s.architecture}

All our numerical  experiments have been carried out
with Tensorflow 2.7  on a 2021 Macbook pro with 64GB unified memory
and Apple M1 Max chip. The code is implemented in
Python and run on CPU (10-core) only.
Throughout the examples, we use the same neural
network architecture and training parameters,  given  below:
\vspace{5pt}

\noindent
{\bf{1.}} We use only one deep neural network
that takes both the time and
the state vector as input.  Then,
the trained boundary is a function of the
continuous time variable, delivering smoother
dependence on time.  This is in contrast to
majority of previous studies that
use one deep neural network for each time
in the discrete set $\cT$.
\vspace{5pt}

\noindent
{\bf{2.}}
The single feedforward neural network that
we employ
consists of $2$ hidden layers
 with leaky rectified linear unit (Leaky ReLU) activation function.
 Each layer has $20+\bar{d}$ many nodes where $\bar{d}$ is the dimension of the latent space $\Xi(\calX)$.
The output layer uses the standard rectified linear unit (ReLU)  activation.
An important parameter
is the  \emph{bias  $\theta_0 \in \R_+$ in the output layer}
that sets  an initial boundary level.
Namely, the output function $g$ is given by,
$$
g(t,\xi; \theta_1,\theta_0,\theta_{in}) = \left(\theta_1 \cdot g_{in}(t,\xi; \theta_{in})+\theta_0\right)^+,
$$
where  $g_{in}(t,\xi; \theta_{in}) \in \R^{20+\bar{d}}$
is the value of the neural network  before entering the output layer,
with parameters $\theta_{in}$.
In our experiments, we have observed that the algorithm
is not sensitive to the initial bias $\theta_0$,
as long as it is set to a reasonable value.  For example,
for options with a given strike $K$,
we  choose $\theta_0=\frac{3K}{2}$,  $\theta_0=\frac{K}{2}$ for  calls
and puts,  respectively.
\vspace{5pt}

\noindent
{\bf{3.}}
We set the number of stochastic gradient
iterations to $M=3,000$ with a
fixed simulation batch size of $B = 512$.
\vspace{5pt}

\noindent
{\bf{4.}} The learning rate
process $\zeta$
is taken from the  Adam optimizer \cite{Kingma}.

\vspace{5pt}

\noindent
{\bf{5.}}
The fuzzy region width $\epsilon$ is chosen to be the strike $K$
times the standard deviation of the  increments of $\alpha(X)$.
While the strike $K$ adjusts to the scale of the payoff,
the second term reflects the typical variation
of the  process $\alpha(X)$ between exercise dates.
\vspace{5pt}

\noindent
{\bf{5.}}
After the training is completed,
the initial price is computed using the sharp boundary and
$J= 2^{22} = 4,194,304$ many Monte Carlo simulations.

\section{One-dimensional Put Option}
\label{s.putBS}
As an initial study, we consider an at the money
 Bermudan put option  in the Black Scholes model.  Hence,
 $\varphi(t,s) = e^{-rt}(K-s)^+$, $\eta=-1$,  and $X_t=S_t$ with $\cX=\R_+$.

 We take the same model and parameters
 as in  Section 4.3.1.2 of \cite{Becker2}. Namely,
 $S_t$ follows the Black-Scholes model with parameters
\begin{equation}
\label{e.p1}
s_0=K=40, \ r = 6\%,\ \sigma = 40\%,\ T=1, \ n=50,
\end{equation}
so that $\cT = \{0,\frac{1}{50},\ldots,1\}$.
For efficient training,
we want the  the stock price to reach low values
to cross the boundary.  Towards this goal, we use
use importance sampling with $\lambda =0.275$,
so that
$S$ becomes a super-martingale under $\Q_{\lambda}$
with $5\%$ negative drift. The fuzzy region width is approximately equal to
$\epsilon = K\sigma \sqrt{T/n}$.

It is classical that the existence of a stopping boundary Assumption
\ref{asm:star} holds with $\alpha(s)=s$ and $\Xi(s)\equiv0$.
Table \ref{tab:resultPutBS} summarizes the
performance values
and compares it to the  values computed by Becker \emph{et.al.}~\cite{Becker1}.
The first column in that table
is the average price of ten experiments,
and the second column the highest price.
In practice, one would of course retain
the trained boundary yielding the highest price.

\begin{table}[H]
\begin{center}
{\bf{Bermudan with parameters as in \eqref{e.p1}}}
\vspace{5pt}

\begin{tabular}{cccc}
\hline \hline
Average Price  & Highest Price & Runtime  & Price in \cite{Becker2}  \\
\hline  \hline
5.308 (0.003)  &  5.313 (0.003) &  57.9 & 5.311 \\
\hline \\[-1em]
\end{tabular}
\caption{Bermudan Put Option with $n=50$ in Black-Scholes model.
First  column is the  average  of ten experiments with its standard deviation in brackets.
The second column is the  highest price  among the realizations and its standard deviation in brackets.
Third column is the  average runtime (in seconds) per experiment for the training phase.}
\label{tab:resultPutBS}
\end{center}
\end{table}

Figure \ref{fig:putBS} displays the initial boundary
and the trained one.  The optimal
boundary, shown in blue, is computed
using a finite-difference scheme.
Figure \ref{fig:putBS2} is the approximation
of the  optimal stopping boundary for the American option
with same parameters.  As the American option
does not restrict the trading dates, we
use $n=250$.  Also,  to capture the known
singular behavior of the optimal boundary,
we use  a
non-uniform discretization that is finer near
maturity.

\begin{figure}[H]
\begin{center}
\includegraphics[height=2in,width=5in]{Figures/Bdry50.pdf}
\caption{Black-Scholes model with parameters
as in \eqref{e.p1}.
Blue trajectories are simulated paths.}%
\label{fig:putBS}
\end{center}
\end{figure}

\begin{figure}[H]
\begin{center}
\includegraphics[height=1.8in,width=5in]{Figures/1DGood.pdf}
\caption{Almost American option in
Black-Scholes model with parameters
 as in \eqref{e.p1} except $n=250$.}%
\label{fig:putBS2}
\end{center}
\end{figure}

\section{Put Option in the Heston Model}
\label{s.putH}

Consider the Bermudan put option of the  Section \ref{s.putBS}
(i.e., $s_0 =K=40$, $T=1$, $n=50$) but
now in the Heston model.
Then, under a risk-neutral measure $\Q$,
the stock and factor dynamics are given by,
\begin{align*}
\frac{dS_t}{S_t} &= r dt + \sqrt{Y_t}\  dW_t, \\
dY_t &= (\kappa(\bar{y}-Y_t) - \gamma_{\Q} Y_t)dt +
\sigma_{Y} \sqrt{Y_t}\  d\tilde{W}_t,

\end{align*}
where  $W, \tilde{W}$ are Brownian motions
with constant correlations  $\rho$,
and $ \gamma_{\Q}$ is the volatility risk premium
related to the risk-neutral measure $\Q$.
We use the parameters
 $$
 r = 6\%,\ \kappa = 1,\
 \sqrt{Y_0} = \sqrt{\bar{y}} = 40\%,  \
 \gamma_{\Q} =0,\  \sigma_{Y} =0.1,\ \rho =-0.5.
 $$
 In particular,
the Feller condition $2 \kappa  \bar{y} \ge  \sigma_{Y}^2$ is satisfied.
To allow for a comparison
with results of Section \ref{s.putBS},
we set the initial value and the long term mean $\sqrt{\bar{y}}$ of
the stochastic volatility $\sqrt{Y}$
equal to the Black-Scholes volatility $\sigma$ from
that section.   Moreover, we use importance sampling with the same
Girsanov parameter, $\lambda =  0.275$.
The  variance process is simulated using the Milstein scheme \cite{KP}
 with $n' =4n=200$ time steps.

 In this example, $X=(S,Y)$, $\cX=\R_+^2$,
 $\alpha(x)=s$, $\beta \equiv 0$, and $\Xi(x)=(1,y)$.  We simplify
and set  $\Xi$  to $y$.
Figure \ref{fig:heston1} shows that the stopping boundary is a function
of time $t$ and the spot variance $Y_t =y$.
As $\varphi$ is convex and bounded,
 the value function $v(t,s,y)$
is  non-decreasing in $y$ for all $t\in [0,T]$, cf.~\cite{LambertonHeston}.
Consequently, the map $y \mapsto f(t,y)$ is non-increasing
and the trained neural network $g(\cdot,\cdot;\theta^M)$
has the same property as seen
 in  \cref{fig:heston1}.
Notice also that $y \mapsto f(t,y)$ becomes less steep as $t$
approaches maturity. This is line with the fact that the
vega of the option  decreases over time until the option is exercised.
Since $s \mapsto v(t,s,y)$ is non-increasing for fixed $(t,y)$,
we therefore expect that the  rectangle $[0,s]\times [0,y]$
is contained in $\calS_t$ whenever $(s,y) \in \calS_t$.
This is also confirmed in  \cref{fig:heston1}.
\vspace{10pt}

\begin{figure}[h]
\centering
\includegraphics[height=1.8in,width=5in]{Figures/2DPlotHestonY.png}
\caption{Stopping and continuation regions for a Bermudan put in the Heston model.}
\label{fig:heston1}
\end{figure}

 Table \ref{tab:resultPutHeston} compares the prices
 of the claims when the stopping decision exploits the spot variance
 (i.e., $f = f(t,y)$) and when it does not (i.e., $f = f(t)$).
 We observe that the price obtained with the ``blind'' strategy
 is only slightly inferior. Interestingly, the neural network coincides with
 the optimal boundary in the Black-Scholes model.
Also note that the longer runtime compared to Table \ref{tab:resultPutBS}
 is due to the thin partition needed to simulate $Y$.

\begin{table}[H]
\centering
{\bf{Bermudan Put Option ($n=50$), Heston model}}
\vspace{3pt}

\begin{tabular}{c c c c }
\hline  \hline
Threshold &  Average Price& Highest Price & Runtime  \\
\hline  \hline
$f(t,y)$ & 5.290 (0.004) & 5.294 (0.003) &  194.3 \\
$f(t)$ & 5.288 (0.004) & 5.293 (0.003) &  193.9  \\
\hline\\[-1em]
\end{tabular}
\caption{Second  column is the
average  of ten experiments with its standard deviation in brackets.
Third column is the highest price
among the realizations and its standard deviation in brackets.
Fourth column is the  average run-time (in seconds) per experiment for the training phase.}
\label{tab:resultPutHeston}
\end{table}

Finally, \cref{fig:heston2} displays a trajectory of $(X,Y)$ together with the
trained stopping boundary
 $g(\cdot;\theta^M)$. As can be seen, the threshold is typically below the
 Black-Scholes boundary when $Y_t$ is above its mean $\bar{y} = 0.16$ and vice versa.

\begin{figure}[h]
\centering
\includegraphics[height=1.8in,width=4in]{Figures/BdryHeston2.pdf}
 \caption{Stock price, variance and threshold process for a Bermudan put in the Heston model.}
\label{fig:heston2}
\end{figure}

\section{Bermudan Max-Call}
\label{s.maxcall}

We now apply the deep empirical risk minimization
to the classical example  of the max-call option
on dividend paying stocks
as studied in the seminal paper by
Broadie \& Detemple \cite{BD},
and also in the book by Detemple  \cite{Dbook}.

We assume that
the stock price process
$S_t =(S^{(1)}_t,\ldots,S^{(d)}_t)\in \R_+^d$
is the state process and solves
$$
dS^{(i)}_t = S^{(i)}_t \Big ( (r-\delta_i)dt + \sigma_i dW^{(i)}_t\Big),
$$
where $\delta_i$ is the dividend rate.
In this example, the reward function  is given by,
$$
\varphi(t,s) = e^{-rt} (\alpha(s) -K)^+,
\quad
\text{where}
\quad\alpha(s) = \max\{ s_1,\ldots,s_d\}, \qquad
s \in \R_+^d.
$$
Hence, $X=S$, $\cX=\R_+^d$, $\eta=1$, and $\beta \equiv 0$.
In view of Theorem \ref{t.fb}, the existence of a
stopping boundary Assumption \ref{asm:star}
holds with $\alpha$
and $\Xi(s)=s/\alpha(s)$.
In two dimensions, the stopping
regions can be visualized effectively
as seen in the Figures \ref{t.2}, \ref{t.3}
reported from \cite{RST}.
These are stopping regions
in two space dimensions
obtained with
initial data $s_0=90$ and $s_0=100$
with parameters  given in \eqref{e.p2}.
Clearly the stopping boundary is independent of
the initial condition
and the below numerical results
verify it.  Also they are very
similar to those obtained in \cite{BD}.

\begin{figure}[H]
\centering
\includegraphics[height=1.75in,width=5in]{Figures/S090Boundary.pdf}
\caption{Stopping Boundaries with $s_0=90$, $K=100$,
$\sigma_i =0.2$, $r =0.05$, $\delta_i =0.1$, reported from \cite{RST}.}
\label{t.2}
\end{figure}

\vspace{-3mm}

\begin{figure}[H]
\begin{center}
\includegraphics[height=1.75in,width=5in]{Figures/S0100Boundary.pdf}
\caption{Stopping Boundaries with $s_0=100$, $K=100$,
$\sigma_i =0.2$, $r =0.05$, $\delta_i =0.1$, reported from \cite{RST}.}
\label{t.3}
\end{center}
\end{figure}

\subsection{Numerical results: Symmetric case}
\label{ss.num-mc-sy}
As in \cite{Becker1,Becker2}, we use the following parameter set:
\begin{equation}
\label{e.p2}
K=100,\ s_0=90,100,110,\ \sigma_i =0.2,\ r =0.05,\ \delta_i =0.1.
\end{equation}
We take the maturity to be $3$ years and $n=9$.  Thus, each
time interval corresponds to four months. Numerical experiments
for these parameters are reported
in the accompanying paper \cite{RST}.
They compare well to the results obtained in \cite{Becker1,Becker2}.

We also compute the initial price of this max option on
$d \in\{2,5,10,20,50\}$ assets.
Since we observe that the drift of the maximum stock
price process is of order $\ln(d)$, we
use important sampling with $\lambda_i=(r-\delta_i-\mu(d))/\sigma$
and $\mu(d)= - 0.01\ln(d)$ to ensure that
the maximum stock value does not cross the stopping boundary too soon.
Utilizing the fact that the stocks
are exchangeable in this example, we order the stock prices and
input the top five assets to the network.
 Table \ref{tab:symMaxOpt} summarizes the results.
 Despite a low cutoff, the  obtained prices are close or above the benchmark.
 Moreover, the runtime remains moderate as $d$ increases.
 The case $d=5$ is a classical example first
 introduced by \citet{BroadieGlasserman} and well-studied later in the literature.
Table \ref{tab:symMaxOpt} also contains the tight confidence intervals, obtained in \cite{AB}, \cite{BC}, and \cite{Becker1}
using a primal-dual method.

\begin{table}[h]
\centering
{\bf{Max-call option with parameters in \eqref{e.p2} and $s_0=100$}}
\vspace{5pt}

\centering
\begin{tabular}{ c c c c c l}
\hline  \hline
$d$ & Average Price& Highest Price & Runtime & Price in \cite{Becker2} &
Confidence Intervals
\\
\hline \hline
2   & 13.883 (0.009)  &  13.898 (0.008) & 29.1  & 13.901 & \ [13.892, 13.934] \; (\cite{AB}) \\
5   & 26.130 (0.010)  & 26.151 (0.009) & 31.8  & 26.147 &  \ [26.115, 26.164] \; (\cite{BC})\\
10   & 38.336 (0.015) &   38.355 (0.011)& 33.1 & 38.272 & \ [38.300, 38.367]  \; (\cite{Becker1})\\
20   & 51.728 (0.018) & 51.753 (0.011)  & 36.0 & 51.572  & \ [51.549, 51.803]
 \; (\cite{Becker1})\\
50  & 69.860 (0.012) & 69.881 (0.011) & 43.5 & 69.572   &  \ [69.560, 69.945]
 \; (\cite{Becker1})\\
\hline\\[-1em]
\end{tabular}
\caption{Max Option on $d \in\{5,10,20,50\}$ symmetric assets.
The second column is the  average  of ten experiments with its standard deviation in brackets. The third column is the highest price  among the realizations and its standard deviation in brackets.
The fourth column is the  average runtime (in seconds) per experiment for the training phase.}
\label{tab:symMaxOpt}
\end{table}

\subsection{Numerical results: Asymmetric case}
\label{ss.num-mc-nons}

We next investigate the Bermudan max option
with  asymmetric assets. We consider two
stocks with same volatility $\sigma$ but different dividend rates.
We use the same parameters for $K,r,\sigma$
as in  the previous subsection, and
$s_0=100$, $\delta = (5 \%,15 \%)$.

 We choose the parameter
 $\lambda$ in the Girsanov theorem,
 so as to make the assets symmetric under an equivalent measure.
 More precisely, we choose
 $\lambda = [(r+\mu(2))(1,1) - \delta]/{\sigma} \in \R^2$
 so that both assets have the same  drift $\mu(2) \in \R$ under $\Q_{\lambda}$,
 and $\mu(2)=-0.01\ln(2)$ is as in the previous subsection.
 The asymmetry is then captured by the Radon-Nikodym derivative
 appearing in the reward function.

The average and highest price after $10$ realizations
are given in Table \ref{tab:asymMaxOpt}. The last column  is a
benchmark price obtained from the Least Square
Monte Carlo (LSMC) algorithm \cite{LS}  with $2^{22}$  simulations.
As can be seen, our method in average computes a price comparable
to the benchmark, but may construct better stopping strategies
reported in the second column.
Notice also that having different dividends gives a premium
 over the symmetric case as shown in the first row of Table \ref{tab:symMaxOpt}.

\begin{table}[h!]
\centering
{\bf{Max-call option with asymmetric assets}
\vspace{5pt}

\centering
\begin{tabular}{ c c c c }
\hline  \hline
 Average Price& Highest Price & Runtime & LSMC Price\\
\hline \hline \\[-1em]
15.551 (0.014)  & 15.575 (0.011) & 30.9 & 15.558 (0.009) \\
\hline\\[-1em]
\end{tabular}
\caption{Max Option on $d=2$ asymmetric assets}
\label{tab:asymMaxOpt}
 \end{table}

\cref{fig:asymCall} displays the stopping and continuation
region over time.
In particular, we clearly see the assymmetry of the problem and the stopping region reflected in the figure.
For $i=1,2$, let the
the connected components of $\calS_t$
be given by,
$\calS^{(i)}_t := \left\{ s \in \cS_t
:  \alpha(s) = s_i \right\}$.
In \cref{fig:asymCall}, we observe that
the light red region $\cI(\calS^{(1)}_t)$
with $\cI(s_1,s_2) := (s_2,s_1)$,
is contained in $\calS^{(2)}_t$.
That is, if it is optimal to stop
at $S_t = (s_1,s_2)$,
then the same is true at   $\cI(S_t)=(s_2,s_1)$.
The structure of the assymmetry is therefore consistent with our expectations.

\begin{figure}[h]
\centering
\includegraphics[height=1.75in,width=5in]{Figures/Asym2MaxCall.pdf}
\caption{Stopping and continuation regions  for a max option with
$d=2$ assets and $T=3$.
The light red region is the reflection of the lower stopping region $\{ s \in \cS_t\
: \ s_1 \ge s_2\}$ through the diagonal.}
\label{fig:asymCall}
\end{figure}

\begin{remark}
\label{r.others}

For
max-call options on two assets,  \cite{Garcia}
also use a two-dimensional parametrization
of the stopping region given by
$$
\cS_t \approx \left\{ (s_1,s_2)\in \R_+^2\ :\
\alpha(s)=\max\{s_1,s_2\} \ge \theta_t^{1}, \ \text{and}\ \
|s_1-s_2| \ge \theta_t^{2}\ \right\}.
$$
While the above hypothesis class of regions
provides satisfactory numerical results,
they are restricted to two dimensions,
and only give a  simple polygonal approximation
of the actual stopping regions
that are more complex as can be seen in the
Figures \ref{t.2}, \ref{t.3}, \ref{fig:asymCall}.
\end{remark}

\section{Look-back Options}
\label{s.look-back}

\emph{Look-back options}  provides exposure to the
minimum or maximum values attained
during the tenure of the option.
There are several types of look-back options
that are commonly traded as summarized
in the Table \ref{tab:lookback}.   In these examples,
when the stock price
$S_t \in \R_+$ is modeled as a Markov process,
then we would have $X=(S, Z)$, where $Z_t$ is either
the \emph{running maximum} $\overline{S}_t$ or
\emph{running minimum} $\underline{S}_t$ given by,
$$
\overline{S}_t:=\max_{u\in[0,t]} S_u ,\qquad
\underline{S}_t:=\min_{u\in[0,t]} S_u .
$$
\begin{table}[H]
\centering
\begin{tabular}{lcrcc}
\hline \hline
\noalign{\vskip 0.5em}
Option & $e^{rt}\varphi(t,X_t)$ & $\eta$ & $\alpha(s,z)$ & $\beta(s,z)$ \\[0.5em] \hline \hline
\noalign{\vskip 0.5em}
Fixed Strike Call & $(\overline{S}_t - K)^{+} $ & $1$ & $z$ & $0$ \\[0.5em] \hline
\noalign{\vskip 0.5em}
Fixed Strike Put & $(K - \underline{S}_t)^{+} $ & $-1$ & $z$ & $0$ \\[0.5em] \hline
\noalign{\vskip 0.5em}
Floating Strike Call & $(S_t - \gamma \underline{S}_t)^+$  & $1$ & $s$ & $\gamma  z$ \\[0.5em] \hline
\noalign{\vskip 0.5em}
Floating Strike Put & $(\overline{S}_t  - S_t)^+$ & $-1$ & $s$ &  $\gamma z$ \\[0.3em] \hline
\end{tabular}
\caption{Different Lookback Options.}
\label{tab:lookback}
\end{table}

Clearly, $X=(S,Z)$ is a  Markov process with
$\cX = \{x=(s,z) \in \R_+^2\ | \ 0< s \le z\}$, when $Z_t=\overline{S}_t$,
or $\cX = \{x=(s,z) \in \R_+^2\ | \ 0< z \le s\}$ when $Z_t=\underline{S}_t$.
The scaling factor $\gamma$ appearing in the payoff of floating strike
options is introduced to reduce the price of these otherwise expensive contracts.
We therefore choose $\gamma \in [1,\infty)$ and $\gamma \in (0,1]$ for call and put options, respectively.
When $\gamma=1$, the floating strike look-back put (call) option delivers precisely the
drawdown (drawup) of the stock.
We refer the reader to \citet{DaiKwok}
for an introduction and discussion of American look-back claims.

\subsection{Fixed strike call}
\label{ss.fsc}

Although our setting
covers all look-backs, as an example we only consider
the fixed strike call.  Then,
$Z_t = \overline{S}_t$,
$\cX = \{x=(s,z) \in \R_+^2 \ | \  0 < s \le z\}$),
$\eta=1$, $\alpha(s,z) = z$.
By Theorem~\ref{t.fb}, a stopping boundary exist with
$\Xi(x) = (s/z,1)$.
Notice that
$1-(s/z) = (z-s)/z$ is
 the relative drawdown of the stock.

As in \cite{DaiKwok}, we use
\begin{equation}\label{e.p3}
T=1/2, \quad  r=2\%, \quad \delta =4\%, \quad  \sigma = 30 \%.
\end{equation}
To approximate the American option,
we consider a regular partition  with $n=200$ many exercise dates.
Moreover, we employ an even thinner  partition with $n'=800$ points
for our simulations.
This allows us to better approximate  the running maximum.
Importance sampling is not needed in this example.
We visualize the the stopping region in the $(s,y)$ plane, similar to Fig.~3 in \cite{DaiKwok},
where the authors use a finite difference scheme to compute the initial price.

\begin{table}[H]
  \centering

  \begin{tabular}{ccccc}
 \hline \hline

    Average Price  & Highest Price & Runtime   & European Price & Upper Bound\\
  \hline  \hline \\[-1em]

  16.827   (0.008)  &  16.844 (0.007) & 216.7  &  16.808 & 16.979 \\%[0.5em]
  \hline
\end{tabular}
  \caption{American look-back call with parameters in \eqref{e.p3} and $K=100$. }
\label{tab:resultLkbk}
  \end{table}

Table \ref{tab:resultLkbk} summarizes the result for $s_0=K=100$.
 First  column in that table is the average of ten experiments
 with the standard deviation in brackets. The second column is the highest price
 among the realizations and its standard deviation in brackets.
 The third column is the  average runtime (in seconds) per experiment for the training phase.
 The last column is the upper bound given by the forward value of the European price.
The European price, denoted by $v_e$ gives a
lower bound on the American price.
The  last  column is the  forward value of European price,
 $e^{rT}v_e$, is an upper bound for the American price.
 Since here $e^{rT} \approx 1.01$, the American option only offers little premium
 over its European counterparts. This is nevertheless  fortunate from a
 performance perspective as it gives a tight interval in which the American price
 must lie (see, e.g., \cite{Conze}).  We indeed obtain prices that are within the bounds.

We display the stopping and continuation region in \cref{fig:lkbkCall}.
Obtaining an overall accurate boundary turns out to be challenging,
as visiting all the pairs $(s,y) \in \calX$ is difficult, especially when $s \ll y$.
We effectively resolve this issue
by randomizing $s_0\in \calX$ and setting the initial running maximum to
$y_0 = K \vee s_0$. As in Fig. $3$ of \cite{DaiKwok}, we observe a flat boundary when $s\ll K$.
That is, the neural network correctly set $g(t,s/y;\theta)$ approximately equal to the strike when
the ratio $s/y$ is small, or equivalently, when the relative drawdown is large.

\begin{figure}[h]
\centering
\includegraphics[scale = 0.365]{Figures/Best.pdf}
\caption{Stopping boundaries for a look-back call  with parameters in \eqref{e.p3} and $K=100$.}
\label{fig:lkbkCall}
\end{figure}

\section{Conclusion}

We have developed an algorithm for the
computation of the ``graph-like'' stopping boundaries
separating the continuation and stopping regions
in optimal stopping.  While the method at the high-level is
empirical risk minimization, a relaxation
based on fuzzy boundaries motivated from phase-field
models for liquid-to-solid phase transitions is used.
Through numerical experiments, the method is shown
to be effective in high-dimensions.  The method has the potential
to incorporate market details like price impact, transaction costs,
and market restrictions.  It is also possible to apply the technique
to other models from financial economics
and other obstacle type problems as long as a
control representation is available.


\section{Introduction}
Coupled-spin systems on the honeycomb lattice have gained tremendous interest with the introduction of Kitaev's model Hamiltonian on this lattice that predicts a spin-liquid state for $S=1/2$ spins \cite{Kitaev2006,Mandal2020b}. Experimental results on the honeycomb lattice \cite{Smirnova2009} have also motivated a plethora of theoretical work on the two-dimensional (2D) Heisenberg $J_1-J_2-J_3$ model on this lattice \cite{Fouet2001,Mulder2010,Oitmaa2011,Bishop2012,Bishop2013b,Bishop2015,Sahoo2020}. It has been found that for $J_3=0$, the classical  $J_1-J_2$ Heisenberg model on the 2D honeycomb lattice exhibits a N\'eel-order ground state for $J_2 < J_1/6$ and various degenerate spiral-like states for $J_2 > J_1/6$  \cite{Mulder2010}.
The realization of the honeycomb model has been suggested  for a wide variety of trigonal $AM_2Pn_2$ ($A=$ Sr, Ca; $M =$ Cr, Mn, Fe, Co, Ni; $Pn$ = P, As, Sb, Bi) compounds with $P\bar{3}m1$ spacegroup. These compounds contain the transition-metal element in a corrugated honeycomb structure [Fig.~\ref{Fig:Structure}(b)] formed by two adjacent trigonal layers (sublattices) that are stacked as A/B-type [Fig.~\ref{Fig:Structure}(a)].  It is interesting to note that the transition-metal bilayer magnetic moments have no intervening binding atoms, as shown in Fig.~\ref{Fig:Structure}(a), so that the major magnetic coupling  between nearest neighbors is likely a direct $M-M$ coupling,  and couplings among further neighbors are likely due to $M-Pn-M$ superexchange.  Recently, it has been concluded that the superexchange within an Mn-$Pn$-Mn moiety becomes stronger as the atomic number of $Pn$ gets  smaller, i.e., stronger for P and weaker for Bi, for similar bond configurations \cite{Islam2020a}.
\begin{figure}[H]
\centering
\includegraphics[width=3 in]{structure_v2.pdf}
\caption{ (a) Chemical structure of $A$Mn$_2Pn_2$ ($A=$ Sr, Ca; $Pn$ = P, As, Sb, Bi) compounds showing the Mn trigonal bilayer without intervening elements. (b) The two trigonal Mn sublattices onto the $ab$-plane are projected, shown with red and green shades. The A and B layers are stacked with two atoms per unit cell (dotted rhombus shows the basal unit cell). The Mn bilayer forms a corrugated honeycomb lattice, where the nearest-neighbor interactions ($J_1$) and the next-nearest-neighbor interactions ($J_2$) are indicated. The magnetic structure shown is typical for the $A$Mn$_2Pn_2$ compounds with $Pn=$ As, Sb, Bi; for which $J_1 >> J_2$, and $J_1$ is antiferromagnetic and $J_2$ is ferromagnetic.}
\label{Fig:Structure}
\end{figure}

Here, we report neutron-diffraction results on single-crystal CaMn$_2$P$_2$, and determine its ground-state magnetic structure. This  compound belongs to the large family of Mn-based trigonal pnictides discussed above that exhibit varying magnetic properties with common characteristics.  The magnetic properties of CaMn$_2$P$_2$ and SrMn$_2$P$_2$ have been reported recently to undergo an unusual first-order antiferromagnetic (AFM) transition
at $T_{\rm N} =70(3)$ and 53(1)~K, respectively \cite{Sangeetha2021}. By contrast, the isostructural CaMn$_2$As$_2$ and SrMn$_2$As$_2$ compounds undergo second order AFM transitions \cite{Sangeetha2016a}.

Furthermore,  $^{31}$P NMR measurements  indicate that the magnetic structure of CaMn$_2$P$_2$ is commensurate with the lattice, whereas SrMn$_2$P$_2$ is incommensurate \cite{Sangeetha2021}. This is consistent with neutron-diffraction measurements of SrMn$_2$P$_2$ that indicate a complex and undetermined magnetic structure \cite{Brock1994}.

On the other hand, all the isostructural compounds that are based on As, Sb, or Bi undergo a simple N\'eel-order on the honeycomb lattice. In particular, neutron-diffraction measurements of SrMn$_2$As$_2$, SrMn$_2$Sb$_2$, CaMn$_2$Sb$_2$, and CaMn$_2$Bi$_2$ reveal N\'eel-ordering (see Fig.~\ref{Fig:Structure}) \cite{Das2017,Ratcliff2009,Bridges2009,Gibson2015}. For  CaMn$_2$Sb$_2$, it has been suggested that the moments are slightly canted towards the c-axis \cite{Bridges2009}.

It is important to note that all these compounds display strong two-dimensional (2D) magnetic fluctuations as manifested in the magnetic susceptibility ($\chi $) measurements that do not show Curie-Weiss behavior at temperatures much higher than $T_{\rm N}$ \cite{Das2017,Ratcliff2009,Bridges2009,Gibson2015,Sangeetha2021}.  In addition, the $\chi (T)$ with applied magnetic field along the $ab$-plane for all these compounds hardly show any anomaly at $T_{\rm N}$. This 2D behavior is also manifested in the magnetic order parameter in neutron-diffraction measurements of SrMn$_2$As$_2$ \cite{Das2017}. These characteristics indicate that the major in-plane coupling among NN ($J_1$) is AFM, common and strong in all these systems. This $J_1$ leads to 2D AFM correlations above $T_{\rm N}$, and a small interplanar coupling switches on the 3D ordering..
Inelastic neutron-scattering measurements analyzed by spin waves using the $J_1-J_2$ Heisenberg model determined  $J_1/J_2 \sim 1/6$ for CaMn$_2$Sb$_2$, suggesting this system is in the proximity of a tricritical point at which the Néel phase and two different spiral-like phases coexist \cite{Mcnally2015}.

\section{\label{Sec:ExpDet} Experimental Details and Methods}
Single crystals of CaMn$_2$P$_2$ were grown in Sn flux, as described previously \cite{Sangeetha2021}, and the crystal used in this study is from the same growth batch.

Single-crystal neutron-diffraction experiments were performed in a zero applied magnetic field using the TRIAX triple-axis spectrometer at the University of Missouri Research Reactor (MURR). An incident neutron beam of energy  14.7 meV was directed at the sample using pyrolytic graphite (PG) monochromator. A PG analyzer was used to reduce the background. Neutron wavelength harmonics were removed from the beam using PG filters placed before the monochromator and in between the sample and analyzer. Beam divergence was limited using collimators before the monochromator; between the monochromator and sample; sample and analyzer; and analyzer and detector of $60^\prime-60^\prime-40^\prime-40^\prime$, respectively. A 40~mg CaMn$_2$P$_2$ crystal was mounted on the cold tip of an Advanced Research Systems closed-cycle refrigerator with a base temperature of approximately 5 K. The crystal was mounted in the $(H,0,L)$ and $(H,H,L)$  scattering planes. The lattice parameters at base temperature are $a=4.096(1)$ and $c = 6.848(2)$~{\AA}. We also note that our sample consists of at least two crystals that are disoriented with respect to each other, as indicated in Fig. \ref{Fig:hh1}. Our diffraction patterns below also show Bragg reflections from the polycrystalline Al sample holder.

\section{\label{Sec:Neutron} Results and Discussion}
\begin{figure}
\centering
\includegraphics[width=3.1 in]{diffraction_data.pdf}
\caption{(a) Diffraction patterns along $(H,H,1)$ at $T=6$ and 100 K showing the emergence of a prominent peak at $H=1/6$ r.l.u.~(reciprocal lattice units). (b) The difference between the $(H,H,1)$ patterns at 6 K and 100 K showing the observed magnetic Bragg reflections in this direction are $(\eta,\eta,1)$ and $(1-\eta,1-\eta,1)$, where $\eta=1/6$. Al peaks (originating from the sample holder) in the difference show both positive and negative signals due to the thermal shift in peak positions. The peaks with asterisk marks originate from a twin of CaMn$_2$P$_2$ oriented in a different direction.}
\label{Fig:hh1}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.0 in]{diffraction_data_subtracted.pdf}
\caption{Difference between scans at low and high temperature, as indicated, along (a) $(-1/6,-1/6, L)$ showing peaks at integer values of $L$ between $-3$ and $3$ (except for $L=0$); (b) $(1/6,1/6, L)$ showing peaks at the integer values of $L$ between 1 and 3 (scans with negative $L$ were not accessible due to the experimental setup); (c) $(H,0,0)$ showing a weak peak at the nuclear $(1,0,0)$ and possibly another at $(2,0,0)$. More details on the temperature dependence of the $(1,0,L)$ peaks are provided in the Supplementary Material (SM) \cite{SM_CaMn2P2}. Signals from the Al sample holder are marked on the figures. As indicated, a minute inclusion of ferromagnetic MnP crystals gives rise to weak peaks.}
\label{Fig:diffraction2}
\end{figure}

Diffraction scans along the $(H,H,1)$ direction at $T=6$ and 100 K in Fig.~\ref{Fig:hh1}(a) show the emergence of a  prominent  peak at $H=1/6$ r.l.u. (reciprocal lattice units) at low temperatures.
As shown Fig.~\ref{Fig:hh1}(b), the difference between these scans at 6 and 100~K displays magnetic Bragg peaks at $(\eta,\eta,1)$ and $(1-\eta,1-\eta,1)$, where $\eta=1/6$.
Figure \ref{Fig:diffraction2}(a) shows the difference between scans at 6 and 100 K along $(-\eta,-\eta, L)$, indicating magnetic Bragg peaks at $L$ = $-3, -2, -1, 1, 2,$ and $3$. Figure \ref{Fig:diffraction2}(b) shows similar observations of magnetic Bragg peaks at $L$ = 1, 2, and 3 in the direction of $(\eta,\eta, L)$. Scans along $(H,H,0)$ do not show any newly-emerging peaks at low temperatures (not shown).
Figure \ref{Fig:diffraction2}(c) shows the difference of scans along $(H,0,0)$ at 6 and 100 K with a weak peak at the nuclear $(1,0,0)$  reflection and possibly another very weak one  at $(2,0,0)$. The other signals that have a negative intensity originate from the Al sample holder. Also, magnetic peaks from a small amount of MnP in the crystal are present in the scan.
As noted below, the temperature dependence of $(1,0,0)$ does not exhibit a transition at $T_{\rm N}$.
\begin{figure}
\centering
\includegraphics[width=2.9 in]{OPFig-eps-converted-to.pdf}
\caption{Integrated intensity as a function of temperature $T$ of the $(1/6,1/6,1)$ magnetic peak showing a sharp transition at $T=70$ K, consistent with specific-heat measurements in Ref.~\cite{Sangeetha2021}.  The dashed line is a guide to the eye. The data near $T_{\rm N}=70$ K indicates a first-order magnetic transition. The inset shows the $(1/6,1/6,1)$ peak at $T = 5$ and 100 K. The weak minimum below $T_{\rm N}$ at $\approx50$ K does not show up in  specific-heat measurements and is currently not understood.}
\label{Fig:OP}
\end{figure}
The temperature dependence of the integrated intensity of the ($\eta$,$\eta$,1) reflections in Fig.~\ref{Fig:OP} shows a very sharp transition at $T=70(1)$ K that coincides with a previous report that indicates a strong first-order magnetic phase transition at this temperature \cite{Sangeetha2021}.
The fact that the peak intensities of the $(\pm \eta,\pm \eta,L)$ reflections fall off for larger $L$, as expected from the magnetic form factor of Mn$^{2+}$, is further evidence that these newly observed Bragg peaks are magnetic in origin.
Below, we propose a magnetic structure that is consistent with all the experimental observations assuming the magnetic propagation vector $\boldsymbol{\tau}$ = $(\eta,\eta,0)$ r.l.u., with $\eta = 1/6$.

 \begin{figure*}
\centering
\includegraphics[width=7 in]{sfandstruct_v2.pdf}
\caption{(a) A compilation of the magnetic reflections observed in the $(H,H,L)$ planes, where the sizes of the spheres roughly reflect observed intensities. The regions in which the neutron-diffraction experiment was performed are shown by solid black lines. (b) The model structure is constructed by creating a $6\times6$ in-plane unit cell consisting of the corrugated honeycomb structure. The sites in green correspond to one trigonal layer (magnetic sublattice) and those in red to the one adjacent to it. The magnetic model has a space-group of $P_Ac$ \cite{bilbao} and is constructed by assigning an arbitrary moment direction at the origin (for instance, the red site at the lower-left corner) and then rotating each nearest neighbor by 60$^{\circ}$ in the same sublattice along the $a$- and $b$-direction.  The other sublattice (sites in green) spins are antiparallel to the one described above (sites in red); more details on the construction of the magnetic structure are provided in the SM \cite{SM_CaMn2P2}. (c) Squared structure factor calculation based on the magnetic structure shown in (b), agreeing with the experimental results shown in (a). Inspection of the model structure in (b) shows a global trend where along the long diagonal of a hexagon, all the spins are collinear, and the nearest neighbors are antiparallel.}
\label{Fig:ModelCalc}
\end{figure*}
To start with, the $(\eta,\eta,0)$ propagation vector indicates that the magnetic structure consists of a $6\times6$ nuclear basal unit cell.  Figure \ref{Fig:ModelCalc}(a) is a compilation of the magnetic reflections observed in the $(H,H,L)$ plane, where the sizes of the peaks roughly reflect observed intensities.
The model structure is constructed by creating a $6\times6$ in-plane nuclear unit cell that spans the corrugated honeycomb structure, i.e., the bilayer magnetic structure stacked along the $c$-axis [Fig.~\ref{Fig:ModelCalc}(b)]. The red sites correspond to one trigonal magnetic sublattice and the green sites to the other magnetic sublattice. The magnetic model is constructed by assigning an arbitrary moment direction at an origin, for instance, at the lower-left corner, and then successively rotating the spin on the nearest neighbors on the same sublattice clockwise by 60$^{\circ}$. The other sublattice is constructed similarly with anti-parallel spins with respect to the first sublattice. Note that along the $(1,0,0)$ and $(0,1,0)$ directions, the magnetic structure of each sublattice is a cycloid with a 60$^{\circ}$ turn angle. Thus, for each sublattice, the overall structure is a cycloid with propagation vector $(\eta, \eta, 0)$, with $\eta=1/6$.

To model the magnetic intensity, $I$, we use the following equation:
\begin{equation}
    I =C|f(Q)|^2\left|{\sum_{j=1}^{k}}{\rm e^{i{\bf Q\cdot r_j}}{\bf {\hat Q}}\times({\bf {\hat m_j}}\times{\bf \hat{Q}}})\right|^2,
\label{sfequation}
\end{equation}
where $C$ is a normalization factor, ${\bf Q}$ is the scattering vector, ${\bf r}_j$ and ${\bf \hat{m}}_j$ are the position of Mn moment and the unit vector of the magnetic moment, respectively.  $f(Q)$ is the magnetic form factor of Mn$^{2+}$. Using Eq.~(\ref{sfequation}), the calculated magnetic intensities are are shown in Fig.~\ref{Fig:ModelCalc}(c), in good agreement with the experimental results shown in \ref{Fig:ModelCalc}(a).

In the Supplementary Material (SM) \cite{SM_CaMn2P2}, we provide more details on the construction of the magnetic structure. As described in the SM, we construct the structure by stacking two trigonal sublattices to form the honeycomb structure. We also note that there are six possible stackings, three of which are inequivalent. We also provide structure-factor calculations for all six possible stackings and conclude that the structure shown in Fig \ref{Fig:ModelCalc}(b) with a magnetic space-group $P_Ac$ \cite{bilbao} is the correct one.

Our intensity calculation [Eq.~(\ref{sfequation})] allows us to estimate the average ordered magnetic moment, $\langle gS \rangle$,  in the model structure shown in Fig.~\ref{Fig:MagStruct}, where $g$ is the spectroscopic-splitting factor and $S$ is the spin quantum number. By comparing nuclear-peak intensities and their structure factors to the observed magnetic-peak intensities, we estimate $\langle gS \rangle = 4.2(5)$ ${\rm \mu_B}$, typical for Mn$^{2+}$ moments.

\begin{figure}[hb]
\centering
\includegraphics[width=3.5 in]{magstruct.pdf}
\caption{Comparison of our proposed magnetic structure with the two extreme cases of $J_2/J_1$. (a) N\'eel structure for $J_1 >> J_2$; (b) Cycloidal structure for (proposed) $\frac{1}{5}<J_2/J_1<\frac{1}{2}$; (c) Trigonal structure for $J_1 << J_2$. The spins on the red sites are not shown as we are assuming the correlation between the layers are insignificant, i.e. $J_1=0$.}
\label{Fig:MagStruct}
\end{figure}

To estimate the ratio of $J_2/J_1$, we examine various magnetic structures, as shown in Fig.~\ref{Fig:MagStruct} for (a) N\'eel ($J_1 >> J_2$), (b) Cycloidal, and (c) Trigonal ($J_1 << J_2$) structures. The ground state energy of these structures is calculated by assuming Heisenberg interactions acting between the nearest-neighbor and next-nearest-neighbor in-plane sites using the following formula:
\begin{equation}
    E = J_1 \sum_{NN} S_i \cdot S_j + J_2 \sum_{NNN} S_i \cdot S_j
    \label{Eq:structureEnergy}
\end{equation}
Consistent with our neutron-diffraction data, the interlayer coupling is ferromagnetic; this term is ignored in Eq.~(\ref{Eq:structureEnergy}). Using Eq.~(\ref{Eq:structureEnergy}) and assuming classical spins ($S=1$), we find the ground state energies of each structure to be:
\begin{equation}
    \begin{aligned}
        E_{N\acute{e}el} &= -3J_1 &+ 6J_2 \\
        E_{Cycloidal} &= -2J_1 &+ J_2 \\
        E_{Trigonal} &= &-3J_2
    \end{aligned}
    \label{Eq:eneries}
\end{equation}
The energy of the cycloidal structure is compared with energies of the other two extreme structures to determine the range of $J_2/J_1$ that give rise to the cycloidal structure. By setting up the two inequalities: (1) $E_{Cycloidal} - E_{N\acute{e}el} < 0$ and (2) $E_{Cycloidal} - E_{Trigonal} < 0$, we find the ratio of $J_2/J_1$ needed for the cycloidal structure is $\frac{1}{5}<J_2/J_1<\frac{1}{2}$. This calculation is not definitive and does not include the strong influence of magnetic fluctuations. Mean-filed calculations show that $J_2/J_1=1.73$ based on high-temperature measurements of the Curie-Weiss temperature \cite{Sangeetha2021} of the isostructural SrMn$_2$P$_2$ (see more details in SM).

\begin{figure}[ht]
\centering
\includegraphics[width=3 in]{otherStructures.png}
\caption{Two possible magnetic structures consistent with the experimentally-observed magnetic peaks. The two moments at the origin at the lower-left corner have the same spin orientation for both structures where the red and green sites are pointed along the $b$- and $a$- directions respectively. (a) shows the magnetic structure with space-group $P_C2$ \cite{bilbao} and is constructed by rotating the nearest neighbors in the red sublattice counterclockwise and those in the green sublattice clockwise by 60 degrees successively. (b) shows the magnetic structure with space-group $P_S\bar{1}$ \cite{bilbao} and is constructed by rotating the nearest neighbors in the red sublattice clockwise and green sublattice counterclockwise by 60 degrees successively. }
\label{Fig:OtherMagStructure}
\end{figure}

We emphasize that Fig.~\ref{Fig:ModelCalc}(b) with magnetic spacegroup $P_Ac$ \cite{bilbao} is not unique with our neutron-diffraction observation shown in Fig.~\ref{Fig:ModelCalc}(a). In fact, systematically searching through the Symmetry-Based Computational Tools for Magnetic Crystallography \cite{bilbao} we find two more magnetic structures, shown in Fig.~\ref{Fig:OtherMagStructure}, that are consistent with the diffraction measurements. In both structures, we start with an origin at the lower-left corner, where the spins on the red site and green site are oriented along the $b$- and $a$-directions, respectively. The magnetic structure in Fig.~\ref{Fig:OtherMagStructure}(a) has a magnetic space-group $P_C2$ \cite{bilbao} and is constructed by rotating the nearest neighbors in the red sublattice counterclockwise and green sublattice clockwise by 60 degrees successively. The other magnetic structure in Fig.~\ref{Fig:OtherMagStructure}(b) has a magnetic space-group of $P_S\bar{1}$ \cite{bilbao} and is constructed similarly by rotating the nearest neighbors in the red sublattice clockwise and green sublattice counterclockwise by 60 degrees successively. The energy associated with these two structures, according to Eq.~(\ref{Eq:structureEnergy}), is $-J_1+J_2$. Comparing with the energy associated with the structure mentioned in Fig.~\ref{Fig:ModelCalc}(b) and assuming the $J_1$ interaction to be AFM, the structures described in Fig.~\ref{Fig:OtherMagStructure} are energetically less favorable than the one in Fig.~\ref{Fig:ModelCalc}(b). Therefore, we propose that the magnetic structure with magnetic space-group $P_Ac$ shown in Fig.~\ref{Fig:ModelCalc}(b) is more favorable than the $P_C2$ and $P_S\bar{1}$ space-groups shown in Fig.~\ref{Fig:OtherMagStructure}. In the SM, we provide more temperature-dependent diffraction results from the crystal that are not correlated with the first-order phase transition at $T_{\rm N}=70$ K.

\section{Conclusions}
Using neutron-diffraction measurements, we find that below $T_{\rm N}=70$ K, CaMn$_2$P$_2$ undergoes a first-order antiferromagnetic transition with a magnetic unit cell that is a $6\times 6$ in-plane chemical unit cell. The average ordered magnetic moment $\langle gS \rangle = 4.2(5)$ $\mu_{\rm B}$. The integrated intensity of the major (1/6,1/6,1) magnetic peak vs. temperature shows an abrupt decrease in intensity at $T_{\rm N}$ that is a characteristic of a first-order phase transition. Our estimate of $\frac{1}{5}<J_2/J_1<\frac{1}{2}$ suggesting this system is very close to a spiral magnetic structure \cite{Mulder2010}. In general, although this and other isostructural pnictide systems possess corrugated-honeycomb structures, our results show that the corrugation may not affect the magnetic structure of these systems.

\acknowledgments
We thank Dr. Simon X M Riberolles for providing us with insights into the Bilbao Crystallographic Server. This research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering.  S.X.M.R.\ and B.U.\ are supported by the Center for Advancement of Topological Semimetals, an Energy Frontier Research Center funded by the U.S. Department of Energy Office of Science, Office of Basic Energy Sciences, through Ames Laboratory. Ames Laboratory is operated for the U.S. Department of Energy by Iowa State University under Contract No.~DE-AC02-07CH11358.