Random neural networks for rough volatility
We construct a deep learning-based numerical algorithm to solve
path-dependent partial differential equations arising in the context of rough
volatility. Our approach is based on interpreting the PDE as a solution to an
SPDE, building upon recent insights by Bayer, Qiu and Yao, and on constructing
a neural network of reservoir type as originally developed by Gonon,
Grigoryeva, Ortega. The reservoir approach allows us to formulate the
optimisation problem as a simple least-square regression for which we prove
theoretical convergence properties.
Antoine Jacquier, Zan Zuric
# Random neural networks for rough volatility

###### Abstract.

We construct a deep learning-based numerical algorithm to solve path-dependent partial differential equations arising in the context of rough volatility. Our approach is based on interpreting the PDE as a solution to an SPDE, building upon recent insights by Bayer, Qiu and Yao, and on constructing a neural network of reservoir type as originally developed by Gonon, Grigoryeva, Ortega. The reservoir approach allows us to formulate the optimisation problem as a simple least-square regression for which we prove theoretical convergence properties.

Key words and phrases:Rough volatility, SPDEs, neural networks, reservoir computing 2020 Mathematics Subject Classification: 60G22, 35K10, 65C20, 68T07, 91G60 AJ acknowledges financial support from the EPSRC/T032146 grant. ZZ is supported by the EPSRC/S023925 CDT in Mathematics of Random Systems: Analysis, Modelling and Simulation. We would like to thank Lukas Gonon, Christian Bayer and Jinniao Qiu for helpful discussions. The python code is available at ZuricZ/RWIN_PDE_solver.
## 1. Introduction

In recent years, a fundamental shift from classical modelling towards so-called rough stochastic volatility models has happened. These "rough" models were first proposed by Gatheral, Jusselin, Rosenbaum [28] and by Bayer, Gatheral, Friz [4], and have since sparked a great deal of research, because of their ability to capture stylised facts of volatility time series and of option prices more accurately, while remaining parsimonious. In essence, they are a class of continuous-path stochastic volatility models, where the instantaneous volatility is driven by a stochastic process with paths rougher than those of Brownian Motion, typically modelled by a fractional Brownian motion [51] with Hurst parameter \(H\in(0,1)\). The reason for this drastic paradigm shift can be found not only under the historical measure, where the roughness of the time series of daily log-realised variance estimates suggests Holder regularity of \(H\approx 0.1\), but also under the pricing measure, where rough volatility models are able to reproduce the power-law behaviour of the ATM volatility skew. Since then, a slew of papers have appeared, providing closed-form expressions for the characteristic functions of rough Heston models [22], machine learning techniques for calibration [40], microstructural foundations [21], option pricing partial differential equations (PDEs) solvers [46, 5], among others.

Dating back to the Black-Scholes equation [11], PDEs have been used to model the evolution of the prices of European-style options. However, rough volatility models give rise to a non-Markovian framework, where the value function for a European option is not deterministic anymore, but is instead random and satisfies a backward stochastic partial differential equation (BSPDE) as was shown in [5]. Moreover, even in classical diffusive models, the so-called curse of dimensionality poses a challenge when solving PDEs in high dimension; until recently, only the backward stochastic differential equation (BSDE) approach by [54] was available to tackle this, which is not really feasible in dimension beyond six.

On a positive note, machine learning methods have spread inside quantitative finance in recent years, and neural networks in particular have become a powerful tool to overcome problems in high-dimensional situations, because of their superior computational performance across a wide range of applications [15, 31, 59]; more precisely in the context of PDEs, examples of applications thereof can be found in [19, 36, 61, 46, 60, 6]. For a more thorough literature review on the use of neural networks in finance and finance-related PDEs, we refer the reader to the surveys in [7, 30].

In this paper, we focus on the works by Hure, Pham and Warin [42], and by Bayer, Qiu and Yao [5], where the classical backward resolution technique is combined with neural networks to estimate both the value function and its gradient. Not only does this approach successfully reduce the curse of dimensionality, but also appears more effective in both accuracy and computational efficiency than existing Euler-based approaches.

Besides research on numerical aspects, a lot of progress has been made on the theoretical foundations for neural network-based methods, in particular showing that they are able to approximate solutions of certain types of PDEs [20, 43, 57, 34]. These results are significant as they show that deep neural networks can be used to solve complex problems that were previously thought intractable. However, in practice, optimal parameters of any given neural network minimising a loss function ultimately have to be calculated approximately. This is usually done through some kind of stochastic gradient descent (SGD) algorithm, which inadvertently introduces an optimisation error. Because of the
non-convexity of the network's loss surface and the stochastic nature of the SGD, the optimisation error is notoriously hard to treat rigorously. One such attempt by Gonon [32] instead involves the use of neural networks in which only certain weights are trainable and the remaining are randomly fixed. This suggests that these random-weight neural networks are, in fact, capable of learning non-degenerate Black-Scholes-type PDEs without succumbing to the curse of dimensionality. Following this, we combine the classical BSDE approach [54, 14] with random-weight neural networks (RWNNs) [41, 55, 56].

Our final algorithm then reduces to a least-square Monte-Carlo, as introduced by Longstaff and Schwartz [49] (see also [1] for related applications), where the usually arbitrary choice of basis is 'outsourced' to the reservoir of the corresponding RWNN. The basis is computationally efficient and ultimately allows us to express the approximation error in terms of the number of stochastic nodes in the network. Moreover, vectorisation of the Randomised Least Square along the sampling direction allows us to evaluate the sum of outer tensor products using the einsum function (available in NumPy and PyTorch) and achieve an even greater speed-up.

To summarise, in contrast with Bayer-Qiu-Yao [5], our numerical scheme employs RWNNs as opposed to the conventional feed-forward neural networks, resulting in significantly faster training times without sacrificing the accuracy of the scheme. Moreover, this structure allows us to provide error bounds in terms of the number of hidden nodes, granting additional insights into the network's performance. Given the comparable performance of RWNNs and conventional feed-forward neural networks, we argue that this paper illuminates an essential lesson, namely that the additional complexity of deep neural networks can sometimes be redundant at the cost of precise error bounds. We note in passing that RWNNs have already been used in Finance to price American options [37], for financial data forecasting ][48], for PIDEs [34], and we refer the interested reader to [16] for a general overview of their applications in data science.

The paper is structured as follows: Section 2 provides a brief overview of Random-weight Neural Networks (RWNNs), including their key features and characteristics. In Section 3, we outline the scheme for the Markovian case and discuss the non-Markovian case in Section 4. The convergence analysis is presented in Section 5. Additionally, Section 6 presents numerical results, which highlight the practical relevance of the scheme and its performance for different models. Some of the technical proofs are postponed to Appendix C to ease the flow of the paper.

**Notations:**\(\mathbb{R}^{+}=[0,\infty)\) represents the non-negative real numbers; \(\aleph\) refers to a random neural network, defined in Section 2; for an open subset \(E\subset\mathbb{R}^{d}\), \(1\leq p\leq\infty\) and \(s\in\mathbb{N}\) we define the Sobolev space

\[\mathcal{W}^{\aleph,p}(E,\mathbb{R}^{m})\coloneqq\Big{\{}f\in L^{p}(E,\mathbb{ R}^{m}):\;\partial_{\aleph}^{\boldsymbol{\alpha}}f\in L^{p}(E,\mathbb{R}^{m}), \text{for all }|\boldsymbol{\alpha}|\leq s\Big{\}},\]

where \(\boldsymbol{\alpha}=(\alpha_{1},\ldots,\alpha_{d})\), \(|\boldsymbol{\alpha}|=\alpha_{1}+\ldots+\alpha_{d}\), and the derivatives \(\partial_{\aleph}^{\boldsymbol{\alpha}}f=\partial_{x_{1}}^{\alpha_{1}}\ldots \partial_{x_{d}}^{\alpha_{d}}f\) are taken in a weak sense.

## 2. Random-weight neural network (Rwnn)

Neural networks with random weights first appeared in the seminal works by Barron [2, 3], but a more modern version was proposed by Huang [41] under the name _Extreme learning machine_. Today these networks are known under different names: reservoir networks, random feature or random-weight networks; we choose to follow the latter as it sounds more explicit to us.
**Definition 2.1** (Neural network).: Let \(L,N_{0},\ldots,N_{L}\in\mathbb{N},\varrho:\mathbb{R}\to\mathbb{R}\) and for \(l=1,\ldots,L\) let \(w_{l}:\mathbb{R}^{N_{l-1}}\to\mathbb{R}^{N_{l}}\) an affine function. A function \(F:\mathbb{R}^{N_{0}}\to\mathbb{R}^{N_{L}}\) defined as

\[F=w_{L}\circ F_{L-1}\circ\cdots\circ F_{1},\quad\text{with }F_{l}=\varrho\circ w _{l}\quad\text{ for }l=1,\ldots,L-1,\]

is called a _neural network_, with activation function \(\varrho\) applied component-wise. \(L\) denotes the total number of layers, \(N_{1},\ldots,N_{L-1}\) denote the dimensions of the hidden layers and \(N_{0}and\)\(N_{L}\) those of the input and output layers respectively. For each \(l\in\{1,\ldots,L\}\) the affine function \(w_{l}:\mathbb{R}^{N_{l-1}}\to\mathbb{R}^{N_{l}}\) is given as \(w_{l}(\mathbf{x})=\mathrm{A}^{(l)}\mathbf{x}+\mathrm{b}^{(l)}\), for \(\mathbf{x}\in\mathbb{R}^{N_{l-1}}\), with \(\mathrm{A}^{(l)}\in\mathbb{R}^{N_{l}\times N_{l-1}}\) and \(\mathrm{b}^{(l)}\in\mathbb{R}^{N_{l}}\). For any \(i\in\{1,\ldots N_{l}\}\) and \(j\in\{1,\ldots,N_{l-1}\}\) the number \(A^{(l)}_{ij}\) is interpreted as the weight of the edge connecting node \(i\) of layer \(l-1\) to node \(j\) of layer \(l\).

A _random-weight neural network_ (RWNN) is a neural network where the hidden layers are randomly sampled from a given distribution and then fixed; consequently, only the last layer is trained: out of all the parameters \((\mathrm{A}^{(l)},\mathrm{b}^{(l)})_{l=0,\ldots,L}\) of the \(L\)-layered neural network, the parameters \((\mathrm{A}^{(0)},\mathrm{b}^{(0)},\ldots,A^{(L-1)},\mathrm{b}^{(L-1)})\) are randomly sampled and frozen and only \((\mathrm{A}^{(L)},\mathrm{b}^{(L)})\) from the last layer are trained.

The training of such an RWNN can then be simplified into a convex optimisation problem. This makes the training easier to manage and understand both practically and theoretically. However, by only allowing certain parts of the parameters to be trained, the overall capacity and expressivity are possibly reduced. Although it is still unclear if random neural networks still maintain any of the powerful approximation properties of general deep neural networks, these questions have been addressed to some extent in e.g. [33, 53], where learning error bounds for RWNNs have been proved.

Denote now \(\aleph_{\infty}^{\varrho}(d_{0},d_{1})\) the set of random neural networks from \(\mathbb{R}^{d_{0}}\) to \(\mathbb{R}^{d_{1}}\), with activation function \(\varrho\)-and we shall drop the explicit reference to input and output dimensions in the notation whenever they are clear from the context. Moreover, for any \(L,K\in\mathbb{N}\), \(\aleph_{L,K}^{\varrho}\) represents a random neural network with a fixed number of hidden layers \(L\) and fixed input and output dimension \(K\) for each hidden layer. We now give a precise definition of a single layer \(\aleph_{K}^{\varrho}\coloneqq\aleph_{1,K}^{\varrho}\), which we will use for our approximation.

**Definition 2.2** (Single layer RWNN).: Let \((\widetilde{\Omega},\widetilde{\mathcal{F}},\widetilde{\mathbb{P}})\) be a probability space on which the iid random variables \(\mathrm{A}_{k}:\widetilde{\Omega}\to\mathbb{R}^{d}\) and \(b_{k}:\widetilde{\Omega}\to\mathbb{R}\), respectively corresponding to weights and biases, are defined. Let \(\boldsymbol{\phi}=\left\{\phi_{k}\right\}_{k\geq 1}\) denote a sequence of random basis functions, where each \(\phi_{k}:\mathbb{R}^{d}\to\mathbb{R}\) is of the form

\[\phi_{k}(\mathbf{x})\coloneqq\varrho\left(\mathrm{A}_{k}^{\top}\mathbf{x}+b_{ k}\right),\qquad x\in\mathbb{R}^{d},\]

with \(\varrho:\mathbb{R}\to\mathbb{R}\) a Lipschitz continuous activation function. For an output dimension \(m\) and \(K\) hidden units, we define the _reservoir_ or _random basis_ as \(\Phi_{K}\coloneqq\phi_{1:K}=(\phi_{1},\ldots,\phi_{K})\) and the random network \(\aleph_{K}^{\varrho}\) with parameter \(\Theta=\left(\theta_{1},\ldots,\theta_{m}\right)^{\top}\in\mathbb{R}^{m\times K}\) as the map

\[\aleph_{K}^{\varrho}:\mathbf{x}\mapsto\Psi_{K}(\mathbf{x};\Theta)\coloneqq \Theta\Phi_{K}(\mathbf{x}).\]

Thus, for each output dimension \(j\in\{1,\ldots,m\}\), \(\aleph_{K}^{\varrho}\) produces a linear combination of the first \(K\) random basis functions \(\theta_{j}^{\top}\phi_{1:K}\coloneqq\sum_{k=1}^{K}\theta_{j,k}\phi_{k}\).

**Remark 2.3**.: In this paper, we will make use of the more compact vector notation

\[\Phi_{K}:\mathbb{R}^{d}\ni\mathbf{x}\mapsto\boldsymbol{\varrho}(\mathrm{A} \mathbf{x}+\mathrm{b})\in\mathbb{R}^{K},\]
where \(\boldsymbol{\varrho}:\mathbb{R}^{K}\to\mathbb{R}^{K}\) acts component-wise \(\boldsymbol{\varrho}(\mathbf{y})\coloneqq(\varrho(y_{1}),\ldots\varrho(y_{K}))\) and \(\mathrm{A}:\widetilde{\Omega}\to\mathbb{R}^{K\times d}\) and \(\mathrm{b}:\widetilde{\Omega}\to\mathbb{R}^{K}\) are the random matrix and bias respectively.

### Derivatives of ReLu-RWNN

In recent years ReLu neural networks have been predominately used in deep learning, because of their simplicity, efficiency and ability to address the so-called vanishing gradient problem [47]. In many ways, ReLu networks also give a more tractable structure to the optimisation problem compared to their smooth counterparts such as \(\tanh\) and sigmoid. Gonon, Grigoryeva and Ortega [33] derived error bounds to the convergence of a single layer RWNN with ReLu activations. Now, while \(\varsigma(y)\coloneqq\max\{y,0\}\) is performing well numerically, it is, however, not differentiable at zero (see [10] for a short exposition on the chain-rule in ReLu networks). As ReLu-RWNNs will be used in our approach to approximate solutions of partial differential equations, a discussion on its derivatives is in order. To that end we let \(\boldsymbol{\varsigma}(\mathbf{y})\coloneqq(\varsigma(y_{1}),\ldots, \varsigma(y_{K}))\) and \(\boldsymbol{H}(y)=\mathds{1}_{(0,\infty)}(\mathbf{y})\in\mathbb{R}^{K}\) for \(\mathbf{y}\in\mathbb{R}^{K}\), where the indicator function is again applied component-wise.

