# Reinforced optimal control

C. Bayer  Christian Bayer
Weierstrass Institute
Berlin, Germany<EMAIL_ADDRESS>,  D. Belomestny  Denis
Belomestny
Duisburg-Essen University, Germany and National University Higher School of
Economics, Russia<EMAIL_ADDRESS>,  P. Hager  Paul Hager
Institut für Mathematik, Humboldt Universität zu Berlin
Germany<EMAIL_ADDRESS>,  P. Pigato  Paolo Pigato
Department of Economics and Finance
University of Rome Tor Vergata
Italy<EMAIL_ADDRESS>,  J. Schoenmakers  John Schoenmakers
Weierstrass Institute
Berlin, Germany<EMAIL_ADDRESS>and  V. Spokoiny  Vladimir
Spokoiny
Weierstrass Institute
Berlin, Germany<EMAIL_ADDRESS>
###### Abstract.

Least squares Monte Carlo methods are a popular numerical approximation method
for solving stochastic control problems. Based on dynamic programming, their
key feature is the approximation of the conditional expectation of future
rewards by linear least squares regression. Hence, the choice of basis
functions is crucial for the accuracy of the method. Earlier work by some of
us [Belomestny, Schoenmakers, Spokoiny, Zharkynbay. Commun. Math. Sci.,
18(1):109–121, 2020] proposes to _reinforce_ the basis functions in the case
of optimal stopping problems by already computed value functions for later
times, thereby considerably improving the accuracy with limited additional
computational cost. We extend the reinforced regression method to a general
class of stochastic control problems, while considerably improving the
method’s efficiency, as demonstrated by substantial numerical examples as well
as theoretical analysis.

###### Key words and phrases:

Reinforced regression, least squares Monte Carlo, stochastic optimal control

###### 2020 Mathematics Subject Classification:

91G20,93E24

C.B., P.H, P.P. J.S. and V.S. were supported by the MATH+ project AA4-2
Optimal control in energy markets using rough analysis and deep networks. The
authors are also thankful to the associate editor and anonymous referees for
there feedback and suggestions.

##  1\. Introduction

Stochastic control problems form an important class of stochastic optimization
problems that find applications in a wide variety of fields, see [Pha09] for
an overview. The general problem can be formulated as follows: How should a
decision-maker control a system with a stochastic component to maximize the
expected reward? In the theory of stochastic control, one distinguishes
between problems with continuous and discrete sets of possible control values.
While the first class of control problems contains, for example, energy
storage problems, the second one includes stopping and multiple stopping
problems. Furthermore one differentiates between discrete-time and continuous-
time optimal control problems. (Neither of these distinctions is fundamental:
for instance, many numerical methods will replace optimal control problems
with a continuous set of control values in continuous time by a surrogate
problem with discrete control values in discrete time. Moreover, many discrete
optimal control problems may well be analyzed as continuous ones, if the
number of possible control values or time-steps is finite but very high).

The range of applications of stochastic control problems is very wide.
Originally, optimal stochastic continuous control problems were inspired by
engineering problems in the continuous control of a dynamic system in the
presence of random noise, see [Åst12] and references therein. In the last
decades, problems in mathematical finance (portfolio optimization, options
with variable exercise possibilities) and economics inspired many new
developments, see [BR11] for some recent developments. Let us also mention a
closely connected area of reinforcement learning with plethora of applications
in robotics, data science, and engineering, see [SB18].

As a canonical general approach for solving a discrete-time optimal control
problem one may consider all possible future evolutions of the process at each
time that a control choice is to be made, see [Åst12]. This method is well
developed and may be effective in some special cases but for more general
problems such as optimal control of diffusion in high dimensions, this
approach is impractical. In [BKS09] a generic Monte Carlo approach combined
with linear regression was proposed and studied, see also [BS18] for an
overview. However, as an important disadvantage, there may be not enough
flexibility when modeling highly non-linear behavior of optimal value
functions. For instance, a regression based on higher-degree polynomials or
local polynomials (splines) may contain too many parameters and, therefore,
may over-fit the Monte Carlo sample or even prohibit parameter estimation
because the number of parameters is too large. As an alternative to the
polynomial bases, nonlinear approximation structures (e.g., artificial neural
networks) can be used instead (see, e.g. [HE16],[BCJ19] and [BCJW21]).

In [BSSZ20] a Monte Carlo based reinforced regression approach is developed
for building sparse regression models at each backward step of the dynamic
programming algorithm in the case of optimal stopping problems. In a nutshell,
the idea is to start with a generic set of basis functions, which is
systematically enlarged with highly problem-dependent additional functions.
The additional basis functions are constructed for the optimal stopping
problem at hand without using a fixed predefined finite dictionary. The new
basis functions are learned during the backward induction via incorporating
information from the preceding backward induction step. More specifically, the
(computed, hence approximate) value function at time $t_{i+1}$ is used as an
additional basis function at time $t_{i}$. Thereby, basis functions highly
specific to the problem at hand are constructed in a completely automatic way.
Indeed, the continuation function at time $t_{i}$ can often be observed to be
very close to the value function at time $t_{i+1}$, especially when the time-
step $t_{i+1}-t_{i}$ is small – alluding to continuity in time of the solution
to some continuous time version of the optimal stopping problem. [BSSZ20]
report that the reinforced basis leads to increased precision over the
starting set of basis functions, comparable to the standard regression
algorithm based on a substantially increased set of basis functions. This
improvement is obtained with a limited increase of the computational cost.

In this work, we carry over the approach of [BSSZ20] to a general class of
discrete-time optimal control problems including multiple stopping problems
(thus allowing pricing of swing options) and a gas storage problem. This
generalization turns out to be rather challenging as the complexity of using
the previously constructed value function in regression basis at each step of
the backward procedure becomes prohibitive when applying the original approach
of [BSSZ20]. We overcome this computational bottleneck by introducing a novel
version of the original reinforced regression algorithm where one uses a
hierarchy of fixed time-depth approximation of the optimal value function
instead of a full-depth approximation employed in [BSSZ20]. As a result, we
regain efficiency and are able to improve upon the standard linear regression
algorithm in terms of achievable precision for a given computational budget.

More precisely, we construct a hierarchy $v^{(i)}$, $i=0,\ldots,I$, of
(approximate) value functions with depth $I>0$. Here, $v^{(0)}$ denotes the
value functions obtained from the classical Monte Carlo regression algorithm.
The higher levels $v^{(i)}$ are computed by regression based on a set of basis
functions reinforced by the value function $v^{(i-1)}$ one level lower. This
way, the added computational cost incurred from reinforcing the basis can be
further decreased with minimal sacrifices of accuracy already for small values
of $I$. In fact, we propose two versions of the algorithm. In the first
version, the levels of the hierarchy of value functions are trained
consecutively, allowing for an adaptive choice of the depth $I$ of the
hierarchy. In the second version, all the levels are trained concurrently,
thereby improving the accuracy at each individual level. As a consequence, $I$
needs to be fixed in advance and cannot be chosen adaptively in the second
variant.

### Outline of the paper

In Section 2 we describe a rather general setting for discrete stochastic
control problems which we are going to use in this paper. The setting is based
on [GHW11]. We recall the reinforced regression algorithm for optimal stopping
problems by [BSSZ20] in detail in Section 3. There we also motivate the
hierarchical construction of the new reinforced regression algorithm as
restricted to the optimal stopping problem. The full algorithm – including
both variants – is introduced in Section 4. A detailed analysis of
computational costs is provided in Section 5. The next Section 6 provides a
detailed convergence analysis for the standard and reinforced regression
algorithms in the current setting. Extensive numerical examples including
optimal stopping problems, multiple stopping problems and a gas storage
optimization problem are provided in Section 7.

##  2\. Setting

First, we present a proper setting for the construction and analysis of
reinforced regression algorithms. The setting will be largely based on
[GHW11]. We will consider stochastic control problems in discrete time with
finite action sets. We note that extensions to continuous action sets are
certainly possible, but are left to future research.

We consider a filtration $\mathcal{F}_{j}$, $j=0,\ldots,J$, which is extended
by $\mathcal{F}_{-1}\coloneqq\set{\emptyset,\Omega}$,
$\mathcal{F}_{J+1}\coloneqq\mathcal{F}_{J}$ for convenience. Let $X$ be a
Markov process with values in $\mathcal{X}$ adapted to
$(\mathcal{F}_{j})_{j=0,\ldots,J}$. Note that we assume that the dynamics of
the underlying process $X$ does not depend on the control.

At time $0\leq j\leq J$ we are given a _control_ $Y_{j}$, which is
$\mathcal{F}_{j-1}$-measurable, and an $\mathcal{F}_{j}$-measurable _cash-
flow_ $Z_{j}=H_{j}(a,Y_{j},X_{j})$ for some deterministic, measurable function
$H_{j}$, where $a$ is an _action_ that we may choose at time $j$ in some
finite action space $\mathcal{K}$. Note that cash-flows may be positive or
negative. We assume that the control $Y_{j}$ takes values in a finite set
$\mathcal{L}$.

######  Remark 2.1.

The assumption that actions $a$ and controls $y$ take values in finite sets
$\mathcal{K}$ and $\mathcal{L}$, respectively, is a weaker assumption than it
may seem at first sight. Many important control problems naturally fall into
this class, see examples below. Even more importantly, it is a well-known fact
that many optimal control problems with genuinely continuous action and
control spaces have solutions of _bang-bang_ type, i.e., all optimal controls
consist of actions taken from a finite set, usually at the boundaries of the
(continuous) action sets. Hence, such control problems can effectively be
reduced to control problems with finite actions sets. Extensions of the
reinforced regression algorithm to infinite action spaces will be studied in
future work.

For a given value of the control $y\in\mathcal{L}$ and a given value $x$ of
the underlying process $X_{j}$, we are given a set of admissible actions

(2.1) $K_{j}(y,x)\subset\mathcal{K},\quad j=0,\ldots,J,$

i.e., $a$ is admissible iff $a\in K_{j}(x,y)$. Finally, if we apply $a\in
K_{j}(Y_{j},X_{j})$, then the control is updated by

(2.2)
$Y_{j+1}\coloneqq\varphi_{j+1}(a,Y_{j}),\quad\varphi_{j+1}:\mathcal{K}\times\mathcal{L}\to\mathcal{L}.$

Suppose that the control and the underlying state process take values $Y_{j}$
and $X_{j}$ at time $0\leq j\leq J$, respectively. For
$\mathbf{a}\coloneqq(a_{j},\ldots,a_{J})\in\mathcal{K}^{J-j+1}$ and
$j\leq\ell\leq J-1$, we define

(2.3)
$Y_{\ell+1}(\mathbf{a};j,Y_{j})\coloneqq\varphi_{\ell+1}(a_{\ell},Y_{\ell}(\mathbf{a};j,Y_{j})),\quad
Y_{j}(\mathbf{a};j,Y_{j})\coloneqq Y_{j},$

noting that $Y_{\ell}(\mathbf{a};j,Y_{j})$ only depends on
$a_{j},\ldots,a_{\ell-1}$. Additionally, we define
$\mathbb{F}_{j,J}(\mathcal{K})$ to be the set of
$(\mathcal{F}_{\ell})_{\ell=j}^{J}$-adapted processes taking values in
$\mathcal{K}$ indexed by $j,\ldots,J$. Clearly, if
$\mathbf{A}\coloneqq(A_{\ell})_{\ell=j}^{J}\in\mathbb{F}_{j,J}(\mathcal{K})$
and $Y_{j}\in\mathcal{F}_{j-1}$, then the process
$Y_{\cdot}(\mathbf{A};j,Y_{j})$ is previsible. The set of _admissible
strategies or admissible policies_ $\mathcal{A}_{j}$ is defined as follows:

(2.4) $\mathcal{A}_{j}(Y_{j},X_{\geq
j})\coloneqq\Bigl{\\{}\mathbf{A}=(A_{\ell})_{\ell=j}^{J}\in\mathbb{F}_{j,J}(\mathcal{K})\,\Big{|}\,A_{\ell}\in
K_{\ell}(Y_{\ell}(\mathbf{A};j,Y_{j}),X_{\ell}),\quad\ell=j,\ldots,J\Bigr{\\}}.$

Now the central issue is the optimal control problem

(2.5)
$V_{j}\coloneqq\sup_{\mathbf{A}=(A_{\ell})_{\ell=j}^{J}\in\mathcal{A}_{j}(Y_{j},X_{\geq
j})}\mathsf{E}_{j}\left[\sum_{\ell=j}^{J}H_{\ell}(A_{\ell},Y_{\ell}(\mathbf{A};j,Y_{j}),X_{\ell})\right],$

at a generic time $0\leq j\leq J,$ where $\mathsf{E}_{j}$ denotes the
conditional expectation w.r.t. $\mathcal{F}_{j}$.

Taking advantage of the Markov property, we introduce the notation
$\mathcal{A}_{j}(y,x):=\mathcal{A}_{j}(y,X^{x}_{\geq j}),$ where $X^{j,x}$
denotes the Markov process $X$ conditioned on $X_{j}=x$, and is defined for
$j\leq\ell\leq J$. We may then define the _value function_ as

