Pricing American options under rough volatility using deep-signatures and signature-kernels

Replacing the continuous time interval $[0,T]$ by some finite grid $0=t_{0}<t_{1}<\dots<t_{N}=T$, the dynammic programming principle [PS06, Page eq. (2.1.7)] for the discrete optimal stopping reads

$$Y_{t_{N}}=Z_{t_{N}},\quad Y_{t_{n}}=\max(Z_{t_{n}},\mathbb{E}[Y_{t_{n+1}}|\mathcal{F}^{X}_{t_{n}}]),\quad 0\leq n\leq N-1,$$

(2)

and $\tau^{\star}=\inf\{t_{n}:Z_{t_{n}}\geq\mathbb{E}[Y_{t_{n+1}}|\mathcal{F}^{X}_{t_{n}}]\}$ is optimal. Motivated by the famous Longstaff and Schwartz algorithm [LS01], the first optimal stopping time can be obtained as $\tau_{0}$ of the following recursion

$$\tau_{N}=t_{N},\quad\tau_{n}=t_{n}1_{\{Z_{t_{n}}\geq\mathbb{E}[Z_{\tau_{n+1}}|\mathcal{F}^{X}_{t_{n}}]\}}+\tau_{n+1}1_{\{Z_{t_{n}}<\mathbb{E}[Z_{\tau_{n+1}}|\mathcal{F}^{X}_{t_{n}}]\}},\quad 0\leq n\leq N-1.$$

(3)

The remaining question is, how to compute the continuation values $\mathbb{E}[Z_{\tau_{n+1}}|\mathcal{F}^{X}_{t_{n}}]$ for possibly non-Markovian state-processes $X$, which will be the main topic of the forthcoming Section 3. For now, assume we are given a suitable family of basis-functions $(\psi_{k})$, such that $\mathbb{E}[Z_{\tau_{n+1}}|\mathcal{F}^{X}_{t_{n}}]\approx\sum_{k}\alpha_{k}^{n}\psi_{k}(X|_{[0,t_{n}]})$ with some coefficients $(\alpha^{n}_{k})$. A general version of the Longstaff-Schwartz algorithm goes as follows:

1. 1.

   Draw $i=1,\dots,M$ samples $X^{(i)}$ and $Z^{(i)}$ and initialize $\tau_{N}^{(i)}\equiv t_{N}$.

2. 2.

   Then, for $n=N-1,\dots,1$

   • •

     solve the minimization problem

     $$\alpha^{\star,n}:=\mathop{\mathrm{argmin}}_{\alpha}\frac{1}{M}\sum_{i=1}^{M}\left(Z^{(i)}_{\tau_{n+1}^{(i)}}-\sum_{k}\alpha_{k}\psi_{k}^{n}(X|_{[0,t_{n}]}^{(i)})\right)^{2}$$

   • •

     Set $\tau_{n}^{(i)}=t_{n}$ if $Z^{(i)}_{t_{n}}\geq\sum_{k}\alpha^{\star,n}_{k}\psi_{k}^{n}(X|_{[0,t_{n}]}^{(i)})$, and $\tau_{n}^{(i)}=\tau_{n+1}^{(i)}$ otherwise.

3. 3.

   Finally, draw $i=1,\dots,M^{\prime}$ independent sample-paths $\tilde{X}^{(i)}$ and $\tilde{Z}^{(i)}$, and compute the stopping times $\tilde{\tau}^{(i)}=\inf\left\{t_{n}:\tilde{Z}^{(i)}_{t_{n}}\geq\sum_{k}\alpha^{\star,n}_{k}\psi_{k}^{n}(\tilde{X}|_{[0,t_{n}]}^{(i)})\right\}$ and the lower-biased estimator

   $$\tilde{y}_{0}=\max\left(Z_{t_{0}},\frac{1}{M^{\prime}}\sum_{i=1}^{M^{\prime}}\tilde{Z}^{(i)}_{\tilde{\tau}^{(i)}}\right)$$

