Recurrent Neural Networks with more flexible memory: better predictions
  than rough volatility
We extend recurrent neural networks to include several flexible timescales
for each dimension of their output, which mechanically improves their abilities
to account for processes with long memory or with highly disparate time scales.
We compare the ability of vanilla and extended long short term memory networks
(LSTMs) to predict asset price volatility, known to have a long memory.
Generally, the number of epochs needed to train extended LSTMs is divided by
two, while the variation of validation and test losses among models with the
same hyperparameters is much smaller. We also show that the model with the
smallest validation loss systemically outperforms rough volatility predictions
by about 20% when trained and tested on a dataset with multiple time series.
Damien Challet, Vincent Ragel
# Recurrent neural networks with more flexible memory: better predictions than rough volatility

###### Abstract

We extend recurrent neural networks to include several flexible timescales for each dimension of their output, which mechanically improves their abilities to account for processes with long memory or with highly disparate time scales. We compare the ability of vanilla and extended long short term memory networks (LSTMs) to predict asset price volatility, known to have a long memory. Generally, the number of epochs needed to train extended LSTMs is divided by two, while the variation of validation and test losses among models with the same hyperparameters is much smaller. We also show that the model with the smallest validation loss systemically outperforms rough volatility predictions by about 20% when trained and tested on a dataset with multiple time series.

Time series Long memory Recurrent Neural Networks Rough Volatility Volatility modelling

## 1 Introduction

Some time series in Nature have a very long memory (Robinson, 2003): fluid turbulence (Resagk et al., 2006), asset price volatility (Cont, 2001) and tick-by-tick events in financial markets (Challet and Stinchcombe, 2001; Lillo and Farmer, 2004). From a modelling point of view, this means that the current value of an observable of interest depends on the past by a convolution of itself with a long-tailed kernel.

Deep learning tackles past dependence in time series with recurrent neural networks (RNNs). These networks are in essence moving averages of nonlinear functions of the inputs and learn the parameters of these averages and functions. Provided that they are sufficiently large, these networks can approximate long-tailed kernels in a satisfactory way, and are of course able to account for more complex problems than a simple linear convolution. Yet, their flexibility may prevent them to learn quickly and efficiently the long memory of time series. Several solutions exist: either one pre-filters the data by computing statistics at with various time scales and use them as inputs to RNNs in the same spirit as multi-time scale volatility modelling (Zumbach and Lynch, 2001; Corsi, 2009), see e.g. Kim and Won (2018), or one extends the neural networks so as to improve their abilities. For example, Zhao et al. (2020) adds delay operators, taking inspiration from the ARIMA processes, to the states of recurrent neural networks, while Ohno and Kumagai (2021) modifies the update equation of the network output so that its dynamics mimics that of a variable with a long memory. In both cases, the time dependence structure is enforced by hand in the dynamics of such networks.

Here, we propose a flexible and parsimonious way to extend the long-memory abilities of recurrent neural networks by using an old trick for approximating long-memory kernels with exponential functions, which helps recurrent neural networks learn faster and better time series with long memory.

Our main contributions are: (i) we introduce RNNs with several multiple flexible time scales for each dimension of the output; (ii) we show that learning to predict time series with long-memory (asset price volatility) is faster and more reliably good with more flexible time scales (iii) rough volatility predictions can be beaten by training a fair number of recurrent neural networks and only using the one with the best validation loss.
Recurrent neural networks with more flexible memory: better predictions than rough volatility

## 2 Methods

Let time series \(y_{t}\) be of interest. Its moving average can be written

\[\tilde{y}_{t}=\int_{-\infty}^{t}K(t-t^{\prime})y_{t}^{\prime}dt^{\prime}, \tag{1}\]

where \(K\) is a kernel. In a discrete time context,

\[\tilde{y}_{t}=\sum_{-\infty}^{t}K(t-t^{\prime})y_{t}^{\prime}. \tag{2}\]