(2.6)
$v^{*}_{j}(y,x)\coloneqq\sup_{\mathbf{A}=(A_{\ell})_{\ell=j}^{J}\in\mathcal{A}_{j}(y,x)}\mathsf{E}\left[\sum_{\ell=j}^{J}H_{\ell}\left(A_{\ell},Y_{\ell}(\mathbf{A};j,y),X_{\ell}^{j,x}\right)\right],$

which satisfies the _dynamic programming principle_ :

(2.7) $v^{*}_{j}(y,x)=\sup_{a\in
K_{j}(y,x)}\left(H_{j}(a,y,x)+\mathsf{E}\left[v^{*}_{j+1}(\varphi_{j+1}(a,y),X_{j+1}^{j,x})\right]\right).$

for $j=0,\ldots,J$ (with $v^{*}_{J+1}(y,x)\coloneqq 0$).

Let us now give a few examples for classical stopping and control problems
which fall into the above setup.

######  Example 2.2.

For a single optimal stopping problem with payoff $g_{j}\geq 0$ at time $j$,
the set of possible control values is $\mathcal{L}=\set{0,1}$, where a control
state $y$ denotes the number of remaining exercise opportunities. The action
$a$ takes the value $1$ if we stop at the current time and $0$ otherwise.
Hence, we have

$K_{j}(y,x)=K(y)\coloneqq\begin{cases}\set{0,1},&y=1,\\\
\set{0},&y=0,\end{cases}$

implying that $\mathcal{K}=\set{0,1}$. The cash-flow is defined by

$H_{j}(a,y,x)\coloneqq a\,g_{j}(x),$

independent of the value of the control $y$. Finally, the update function of
the control is defined by $\varphi_{j+1}(a,y)\coloneqq\max(y-a,0).$ Note that
the value function $v_{j}^{\ast}(0,\cdot)\equiv 0$, and, hence, the optimal
stopping literature usually only considers $(j,x)\mapsto v_{j}^{\ast}(1,x)$.

######  Example 2.3.

Let us now suppose that we have a multiple stopping problem with
$L\in\mathbb{N}$ exercise rights. Again, the control state $y$ signifies the
remaining exercise opportunities, leading to $\mathcal{L}=\set{0,1,\ldots,L}$.
The admissible action set is now defined as

$K_{j}(y,x)=K(y)\coloneqq\begin{cases}\set{0,1},&y\geq 1,\\\
\set{0},&y=0.\end{cases}$

Again, $\mathcal{K}=\set{0,1}$. The cash-flow $H_{j}$ and the update function
$\varphi_{j+1}$ are defined as in Example 2.2.

######  Example 2.4.

Consider a simple gas storage problem: given $N\in\mathbb{N}$ and
$\Delta=1/N$, we assume that the volume of gas in a storage can only be
increased and decreased by a fraction $\Delta$ over a given time increment.
Let the control $y$ denote the status (fill level) of the storage at time $j$.
Hence, we define $\mathcal{L}\coloneqq\set{0,\Delta,2\Delta,\ldots,1}.$ At
time $j$, we may either sell $\Delta$ (volume of gas; $a=-1$), buy $\Delta$
($a=+1$) – at the current market price $X_{j}$ – or do nothing ($a=0$). Hence,
the admissible policy set is

$K_{j}(y)\coloneqq\begin{cases}\set{0,1},&y=0,\\\ \set{-1,0,1},&\Delta\leq
y\leq 1-\Delta,\\\ \set{-1,0},&y=1,\end{cases}$

with $\mathcal{K}=\set{-1,0,1}$, while the cash-flow is given by

$H_{j}(a,y,x)\coloneqq-a\Delta x.$

The update function in given by $\varphi_{j+1}(a,y)\coloneqq((y+a\Delta)\wedge
1)\vee 0.$

######  Remark 2.5.

While we do not allow the actions to have an effect on the dynamics of the
state process $X$, a large class of more general control problems could be
incorporated by a simple modification of our setting. If we allow updates of
the control variable $y$ to depend on the state $x$ as well as on the previous
control and the action, i.e., $Y_{j+1}=\varphi_{j}(a,Y_{j},X_{j})$, then our
theoretical analysis remains intact. However, it now becomes possible to
control the dynamics of the state process $X$, provided that the law of the
controlled process remains absolutely continuous w.r.t. the law of the
original process $(X_{j})_{j=0,\ldots,J}$. We refer to the discussion of the
optimal liquidation example in [GHW11, Section 2] for more details. Note,
however, allowing $Y$ to depend on $X$ in such a way might require us to use
regression in $(x,y)$ rather than just $x$ for most practical problems.

##  3\. Reinforced regression for optimal stopping

In this section, we recall the standard regression algorithm as well as the
reinforced regression algorithm introduced in [BSSZ20] for optimal stopping
problems. We will point out the drawbacks of the latter algorithm for more
general control problems, and propose and motivate several modifications.
However, for the purpose of a clear illustration, we will restrict ourselves
in this section to the optimal stopping case.

Let us recall the optimal stopping setup from Example 2.2 and denote by
$v^{*}_{j}(x)$ the value function at $j\in\\{0,...,J\\}$ evaluated in
$x\in\mathbb{R}^{d}$ and $y=1$. Further recall that the dynamic programming
principle is given by

$v^{*}_{j}(x)=\max(g_{j}(x),c^{*}_{j}(x)),\;\;0\leq j\leq J-1,\quad
v^{*}_{J}(x)=g_{J}(x),\quad x\in\mathbb{R}^{d},$

where the continuation function is given by
$c^{*}_{j}(x)=\mathsf{E}_{j}[v^{*}_{j+1}(X^{j,x}_{j+1})]$. Fix a set of basis
function $\\{\psi_{1},...,\psi_{K}\\}$ with
$\psi_{k}:\mathbb{R}^{d}\to\mathbb{R}$, $k=1,\ldots,K$, and sample
trajectories $(X^{(m)}_{j})_{0\leq j\leq J,1\leq m\leq M}$ from the underlying
Markov chain, i.e., $(X^{(m)}_{j})_{0\leq j\leq J}$ are i.i.d. samples from
the distribution of $(X_{j})_{0\leq j\leq J}$, $m=1,\ldots,M$. Then the
regression method due to Tsitsiklis-van Roy [TVR01], which we will refer to as
the _standard regression method_ , inductively constructs an approximation
$v=(v_{j})_{j=0,...,J}$ to the value function $v^{*}$ as follows: For $j=J$
initialize $v_{J}=g_{J}$. For $j\in\\{J-1,...,0\\}$ set

(3.1) $\displaystyle v_{j}(x):=\max(g_{j}(x),c_{j}(x)),\quad
c_{j}(x)=\sum_{k=1}^{K}\gamma_{j,k}\psi_{k}(x),$

where the regression coefficients are given by the solution to the least
squares problem

(3.2)
$\displaystyle\gamma_{j,1},...,\gamma_{j,K}:=\underset{\gamma_{1},...,\gamma_{K}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}v_{j+1}(X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})\bigg{|}^{2}.$

The procedure is illustrated in Figure 1. Note that the costs of this
algorithm are of the order $M\cdot J\cdot K^{2}$ (see, e.g., [BSSZ20] or
Section 5).

Figure 1. Illustration of standard regression approach due to Tsitsiklis-van
Roy [TVR01]. The solid arrows indicate the dependencies in the (feed forward)
evaluation of $c_{j}$ and $v_{j}$ in (3.1). The dashed arrows start from the
regression data $v_{j+1}$ and symbolize the regression procedure (3.2).

One problem of the standard regression algorithm is that its performance
strongly depends on the choice of basis functions. Indeed, while standard
classes such as polynomials or splines usually form the backbone of the
construction of basis functions, practitioners usually add customized basis
functions, for instance the payoff function $g_{j}$ and some functionals
applied to it.

As a more systematic approach, the authors of [BSSZ20] proposed a _reinforced
regression algorithm_. In this procedure the regression basis at each step of
the backward induction is reinforced with the approximate value function from
the previous step of the induction. The approximate continuation function at
$j\in\\{0,...,J-1\\}$ is then given by

$c_{j}(x)\coloneqq\sum_{k=1}^{K}\gamma_{j,k}\psi_{k}(x)+\gamma_{j,K+1}v_{j+1}(x),$

where the regression coefficient are the solution to the least squares problem

$\gamma_{j,1},...,\gamma_{j,K+1}\coloneqq\underset{\gamma_{1},...,\gamma_{K+1}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}v_{j+1}(X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})-\gamma_{K+1}v_{j+1}(X^{(m)}_{j})\bigg{|}^{2}.$

Note that this procedure induces a recursion whenever an approximate value
function is evaluated: in order to evaluate $v_{j}(x)$ we need to evaluate
$c_{j}(x)$, which in turn requires an evaluation of $v_{j+1}(x)$ and so forth,
until $v_{J}(x)=g_{J}(x)$ terminates the recursion. Figure 2 illustrates this
procedure. The costs of the reinforced regression method are proportional to
$M\cdot J\cdot K^{2}+M\cdot J^{2}\cdot K$ (see [BSSZ20]).

Figure 2. Illustration of the reinforced regression approach. Evaluation of
$v_{j}$ in the reinforced regression algorithm leads to a recursion with $J-j$
steps.

Despite the increased computational cost compared to the standard regression
algorithm with the same set of basis functions $\psi_{1},\ldots,\psi_{K}$, the
reinforced regression algorithm can improve the overall computational cost for
a fixed error tolerance drastically. As a rule of thumb, [BSSZ20] report that
the reinforced regression algorithm with a standard basis consisting of
polynomials of a given degree leads to similar accuracy as the standard
regression algorithm based on polynomials of one degree higher. In particular,
the reinforced regression algorithm already outperforms the standard
regression algorithm for small dimensions $d>1$, as long as the number $J$ of
time-steps is not too large.

A direct generalization of the reinforced regression algorithm to more general
control problems is certainly possible. The main difference to the optimal
stopping problem is that for fixed time $j$ we have to choose from many
potential candidates to reinforce with, namely any $v_{j+1}(y,\cdot),$
$y\in\mathcal{L}$ is a candidate. Additionally, the dynamic programming
principle (2.7) now entails a possibly non-trivial optimization problem in
terms of the policy $a$. Especially the second point makes the recursion at
the heart of the reinforced regression algorithm untenable for general control
problems.

One solutions immediately comes to mind: If performing the recursion all the
way to terminal time $J$ is too costly, why not truncate at a certain
recursion depth? This idea is, in principle, sound, and is the basis of the
adaptations suggested below. However, some care is needed in the
implementation of this idea. Indeed, if “truncation” simply were to mean
“replace the reinforcing basis functions by $0$ after a certain truncation
step”, this would introduce a structural error in the procedure, as regression
coefficients formerly computed in the presence of these basis functions would
suddenly be incorrect. Instead, we propose to compute a _hierarchy_ of
reinforced regression solutions, corresponding to different “cut-off depths”
of the recursion. This way, we can make sure that the coefficients are always
consistent, that is, an error as mentioned above can be avoided. We introduce
two versions, which both adhere to the same general idea, but differ in an
important implementation detail.

The _hierarchical reinforced regression algorithm A_ iteratively constructs
approximations $(v^{(i)})_{i=0,1,...}$ to the true value function as follows:
For $i=0$ we construct $(v^{(0)}_{j})_{0\leq j\leq J}$ using the standard
regression method described above. Then for any $i\geq 1$, given that
$v^{(l)}$ is already constructed for $0\leq l\leq i-1$, define $v^{(i)}$ with
the usual backwards induction, where the regression basis at step
$j\in\\{J-1,...,0\\}$ is reinforced with $v^{(i-1)}_{j}$. The approximate
continuation function of the $i^{th}$ iteration is given by

(3.3)
$c^{(i)}_{j}(x)\coloneqq\sum_{k=1}^{K}\gamma^{(i)}_{j,k}\psi_{k}(x)+\gamma_{j,K+1}^{(i)}v^{(i-1)}_{j+1}(x),$

where the regression coefficients are the solutions to the least squares
problem

$\gamma^{(i)}_{j,1},...,\gamma^{(i)}_{j,K+1}\coloneqq\underset{\gamma_{1},...,\gamma_{K+1}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}v^{(i)}_{j+1}(X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})-\gamma_{K+1}v^{(i-1)}_{j+1}(X^{(m)}_{j})\bigg{|}^{2}.$

The procedure may be stopped after a fixed number of iterations, or using an
adaptive criterion. An illustration of the method can be found in Figure 3.
Note that the recursion that is started when evaluating $v^{(i)}_{j}(x)$
always terminates after at most $i$ steps in the evaluation of
$v^{(0)}_{j+i}(x)$ for $i\leq J-j$ or in $v^{(i-J-j)}_{J}(x)=g_{J}(x)$ for
$J-j\leq i$.

Figure 3. Illustration of the hierarchical reinforced regression algorithm A,
for three iterations. In the lower right part of the diagram, the vertical
lines indicate the equality $v^{(i)}_{j}\equiv v^{(l)}_{j}$ for $J-j\leq i$.