**Lemma 2.4**.: _For any linear function \(\ell(\mathbf{x})=\mathrm{A}\mathbf{x}+\mathrm{b}\), with \(\mathrm{A}\in\mathbb{R}^{K\times d}\) and \(\mathrm{b}\in\mathbb{R}^{K}\), then_

\[\nabla_{x}(\boldsymbol{\varsigma}\circ\ell)(\mathbf{x})=\mathrm{diag}( \boldsymbol{H}(\mathrm{A}\mathbf{x}+\mathrm{b}))\mathrm{A},\qquad\text{for a.e. }\mathbf{x}\in\mathbb{R}^{d}.\]

Proof.: Let \(\mathcal{A}\coloneqq\big{\{}\mathbf{x}\in\mathbb{R}^{d}:(\boldsymbol{\varsigma }\circ\ell)(\mathbf{x})=0\big{\}}=\big{\{}\mathbf{x}\in\mathbb{R}^{d}:\ell( \mathbf{x})\leq 0\big{\}}\). Then \((\boldsymbol{\varsigma}\circ\ell)(\mathbf{x})=\ell(\mathbf{x})\) for all \(\mathbf{x}\in\mathbb{R}^{d}\backslash\mathcal{A}\). Since \(\ell\) is Lipschitz, differentiability on level sets [23, Section 3.1.2, Corollary I] implies that \(\nabla_{\mathbf{x}}\left(\boldsymbol{\varsigma}\circ\ell\right)(\mathbf{x})= \mathbf{0}\in\mathbb{R}^{d}\) for almost every \(\mathbf{x}\in\mathcal{A}\), and hence

\[\nabla_{\mathbf{x}}(\boldsymbol{\varsigma}\circ\ell)(\mathbf{x})=\mathrm{diag }\left(\mathds{1}_{\{\ell(\mathbf{x})\in\mathbb{R}^{d}\backslash\mathcal{A} \}}\right)\nabla_{\mathbf{x}}\ell(\mathbf{x})=\mathrm{diag}\left(\mathds{1}_ {(0,\infty)}(\ell(\mathbf{x}))\right)\nabla_{\mathbf{x}}\ell(\mathbf{x})= \mathrm{diag}(\boldsymbol{H}(\mathrm{A}\mathbf{x}+\mathrm{b}))\mathrm{A}.\]

Thus by Lemma 2.4, the first derivative of \(\Psi(\cdot;\Theta)\in\aleph_{K}^{\varsigma}\) is equal to

\[\nabla_{\mathbf{x}}\Psi_{K}(\mathbf{x};\Theta)=\Theta\,\mathrm{diag}( \boldsymbol{H}(\mathrm{A}\mathbf{x}+\mathrm{b}))\mathrm{A}\qquad\text{for a.e. }\mathbf{x}\in\mathbb{R}^{d}. \tag{2.1}\]

The above statements hold almost everywhere, it is thus appropriate we introduce a notion of _approximate differentiability_.

**Definition 2.5** (Approximate limit, [23, Section 1.7.2]).: Consider a Lebesgue-measurable set \(E\subset\mathbb{R}^{d}\), a measurable function \(f:E\to\mathbb{R}^{m}\) and a point \(\mathbf{x}_{0}\in E\). We say \(l\in\mathbb{R}^{m}\) is the approximate limit of \(f\) at \(\mathbf{x}_{0}\), and write \(\mathrm{ap}\lim_{\mathbf{x}\to x_{0}}f(\mathbf{x})=l\), if for each \(\varepsilon>0\),

\[\lim_{r\downarrow 0}\frac{\lambda\left(\mathcal{B}_{r}(\mathbf{x}_{0})\cap\{ \mathbf{x}\in E:\;|f(\mathbf{x})-l|\geq\varepsilon\}\right)}{\lambda(\mathcal{B }_{r}(\mathbf{x}_{0}))}=0,\]

with \(\lambda\) the Lebesgue measure and \(\mathcal{B}_{r}(\mathbf{x}_{0})\) the closed ball with radius \(r>0\) and center \(\mathbf{x}_{0}\).

**Definition 2.6** (Approximate differentiability, [23, Section 6.1.3]).: Consider a measurable set \(E\subset\mathbb{R}^{d}\), a measurable map \(f:E\to\mathbb{R}^{m}\) and a point \(\mathbf{x}_{0}\in E\). The map \(f\) is approximately differentiable at \(\mathbf{x}_{0}\) if there exists a linear map \(\mathrm{D}_{\mathbf{x}}:\mathbb{R}^{d}\to\mathbb{R}^{m}\) such that

\[\mathrm{ap}\lim_{\mathbf{x}\to x_{0}}\frac{f(\mathbf{x})-f(\mathbf{x}_{0})- \mathrm{D}_{\mathbf{x}}(\mathbf{x}-\mathbf{x}_{0})}{|\mathbf{x}-\mathbf{x}_{0}| }=0.\]

Then \(\mathrm{D}_{\mathbf{x}}\) is called the approximate differential of \(f\) at \(\mathbf{x}_{0}\).

**Remark 2.7**.: The usual rules from classical derivatives, such as the uniqueness of the differential, and differentiability of sums, products and quotients, apply to approximately differentiable functions. Moreover, the chain rule applies to compositions \(\varphi\circ f\) when \(f\) is approximately differentiable at \(\mathbf{x}_{0}\) and \(\varphi\) is classically differentiable at \(f(\mathbf{x}_{0})\).
**Remark 2.8** ([23, Theorem 4, Section 6.1.3]).: For \(f\in\mathcal{W}^{1,p}_{\mathrm{loc}}\left(\mathbb{R}^{d}\right)\) and \(1\leq p\leq\infty\), \(f\) is approximately differentiable almost everywhere and its approximate derivative equals its weak derivative almost everywhere. We will thus use the operator \(\mathrm{D}_{\mathbf{x}}\) to denote the weak derivative and approximate derivative interchangeably, to distinguish them from the classical derivative denoted by \(\nabla\).

**Lemma 2.9**.: _Let \(E\subset\mathbb{R}^{d}\) be a measurable set with finite measure, \(X:\Omega\to E\) a continuous random variable on some probability space \((\Omega,\mathcal{F},\mathbb{P})\), \(\varphi\in\mathcal{C}^{1}(\mathbb{R}^{m})\), \(\Phi_{\mathrm{ap}}:E\to\mathbb{R}^{m}\) an approximately differentiable function, and \(\Phi\in\mathcal{C}^{1}(\mathbb{R}^{d};\mathbb{R}^{m})\) its continuously differentiable extension to \(\mathbb{R}^{d}\). Then \(\mathbb{E}[\varphi(\mathrm{D}_{x}\Phi_{\mathrm{ap}}(X))]=\mathbb{E}[\varphi( \nabla_{x}\Phi(X))]\)._

Proof.: By [24, Theorem 3.1.6] a function \(\Phi_{\mathrm{ap}}:E\to\mathbb{R}^{m}\) is approximately differentiable almost everywhere if for every \(\varepsilon>0\) there is a compact set \(F\subset E\) such that the Lebesgue measure \(\lambda(E\backslash F)<\varepsilon\) and \(\left.\Phi_{\mathrm{ap}}\right|_{F}\) is \(\mathcal{C}^{1}\). Since \(\varphi\) is everywhere differentiable, it maps null-sets to null-sets [58, Lemma 7.25]. The claim follows since \(\mathbb{P}\) is absolutely continuous with respect to the Lebesgue measure \(\lambda\), \(X\) being a continuous random variable.

**Corollary 2.10**.: _Let \(E\subset\mathbb{R}^{d}\) be a measurable set with finite measure, \(X:\Omega\to E\) a continuous random variable on some probability space \((\Omega,\mathcal{F},\mathbb{P})\), \(\varphi\in\mathcal{C}^{1}(\mathbb{R}^{m})\), \(\Phi:E\to\mathbb{R}^{m}\) an approximately differentiable function, and \(\Psi\in\mathcal{W}^{1,p}(E,\mathbb{R}^{m})\) for \(p\geq 1\) such that \(\Phi=\Psi\) almost everywhere. Then \(\mathbb{E}[\varphi(\mathrm{D}_{\mathbf{x}}\Phi(X))]=\mathbb{E}[\varphi( \mathrm{D}_{\mathbf{x}}\Psi(X))]\)._

Proof.: This is a direct consequence of Lemma 2.9, after noting that the two notions of derivatives are the same on \(\mathcal{W}^{1,p}(E,\mathbb{R}^{m})\) (see Remark 2.8). 

From a practical perspective, the second-order derivative of the network with respect to the input will be zero for all intents and purposes. However, as will become apparent in Lemma 5.6, we need to investigate it further, in particular the measure zero set of points where ReLu-RWNN is not differentiable. Rewriting the diagonal operator in terms of the natural basis \(\{e_{i}\}\) and evaluating the function \(\boldsymbol{H}\) component-wise yields

\[\nabla_{\mathbf{x}}\Psi_{K}(\mathbf{x};\Theta)=\Theta\left(\sum_{j=1}^{K}e_{j }e_{j}^{\top}H\left(e_{j}^{\top}\mathrm{A}\mathbf{x}+b_{j}\right)\right) \mathrm{A}.\]

The \(i\)-th component of the second derivative is thus

\[\left[\nabla_{\mathbf{x}}^{2}\Psi_{K}(\mathbf{x};\Theta)\right]_{i}=\Theta \left(\sum_{j=1}^{K}e_{j}e_{j}^{\top}a_{ji}H^{\prime}\left(e_{j}^{\top} \mathrm{A}\mathbf{x}+b_{j}\right)\right)\mathrm{A}=\Theta\operatorname{diag} \left(a_{i}\right)\operatorname{diag}\left(\boldsymbol{H}^{\prime}(\mathrm{A} \mathbf{x}+\mathrm{b})\right)\mathrm{A},\]

where \(a_{i}\) denotes the \(i\)-th column of the matrix \(\mathrm{A}\). Next, we let

\[\delta_{0}^{\varepsilon}(\mathbf{x})\coloneqq\frac{H(\mathbf{x})-H(\mathbf{x }-\varepsilon)}{\varepsilon}\]

for \(\mathbf{x}\in\mathbb{R}\) and define the left derivative of \(H\) as \(H^{\prime}=\lim_{\varepsilon\downarrow 0}\delta_{0}^{\varepsilon}=\delta_{0}\) in the distributional sense. This finally gives the second derivative of the network:

\[\left[\nabla_{\mathbf{x}}^{2}\Psi_{K}(\mathbf{x};\Theta)\right]_{i}=\Theta \operatorname{diag}\left(a_{i}\right)\operatorname{diag}\left(\boldsymbol{ \delta}_{0}(\mathrm{A}\mathbf{x}+\mathrm{b})\right)\mathrm{A},\]

where \(\boldsymbol{\delta}_{0}\) denotes the vector function applying \(\delta_{0}\) component-wise.
### Randomised least squares (RLS)

Let \(Y\in\mathbb{R}^{d}\) and \(X\in\mathbb{R}^{k}\) random variables defined on some probability space \((\Omega,\mathcal{F},\mathbb{P})\) and let \(\boldsymbol{\beta}\in\mathbb{R}^{d\times k}\) be a deterministic coefficient matrix. If the loss function is the mean square error (MSE), we can derive the estimator for the randomised least square:

\[\nabla_{\boldsymbol{\beta}}\mathbb{E}\left[\|Y-\boldsymbol{\beta}X \|^{2}\right] =\nabla_{\boldsymbol{\beta}}\mathbb{E}[(Y-\boldsymbol{\beta}X)^{ \top}(Y-\boldsymbol{\beta}X)]\] \[=\mathbb{E}\left[\nabla_{\boldsymbol{\beta}}(Y^{\top}Y-Y^{\top} \boldsymbol{\beta}X-X^{\top}\boldsymbol{\beta}^{\top}Y+X^{\top}\boldsymbol{ \beta}^{\top}\boldsymbol{\beta}X)\right]\] \[=\mathbb{E}\left[2\boldsymbol{\beta}XX^{\top}-2YX^{\top}\right],\]

which gives the minimiser1\(\boldsymbol{\beta}=\mathbb{E}\left[YX^{\top}\right]\mathbb{E}\left[XX^{\top} \right]^{-1}\), and its estimator

Footnote 1: The matrix \(\mathbb{E}[XX^{\top}]\) may not be invertible, but its generalised Moore-Penrose inverse always exists.

\[\widehat{\boldsymbol{\beta}}\coloneqq\left(\sum_{j=1}^{n}Y_{j}X_{j}^{\top} \right)\left(\sum_{j=1}^{n}X_{j}X_{j}^{\top}\right)^{-1}.\]

Depending on the realisation of the reservoir of the RWNN, the covariates of \(X\) may be collinear, so that \(X\) is close to rank deficient. A standard remedy is to use the Ridge regularised version [38] of the estimator

\[\widehat{\boldsymbol{\beta}}_{R}=\left(\sum_{j=1}^{n}Y_{j}X_{j}^{\top}\right) \left(\sum_{j=1}^{n}X_{j}X_{j}^{\top}+\lambda I\right)^{-1},\quad\text{for }\lambda>0,\]

which results in a superior, more robust performance in our experiments.

**Remark 2.11**.: The above derivation holds true for the approximate derivative \(\mathrm{D}_{x}\) as well because all operations above hold for approximately differentiable functions (Remark 2.7).

**Remark 2.12**.: One of the advantages of RLS estimators is the ability to vectorise them over the sampling dimension (i.e. over \(n\) samples) in order to use tensor functionalities of packages such as NumPy and PyTorch to efficiently evaluate the sum of outer products using the einsum function. Details are provided in the code available at ZuricZ/RWNN_PDE_solver.

## 3. The Markovian case

Let the process \(\mathbf{X}\) of the traded and non-traded components of the underlying under the risk-neutral measure \(\mathbb{Q}\) be given by the following \(d\)-dimensional SDE:

\[\mathbf{X}_{s}^{t,\mathbf{x}}=\mathbf{x}+\int_{t}^{s}\mu(r,\mathbf{X}_{r}^{t, \mathbf{x}})\mathrm{d}r+\int_{t}^{s}\Sigma\left(r,\mathbf{X}_{r}^{t,\mathbf{x }}\right)\mathrm{d}W_{r}, \tag{3.1}\]

where \(\mu:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d}\) and \(\Sigma:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}^{d\times d}\) adhere to Assumption A.1, and \(W\) is a standard \(d\)-dimensional Brownian motion on the probability space \((\Omega,\mathcal{F},\mathbb{Q})\) equipped with the natural filtration \(\mathbb{F}=\{\mathcal{F}_{t}\}_{0\leq t\leq T}\) of \(W\). By the Feynman-Kac formula, options whose discounted expected payoff under \(\mathbb{Q}\) can be represented as