Notice that the resimulation in the final step ensures, that the $\tilde{\tau}^{(i)}$ are indeed stopping times, since the computed coefficients $\alpha$ are now independent of the samples, and thus $\max(Z_{t_{0}},\mathbb{E}[Z_{\tilde{\tau}_{1}}])$ is a true lower-bound. The Monte-Carlo approximation $\tilde{y}_{0}$ is therefore lower-biased and strictly speaking not a true lower-bound, due to the Monte-Carlo error with respect to $M^{\prime}$. Nevertheless, this error is typically chosen to be very small, and we will refer to $\tilde{y}_{0}$ as lower-bound anyways.

While the Longstaff-Schwartz procedure described in the last section provides us with lower-bounds for the optimal value, it is desirable to additionally have an algorithm producing upper-bounds. To this end, it is useful to consider the (discretized) dual formulation of (1) due to [Rog02]

$$\inf_{M\in\mathcal{M}^{2}}\mathbb{E}[\max_{0\leq k\leq N}(Z_{t_{k}}-M_{t_{k}})],$$

(4)

where the minimization is over discrete $L^{2}$-martingales with respect to the filtration $\mathcal{G}_{n}=\mathcal{F}^{X}_{t_{n}},n=0,\dots,N$. By the Martingale Representation Theorem [KS91, Theorem 4.5], if the underlying filtration $(\mathcal{F}^{X}_{t})$ is generated by a Brownian motion $W$, such martingales are of the form $M_{t_{k}}^{\alpha}=\int_{0}^{t_{k}}\alpha_{s}dW_{s}$ for some $(\mathcal{F}_{t}^{X})-$progressive process $\alpha$. On the other hand, due to the progressive measurability, the integrands are of the form $\alpha_{t}=f(X|_{[0,t]})$ for some (deterministic) function $f$. Assume again that we are given some family of basis-functions $(\psi_{k})$ such that $f(X|_{[0,t]})\approx\sum_{k}\beta_{k}\psi_{k}(X|_{[0,t]})$ for some coefficients $\beta$. If the family $(\psi_{k})$ is rich enough, it is reasonable to parametrize the space $\mathcal{M}^{2}$ by the span of basis-martingales $\{M^{k}_{t}=\int_{0}^{t}\psi_{k}(X|_{[0,s]})dW_{s}:k\geq 0\}$, which leads to the following dual algorithm:

1. 1.

   Draw $i=1,\dots,M$ sample-paths of $Z^{(i)},X^{(i)},W^{(i)}$, and compute the basis martingales $M^{k,(i)}$ (e.g. using an Euler-scheme)

2. 2.

   Solve the minimization problem (see, e.g., [DFM12, Bel13])

   $$\beta^{\star}=\mathop{\mathrm{argmin}}_{\beta}\frac{1}{M}\sum_{i=1}^{M}\max_{0\leq n\leq N}\left(Z^{(i)}_{t_{n}}-\sum_{k}\beta_{k}M^{k}_{t_{n}}\right)$$

3. 3.

   Finally, draw $i=1,\dots,M^{\prime}$ independent sample-paths $\tilde{Z}^{(i)},\tilde{X}^{(i)},\tilde{W}^{(i)}$ and compute the martingales $\tilde{M}^{k,(i)}$ and the upper-biased estimator

   $$\tilde{y}_{0}=\frac{1}{M^{\prime}}\sum_{i=1}^{M}\max_{0\leq n\leq N}\left(\tilde{Z}^{(i)}_{t_{n}}-\sum_{k}\beta^{\star}_{k}\tilde{M}^{k,(i)}_{t_{n}}\right)$$

Similar as in the primal case, the independent resimulations ensure that we get true martingales $\tilde{M}^{\star}:=\sum_{k}\beta^{\star}_{k}\tilde{M}^{k}$, since $\beta^{\star}$ is independent of the samples, and therefore $\mathbb{E}[\max_{k}(Z_{t_{k}}-\tilde{M}^{\star}_{t_{k}})]$ is a true upper-bound. By the same abuse of terminology as in the primal case, we call the Monte-Carlo approximation $\tilde{y}_{0}$ upper-bound, keeping in mind that it is strictly speaking only upper-biased due to the Monte-Carlo error with respect to $M^{\prime}$.