For a fixed number of iterations $i\in\\{0,...,I\\}$ we can modify the
structure of the previous method so that the primary iteration is the
backwards induction over $j\in\\{J,J-1,...,0\\}$ and the secondary iteration
is over $i\in\\{0,...,I\\}$. In this case we can further modify the algorithm
by using $v^{(I)}_{j+1}$ as the regression target for the continuation
functions $c^{(i)}_{j}$ for all $i\in\\{0,...,I\\}$. We name the resulting
algorithm the _hierarchical reinforced regression algorithm B_. The
approximate continuation function at step $j$ and iteration $i$ is then still
given by (3.3) and the least squares problem is given by

$\gamma^{(i)}_{j,1},...,\gamma^{(i)}_{j,K+1}\coloneqq\underset{\gamma_{1},...,\gamma_{K+1}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}v^{(I)}_{j+1}(X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})-\gamma_{K+1}v^{(i-1)}_{j+1}(X^{(m)}_{j})\bigg{|}^{2}.$

Also in this algorithm, the recursion that is started when evaluating
$v^{(I)}_{j}$ stops after at most $I$ steps. The costs of the algorithm are
discussed in Section 5.

##  4\. Iterated reinforced regression for optimal control

Following the ideas and motivations of Section 3 we now present hierarchical
reinforced regression algorithms for optimal control based on the Bellman
equation (2.7). The algorithms are based on $M$ sample trajectories
$(X^{(m)}_{j})_{j=0,...,J,m=1,...,M}$ from the underlying Markov chain $X$.
and some initial set $\\{\psi_{1},...,\psi_{K}\\}$ of basis functions
$\psi_{i}:\mathbb{R}^{d}\to\mathbb{R}$. For each $y\in\mathcal{L}$ we will
define a subset $\mathcal{L}^{y}\subset\mathcal{L}$ of cardinality
$R^{y}\coloneqq|\mathcal{L}^{y}|$ and reinforce the basis
$\\{\psi_{1},\ldots,\psi_{K}\\}$ by
$\\{v_{j+1}(z,\cdot)|z\in\mathcal{L}^{y}\\}$. The respective algorithms
iteratively construct sequences of approximations to the value function

$v^{(i)}=(v^{(i)}_{j})_{j=0,...,J}\quad\text{with}\quad
v^{(i)}_{j}:\mathcal{L}\times\mathbb{R}^{d}\to\mathbb{R},$

for $i=\\{0,1,...\\}$ until the iteration is terminated.

###  4.1. Hierarchical reinforced regression algorithm A

For $i=0$ construct $v^{(0)}$ using the standard regression method inductively
as follows: At the terminal time $J$ initialize $v^{(0)}_{J}\coloneqq v_{J}$
where

(4.1) $v_{J}(y,x)=\max_{a\in K_{J}(y,x)}H_{J}(a,y,x),\quad\text{for all}\quad
y\in\mathcal{L},\;x\in\mathbb{R}^{d}.$

For a $j\in\\{0,...,J-1\\}$, assume that $v^{(0)}_{l}$ is already constructed
for all $l\in\\{j+1,...,J\\}$. Then for each $y\in\mathcal{L}$ define the
regression coefficients by solving the following least squares problem

(4.2)
$\gamma^{(0),y}_{j,1},...,\gamma^{(0),y}_{j,K}\;\coloneqq\;\underset{\gamma_{1},...,\gamma_{K}}{\operatornamewithlimits{arg\,min}}{\sum_{m=1}^{M}\left|v^{(0)}_{j+1}(y,X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})\right|^{2}}.$

Next define the continuation function by

(4.3)
$c^{(0)}_{j}(y,x)\coloneqq\sum_{k=1}^{K}\gamma_{j,k}^{(0),y}\psi_{k}(x),\quad\text{for
all}\quad y\in\mathcal{L},\;x\in\mathbb{R}^{d}$

and the approximate value function $v^{(0)}_{j}$ through the dynamic
programming principle

(4.4) $v^{(0)}_{j}(y,x)\coloneqq\max_{a\in
K_{j}(y,x)}\Big{(}H_{j}(a,y,x)+c^{(0)}_{j}(\varphi_{j}(a,y),x)\Big{)}\quad\text{for
all}\quad y\in\mathcal{L},\;x\in\mathbb{R}^{d}.$

Given the approximation $v^{(i)}$ for some $i\geq 0$ we construct a new
approximation $v^{(i+1)}$ using reinforced regression inductively as follows:
Initialize at the terminal time $v^{(i+1)}_{J}\coloneqq v_{J}$. For
$j\in\\{0,...,J-1\\}$ assume that $v^{(i+1)}_{l}$ is already constructed for
$l\in\\{j+1,...,J\\}$. Then for each $y\in\mathcal{L}$ define the regression
coefficients by solving the following least squares problem

(4.5)
$\begin{split}\gamma^{(i+1),y}_{j,1},...,\gamma^{(i+1),y}_{j,K+R^{y}}\coloneqq\underset{\gamma_{1},...,\gamma_{K+R^{y}}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}\;v^{(i+1)}_{j+1}(y,X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})\\\
-\sum_{k=1}^{R^{y}}\gamma_{K+k}v^{(i)}_{j+1}(y_{k},X^{(m)}_{j})\;\bigg{|}^{2},\end{split}$

where $\\{y_{k}\\}_{k=1,...,R^{y}}=\mathcal{L}^{y}$, and define the
continuation function $c^{(i+1)}_{j}$ by

(4.6)
$c^{(i+1)}_{j}(y,x)\coloneqq\sum_{k=1}^{K}\gamma^{(i+1),y}_{j,k}\psi_{k}(x)+\sum_{k=1}^{R^{y}}\gamma^{(i+1),y}_{j,K+k}v^{(i)}_{j+1}(y_{k},x),$

for all $y\in\mathcal{L}$ and $x\in\mathbb{R}^{d}$. Finally define the
approximation $v_{j}^{(i+1)}$ through the dynamic programming principle by

(4.7) $v^{(i+1)}_{j}(y,x)\coloneqq\max_{a\in
K_{j}(y,x)}\Big{(}H_{j}(a,y,x)+c^{(i+1)}_{j}(\varphi_{j}(a,y),x)\Big{)},$

for all $y\in\mathcal{L}$ and $x\in\mathbb{R}^{d}$.

The iteration over $i\in\\{0,1,...\\}$ can be terminated after
$I\in\mathbb{N}$ steps, yielding $v^{(I)}$ as an approximation to the true
value function. Alternatively one can introduce an adaptive termination
criterion, for example by comparing the relative change in the error of the
least squares problem (4.5), terminating after the change falls under a given
threshold.

######  Remark 4.1.

Recall that in the initialization we have $v^{(i)}_{J}=v^{(0)}_{J}$ for all
$i\in\\{1,...,I\\}$. It then follows inductively that

(4.8) $v^{(i)}_{j}\equiv v^{(l)}_{j},\quad\text{ for all }\quad J-j\leq i\leq
I,\quad l\geq i.$

This identity can be used to reduce the costs of the algorithm, since the
regression problem only needs to be solved for all $(j,i)$ with $0\leq j\leq
J-1$ and $0\leq i\leq(J-j)\wedge I$.

######  Remark 4.2.

More general or other forms of reinforced basis functions are certainly
possible. The essential point is that they are based on the regression result
from the preceding step in the backwards induction and the preceding
iteration. Our specific choice may be seen as a natural primal choice. We left
flexibility in the choice of the sets $\mathcal{L}^{y}$, for which depending
on the cardinality of the set $\mathcal{L}$, possible choices are the trivial
$\mathcal{L}^{y}=\mathcal{L}$ and $\mathcal{L}^{y}=\\{y\\}$, or
$\mathcal{L}^{y}=\mathcal{L}^{\prime}$ for some set $\mathcal{L}^{\prime}$
independent of $y$, or more elaborately
$\mathcal{L}^{y}_{j}=\\{\varphi(a,y)\;|\;a\in K_{j}(y,x_{j})\\}$ for some
$x_{j}\in\mathbb{R}^{d}$. Note that the use of a step dependent set
$\mathcal{L}_{j}^{y}$ in the above method is straightforward.

###  4.2. Hierarchical reinforced regression algorithm B

Note that (4.2) and (4.5) are based on the approximate value functions
$v^{(0)}_{j+1}$ and $v^{(i+1)}_{j+1}$, respectively, even though the more
accurate approximation $v^{(I)}_{j+1}$ is already available at this point.
Hence, we can potentially improve the algorithm’s accuracy by always
considering the most accurate approximation of the value function $v_{j+1}$ in
the Bellman equation.

Fix a number of iterations $I\in\mathbb{N}$ and initialize the approximate
value functions at the terminal time by $v^{(i)}_{J}\equiv v_{J}$ for all
$i\in\\{0,...,I\\}$, where $v_{J}$ is given by (4.1). The approximate value
functions at times previous to $J$ are defined inductively as follows:

Let $j\in\\{0,...,J-1\\}$ and assume that $v^{(i)}_{l}$ is already defined for
all $l\in\\{j+1,...,J\\}$ and $i\in\\{0,...,I\\}$. For $i=0$ and each
$y\in\mathcal{L}$ determine the coefficients for the regression basis by
solving the least squares problem

(4.9)
$\gamma^{(0),y}_{j,1},...,\gamma^{(0),y}_{j,K}\;\coloneqq\;\underset{\gamma_{1},...,\gamma_{K}}{\operatornamewithlimits{arg\,min}}{\sum_{m=1}^{M}\left|v^{(I)}_{j+1}(y,X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})\right|^{2}}$

and define the approximate continuation function $c^{(0)}$ by (4.3). For
$i\in\\{1,...,I\\}$ and each $y\in\mathcal{L}$ determine the regression
coefficients by solving the least squares problem

(4.10)
$\gamma^{(i),y}_{j,1},...,\gamma^{(i),y}_{j,K+R^{y}}\coloneqq\underset{\gamma_{1},...,\gamma_{K+R^{y}}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}\;v^{(I)}_{j+1}(y,X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})\\\
-\sum_{k=1}^{R^{y}}\gamma_{K+k}v^{(i-1)}_{j+1}(y_{k},X^{(m)}_{j})\;\bigg{|}^{2},$

where $\\{y_{k}\\}_{k=1,...,R^{y}}=\mathcal{L}^{y}$, and define the
continuation function $c^{(i)}_{j}$ by (4.6).

Finally, define the approximation to the value function $v^{(i)}_{j}$ for all
$i=\\{0,...,I\\}$ by (4.7). After ending the backwards induction use
$(v^{(I)}_{j})_{j=0,...,J}$ as an approximation to the true value function.

######  Remark 4.3.

Note that the identity (4.8) also holds for the above algorithm. Moreover,
since we are only interested in $v^{(I)}$, we can discard the computation of
$c^{(i)}_{j}$ and $v^{(i)}_{j}$ for all $0\leq j+i\leq I-1$, since they do not
contribute to the construction of $v^{(I)}$. The least squares problem then
only needs to be solved for $(j,i)\in\\{0,...,J-1\\}\times\\{0,...,I\\}$ with
$0\leq j+i\leq I-1$ and $0\leq i\leq(J-j)\wedge I.$

######  Remark 4.4.

Choosing the number of iterations $I=J$ we then have from the previous remark
that only the value functions on the diagonal $j=i$ need to be constructed. In
this case, denote $v_{j}=v^{(j)}_{j}$, $c_{j}=c^{(j)}_{j}$ etc., and observe
that the least squares problem which is solved in each step
$j\in\\{J-1,...,0\\}$ of the backwards induction is given by

$\gamma^{y}_{j,1},...,\gamma^{y}_{j,K+R^{y}}\coloneqq\underset{\gamma_{1},...,\gamma_{K+1}}{\operatornamewithlimits{arg\,min}}\sum_{m=1}^{M}\bigg{|}\;v_{j+1}(y,X^{(m)}_{j+1})-\sum_{k=1}^{K}\gamma_{k}\psi_{k}(X^{(m)}_{j})-\sum_{k=1}^{R^{y}}\gamma_{K+k}v_{j+1}(y_{k},X^{(m)}_{j})\;\bigg{|}^{2}.$

where $\\{y_{k}\\}_{k=1,...,R^{y}}=\mathcal{L}^{y}$. Hence, for $I=J$ the
above algorithm represents a direct extension of the _reinforced regression
algorithm_ in [BSSZ20] from optimal stopping to optimal control problems.

##  5\. Computational cost

We study the computational work of the modified reinforced regression
algorithm of Section 4.2. In what follows, the following operations are
considered to be performed at _constant_ cost:

  * •

Multiplications, additions and other primitive operations at cost $c_{\ast}$;

  * •

Simulation from the distribution of the Markov process $X_{j}$ at cost
$c_{X}$;

  * •

Evaluation of the standard basis functions $\psi_{i}$ or of the payoff $H_{j}$
at cost $c_{f}$;

We furthermore introduce the following notations:

  * •

We set $R\coloneqq\max_{y\in\mathcal{L}}R^{y}$.

  * •

The cost of evaluating other non-trivial, but known functions $\varphi$ (think
of the value function when all the required regression coefficients are
already known) will be denoted by $\operatorname{cost}(\varphi)$.