\[u(t,\mathbf{x})=\mathbb{E}\left[\int_{t}^{T}\mathrm{e}^{-r(s-t)}f\left(s, \mathbf{X}_{s}^{t,\mathbf{x}}\right)\mathrm{d}s+\mathrm{e}^{-r(T-t)}g\left( \mathbf{X}_{T}^{t,\mathbf{x}}\right)\right]\quad\text{ for all }(t,\mathbf{x})\in[0,T]\times \mathcal{A},\]
for \(\mathcal{A}\subset\mathbb{R}^{d}\) with interest rate \(r\geq 0\) and continuous functions \(f:[0,T]\times\mathbb{R}^{d}\to\mathbb{R}\) and \(g:\mathbb{R}^{d}\to\mathbb{R}\) can be viewed as solutions to the Cauchy linear parabolic PDE

\[\left\{\begin{array}{rl}\partial_{t}u+\mathcal{L}u+f-ru=0,&\text{ on }[0,T)\times\mathcal{A},\\ u(T,\cdot)=g,&\text{ on }\mathcal{A},\end{array}\right.\]

where

\[\mathcal{L}u\coloneqq\frac{1}{2}\operatorname{Tr}\left(\Sigma\Sigma^{\top} \nabla_{\mathbf{x}}^{2}u\right)+(\nabla_{\mathbf{x}}u)\mu,\qquad\text{on }[0,T)\times\mathcal{A}, \tag{3.2}\]

is the infinitesimal generator associated with diffusion (3.1). In this Markovian setting, we thus adopt a set-up similar to [42] and consider a slightly more general PDE

\[\left\{\begin{array}{rl}\partial_{t}u(t,\mathbf{x})+\mathcal{L}u(t,\mathbf{ x})+f\Big{(}t,\mathbf{x},u(t,\mathbf{x}),\nabla_{\mathbf{x}}u(t,\mathbf{x}) \cdot\Sigma(t,\mathbf{x})\Big{)}&=0,\qquad\text{on }[0,T)\times\mathcal{A},\\ u(T,\mathbf{x})&=g(\mathbf{x}),&\text{on }\mathcal{A},\end{array}\right. \tag{3.3}\]

with \(f:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\to\mathbb{R}\) such that Assumption A.1 is satisfied, which guarantees the existence and uniqueness of the solution to the corresponding backward stochastic differential equation (BSDE) [54, Section 4]. The corresponding second-order generator is again given by (3.2).

The following assumption is only required to cast the problem into a regression. Otherwise, the optimisation (below) can still be solved using other methods, such as stochastic gradient descent. Another solution would be to use the so-called splitting method to linearise the PDE (as in [6] and the references therein for example).

**Assumption 3.1**.: The function \(f:[0,T]\times\mathbb{R}^{d}\times\mathbb{R}\times\mathbb{R}^{d}\times\mathbb{R }^{d}\to\mathbb{R}\) has an affine structure in \(\mathbf{y}\in\mathbb{R}^{m}\) and in \(\mathbf{z}^{1},\mathbf{z}^{2}\in\mathbb{R}^{d}\):

\[f\left(t,\mathbf{x},\mathbf{y},\mathbf{z}^{1},\mathbf{z}^{2}\right)=a(t, \mathbf{x})\mathbf{y}+b(t,\mathbf{x})\mathbf{z}^{1}+c(t,\mathbf{x})\mathbf{z} ^{2}+\widetilde{f}(t,\mathbf{x}),\]

for some real-valued functions \(a,b,c,\widetilde{f}\) on \([0,T]\times\mathbb{R}^{d}\) that map to conformable dimensions.

### Random weighted neural network scheme

The first step in so-called deep BSDE schemes [19, 36, 42] is to establish the BSDE associated with the PDE (3.3) and the process (3.1) through the nonlinear Feynman-Kac formula. By [54] there exist \(\mathbb{F}\)-adapted processes \((Y,Z)\), which are unique solutions to the BSDE

\[Y_{t}=g\left(\mathbf{X}_{T}\right)+\int_{t}^{T}f\left(s,\mathbf{X}_{s},Y_{s},Z _{s}\right)\mathrm{d}s-\int_{t}^{T}Z_{s}\mathrm{d}W_{s},\qquad\text{for any }t\in[0,T], \tag{3.4}\]

and which are connected to the PDE (3.3) via

\[Y_{t}=u(t,\mathbf{X}_{t})\qquad\text{and}\qquad Z_{t}=\nabla_{x}u(t,\mathbf{ X}_{t})\cdot\Sigma(t,\mathbf{X}_{t}).\]

with terminal condition \(u(T,\cdot)=g\). Next, the BSDE (3.4) is rewritten in forward form

\[Y_{t}=Y_{0}-\int_{0}^{t}f\left(s,\mathbf{X}_{s},Y_{s},Z_{s}\right)\mathrm{d}s +\int_{0}^{t}Z_{s}\mathrm{d}W_{s},\qquad\text{for any }t\in[0,T],\]

and both processes are discretised according to the Euler-Maruyama scheme. To this end let \(\pi\coloneqq\{0=t_{0}<t_{1}<\ldots<t_{N}=T\}\) be a partition of the time interval \([0,T]\) with modulus \(|\pi|=\max_{i=\{0,1,\ldots,N-1\}}\delta_{i}\) and \(\delta_{i}\coloneqq t_{i+1}-t_{i}\). Then the scheme is given by

\[\left\{\begin{array}{rl}\mathbf{X}_{t_{i+1}}&=\mathbf{X}_{t_{i}}+\mu(t_{i}, \mathbf{X}_{t_{i}})\delta_{i}+\Sigma(t_{i},\mathbf{X}_{t_{i}})\Delta_{i}^{W}, \\ Y_{t_{i+1}}&=Y_{t_{i}}-f\left(t_{i},\mathbf{X}_{t_{i}},Y_{t_{i}},Z_{t_{i}} \right)\delta_{i}+Z_{t_{i}}\Delta_{i}^{W},\end{array}\right. \tag{3.5}\]
where naturally \(\Delta_{i}^{W}\coloneqq W_{t_{i+1}}-W_{t_{i}}\). Then for all \(i\in\{N-1,\ldots,0\}\) we approximate \(u(t_{i},\cdot)\) with \(\mathfrak{U}_{i}(\cdot;\Theta^{i})\in\aleph_{K}^{\rho}\) and \(Z_{t_{i}}\) as

\[\begin{array}{llll}u(t_{i},\mathbf{X}_{t_{i}})=Y_{t_{i}}&\approx\mathfrak{U} _{i}(\mathbf{X}_{t_{i}};\Theta^{i})&\coloneqq\Theta^{i}\Phi_{K}^{i}(\mathbf{X} _{t_{i}}),\\ Z_{t_{i}}&\approx\mathfrak{J}_{i}(\mathbf{X}_{t_{i}})&\coloneqq\mathrm{D}_{ \mathbf{x}}\mathfrak{U}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})\cdot\Sigma(t_{i}, \mathbf{X}_{t_{i}})=\Theta^{i}\mathrm{D}_{\mathbf{x}}\Phi_{K}^{i}(\mathbf{X}_ {t_{i}})\cdot\Sigma(t_{i},\mathbf{X}_{t_{i}}).\end{array}\]

Recall that the derivative \(\mathfrak{J}_{i}(\mathbf{X}_{t_{i}})\) is the approximate derivative from Definition 2.6. The following formulation of the loss function \(\ell\) using the approximate derivative is sensible by Lemma 2.9: notice that for the optimal parameter \(\Theta^{i+1,*}\) in step \((i+1)\), the optimal approximation \(\widehat{\mathfrak{U}}_{i+1}(\mathbf{X}_{t_{i+1}})\coloneqq\mathfrak{U}_{i+1} (\mathbf{X}_{t_{i+1}};\Theta^{i+1,*})\) does not depend on \(\Theta^{i}\), hence under Assumption 3.1 with \(c=0\) the loss function at the \(i\)-th discretisation step reads

\[\begin{array}{l}\ell(\Theta^{i})\coloneqq\mathbb{E}^{\Phi}\left[\left\| \widehat{\mathfrak{U}}_{i+1}(\mathbf{X}_{t_{i+1}})-\left[\mathfrak{U}_{i}( \mathbf{X}_{t_{i}};\Theta^{i})-f\Big{(}t_{i},\mathbf{X}_{t_{i}},\mathfrak{U}_{ i}(\mathbf{X}_{t_{i}};\Theta^{i}),\mathfrak{J}_{i}(\mathbf{X}_{t_{i}};\Theta^{i}) \Big{)}\delta_{i}+\mathfrak{J}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})\Delta_{i}^{W }\right]\right\|^{2}\right]\\ \phantom{\ell(\Theta^{i})}=\mathbb{E}^{\Phi}\left[\left\|\widehat{\mathfrak{U }}_{i+1}(\mathbf{X}_{t_{i+1}})-\left[(\mathfrak{U}_{i}(\mathbf{X}_{t_{i}}; \Theta^{i})-\left(a_{i}\mathfrak{U}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})+b_{i} \mathfrak{J}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})+\widetilde{f}_{i}\right)\delta _{i}+\mathfrak{J}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})\Delta_{i}^{W}\right] \right\|^{2}\right]\\ \phantom{\ell(\Theta^{i})}=\mathbb{E}^{\Phi}\left[\left\|\widehat{\mathfrak{U }}_{i+1}(\mathbf{X}_{t_{i+1}})+\widetilde{f}_{i}\delta_{i}-\Theta^{i}\Big{\{}( 1-a_{i}\delta_{i})\Phi_{K}^{i}(\mathbf{X}_{t_{i}})+\mathrm{D}_{x}\Phi_{K}^{i}( \mathbf{X}_{t_{i}})\Sigma_{i}\left(b_{i}\delta_{i}+\Delta_{i}^{W}\right) \Big{\}}\right\|^{2}\right]\\ \phantom{\ell(\Theta^{i})}=\mathbb{E}^{\Phi}\left[\left\|\mathbf{Y}^{i}-\Theta ^{i}\mathrm{X}^{i}\right\|^{2}\right]\end{array}\]

where \(p_{i}\coloneqq p(t_{i},\mathbf{X}_{t_{i}})\) for \(p\in\{a,b,\widetilde{f},\Sigma\}\), and the expectation \(\mathbb{E}^{\Phi}\) is of course conditional on the realisation of the random basis \(\Phi_{K}^{i}\), i.e., conditional on the random weights and biases of the RWNN. Furthermore, we used the notations

\[\mathrm{Y}^{i}\coloneqq\widehat{\mathfrak{U}}_{i+1}(\mathbf{X}_{t_{i+1}})+ \widetilde{f}_{i}\delta_{i}\qquad\text{and}\qquad\mathrm{X}^{i}\coloneqq(1-a_ {i}\delta_{i})\Phi_{K}(\mathbf{X}_{t_{i}})+\mathrm{D}_{\mathbf{x}}\Phi_{K}( \mathbf{X}_{t_{i}})\cdot\Sigma_{i}\left(b_{i}\delta_{i}+\Delta_{i}^{W}\right).\]

The problem can now be solved via least squares from Section 2.2, yielding the estimator

\[\Theta^{i,*}=\mathbb{E}^{\Phi}\left[\mathrm{Y}^{i}\mathrm{X}^{i\top}\right] \mathbb{E}^{\Phi}\left[\mathrm{X}^{i}\mathrm{X}^{i\top}\right]^{-1}.\]

### Algorithm

We now summarise the algorithmic procedure of our RWNN scheme. See how the algorithm resembles the Least-Square Monte-Carlo method of [49] after considering sample estimator version of RLS from Section 2.2:
```
Input: \(\pi=\{0=t_{0}<t_{1}<\ldots<t_{N}=T\}\) time grid Initialisation: \(\Phi_{K}^{i}\) reservoirs with \(K\in\mathbb{N}\) hidden nodes with weights and biases distributed as \(\mathcal{U}_{[-R,R]}\) with \(R>0\) for all \(i\in\{0,\ldots,N-1\}\)do:  Generate \(n\in\mathbb{N}\) paths of \(\{\mathbf{X}_{t_{i}}^{\pi,j}\}_{i=0}^{N}\) for \(j\in\{1,\ldots,n\}\) with the Euler-Maruyama scheme (3.5)  Set \(\widehat{\mathfrak{U}}_{N}(\mathbf{x})=g(\mathbf{x})\) for all \(\mathbf{x}\in\mathbb{R}^{d}\) for\(i\in\{N-1,\ldots,0\}\)do  Approximate \(u(t_{i},\cdot)\) with \(\mathfrak{U}(\cdot;\Theta^{i})\in\aleph_{K}^{\varsigma}\) based on reservoir \(\Phi_{K}^{i}\)  Evaluate the derivative of \(\mathfrak{U}(\cdot;\Theta^{i})\) according to (2.1)  Solve the regression problem (possibly using the Ridge estimator, see Section 2.2) \[\Theta^{i,*}=\operatorname*{arg\,min}_{\Theta^{i}}\ell(\Theta^{i})= \operatorname*{arg\,min}_{\Theta^{i}}\mathbb{E}^{\Phi,n}\left[\left\|\mathrm{ Y}^{i}-\Theta^{i}\mathrm{X}^{i}\right\|^{2}\right]\] where \[\mathrm{Y}^{i} \coloneqq\mathfrak{U}_{i}(\mathbf{X}_{t_{i}};\Theta^{i})+ \widetilde{f}\delta_{i}\] \[\mathrm{X}^{i} \coloneqq(1-a)\Phi_{K}^{i}(\mathbf{X}_{t_{i}})+(\nabla_{x} \Phi_{K}^{i}(\mathbf{X}_{t_{i}}))\Sigma_{i}\left(b\delta_{i}+\Delta_{i}^{W} \right).\] and \(\mathbb{E}^{\Phi,n}\) is evaluated over the empirical measure of \(\{\mathbf{X}_{t_{i}}^{\pi,j}\}_{i=0}^{N}\) for \(j\in\{1,\ldots,n\}\)  Update \(\widehat{\mathfrak{U}}_{i}=\mathfrak{U}_{i}\left(\cdot,\Theta^{i,*}\right)\) endfor return\(\mathfrak{U}=\{\mathfrak{U}(\cdot;\Theta^{i,*})\}_{i=0}^{N}\).
```

**Algorithm 1** RWNN scheme

## 4. The non-Markovian case

We now consider a stochastic volatility model under a risk-neutral measure so that \(\mathbf{X}=(X,V)\), where the dynamics of log-price process \(X\) are given by

\[\mathrm{d}X_{s}^{t,x}=\bigg{(}r-\frac{V_{s}}{2}\bigg{)}\,\mathrm{d}s+\sqrt{V_{s }}\Big{(}\rho_{1}\mathrm{d}W_{s}^{1}+\rho_{2}\mathrm{d}W_{s}^{2}\Big{)},\qquad 0 \leq t\leq s\leq T, \tag{4.1}\]