When the process is Markovian, its kernel \(K(x)\simeq e^{-x/\tau_{0}}\) for large \(x\), where \(\tau_{0}\) is the slowest timescale at which the process forgets its past. In this case, one can write \(y_{t}\) in a recursive way

\[\tilde{y}_{t}=\tilde{y}_{t-1}(1-\lambda)+\lambda y_{t}, \tag{3}\]

where \(\lambda\simeq 1/\tau_{0}\); \(\tilde{y}_{t}\) is then an exponentially moving average (EMA) of \(y_{t}\).

Long memory processes however, have a kernel that decreases at least as slowly as a power-law. In turn, power-laws can be approximated by a sum of exponential functions: naively, if \(K(x)=x^{-\alpha}\), one writes

\[K(x)\propto\sum_{i=1}^{\infty}w_{i}\exp(-x/\tau_{i}) \tag{4}\]

with \(w_{i}\propto(1/c^{\alpha})^{i}\) and \(\tau_{i}=c^{i}\) for a well-chosen constant \(c\): one covers the \(x\) space in a geometric way and the weights \(w_{i}\) account for the power-law decreasing nature of \(K(x)\). This rough approach works well and is widespread. Bochud and Challet (2007) propose a more refined method to determine how many exponential functions one needs to approximate optimally \(K\) and how to compute \(w_{i}\) for a given \(\alpha\) and for a given range of \(x\) over which the kernel has to be approximated by a sum of exponential functions (e.g. \(x\in[1,1000]\)). For example, one needs about 4 exponential functions to approximate 3 decades).

Writing down the update equations of well-known recurrent neural network architectures makes it clear that they use exponentially moving averages with a single time scale for each output dimension. For example, Gate Recurrent Units (GRU) Cho et al. (2014) transform the input vector \(x_{t}\) into a vector of timescales \(\lambda_{t}\) defined as

\[\lambda_{t}=\sigma(W_{\lambda}x_{t}+U_{\lambda}c_{t-1}+b_{\lambda}) \tag{5}\]

which is then used in the update of the output \(c_{t}\)

\[c_{t}=c_{t-1}\odot(1-\lambda)+\lambda\odot\tilde{c}_{t}, \tag{6}\]

where the update \(\tilde{c}_{t}\) is also computed from the input with learned weights, i.e.

\[\tilde{c}_{t} =\sigma_{c}(W_{c}x_{t}+U_{c}(c_{t-1}\odot r_{t})+b_{c}) \tag{7}\] \[r_{t} =\sigma_{r}(W_{r}x_{t}+U_{r}c_{t-1}+b_{r}) \tag{8}\]

for a non-linear functions \(\sigma_{c}\) and \(\sigma_{r_{t}}\odot\) is the element-wise (Hadamard) product and \(r_{t}\) the reset gate which modifies the value of \(c_{t-1}\) when computing \(\tilde{c}_{t}\). By design, GRUs can only compute exponentially moving averages of \(\tilde{c}_{t}\), although they possess the interesting ability to learn both \(\lambda_{t}\) and the update \(\tilde{c}_{t}\) as a function of their inputs. It is straightforward to extend GRUs to an arbitrary number of timescales \(n\) by using \(n\)\(\tilde{c}_{t}^{(k)}\), \(k=1,\cdots,n\) and

\[\lambda_{t}^{(k)} =\sigma(W_{\lambda}^{(k)}x_{t}+U_{\lambda}^{(k)}c_{t-1}^{(k)}+b_{ \lambda}^{(k)}) \tag{9}\] \[c_{t}^{(k)} =c_{t-1}^{(k)}\odot(1-\lambda_{t}^{(k)})+\lambda_{t}^{(k)}\tilde {c}_{t}, \tag{10}\]

where each \(c_{t}^{(k)}\) is an exponentially moving average at time scale \(\sim 1/\lambda_{t}^{(k)}\). Finally, the output will be