If an expression involves several such operations, then only the most
expensive constant is reported. (E.g., evaluating a basis function and
multiplying the value by a scalar constant is considered to incur a cost
$c_{f}$.) We may also use constants $c$ which do not depend on the specifics
of the algorithm. We now go through the individual stages of the algorithm.

  1. (1)

Simulating trajectories at cost $\operatorname{cost}_{1}=c_{X}M(J+1)$.

  2. (2)

Computing the terminal value function as in (4.1) for a given
$x\in\mathbb{R}^{d}$ and all $y\in\mathcal{L}$ at cost
$\operatorname{cost}_{2}=c_{f}\left\lvert\mathcal{L}\right\rvert\left\lvert\mathcal{K}\right\rvert$.

  3. (3)

For fixed $0\leq j\leq J-1$ and $y\in\mathcal{L}$ set up the least squares
problem (4.9) at cost
$M\left(c_{f}K+\operatorname{cost}\left(v^{(I)}_{j+1}\right)\right)$.

  4. (4)

For fixed $0\leq j\leq J-1$ and $y\in\mathcal{L}$, we solve the least squares
problem (4.9) at cost $c_{\ast}MK^{2}$. The total cost is
$\operatorname{cost}_{4}=c_{\ast}JMK^{2}\left\lvert\mathcal{L}\right\rvert$.

  5. (5)

For fixed $0\leq j\leq J-1$, $y\in\mathcal{L}$, and $1\leq i\leq I$ set up the
least squares problem (4.10) at cost
$M\left(c_{f}K+\operatorname{cost}\left(v^{(I)}_{j+1}\right)+R\operatorname{cost}\left(v_{j+1}^{(i-1)}\right)\right)$.

  6. (6)

For fixed $0\leq j\leq J-1$, $y\in\mathcal{L}$, and $1\leq i\leq I$ solve the
least squares problem (4.10) at cost $c_{\ast}M(K+R)^{2}$, leading to a total
cost of
$\operatorname{cost}_{6}=c_{\ast}M(K+R)^{2}J\left\lvert\mathcal{L}\right\rvert$.

List 5.1:  Stages of the algorithm

For simplicity of the presentation, we shall only consider the following
scenario:

######  Assumption 5.1.

The total set of reinforced basis functions contains all available value
functions, i.e., $\bigcup_{y\in\mathcal{L}}\mathcal{L}^{y}=\mathcal{L}.$

For fixed $0\leq i\leq I$ and $0\leq j\leq J$ let

(5.1)
$\mathbf{v}^{(i)}_{j}\coloneqq\left(v^{(i)}_{j}(y,\cdot)\right)_{y\in\mathcal{L}},\quad\mathbf{c}^{(i)}_{j}\coloneqq\left(c^{(i)}_{j}(y,\cdot)\right)_{y\in\mathcal{L}}.$

The key step of the cost analysis is understanding the cost of evaluating the
reinforced basis functions, which are, in turn, given in terms of reinforced
basis functions at later time steps. We note that it is essential to analyze
the cost of evaluating the full set of reinforced basis functions
$\mathbf{v}^{(i)}_{j}$ rather than individual ones $v^{(i)}_{j}(y,\cdot)$, as
the latter method would show us an apparent explosion of basis functions as we
increase time.111Suppose that each reinforced basis function
$v^{(i)}_{j}(y,\cdot)$ depends on two reinforced basis functions
$v^{(i-1)}_{j+1}(y^{\prime},\cdot)$ and
$v^{(i-1)}_{j+1}(y^{\prime\prime},\cdot)$. If we follow this recursion for
$l\leq i$ steps, we arrive at a total set of $2^{l}$ basis functions. The
catch is that many, if not all, of these basis functions overlap with basis
functions for other reinforced basis functions $v^{(i)}_{j}(\tilde{y},\cdot)$.
By (4.7), evaluating $\mathbf{v}^{(i)}_{j}$ requires evaluating the payoff
functions for all combinations of controls $y\in\mathcal{L}$ and policies
$a\in\mathcal{K}$, then evaluating $\mathbf{c}^{(i)}_{j}$, and taking the
corresponding maxima. In total, this means

$\operatorname{cost}\left(\mathbf{v}^{(i)}_{j}\right)\leq\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+\operatorname{cost}\left(\mathbf{c}^{(i)}_{j}\right).$

On the other hand, by (4.6) evaluating $\mathbf{c}^{(i)}_{j}$ requires $K$
evaluations of standard basis functions, $K\left\lvert\mathcal{L}\right\rvert$
elementary operations for summing them, one evaluation of
$\mathbf{v}^{(i-1)}_{j+1}$, and $\left\lvert\mathcal{L}\right\rvert^{2}$
elementary operations for their summation. In total, this means that

$\operatorname{cost}\left(\mathbf{c}^{(i)}_{j}\right)\leq
Kc_{f}+K\left\lvert\mathcal{L}\right\rvert
c_{\ast}+\mathds{1}_{i>0}\left(\left\lvert\mathcal{L}\right\rvert^{2}c_{\ast}+\operatorname{cost}\left(\mathbf{v}^{(i-1)}_{j+1}\right)\right).$

This implies the cost estimate

(5.2)
$\operatorname{cost}\left(\mathbf{v}^{(i)}_{j}\right)\leq\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+K\left\lvert\mathcal{L}\right\rvert
c_{\ast}+\mathds{1}_{i>0}\left(\left\lvert\mathcal{L}\right\rvert^{2}c_{\ast}+\operatorname{cost}\left(\mathbf{v}^{(i-1)}_{j+1}\right)\right).$

######  Lemma 5.2.

The cost of evaluating $\mathbf{v}^{(i)}_{j}$, $i=0,\ldots,I$, $j=0,\ldots,J$
can be bounded by

$\operatorname{cost}\left(\mathbf{v}^{(i)}_{j}\right)\leq\begin{cases}(i+1)\left(\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+K\left\lvert\mathcal{L}\right\rvert
c_{\ast}\right)+i\left\lvert\mathcal{L}\right\rvert^{2}c_{\ast},&j+i\leq J,\\\
\left\lvert\mathcal{L}\right\rvert\left\lvert\mathcal{K}\right\rvert(c_{f}+c_{\ast})+(J-j)\left(\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+(K+1)\left\lvert\mathcal{L}\right\rvert
c_{\ast}\right),&j+i>J.\end{cases}$

###### Proof.

For
$a\coloneqq\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+K\left\lvert\mathcal{L}\right\rvert
c_{\ast}$, consider the cost recursion

$c(k+1)\leq a+\left\lvert\mathcal{L}\right\rvert c_{\ast}c(k),\quad k\geq 0.$

Assuming that the recursion hits $i=0$ before $j=J$, i.e., $i+j\leq J$, the
cost $c(k)\coloneqq\operatorname{cost}\left(\mathbf{v}_{j+i-k}^{(k)}\right)$
satisfies the recursion with $c(0)\leq a$, and, hence, we obtain

$c(k)\leq(k+1)a+k\left\lvert\mathcal{L}\right\rvert^{2}c_{\ast}.$

This gives the first expression in the statement of the lemma with $k=i$.

On the other hand, if $i+j>J$, we hit $j=J$ before $i=0$. In this case,
$c(k)\coloneqq\operatorname{cost}\left(\mathbf{v}^{(i+j-J+k)}_{J-k}\right)$
satisfies the same recursion, but with initial value
$c(0)\leq\left\lvert\mathcal{L}\right\rvert\left\lvert\mathcal{K}\right\rvert(c_{f}+c_{\ast})$.
∎

In order to shorten notation, we introduce

$\displaystyle a$ $\displaystyle\coloneqq$
$\displaystyle\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+K\left\lvert\mathcal{L}\right\rvert
c_{\ast},$ $\displaystyle b$ $\displaystyle\coloneqq$
$\displaystyle\left\lvert\mathcal{L}\right\rvert^{2}c_{\ast},$ $\displaystyle
d$ $\displaystyle\coloneqq$
$\displaystyle\left\lvert\mathcal{L}\right\rvert\left\lvert\mathcal{K}\right\rvert(c_{f}+c_{\ast}),$
$\displaystyle e$ $\displaystyle\coloneqq$
$\displaystyle\left\lvert\mathcal{K}\right\rvert\left\lvert\mathcal{L}\right\rvert(c_{f}+c_{\ast})+Kc_{f}+(K+1)\left\lvert\mathcal{L}\right\rvert
c_{\ast},$

so that the estimate of Lemma 5.2 shortens to

$\operatorname{cost}\left(\mathbf{v}^{(i)}_{j}\right)\leq\begin{cases}(i+1)a+ib,&j+i\leq
J,\\\ d+(J-j)e,&j+i>J.\end{cases}$

We next estimate the cost of setting up the regression problem (4.9), which is
proved similarly.

######  Lemma 5.3.

The cost of setting up the regression problem for $c^{(0)}_{j}(y,\cdot)$,
$j=0,\ldots,J-1$, $y\in\mathcal{L}$, can be bounded by

$\operatorname{cost}_{3}\leq
JMKc_{f}+M(J-I)\left((I+1)a+Ib\right)+MId+\frac{1}{2}MI(I+1)e.$

The cost for setting up the least squares problem (4.10) is computed in a
similar way.

######  Lemma 5.4.

The cost of setting up the regression problem for $c^{(i)}_{j}(y,\cdot)$,
$i=1,\ldots,I$, $j=0,\ldots,J-1$, $y\in\mathcal{L}$, can be bounded by

$\operatorname{cost}_{5}\leq
JMKc_{f}+\frac{M}{2}I\left[(I+1)a+(I-1)b\right](J-I+2)+\\\
+\frac{M}{6}I\left[11a+2b-9d+5e+I(I+6)a+3I(I+b)+3I^{2}d+I(I+6)e\right].$

###### Proof.

A closer look at (4.10) reveals that the total cost of setting up all these
least squares problems can be bounded by

(5.3)
$\operatorname{cost}_{5}\leq\sum_{j=0}^{J-1}M\left(Kc_{f}+\sum_{i=1}^{I}\operatorname{cost}\left(\mathbf{v}^{(i-1)}_{j+1}\right)\right),$

taking into account that $\mathbf{v}^{(I)}_{j+1}$ was already evaluated during
the set-up of the least squares problem (4.9) and, hence, does not need to be
evaluated again. Using Lemma 5.2, we obtain

$\displaystyle\operatorname{cost}_{5}$ $\displaystyle\leq
JMKc_{f}+M\sum_{j=0}^{J-1}\left[\sum_{i=0}^{(J-j)\wedge(I-1)}\left((i+1)a+ib\right)+\sum_{i=1+(J-j)\wedge(I-1)}^{I-1}\left(d+(J-j-1)e\right)\right]$
$\displaystyle=JMKc_{f}+M\sum_{j=0}^{J-I+1}\left[\sum_{i=0}^{I-1}\left((i+1)a+ib\right)\right]+$
$\displaystyle\quad\quad+M\sum_{j=J-I}^{J-1}\left[\sum_{i=0}^{J-j}\left((i+1)a+ib\right)+\sum_{i=J-j+1}^{I-1}\left(d+(J-j-1)e\right)\right].$

Evaluating the double sums gives the estimate from the statement of the lemma.
∎

Abandoning the difference between $c_{f}$ and $c_{\ast}$ using the trivial
bounds
$\operatorname{cost}_{3}\leq\operatorname{const}\operatorname{cost}_{5}$,
$\operatorname{cost}_{4}\leq\operatorname{const}\operatorname{cost}_{6}$, we
obtain

######  Theorem 5.5.

The computational cost of the algorithm presented in Section 4.2 can be
bounded by

$\operatorname{cost}\leq\operatorname{const}MJ\left(c_{X}+I^{2}(K+\left\lvert\mathcal{K}\right\rvert+\left\lvert\mathcal{L}\right\rvert)\left\lvert\mathcal{L}\right\rvert+(K+R)^{2}\left\lvert\mathcal{L}\right\rvert\right),$

where $\operatorname{const}$ is a positive number independent of
$\left\lvert\mathcal{K}\right\rvert$, $\left\lvert\mathcal{L}\right\rvert$,
$K$, $J$, and $I$.

######  Remark 5.6.

Recall that Remark 4.4 introduced a significantly cheaper variant of algorithm
B for the case $I=J$. It is easy to see that the computational cost of this
variant is bounded by

$\operatorname{cost}\leq\operatorname{const}MJ\left(c_{X}+J(K+\left\lvert\mathcal{K}\right\rvert+\left\lvert\mathcal{L}\right\rvert)\left\lvert\mathcal{L}\right\rvert+(K+R)^{2}\left\lvert\mathcal{L}\right\rvert\right),$

i.e., the total cost is proportional to $J^{2}$ rather than $J^{3}$. Indeed,
the main difference in the cost analysis as compared to the full modified
algorithm is that (5.3) can be replaced by

$\operatorname{cost}_{5}\leq\sum_{j=0}^{J-1}M\left(Kc_{f}+\operatorname{cost}\left(\mathbf{v}^{(J-j-1)}_{j+1}\right)\right).$

Note that this essentially corresponds to the algorithm of [BSSZ20] directly
generalized to optimal control problems.

##  6\. Convergence analysis