starting from \(X_{t}^{t,x}=x\in\mathbb{R}\), with interest rate \(r\in\mathbb{R}\), correlation \(\rho_{1}\in[-1,1]\), and denote \(\rho_{2}\coloneqq\sqrt{1-\rho_{1}^{2}}\), and \(W^{1},W^{2}\) are two independent Brownian motions. We allow for a general variance process process \(V\), satisfying the following:

**Assumption 4.1**.: The process \(V\) has continuous trajectories, is non-negative, adapted to the natural filtration of \(W^{1}\) and integrable, i.e. \(\mathbb{E}\left[\int_{0}^{t}V_{s}\mathrm{d}s\right]<\infty\), for any \(t\geq 0\).

By no-arbitrage, the fair price of a European option with payoff \(h:\mathbb{R}^{+}\to\mathbb{R}^{+}\) reads

\[u(t,x)\coloneqq\mathbb{E}\left[\mathrm{e}^{-r(T-t)}h\left(\mathrm{e}^{X_{T}^{t,x}+rT}\right)\middle|\mathcal{F}_{t}\right],\quad\text{for all }(t,x)\in[0,T]\times\mathbb{R},\]

subject to (4.1). Since \(\mathbf{X}\) is not Markovian, one cannot characterise the value function \(u(t,x)\) via a deterministic PDE. Bayer, Qiu and Yao [5] proved that \(u\) can be viewed as a random field which, together with another random field \(\psi\), satisfies the backward
stochastic partial differential equation (BSPDE)

\[-{\rm d}u(t,x)=\left[\frac{V_{t}}{2}\partial_{x}^{2}u(t,x)+\rho\sqrt{V_{t}} \partial_{x}\psi(t,x)-\frac{V_{t}}{2}\partial_{x}u(t,x)-ru(t,x)\right]{\rm d}t- \psi(t,x){\rm d}W_{t}^{1}, \tag{4.2}\]

in a distributional sense, for \((t,x)\in[0,T]\times\mathbb{R}\), with boundary condition \(u(T,x)=h\left({\rm e}^{x+rT}\right)\). where the variance process \(\{V_{t}\}_{t\geq 0}\) is defined exogenously under Assumption 4.1. We in fact consider the slightly more general BSPDEs

\[\left\{\begin{array}{rl}-{\rm d}u(t,x)&=\left\{\frac{V_{t}}{2} \mathrm{D}_{x}^{2}u(t,x)+\rho\sqrt{V_{t}}\mathrm{D}_{x}\psi(t,x)-\frac{V_{t}}{2 }\mathrm{D}_{x}u(t,x)\right.\\ &\left.+f\left(t,{\rm e}^{x},u(t,x),\rho_{2}\sqrt{V_{t}}\mathrm{D}_{x}u(t,x), \psi(t,x)+\rho_{1}\sqrt{V_{t}}\mathrm{D}_{x}u(t,x)\right)\,\right\}{\rm d}t\\ &-\psi(t,x){\rm d}W_{t}^{1},\quad(t,x)\in[0,T)\times\mathbb{R},\\ u(T,x)&=g\left({\rm e}^{x}\right),\quad x\in\mathbb{R},\end{array}\right. \tag{4.3}\]

for \(g:\mathbb{R}^{+}\to\mathbb{R}\). We present exact conditions on \(f\) and \(g\) ensuring well-posedness in Assumption A.2, but additionally require the existence of a weak-Sobolev solution (see Assumption 5.1 for more details). Note that (4.2) is just a particular case of the general BSPDE (4.3) for the choice \(f(t,x,y,z,\widetilde{z})\equiv-ry\) and \(g({\rm e}^{x})\equiv h({\rm e}^{x+rT})\). Again, this general form is shown to be well-posed in the distributional sense under Assumption A.2 (borrowed from [5]). By [14] the corresponding BSDE is then, for \(0\leq t\leq s\leq T\),

\[\left\{\begin{aligned} -{\rm d}Y_{s}^{t,x}&=f\left(s,{ \rm e}^{X_{s}^{t,x}},Y_{s}^{t,x},{Z_{s}^{1}}^{t,x},{Z_{s}^{2}}^{t,x}\right){ \rm d}s-{Z_{s}^{1}}^{t,x}{\rm d}W_{s}^{1}-{Z_{s}^{2}}^{t,x}{\rm d}W_{s}^{2},\\ Y_{T}^{t,x}&=g\left({\rm e}^{X_{T}^{t,x}}\right), \end{aligned}\right. \tag{4.4}\]

where \((Y_{s}^{t,x},{Z_{s}^{1}}^{t,x},{Z_{s}^{2}}^{t,x})\) is defined as the solution to BSDE (4.4) in the weak sense.

### Random neural network scheme

Let the quadruple \(\left(X_{s},Y_{s},Z_{s}^{1},Z_{s}^{2}\right)\) be the solution to the forward BSDE (FBSDE)

\[\left\{\begin{aligned} -{\rm d}Y_{s}&=f\left(s,{ \rm e}^{X_{s}},Y_{s},Z_{s}^{1},Z_{s}^{2}\right){\rm d}s-{Z_{s}^{1}}{\rm d}W_{s }^{1}-{Z_{s}^{2}}{\rm d}W_{s}^{2},\\ {\rm d}X_{s}&=-\frac{V_{s}}{2}{\rm d}s+\sqrt{V_{s}} \left(\rho_{1}{\rm d}W_{s}^{1}+\rho_{2}{\rm d}W_{s}^{2}\right),\\ V_{s}&=\xi_{s}\mathcal{E}\left(\eta\widehat{W}_{s} \right),\quad\text{ with }\quad\widehat{W}_{s}=\int_{0}^{s}\mathcal{K}(s,r){\rm d}W_{r}^{1}, \quad,\end{aligned}\right. \tag{4.5}\]

for \(s\in[0,T]\), with terminal condition \(Y_{T}=g\left({\rm e}^{X_{T}}\right)\), initial condition \(X_{0}=x\) and \(\mathcal{K}\) a locally square-integrable kernel. For notational convenience below, we use \(\rho_{2}\coloneqq\sqrt{1-\rho_{1}^{2}}\), with \(\rho_{1}\in[-1,1]\). Here \(\mathcal{E}(\cdot)\) denotes the Wick stochastic exponential and is defined as \(\mathcal{E}(\zeta)\coloneqq\exp\left(\zeta-\frac{1}{2}\mathbb{E}[|\zeta|^{2}]\right)\) for a zero-mean Gaussian variable \(\zeta\). Then by [5, Theorem 2.4],

\[Y_{t} =u\left(t,X_{t}\right), \text{ for }t\in[0,T],\] \[Z_{t}^{1} =\psi\left(t,X_{t}\right)+\rho_{1}\sqrt{V_{t}}\mathrm{D}_{x}u \left(t,X_{t}\right), \text{ for }t\in[0,T),\] \[Z_{t}^{2} =\rho_{2}\sqrt{V_{t}}\mathrm{D}_{x}u\left(t,X_{t}\right), \text{ for }t\in[0,T),\]

where \((u,\psi)\) is the unique weak solution to (4.4). Accordingly, the forward equation reads

\[Y_{t}=Y_{0}-\int_{0}^{t}f\left(s,{\rm e}^{X_{s}},Y_{s},Z_{s}^{1},Z_{s}^{2} \right){\rm d}s+\int_{0}^{t}Z_{s}^{1}{\rm d}W_{s}^{1}+\int_{0}^{t}Z_{s}^{2}{ \rm d}W_{s}^{2},\qquad\text{ for }t\in[0,T].\]
By simulating \((W^{1},W^{2},V)\), the forward process \(X\) may be approximated by an Euler scheme-with the same notations as in the Markovian case-and the forward representation above yields the approximation

\[u\left(t_{i+1},X_{t_{i+1}}\right)\approx u(t_{i},X_{t_{i}})-f\left(t_{i},\mathrm{ e}^{X_{t_{i}}},u\left(t_{i},X_{t_{i}}\right),Z^{1}_{t_{i}},Z^{2}_{t_{i}} \right)\delta_{i}+Z^{1}_{t_{i}}\Delta_{i}^{W^{1}}+Z^{2}_{t_{i}}\Delta_{i}^{W^{ 2}},\]

with

\[Z^{1}_{t_{i}}=\rho_{1}\sqrt{V_{t_{i}}}\mathrm{D}_{x}u\left(t_{i},X_{t_{i}} \right)+\psi\left(t_{i},X_{t_{i}}\right)\qquad\text{and}\qquad Z^{2}_{t_{i}}= \rho_{2}\sqrt{V_{t_{i}}}\mathrm{D}_{x}u\left(t_{i},X_{t_{i}}\right).\]

By Lemma B.2 we can, for each time step \(i\in\{0,\ldots,N-1\}\), approximate the solutions \(u(t_{i},\cdot)\) and \(\psi(t_{i},\cdot)\) by two separate networks \(\mathfrak{U}_{i}\) and \(\Psi_{i}\) in \(\aleph_{K}^{\varsigma}\):

\[Y_{t_{i}} \approx\mathfrak{U}_{i}(X_{t_{i}};\Theta^{i})\qquad=\Theta^{i} \Phi^{\Theta,i}_{K}(X_{t_{i}}),\] \[Z^{1}_{t_{i}} \approx\mathfrak{Z}^{1}_{i}(X_{t_{i}};\Theta^{i},\Xi^{i})=\Theta ^{i}\left(\mathrm{D}_{x}\Phi^{\Theta,i}_{K}(X_{t_{i}})\right)\rho_{1}\sqrt{V_ {t_{i}}}+\Xi^{i}\Phi^{\Xi,i}_{K}(X_{t_{i}}),\] \[Z^{2}_{t_{i}} \approx\mathfrak{Z}^{2}_{i}(X_{t_{i}};\Theta^{i},\Xi^{i})=\Theta ^{i}\left(\mathrm{D}_{x}\Phi^{\Theta,i}_{K}(X_{t_{i}})\right)\rho_{2}\sqrt{V_ {t_{i}}}.\]

Here \(\Phi^{\Xi}_{K}\) and \(\Phi^{\Theta}_{K}\) are realisations of random bases (reservoirs) of the RWNNs with respective parameters \(\Xi\) and \(\Theta\). The next part relies on Assumption 3.1, namely

\[f\left(t_{i},\mathrm{e}^{X_{t_{i}}},Y_{t_{i}},Z^{1}_{t_{i}},Z^{2}_{t_{i}} \right)=a(t_{i},X_{t_{i}})Y_{t_{i}}+b(t_{i},X_{t_{i}})Z^{1}_{t_{i}}+c(t_{i},X_ {t_{i}})Z^{2}_{t_{i}}+\widetilde{f}(t_{i},X_{t_{i}}),\]

for some functions \(a,b,c,\widetilde{f}\) mapping to \(\mathbb{R}\), so that, as in the Markovian case, the minimisation of the expected quadratic loss at every time step \(i\in\{N-1,\ldots,0\}\) reads

\[\ell(\Theta^{i},\Xi^{i})\] \[\coloneqq\mathbb{E}^{\Phi}\Bigg{[}\Bigg{|}\widehat{\mathfrak{U}}_ {i+1}(X_{t_{i+1}})-\bigg{\{}\mathfrak{U}_{i}(X_{t_{i}};\Theta^{i})-f\Big{(}a \mathfrak{U}_{i}(X_{t_{i}};\Theta^{i})+b\mathfrak{Z}^{1}_{i}(X_{t_{i}};\Theta^ {i},\Xi^{i})+c\mathfrak{Z}^{2}_{i}(X_{t_{i}};\Theta^{i},\Xi^{i})+\widetilde{f }\Big{)}\,\delta_{i}\] \[\qquad\qquad\qquad\qquad+\sum_{k=1}^{2}\mathfrak{Z}^{k}_{i}(X_{t_ {i}};\Theta^{i},\Xi^{i})\Delta_{i}^{W^{k}}\bigg{\}}\Bigg{|}^{2}\Bigg{]}\] \[=\mathbb{E}^{\Phi}\Bigg{[}\Bigg{|}\widehat{\mathfrak{U}}_{i+1}(X _{t_{i+1}})-\bigg{\{}\mathfrak{U}_{i}(X_{t_{i}};\Theta^{i})-\Big{(}a\mathfrak{ U}_{i}(X_{t_{i}};\Theta^{i})+b\mathfrak{Z}^{1}_{i}(X_{t_{i}};\Theta^{i},\Xi^{i})+c \mathfrak{Z}^{2}_{i}(X_{t_{i}};\Theta^{i},\Xi^{i})+\widetilde{f}\Big{)}\, \delta_{i}\] \[\qquad\qquad\qquad\qquad+\sum_{k=1}^{2}\mathfrak{Z}^{k}_{i}(X_{t _{i}};\Theta^{i},\Xi^{i})\Delta_{i}^{W^{k}}\bigg{\}}\Bigg{|}^{2}\Bigg{]}\] \[=\mathbb{E}^{\Phi}\Bigg{[}\Bigg{|}\widehat{\mathfrak{U}}_{i+1}(X _{t_{i+1}})+\widetilde{f}\delta_{i}-\bigg{\{}\Xi^{i}\Phi^{\Xi,i}_{K}(X_{t_{i}}) \left(\Delta_{i}^{W^{1}}-b\delta_{i}\right)\] \[\qquad\qquad\qquad\qquad\qquad+\Theta^{i}\left((1-a\delta_{i}) \Phi^{\Theta,i}_{K}(X_{t_{i}})+\mathrm{D}_{x}\Phi^{\Theta,i}_{K}(X_{t_{i}}) \sqrt{V_{t_{i}}}\left(\Delta_{i}^{B}-(b\rho_{1}+c\rho_{2})\delta_{i}\right) \bigg{)}\,\Bigg{\}}\Bigg{|}^{2}\Bigg{]}\] \[=\mathbb{E}^{\Phi}\left[\big{|}\mathrm{Y}^{i}-\Xi^{i}\mathrm{X}^{ i}_{1}-\Theta^{i}\mathrm{X}^{i}_{2}\big{|}^{2}\right],\]
with \(\Delta_{i}^{B}=(\rho_{1}\Delta_{i}^{W^{1}}+\rho_{2}\Delta_{i}^{W^{2}})\) and where \(\widehat{\mathfrak{U}}_{i+1}(X_{t_{i+1}})\coloneqq\mathfrak{U}_{i+1}(X_{t_{i+1}}; \Theta^{i+1,*})\) was set in the previous time step and is now constant (without dependence on \(\Theta^{i}\)). We defined

\[\left\{\begin{array}{ll}\mathrm{Y}^{i}&\coloneqq\widehat{\mathfrak{U}}_{i+1} (X_{t_{i+1}})+\widetilde{f}\delta_{i},\\ \mathrm{X}_{1}^{i}&\coloneqq\Phi_{K}^{\Xi}(X_{t_{i}})\left(\Delta_{i}^{W^{1}} -b\delta_{i}\right),\\ \mathrm{X}_{2}^{i}&\coloneqq(1-a\delta_{i})\Phi_{K}^{\Theta}(X_{t_{i}})+D_{x }\Phi_{K}^{\Theta}(X_{t_{i}})\sqrt{V_{t_{i}}}\left(\Delta_{i}^{B}-(b\rho_{1}+c \rho_{2})\delta_{i}\right).\end{array}\right. \tag{4.6}\]

In matrix form, this yields \(\ell(\Theta^{i},\Xi^{i})=\mathbb{E}^{\Phi}[\|\mathrm{Y}^{i}-\boldsymbol{ \beta}^{i}\mathrm{X}^{i}\|^{2}]\), with \(\boldsymbol{\beta}^{i}=\left[\Xi^{i},\Theta^{i}\right]\) and \(\mathrm{X}^{i}=\left[\mathrm{X}_{1}^{i},\mathrm{X}_{2}^{i}\right]^{\top}\), for which the RLS from Section 2.2 yields the solution

\[\boldsymbol{\beta}^{i}=\mathbb{E}^{\Phi}\left[\left[\mathrm{Y}^{i}\mathrm{X} _{1}^{i\top}\quad\mathrm{Y}^{i}\mathrm{X}_{2}^{i\top}\right]\right]\mathbb{E} ^{\Phi}\left[\left[\begin{array}{cc}\mathrm{X}_{1}^{i}\mathrm{X}_{1}^{i\top }&\mathrm{X}_{1}^{i}\mathrm{X}_{2}^{i\top}\\ \mathrm{X}_{2}^{i}\mathrm{X}_{1}^{i\top}&\mathrm{X}_{2}^{i}\mathrm{X}_{2}^{i \top}\end{array}\right]\right]^{-1}. \tag{4.7}\]

### Algorithm

We summarise the steps of the algorithm below:

_Inputs_: \(\mathrm{time\,grid}\,\pi=\{0=t_{0}<t_{1}<\ldots<t_{N}=T\}\); number \(K\) of hidden nodes; \(R>0\);

_Initialisation_:

\(\left(\{\Phi_{K}^{\Theta,i}\}_{i=0}^{N-1},\{\Phi_{K}^{\Xi,i}\}_{i=0}^{N-1}\right)\) reservoirs with weights and biases distributed as \(\mathcal{U}_{[-R,R]}\);

**do**:

Generate \(n\) paths of \((\{X_{t_{i}}^{\pi,j}\}_{i=0}^{N},\{V_{t_{i}}^{\pi,j}\}_{i=0}^{N})\) for \(j\in\{1,\ldots,n\}\) with Euler-Maruyama

Set \(\widehat{Y}_{N}(x)=g(x)\) for all \(x\in\mathbb{R}\)

**for \(i\in\{N-1,\ldots,0\}\) do**:

Approximate \(u(t_{i},\cdot)\) with \(\mathfrak{U}(\cdot;\Theta^{i})\in\aleph_{K}^{\varsigma}\) based on reservoir \(\Phi_{K}^{\Theta,i}\)

Approximate \(\psi(t_{i},\cdot)\) with \(\Psi(\cdot;\Xi^{i})\in\aleph_{K}^{\varsigma}\) based on reservoir \(\Phi_{K}^{\Xi,i}\)

Evaluate derivatives of \(\left(\mathfrak{U}(\cdot;\Theta^{i}),\Psi(\cdot;\Xi^{i})\right)\) according to (2.1)

Solve the regression problem (possibly using the Ridge estimator from Section 2.2)

\[\Theta^{i,*}=\operatorname*{arg\,min}_{\Theta^{i}}\ell(\Theta^{i},\Xi^{i})= \operatorname*{arg\,min}_{\Theta^{i}}\mathbb{E}^{\Phi,n}\left[\left\|\mathrm{Y }^{i}-\boldsymbol{\beta}^{i}\mathrm{X}^{i}\right\|^{2}\right]\]

with \(\boldsymbol{\beta}^{i}=\left[\Xi^{i},\Theta^{i}\right]\), \(\mathrm{X}^{i}=\left[\mathrm{X}_{1}^{i},\mathrm{X}_{2}^{i}\right]^{\top}\), where \(\mathrm{Y}^{i},\mathrm{X}_{1}^{i},\mathrm{X}_{2}^{i}\) are given in (4.6), and \(\mathbb{E}^{\Phi,n}\) is computed with the empirical measure of \(\left(\{X_{t_{i}}^{\pi,j}\}_{i=0}^{N},\{V_{t_{i}}^{\pi,j}\}_{i=0}^{N}\right)_{j \in\{1,\ldots,n\}}\);

Update \(\widehat{\mathfrak{U}}_{i}(X_{t_{i}})=\mathfrak{U}_{i}\left(X_{t_{i}},\Theta^{ i,*}\right)\)

**end for**:

**return \(\{\mathfrak{U}_{i}(\cdot;\Theta^{i,*})\}_{i=0}^{N}\).**

## 5. Convergence analysis

In this section, whenever there is any ambiguity, we use the notation \(X^{\pi}\) to denote the discretised version of the solution process of (4.5) over the partition \(\pi=\{0=t_{0}<t_{1}<\ldots<t_{N}=T\}\) of the interval \([0,T]\), with modulus \(|\pi|=\max_{i\in\{0,1,\ldots,N-1\}}\delta_{i}\) with \(\delta_{i}=t_{i+1}-t_{i}\). As mentioned just before Assumption 3.1, the linearity of \(f\) assumed before was only required to cast the optimisation in Algorithm 2 into a regression problem. In the forthcoming convergence analysis, this does not play any role, and we, therefore, allow for a more general function \(f\).
**Assumption 5.1**.:
1. There exists a unique weak solution to the BSPDE system (4.3) with \(u,\psi\in\mathcal{W}^{3,2}\);
2. There is an increasing continuous function \(\omega:\mathbb{R}^{+}\to\mathbb{R}^{+}\) with \(\omega(0)=0\) such that \[\mathbb{E}\left[\int_{t_{1}}^{t_{2}}V_{s}\mathrm{d}s\right]+\mathbb{E}\left[ \left|\int_{t_{1}}^{t_{2}}V_{s}\mathrm{d}s\right|^{2}\right]\leq\omega(|t_{2}- t_{1}|),\quad\text{for any }0\leq t_{1}\leq t_{2}\leq T;\]
3. There exists \(L_{f}>0\) such that, for all \((t,x,z^{1},z^{2})\) and \((\widetilde{t},\widetilde{x},\widetilde{z}^{1},\widetilde{z}^{2})\), \[\left|f\left(t,\mathrm{e}^{x},y,z^{1},z^{2}\right)-f\left( \widetilde{t},\mathrm{e}^{\widetilde{x}},\widetilde{y},\widetilde{z}^{1}, \widetilde{z}^{2}\right)\right|\] \[\qquad\qquad\leq L_{f}\left\{\omega(|t-\widetilde{t}|)^{\frac{1} {2}}+|x-\widetilde{x}|+|y-\widetilde{y}|+|z^{1}-\widetilde{z}^{1}|+|z^{2}- \widetilde{z}^{2}|\right\}.\]

**Assumption 5.2**.: The first absolute moment of the discretisation scheme over the partition \(\pi\) for \(\{V_{s}\}_{s\in[0,T]}\) is bounded, namely \(\mathbb{E}\left[\left|V_{t_{i}}^{\pi}\right|\right]<\infty\) for all \(i\in\{0,\dots,N-1\}\).

Under Assumptions 4.1-5.1 Briand, Delyon, Hu, Pardoux, Stoica [14] established that

\[\mathbb{E}\left[\sup_{0\leq t\leq T}|X_{t}|^{2}\right] \leq C\left(1+|x_{0}|^{2}\right), \tag{5.1}\] \[\max_{i\in\{0,\dots,N-1\}}\mathbb{E}\left[\left|X_{t_{i+1}}-X_{t_{i+1}}^{ \pi}\right|^{2}+\sup_{t\in[t_{i},t_{i+1}]}\left|X_{t}-X_{t_{i}}^{\pi}\right|^{ 2}\right] \leq C\omega(|\pi|),\]

for some \(C>0\) independent of \(|\pi|\), and we furthermore have [14]

\[\mathbb{E}\left[\int_{0}^{T}\left|f\left(t,\mathrm{e}^{X_{t}},Y_{t},Z_{t}^{1},Z_{t}^{2}\right)\right|^{2}\mathrm{d}t\right]<\infty, \tag{5.2}\]

as well as the standard \(L^{2}\)-regularity result on \(Y\):

\[\max_{i\in\{0,\dots,N-1\}}\mathbb{E}\left[\sup_{t\in[t_{i},t_{i+1}]}\left|Y_{ t}-Y_{t_{i}}^{\pi}\right|^{2}\right]=\mathcal{O}(|\pi|). \tag{5.3}\]

For \(k\in\{1,2\}\), define the errors

\[\varepsilon^{Z^{k}}(\pi)\coloneqq\mathbb{E}\left[\sum_{i=0}^{N-1}\int_{t_{i}} ^{t_{i+1}}\left|Z_{t}^{k}-\overline{Z}_{t_{i}}^{k}\right|^{2}\mathrm{d}t\right],\quad\text{ with }\quad\overline{Z}_{t_{i}}^{k}\coloneqq\frac{1}{\delta_{i}} \mathbb{E}_{i}\left[\int_{t_{i}}^{t_{i+1}}Z_{t}^{k}\mathrm{d}t\right], \tag{5.4}\]

where \(\mathbb{E}_{i}\) denotes the conditional expectation given \(\mathcal{F}_{t_{i}}\). We furthermore define the auxiliary processes, for \(i\in\{0,\dots,N-1\}\),

\[\widehat{\mathcal{V}}_{t_{i}} \coloneqq\mathbb{E}_{i}\left[\widehat{\mathfrak{U}}_{i+1}\left(X_ {t_{i+1}}^{\pi}\right)\right]+f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}}, \widehat{\mathcal{V}}_{t_{i}},\overline{\widehat{Z}}_{t_{i}}^{1},\overline{ \widehat{Z}}_{t_{i}}^{2}\right)\delta_{i},\] \[\overline{\widehat{Z}}_{t_{i}}^{1} \coloneqq\widehat{\Psi}_{i}(X_{t_{i}}^{\pi})+\frac{1}{\delta_{i}} \mathbb{E}_{i}\left[\widehat{\mathfrak{U}}_{i+1}\left(X_{t_{i+1}}^{\pi}\right) \Delta_{i}^{W^{1}}\right],\] \[\overline{\widehat{Z}}_{t_{i}}^{2} \coloneqq\frac{1}{\delta_{i}}\mathbb{E}_{i}\left[\widehat{ \mathfrak{U}}_{i+1}\left(X_{t_{i+1}}^{\pi}\right)\Delta_{i}^{W^{2}}\right], \tag{5.5}\]

with \(\widehat{\mathfrak{U}}_{i}(\mathbf{x})\coloneqq\mathfrak{U}_{i}(x;\Theta^{i, *})\) and \(\widehat{\Psi}_{i}(x)\coloneqq\Psi(x;\Xi^{i,*})\) as before. Observe that \(\widehat{\Psi}_{i+1}\) and \(\widehat{\mathfrak{U}}_{i+1}\) do not depend on \(\Theta^{i}\) because the parameters were fixed at \((i+1)\) time step and are held
constant at step \(i\) (see Algorithm 2). Next, notice that \(\widehat{\mathcal{V}}\) is well defined by a fixed-point argument since \(f\) is Lipschitz. By Assumption 5.1(i), there exist \(\widehat{v}_{i},\widetilde{\widetilde{z}}_{t_{i}}^{k}\) for which

\[\widehat{\mathcal{V}}_{t_{i}}=\widehat{v}_{i}(X_{t_{i}}^{\pi})\qquad\text{and} \qquad\widetilde{\widetilde{Z}}_{t_{i}}^{k}=\widetilde{\widetilde{z}}_{t_{i}}^ {k}(X_{t_{i}}^{\pi})\quad\text{for }i\in\{0,\dots,N-1\},k\in\{1,2\}. \tag{5.6}\]

Then by the martingale representation theorem [17, Theorem 14.5.1], there exist integrable processes \(\overline{\widetilde{Z}}^{1}\) and \(\overline{\widetilde{Z}}^{2}\) such that

\[\widehat{\Omega}_{i+1}\left(X_{t_{i+1}}^{\pi}\right)=\widehat{\mathcal{V}}_{t _{i}}-f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\widehat{\mathcal{V}}_{t_{i}},\widetilde{\widetilde{Z}}_{t_{i}}^{1},\widetilde{\widetilde{Z}}_{t_{i}}^{2} \right)\delta_{i}+\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{t}^{1}\mathrm{d}W_{t}^{1 }+\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{t}^{2}\mathrm{d}W_{t}^{2}, \tag{5.7}\]