\[c_{t}=\sum_{k=1}^{n}w_{k}c_{t}^{(k)}, \tag{11}\]
The simple \(\alpha-\)RNN (Dixon and London, 2021), which are simplified GRUs, share the same assumption of a single time scale per output dimension and thus can be generalized in the same way. Let us show now how extend LSTMs with a forget gate (Gers et al., 2000). Starting from their output \(h_{t}\), one has

\[h_{t} =o_{t}\odot\sigma_{h}(c_{t}) \tag{12}\] \[c_{t} =f_{t}\odot c_{t-1}+i_{t}\odot\tilde{c}_{t}, \tag{13}\]

where \(o_{t}\), \(i_{t}\), and \(\tilde{c}_{t}\) are determined from the input \(x_{t}\) and the previous output \(h_{t-1}\) with learned weights and \(\sigma_{h}\) is a nonlinear function. Writing \(f_{t}=1-\lambda_{t}\) makes it obvious that the cell vector \(c_{t}\) evolves in the same way as \(y_{t}\) in Eq. (6) if \(i_{t}\simeq\lambda_{t}\).

Extending LSTMs to include \(n\) time scales by cell state dimension is therefore straightforward: one needs to compute \(n\) EMAs and their associated \(\lambda\)s as follows

\[f_{t}^{(k)} =\sigma(W_{f}^{(k)}x_{t}+U_{f}^{(k)}h_{t-1}+b_{f}^{(k)}) \tag{14}\] \[i_{t}^{(k)} =\sigma(W_{i}^{(k)}x_{t}+U_{i}^{(k)}h_{t-1}+b_{i}^{(k)})\] (15) \[c_{t}^{(k)} =f_{t}^{(k)}\odot c_{t-1}^{(k)}+i_{t}^{(k)}\odot\tilde{c}_{t}, \tag{16}\]

where \(\tilde{c}_{t}\) follows Eq. (11). Note that one could set \(i_{t}^{(k)}=1-f_{t}^{(k)}=\lambda_{t}^{(k)}\) and not learn the weights associated to \(i_{t}\). Learning as well \(i^{(k)}\) is equivalent to modulate the importance of the update, which is known as to as cognitive bias (Palmitneri et al., 2017): this is made clear by writing \(i_{t}^{(k)}=v_{t}^{(k)}(1-f_{t}^{(k)})=v_{t}^{(k)}\lambda_{t}^{(k)}\), where \(v_{t}^{(k)}\) is the modulation of learning speed.

We will focus on the case \(n=2\): (11) amounts to

\[c_{t}=c_{t}^{(1)}\odot\alpha+(1-\alpha)\odot c_{t}^{(2)}, \tag{17}\]

where the vector \(\alpha\) is learnable and its components are bounded to the \([0,1]\) interval. We call LSTMs with several time scales (\(n>1\)) per dimension VLSTMs, which stands for very long short term memory.

Note that LSTMs with a sufficiently large cell dimension \(N_{h}\) can in principle learn to superpose timescales in the same way as Eq. (11) by learning one time scale per dimension and using final dense layer to learn how to combine them. However, imposing constraints (or equivalently, injecting some known structure) is known to lead to faster learning and better results (e.g. Physics-guided deep learning, see Thuerey et al. (2021) for a review).

Naively, when learning to predict a process that is not too noisy, we expect the difference between VLSTM and LSTM to be the highest when \(N_{h}=1\), i.e. precisely when LSTMs do not have the possibility to compute long-term averages and to decrease when \(N_{h}\) increases.

### Volatility prediction

Given an asset price \(P_{t}\), its and its log return \(r_{t}=\log P_{t}-\log P_{t-1}\), the asset price volatility \(\sigma\) is defined as \(\sigma^{2}=E(r^{2})\). The dynamics of financial markets is ever-changing, which results in a temporal dependence of \(\sigma\) with clear patterns of long-term dependence (Cont, 2001) (see Fig. 1).

