MOrdReD: Memory-based Ordinal Regression Deep Neural Networks for Time Series Forecasting
Bernardo Pérez Orozco Gabriele Abbati Stephen Roberts




Corresponding author: B. Pérez Orozco 
Information Engineering, Dept. of Engineering Science, University of Oxford, UK

Tel: +44 (0)1865 616600

ORCID: 0000-0001-9741-7950

ber@robots.ox.ac.uk

 
G. Abbati Information Engineering, Dept. of Engineering Science, University of Oxford,  UK
gabb@robots.ox.ac.uk


S.J. Roberts Information Engineering, Dept. of Engineering Science & Oxford-Man Institute of Quantitative Finance, University of Oxford, UK

& Mind Foundry Ltd.
sjrob@robots.ox.ac.uk


MOrdReD: Memory-based Ordinal Regression Deep Neural Networks for Time Series Forecasting
    Bernardo Pérez Orozco Gabriele Abbati Stephen J. Roberts
    Received: 24 October 2018 / Accepted: TBD
=========================================================================================




Time series forecasting is ubiquitous in the modern world. Applications range from health care, astronomy and include climate modelling, to financial trading and monitoring of critical engineering equipment. To offer robust value over this wide range of activities, models must not only provide accurate forecasts, but also quantify and adjust their uncertainty over time. In this work, we directly tackle both tasks with a novel, fully end-to-end deep learning method. By recasting time series prediction as an ordinal regression task, we develop a principled methodology to assess long-term predictive uncertainty and describe the rich, multi-modal, non-Gaussian behaviour which arises regularly in many problem domains. 

Notably, our framework is a wholly general-purpose approach that requires little to no user intervention to be used. We showcase this key feature in a large-scale benchmark test with 45 datasets drawn from both a wide range of real-world applications and a comprehensive list of synthetic data. This wide comparison uses, as benchmark methods, state-of-the-art models from both Machine Learning and Statistics literature, such as the Gaussian Process. We find that our approach not only provides excellent predictive forecasts and associated uncertainty bounds, but also allows us to infer derived information, such as the predictive distribution of critical events of interest. We show that the method we propose performs accurately and reliably, even over long time horizons. 
 


    
	
	

§ INTRODUCTION








Performing long-term, reliable time series forecasting has wide interest in both academia and industry. Application domains are great in number, including meteorology, energy systems, astronomy, finance, dynamical systems, physiology and many others. Crucially, forecasts about a system's future should not be restricted to merely investigating the most likely, or expected, values of the system at a specific time instant. Information regarding the likely timing of an event of interest is important, and so are beliefs regarding the value at a given time. Such questions can naturally be answered by providing full predictive probability distributions at each time step. 



In the main, forming a rich probability distribution over future states has been addressed by statistical models. The latter are often parametric, such as polynomial regressors or the family of Markovian state-space models. More recently, non-parametric approaches such as Gaussian Processes (GPs) <cit.> allow avoidance of pre-specified parametric forms. In part, no doubt, due to their native ability to infer full predictive distributions, GPs (and related models) have increasingly dominated the time series modelling literature over the last decade, especially as  classic Markov models (including autoregressive processes and the like) can be seen as special cases of the GP. We simply note that, powerful though the GP approach is, the standard GP model makes the tacit assumption that the residuals are normally distributed and hence the conditional predictive posterior is also Gaussian. This could be restrictive in the long-term forecasting scenario, where short-term forecast uncertainty may induce distant scenarios whose nature is inherently multi-modal. Furthermore, in its naive form, the GP also suffers from poor scalability, which can hinder the use of the large amounts of data readily available in some tasks. Lastly, to obtain best performance, expert knowledge, or extensive searches may be needed to refine the kernel function at the heart of the GP.

In parallel to these developments, deep and recurrent neural network models have gained significant attention in the domain of classification tasks, such as those that arise in the Natural Language Processing community. Neural networks, despite their complexity, scale well with dataset size and the increasing availability and development of algorithms, data, hardware and software have enabled these models to become the state-of-the-art in sequential classification tasks such as speech recognition, machine translation, scene labelling, and others. In the context of regression tasks such as time series forecasting, however, these models have not achieved the same success. Even though time series forecasting with neural networks has been investigated extensively in the literature, to the authors' best knowledge, this work only deals with point forecasts. A potential reason is that fully probabilistic models such as Bayesian neural networks are analytically intractable, and approximate Bayesian inference methods can be computationally expensive even for networks with relatively few units. This leaves little room to justify their usage in comparison with probabilistic methods.

In this work, we develop a recurrent neural network framework for time series prediction that is capable of producing fully (though approximate) probabilistic forecasts. We achieve this by recasting the time series forecasting task as an ordinal (auto-)regression one, which consists of passing the time series through a quantisation step, thus transforming the prediction task into a sequential classification problem. This enables us to incorporate state-of-the-art developments from the Machine Learning and Natural Language Processing communities. Such ordinal forecasts can further describe rich multi-modal, non-Gaussian behaviour, akin to that which may describe long-term forecasts more faithfully. The number of modes and other shape features of the predictive posterior distribution can thence adapt dynamically for different predictive horizons.

Moreover, our ordinal framework is a wholly general-purpose approach that is scalable with dataset size and requires little to no human intervention or expert craftsmanship for model specification. We showcase this key feature by means of an exhaustive comparison over 45 datasets from a wide range application domains, encompassing a variety of state-of-the-art methods from both the Machine Learning and Statistical modelling communities. Crucially, we evaluate both the accuracy and calibration of our probabilistic forecasts, in addition to the deviation of point forecasts, such as the mean and the median of the predictive distributions. We provide evidence that our method achieves superior performance in the majority of datasets, and furthermore, we show that even when it does not, on average we still expect it to surpass most of its competitors.

To further showcase how our probabilistic forecasts can be used to address other queries of interest, we also evaluate our model in the context of event timing prediction. In particular, we demonstrate how to construct predictive densities for the occurrence of local maxima with Kernel Density Estimation, and show that our framework is able to accurately form beliefs over the timing of such events.


















In summary, the contributions of our work are as follows:

	
  * We develop a fully end-to-end framework for long-term time series forecasting that requires little to no manual model tuning and outperforms other state-of-the-art techniques;
	
  * we incorporate a sequential procedure to infer reliable long-term uncertainty bounds in an autoregressive fashion;
	
  * we perform an exhaustive, large-scale comparison study of this and other state-of-the-art forecasting techniques;
	
  * we demonstrate how our probabilistic forecasts can be employed to predict the occurrence of (e.g.) local maxima in a timely and reliable manner.


The remainder of this paper is organised as follows: in Section <ref> we offer a review of the relevant literature; then, in Section <ref> we present an overview of our ordinal autoregression framework, as well as a succinct review of the methods we used as benchmarks; in Section <ref> the experiments that support the paper and the data used are detailed. We lastly present our conclusions in Section <ref>.











































	
	

§ LITERATURE REVIEW
 



 §.§ MACHINE LEARNING METHODS FOR TIME SERIES FORECASTING
 