One of the most important algebraic properties of the signature is the following: For any two linear functionals on the tensor-algbera, that is $\ell_{1},\ell_{2}\in T((\mathbb{R}^{d})^{\star})$, we can find $\ell_{3}$ (in fact, a closed-form expression can be given), such that $\big{\langle}\ell_{1},\mathbbm{X}^{<\infty}_{0,T}\big{\rangle}\big{\langle}\ell_{2},\mathbbm{X}^{<\infty}_{0,T}\big{\rangle}=\big{\langle}\ell_{3},\mathbbm{X}^{<\infty}_{0,T}\big{\rangle}$, for all signatures $\mathbbm{X}^{<\infty}_{0,T}$. In words, the space of linear functionals of the signature forms an algbera, which is point-seperating (since we assume to have a unique signature, see Remark 3.1). An application of a general Stone-Weierstrass result, see [CO22, Section 2.1], tells us that the set of linear functional of the robust signature is dense in the set of bounded continuous functions on the corresponding path space. This fact can be used to deduce a global $L^{p}$-approximation for linear functionals of the signature, see [BPS23, Theorem 2.8], which then yields the following convergence results. We will always assume to be in the framework described in the beginning of Section [BPS23, Section 3.1], where certain assumptions on the path space, as well as on the payoff process are made precise. The following proposition summarizes the results [BPS23, Proposition 3.3 and 3.8], to which we refer for more precise statements and detailed proofs.

While the last result ensures converges for the primal and dual procedures, is does not come with quantitative statements about convergence rates. Additionally, these algorithms can become computationally expensive, as the size of the signature grows exponentially with $K$. In [BPS23] several numerical experiments were performed for these two algorithms in non-Markovian frameworks, in particular for models driven by fractional Brownian motion with small Hurst parameters. The accuracy of the method can be measured by the duality gap between lower and upper bounds, and in some examples these gaps were observed to be quite significant even for high truncation levels. Two important sources for this gap can be described as follows: First, especially in very rough regimes, the considered signature levels might not suffice to capture the relevant past of the non-Markovian processes, so that information is lost when truncating the signature. Secondly, the functionals on the path spaces we try to learn are highly non-linear, so that more sophisticated learning technologies are required to improve the performance. The goal of the next sections is to introduce non-linear extensions of the primal and dual algorithm, based on deep-, resp. kernel-learning methodologies, with the objective of improving the duality gap.

Applying Deep Neural Networks (DNNs) on top of the signature was considered several times in the literature, see for instance [KBPA⁺19], where even more generally, the signature can act as a layer in a network. Here, we simply consider the truncated signature as the input for a DNN of the form

$$\theta=\beta\circ\sigma\circ A_{1}\circ\sigma\cdots\circ A_{I},$$

where $\beta,A_{0},\dots,A_{I}$ are affine, $\beta:\mathbb{R}^{m}\rightarrow\mathbb{R},A_{I}:T(\mathbb{R}^{d})\rightarrow\mathbb{R}^{m},A_{I-i}:\mathbb{R}^{m}\to\mathbb{R}^{m}$, for $i=1,\dots,I-1$, and $\sigma$ is an activation function. We call $\theta$ a signature DNN with $I$ hidden layers, each of which having $m$ neurons. The signature map $\mathbf{x}|_{[0,t]}\mapsto\mathbbm{x}^{<\infty}_{0,t}$ represents the input layer, and the first hidden layer acting on the signature, can be seen as a vector of linear functionals $A_{I}=(\ell_{1},\dots,\ell_{m})$, and $A_{I}(\mathbbm{x}^{<\infty})=(\big{\langle}\ell_{1},\mathbbm{x}^{<\infty}\big{\rangle},\dots,\big{\langle}\ell_{m},\mathbbm{x}^{<\infty}\big{\rangle})^{\top}\in\mathbb{R}^{m}$. The remaining hidden layers are all of the form $A_{j}(x)=A_{j}x+b_{j}$, where $A_{j}\in\mathbb{R}^{m\times m},b_{j}\in\mathbb{R}^{m}$, and the output layer is given by $\beta\in\mathbb{R}^{m}$. We denote by $\mathrm{DNN}^{\sigma,I}_{\mathrm{sig}}$ the set of all such signature DNNs $\theta$ with $I$ hidden layers, and activation function $\sigma$.