Risk management, portfolio optimization. and option pricing benefit from the ability to predict the evolution of \(\sigma_{t}\). Fortunately, \(\sigma_{t}\) is relatively easy to predict, owing to its long memory (Cont, 2001): for example, its auto-correlation decreases very slowly, presumably as a power-law over more than a year. Econometric models include GARCH, whose simplest version involves only one timescale, while many variations use several time scales (Zumbach and Lynch, 2001; Corsi, 2009; Zumbach, 2015). Rough volatility (Gatheral et al., 2018), on the other hand, considers \(\log\sigma_{t}\) as fractional Brownian motion and thus includes all the time scales. As can be expected, rough volatility models outperform GARCH-based models for volatility prediction. Using LSTMs for volatility prediction is found e.g. in Kim and Won (2018); Filipovic and Khalilzadeh (2021); Rosenbaum and Zhang (2022) that use various types of predictors (including GARCH models) and architectures. Notably, Rosenbaum and Zhang (2022) show that the average prediction of 10 stacked LSTMs with past volatility and price return as predictors match the performance of rough volatility.

### Architecture and hyperparameters

Our first aim is to characterise the effects of multiple time scales per cell dimension. Therefore, we compare simple non-stacked LSTMs with or without the proposed modification. Stacked LSTMs can learn additional time scales at
the cost of doubling the number of parameters, which we precisely wish to avoid here. We pass the outputs \(h_{t}\) of the LSTMs and VLSTMs through a dense layer of size \(N_{h}\) with sigmoid activation functions, so as to combine the outputs in a non-trivial way, and a final dense layer with linear activation. Both final layers have a bias term, which allows the model to learn a baseline volatility level.

We report a systematic study of the relative performance of LSTMs vs VLSTMs. We vary the sequence length \(T_{\text{seq}}\) from 10 to 100 by steps of 15, and the dimension of the hidden state \(N_{h}\in\{1,\cdots,5\}\). Finally, we train models with and without biases (except for the final two dense layers which always have biases). There are thus 70 variations of hyperparameters per architecture choice.

For each hyperparameter and architecture couple, we train 20 networks which yields 2800 models altogether. We use a standard 60/20/20 train/validation/test splits and apply early stopping criterion of the minimum validation loss over the 5 last epochs, with a maximum of 1000 epochs. Batch size is set to 128. We train the networks to predict \(\log\sigma_{t+1}\) with an MSE loss function.

Data comes from volatility computed by Oxford-Man Institute from intraday data with the two-scale realized kernel estimation method (Barndorff-Nielsen et al., 2008), which contain volatility time series for 31 indices and 2117 to 5385 data points per index (Heber et al., 2022). Since the volatility individual time series start and end at heterogeneous dates, we used the dates to define the train/validation/test splits: the train set ranges from 2000-01-04 to 2012-09-06, validation set from 2012-09-07 to 2016-11-23 and test set from 2016-11-24 to 2021-02-17. This is necessary as the time series are cross-correlated, hence, splitting them according to their respective length would cause information leakage from the future and thus overfitting.

## 3 Results

### Average loss

Let us plot the average test loss of LSTMs and VLSTMs as a function of \(N_{h}\) at fixed \(T_{\text{seq}}\), the dimension of the memory cell, and of \(T_{seq}\) at fixed \(N_{h}\). This approach is taken by Rosenbaum and Zhang (2022) who trained 10 LSTMs instead of 20 here. Figure 2 shows that VLSTMs enjoy a sizeable advantage on average. We note that when \(N_{h}=1\), our initial intuition was correct: VLSTMs have a smaller average test loss for all variations of hyperparameters (\(T_{\text{seq}}\) and bias)