In this section we analyze the convergence properties of the standard and
reinforced regression algorithms introduced in the previous sections. For
related convergence analysis in the case of optimal stopping problems we refer
the interested reader to [Zan13], [Zan18], and [BS20], see also [BRS21].
Henceforth we assume that

(6.1)
$\max_{j=0,\ldots,J}\sup_{y\in\mathcal{L}}\sup_{a\in\mathcal{K}}\sup_{x\in\mathcal{X}}\left|H_{j}(a,y,x)\right|\leq
C_{H},$

then all the value functions

$v_{j}^{\ast}(y,x):=\sup_{\mathbf{A}=(A_{\ell})_{\ell=j}^{J}\in\mathcal{A}_{j}(y,x)}\mathsf{E}\left[\sum_{\ell=j}^{J}H_{\ell}\left(A_{\ell},Y_{\ell}(\mathbf{A};j,y),X_{\ell}^{j,x}\right)\right]$

are uniformly bounded by $JC_{H}.$ Fix a sequence of spaces $\Psi_{j},$
$j=0,\ldots,J,$ of functions defined on $\mathcal{X}.$ We stress that these
spaces are not necessarily linear at this point. Construct the corresponding
sequence of estimates $(v_{j,M}(y,x))_{j=0}^{J}$ via

(6.2) $\displaystyle v_{J,M}(y,x)$ $\displaystyle=\sup_{a\in
K_{j}(y,x)}H_{J}(a,y,x)\text{ \ \ }\mathrm{and}$ $\displaystyle v_{j,M}(y,x)$
$\displaystyle=\sup_{a\in
K_{j}(y,x)}\left(H_{j}(a,y,x)+T_{W}\mathcal{P}_{j,M}[v_{j+1,M}](\varphi_{j+1}(a,y),x)\right),\text{
\ \ }j<J,$

where $\mathcal{P}_{j,M}[g](z,x)$ stands for the empirical projection of the
conditional expectation $\mathsf{E}[g(z,X_{j+1}^{j,x})]$ on $\Psi_{j},$ based
on a sample

(6.3) $\mathcal{D}_{M,j}=\Bigl{\\{}(X_{j}^{(m)},X_{j+1}^{(m)}),\text{ \ \
}m=1,\ldots,M\Bigr{\\}}$

from the joint distribution of $(X_{j},X_{j+1}),$ that is,

${\mathcal{P}}_{j,M}[g](z,\cdot)\in\arg\inf_{\psi\in\Psi_{j}}\sum_{m=1}^{M}\left[\left|g(z,X_{j+1}^{(m)})-\psi(X_{j}^{(m)})\right|^{2}\right].$

In (6.2) $T_{W}$ is a truncation operator at level $W=JC_{H}$ defined by

$T_{W}f(x)=\begin{cases}f(x),&|f(x)|\leq W,\\\
W\mathrm{sign}(f(x)),&\mathrm{otherwise.}\end{cases}$

Due to Theorem 11.5 in [GKKW02], one has for all $g$ with
$\left\|g\right\|_{\infty}\leq W,$ $j=0,\ldots,J-1,$ and all
$z\in\mathcal{L},$ that

(6.4)
$\mathsf{E}\left[\left\|T_{W}\mathcal{P}_{j,M}[g](z,\cdot)-\mathsf{E}\left[g(z,X_{j+1}^{j,\cdot})\right]\right\|_{L_{2}(\mu_{j})}^{2}\right]\\\
\leq\varepsilon_{j,M}^{2}+2\underset{w\,\in\,\Psi_{j}}{\inf}\left\|g(z,\cdot)-w\right\|_{L_{2}(\mu_{j})}^{2}\text{
\ \ with \ \ }\varepsilon_{j,M}^{2}:=cW^{4}\frac{1+\log
M}{M}\mathtt{VC}(\Psi_{j}),$

where $\mathtt{VC}(\Psi_{j})$ is the Vapnik-Chervonenkis dimension of
$\Psi_{j}$ (see Definition 9.6 in [GKKW02]$),\mu_{j}$ is the distribution of
$X_{j},$ and $c$ is an absolute constant. In order to keep the analysis
tractable, we assume that the sets $\mathcal{D}_{M,j}$ are independent for
different $j$, see Remark 6.5 below. More specifically, we consider an
algorithmic framework based on (6.2), where for every exercise date the
samples (6.3) are simulated independently, and consider the information sets

$\mathcal{G}_{j,M}:=\sigma\left\\{\mathbf{X}^{j;M},\ldots,\mathbf{X}^{J;M}\right\\}\text{
with }\mathbf{X}^{j;M}:=\bigl{(}X_{j}^{(m)},\,m=\ 1,\ldots,M\bigr{)}.$

Let us define for $j<J,$ $z\in\mathcal{L},$ $x\in\mathcal{X},$

(6.5) $\widehat{C}_{j}(z,x):=T_{W}\mathcal{P}_{j,M}[v_{j+1,M}](z,x),$

and for a generic (exact) dummy trajectory $\left(X_{l}\right)_{l=0,\ldots,J}$
independent of $\mathcal{G}_{j,M},$ let

(6.6)
$\widetilde{C}_{j}(z,x):=\mathsf{E}_{\mathcal{G}_{j+1,M}}\left[v_{j+1,M}(z,X^{j,x}_{j+1})\right].$

Note that $\widetilde{C}_{j}\left(\cdot,\cdot\right)$ is a
$\mathcal{G}_{j+1,M}$-measurable random function while the estimate
$\widehat{C}_{j}\left(\cdot,\cdot\right)$ is a $\mathcal{G}_{j}$-measurable
one. We further define

(6.7)
$C_{j}^{\ast}(z,x)=\mathsf{E}\left[v_{j+1}^{\ast}(z,X^{j,x}_{j+1})\right],\quad
j<J.$

The following lemma holds.

######  Lemma 6.1.

We have that,

(6.8)
$\mathsf{E}\left[\Bigl{\|}\sup_{z\in\mathcal{L}}\Bigl{|}\widetilde{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\Bigr{|}\Bigr{\|}_{L_{2}(\mu_{j})}^{2}\right]\leq\mathsf{E}\left[\Bigl{\|}\sup_{z\in\mathcal{L}}\Bigl{|}\widehat{C}_{j+1}(z,\cdot)-C_{j+1}^{\ast}(z,\cdot)\Bigr{|}\Bigr{\|}_{L_{2}(\mu_{j+1})}^{2}\right].$

###### Proof.

Let $X$ be a generic (exact) dummy trajectory independent of
$\mathcal{G}_{j+1,M}.$ Then from (6.6), and (6.7) we see that for $j<J,$

(6.9)
$\left|\widetilde{C}_{j}(z,X_{j})-C_{j}^{\ast}(z,X_{j})\right|\leq\mathsf{E}_{\mathcal{G}_{j+1,M}}\left[\left.\left|v_{j+1,M}(z,X_{j+1})-v_{j+1}^{\ast}(z,X_{j+1})\right|\right|X_{j}\right]$

Next, by (2.7), (6.2), (6.5), and (6.7) we have that

(6.10) $\displaystyle\left|v_{j+1,M}(z,x)-v_{j+1}^{\ast}(z,x)\right|$
$\displaystyle\leq$ $\displaystyle\sup_{a\in
K_{j+1}(z,x)}\left|\widehat{C}_{j+1}(\varphi_{j+2}(a,z),x)-C_{j+1}^{\ast}(\varphi_{j+2}(a,z),x)\right|$
$\displaystyle\leq$
$\displaystyle\sup_{z^{\prime}\in\mathcal{L}}\left|\widehat{C}_{j+1}(z^{\prime},x)-C_{j+1}^{\ast}(z^{\prime},x)\right|.$

Hence, by (6.9) one has that

$\sup_{z\in\mathcal{L}}\left|\widetilde{C}_{j}(z,X_{j})-C_{j}^{\ast}(z,X_{j})\right|\leq\mathsf{E}_{\mathcal{G}_{j+1,M}}\left[\left.\sup_{z\in\mathcal{L}}\left|\widehat{C}_{j+1}(z,X_{j+1})-C_{j+1}^{\ast}(z,X_{j+1})\right|\right|X_{j}\right].$

Finally, by taking the “all-in expectation” w.r.t. the law
$\mathcal{\mu}_{j}\otimes\mathbb{P}_{M}$, we observe that

$\displaystyle\mathsf{E}\left[\sup_{z\in\mathcal{L}}\left|\widetilde{C}_{j}(z,X_{j})-C_{j}^{\ast}(z,X_{j})\right|^{2}\right]$
$\displaystyle\leq$
$\displaystyle\mathsf{E}\left\\{\mathsf{E}_{\mathcal{G}_{j+1,M}}\left[\left.\sup_{z\in\mathcal{L}}\left|\widehat{C}_{j+1}(z,X_{j+1})-C_{j+1}^{\ast}(z,X_{j+1})\right|\right|X_{j}\right]\right\\}^{2}$
$\displaystyle\leq$
$\displaystyle\mathsf{E}\left[\sup_{z\in\mathcal{L}}\left|\widehat{C}_{j+1}(z,X_{j+1})-C_{j+1}^{\ast}(z,X_{j+1})\right|^{2}\right]$

by Jensen’s inequality and the tower property. ∎

In fact, Lemma 6.1 is the key to the next proposition.

######  Proposition 6.2.

Set

$\mathcal{E}_{j}:=\Bigl{\|}\sup_{z\in\mathcal{L}}\left|\widehat{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\right|\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})},\quad
j=0,\ldots,J-1,$

with $\mathbb{P}_{M}$ being the law of the sample $X_{j}^{(m)},\,m=\
1,\ldots,M,\,j=1,\ldots,J.$ Then it holds

(6.11)
$\mathcal{E}_{j}\leq\left|\mathcal{L}\right|\Bigl{(}\varepsilon_{j,M}+\sqrt{2}\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\bigl{\|}\widetilde{C}_{j}(z,\cdot)-w\bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\Bigr{)}+\left|\mathcal{L}\right|\mathcal{E}_{j+1}$

for all $j=0,\ldots,J-1,$ with $\mathcal{E}_{J}=0$ by definition.

###### Proof.

The case $j=J-1$ follows from (6.4) and the fact that
$\widetilde{C}_{J-1}=C^{\ast}_{J-1}.$ Set
$r_{j,M}(z)=\underset{w\,\in\,\Psi_{j}}{\inf}\|\widetilde{C}_{j}(z,\cdot)-w\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}.$
Due to (6.4) we have with probability $1,$

(6.12)
$\mathsf{E}_{\mathcal{G}_{j+1,M}}\left[\left\|\widehat{C}_{j}(z,\cdot)-\widetilde{C}_{j}(z,\cdot)\right\|_{L_{2}(\mu_{j})}^{2}\right]\leq\varepsilon_{j,M}^{2}+2r^{2}_{j,M}(z).$

Hence

(6.13)
$\Bigl{\|}\widehat{C}_{j}(z,\cdot)-\widetilde{C}_{j}(z,\cdot)\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\varepsilon_{j,M}+\sqrt{2}r_{j,M}(z).$

By applying (6.13) it follows that

$\displaystyle\Bigl{\|}\widehat{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}$
$\displaystyle\leq$
$\displaystyle\varepsilon_{j,M}+\sqrt{2}r_{j,M}(z)+\left\|\widetilde{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}.$

From this and Lemma 6.1 we imply

$\displaystyle\sup_{z\in\mathcal{L}}\Bigl{\|}\widehat{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}$
$\displaystyle\leq$
$\displaystyle\varepsilon_{j,M}+\sqrt{2}\sup_{z\in\mathcal{L}}r_{j,M}(z)$
$\displaystyle+\Bigl{\|}\sup_{z\in\mathcal{L}}\left|\widehat{C}_{j+1}(z,\cdot)-C_{j+1}^{\ast}(z,\cdot)\right|\Bigr{\|}_{L_{2}(\mu_{j+1}\otimes\mathbb{P}_{M})}$

and then (6.11) follows. ∎

######  Corollary 6.3.

Suppose that

$\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\bigl{\|}\widetilde{C}_{j}(z,\cdot)-w\bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\delta,\quad\mathtt{VC}(\Psi_{j})\leq
D,\quad 0\leq j\leq J-1,$

for some $\delta>0$ and $D>0.$ Proposition 6.2 then yields for
$j=0,\ldots,J-1,$ by using (6.10),

(6.14)
$\Bigl{\|}\sup_{z\in\mathcal{L}}\left|v_{j,M}(z,\cdot)-v_{j}^{\ast}(z,\cdot)\right|\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\left(cW^{4}\frac{1+\log
M}{M}D+\sqrt{2}\delta\right)\frac{\left|\mathcal{L}\right|^{J-j+1}-\left|\mathcal{L}\right|}{\left|\mathcal{L}\right|-1}.$

######  Corollary 6.4.

By inserting the estimate

$\underset{w\,\in\,\,\Psi_{j}}{\inf}\left\|\widetilde{C}_{j}(z,\cdot)-w\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\left\|\widetilde{C}_{j}(z,\cdot)-C_{j}^{\ast}(z,\cdot)\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\\!\\!\\!+\underset{w\,\in\,\,\Psi_{j}}{\inf}\left\|C_{j}^{\ast}(z,\cdot)-w\right\|_{L_{2}(\mu_{j})}$