since \(\overline{\widetilde{Z}}_{t}^{k}\) are \(\mathcal{F}_{t}^{W}\)-adapted as asserted by the martingale representation theorem. From here, Ito's isometry yields

\[\overline{\widetilde{Z}}_{t_{i}}^{1} =\widehat{\Psi}_{i}(X_{t_{i}}^{\pi})+\frac{1}{\delta_{i}}\mathbb{ E}_{i}\left[\widehat{\Omega}_{i+1}(X_{t_{i+1}}^{\pi})\Delta_{i}^{W^{k}}\right]\] \[=\frac{1}{\delta_{i}}\int_{t_{i}}^{t_{i+1}}\widehat{\Psi}_{i}(X_ {t_{i}}^{\pi})\mathrm{d}t+\frac{1}{\delta_{i}}\mathbb{E}_{i}\left[\left( \widehat{\mathcal{V}}_{t_{i}}+\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{t}^{1} \mathrm{d}W_{t}^{1}+\int_{t_{i}}^{t_{i+1}}\widehat{Z}_{t}^{2}\mathrm{d}W_{t}^ {2}\right)\int_{t_{i}}^{t_{i+1}}\mathrm{d}W_{t}^{1}\right]\] \[=\frac{1}{\delta_{i}}\mathbb{E}_{i}\left[\int_{t_{i}}^{t_{i+1}} \left(\widehat{\Psi}_{i}(X_{t_{i}}^{\pi})+\widehat{Z}_{t}^{1}\right)\mathrm{d }t\right],\]

and similarly,

\[\overline{\widetilde{Z}}_{t_{i}}^{2}=\frac{1}{\delta_{i}}\mathbb{E}_{i}\left[ \int_{t_{i}}^{t_{i+1}}\widehat{Z}_{t}^{2}\mathrm{d}t\right].\]

We consider convergence in terms of the following error:

\[\mathscr{E}\left(\widehat{\mathfrak{U}},\widehat{\Psi}\right)\coloneqq\max_{i \in\{0,\dots,N-1\}}\mathbb{E}^{\Phi}\left[\left|Y_{t_{i}}-\widehat{\mathfrak{U }}_{i}\left(X_{t_{i}}^{\pi}\right)\right|^{2}\right]+\mathbb{E}^{\Phi}\left[ \sum_{i=0}^{N-1}\int_{t_{i}}^{t_{i+1}}\sum_{k=1}^{2}\left|Z_{t}^{k}-\widehat{ \mathcal{Z}}_{i}^{k}\left(X_{t_{i}}^{\pi}\right)\right|^{2}\mathrm{d}t\right],\]