Machine Learning (ML) has a rich history in the time series literature, covering a broad range of methods such as random forest regressors <cit.>, quantile random forests <cit.> and artificial neural networks (ANNs) <cit.>. However, many of these do not readily provide probabilistic forecasts, and some others only do so in a limited fashion. 

Consequently, the time series literature has increasingly been dominated by models such as the Gaussian Process <cit.> which enable the inference of full predictive distributions in a principled manner, even in the presence of little data. Indeed, even summarised predictions, such as those given by the predictive posterior mean, have been shown to produce smaller deviations than those provided by competitors such as neural networks in a variety of application domains, for example in <cit.>. Their standard form, however, assumes that residuals distribute normally and thence the conditional predictive posterior is also a Gaussian. Furthermore, this approach (in its most general form) suffers from poor scalability, which can hinder the use of large datasets.

Despite not being able to produce probabilistic forecasts by design, ANNs are popular ML models that are known to scale well with dataset size. As a result, there is a rich literature covering applications of time series forecasting with feed-forward neural network models such as the multi-layer perceptron (MLP) and the radial basis function network (RBFNs), for example: electricity pricing <cit.>, stock pricing <cit.>, water quality <cit.> and solar cell energy output <cit.>. Some neural network architectures are also endowed with a memory mechanism, as is the case for the recurrent neural networks (RNNs). This mechanism refers to an encoding scheme in which the correlation between successive elements in sequential data is exploited to extract characteristic temporal features. RNNs however suffer from the vanishing and exploding gradient problems, which renders them more difficult to optimise by gradient descent than their memory-less counterparts. This has been addressed by introducing learnable gating mechanisms that allow for error signals to be propagated through time, thus enabling efficient learning via gradient descent, as in the case of the long short-term memory (LSTM) <cit.>. In consequence, LSTMs have shown themselves to offer excellent, end-to-end performance in predictions over symbol sequences, especially in the Natural Language Processing literature <cit.>.

Gating mechanisms, however, come at the expense of introducing more model parameters, which in turn may require more data to be inferred without overfitting. Indeed, such networks could only be trained at small scale until recently, and this could be a potential reason why early findings showed that they were easily outperformed by simpler, memory-less methods <cit.> in the context of time series. With the development of specialised hardware and software however, as well as increased availability of data for some tasks, recurrent models with a greater numbers of units have now been successfully used in several application domains, such as traffic flow prediction <cit.>, precipitation nowcasting <cit.>, chaotic dynamical systems <cit.>, electricity market price <cit.>, financial market predictions <cit.> and energy output forecasting <cit.>.



 §.§ COMPUTING PREDICTIVE UNCERTAINTY WITH RECURRENT NEURAL NETWORKS
 

In spite of its extent, the time series literature addressing the quantification and evaluation of predictive probability distributions with (recurrent or otherwise) neural networks remains scarce, with much being historic research from decades ago. Further, assessments are often constrained to benchmarking average deviation metrics, such as the mean squared or absolute errors (MSE, MAE) and the like. Some portion of the literature has attempted to quantify uncertainty, either in a limited fashion (for instance, by estimating upper and lower bounds <cit.>), or fully, by performing approximate or sample-based Bayesian inference to directly learn a Bayesian Neural Network <cit.>. The latter approach however introduces scalability issues that make it difficult to justify their usage for time series prediction in lieu of models such as the Gaussian Process.

Nonetheless, recent findings linking the dropout regularisation technique to approximate Bayesian inference <cit.> have enabled an efficient and principled, though still approximate, quantification of uncertainty in a variety of settings. Indeed, recent work looks to apply these approaches to time series forecasting <cit.>, however still making the tacit assumption that the predictive distribution is uni-modal Gaussian.

Dropout <cit.> is a stochastic regularisation technique that has been successfully employed to prevent overfitting and approximately combine different network architectures. At training time, before each feed-forward operation (necessary to compute the gradients of the loss function), each unit of the neural network is dropped with probability p_drop. In this way, a mask computed accordingly to this Bernoulli distribution is placed on the network, which in turn gets “thinned”. For gradient-based optimisers, gradients are computed with these new approximate outputs before updating the learnable parameters. A new mask is then drawn for each subsequent iteration. 

With this approach, a neural network with n units can be seen as a collection 2^n thinned networks sharing a substantial number of weights. As shown in <cit.>, a Monte Carlo approach can then be followed to approximate the predictive posterior distribution. This derives from the fact that, under certain assumptions, a neural network of arbitrary depth trained with dropout before each layer is mathematically equivalent to a probabilistic Deep Gaussian Process <cit.>.








	
	

§ METHODOLOGY



In this Section, we develop  the methods  we use for our large-scale benchmarking task. We start by introducing our own MOrdReD[In Arthurian legend, Mordred was the illegitimate son of King Arthur. Mordred fell in battle against his own father.] methodology: a Memory-endowed Ordinal Regression Deep neural network for time series forecasting. We then give a succinct summary of the baselines used to compare our model with, namely Gaussian Processes and the state-space formulation of autoregressive (AR) forecasting. 



 §.§ ORDINAL REGRESSION


We start by defining what an ordinal regression task is. Consider a bounded time series with range I ⊂ℝ and let C = {C_i}_i=1^M be a partition of I with cardinality M. Without loss of generality, we assume the C_i have all the same measure on ℝ (e.g. the C_i are non-overlapping sub-intervals of I with equal sizes). 


We can then define the time series forecasting problem in an ordinal and autoregressive fashion. Assume the P-sample-long observed sequence 𝐗^(t)∈ℝ^P × M, with 𝐗^(t) = (𝐱_1^(t), …, 𝐱_P^(t)) = (𝐱_t - P + 1, …, 𝐱_t), where P is the lookback window horizon and each 𝐱_i∈ℝ^M is the quantised one-hot encoded representation of a time series observation, originally assumed to be real-valued. The task of ordinal regression, which lies between regression and classification, consists in learning a map from the latter sequence to a symbol 𝐲^(t). Autoregressive time series forecasting can then be enabled by letting 𝐲^(t) = 𝐱^(t)_P+1=𝐱_t+1.

The procedure described above is indeed equivalent to framing the quantised one-hot encoded time series 𝐱  with a sliding window of length P. Given the symbol set cardinality M and lookback P, our time series datasets can be now defined as 𝒳 = {𝐗^(P), …, 𝐗^(N)} and 𝒴 = {𝐲^(P), …, 𝐲^(N)}.

We remark that such autoregression may induce a potential loss of accuracy during the quantisation phase. This can be mitigated by choosing a relatively large M that allows for a sufficiently fine-grained partition C, and we note that, within reason, the choice of M (problem dependent and a trade-off between resolution and computational simplicity) does not affect results. LSTMs (introduced in the next Section) thence naturally arise as a befitting methodology, due to their empirically proven ability to model categorical sequences over large symbol sets, such as those that frequently arise in the Natural Language Processing domain <cit.>. We further note that, within reason, quantisation has little impact in other queries of interest, such as long-term trend changes or event occurrence forecasting.



 §.§ MEMORY-ENDOWED ORDINAL REGRESSION DEEP NEURAL NETWORKS