in Proposition 6.2, we get the alternative recursion

$\mathcal{E}_{j}\leq\left|\mathcal{L}\right|\left(\varepsilon_{K,M}+\sqrt{2}\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\left\|C_{j}^{\ast}(z,\cdot)-w\right\|_{L_{2}(\mu_{j})}\right)+\left|\mathcal{L}\right|(1+\sqrt{2})\mathcal{E}_{j+1},$

and under the alternative assumption

$\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\bigl{\|}{C}^{\ast}_{j}(z,\cdot)-w\bigr{\|}_{L_{2}(\mu_{j})}\leq\delta,\quad\mathtt{VC}(\Psi_{j})\leq
D,\quad 0\leq j\leq J-1,$

for some $\delta>0$ and $D>0,$ we obtain for $j=0,\ldots,J,$ the bounds

$\Bigl{\|}\sup_{z\in\mathcal{L}}\left|v_{j,M}(z,\cdot)-v_{j}^{\ast}(z,\cdot)\right|\Bigr{\|}_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\left(cW^{4}\frac{1+\log
M}{M}D+\sqrt{2}\delta\right)\left|\mathcal{L}\right|\frac{\left((1+\sqrt{2})\left|\mathcal{L}\right|\right)^{J-j}-1}{(1+\sqrt{2})\left|\mathcal{L}\right|-1}.$

######  Remark 6.5.

In [Zan13] the convergence of an “independent sample version” of the
Longstaff-Schwartz algorithm was studied based on an assumption similar to the
independence $\mathcal{D}_{M,j}$ for different $j$ here. However, in a later
paper [Zan18] it was shown (by more involved analysis) that the convergence
rates based on one and the same sample are basically the same as in [Zan13] up
to certain constants. One therefore may naturally expect that similar
conclusions apply in our context. Therefore the numerical examples in Section
7 are based on a single sample of $M$ trajectories.

The proposed reinforced regression algorithm with $I=J$ uses linear
approximation spaces of the form

(6.15)
$\displaystyle\Psi_{j}=\mathrm{span}\\{\psi_{1}(x),\ldots,\psi_{K}(x),{v}_{j+1,M}(y_{1},x),\ldots,{v}_{j+1,M}(y_{R},x)\\},\quad
j=0,\ldots,J-1,$

for $\mathcal{L}=\\{y_{1},\ldots,y_{R}\\},$ where
$\psi_{1}(x),\ldots,\psi_{K}(x)$ are some fixed basis functions (e.g.
polynomials) on $\mathcal{X}.$ In this case $\mathtt{VC}(\Psi_{j})\leq
K+R,\quad 0\leq j\leq J-1.$ In order to see the advantage of adding additional
basis functions more clearly, we prove the following proposition.

######  Proposition 6.6.

Assume additionally that

(6.16)
$\max_{j=1,\ldots,J}\sup_{y\in\mathcal{L}}\sup_{a\in\mathcal{K}}\sup_{x\in\mathcal{X}}\left|H_{j}(a,y,x_{1})-H_{j}(a,y,x_{2})\right|\leq
L_{H}|x_{1}-x_{2}|,$

for all $x_{1},x_{2}\in\mathcal{X}$ and

(6.17)
$\displaystyle\max_{j=1,\ldots,J}\max_{\ell=j,\ldots,J}\mathsf{E}[|X_{\ell}^{j,x_{1}}-X_{\ell}^{j,x_{2}}|]\leq
L_{X}|x_{1}-x_{2}|,\quad\forall x_{1},x_{2}\in\mathcal{X}$

for some constants $L_{H}>0,$ $L_{X}>0.$ Then it holds for reinforced spaces
$(\Psi_{j})$ from (6.15)

$\displaystyle\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\left\|\widetilde{C}_{j}(z,\cdot)-w\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq
JL_{X}L_{H}\left[\mathsf{E}\int|X^{j,x}_{j+1}-x|^{2}\mu_{j}(dx)\right]^{1/2}.$

###### Proof.

We have

$\displaystyle\sup_{z\in\mathcal{L}}\underset{w\,\in\,\,\Psi_{j}}{\inf}\left\|\widetilde{C}_{j}(z,\cdot)-w\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}\leq\sup_{z\in\mathcal{L}}\left\|\mathsf{E}_{\mathcal{G}_{j+1,M}}\bigl{[}v_{j+1,M}(z,X^{j,\cdot}_{j+1})-v_{j+1,M}(z,\cdot)\bigr{]}\right\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}.$

Under assumptions (6.16) and (6.17), we then have

$\displaystyle\max_{j=1,\ldots,J}\sup_{y\in\mathcal{L}}|v^{\ast}_{j}(y,x_{1})-v^{\ast}_{j}(y,x_{2})|\leq
JL_{X}L_{H}|x_{1}-x_{2}|,\quad\forall x_{1},x_{2}\in\mathcal{X}.$

By using an additional truncation, one can achieve that the Lipschitz
constants of the estimates $v_{j,M}(z,\cdot),$ $j=0,\ldots,J-1,$ are all
uniformly bounded by a constant $JL_{X}L_{H}$ with probability $1.$ ∎

The above proposition implies that if $JL_{H}$ stays bounded for $J\to\infty,$
(for example if $H$ scales as $1/J$), then the approximation error
$\underset{w\,\in\,\,\Psi_{j}}{\inf}\|\widetilde{C}_{j}(z,\cdot)-w\|_{L_{2}(\mu_{j}\otimes\mathbb{P}_{M})}$
becomes small as $J\to\infty.$ Note that the latter property can not be
guaranteed when using fixed (nonadaptive) linear spaces $\Psi_{j}.$ Of course,
the exponential in $J$ factor in (6.14) will lead to explosion of the overall
error as $J\to\infty,$ but the above observation still indicates that the
inclusion of the functions ${v}_{j+1,M}(y_{1},x),\ldots,{v}_{j+1,M}(y_{R},x)$
into $\Psi_{j}$ can significantly improve the quality of the estimates
$v_{j,M}(z,\cdot)$ especially in the case of large $J$. Concerning the
dependence of the bound (6.14) on $J$, we note that this estimate is likely to
be too pessimistic, see also a discussion in [Zan13].

##  7\. Numerical examples

We now present various numerical examples which demonstrate the accuracy of
the reinforced regression algorithm in practice. To allow for a direct
comparison with the reinforced regression algorithm of [BSSZ20], we first
consider a (single) optimal stopping problem, more particularly a Bermudan
max-call option. Our second example is a multiple stopping problem, for which
the hyperparameters already become crucial. Finally, our last example is an
optimal control of a gas storage.

We have tested both algorithms A and B for the HRR (hierarchical reinforced
regression) method, however, the latter version always gave slightly better
results and therefore we have only included the value obtained with the
algorithm B. Intuitively this is to be expected, as the algorithm B uses more
accurate regression targets in the backwards induction. The algorithm A may be
of use in situations where one is interested in improving the approximation
until a certain accuracy threshold is reached, however, properly done this
approach should also include a calculation of upper bounds, which we do not
discuss in our paper.

Before, let us also mention how a lower biased estimate to the value of a
control problem in a Markovian setting is calculated using the result of a
regression procedure. Let $c$ be an approximation to the function $c^{\ast}$
given by

$c^{*}_{j}(x,y)=\mathsf{E}\left[v_{j+1}^{\ast}(y,X_{j+1}^{j,x})\right],\quad
x\in\mathbb{R}^{d},\;y\in\mathcal{L},\;j\in\\{0,...,J-1\\},$

with $c_{J}\equiv c^{*}_{J}\equiv 0$ by convention. Using the hierarchical
reinforced regression method, such an approximation is given by $c^{(I)}$
defined in (4.6). Further let $(X^{(m)})_{1\leq m\leq M_{\mathrm{test}}}$ be
sample trajectories from the underlying Markov chain, generated independently
from the samples used in the regression procedure. Then we can iteratively
define a sequence of polices $(\mathbf{A}^{(m)})_{1\leq m\leq
M_{\mathrm{test}}}$ with $\mathbf{A}^{(m)}=(A^{(m)}_{0},...,A^{(m)}_{J})$ by

$A^{(m)}_{j}:=\operatornamewithlimits{arg\,max}_{A\in
K_{j}(Y^{(m)}_{j},X_{j}^{(m)})}\left(H_{j}(A,Y_{j}^{(m)},X^{(m)}_{j})+c_{j}(\varphi(A,Y^{(m)}_{j}),X^{(m)}_{j})\right),$

or all $m=1,...,M_{\mathrm{test}}$ and $j=0,...,J$, where
$Y^{(m)}_{0}:=y_{0}\in\mathcal{L}$ and
$Y^{(m)}_{j+1}:=\varphi_{j+1}(A_{j}^{m},Y^{(m)}_{j})$. It then follows from
the definition, that each $\mathbf{A}^{(m)}$ is an admissible sequence of
policies, i.e. $\mathbf{A}^{(m)}\in\mathcal{A}_{0}(y_{0},X^{(m)})$. Therefore,
a lower estimate to the value $\mathsf{E}(v(y_{0},X_{0}))$ is given by

$\frac{1}{M_{\mathrm{test}}}\sum_{m=1}^{M_{\mathrm{test}}}\sum_{j=0}^{J}H_{j}(A^{(m)}_{j},Y_{j}^{(m)},X_{j}^{(m)}).$

Lower bounds allow a direct comparison of the performance of different
methods, in the sense that the method yielding the highest lower bound (up to
Monte Carlo errors) performed best, since this value must be closest to the
true value of the control problem. This direct comparison is, however, not
possible for the approximate value $v_{0}$ which may lie above or below the
true value. Our main premise in the following is that the HRR algorithm is a
more efficient way to improve the performance of the regression algorithm as
compared with increasing the complexity of the regression basis. Hence, for
this relative comparison it is also sufficient to study lower bounds.

###  7.1. Bermudan max-call option

In this section we evaluate the performance of the hierarchical reinforced
regression (HRR) method from Section 4 on the valuation of a Bermudan max-call
option. Let $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{0\leq t\leq T},P)$ be a
filtered probability space on which a $d$-dimensional Brownian motion
$W=(W(t))_{0\leq t\leq T}$ is defined. Further let $X=(X(t))_{0\leq t\leq T}$
be the geometric Brownian motion defined by

$dX^{k}(t)=(r-\delta)\,X^{k}(t)\,dt+X^{k}(t)\,\sigma\,dW^{k}(t),\quad
X^{k}(0)=x_{0},\quad 0\leq t\leq T,\quad k\in\\{1,...,d\\},$

where $x_{0},r,\delta,\sigma>0$. Option rights can be exercised on a
predefined set of possible exercise dates $\\{t_{0},t_{1},...,t_{J}\\}$, where
at most one right can be exercised on any given date. Assume that the exercise
dates are equidistant $t_{j}:=j\cdot\Delta t$ for all $j\in\\{0,...,J\\}$ with
$\Delta t:=T/J$ and define the underlying Markov chain $(X_{j})_{j=0,...,J}$
by $X_{j}=X({j\Delta t})$. Recall from Example 2.3 that in order to model a
multiple stopping problem in the optimal control framework we define the set
of policies by $\mathcal{K}=\\{0,1\\}$ and the set of controls by
$\mathcal{L}=\\{0,...,y_{max}\\}$, where $y_{max}$ is the number of exercise
rights. Further, we define

$\varphi_{j}(a,y)\coloneqq(y-a)_{+},\quad K_{j}(x,y)\coloneqq\\{0,1\wedge
y\\},\quad H_{j}(a,y,x)\coloneqq a\cdot g_{j}(x)\,e^{-t_{j}r},$

for all $y\in\mathcal{L}$, $a\in\mathcal{K}$ and $j\in\\{0,...,J\\}$, where
$g$ is the max-call pay-off function defined by

$g(x)\coloneqq(\max\\{x^{1},....,x^{d}\\}-C)_{+},\quad x\in\mathbb{R}^{d},$

where $C\in\mathbb{R}_{+}$ is the option strike. Then the value function
$v^{*}_{0}(y_{max},x_{0})$ defined in (2.6) yields the value of the Bermudan
max-call option with underlying $X$ and data
$(d,J,T,y_{max},x_{0},C,r,\delta,\sigma)$.

####  7.1.1. Single exercise right

We will first consider the case of a single exercise right $y_{max}=1$. This
is a standard example in the literature, see for example [BG04, Rog02, AB04]
and more recently [BCJ19]. The performance of the reinforced regression method
for this example was already analyzed in [BSSZ20]. We revisit this example in
order to demonstrate that even in the optimal stopping case our novel HRR
method allows for improvements in computational costs without sacrificing the
quality of the estimations.