Large loss fluctuations among models are associated with large average test losses for both VLSTMs and LSTMs; however test losses of VLSTMs are more likely to be small (and have accordingly small fluctuations). This is explained by a large difference in training convergence time, as shown below. We also note that, at least for volatility prediction, keeping bias terms in the computation of \(i\), \(f\), \(\bar{c}\), and \(c\) (referred to as internal biases henceforth) is manifestly problematic; it turns out to be the default option both for PyTorch and Keras and is probably implicitly used in other papers. On the whole, we note that a simple average of the outputs of an ensemble of models leads to quite large fluctuations, hence that the question of the convergence of the models must be investigated and a way to select the good models would much improve the usefulness of LSTMs in that context.

Figure 1: Volatility of the SPX index (annualized) as a function of time, showing traces of long-term dependence. Data source: Oxford-Man Institute (Heber et al., 2022).
Train convergence, it turns out, is a hit and miss process: some models are stuck in a high loss regime, while some models do learn a more realistic dynamical process and reach much lower losses. This yields a bi-modal density of losses (see Fig. 3). It is noteworthy that the fraction of VLSTMs that learn better is much larger. This is linked to the fact that VLSTMs learn much faster (see below) and that VLSTMs without internal biases are less likely to be stuck in a high loss regime.

Since the volatility process is well approximated e.g. by a rough volatility model (Gatheral et al., 2018), the test loss is commensurate with the validation loss as expected, itself commensurate with the train loss. We plot in Fig. 3 the test loss versus the validation loss, which shows that test losses are accordingly also bimodal, with a majority of models not stuck in the high loss regime, some having a test loss smaller than rough volatility models.

### Keeping the better models

This result suggests a way to select the good models, since the validation loss distribution is bimodal and since the test losses are roughly proportional to validation losses. To select models whose validation loss belongs in the lower peak, we compute 9 quantiles \(q(p)\) with regular sequence of probabilities \(p=0.1,\cdots,0.9\), and keep the models whose validation loss is smaller than the quantile corresponding to the maximum change between quantiles, a simple yet effective way to find well separated peaks. We call these models the better ones in the following. This procedure allows a fairer comparison between LSTMs and VLSTMs.

VLSTMs are still better than LSTMs, even for larger \(N_{h}\). Figure 4 plots the average test loss of the models with below-average validation loss versus \(N_{h}\) and the sequence length. The test losses are now much closer, but VLSTMs still retain a sizable advantage: their test losses are both lower on average and their fluctuations are much smaller.

Both the variability of results and the strange results for \(N_{h}=1\) when biases are allowed in the computation of the internal states of (V)LSTMs can be traced back to training convergence problems. A simple way to ascertain the main

Figure 2: Volatility prediction. Upper plots: mean test loss vs the memory cell dimension \(N_{h}\) (\(T_{\text{seq}}=40\)); lower plots: mean test loss vs the sequence length \(T_{\text{seq}}\) (\(N_{h}=2\)). Left plots: (V)LSTMs with bias weights; right plots: (V)LSTMs with no bias weights. The dashed line is the average MSE of predictions made with rough volatility models.
Recurrent neural networks with more flexible memory: better predictions than rough volatility

Figure 4: Multiple volatility time series prediction test losses of the models with below-average validation losses. Upper plots: mean test loss vs the memory cell dimension \(N_{h}\) (\(T_{\text{seq}}=100\)); lower plots: mean test loss vs the sequence length \(T_{\text{seq}}\) (\(N_{h}=3\)). Left plots: (V)LSTMs with bias weights; right plots: (V)LSTMs with no bias weights. The dashed line is the average MSE of predictions made with rough volatility models.

Figure 3: Left plot: density of validation losses by architecture. Right plot: test loss vs validation loss. Multiple volatility time series prediction. The dashed line is the average MSE of predictions made with rough volatility models.
difference between VLSTMs and LSTMs is to measure the time it takes for their training to converge, i.e., to reach the early stopping criterion. Figure 5 reports the fraction of models that have converged as a function of the number of epochs (limited to 1000). LSTMs need more epochs to be trained. We also found that the case \(N_{h}=1\) and small \(T_{\text{seq}}\) is hard to learn for this kind of architecture, the training of many models requiring more than 1000 epochs to reach the early stopping criterion.