In this Section our forecasting method is introduced. We consider a Neural Network parameterised by a vector of parameters θ, whose output is f^θ(𝐗). In particular, we consider LSTM-based architectures trained on a dataset 𝒳, 𝒴, with model likelihood p(𝐲|𝐗, θ).



  §.§.§ The Long Short Term Memory


The Long Short-Term Memory (LSTM) <cit.> is a particular type of recurrent neural network (RNN) endowed with a gating mechanism that enables efficient propagation of error signals in time. Such mechanisms prevent well-known gradient instability issues that arise in other simple recurrent neural network models <cit.>. Gating thus allows for efficient gradient-based learning of temporal features, which are encoded internally in memory cells that are updated with every new observation in the sequence. 

Consider an observed sequence 𝐗^(n) = (𝐱_1^(n), …, 𝐱_P^(n)), and let ⊙ be the element-wise product operator. For each t = 1,…,P, the LSTM outputs a vector of temporal features 𝐡_t, which is given by:

    𝐡_t    = 𝐨_t ⊙tanh(𝐂_t).

𝐂_t is the memory cell at time t:

    𝐂_t = 𝐢_t⊙𝐒_t + 𝐟_t⊙𝐂_t-1,

where

    𝐒_t = tanh(𝐖_𝐒𝐱'_t + 𝐛_𝐒)

and 𝐢_t, 𝐨_t, 𝐟_t are the input, output and forget gates respectively:

    𝐢_t =   σ(𝐖_𝐢𝐱'_t + 𝐛_𝐢),
    	
    𝐨_t =   σ(𝐖_𝐨𝐱'_t + 𝐛_𝐨),
    	
    𝐟_t =   σ(𝐖_𝐟𝐱'_t + 𝐛_𝐟).

Here σ is the element-wise sigmoid function and 𝐱'_t = (𝐱_t, 𝐡_t-1) is the concatenation of the observation 𝐱_t and the LSTM output at the previous time step 𝐡_t-1. In the ordinal regression setting, we optimise the categorical cross-entropy loss function with respect to the LSTM's learnable parameters  θ = {𝐖_𝐢, 𝐖_𝐨, 𝐖_𝐟, 𝐖_𝐒, 𝐛_𝐢, 𝐛_𝐨, 𝐛_𝐟, 𝐛_𝐒}.




  §.§.§ Sequence-to-Sequence


LSTMs can be further enriched by incorporating learnt features from the backwards temporal dynamics of the time series, in a bidirectional fashion. This is simply given by a second LSTM that scans the input sequence 𝐗 in reverse order, and whose outputs are then averaged with those of the initial LSTM. This bidirectional model has shown remarkable success in the NLP literature <cit.>, including extensions to LSTM-based models such as the sequence-to-sequence architecture <cit.>. 

We now focus our attention to the latter. The sequence-to-sequence model consists of two recurrent neural network models: an encoder f^(enc), which maps the observed sequence, 𝐗, into a fixed-dimensional summary, 𝐡^(dec)_0, 𝐂^(dec)_0; and a decoder f^(dec), which uses 𝐡^(dec)_0, 𝐂^(dec)_0 as an informed initial state to predict future observations iteratively. More precisely:

    (𝐡^(dec)_0, 𝐂^(dec)_0)    = f^(enc)(𝐗),
    𝐲̂_t    = Softmax( f^(dec)(𝐡^(dec)_t-1, 𝐂^(dec)_t-1, 𝐱_t))

We summarise this sequence-to-sequence model in Figure <ref>.





  §.§.§ MOrdReD: an Ordinal Regression Sequence-to-Sequence framework


In our ordinal (auto-)regression setting, a sequence-to-sequence's encoder firstly produces a summary 𝐡^(dec)_0, 𝐂^(dec)_0 of the last P observed samples of a one-hot encoded, quantised time series, 𝐗^(enc) = (𝐱_1, …, 𝐱_P). Then, the informed summary 𝐡^(dec)_0, 𝐂^(dec)_0 is fed into the model's decoder, which observes the last available sample 𝐱_P and finally outputs a class probability density 𝐲̂_P∈ℝ^M over the bins C, with 𝐱_P+1∼𝐲̂_P. 

We note that 𝐲̂_P has two possible interpretations: on the one hand, it is the categorical distribution that governs the behaviour of 𝐱_P+1; on the other hand, 𝐲̂_P itself can be interpreted as a corrupted or otherwise noisy representation of 𝐱_P+1. In our model, we enable direct autoregression by feeding back 𝐲̂_P into the decoder to forecast 𝐱_P+2∼𝐲̂_P+1. We have observed empirically that such autoregression allows our model to perform better than by one-hot encoding either a summary statistic or a sample drawn from 𝐲̂_P. We argue that this is the case because the model has full information about its own previous decisions, including partial bin allocations.

This loopback process can be repeated an arbitrary number of times to forecast over long horizons, as we show in Section <ref>. During training time, we enable teacher forcing by allowing the decoder to see the true, uncorrupted observation 𝐗_P+k, k>1. In other words, we jointly train the encoder and the decoder by allowing them to see 𝐗^(i), 𝐗^(i+1) respectively. At testing time however, the decoder only has access to at most 1 observed sample, which it uses as a seed to perform iterative one-step ahead prediction.






 §.§ QUANTIFYING UNCERTAINTY WITH MORDRED
  
As described in Section <ref>, Monte Carlo dropout can be used to approximate the predictive posterior distribution of a deep neural network. The predictive distribution of the model presented in Section <ref> is given by:

    p(𝐲^* | 𝐗^*) = ∫ p(𝐲^* | 𝐗^*, θ) p_drop(θ) dθ,

where p_drop(θ) denotes the dropout (Bernoulli) distribution acting as described in Section <ref>. It is then possible to approximate the integral via Monte Carlo:

    p(𝐲^* | 𝐗^*) ≈1N_s∑_n=1^N_sSoftmax( f^θ̂^(n)(𝐗^*) ),

where θ̂^(n)∼ p_drop(θ).
Thanks to the theoretical results that connect Deep Gaussian Processes with neural networks, the use of the above mean for the predictive distribution is grounded and justified for feed-forward neural networks of arbitrary depth. These results are further extended in <cit.>, where the authors make a theoretically-sound argument for the application of dropout to the case of recurrent neural networks. In essence, the masks need to be fixed for all the recurrent passes. This is the approach we follow when performing the experiments illustrated in the remainder of this article.

We note that it is straightforward to reinterpret our ordinal predictive distribution as one over the original, pre-quantised time series. Consider a categorical output ŷ_t over the partition C of intervals C_i, each of length |C_i|. Then we can straightforwardly construct a piece-wise uniform probability density function over the range of the time series I:

    p(𝐱_t+1|𝐱_t,...,𝐱_t-P+1)    = p_MOR(𝐱_t+1|ŷ_t)
       = ∏_i=1^Mŷ_t(i)/|C_i|^𝕀[𝐱_t+1∈ C_i]
    logp_MOR(𝐱_t+1|ŷ_t)   = ∑_i=1^Mπ_t+1^ilogŷ_t(i)/|C_i|
    π_t+1^i    = 𝕀[𝐱_t+1∈ C_i]

Quantities such as the negative sequence log-likelihood for a sequence of P_h samples drawn from the pre-quantised time series can now be calculated as:

    NLL_MOR   = -log∏_k=1^P_h∏_i=1^Mŷ_N+k-1(i)/|C_i|^π_N+k^i
       =-∑_k=1^P_h∑_i=1^Mπ_N+k^ilogŷ_N+k-1(i)/|C_i|





 §.§ REMARKS

Two crucial aspects to highlight about our framework are the following. Firstly, we achieve end-to-end time series forecasting by allowing our model to learn directly from data. Importantly, we do not make of use of handcrafted features, nor require expert knowledge to specify the core parts of the model (as often required in GP kernel design for instance).

Secondly, our model learns a non-parametric predictive posterior at each time step. This allows for rich, multi-modal behaviour to be encoded and fed back in an autoregressive fashion. We have observed empirically that providing our model with full knowledge of its own previous decisions enables it to achieve more reliable, longer-term forecasting than by one-hot encoding some statistic (such as the empirical mean or the mode) at every time step.









 §.§ BASELINE MODELS


In this Section, we review the baseline models we use to compare our framework with. These largely belong to three different families:


    
  * RNN regression. This baseline is a simple counterpart of our model that motivates the transition from direct to ordinal regression with RNNs. The same sequence-to-sequence architecture described in Section <ref>, but learnt by minimising the mean squared error of the model output with respect to the raw time series.
    
  * AR(p) modelling formulated as a state-space model. Autoregressive modelling of order p is a well-established statistical method, and widely used as a benchmark for time series forecasting tasks. We simply recapitulate that AR(p) models with time-varying coefficients can be recast as linear Gaussian state-space models. Such systems admit exact parameter inference and prediction in an efficient and principled manner by means of the Kalman filter. We refer the reader to <cit.> for a more detailed presentation of this model.
    
  * Gaussian Process autoregression. GPs, and in particular autoregressive GPs, have dominated the time series literature over the last decade. One potential limitation of this approach is that iterative one-step forecasting indeed requires GPs able to handle noisy inputs, which is still an active topic of research <cit.>. We thus propose two autoregressive GP baselines: one in which the one-step-ahead predictive posterior distributions are approximated via Monte-Carlo integration (described in Section <ref>), and an alternative approach in which each predictive posterior is approximated with a Bayesian Gaussian Mixture Model (described in Section <ref>). The latter approach requires more computational resources, but also enables the model to represent rich, multi-modal behaviour and adapt it over time.




  §.§.§ Autoregressive Gaussian Processes
 

Gaussian processes describe distributions over functions. They are defined as a collection of random variables, any finite number of which have a joint Gaussian distribution. A GP f is fully specified by a mean function m(𝐱) = 𝔼[f(𝐱)] and a covariance function k(𝐱, 𝐱) = 𝔼[(f(𝐱) - m(𝐱))(f(𝐱') - m(𝐱'))]. The covariance function (or kernel) describes correlations between data points and can be used to encapsulate prior belief about the problem. We refer the reader to <cit.> for a more detailed discussion about Gaussian Processes.

Consider a time series of length N. Long-term forecasting can be done in an autoregressive manner by learning a Gaussian Process f that maps a sequence of length P, 𝐗 = (𝐱_N-P-1, …, 𝐱_N), to a normal distribution over its next sample 𝒩(𝐱_N+1|μ_N+1, σ^2_N+1). However, sequentially inferring μ_t, σ^2_t for any t>N+1 requires adequate forward-propagation of the uncertainty computed in all previous forecasts. This effectively turns the task into one of GP regression with noisy inputs.





Uncertain, or noisy input, GPs are analytically intractable and therefore developing approximations remains an active area of research. In our work, we use a Monte Carlo-based approach <cit.>. S_GP sample trajectories are drawn up to time N+P_h for a predictive horizon P_h, i.e. 𝐗̂^[s] = (𝐱̂^[s]_N+1, …, 𝐱̂^[s]_N+P_h), and used to estimate the moments μ̃_N+k, σ̃^2_N+k of the corrected normal distribution p(𝐱_N+k|μ̃_N+k, σ̃^2_N+k).

Each sample trajectory is built iteratively ∀ k, 2≤ k≤ P_h in the following fashion: 

	
  * at time N+k-1, compute 𝒩(𝐱_N+k-1| μ_N+k-1, σ^2_N+k-1) from the sequence 𝐗̂^(N+k-2),
	
  * draw a sample 𝐱̂^[s]_N+k-1 from 𝒩(𝐱_N+k-1|  μ_N+k-1, σ^2_N+k-1), and build the next input seed sequence 
	
    𝐗̂^(N+k-1) = (𝐱_N+k-P-1, …, 𝐱_N+k-1)

	
  * repeat from step 1 until k>P_h.

After repeating the procedure above S_GP times, the corrected mean and variance ∀ k, 2≤ k ≤ P_h are thence given by:

    μ̃_N+k   = 𝔼[𝐱_N+k]
    	≈1/S_GP∑_s=1^S_GP𝐱̂^[s]_N+k
    σ̃^2_N+k   = 𝔼[ (𝐱_N+k)^2 ] - 𝔼[𝐱_N+k]^2 
       ≈1/S_GP∑_s=1^S_GP(𝐱̂^[s]_N+k - μ̃_N+k)^2




  §.§.§ Bayesian Gaussian Mixture Models
  
The baseline described above assumes each one-step-ahead predictive distribution is uni-modal Gaussian. We propose an alternative baseline in which we characterise the sample trajectories just described as a collection of Gaussian Mixture Models (GMM), which allows for one-step-ahead predictive distributions with a variable number of modes at each time step. In other words, the prediction 𝐱_t+1∼GMM(𝐱_t+1|π_t+1,μ_t+1,σ^2_t+1) with

    GMM(𝐱_t+1|π_t+1,μ_t+1,σ^2_t+1) =
    ∑_k=1^K π_k, t+1𝒩(𝐱_t+1|μ_k, t+1,σ^2_k, t+1)

where each π_k, t+1 is the mixture weight of the individual Gaussians and K is a hyperparameter for the number of mixtures in the model. Indeed this enables the description of richer behaviour than the one given by uni-modal Gaussians, but comes at the expense of the potentially computationally expensive process of learning up to P_h distinct GMMs.

There is extensive literature available on the inference of the parameters (π_t+1,μ_t+1,σ^2_t+1), and indeed on model selection for the appropriate number of mixtures K. State-of-the-art methods however aim to perform Bayesian inference over the learnable parameters, i.e. inferring a posterior distribution over each learnable parameter that best explains the observed evidence or data. This enables a principled prediction framework for other quantities of interest and further guards against overfitting and other well-documented pathologies that occur in point-estimate inference mechanisms, such as Maximum a Posteriori estimation <cit.>. Furthermore, such Bayesian methods also allow for natural shrinkage of the number of mixtures K, i.e. the weights of less informative mixtures naturally shrink to 0 and thus one need only provide an upper bound on K. We note that, although full Bayesian inference is intractable, efficient approximate Bayesian inference can be performed via Variational Bayes techniques, which is the approach we take in this article. We refer the reader to <cit.>, <cit.>, <cit.> and <cit.> for a more detailed discussion of Variational Bayes GMMs.













	
	

§ EXPERIMENTS



In this Section, we discuss our model's performance in two different long-term different tasks following the presentation of our baselines in Section <ref>.



 §.§ DATA




Our large-scale experiment was performed on 45 datasets compiled from different sources <cit.>. We used the following criteria to obtain this compilation:


    
  * Length. The compilation includes time series of different lengths, mostly ranging between 1,000 and 20,000 samples. In this work, we aim to compare our model in a setting in which poor performance should not be attributed to lack of data. We therefore only use time series with at least 10,000 observations, and further truncate those that contain over 30,000 so as to retain computational tractability during model fitting.
    
  * Structure. All time series were drawn from either the provided real-world references, or are a synthetic dynamical system.
    
  * Domain variety. Time series were queried and grouped by application domain, and at most two from each resulting group were subsequently drawn at random.


The selected datasets correspond with a wide range of dynamical systems, chaotic maps, and data obtained from varied application domains, including meteorology, astronomy and anatomy. A summary list and short descriptions of the selected datasets is given in Appendix A, and full descriptions can be found in the supplementary material of <cit.>. All the time series in the library were linearly and seasonally de-trended as detailed in the corresponding sources. Furthermore, all time series were standardised to zero mean, unit variance.




 §.§ TASK 1: LONG-TERM FORECASTING

In this experiment, we perform a large-scale comparison of our MOrdReD framework against the 4 baselines described in Section <ref>. Our long-term forecasting task consists in extrapolating P_h emissions (𝐱̂_N+1, …, 𝐱̂_N+P_h) forward from the last P observations (𝐱_N-P, ..., 𝐱_N) in the time series of length N, with P<<P_h. In our setting, we chose a lookback P=100, and a predictive horizon P_h=1000. Such iterative prediction is obtained for each of the 45 datasets described in Section <ref> without looking at real testing data or re-estimating model parameters during the process. 

We note that the aim of this forecasting task is to showcase the performance of our method as an “out-of-the-box" methodology, and thus we perform our comparison prior to any application domain crafting takes place (such as tailoring a task-optimal kernel, or appending handcrafted features and exogenous or correlated time series). We now give specific details about data preparation, model selection, model implementation and evaluation metrics.




  §.§.§ Data pre-processing
  
We used 45 different datasets drawn from both synthetic and real-world sources. Synthetic data includes an ample number of dynamical systems and chaotic maps, and real-world data was drawn from a number of sources and application domains, which encompass Meteorology, Astronomy, Physiology, Acoustics, and others. A short list is provided in Appendix A, and a more detailed description is given in <cit.>. All datasets were linearly and seasonally de-trended as detailed in Section 2.6 of the complementary material of <cit.>. Furthermore, all time series were standardised to zero mean, unit variance.  In the case of MOrdReD, all time series were quantised into a maximum of M=300 bins. For the rest of the models, which perform direct regression, regularising white noise with σ=10^-3 was added to the standardised time series. Datasets were further split into training, validation and testing time series. The training sets correspond to the first 70% of the observations in each time series, validation sets correspond to the subsequent 15%, and test sets to the last 15%. 



  §.§.§ Model learning, selection and implementation
  
We now provide details for the optimisation of each model's parameters and hyperparameters. The optimal hyperparameter configurations found for our neural network and AR(p) models are given in Appendix C.


    
  * MOrdReD. All our MOrdReD models were built in Python 2.7 and using Keras 2.1.2 with the Tensorflow backend <cit.>. Parameter optimisation was performed by minimising the categorical cross-entropy with Nesterov momentum Adam using α=0.002, β_1=0.9, β_2=0.999 and a schedule decay of 0.004, and parameters were initialised using the Glorot Uniform method. At training time, all models incorporate Early Stopping to guard against overfitting, while allowing for a maximum of 50 epochs using mini-batches of 256 training examples. Both the bidirectional encoder and the decoder share the same number of hidden units n_u∈{64, 128, 256, 320}, and the output fully-connected layers have M softmax output units corresponding to each time series' ordinal classes. All models further incorporate dropout and L2 regularisation with hyperparameters p_dout∈{0.25, 0.35, 0.5} and λ_L2∈{1e-6, 1e-7, 1e-8}, which were jointly optimised with the number of hidden units n_u by hypergrid search on the validation set. MC dropout predictive posteriors were estimated using N_s=100 samples.
    
  * Direct regression sequence-to-sequence neural network. These were learnt in a similar fashion as the above, but by minimising the mean squared error loss function of the single output unit (with no activation function) with respect to the ground truth. This baseline serves to motivate the transition into ordinal regression.
    
  * Autoregressive GP. All our GP models were implemented using GPy 1.9.2 <cit.> and using a Matérn 5/2 kernel with additive white noise and Automatic Relevance Determination enabled. Samples drawn from the Matérn 5/2 covariance function are twice-differentiable, which enables the modelling of a wide range of physical systems such as those in our datasets. Time-indexed Gaussian predictive distributions were then estimated following the Monte Carlo-based procedure described in Section <ref>, with N_s=100 samples trajectories.
    
  * ARGP with predictive Gaussian Mixture Model. This baseline follows the implementation described above, but the time-indexed predictive distributions are estimated via a Gaussian Mixture Model as described in Section <ref>. This allows for rich, multi-modal behaviour to be described succinctly at each time step. We performed approximate Bayesian inference to learn these models, with up to K=5 mixtures at each time step, using the open source Variational Bayes GMM implementation available in the Scikit-learn toolbox <cit.>, v. 0.19.1.
    
  * AR(p). The Python library StatsModels 0.9 <cit.> offers an open source implementation of the state-space formulation of AR(p) modelling. We used this and optimised the lookback hyperparameter p∈{16, 32, 64}.


Our experiments are further accompanied by a Python library that provides an implementation of our methodology, as well as an interface for the ample number of libraries that implement our baselines. This code and the scripts used to perform these experiments can be found at <https://www.github.com/bperezorozco/ordinal_tsf>.



  §.§.§ Evaluation metrics
  
Our evaluation is performed across three axes: predictive accuracy, uncertainty quantification and forecast reliability. This exhaustive comparison enables us to speak not only about the performance of the forecasts predicted by our model, but also about the credibility and reliability of the uncertainty bounds they infer. In terms of forecast accuracy, we measure the Symmetric Mean Absolute Percentage Error (SMAPE) and the Root Mean Squared Error (RMSE) of both the mean and the median of the predictive distributions of each method:

    SMAPE   = 2/P_h∑_t'=N^N+P_h|𝐱_t' - 𝐱̂_t'|/|𝐱_t'| + |𝐱̂_t'|
    RMSE   = √(∑_t'=N^N+P_h‖𝐱_t' - 𝐱̂_t'‖^2/P_h)


Crucially, one goal in this experiment is not just to evaluate a model's predictive accuracy, but also to measure its honesty expressed through uncertainty estimations. We prefer models that are able to either produce an accurate confident forecast, or otherwise explicitly confess their ignorance - for instance, through a mean-reversal process. 

We quantify this via the negative sequence log-likelihood (NLL), which we rewrite using the chain rule for joint distributions:

    -logp(𝐗^(N+P_h)) =    -logp(𝐱_N+1)
       -∑_k=2^P_hlogp(𝐱_N+k|⋂_k'=1^k-1𝐱_N+k')

	

For MOrdReD, each term in the summation acquires the piece-wise uniform shape described in Section <ref>; for ARGPs with a GMM predictive distribution, it acquires the shape described in Section <ref>; for all the other baselines, it is given by a standard Gaussian density. 

In addition to the NLL, we also provide the integrated or cumulative NLL:


    CNLL = ∑_i=1^P_h-logp(𝐗^(N+i))


This metric summarises information about how the NLL changes with time, e.g. greater penalties are incurred by models that make inaccurate short-term predictions, whereas erroneous long-term forecasts have a lesser contribution.

Lastly, we benchmark our methodology in terms of the calibration of its uncertainty bound predictions. That is, we measure its ability to produce output densities that can be interpreted as real-world probabilities. Graphics tools such as quantile-quantile (QQ) plots and reliability diagrams are widespread in the time series literature to measure uncertainty calibration. However, due to the large-scale nature of this task, we instead propose a related summary metric of these which we now introduce. 

Consider the time series 𝐱̂_α whose t-th entry is given by the α-quantile of the predictive distribution at time t, i.e.:

    𝐱̂_α,t = q_t    s.t._-∞^q_tp(𝐱_t|⋂_t'=1^t-1𝐱_t')d𝐱_t = α


For a given calibrated model and α, 0 < α < 1, we expect that exactly α% of the observations in the ground truth 𝐱 will lie below 𝐱̂_α. Let r_α be this proportion. Then for better calibrated models, r_α→α with:


    r_α = 1/P_h∑_t=N+1^N+P_h𝕀[𝐱_t < 𝐱̂_α, t]


We thus define the calibration metric QQDist as the (Euclidean) distance between r_α and α:

    QQDist = ∫_0^1(r_α - α)^2dα


We provide two instances of this metric. One in which the r_α for each method are computed for the full predictive horizon P_h=1000, and one up to a truncated horizon P'_h=250. We provide both instances so as to distinguish calibration performance in both the medium and long term.

Performance is measured over the following 6 statistics: SMAPE and MASE for the predictive median's accuracy; NLL and CNLL to evaluate the predictive distributions' accuracy; and QQDist and QQDist-250 for uncertainty calibration.



  §.§.§ Results

We summarise our results in Tables <ref>, <ref> and <ref>. We note that our model achieves state-of-the-art performance in an wholly end-to-end fashion, with no human intervention to specify the core architecture of the model. Unlike GP Regression, where knowledge of kernels is required and often inaccessible to unfamiliar users, our framework is end-to-end and  requires almost no user intervention. Additionally, other advantages of neural networks such as scalability are now readily accessible. Our framework can make the most out of big datasets and produce reliable forecasts in the long term that outperform other state-of-the-art techniques.

From Table <ref> we note that our model is the one that performs best in the largest number of datasets across all the proposed metrics. Interestingly, we see that the worst performing baseline is the one given by direct regression neural networks. This gives further empirical motivation to transition from direct into ordinal regression for the family of neural network models. We also note that classical statistical models perform very well in the calibration metrics, whereas GPs are the worst performing in this axis. We argue that this could be the case due to the exact sequential inference method of the predictive variance proposed in the state-space formulation of AR(p) models. Finally, GPs are the closest competitors in terms of forecasting accuracy.

We note that, for some time series, our model is the only one capable of learning any structure for long-term forecasting. One such example is the electrocardiogram time series, depicted in Figure <ref>. This complex-shaped signal is consistently predicted in a timely and accurate fashion by MOrdReD in the long-term, whereas no baseline is capable of reproducing the characteristic shape of the QRS complex beyond its first occurence.

We now focus on Table <ref>, which provides the mean performance rank for each baseline and metric. This provides further information about how our method performs when it is not the best-performing one, and we observe that MOrdReD still achieves the best mean rank for all metrics. In order to produce this table, two considerations were taken into account: on one hand, both GP baselines could propose very similar predictive distributions for some datasets (e.g. those in which multi-modality is not required and therefore both propose a uni-modal Gaussian), and this induced twice the penalty for those methods that underperformed against GPs. We solve this by merging the results of both GPs and retaining the better performing one. This avoids allocating a double penalty on models that underperform with respect to GPs if both GPs propose equivalent predictive distributions (e.g. if the optimal number of Gaussian mixtures is 1).

The information in Table <ref> crucially summarises that whenever our model is not the best performing, it is still likely to be the second best-performing, etc. An example of this scenario is given in Figure <ref>, which provides our forecasts for the tide height dataset. The GP baseline is the best performing, as it achieves better timing accuracy in the longer term in comparison with MOrdReD. However, the remaining baselines are not able to capture the basic structure of the signal at all, and therefore our model comes out as a clear runner-up.

Nevertheless, we also observed that MOrdReD is not always so close to the best-performing model. For instance, we noted that GPs are often the best-performing in chaotic maps and dynamical systems. An example is given in Figure <ref>, where we provide the forecasts for a Lorenz map. Whereas both MOrdReD and GPs are able to model multi-modal behaviour in the long-term, MOrdReD syncs out with respect to the ground truth considerably faster than GPs. Indeed, GPs only make use of their multi-modal predictive capability relatively late in comparison with MOrdReD.

In a similar fashion to Table <ref>, in Table <ref> we introduce the reciprocal metric: mean worst-performance rank. This metric speaks about the frequency in which a method is the worst performing for a given task. Rank 3 for each task in this case is given to the best-performing model, and thus larger overall metric values indicate better performance. These results indicate that there are relatively few cases in which all other models outperform it for any given metric in the assessment.


 

 







 §.§ TASK 2: QUANTIFYING EVENT TIMING UNCERTAINTY

As motivated earlier, a crucial component of long-term prediction is timely forecasting. Decision-makers often prioritise knowing the time of occurrence of an event over the magnitude of the event itself. For instance, consider the tide height dataset of Section <ref>. Experts are often interested in knowing when the tide levels will reach their maximum, and not just the actual height they will reach.[Despite this being a periodic phenomenon, sensor quantisation leads to representation errors that make the time series quasi-periodic in practice.]

In a similar vein, cardiologists are interested, for instance, in forecasting the timing of events such as the QRS complex. Predicting these timings further enables them to compute other metrics of interest, such as the RR interval – which has been linked, for example, to Parkinson's disease. Accurately quantifying these metrics is therefore of utmost importance, and fully relies on models that can speak honestly about their forecasts and their uncertainties.

 The task we now focus on is to predict when a certain event will happen, alongside reliable uncertainty bounds. In this setting, we compare our MOrdReD framework with the same baselines described in Section <ref>. However, in contrast to Task 1 in Section <ref>, we focus our attention on only 8 datasets drawn from meteorology and physiology. We motivate our decision as follows: on one hand, evaluating the occurrence of local maxima is not relevant for some data (e.g. acoustic signals); on the other, we wish to focus our attention on real-world data, subject to realistic noise. 



  §.§.§ Constructing densities for critical event occurrence


In order to adequately quantify event timings, we construct a non-parametric probability distribution p(t) that describes the probability of an event happening at time t. In our setting, we take, as example, this event to be the time series hitting a local maximum. To compute such timing forecasts, we employ Kernel Density Estimation (KDE) to construct a probability density function p(t) that describes the probability that a certain event (hitting a local maximum in this example) will happen at time t. KDE is a smoothing technique that builds a non-parametric distribution from a sample set X. The estimator is given by


    p(t) = 1/n∑_i=1^|X|K_h(t-t_i),

where t_i∈ X, K_h is a non-negative function that integrates to one, and h is the bandwidth hyperparameter. The reader is referred to <cit.> for a more detailed discussion. In our work, we used the freely available implementation of Scikit-learn <cit.>. Now consider the true event timings X^(true) = (t_1^(true), …, t_L'^(true)). The constructed densities are then evaluated in terms of the negative log-likelihood for effectively predicting such correct timings:

    NLL_X^true = -log∏_i=1^L'p(t_i^(true)) = -∑_i=1^L'logp(t_i^(true))


We now briefly describe how to construct the sets X, X^(true). In the case of X, sample trajectories are drawn from the predictive densities of each framework, and the timing of the desired event in each sample is recorded. An estimator p(t) is subsequently built for each method. In the case of X^(true), the true event timings are recorded directly from the ground truth. Maxima detection was achieved in our experimental setting with the open source PeakUtils Python library. Recording the timing of events can be inhibited by the presence of noise, even in the ground truth data. To construct the set X^(true), we decomposed each dataset's ground truth using Empirical Mode Decomposition <cit.> and then recorded the events from the (noiseless) intrinsic mode function that captured long-term quasi-periodic nature. To this end, we used an additional Python open source library, PyEMD.



  §.§.§ Results





In Table <ref>, the negative log-likelihood values for drawing the true timings X^(true) from the constructed densities p_ARGP(t), p_MOR(t) are shown for predicting the occurrence of a local maximum in 8 different datasets. We observed that our framework does best in more dataset examples than the rest of the baselines, closely followed by the GP baseline, which dominates in the remaining datasets. 

We further provide a visual example of these densities in Figures <ref> and <ref>. In the former, we observe that MOrdReD accurately predicts the timing of all maxima in the long-term, clearly outperforming all the other baselines in this task. On the contrary, in the case of the tide height dataset, we observe that the GP baseline is the one that accurately forecasts the timing of all maxima, since MOrdReD syncs out with respect to the ground truth towards the end of the forecast.








	
	

§ CONCLUDING REMARKS



In this article we introduce a novel, fully end-to-end ordinal time series forecasting method based on recurrent neural networks. By recasting time series forecasting as an ordinal regression task, we demonstrate how state-of-the-art Machine Learning methods can be employed to perform reliable long-term time series prediction in a scalable and general-purpose fashion. We show how recent developments in the Bayesian neural networks literature can be used to quantify predictive uncertainty in time series forecasting, a major task that the literature of time series prediction with ANNs has paid relatively little attention to. Furthermore, our ordinal framework enables the description of adaptable multi-modal, non-Gaussian behaviour through time. 

Crucially, we highlight that these results are achieved in a fully end-to-end fashion. Little to no human intervention is needed to use MOrdReD for model design decisions. Direct learning from data is thence enabled without compromising scalability, assuming restrictive shapes of predictive distributions, or requiring expert craftsmanship for model design, as is often the case in designing a kernel for Gaussian Processes. 

In order to assess the performance of our framework, we provide a large-scale benchmark test over 45 different datasets drawn from an ample range of application domains and synthetic maps. Our 45 datasets were drawn from a variety of sources and include time series with quasi-periodic behaviour and complex shapes that in some cases are only modelled accurately by MOrdReD, as we exemplify with the electrocardiogram case. We compare our method with state-of-the-art baselines in the time series literature from both the Statistics and the Machine Learning perspectives, and evaluate a number of metrics to do with predictive accuracy, uncertainty quantification and forecast reliability.

We find that our model is empirically capable of outperforming other state-of-the-art competitors in terms of long-term forecasting uncertainty estimation, while also inheriting all the advantages of neural network models, such as scalability. As a result, our framework is consistently top-ranked in all metrics. We additionally provide evidence that it would also be the second-best performing model whenever another baseline achieved the best rank, and in a similar vein, we show that it is unlikely that it will be outperformed by all other baselines simultaneously in any given task, for any metric considered in the assessment.

Finally, we demonstrate how our model can be used to construct the predictive distribution of the occurrence of critical events of interest in applied settings. We show that our framework yields long-term, reliable confidence intervals even in safety-critical environments such as cardiology and meteorology.

This study, though comprehensive, still assumes the existence of enough data to fit our model, and only considers the case of bounded univariate time series autoregression. Future directions include extending this framework to the univariate case with exogenous observations, and more generally, to the multivariate case; as well as extensions that can achieve competitive performance even in presence of little data.
	
	

§ ACKNOWLEDGEMENTS


The authors would like to thank Xiaowen Dong and other members of the Machine Learning Research Group at the University of Oxford for their feedback and edits on the paper clarity and structure.



§ CONFLICT OF INTEREST


The authors of this manuscript declare that there is no conflict of interest between this manuscript and other published works.
	
	

Bernardo Pérez Orozco is funded by Consejo Nacional de Ciencia y Tecnología (CONACYT) in Mexico CVU ID #598304. Gabriele Abbati was supported by University of Oxford,
Google DeepMind. We also gratefully acknowledge the support provided by the Royal Academy of Engineering and the Oxford-Man Institute of Quantitative Finance.
	
	
	
	
	

§ APPENDIX A

In this appendix we provide a list of the 45 datasets used throughout the article. Those marked with an (*) were further used for the event occurrence forecasting task. Unless otherwise stated, these were all drawn from the compilation made by <cit.>.


    
  * {MACKEY} Mackey-Glass chaotic attractor, synthesised by the authors according to the system:
	
    x_t+1 = (1 - b)x_t + ax_t-τ/1 + x_t-τ^n

	
	with a=0.2, b=0.1, τ=17, n=10 and discarding the first 1,000 burnout samples.
	
	
  * (*) {ECG} a recording of the electrical activity of the human heart <cit.>, which highlights forecasting of a quasi-periodic, complex-shaped signal where accurate timing is essential;
	
  * (*) {AIRFLOW} a recording of the breathing activity of a human<cit.>, which highlights forecasting of a quasi-periodic signal where each quasi-cycle has a relatively long time period;
	
  * (*) {TIDE} data taken from the Sotonmet <cit.> environmental sensors that record tide heights in an off-shore weather station. 
    
  * {AS_s3.2_} birdsong excerpt from the Macaulay Animal Sounds Library,
    
  * (*) {CM_air.s} air temperature time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * {CM_air.2} a second air temperature time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * {CM_prate} precipitation rate time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * (*){CM_slp19} a sea level pressure time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * {CM_SLP.2} a second sea level pressure time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * (*) {CM_rhum1} relative humidity time series taken by different sensors in the Eurasia region between 1948 and 2007,
    
  * {CM_lwtla} Lamb/Jenkins weather type series measured from 1861 to 1997,
    
  * {EMexptqp} Quasi-periodic output of Eric Weeks' Annulus experiment,
    
  * {EM_henon} Henon map 
    
    x' =    a + by - x^2 
    
            y' =    x

    with a=1.4, b=0.3
    
    
  * {EMlorenz} Lorenz map 
    
    x' =   σ (y-x)
    
            y' =    rx -y - xz
     
            z' =    xy - bz

    with σ=10, r=28, b=8/3,
    
    
  * {EM_rossl} Rössler attractor
    
    x' =    -z-y 
     y' =    x + ay 
     z' =    b + z(x-c)

    with a=0.15, b=0.20, c=10.0,
    
    
  * {FI_yahoo} Log returns of GPSC stock chart,
    
  * {FL_ACT_L} z-channel of the ACT attractor:
    
    x' =   α(x-y)
     y' =    -4α y + xz + μ x^3 
     z' =    -δα z + xy + β z^2

    computed with α=1.8, β=-0.07, δ=1.5, μ=0.02,
    
    
  * {FL_chen_} x-channel of Chen's system 
    
    x' =    a(y-x)
     y' =    (c-a)x -xz + cy 
     z' =    xy - bz

    computed with a=35, b=3, c=28,
    
    
  * {FL_dblsc} y-channel of the double-scroll system 
    
    x' =    y 
     y' =    z 
     z' =    -a[z + y + x - sgn(x)]

    computed with a=0.8 and initial conditions x_0=0.01, y_0=0.01, z_0=0,
    
    
  * {FL_hadle} x-channel of the Hadley circulation system 
    
    x' =    -y^2 - z^2 -ax + aF
     y' =   xy -bxz -y + G 
     z' =    bxy + xz - z

    computed with a=0.25, b=4, F=8, G=1 and initial conditions x_0=0, y_0=0, z_0=1.3,
    
    
  * {FL_labyr} y-channel of the Labyrinth Chaos system 
    
    x' =   sin(y) 
     y' =    -sin(z) 
     z' =   sin(x)

    computed with initial conditions x_0=0.1, y_0=0, z_0=0,
    
    
  * {FL_moore} z-channel of the Moore-Spiegel oscillator 
    
    x' =    y 
     y' =    z 
     z' =    -z - (T-R + Rx^2)y - Tx

    computed with T=6, R=20,
    
    
  * {FL_noseh} z-channel of the Nosé-Hoover oscillator 
    
    x' =    y 
     y' =    -x + yz 
     z' =    a - y^2

    computed with a=1 and initial conditions x_0=0, y_0=5, z_0=0,
    
    
  * {FL_ruckl} z-channel of the Rucklidge attractor 
    
    x' =    -κ x + λ y - yz 
     y' =    x 
     z' =    -z + y^2

    computed with κ=2, λ=6.7 and initial conditions x_0=1, y_0=0, z_0=4.5,
    
    
  * {FL_simpq} y-channel of the simplest quadratic flow 
    
    x' =    y 
     y' =    z 
     z' =    -az + y^2 - x

    computed with a=2.028 and initial conditions x_0=0.9, y_0=0, z_0=0.5,
    
    
  * {FL_thoma} y-channel of Thomas cyclically symmetric attractor 
    
    x' =    -bx + sin(y) 
     y' =    -by + sin(z) 
     z'  -bz + sin(x)

    computed with b=0.18 and initial conditions x_0=0.1, y_0=0, z_0=0,
    
    
  * {FL_windm} y-channel of the WINDMI attractor 
    
    x' =    y 
     y' =    z 
     z' =    -az - y + b - exp(x)

    for a=0.7, b=2.5 and initial conditions  x_0=0, y_0=0.8, z_0=0,
    
    
  * {K_standa} Log returns of the Standard & Poor index,
    
  * {MC_inttr} Internet traffic data from an ISP, provided by the Time Series Data Library,
    
  * {MP_Lozi_} Lozi map x_t+1 = 1 - a|x_t| + bx_n-1 with a=1.7, b=0.5 and initial conditions x_1=-0.1, x_0=0.15, 
    
  * {MP_freit} Freitas' stochastic sine map x_t+1=μsin(x_t) + Y_tη_t, with parameters μ=2.4, b=3, q=0.2, η_t∼Uniform(-b, b), Y_t∼Bernoulli(q),
    
  * {MP_logis} Logistic map x_t+1 = Ax_t(1-x_t) with A=3.2 and initial condition x_0=0.91,
    
  * {MUS_Si_l}  a music excerpt from B. Fulcher's personal collection, entitled Si loin de vous,
    
  * {MUS.3_78} a second music excerpt from B. Fulcher's personal collection, entitled 3,
    
  * {SFX_mach} Sound Jay Mach Electric Drill sound effect,
    
  * {SF_Acont}  Santa Fe laser generated data excerpt,
    
  * {SF_B1_1} Santa Fe heart rate data excerpt,
    
  * {SF_D1}  Santa Fe synthetically generated snippet,
    
  * {SL_perci} StatLib demeaned water level measurements,
    
  * {SPIDR_hp} Hemispheric Power Index excerpt provided by the Space Physics Interactive Data Resource,
    
  * {SY_AR2_T} Timmer nonstationary autoregressive process with τ=20, T_mean=20, T_mod=20, σ=1, ℳ_T=5, η=1000,
    
  * {SY_NLAR2} Faes nonlinear autoregressive process with a_1=3.6, a_2=0.8,
    
  * {TSAR_eqe} Earthquake and explosion seismic series provided in Stoffer's Time Series Analysis and its applications, with R examples.
    
  * {TXT_slc_} Project Gutenberg excerpt of Dickens' Oliver Twist.


	

§ APPENDIX B

In this appendix we provide the full results for the performance of each model on each dataset for all metrics. Full details of the benchmarking task and the definition of our metrics are given in Section <ref>.



 

 

 

 

 

 

 


	

§ APPENDIX C

In this appendix we provide the optimised hyperparameters for our MOrdReD and AR(p) models. 








	
	