$d$ | Basis | Lower bounds | CI from [BCJ19]
---|---|---|---
Regression | HRR | Reinf. Reg.
$I=0$ | $I=1$ | $I=9$
2 | $\Psi_{1}$ | 13.015 (0.022) | 13.772 (0.015) | 13.794 (0.015) | [13.880,13.910]
$\Psi_{1,g}$ | 13.679 (0.019) | - | -
$\Psi_{2}$ | 13.775 (0.016) | 13.871 (0.015) | 13.882 (0.014)
$\Psi_{3}$ | 13.874 (0.016) | - | -
3 | $\Psi_{1}$ | 17.764 (0.029) | 18.526 (0.017) | 18.540 (0.018) | [18.673,18.699]
$\Psi_{1,g}$ | 18.404 (0.022) | - | -
$\Psi_{2}$ | 18.519 (0.020) | 18.639 (0.017) | 18.653 (0.017)
$\Psi_{3}$ | 18.655 (0.021) | - | -
5 | $\Psi_{1}$ | 25.463 (0.024) | 25.998 (0.021) | 25.990 (0.019) | [26.138, 26.174]
$\Psi_{1,g}$ | 25.823 (0.026) | - | -
$\Psi_{2}$ | 25.990 (0.023) | 26.097 (0.019) | 26.109 (0.020)
| $\Psi_{3}$ | 26.111 (0.022) | - | - |
10 | $\Psi_{1}$ | 38.022 (0.025) | 38.234 (0.022) | 38.225 (0.024) | [38.300,38.367]
$\Psi_{1,g}$ | 38.058 (0.024) | - | -
$\Psi_{2}$ | 38.299 (0.023) | 38.316 (0.020) | 38.331 (0.021)
| $\Psi_{3}$ | 38.349 (0.021) | - | - |
Table 1. Lower bounds ($\pm$ 99.7% Monte-Carlo error) for the value of the
Bermudan max-call option with data $J=9$, $T=1$, $y_{max}=1$, $x_{0}=C=100$,
$r=0.05$, $\delta=0.1$, $\sigma=0.2$ and different numbers of underlying
assets $d\in\\{2,3,5,10\\}$. For all methods we used $M=10^{6}$ training
sample paths and $M_{\mathrm{test}}=10^{7}$ paths for calculating the lower
bound. The last column presents the 95% confidence intervals for the value of
the Bermuda option from [BCJ19].

Define the functions $f_{i}:\mathbb{R}^{d}\to\mathbb{R}$,
$x\mapsto\operatorname{sort}(x^{1},\ldots,x^{d})^{i}$, the $i$th largest entry
of $x$, for $i\in\\{1,\dots,d\\}$ and consider the following three sets of
regression basis functions:

$\begin{array}[]{c}\Psi_{1}\coloneqq\\{1,f_{1},...,f_{d}\\},\quad\Psi_{1,g}\coloneqq\Psi_{1}\cup\\{g\\},\quad\Psi_{2}\coloneqq\Psi_{1}\cup\\{f_{i}\cdot
f_{j}\;|\;1\leq i\leq j\leq d\\}\\\
\Big{.}\Psi_{3}\coloneqq\Psi_{1}\cup\Psi_{2}\cup\\{f_{i}\cdot f_{j}\cdot
f_{k}\;|\;1\leq i\leq j\leq k\leq d\\}.\end{array}$

Note that the cardinalities of these sets are given by $|\Psi_{1}|=d+1$,
$|\Psi_{1,p}|=d+2$, $|\Psi_{2}|=\frac{1}{2}d^{2}+\frac{3}{2}d+1$, and
$|\Psi_{3}|=\frac{1}{6}d^{3}+d^{2}+\frac{11}{6}d+1$, respectively. Regarding
the HRR method, we use the algorithm of the second type described in Section
4.2. Further note that in the optimal stopping case there is only one choice
for the set of reinforced value function for the HRR method since
$\mathcal{L}=\\{1\\}$ and therefore we always set $\mathcal{L}^{1}=\\{1\\}$.

Figure 4. A visualization of the lower bounds from Table 1.

We considered two different set-ups for the comparison of the different
methods:

  * •

First we keep the number of exercises dates $J$ fixed and vary the number of
underlying assets $d$;

  * •

Second we keep $d$ fixed and vary $J$ (while also keeping $T$ fixed).

In Table 1 we present lower estimates to the value of a Bermudan max-call
option with a single exercise right for $J=9$ and $d\in\\{2,3,5,10\\}$. In the
corresponding Figure 4 we have visualized the lower bounds for the comparison
between the different regression methods. For each of the considered methods
we used $M=10^{6}$ simulated training samples paths to determine the
regression coefficients and $M_{\text{test}}=10^{7}$ paths for calculating the
lower bounds.

In order to give the reader an easy reference point, in Table 1 we have also
included the confidence intervals for the value of the Bermudan max-call
option from [BCJ19]. Note however that we do not aim for an improvement of the
latter values in terms of benchmarking. In fact the method used in [BCJ19] is
quite different from ours, as it uses deep neural networks for approximating
the optimal stopping policies at each time step, and a direct comparison would
require the usage of higher order polynomials for our method.

We first observe that across all numbers of assets $d$ the HRR method with the
set of basis functions $\Psi_{1}$ performs significantly better then the
standard regression method with the basis $\Psi_{1}$ and $\Psi_{1,g}$. The
same holds true when comparing the methods using the regression basis
$\Psi_{2}$. More importantly however, we observe that for $d\leq 5$ the HRR
method with basis $\Psi_{1}$ yields lower bounds of the same quality as
obtained with the standard method and the larger basis $\Psi_{2}$. The same
holds true when comparing the HRR method with basis $\Psi_{2}$ against the
standard method with the basis $\Psi_{3}$. In the case $d=10$ assets, the
lower bounds obtained with the HRR method and basis $\Psi_{1}$ lie just
slightly below the values of the lower bounds obtained with standard method
and the basis $\Psi_{2}$, however one has to keep in mind that in this case
$|\Psi_{1}|=11$ and $|\Psi_{2}|=286$. Moreover, we see that the HRR method
with a recursion depth $I=1$ performs just as well as the (full depth)
reinforced regression method ($I=J=9$).

Figure 5. The elapsed CPU times during the backwards induction and calculation
of the lower bounds from Table 1, plotted with respect to the number of
underlying assets $d$.

Furthermore, in Figure 5 we have visualized the corresponding elapsed CPU
times during the backwards induction and the calculation of the lower bounds.
As foreshadowed in Section 5, we see that the computational costs of the HRR
method are significantly reduced by choosing a small recursion depth $I$. In
particular, we are able to state that for sufficiently large $d$ ($d\geq 5$
respectively $d\geq 3$) the HRR method with recursions depth $I=1$ and the
basis $\Psi_{1}$ respectively $\Psi_{2}$ is more efficient then the standard
method with the basis $\Psi_{2}$ respectively $\Psi_{3}$.

Figure 6. Visualization of the lower bounds ($\pm$ 99.7% Monte-Carlo error) of
the values of Bermuda max-call options with $J=9,18,36,72$ exercise dates and
$d=4$, $T=1$, $y_{max}=1$, $x_{0}=C=100$, $r=0.05$, $\delta=0.1$,
$\sigma=0.2$. The values are obtained with the standard regression method
$I=0$ and the HRR method $I=1,2,9$. For all methods we used $M=10^{6}$
training sample paths and $M_{\mathrm{test}}=10^{6}$ sample paths for the
calculation for the lower bounds.

The results of the second set-up are visualized in Figure 6. In this case, we
have approximated the value of Bermudan max-call options with a fixed number
of assets $d=4$ and different numbers of exercise dates
$J\in\\{9,18,36,72\\}$, while also keeping $T=1$ fixed. When keeping all other
parameters fixed, the value of the option clearly is non-decreasing in the
number of exercise dates $J$. Our first observation is that for all considered
methods there exists a threshold for $J$ at which the performance worsens.
Indeed, Corollary 6.3 shows that the approximation error depends exponentially
on $J$.

However, in practice, some methods are less vulnerable to the error explosion
in $J$ than others. In this example, we see that the standard regression
method with the basis $\Psi_{1,g}$ and $\Psi_{2}$, respectively, starts to
perform worse for $J\geq 36$. The lower bounds calculated with the HRR method
with the basis $\Psi_{1}$ and $I=1$ stay approximately on the same level as
the lower bounds calculated with the standard method and the basis $\Psi_{2}$,
for all numbers of exercise dates. The lower bounds calculated with the
standard method and the basis $\Psi_{3}$ first increase when moving from $9$
to $18$ exercise dates and decrease at last when moving from $36$ to $72$
exercise dates.

The main observation is that when we increase $J$, the HRR methods with $I=2$
and $I=9$ come closer to the lower bounds calculated with the standard method
and the basis $\Psi_{3}$. This underlines the theoretical discussion for
$J\to\infty$ from Section 6. Moreover, we see that the HRR method performs at
least as well with $I=2$ as with $I=9$. Regarding CPU time, the HRR method
with $I=2$ is more efficient than then the standard regression method with the
basis $\Psi_{3}$ for $J\geq 9$. Comparing the HRR methods with $I=2$ and
$I=9$, we see that choosing the parameter $I$ small is necessary to in order
to obtain desirable efficiency. Since on the other hand, the HRR method with
$I=1$ performed significantly worse then with $I=2$, we also see that in this
case it was necessary to choose $I>1$. These observations underline the
relevance of the HRR in its full complexity even in the case of optimal
stopping problems.

####  7.1.2. Multiple exercise rights

Next we consider a Bermudan max-call option with $y_{max}=4$ exercise rights.
In this case the HRR method allows for different possibilities of reinforced
value functions depending on the choice of the sets $\mathcal{L}^{y}$ (recall
Remark 4.2). Since $y_{max}$ is small, we choose
$\mathcal{L}^{y}\equiv\\{1,2,3,4\\}$ for simplicity.

Basis | Regression | Hierarchical Reinforced Regression
---|---|---
$I=0$ | $I=1$ | $I=2$ | $I=3$ | $I=5$
$\Psi_{1}$ | 90.863 (0.072) | 92.038 (0.070) | 92.287 (0.070) | 92.311 (0.067) | 92.357 (0.061)
$\Psi_{1,g}$ | 91.837 (0.082) | - | - | - | -
$\Psi_{2}$ | 92.140 (0.070) | 92.418 (0.064) | 92.548 (0.060) | 92.631 (0.061) | 92.625 (0.061)
$\Psi_{3}$ | 92.571 (0.069) | - | - | - | -
Table 2. Lower bounds ($\pm$ 99.7% Monte-Carlo error) for the value of the
Bermudan max-call option with data $J=24$, $T=2$, $y_{max}=4$, $x_{0}=C=100$,
$d=5$, $r=0.05$, $\delta=0.1$, $\sigma=0.2$. For all methods we used
$M=10^{6}$ training sample paths and $M_{\mathrm{test}}=10^{7}$ paths for
calculating the lower bound. An upper bound to the value, calculated with the
dual approach from [Sch12] and [BSZ15], is given by $92.971$ $(0.043)$, with
the HRR method with $I=3$ and basis $\Psi_{2}$ using $10^{5}$ outer and
$10^{3}$ inner sample paths.

In Table 3 we present lower bounds to the value of a Bermuda max-call option
with $y_{max}=4$ exercises rights, obtained with the standard regression
method and HRR method for different choices of regression basis functions and
the parameter $I$, with the implementation of the second type described in
Section 4.2. We first observe that for a fixed set of basis functions
$\Psi_{1}$ or $\Psi_{2}$ increasing the parameter $I$ yields increased, and
thus improved, lower bounds. This improvement is most significant when moving
from $I=0$ (standard regression) to $I=1$ and from $I=1$ to $I=2$ and becomes
less significant when further increasing $I$. Moreover, we observe that the
HRR method with $I=1$ and basis functions $\Psi_{1}$ yields better lower
bounds then the standard regression method with the larger set of basis
functions $\Psi_{1,g}$ and more importantly, for $I\geq 2$ the HRR with basis
functions $\Psi_{1}$ method yields better lower bounds than the standard
regression method with the even larger set of basis functions $\Psi_{2}$. This
observation prevails when comparing the standard regression method with the
basis $\Psi_{3}$ against the HRR method with the basis $\Psi_{2}$. We can
therefore conclude that the HRR method yields results of better quality than
standard regression using fewer regression basis functions. Moreover, we
realize that up to changes that are insignificant with respect to the Monte
Carlo error, the HRR reaches its best performance already for $I=3$, thus
further increasing $I$ is not necessary.

###  7.2. A gas storage problem

In this subsection we consider a gas-storage problem of the kind introduced in
Example 2.4. In contrast to the example in the previous subsection, this
optimal control problem is not of a multiple stopping type, which is a
consequence of the anti-symmetry in the policy set: injection of gas into the
facility $(a=1)$, no action $(a=0)$ and production of gas $(a=-1)$.

For the gas price we use a similar but slightly more elaborate model to the
one proposed in [TDR09] (and also used in [GHW11]). More specifically, we use
the following joint dynamics to model the _price of crude oil_ $X^{1}$ and the
_price of natural gas_ $X^{2}$

(7.1)
$\begin{split}\mathrm{d}X^{1}(t)&=\alpha_{1}(\beta-X^{1}(t))\mathrm{d}t+\sigma_{1}X^{1}(t)\mathrm{d}W^{1}(t)+\left(J^{1}_{N(t-)+1}-X^{1}(t)\right)\mathrm{d}N(t)\\\
\mathrm{d}X^{2}(t)&=\alpha_{2}(X^{1}(t)-X^{2}(t))\mathrm{d}t+\sigma_{2}X^{2}(t)\mathrm{d}W^{2}(t)+\left(J^{2}_{N(t-)+1}-X^{2}(t)\right)\mathrm{d}N(t),\end{split}$