### Best model

Finally, let us investigate the test loss of the model with the best validation loss among the 20 models trained for each of the 140 hyperparameters/architecture choices. It turns out that under these conditions VLSTMs and LSTMs have essentially the same performance. What differentiate them however is the speed at which they learn. Let us plot the test loss versus the time of convergence for LSTMs and VLSTMs with and without biases (left plot of Fig. 6): there is slight negative dependence between test loss and convergence times, the longer one learns, the better. Notably convergence times of LSTMs are spread all over the whole \([1,1000]\) interval, while VLSTMs converge before 400 epochs. The right plot of Fig. 6 displays the ECDF of the convergence times, which shows a sizable difference between LSTMs and VLSTMs: whereas 20% of LSTMs models do not manage to converge before 1000 epochs, all VLSTMs do before 400, except one, when biases are allowed.

Thus, training a given number of models is significantly shorter with VLSTMs because they do not need to learn how to approximate the kernel \(K(x)\). One also sees that models with internal biases converge more slowly than those without them. We also wish to point out that because the fluctuation of validation losses among the trained models is much smaller for VLSTMs than for LSTMs, hence, that in practice, one needs to train fewer VLSTMs than LSTMs before finding a good one.

## 4 Conclusion

Adding an explicit but flexible kernel structure to LSTMs brings significant improvements in every metric: number of epochs needed to reach convergence, overall prediction accuracy, and accuracy variation between models at fixed

\begin{table}
\begin{tabular}{l c c c} Architecture & bias & test loss & test loss \\  & & average & std dev. \\ \hline rough vol. & & 0.288 & \\ LSTM & yes & 0.241 & 0.032 \\ LSTM & no & 0.245 & 0.057 \\ VLSTM & yes & 0.232 & 0.017 \\ VLSTM & no & 0.230 & 0.015 \\ \end{tabular}
\end{table}
Table 1: Standard deviation of the test losses of the better models computed over all the values of \(N_{h}\) and \(T_{\text{seq}}\). Multiple time series volatility prediction

Figure 5: Fraction of models having converged before a given number of epochs. Left plot: \(N_{h}=1\), right plot: all \(N_{h}\); multiple time series volatility prediction, \(T_{\text{seq}}=100\).
parameters. There is a cost as the number of trainable parameters of VLSTMs is larger than LSTMs are fixed hyperparameters, but doubling the number of time scales does not require to double the number of trainable parameters, thanks to the explicit kernel approximation structure. Although this paper focuses on LSTMs, the same idea can be applied to GRUs and \(\alpha-\)RNNs in a straightforward way.

Our results mirror those of Rosenbaum and Zhang (2022): we also succeeded in training a single model at a time on many volatility time series of various underlying asset types. This reflects the universality of volatility dynamics, a fact also hinted at by rough volatility and multi-scale GARCH-like models (Zumbach, 2015).

While it is hard to beat rough volatility for volatility prediction, we found that even simple LSTMs can beat it, provided that one trains several models and selects the best one according to its validation loss. Using LSTMs for that purpose requires to train more models over more epochs than using VLSTMs. Volatility prediction can be further improved by adding some more features, such as prior knowledge of predictable special events, and possibly by using more complex neural architectures.

## 5 Code and data availability

Full code, including the Keras VLSTM class, and data, are available at [https://github.com/damienchallet/VLSTM](https://github.com/damienchallet/VLSTM).

## Acknowledgments

This work used HPC resources from the "Mesocentre" computing center of CentraleSupelec and Ecole Normale Superieure Paris-Saclay supported by CNRS and Region Ile-de-France.