with \(\widehat{\mathcal{Z}}_{i}^{1},\widehat{\mathcal{Z}}_{i}^{2}\) introduced before Lemma 5.7. We now state the main convergence result:

**Theorem 5.3**.: _Under Assumptions 4.1-A.2-5.1, there exists \(C>0\) such that_

\[\mathscr{E}\left(\widehat{\mathfrak{U}},\widehat{\Psi}\right)\leq C\left\{ \omega(|\pi|)+\mathbb{E}\left[\left|g(X_{T})-g(X_{T}^{\pi})\right|^{2}\right]+ \sum_{k=1}^{2}\varepsilon^{Z^{k}}(\pi)+\frac{C^{*}}{K}N+M|\pi|^{2}\right\}\]

_over a compact \(\mathcal{Q}\subset\mathbb{R}\), with \(C^{*},M>0\) given in Lemma 5.8._

The following corollary is immediate from (C.5), established in Part II of the proof of Theorem 5.3:

**Corollary 5.4**.: _Under Assumptions 4.1-A.2-5.1, there exists \(C>0\) such that_

\[\max_{i\in\{0,\dots,N-1\}} \mathbb{E}^{\Phi}\left[\left|Y_{t_{i}}-\widehat{\mathfrak{U}}_{i} \left(X_{t_{i}}^{\pi}\right)\right|^{2}\right]\leq\] \[C\left\{\omega(|\pi|)+\mathbb{E}\left[\left|g(X_{T})-g(X_{T}^{ \pi})\right|^{2}\right]+\sum_{k=1}^{2}\varepsilon^{Z^{k}}(\pi)+\frac{C^{*}}{K} N+M|\pi|^{2}\right\}\]

_over a compact \(\mathcal{Q}\subset\mathbb{R}\), with \(C^{*},M>0\) given in Lemma 5.8._

**Remark 5.5**.: The second error term is the strong \(L^{2}\)-Monte-Carlo error and is \(\mathcal{O}(N^{-H})\) for processes driven by an fBm with Hurst parameter \(H\in(0,1)\). We refer the reader to [13, 27] for an exposition on strong versus weak error rates in rough volatility models.
To prove Theorem 5.3, the following bounds on the derivatives are key.

**Lemma 5.6**.: _Let \(\Psi_{K}(\cdot;\Theta)\in\mathbb{N}_{K}^{\varsigma}\) and \((X_{t_{i}}^{\pi},V_{t_{i}}^{\pi})_{i}\) denote the discretised versions of (4.5) over the partition \(\pi\), then there exist \(L_{1},L_{2}>0\) such that, for all \(i\in\{0,\ldots,N-1\}\),_

\[\left\|\mathbb{E}_{i}^{\Phi}\left[\mathrm{D}_{x}\Psi_{K}(X_{t_{i+1}}^{\pi}; \Theta)\right]\right\|\leq L_{1}\qquad\text{and}\qquad\left\|\mathbb{E}_{i}^{ \Phi}\left[\mathrm{D}_{x}^{2}\Psi_{K}(X_{t_{i+1}}^{\pi};\Theta)\right]\right\| \leq L_{2}.\]

Proof.: We start with the first derivative. For all \(x,y\in\mathbb{R}^{d}\),

\[\|\Psi_{K}(x;\Theta)-\Psi_{K}(y;\Theta)\| =\|\Theta\left(\upvarsigma(\mathrm{A}x+b)-\upvarsigma(Ay+b) \right)\|\leq\|\Theta\|_{F}\|\upvarsigma(\mathrm{A}x+b)-\upvarsigma(Ay+b)\|\] \[\leq\|\Theta\|_{F}\|\mathrm{A}x-Ay\|\leq\|\Theta\|_{F}\|A\|_{F}\| x-y\|\leq L_{1}\|x-y\|,\]

since \(\varsigma\) is \(1\)-Lipschitz. The estimator \(\Theta\) has an explicit form (4.7) and its norm is finite, therefore \(\Psi_{K}(\cdot;\Theta)\) is globally Lipschitz and its first derivative is bounded by \(L_{1}>0\). Next, without loss of generality, we can set \(\mathrm{A}=\mathrm{I}\) and \(\mathrm{b}=0\), since their support is bounded. As in Section 2.1, for the component \(j\in\{1,\ldots,m\}\),

\[\mathbb{E}_{i}^{\Phi}\left[\mathrm{D}_{x}^{2}\Psi_{K}(X_{t_{i+1}} ^{\pi};\Theta)\right]_{j} =\int\Theta\operatorname{diag}\left(e_{j}\right)\operatorname{ diag}\left(\mathbf{\delta}_{0}\left(x-\frac{1}{2}V_{t_{i}}^{\pi}\delta_{i}+\sqrt{V_{t_ {i}}}w\right)\right)\mathbf{p}_{\mathcal{N}}(w)\mathrm{d}w\] \[=\Theta\operatorname{diag}\left(e_{j}\right)\operatorname{diag} \left(\mathbf{p}_{\mathcal{N}}\left(0;x-\frac{1}{2}V_{t_{i}}^{\pi}\delta_{i},V_{t_ {i}}^{\pi}\delta_{i}\right)\right),\]

since \(\Delta_{i}^{B}\sim\mathcal{N}(0,\delta_{i})\) and \(\mathbf{p}_{\mathcal{N}}\) is the Gaussian density applied component-wise. Since the weights are sampled on a compact and \(\|\Theta\|\) is finite, then there exists \(C>0\) such that

\[\left\|\mathbb{E}_{i}^{\Phi}\left[\mathrm{D}_{x}^{2}\Psi_{K}(X_{t_{i+1}}^{\pi} ;\Theta)\right]\right\|\leq C\|\Theta\|_{F}=L_{2}.\]

From here the error bound of approximating \(\widehat{\mathcal{V}}_{t_{i}}\), \(\overline{\widehat{Z}}_{t_{i}}^{1}\) and \(\overline{\widehat{Z}}_{t_{i}}^{2}\) with their RWNN approximators \(\widehat{\mathfrak{U}}_{i}\), \(\widehat{\mathcal{Z}}_{i}^{1}\) and \(\widehat{\mathcal{Z}}_{i}^{2}\) (defined in the lemma below) can be obtained. For \(i\in\{0,\ldots,N-1\}\), \((\mathfrak{U}_{i},\Psi_{i})\in\mathbb{N}_{K}^{\varsigma}\), introduce

\[\mathcal{Z}_{i}^{1}(x) \coloneqq\Psi_{i}(x)+\rho_{1}\sqrt{V_{t_{i}}}\,\mathrm{D}_{x} \mathfrak{U}_{i}(x),\ \ \mathcal{Z}_{i}^{2}(x) \coloneqq\rho_{2}\sqrt{V_{t_{i}}}\,\mathrm{D}_{x} \mathfrak{U}_{i}(x),\] \[\widehat{\mathcal{Z}}_{i}^{1}(x) \coloneqq\widehat{\Psi}_{i}(x)+\rho_{1}\sqrt{V_{t_{i}}}\,\mathrm{ D}_{x}\widehat{\mathfrak{U}}_{i}(x),\ \ \widehat{\mathcal{Z}}_{i}^{2}(x) \coloneqq\rho_{2}\sqrt{V_{t_{i}}}\,\mathrm{D}_{x} \widehat{\mathfrak{U}}_{i}(x).\]

**Lemma 5.7**.: _There exists \(M>0\) such that, for any \(i\in\{0,\ldots,N-1\}\), \(k=1,2\),_

\[\mathbb{E}^{\Phi}\left[\left|\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi})-\overline{ \widehat{Z}}_{t_{i}}^{k}\right|^{2}\right]\leq\rho_{k}^{2}|\pi|^{2}M.\]

Proof.: From (5.5) and (5.6), we have, for \(i\in\{0,\ldots,N-1\}\) and \(k\in\{1,2\}\),

\[\widehat{v}_{i}(x) =\mathbb{E}_{i}^{\Phi}\left[\widehat{\mathfrak{U}}_{i+1}\left(X_ {t_{i+1}}^{x,\pi}\right)\right]+f\left(t_{i},\mathrm{e}^{x},\widehat{v}_{i}(x),\overline{\widehat{z}}_{i}^{1}(x),\overline{\widehat{z}}_{i}^{2}(x)\right) \delta_{i},\] \[\overline{\widehat{z}}_{i}^{k}(x) =\Psi_{i}\left(X_{t_{i+1}}^{x,\pi}\right)\!\!1\!\!1_{\{k=1\}}+ \frac{1}{\delta_{i}}\mathbb{E}_{i}^{\Phi}\left[\widehat{\mathfrak{U}}_{i+1} \left(X_{t_{i+1}}^{x,\pi}\right)\Delta_{i}^{W^{k}}\right],\]

where \(X_{t_{i+1}}^{x,\pi}=x+\left(r-\frac{1}{2}V_{t_{i}}\right)\delta_{i}+\sqrt{V_{t _{i}}}\,\Delta_{i}^{B}\) is the Euler discretisation of \(\{X_{t}\}_{t\in[0,T]}\) over \(\pi\) and \(\{V_{t_{i}}^{\pi}\}_{i=0}^{N}\) is the appropriate discretisation of the volatility process over the same partition. For \(\{\mathcal{R}^{k}\}\stackrel{{\text{iid}}}{{\sim}}\mathcal{N}(0,1)\), the two auxiliary processes can be written as

\[\overline{\widehat{z}}_{i}^{k}(x)=\Psi_{i}\left(X_{t_{i+1}}^{x,\pi}\right)\!\!1 \!\!1_{\{k=1\}}+\frac{1}{\delta_{i}}\mathbb{E}_{i}^{\Phi}\left[\widehat{ \mathfrak{U}}_{i+1}\left(x-\frac{1}{2}V_{t_{i}}^{x,\pi}\delta_{i}+\sqrt{V_{t_{i} }^{x,\pi}}\sqrt{\delta_{i}}\left(\rho_{1}\mathcal{R}^{1}+\rho_{2}\mathcal{R}^{2} \right)\right)\sqrt{\delta_{i}}\mathcal{R}^{k}\right].\]
Notice that, while any sensible forward scheme for \(\{V_{t}\}_{t\in[0,T]}\) does depend on a series of Brownian increments, \(V_{t_{i}}^{x,\pi}\) only depends on \(\left(\Delta_{i-1}^{W},\ldots,\Delta_{0}^{W}\right)\), which are known at time \(t_{i}\). Thus, since usual derivative operations are available for approximately differentiable functions (Remark 2.7) multivariate integration by parts for Gaussian measures (a formulation of Isserlis' Theorem [44]) yields

\[\widetilde{\widetilde{z}}_{i}^{\,k}(x)=\Psi_{i}\left(X_{t_{i+1}}^{x,\pi}\right) \mathbf{1}_{\{k=1\}}+\rho_{k}\sqrt{V_{t_{i}}}\mathbb{E}^{\Phi}\left[\mathrm{D} _{x}\widehat{\Omega}_{i+1}\left(X_{t_{i+1}}^{x,\pi}\right)\right],\]

with corresponding derivatives

\[\mathrm{D}_{x}\widetilde{\widetilde{z}}_{i}^{\,k}(x)=\mathrm{D}_{x}\Psi_{i} \left(X_{t_{i+1}}^{x,\pi}\right)\mathbf{1}_{\{k=1\}}+\rho_{k}\sqrt{V_{t_{i}}} \mathbb{E}^{\Phi}\left[\mathrm{D}_{x}^{2}\widehat{\Omega}_{i+1}\left(X_{t_{i +1}}^{x,\pi}\right)\right].\]

An application of the implicit function theorem then implies

\[\mathrm{D}_{x}\widehat{v}_{i}(x)=\mathbb{E}_{i}^{\Phi}\left[\mathrm{D}_{x} \widehat{\Omega}_{i+1}\left(X_{t_{i+1}}^{x,\pi}\right)\right]+\delta_{i} \left\{\mathrm{D}_{x}\widehat{f}_{i}(x)+\mathrm{D}_{y}\widehat{f}_{i}(x) \mathrm{D}_{x}\widehat{v}_{i}(x)+\sum_{k=1}^{2}\mathrm{D}_{z^{k}}\widehat{f}_{ i}(x)\mathrm{D}_{x}\widetilde{\widetilde{z}}_{i}^{\,k}(x)\right\},\]

where \(\widehat{f}_{i}(x)\coloneqq f\left(t_{i},\mathrm{e}^{x},\widehat{v}_{i}(x), \widetilde{\widetilde{z}}_{i}^{\,1}(x),\widetilde{\widetilde{z}}_{i}^{\,2}(x )\right)\) and

\[\Psi_{i}\left(X_{t_{i+1}}^{x,\pi}\right)\mathbf{1}_{\{k=1\}}+ \left(1-\delta_{i}\mathrm{D}_{y}\widehat{f}_{i}(x)\right)\rho_{k}\sqrt{V_{t_{i }}^{\pi}}\mathrm{D}_{x}\widehat{v}_{i}(x)\] \[=\widetilde{\widetilde{z}}_{i}^{k}(x)+\rho_{k}\sqrt{V_{t_{i}}^{ \pi}}\delta_{i}\left(\mathrm{D}_{x}\widehat{f}_{i}(x)+\mathrm{D}_{z^{1}} \widehat{f}_{i}(x)\mathrm{D}_{x}\widetilde{\widetilde{z}}_{i}^{\,1}(x)+ \mathrm{D}_{z^{2}}\widehat{f}_{i}(x)\mathrm{D}_{x}\widetilde{\widetilde{z}}_{ i}^{\,2}(x)\right).\]

Thus, for small enough \(|\pi|\), since \(f\) is Lipschitz (Assumption 5.1) and Lemma 5.6,

\[\Psi_{i}\left(X_{t_{i+1}}^{x,\pi}\right)\mathbf{1}_{\{k=1\}}+ \rho_{k}\sqrt{V_{t_{i}}^{\pi}}\mathrm{D}_{x}\widehat{v}_{i}(x) \leq\overline{\widetilde{z}}_{i}^{k}(x)+\rho_{k}\delta_{i}\sqrt{V_ {t_{i}}^{\pi}}\left(1+\mathrm{D}_{x}\widetilde{\widetilde{z}}_{i}^{\,1}(x)+ \mathrm{D}_{x}\widetilde{\widetilde{z}}_{i}^{\,2}(x)\right)\] \[\leq\overline{\widetilde{z}}_{i}^{k}(x)+\rho_{k}\delta_{i}\sqrt{V _{t_{i}}^{\pi}}\left(1+L_{1}+\sqrt{2}L_{2}\sqrt{V_{t_{i}}^{\pi}}\right).\]

Therefore

\[\mathbb{E}^{\Phi}\left[\left|\Psi_{i}(X_{t_{i}}^{\pi})+\rho_{1} \sqrt{V_{t_{i}}^{\pi}}\mathrm{D}_{x}\mathfrak{U}_{i}(X_{t_{i}}^{\pi})- \widetilde{\widetilde{Z}}_{t_{i}}^{\,1}\right|^{2}\right]\] \[\leq\mathbb{E}^{\Phi}\left[\left|\overline{\widetilde{z}}_{i}^{\, 1}(X_{t_{i}}^{\pi})-\overline{\widetilde{Z}}_{t_{i}}^{\,1}+\rho_{1}\delta_{i} \sqrt{V_{t_{i}}^{\pi}}\left[1+L_{1}+\sqrt{2}L_{2}\sqrt{V_{t_{i}}^{\pi}}\right] \right|^{2}\right]\] \[\leq\rho_{1}^{2}|\pi|^{2}\left\{\mathbb{E}\left[\left|V_{t_{i}}^{ \pi}\right|\right]\left(1+L_{1}+\sqrt{2}L_{2}\mathbb{E}\left[\left|V_{t_{i}}^{ \pi}\right|\right]\right)\right\}\leq\rho_{1}^{2}|\pi|^{2}M,\]

using Corollary 2.10 in the second inequality and the boundedness of \(\mathbb{E}[|V^{\pi}|]\) from Assumption 5.2 in the last line. The proof of the other bound is analogous. 

**Lemma 5.8**.: _Under Assumptions 4.1-A.2-5.1, for sufficiently small \(|\pi|\) we have_

\[\mathbb{E}^{\Phi}\left[\left|\widehat{\mathcal{V}}_{t_{i}}-\widehat{\Omega}_{i} \left(X_{t_{i}}^{\pi}\right)\right|^{2}\right]+\delta_{i}\mathbb{E}^{\Phi} \left[\sum_{k=1}^{2}\left|\widehat{Z}_{t_{i}}^{k}-\widehat{Z}_{i}^{k}\left(X_{t _{i}}^{\pi}\right)\right|^{2}\right]\leq C\left\{\frac{C^{*}}{K}+M|\pi|^{3}\right\}\]

_on a compact \(\mathcal{Q}\), for all \(i\in\{0,\ldots,N-1\}\) and \(K\) hidden units, for some \(C>0\), where \(C^{*}\) is as in Proposition B.1 and \(M\) in Lemma 5.7._
Proof of Lemma 5.8.: Fix \(i\in\{0,\ldots,N-1\}\). Relying on the martingale representation in (5.7) and Lemma B.2, we can define the following loss function for the pair \((\mathfrak{U}_{i}(\cdot;\Theta),\Psi_{i}(\cdot;\Xi))\in\aleph_{K}^{\kappa}\) and their corresponding parameters \(\Theta\) and \(\Xi\):

\[\widehat{L}_{i}(\Theta,\Xi)\coloneqq\widetilde{L}_{i}(\Theta,\Xi)+\mathbb{E}^ {\Phi}\left[\int_{t_{i}}^{t_{i+1}}\sum_{k=1}^{2}\left|\widehat{Z}_{t}^{k}- \overline{\widetilde{Z}}_{t_{i}}^{k}\right|^{2}\right], \tag{5.8}\]

with

\[\widetilde{L}_{i}(\Theta,\Xi)\coloneqq\mathbb{E}^{\Phi}\bigg{[} \bigg{|}\widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta) +\delta_{i}\Big{\{}f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\mathfrak{U}(X_{ t_{i}}^{\pi};\Theta),\mathcal{Z}_{i}^{1}(X_{t_{i}}^{\pi};\Theta,\Xi), \mathcal{Z}_{i}^{2}(X_{t_{i}}^{\pi};\Theta,\Xi)\right)\] \[-f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\mathcal{V}_{t_{i}}, \overline{\widetilde{Z}}_{t_{i}}^{1},\overline{\widetilde{Z}}_{t_{i}}^{2} \right)\Big{\}}\Big{|}^{2}\bigg{]}+\delta_{i}\sum_{k=1}^{2}\mathbb{E}^{\Phi} \left[\left|\overline{\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i} }^{\pi};\Theta,\Xi)\right|^{2}\right].\]