It is well known, that already $I=1$ hidden layer DNNs are universal, see, e.g. [LLPS93]. For this reason, and to ease the notation for the theoretical results, in the remainder of this section we fix $I=1$. Nevertheless, all the results trivially extend to multiple hidden layers $I>1$, and in the numerical experiments we mostly choose more than one hidden layer. Setting $\mathrm{DNN}^{\sigma}_{\mathrm{sig}}:=\mathrm{DNN}^{\sigma,1}_{\mathrm{sig}}$, we can explicitly write the set as

$$\mathrm{DNN}^{\sigma}_{\mathrm{sig}}=\left\{\mathbf{x}|_{[0,t]}\mapsto\sum_{j=1}^{m}\beta_{j}\sigma\Big{(}\big{\langle}\ell_{j},\mathbbm{x}^{<\infty}_{0,t}\big{\rangle}\Big{)}:(\beta,\ell)\in\mathbb{R}^{m}\times\mathcal{W}^{m},m\geq 1\right\},$$

(7)

where $\mathcal{W}$ denotes the set of all linear functionals, and $\mathbf{x}|_{[0,t]}$ are stopped $\alpha$-Hölder rough paths, see [BPS23, Section 2] for details. For fixed $(m,k)\in\mathbb{N}^{2}$ and any pair $(\beta,\ell)\in\mathbb{R}^{m}\times(\mathcal{W}^{\leq k})^{m}$, where $\mathcal{W}^{\leq k}$ are linear functionals on the $k$-step signature $\mathbbm{x}^{\leq k}$, denote by $\theta^{(\beta,\ell)}$ the corresponding functional in the set $\mathrm{DNN}^{\sigma}_{\mathrm{sig}}$. Revisiting the Longstaff & Schwartz algorithm presented in Section 2.1, it is tempting to define a sequence of deep stopping times $(\tau^{(m,k)}_{n}:1\leq n\leq N)$, recursively defined similar to (3), but with conditional expectations replaced by $\theta^{(\beta_{P}^{\star,n},\ell_{P}^{\star,n})}(\mathbf{X}|_{[0,t_{n}]})$, where for $0\leq n\leq N$

$$(\beta_{P}^{\star,n},\ell_{P}^{\star,n})=\mathop{\mathrm{argmin}}_{(\beta,\ell)\in\mathbb{R}^{m}\times(\mathcal{W}^{\leq k})^{m}}\mathbb{E}\left[\big{(}Z_{\tau^{(k,m)}_{n+1}}-\theta^{(\beta,\ell)}(\mathbf{X}|_{[0,t_{n}]})\big{)}^{2}\right].$$

(8)

In words, we recursively learn the continuation values as neural networks of truncated signatures. Similarly, we approximate the Doob-martingale from Section 2.2 by deep martingales $M^{(k,m)}_{t}=\int_{0}^{t}\theta^{(\beta_{D}^{\star},\ell_{D}^{\star})}(\mathbf{X}|_{[0,s]})dW_{s}$, where

$$(\beta_{D}^{\star},\ell_{D}^{\star})=\mathop{\mathrm{argmin}}_{(\beta,\ell)\in\mathbb{R}^{m}\times(\mathcal{W}^{\leq k})^{m}}\mathbb{E}\left[\max_{0\leq n\leq N}\bigg{(}Z_{t_{n}}-\int_{0}^{t_{n}}\theta^{(\beta,\ell)}(\mathbf{X}|_{[0,s]})dW_{s}\bigg{)}\right].$$

(9)

Thanks to the universality of both signatures and DNNs, the following result should not come as a surprise. Similar as in the last section, we adopt the assumptions described in the beginning of [BPS23, Section 3.1].

An outline of the proof can be found in Appendix A. Similar to the linear case, in practice we solve the sample average approximation of (8)-(9), that is

	$\displaystyle(\beta_{P}^{\star,n},\ell_{P}^{\star,n})$	$\displaystyle=\mathop{\mathrm{argmin}}_{(\beta,\ell)\in\mathbb{R}^{m}\times(\mathcal{W}^{\leq k})^{m}}\frac{1}{M}\sum_{i=1}^{M}\left(Z^{(i)}_{\tau_{n+1}^{(i)}}-\theta^{(\beta,\ell)}(\mathbf{X}^{(i)}|_{[0,t_{n}]})\right)^{2},$		(10)
	$\displaystyle(\beta_{D}^{\star},\ell_{D}^{\star})$	$\displaystyle=\mathop{\mathrm{argmin}}_{(\beta,\ell)\in\mathbb{R}^{m}\times(\mathcal{W}^{\leq k})^{m}}\frac{1}{M}\sum_{i=1}^{M}\max_{0\leq n\leq N}\left(Z^{(i)}_{t_{n}}-M^{(\beta,\ell),(i)}_{t_{n}}\right),$