for $0\leq t\leq T$, where $\beta,\alpha_{i},\sigma_{i}>0$ for $i=1,2$,
$W^{1}$ and $W^{2}$ are Brownian motions with correlation $\rho_{W}\in[0,1]$,
$N$ is a Poisson process with intensity $\lambda>0$ and $(J_{k})_{k=1,...}$
are i.i.d. normal distributed random vectors with
$J^{i}_{1}\sim\mathcal{N}(\mu_{i},\eta_{i}^{2})$, $\mu_{i},\eta_{i}>0$ and
$\rho_{J}=\mathrm{Cor}(J_{1}^{1},J_{1}^{2})\in[0,1]$. Moreover we assume that
$(W^{1},W^{2})$, $N$ and $(J^{1},J^{2})$ are independent. Note that both
$X^{1}$ and $X^{2}$ are mean reverting processes with jump contributions. The
oil price process $X^{1}$ reverts to the long-term constant mean $\beta$ and
the gas price process $X^{2}$ reverts towards the oil price $X^{1}$, which is
aiming to model the well known strong correlation between crude oil and
natural gas prices. Note also that we have assumed for simplicity that the
jump signal, which has the purpose of modeling price peaks, is the same
Poisson process for both oil and gas prices, however the magnitude of the
jumps is given by different (but correlated) normal distributed random
variables.

Denote by $(\widetilde{X}_{j})_{j=1,...,365}$ the 2-dimensional Markov chain
that is obtained by discretizing the above SDE (7.1) with an Euler-scheme on
the time interval interval $[0,1]$. We assume that the manager of the gas
storage facility has the possibility to buy and sell gas on a predefined set
of dates in the year $\\{t_{j}\\}_{j=1,...,J}\subset\\{1,...,365\\}$ with
$t_{j}=j\cdot\delta t$ and some $\delta t,J\in\mathbb{N}$ such that $\delta
t\cdot J\leq 365$. The 2-dimensional Markov chain underlying the optimal
control problem is then given by $X=(X_{j})_{j=0,...,J}$ with
$X_{j}:=\widetilde{X}_{t_{j}}$.

Recall from Example 2.4 that we assume that the volume of gas in the storage
can only be increased or decreased by a fraction $\Delta=1/N$ for some
$N\in\mathbb{N}$ over the time interval of $\delta t$ days. The state space of
the control variable is then given by
$\mathcal{L}=\\{0,\Delta,2\Delta,...,1\\}$. Also recall the definition of the
space of policies $\mathcal{K}$, the constraint sets $K_{j}$ and the function
$\varphi_{j}$ from Example 2.4. We assume that there is no trading at $j=0$
hence $K_{0}\equiv\\{0\\}$. The cash-flow underlying to the optimal control
problem only depends on the second component of the Markov chain $X_{j}$ and
is given by

$H_{j}(a,y,x)=-a\cdot\Delta\cdot x^{2}\cdot e^{-rj(\delta t/365)},\quad
a\in\mathcal{K},\,y\in\mathcal{L},\;j=1,...,J,$

where $r>0$ is the interest rate.

We do not pay attention to the physical units of the parameters quantifying
the gas storage capacity and the quotation of the gas price, since the
linearity of the pay-off with respect to the parameter $\Delta$ and the gas
price $X^{2}_{j}$ allows to properly scale the resulting value of the optimal
control problem. The following specific choice of the price model parameters
are oriented at the values in [TDR09]

(7.2)
$\begin{array}[]{c}\beta=45,\quad\alpha_{1}=0.25,\quad\alpha_{2}=0.5,\quad\sigma_{1}=\sigma_{2}=0.2,\quad\rho_{W}=0.6,\\\
\quad\lambda=2,\quad\mu_{1}=\mu_{2}=100,\quad\eta_{1}=\eta_{2}=30,\quad\rho_{J}=0.6.\end{array}$

Figure 7 shows a sample trajectory of the Markov chain $X$ with the above
parameters.

Figure 7. A sample path of the Markov chain
$(X_{j})_{j=0,...,J}=(\widetilde{X}_{t_{j}})_{j=0,...,J}$ where $t_{j}=j\cdot
7$ and $J=52$. The approximation
$\widetilde{X}=(\widetilde{X}_{j})_{j=1,...,365}$ to the SDE (7.1) is
simulated with the parameters given in (7.2) and
$\widetilde{X}_{0}=(100,100)$. $X^{1}$ and $X^{2}$ serve as models for the
prices of crude oil and natural gas.

Further we define the following sets of polynomial regression basis functions

(7.3)
$\begin{split}P_{i}(X^{2})&:=\left\\{(x_{1},x_{2})\mapsto(x_{2})^{p}\;\big{|}\;p=0,...,i\right\\}\\\
P_{i}(X^{1},X^{2})&:=\left\\{(x_{1},x_{2})\mapsto(x_{1})^{p}(x_{2})^{q}\;\big{|}\;p,q=0,...,i,\;\;p+q\leq
i\right\\}.\end{split}$

I | Basis | $v_{0}(Y_{0},X_{0})$ | Lower bounds
---|---|---|---
0 | $P_{1}(X^{2})$ | 78.381 | 70.489 (0.066)
$P_{1}(X^{1},X^{2})$ | 78.575 | 70.635 (0.068)
$P_{2}(X^{2})$ | 73.072 | 71.253 (0.068)
$P_{2}(X^{1},X^{2})$ | 73.207 | 71.402 (0.068)
$P_{3}(X^{1},X^{2})$ | 72.929 | 71.333 (0.081)
| $P_{4}(X^{1},X^{2})$ | 72.595 | 71.498 (0.068)
1 | $P_{1}(X^{1},X^{2})$ | 71.991 | 71.579 (0.070)

Table 3. Approximate values and lower bounds for the gas storage problem with
parameters given in (7.4) and price model parameters given in (7.2). Note that
- although seemingly so - the estimate $v_{0}$ not necessarily presents an
upper bound to the true value and is included in the table only for
verification purposes.

The quantities were obtained with the standard regression method ($I=0$) and
the HRR method ($I=1$), the different sets of basis functions (7.3),
$M=10^{5}$ training sample paths and $M_{\mathrm{test}}=10^{6}$ paths for
calculating the lower bounds.

We have approximated the value of the gas storage problem with the following
parameters

(7.4) $\delta t=7,\quad J=52,\quad\Delta=1/8,\quad X_{0}=(100,100),\quad
Y_{0}=4/8,\quad r=0.1.$

In this configuration the gas storage facility is initially loaded with half
its capacity and the gas storage manager has the possibility to trade gas
every seven days, and the amount by which the manager can inject or produce
gas is one height of the total capacity. In Table 3 we present the numerical
results that were obtained with the standard regression method and the HRR
method.

We used $M=10^{5}$ training sample paths and $M_{test}=10^{6}$ sample paths
for the calculation of the lower bounds. For the set of reinforced basis
functions in the HRR method we have chosen $\mathcal{L}^{y}\equiv\\{Y_{0}\\}$,
i.e. in each step of the backwards induction, the regression basis was
reinforced with only one function.

We observe at first, that the lower bounds obtained with the standard
regression method are improved when using polynomials in both variables
$(X^{1},X^{2})$ instead of just in the second variable $X^{2}$ (gas price) and
are also improved when using polynomials of increasing order (with the only
exception of the third degree polynomials). Moreover, we observe that the
lower bound obtained with the HRR method, using the set of basis functions
$P_{1}(X^{1},X^{2})$ and $I=1$, lies above all lower bounds that were obtained
with the standard regression, in particular the bound obtained with the
regression basis $P_{4}(X^{1},X^{2})$ (up to Monte Carlo errors). Hence the
HRR method based on polynomials of degree one performed at least as well as
the standard method with polynomials of degree four.

###  7.3. Conclusions

Let us summarize the findings of the numerical experiments. We observe that
the hierarchical reinforced regression algorithm (HRR) based on polynomial
basis functions of a certain degree $\deg$ tend to produce results comparable
to standard regression (SR) based on polynomial basis functions of degree
$\deg+1$ or even higher, see Figures 4, 6, Tables 2, and, most impressively,
3.

The numerical results also indicate that, indeed, HRR with low depth of the
hierarchy $I$ already performs very well, even if $I\ll J$, see Figures 4, 6
and Table 2. Hence, HRR performs with similar accuracy to the reinforced
regression algorithm (RR) of [BSSZ20], but at much improved cost.
Additionally, when comparing HRR with SR at fixed accuracy, the computational
cost of HRR is usually much smaller, especially for $d$ large, see Figure 5.

Finally, we note that the accuracy of the HRR method increases substantially
when the time discretization is refined, i.e., when $J$ is increased for fixed
time horizon $T$. This theoretically very plausible observation (see Section
6) is backed up by numerical experiments, see Figure 6.

## References

  * [AB04] L. Andersen and M. Broadie.  A Primal-Dual Simulation Algorithm for Pricing Multi-Dimensional American Options.  Management Science, 50(9):1222–1234, 2004.
  * [Åst12] Karl J Åström.  Introduction to stochastic control theory.  Courier Corporation, 2012.
  * [BCJ19] Sebastian Becker, Patrick Cheridito, and Arnulf Jentzen.  Deep optimal stopping.  Journal of Machine Learning Research, 20:74, 2019.
  * [BCJW21] Sebastian Becker, Patrick Cheridito, Arnulf Jentzen, and Timo Welti.  Solving high-dimensional optimal stopping problems using deep learning.  European Journal of Applied Mathematics, 32(3):470–514, 2021.
  * [BG04] Mark Broadie and Paul Glasserman.  A stochastic mesh method for pricing high-dimensional American options.  Journal of Computational Finance, 7:35–72, 2004.
  * [BKS09] D. Belomestny, A. Kolodko, and J. Schoenmakers.  Regression methods for stochastic control problems and their convergence analysis.  SIAM Journal on Control and Optimization, 48:3562–3588, 01 2009\.
  * [BR11] Nicole Bäuerle and Ulrich Rieder.  Markov decision processes with applications to finance.  Springer Science & Business Media, 2011.
  * [BRS21] Christian Bayer, Martin Redmann, and John Schoenmakers.  Dynamic programming for optimal stopping via pseudo-regression.  Quantitative Finance, 21(1):29–44, 2021.
  * [BS18] Denis Belomestny and John Schoenmakers.  Advanced simulation-based methods for optimal stopping and control.  Palgrave Macmillan, London, 2018.  With applications in finance.
  * [BS20] Denis Belomestny and John Schoenmakers.  Optimal stopping of McKean–Vlasov diffusions via regression on particle systems.  SIAM Journal on Control and Optimization, 58(1):529–550, 2020.
  * [BSSZ20] Denis Belomestny, John Schoenmakers, Vladimir Spokoiny, and Bakhyt Zharkynbay.  Optimal stopping via reinforced regression.  Communications in Mathematical Sciences, 16(1):109–121, 2020.
  * [BSZ15] Christian Bender, John Schoenmakers, and Jianing Zhang.  Dual representations for general multiple stopping problems.  Math. Finance, 25(2):339–370, 2015.
  * [GHW11] Lajos Gyurko, B. Hambly, and Jan Witte.  Monte Carlo methods via a dual approach for some discrete time stochastic control problems.  Mathematical Methods of Operations Research, 81, 12 2011.
  * [GKKW02] László Györfi, Michael Kohler, Adam Krzyżak, and Harro Walk.  A distribution-free theory of nonparametric regression.  Springer Series in Statistics. Springer-Verlag, New York, 2002.
  * [HE16] Jiequn Han and Weinan E.  Deep learning approximation for stochastic control problems, 2016.
  * [Pha09] Huyên Pham.  Continuous-time stochastic control and optimization with financial applications, volume 61.  Springer Science & Business Media, 2009.
  * [Rog02] L. C. G. Rogers.  Monte Carlo valuation of American options.  Math. Finance, 12(3):271–286, 2002.
  * [SB18] Richard S Sutton and Andrew G Barto.  Reinforcement learning: An introduction.  MIT press, 2018.
  * [Sch12] J. Schoenmakers.  A pure martingale dual for multiple stopping.  Finance Stoch., 16:319–334, 2012.
  * [TDR09] Matt Thompson, Matt Davison, and Henning Rasmussen.  Natural gas storage valuation and optimization: A real options application.  Naval Research Logistics (NRL), 56(3):226–238, 2009.
  * [TVR01] John Tsitsiklis and Benjamin Van Roy.  Regression methods for pricing complex American style options.  IEEE transactions on neural networks / a publication of the IEEE Neural Networks Council, 12:694–703, 02 2001.
  * [Zan13] Daniel Z. Zanger.  Quantitative error estimates for a least-squares Monte Carlo algorithm for American option pricing.  Finance Stoch., 17(3):503–534, 2013.
  * [Zan18] Daniel Z. Zanger.  Convergence of a least-squares Monte Carlo algorithm for American option pricing with dependent sample data.  Math. Finance, 28(1):447–479, 2018.