Recall the following useful identity, valid for any \(a,b\in\mathbb{R}\) and \(\chi>0\):

\[(a+b)^{2}\leq(1+\chi)\,a^{2}+\left(1+\frac{1}{\chi}\right)b^{2}. \tag{5.9}\]

Applying (5.9) with \(\chi=\gamma\delta_{i}\) yields

\[\widetilde{L}_{i}(\Theta,\Xi)\leq\delta_{i}\sum_{k=1}^{2}\mathbb{ E}^{\Phi}\left[\left|\overline{\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_ {i}}^{\pi};\Theta,\Xi)\right|^{2}\right]+(1+C\delta_{i})\mathbb{E}^{\Phi}\left[ \left|\widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta) \right|^{2}\right]\] \[+\left(1+\frac{1}{C\delta_{i}}\right)\mathbb{E}^{\Phi}\left[ \left|f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\mathfrak{U}(X_{t_{i}}^{\pi}; \Theta),\mathcal{Z}_{i}^{1}(X_{t_{i}}^{\pi};\Theta,\Xi),\mathcal{Z}_{i}^{2}(X _{t_{i}}^{\pi};\Theta,\Xi)\right)-f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}}, \mathcal{V}_{t_{i}},\overline{\widetilde{Z}}_{t_{i}}^{1},\overline{\widetilde{ Z}}_{t_{i}}^{2}\right)\right|^{2}\right].\]

Now by the Lipschitz condition on \(f\) from Assumption 5.1,

\[\widetilde{L}_{i}(\Theta,\Xi)\leq(1+C\delta_{i})\mathbb{E}^{\Phi}\left[\left| \widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta)\right| ^{2}\right]+C\delta_{i}\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left|\overline{ \widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi};\Theta,\Xi) \right|^{2}\right].\]

For any \(a,b\in\mathbb{R}\), \(\chi<0\), Inequality (5.9) holds with the reverse sign, hence, with \(\chi=-\gamma\delta_{i}\),

\[\widetilde{L}_{i}(\Theta,\Xi)\geq\delta_{i}\sum_{k=1}^{2}\mathbb{ E}^{\Phi}\left[\left|\overline{\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_ {i}}^{\pi};\Theta,\Xi)\right|^{2}\right]+(1-\gamma\delta_{i})\mathbb{E}^{ \Phi}\left[\left|\widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{ \pi};\Theta)\right|^{2}\right]\] \[-\frac{1}{\gamma\delta_{i}}\mathbb{E}^{\Phi}\left[\left|f\Big{(}t _{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\mathfrak{U}(X_{t_{i}}^{\pi};\Theta), \mathcal{Z}_{i}^{1}(X_{t_{i}}^{\pi};\Theta,\Xi),\mathcal{Z}_{i}^{2}(X_{t_{i}}^ {\pi};\Theta,\Xi)\right)-f\left(t_{i},\mathrm{e}^{X_{t_{i}}^{\pi}},\mathcal{V }_{t_{i}},\overline{\widetilde{Z}}_{t_{i}}^{1},\overline{\widetilde{Z}}_{t_{i} }^{2}\Big{)}\right|^{2}\right].\]
Again since \(f\) is Lipschitz, the arithmetic-geometric inequality implies

\[\widetilde{L}_{i}(\Theta,\Xi) \geq(1-\gamma\delta_{i})\mathbb{E}^{\Phi}\left[\left|\widehat{ \mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta)\right|^{2} \right]+\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left|\widetilde{\widetilde{Z}}_{t _{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi};\Theta,\Xi)\right|^{2}\right]\] \[\geq(1-\gamma\delta_{i})\mathbb{E}^{\Phi}\left[\left|\widehat{ \mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta)\right|^{2} \right]+\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left|\widetilde{\widetilde{Z}}_{t _{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi};\Theta,\Xi)\right|^{2}\right]\] \[\qquad-\frac{3\delta_{i}L_{f}^{2}}{\gamma}\left(\mathbb{E}^{\Phi }\left[\left|\widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi}; \Theta)\right|^{2}\right]+\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left|\widetilde {\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi};\Theta,\Xi) \right|^{2}\right]\right).\]

Taking \(\gamma=6L_{f}^{2}\) gives

\[\widetilde{L}_{i}(\Theta,\Xi)\geq(1-C\delta_{i})\mathbb{E}^{\Phi}\left[ \left|\widehat{\mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta) \right|^{2}\right]+\frac{\delta_{i}}{2}\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[ \left|\widetilde{\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{ \pi};\Theta,\Xi)\right|^{2}\right].\]

For a given \(i\in\{0,\ldots,N-1\}\) take \((\Theta^{*},\Xi^{*})\in\arg\min_{\Theta,\Xi}\widehat{L}_{i}(\Theta,\Xi)\) so that \(\widehat{\mathfrak{U}}_{i}=\mathfrak{U}_{i}(\cdot;\Theta^{*})\) and \(\widehat{\mathcal{Z}}_{i}^{k}(\cdot)\coloneqq\mathcal{Z}_{i}^{k}(\cdot;\Theta ^{*},\Xi^{*})\). From (5.8), \(\widetilde{L}_{i}\) and \(\widetilde{L}_{i}\) have the same minimisers, thus combining both bounds gives for all \((\Theta,\Xi)\in\mathbb{R}^{m\times K}\times\mathbb{R}^{m\times K}\),

\[(1-C\delta_{i})\,\mathbb{E}^{\Phi}\left[\left|\widehat{\mathcal{ V}}_{t_{i}}-\widehat{\mathfrak{U}}_{i}\left(X_{t_{i}}^{\pi}\right)\right|^{2} \right]+\frac{\delta_{i}}{2}\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left| \widetilde{\widetilde{Z}}_{t_{i}}^{k}-\widehat{\mathcal{Z}}_{i}^{k}\right|^{2 }\right]\leq\widetilde{L}_{i}(\Theta^{*},\Xi^{*})\leq\widetilde{L}_{i}(\Theta,\Xi)\] \[\leq(1+C\delta_{i})\,\mathbb{E}^{\Phi}\left[\left|\widehat{ \mathcal{V}}_{t_{i}}-\mathfrak{U}_{i}\left(X_{t_{i}}^{\pi};\Theta\right) \right|^{2}\right]+C\delta_{i}\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left| \widetilde{\widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi}; \Theta,\Xi)\right|^{2}\right].\]

Letting \(|\pi|\) be sufficiently small gives together with Lemma 5.7

\[\mathbb{E}^{\Phi}\left[\left|\widehat{\mathcal{V}}_{t_{i}}- \widehat{\mathfrak{U}}_{i}\left(X_{t_{i}}^{\pi}\right)\right|^{2}\right] +\delta_{i}\sum_{k=1}^{2}\mathbb{E}^{\Phi}\left[\left|\widetilde{ \widetilde{Z}}_{t_{i}}^{k}-\mathcal{Z}_{i}^{k}(X_{t_{i}}^{\pi};\Theta,\Xi) \right|^{2}\right]\] \[\leq C\left\{\inf_{\Theta}\mathbb{E}^{\Phi}\left[\left|\widehat{ v}_{i}(X_{t_{i}}^{\pi})-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta)\right|^{2}\right]+|\pi|^{3} \left(\rho_{1}^{2}+\rho_{2}^{2}\right)M\right\},\]

therefore

\[\mathbb{E}^{\Phi}\left[\left|\widehat{\mathcal{V}}_{t_{i}}- \widehat{\mathfrak{U}}_{i}\left(X_{t_{i}}^{\pi}\right)\right|^{2}\right]+ \delta_{i}\mathbb{E}^{\Phi}\left[\sum_{k=1}^{2}\left|\widetilde{\widetilde{Z}}_ {t_{i}}^{k}-\widehat{\mathcal{Z}}_{i}^{k}\left(X_{t_{i}}^{\pi}\right)\right|^{ 2}\right]\] \[\leq C\left\{\inf_{\Theta}\mathbb{E}^{\Phi}\left[\left|\widehat{ v}_{i}(X_{t_{i}}^{\pi})-\mathfrak{U}_{i}(X_{t_{i}}^{\pi};\Theta)\right|^{2}\right]+M|\pi|^{3} \right\}\leq C\left\{\frac{C^{*}}{K}+M|\pi|^{3}\right\},\]

by Proposition B.1 over any compact \(\mathcal{Q}=\{x\in\mathbb{R}:|x|\leq Q\}\), \(Q>0\). 

The rest of the proof is similar to the proofs of [5, Theorem A.2] and [42, Theorem 4.1], but we include it in Appendix C for the sake of completeness.
## 6. Numerical results

We now showcase the performance of the RWNN scheme on a representative model from each-Markovian and non-Markovian-class. We first test our scheme in the multi-dimensional Black-Scholes (BS) setting [11] and then move to the non-Markovian setup with the rough Bergomi (rBergomi) model [4]. We develop numerical approximations to European option prices given in (3.3) and (4.2), choosing

\[f(t,x,y,z^{1},z^{2})=-ry\qquad\text{and}\qquad g_{\text{call}}\left(x\right)= \left(\text{e}^{x}-\mathscr{K}\right)^{+},\]

and discretising over the partition \(\pi=\{0=t_{0},t_{1},\ldots t_{N}=T\}\) for some \(N\in\mathbb{N}\). The precise discretisation schemes of the individual processes are given in their corresponding sections below. We remark, however, that the approximated option price for a given Monte-Carlo sample can become (slightly) negative by construction so we add an absorption feature for both models:

\[Y_{t_{i}}^{\pi}\coloneqq\max\left\{0,\widetilde{Y}_{t_{i}}^{\pi}\right\}, \qquad\text{for }i\in\{0,\ldots,N-1\}, \tag{6.1}\]

where \(\left\{\widetilde{Y}_{t_{i}}^{\pi}\right\}_{i=0}^{N}\) denotes the approximation obtained through the RWNN scheme.

**Remark 6.1**.: This is a well-studied problem and is especially prevalent in the simulation of square-root diffusions. We acknowledge that the absorption scheme possibly creates additional bias (see [50] for the case of the Heston model), however, a theoretical study in the case of the PDE-RWWN scheme is out of the scope of this paper.

The reservoir used as a random basis of RWNNs here is the classical linear reservoir from Definition 2.2. For numerical purposes, we introduce a so-called _connectivity_ parameter, a measure of how interconnected the neurons in a network are: the higher the connectivity, the more inter-dependence between the neurons (see [18] for effects of connectivity in different reservoir topologies). In practice, however, too high a connectivity can lead to overfitting and poor generalisation. Recall that our reservoir is given by

\[\Phi_{K}:\mathbb{R}^{d}\to\mathbb{R}^{K},\qquad x\mapsto\Phi_{K}(x)\coloneqq \boldsymbol{\varrho}(\text{A}x+\text{b}),\]

where only \(\text{A}\in\mathbb{R}^{K\times d}\) is affected by the connectivity parameter \(c\in(0,1]\). A value \(c=1\) means that A is dense and fully determined by sampled weights. We find that the choice \(c\approx 0.5\) results in superior performance.

All the experiments below were run on a standard laptop with an AMD Ryzen 9 5900HX processor without any use of GPU, which would most certainly speed up the algorithms further. The code for both models is available on GitHub at ZuricZ/RWNN_PDE_solver.

### Example: Black-Scholes