for i.i.d sample-paths $Z^{(i)},\mathbf{X}^{(i)},W^{(i)}$. The latter (non-convex) optimization problems can then be solved using stochastic gradient descent. More details about this and the DNN architecture will be discussed in Section 4. Let us conclude this section with a remark.

Let us start by recalling some notions from general RKHS theory, see, e.g. [Ste08]. Given a feature map $\phi:\mathcal{X}\rightarrow\mathcal{H}$, where $\mathcal{H}$ is a Hilbert space (the so-called feature space), one can define the associated kernel

$$k:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R},\quad k(x,y):=\big{\langle}\phi(x),\phi(y)\big{\rangle}_{\mathcal{H}}.$$

(11)

If the kernel is positive-definite, there exists a unique RKHS with reproducing kernel $k$, see [Ste08, Theorem 4.21], in the sense that $k$ reproduces elements $f\in\mathcal{H}$, that is $f(x)=\big{\langle}k(\cdot,x),f\big{\rangle}_{\mathcal{H}}$, and it generates the Hilbert space, $\mathcal{H}=\overline{\mathrm{span}\{k(\cdot,x):x\in\mathcal{X}\}}$. In our case, the feature map $\phi$, maps a path $x\in\mathcal{X}$ to its path-signature $\mathbbm{x}^{<\infty}$, and therefore signature kernel is naturally defined by

$$k_{s,t}(x,y):=\big{\langle}\mathbbm{x}^{<\infty}_{0,s},\mathbbm{y}^{<\infty}_{0,t}\big{\rangle},\quad 0\leq s,t\leq T$$

where $\big{\langle}\cdot,\cdot\big{\rangle}$ is the natural extension of the inner products $(\mathbb{R}^{d})^{\otimes k}$ to the tensor algebra $T((R^{d}))$, see [CS24, Definition 2.1.1]. It has been shown in [SCF⁺21], that the signature-kernel solves a Goursat-type PDE with respect to the time variables. This kernel trick allows us to evaluate the signature-kernel, representing the whole signature, by numerically solving a second-order PDE. An important extension, designed for rougher inputs, was discussed in the recent work [LL24].

An important observation is that universality of the signature feature map is equivalent to the universality of the signature-kernel, see for instance [CO22, Section 6], but also [CS24, Chapter 2.1]. Returning to optimal stopping, this motivates us to seek for functionals in the primal and dual problems, solving the regularized minimizing problems on the RKHS

$$f^{\lambda}=\mathop{\mathrm{argmin}}_{f\in\mathcal{H}}\mathcal{L}(f)+\lambda\|f\|^{2}_{\mathcal{H}},\quad\lambda>0,$$

(12)

where the loss $\mathcal{L}$ is either the primal or the dual loss function. More precisely, in the primal case, at each exercise date $n$, the loss is simply the mean-square error $\mathcal{L}_{n}(f)=\|Z_{t_{n}}-f(X|_{[0,t_{n}]})\|_{2}^{2}$. In the dual case, it is given by $\mathcal{L}(f)=\|\max_{1\leq k\leq N}(Z_{t_{k}}-M^{f}_{t_{k}})\|_{1}$, where $M_{t}^{f}=\int_{0}^{t}f(X|_{[0,s]})dW_{s}$. Then, replacing the losses by their sampled version, the general representer theorem [SHS01, Theorem 1] reveals that the minimizers are of the form $f^{\alpha}(x)=\sum_{i=1}^{M}\alpha_{i}k(X^{(i)},x)$. Thus, let us define the minimizers

$$\alpha_{P}^{\star,n}:=\mathop{\mathrm{argmin}}_{\alpha\in\mathbb{R}^{M}}\frac{1}{M}\sum_{i=1}^{M}\left(Z^{(i)}_{\tau_{n+1}^{(i)}}-\sum_{j=1}^{M}\alpha_{j}k_{t_{n},t_{n}}(X^{(i)},X^{(j)})\right)^{2}+\lambda\|f^{\alpha}\|^{2}_{\mathcal{H}},\quad\lambda>0,$$

(13)

for all exercise date $n=1,\cdots,N-1$. Similarly,

$\displaystyle\alpha_{D}^{\star}:=\mathop{\mathrm{argmin}}_{\alpha\in\mathbb{R}^{M}}\frac{1}{M}\sum_{i=1}^{M}\max_{0\leq n\leq N}\left(Z^{(i)}_{t_{n}}-\sum_{j=1}^{M}\alpha_{j}M^{(i),(j)}_{t_{n}}\right)+\lambda\|f^{\alpha}\|^{2}_{\mathcal{H}},\quad\lambda>0,$

(14)

for the kernel-martingales $M^{(i),(j)}_{t}=\left(\int_{0}^{t}k_{s,s}(X^{j},X)dW_{s}\right)^{(i)}.$ The algorithms in Section 2 can then easily be translated to the kernel-learning framework, by replacing the minimizers in the second step accordingly.

Having mentioned the theoretical advantages of the kernel-method, evaluating the signature-kernel, which means solving the Goursat-PDE [SCF⁺21], becomes the main difficulty in this procedure. Due to the recursive nature of the regression problems in the Longstaff and Schwartz algorithm in Section 2.1, typically large sample size are required to ensure stability. Moreover, as we are especially interested in paths of low regularity, a fine discretization grid is required to solve the Goursat PDE. Combining these observations with the fact, that the kernel-ridge-regressions involve computing and inverting the Gram-matrices $(k_{t_{n},t_{n}}(X^{(i)},X^{(j)})_{i,j})$, it becomes clear that a reduction of the computational costs is required.

To this end, we propose the following approach, related to the so-called Nyström-method for kernel-learning [DMC05, WS00]. We randomly select $L\ll M$ subsamples, denoting by $I_{L}$ the set of those indices, and in both the primal and dual optimization problems we restrict the minimization to functions of the form $f(\cdot)=\sum_{i\in I^{L}}\beta_{i}k(x^{i},\cdot)$. It follows that we only need to compute the $(M\times L)$-matrices

$$\mathbf{K}_{t_{n}}=(k_{t_{n},t_{n}}(X^{j},X^{i}):i\in I_{L},j=1,\dots,M),\quad n=1,\dots,N.$$

For the primal example, that is, for the kernel ridge regression, we can additionally observe that the explicit solutions is given by

$$\alpha_{P}^{\star,n}=(\mathbf{K}_{t_{n}}^{T}\mathbf{K}_{t_{n}}+M\lambda\mathbf{R}_{t_{n}})^{-1}\mathbf{K}_{t_{n}}^{T}Z_{\tau_{n+1}}$$

where $\mathbf{R}_{t_{n}}=(k_{t_{n},t_{n}}(x^{i},x^{j}):i,j\in I_{L})$. Further details about the implementation will be presented in the next section.

The code accompanying this section can be found in https://github.com/lucapelizzari/Optimal_Stopping_with_signatures. Before discussing numerical experiments, let us specify in detail how the procedures introduced in Section 3 are implemented. All the signatures are computed using the package iisignature, see [RG18]. For the signature kernel, we rely on the PDE solvers from the package sigkernel to obtain the kernel as finite-difference solution to the Goursat PDE derived in [SCF⁺21]. For the simulation of the rough Bergomi, we use https://github.com/ryanmccrickerd/rough_bergomi related to [MP18], and for the simulation of rough Heston https://github.com/SimonBreneis/approximations_to_fractional_stochastic_volterra_equations related to [BB24].