The Black-Scholes model [11] is ubiquitous in mathematical finance, allowing for closed-form pricing and hedging of many financial contracts. Despite its well-known limitations, it remains a reference and is the first model to check before exploring more sophisticated ones. Since it offers an analytical pricing formula as a benchmark for numerical results, it will serve as a proof of concept for our numerical scheme. Under the pricing measure \(\mathbb{Q}\), the underlying assets \(\boldsymbol{S}=(S^{1},\ldots,S^{d})\) satisfy

\[\text{d}S_{t}^{j}=rS_{t}^{j}\text{d}t+\sigma_{j}S_{t}^{j}\text{d}W_{t}^{j}, \quad\text{for }t\in[0,T],j\in\{1,\ldots,d\},\]

where \(\{W_{t}^{j}\}_{t\in[0,T]}\) are standard Brownian motions such that \(\langle W^{i},W^{j}\rangle=\rho_{i,j}\text{d}t\), with \(\rho_{i,j}\in[-1,1]\), \(r\geq 0\) is the risk-free rate and \(\sigma_{j}>0\) is the volatility coefficient. The
corresponding \(d\)-dimensional option pricing PDE is then given by

\[\frac{\partial u(t,\mathbf{S})}{\partial t}+\sum_{j=1}^{d}rS^{j}\frac{\partial u(t, \mathbf{S})}{\partial S^{j}}+\sum_{j=1}^{d}\frac{(\sigma_{i}S^{j})^{2}}{2}\frac{ \partial^{2}u(t,\mathbf{S})}{(\partial S^{j})^{2}}+\sum_{j=1}^{d-1}\sum_{k=j+1}^{d} \rho_{j,k}\sigma_{j}\sigma_{k}S^{j}S^{k}\frac{\partial^{2}u(t,\mathbf{S})}{\partial S ^{j}\partial S^{k}}=ru(t,\mathbf{S}),\]

for \(t\in[0,T)\) with terminal condition \(u(T,\mathbf{S}_{T})=g(S_{T}^{1},\ldots,S_{T}^{d})\).

To use Algorithm 1, the process \(S\) has to be discretised, for example using the standard log-Euler-Maruyama scheme, for each \(j=1,\ldots,d\) and \(i\in\{0,1,\ldots,N-1\}\):

\[\left\{\begin{array}{ll}X_{t_{i+1}}^{\pi,j}&=X_{t_{i}}^{\pi,j}+\left(r- \frac{\sigma_{j}^{2}}{2}\right)\delta_{i}+\sigma_{j}\Delta_{i}^{W^{j}},\\ S_{t_{i+1}}^{\pi,j}&=\exp\left\{X_{t_{i+1}}^{\pi,j}\right\},\end{array}\right.\]

with initial value \(X_{0}^{\pi,j}=\log\left(S_{0}^{\pi,j}\right)\). If not stated otherwise, we let \(\mathscr{K}=S_{0}=1\), \(r=0.01\), \(T=1\), and run the numerical scheme with \(N=21\) discretisation steps and \(n_{\mathrm{MC}}=50,000\) Monte-Carlo samples. The reservoir has \(K\in\{10,100,1000\}\) hidden nodes, in Sections 6.1.3 and 6.1.2 the connectivity parameter is set to \(c=0.5\).

#### 6.1.1. Convergence rate

We empirically analyse the error rate in terms of the number of hidden nodes \(K\) obtained in Corollary 5.4. To isolate the dependence on the number of nodes, we consider a single ATM vanilla Call, fix \(c=1\), \(n_{\mathrm{MC}}=50,000\), \(\sigma=0.1\) and vary \(K\in\{10,100,1000\}\). Due to our vectorised implementation of the algorithm, the reservoir basis tensor cannot fit into the random-access memory of a standard laptop for \(K\geq 10000\). The results are compared to the theoretical price only computed using the Black-Scholes pricing formula. The absorption scheme (6.1) is applied.

#### 6.1.2. ATM Call option

As a proof of concept we first test Algorithm 1 with Call options written on \(d\in\mathbb{N}\) independent assets, i.e. with \(\rho_{j,k}=0\) for \(j\neq k\) and \(V_{T}=\left(g_{\mathrm{call}}(S_{T}^{j})\right)_{j\in\{1,\ldots,d\}}\). This is only done so that the results can be directly compared to the theoretical price computed using the Black-Scholes pricing formula and is, in effect, the same as pricing \(d\) options on \(d\) independent assets, each with their own volatility parameter \(\sigma\). All the listed results in this section are for \(K=100\) hidden nodes.

In Table 1, results and relative errors are shown for \(d=5\) and \(\mathbf{\sigma}\coloneqq(\sigma_{1},\ldots,\sigma_{d})\) uniformly spaced over \([0.05,0.4]\). Next, the effects of the absorption scheme (6.1) are investigated. Curiously, the absorption performs noticeably worse compared to the basic scheme, where one does not adjust for negative paths. This leads us to believe that absorption adds a substantial bias, similar to the Heston case (see Remark 6.1). Therefore, such a scheme should only be used for purposes, when positivity of the option price paths is strictly necessary (e.g. when hedging). Finally, in Table 2, total MSE and computational times are given for different dimensions. The computational times for different dimensions are then plotted in Figure 2. It is important to note that our results do not allow us to make definitive claims about the computational times of the PDE-RWNN scheme across different dimensions. This was not the goal of our experiments, and further detailed theoretical study and experiments would be necessary to draw more definitive conclusions regarding the efficiency of the scheme in various dimensions.
\begin{table}
\begin{tabular}{l|c c c c}  & \multicolumn{3}{c}{Price} \\ \hline \(\sigma\) & True & PDE (w/ abs) & PDE (w/o abs) & MC \\ \hline
0.05 & 0.02521640 & 0.02960256 & 0.02531131 & 0.02574731 \\
0.1 & 0.04485236 & 0.05523114 & 0.04467687 & 0.04547565 \\
0.15 & 0.06459483 & 0.07719949 & 0.06477605 & 0.06520783 \\
0.2 & 0.08433319 & 0.10307868 & 0.08443957 & 0.08484961 \\
0.25 & 0.10403539 & 0.12660871 & 0.10412393 & 0.10513928 \\ \multicolumn{3}{c}{Rel. Error} \\ \hline \(\sigma\) & PDE(w/ abs) & PDE (w/o abs) & MC \\ \hline
0.05 & 1.74e-01 & -3.76e-03 & -2.11e-02 \\
0.1 & 2.31e-01 & 3.91e-03 & -1.39e-02 \\
0.15 & 1.95e-01 & -2.81e-03 & -9.49e-03 \\
0.2 & 2.22e-01 & -1.26e-03 & -6.12e-03 \\
0.25 & 2.17e-01 & -8.51e-04 & -1.06e-02 \\ \end{tabular}
\end{table}
Table 1. A single run for \(d=5\) independent underlyings, where European Calls are compared to the price obtained through PDE-RWNN (_with_ and _without_ absorption) and the Monte Carlo methods along each dimension. Below, the relative errors of both methods are given. The MC method was run using the same paths as in the PDE-RWNN.

Figure 1. Empirical convergence of MSE under Black-Scholes in terms of the number of hidden nodes. Error bars mark 0.1 and 0.9 quantiles of 20 separate runs of the algorithm. The slope coefficient of the dashed line is obtained through regression of the means of individual runs, while the solid line represents \(1/K\) convergence and is shown as a reference.
#### 6.1.3. Basket option

We consider an equally weighted basket Call option with a payoff

\[g_{\mathrm{basket}}(\mathbf{S}_{T})\coloneqq\left(\frac{1}{d}\sum_{j=1}^{d}S_{T}^{ j}-\mathscr{K}\right)^{+},\]

where \(\mathscr{K}>0\) denotes the strike price. For simplicity, we consider an ATM option with \(\mathscr{K}\coloneqq\frac{1}{d}\sum_{j=1}^{d}S_{0}^{j}\) and set all \(S_{0}^{j}=1\) for \(j\in\{1,\dots,5\}\). The volatilities \(\sigma_{j}\) are uniformly spaced between \([0.05,0.25]\) and the correlation matrix is randomly chosen as

\[\mathbf{\rho}\coloneqq\begin{bmatrix}1&0.84&-0.51&-0.70&0.15\\ 0.84&1&-0.66&-0.85&0.41\\ -0.51&-0.66&1&0.55&-0.82\\ -0.70&-0.85&0.55&1&-0.51\\ 0.15&0.41&-0.82&-0.51&1\end{bmatrix},\]

so that \(\Sigma\coloneqq\operatorname{diag}(\mathbf{\sigma})\mathbf{\rho}\operatorname{diag}( \mathbf{\sigma})\). Since the distribution of a sum of Lognormal is not known explicitly, no closed-form expression is available for the option's price. Hence, the reference price is computed using Monte-Carlo with \(100\) time steps and \(400,000\) samples. In Table 3, we compare our scheme with a classical MC estimator in terms of relative error for \(K=100\) hidden nodes.

### Rough Bergomi

The rough Bergomi model belongs to the recently developed class of rough stochastic volatility models, first proposed in [4, 28, 35], where the instantaneous

\begin{table}
\begin{tabular}{l|c c} \(d\) & Total MSE (with abs) & CPU Time (s) \\ \hline
5 & 3.4826e-08 & 10.5 \\
10 & 5.4169e-08 & 16.0 \\
25 & 4.9015e-08 & 34.5 \\
50 & 1.6533e-07 & 65.0 \\
100 & 2.5340e-07 & 135.0 \\ \end{tabular}
\end{table}
Table 2. Total MSE calculated across all \(d\) assets and CPU training times for varying dimension \(d\), where \(\mathbf{\sigma}\) uniformly spaced over \([0.05,0.4]\).

Figure 2. Computational time vs number of dimensions, as in Table 2.
variance is driven by a fractional Brownian motion (or more generally a continuous Gaussian process) with Hurst parameter \(H<\frac{1}{2}\). As highlighted in many papers, they are able to capture many features of (Equities, Commodities,...) data, and clearly seem to outperform most classical models, with fewer parameters. Precise examples with real data can be found in [4] for SPX options, in [12, 29, 45, 39] for joint SPX-VIX options and in [8, 28] for estimation on historical time series, the latter being the state-of-the-art statistical analysis under the \(\mathbb{P}\)-measure. We consider here the price dynamics under \(\mathbb{Q}\) with constant initial forward variance curve \(\xi_{0}(t)>0\) for all \(t\in[0,T]\):

\[\left\{\begin{array}{rl}\frac{\mathrm{d}S_{t}}{S_{t}}&=r\mathrm{d}t+\sqrt{V_ {t}}\mathrm{d}\left(\rho_{1}\mathrm{d}W_{t}^{1}+\rho_{2}W_{t}^{2}\right),\\ V_{t}&=\xi_{0}(t)\mathcal{E}\left(\eta\sqrt{2H}\int_{0}^{t}(t-u)^{H-\frac{1}{2 }}\mathrm{d}W_{u}^{1}\right),\end{array}\right.\]

where \(\eta>0\) and \(H\in(0,1)\) is the Hurst parameter. The corresponding BSPDE reads

\[-\mathrm{d}u(t,x)=\left[\frac{V_{t}}{2}\partial_{x}^{2}u(t,x)+\rho\sqrt{V_{t}} \partial_{x}\psi(t,x)-\frac{V_{t}}{2}\partial_{x}u(t,x)-ru(t,x)\right]\mathrm{ d}t-\psi(t,x)\mathrm{d}W_{t}^{1},\]

with terminal condition \(u(T,x)=g_{\mathrm{call}}\left(\mathrm{e}^{x+rT}\right)\). While the existence of the solution was only proven in the distributional sense [5], we nevertheless apply our RWNN scheme. To test Algorithm 2, both the price and the volatility processes are discretised according to the Hybrid scheme developed in [9, 52]. We set the rBergomi parameters as follows: \(H=0.3\), \(\eta=1.9\), \(\rho=-0.7\), \(r=0.01\), \(T=1\), \(S_{0}=1\) and choose the forward variance curve to be flat with \(\xi_{0}(t)=0.235\times 0.235\). Again, we are pricing an ATM vanilla Call option with \(\mathscr{K}=S_{0}=1\). The number of discretisation steps is again \(N=21\), the number of Monte-Carlo samples \(n_{\mathrm{MC}}=50,000\) and the reservoir has \(K\in\{10,100,1000\}\) nodes with connectivity \(c=0.5\) in Section 6.2.2.

#### 6.2.1. Convergence rate

As in Section 6.1.1, we conduct an empirical analysis of the convergence error for the ATM Call. To isolate the dependence on the number of nodes we fix \(c=1\), \(n_{\mathrm{MC}}=50,000\) and vary \(K\in\{10,100,1000\}\). The reference price is computed by Monte-Carlo with 100 time steps and \(800,000\) samples. The absorption scheme has been applied and the results are displayed in Figure 3. In this section, the same random seed was used as in Section 6.1.1, to ensure consistent results across different simulations.

#### 6.2.2. ATM Call option

We now evaluate the performance of our PDE-RWNN method for option pricing in the rough Bergomi model using the market parameters listed above and compare the results to those obtained with the MC method over the same sample paths. We also investigate the effect of the absorption scheme (Table 4) and find that,

\begin{table}
\begin{tabular}{l|c c c c}  & Reference & PDE (w/ abs) & PDE (w/o abs) & MC \\ \hline Price & 0.016240 & 0.018220 & 0.016131 & 0.016251 \\ \hline Rel. error & - & -1.22e-01 & -6.71e-03 & -6.50e-04 \\ \hline Time & 12.8 & 9.7s & 9.8s & 0.3s \\ \end{tabular}
\end{table}
Table 3. Comparison of prices, relative errors and CPU time of the Monte-Carlo estimator, PDE-RWNN scheme _with_ and _without_ absorption (using same sampled MC paths and \(K=100\)) and the reference price computed with 100 time steps and 400,000 samples.
interestingly, despite keeping the paths positive, the absorption scheme adds noticeable bias. Nevertheless, the relative error of the proposed scheme with absorption is comparable to the results using regular artificial networks found in the literature [5, Table 1], yet, our scheme learns much faster with orders of magnitudes faster training times.

\begin{table}
\begin{tabular}{l|c|c|c|c}  & Reference & PDE (w/ abs) & PDE (w/o abs) & MC \\ \hline Price & 0.079932 & 0.0819236 & 0.079729 & 0.080310 \\ \hline Rel. error & - & 24.9e-3 & 2.54e-03 & -4.73e-03 \\ \hline Time & 10.1s & 7.4s & 7.5 & 0.4s \\ \end{tabular}
\end{table}
Table 4. Comparison of prices, relative errors and CPU time of the Monte-Carlo estimator, PDE-RWNN scheme _with_ absorption, PDE-RWNN scheme _without_ absorption both with \(K=100\) (and both using same sampled MC paths) and the reference price computed with 100 time steps and 800,000 samples.

Figure 3. Empirical convergence of MSE under rBergomi in terms of the number of hidden nodes. Error bars mark 0.1 and 0.9 quantiles of 20 separate runs of the algorithm. The slope coefficient of the dashed line is obtained through regression of the means of individual runs, while the solid line represents \(1/K\) convergence and is shown as a reference.