In this section we compare all the signature-based methods to price Bermudan put options. In Table 1-2, we present the rough Bergomi model with Hurst parameters $H=0.8$, resp. $H=0.07$, and for a range of strikes $K\in\{70,80,\dots,120\}$ and $S_{0}=100$. We choose the contract duration of $T=1$ year and $J=12$ early exercise options and interest rate $r=0.05$ and correlation $\rho=-0.9$. In the first column we report the European price, that is, the price for no early exercise opportunity, $\mathbf{E}=\mathbb{E}(e^{-rT}(K-S_{T})^{+})$. The second columns correspond to lower-bounds obtained in [GMZ20], see also [BTW20]. Finally, in the remaining columns, we compare the point estimate, which is the biased-estimator from the Longstaff and Schwartz algorithm based on the training samples, with the lower- and upper-bounds obtained from independent testing samples. While the intervals in the linear case are taken from [BPS23], we derive the intervals for the deep signature using truncation level $K=4$, discretization for the signature and martingales $N=600$, and both $M=12^{18}$ samples for training and testing, to ensure stability of the procedure. For the signature-kernel, we use $M=2^{17}$ samples, with $N=240$ discretization points for solving the Goursat PDE. For each strike, we highlight the best lower-, resp. upper-bound, and in the last column present the relative duality gap with respect to these two values. In Table 3 we present the same considerations for the rough Heston model with $H=0.1,r=0.06$ and $\rho=-0.7$, a choice motivated in [BB24], with smaller regularization parameter. The trainings were repeated $20$ times and the overall Monte-Carlo error is below $0.01$. Finally, in Table 4 we compare the computational time (in seconds) of one training each, as well as the offline computation of the signature, resp. signature-kernel.

In Figure 1 we present the dependence of the lower-bounds, resp. point estimates with respect to $\lambda\in[0,2]$. The optimal choice is made for highest possible lower bounds (blue lines), so that the distance to the point estimate is reasonable, since larger difference suggest overfitting. According to this remark, we choose $\lambda=10^{-3}$ for rough Bergomi, and $\lambda=10^{-8}$ in the rough Heston model.

First and most visibly, we observe that the the deep signature method delivers significantly better upper-bounds than both other methods, with reasonable training and offline costs, see Table 4. Moreover, the deep method also mostly produces the smallest duality gap, that is the smallest distance between lower- and upper-bound. However, one can also observe that for both the point estimates and lower-bounds, the linear and kernel methods, which are based on conventional convex optimization problems, slightly outperform the deep signature method in general. It should be noted that, however, this comes at the price of offline computational costs, see Table 4, where we recall that the DNNs allow to only lift the augmented volatility process. As discussed in [BPS23], for the linear procedure it is necessary to lift the three-dimensional path $(t,X_{t},\phi(X_{t}))$ and add polynomials of the state, which explains the increased computational time in the linear case. While the signature-kernel method improves the linear approach, the complexity of deriving the Gram-matrix, i.e. solving the Goursat-PDE, see Table 4, leaves the signature-kernel procedure as the most expensive approach to solve the optimal stopping problem. To handle large data sets, it is therefore necessary to solve the complexity issue for the kernel, and a potential direction could be the random Fourier signature features by [TOS23], but this is outside the scope of this paper.

In summary, for our goal of achieving minimal duality gaps for arbitrary large sample sizes, the best possible results can be achieved by combining a Longstaff & Schwartz algorithm based on the linear signature or the signature-kernel, exploiting the simplicity of the training procedure, together with the deep-signature approach for the upper-bounds, to capture the non-linearity of the Doob martingale integrand. If one has to choose one method, we recommend the deep-signature method, simply as it shows the smallest duality gaps and overall computational times.

On the other hand, it is important to note that we observe the primal signature-kernel method to be more stable compared to deep-signatures, with respect to smaller data sets. This unsurprising advantage of the kernel-method can become important when one works with real-word rather than synthetic data, where the number of samples cannot be arbitrary large, in which case the signature-kernel might be favorable. This is illustrated in Figure 2 (a) for the rough Bergomi model with $H=0.07$, where we can observe more severe instability for the neural network training for smaller $M$, both in the sense of small lower-bounds and high training variance. The Monte-Carlo errors reflect this variance, which is obtained after independently repeating the training $20$ times for each sample size. For completeness, we show a similar plot in the dual case in Figure 2 (b), where it is confirmed again that the deep-signature method heavily outperforms the signature-kernel approach, also in every data-size regime. The variance in the dual procedure is much smaller, as we only solve one optimization problem in the training, rather than at each exercise date.