Title: N-BEATS: NEURAL BASIS EXPANSION ANALYSIS FOR INTERPRETABLE TIME SERIES FORECASTING

Abstract: We focus on solving the univariate times series point forecasting problem using deep learning. We propose a deep neural architecture based on backward and forward residual links and a very deep stack of fully-connected layers. The architecture has a number of desirable properties, being interpretable, applicable without modification to a wide array of target domains, and fast to train. We test the proposed architecture on several well-known datasets, including M3, M4 and TOURISM competition datasets containing time series from diverse domains. We demonstrate state-of-the-art performance for two configurations of N-BEATS for all the datasets, improving forecast accuracy by 11% over a statistical benchmark and by 3% over last year's winner of the M4 competition, a domain-adjusted hand-crafted hybrid between neural network and statistical time series models. The first configuration of our model does not employ any time-series-specific components and its performance on heterogeneous datasets strongly suggests that, contrarily to received wisdom, deep learning primitives such as residual blocks are by themselves sufficient to solve a wide range of forecasting problems. Finally, we demonstrate how the proposed architecture can be augmented to provide outputs that are interpretable without considerable loss in accuracy.

1. INTRODUCTION:

Time series (TS) forecasting is an important business problem and a fruitful application area for machine learning (ML). It underlies most aspects of modern business, including such critical areas as inventory control and customer management, as well as business planning going from production and distribution to finance and marketing. As such, it has a considerable financial impact, often ranging in the millions of dollars for every point of forecasting accuracy gained (Jain, 2017; Kahn, 2003) . And yet, unlike areas such as computer vision or natural language processing where deep learning (DL) techniques are now well entrenched, there still exists evidence that ML and DL struggle to outperform classical statistical TS forecasting approaches (Makridakis et al., 2018a; b) . For instance, the rankings of the six "pure" ML methods submitted to M4 competition were 23, 37, 38, 48, 54, and 57 out of a total of 60 entries, and most of the best-ranking methods were ensembles of classical statistical techniques (Makridakis et al., 2018b) .

On the other hand, the M4 competition winner (Smyl, 2020) , was based on a hybrid between neural residual/attention dilated LSTM stack with a classical Holt-Winters statistical model (Holt, 1957; 2004; Winters, 1960) with learnable parameters. Since Smyl's approach heavily depends on this Holt-Winters component, Makridakis et al. (2018b) further argue that "hybrid approaches and combinations of method are the way forward for improving the forecasting accuracy and making forecasting more valuable". In this work we aspire to challenge this conclusion by exploring the potential of pure DL architectures in the context of the TS forecasting. Moreover, in the context of interpretable DL architecture design, we are interested in answering the following question: can we inject a suitable inductive bias in the model to make its internal operations more interpretable, in the sense of extracting some explainable driving factors combining to produce a given forecast?

1.1. SUMMARY OF CONTRIBUTIONS:

Deep Neural Architecture: To the best of our knowledge, this is the first work to empirically demonstrate that pure DL using no time-series specific components outperforms well-established statistical approaches on M3, M4 and TOURISM datasets (on M4, by 11% over statistical benchmark, by 7% over the best statistical entry, and by 3% over the M4 competition winner). In our view, this provides a long-missing proof of concept for the use of pure ML in TS forecasting and strengthens motivation to continue advancing the research in this area.

Interpretable DL for Time Series: In addition to accuracy benefits, we also show that it is feasible to design an architecture with interpretable outputs that can be used by practitioners in very much the same way as traditional decomposition techniques such as the "seasonality-trend-level" approach (Cleveland et al., 1990) .

2. PROBLEM STATEMENT:

We consider the univariate point forecasting problem in discrete time. Given a length-H forecast horizon a length-T observed series history [y 1 , . . . , y T ] ∈ R T , the task is to predict the vector of future values y ∈ R H = [y T +1 , y T +2 , . . . , y T +H ]. For simplicity, we will later consider a lookback window of length t ≤ T ending with the last observed value y T to serve as model input, and denoted x ∈ R t = [y T -t+1 , . . . , y T ]. We denote y the forecast of y. The following metrics are commonly used to evaluate forecasting performance (Hyndman & Koehler, 2006; Makridakis & Hibon, 2000; Makridakis et al., 2018b; Athanasopoulos et al., 2011) : .

Here m is the periodicity of the data (e.g., 12 for monthly series). MAPE (Mean Absolute Percentage Error), sMAPE (symmetric MAPE) and MASE (Mean Absolute Scaled Error) are standard scale-free metrics in the practice of forecasting (Hyndman & Koehler, 2006; Makridakis & Hibon, 2000) : whereas sMAPE scales the error by the average between the forecast and ground truth, the MASE scales by the average error of the naïve predictor that simply copies the observation measured m periods in the past, thereby accounting for seasonality. OWA (overall weighted average) is a M4specific metric used to rank competition entries (M4 Team, 2018b) , where sMAPE and MASE metrics are normalized such that a seasonally-adjusted naïve forecast obtains OWA = 1.0.

3. N-BEATS:

Our architecture design methodology relies on a few key principles. First, the base architecture should be simple and generic, yet expressive (deep). Second, the architecture should not rely on timeseries-specific feature engineering or input scaling. These prerequisites let us explore the potential of pure DL architecture in TS forecasting. Finally, as a prerequisite to explore interpretability, the architecture should be extendable towards making its outputs human interpretable. We now discuss how those principles converge to the proposed architecture.

3.1. BASIC BLOCK:

The proposed basic building block has a fork architecture and is depicted in Fig. 1 (left). We focus on describing the operation of -th block in this section in detail (note that the block index is dropped in Fig. 1 for brevity). The -th block accepts its respective input x and outputs two vectors, x and y . For the very first block in the model, its respective x is the overall model input -a history lookback window of certain length ending with the last measured observation. We set the length of Figure 1 : Proposed architecture. The basic building block is a multi-layer FC network with RELU nonlinearities. It predicts basis expansion coefficients both forward, θ f , (forecast) and backward, θ b , (backcast). Blocks are organized into stacks using doubly residual stacking principle. A stack may have layers with shared g b and g f . Forecasts are aggregated in hierarchical fashion. This enables building a very deep neural network with interpretable outputs. input window to a multiple of the forecast horizon H, and typical lengths of x in our setup range from 2H to 7H. For the rest of the blocks, their inputs x are residual outputs of the previous blocks. Each block has two outputs: y , the block's forward forecast of length H; and x , the block's best estimate of x , also known as the 'backcast', given the constraints on the functional space that the block can use to approximate signals.

Internally, the basic building block consists of two parts. The first part is a fully connected network that produces the forward θ f and the backward θ b predictors of expansion coefficients (again, note that the block index is dropped for θ b , θ f , g b , g f in Fig. 1 for brevity). The second part consists of the backward g b and the forward g f basis layers that accept the respective forward θ f and backward θ b expansion coefficients, project them internally on the set of basis functions and produce the backcast x and the forecast outputs y defined in the previous paragraph.

The operation of the first part of the -th block is described by the following equations:

h ,1 = FC ,1 (x ), h ,2 = FC ,2 (h ,1 ), h ,3 = FC ,3 (h ,2 ), h ,4 = FC ,4 (h ,3 ). θ b = LINEAR b (h ,4 ), θ f = LINEAR f (h ,4 ).

(1)

Here LINEAR layer is simply a linear projection layer, i.e. θ f = W f h ,4 . The FC layer is a standard fully connected layer with RELU non-linearity (Nair & Hinton, 2010) , such that for FC ,1 we have, for example: h ,1 = RELU(W ,1 x + b ,1 ). One task of this part of the architecture is to predict the forward expansion coefficients θ f with the ultimate goal of optimizing the accuracy of the partial forecast y by properly mixing the basis vectors supplied by g f . Additionally, this sub-network predicts backward expansion coefficients θ b used by g b to produce an estimate of x with the ultimate goal of helping the downstream blocks by removing components of their input that are not helpful for forecasting.

The second part of the network maps expansion coefficients θ f and θ b to outputs via basis layers, y = g f (θ f ) and x = g b (θ b ). Its operation is described by the following equations:

y = dim(θ f ) ∑ i=1 θ f ,i v f i , x = dim(θ b ) ∑ i=1 θ b ,i v b i .

Here v f i and v b i are forecast and backcast basis vectors, θ f ,i is the i-th element of θ f . The function of g b and g f is to provide sufficiently rich sets

{v f i } dim(θ f ) i=1 and {v b i } dim(θ b ) i=1

such that their respective outputs can be represented adequately via varying expansion coefficients θ f and θ b . As shown below, g b and g f can either be chosen to be learnable or can be set to specific functional forms to reflect certain problem-specific inductive biases in order to appropriately constrain the structure of outputs. Concrete examples of g b and g f are discussed in Section 3.3.

3.2. DOUBLY RESIDUAL STACKING:

The classical residual network architecture adds the input of the stack of layers to its output before passing the result to the next stack (He et al., 2016) . The DenseNet architecture proposed by Huang et al. (2017) extends this principle by introducing extra connections from the output of each stack to the input of every other stack that follows it. These approaches provide clear advantages in improving the trainability of deep architectures. Their disadvantage in the context of this work is that they result in network structures that are difficult to interpret. We propose a novel hierarchical doubly residual topology depicted in Fig. 1 (middle and right). The proposed architecture has two residual branches, one running over backcast prediction of each layer and the other one is running over the forecast branch of each layer. Its operation is described by the following equations:

x = x -1 -x -1 , y = ∑ y .

As previously mentioned, in the special case of the very first block, its input is the model level input x, x 1 ≡ x. For all other blocks, the backcast residual branch x can be thought of as running a sequential analysis of the input signal. Previous block removes the portion of the signal x -1 that it can approximate well, making the forecast job of the downstream blocks easier. This structure also facilitates more fluid gradient backpropagation. More importantly, each block outputs a partial forecast y that is first aggregated at the stack level and then at the overall network level, providing a hierarchical decomposition. The final forecast y is the sum of all partial forecasts. In a generic model context, when stacks are allowed to have arbitrary g b and g f for each layer, this makes the network more transparent to gradient flows. In a special situation of deliberate structure enforced in g b and g f shared over a stack, explained next, this has the critical importance of enabling interpretability via the aggregation of meaningful partial forecasts.

3.3. INTERPRETABILITY:

We propose two configurations of the architecture, based on the selection of g b and g f . One of them is generic DL, the other one is augmented with certain inductive biases to be interpretable.

The generic architecture does not rely on TS-specific knowledge. We set g b and g f to be a linear projection of the previous layer output. In this case the outputs of block are described as:

y = V f θ f + b f , x = V b θ b + b b .

The interpretation of this model is that the FC layers in the basic building block depicted in Fig. 1 learn the predictive decomposition of the partial forecast y in the basis V f learned by the network. Matrix V f has dimensionality H × dim(θ f ). Therefore, the first dimension of V f has the interpretation of discrete time index in the forecast domain. The second dimension of the matrix has the interpretation of the indices of the basis functions, with θ f being the expansion coefficients for this basis. Thus the columns of V f can be thought of as waveforms in the time domain. Because no additional constraints are imposed on the form of V f , the waveforms learned by the deep model do not have inherent structure (and none is apparent in our experiments). This leads to y not being interpretable.

The interpretable architecture can be constructed by reusing the overall architectural approach in Fig. 1 and by adding structure to basis layers at stack level. Forecasting practitioners often use the decomposition of time series into trend and seasonality, such as those performed by the STL (Cleveland et al., 1990 ) and X13-ARIMA (U.S. Census Bureau, 2013) . We propose to design the trend and seasonality decomposition into the model to make the stack outputs more easily interpretable. Note that for the generic model the notion of stack was not necessary and the stack level indexing was omitted for clarity. Now we will consider both stack level and block level indexing. For example, y s, will denote the partial forecast of block within stack s.

Trend model. A typical characteristic of trend is that most of the time it is a monotonic function, or at least a slowly varying function. In order to mimic this behaviour we propose to constrain g b s, and g f s, to be a polynomial of small degree p, a function slowly varying across forecast window:

y s, = p ∑ i=0 θ f s, ,i t i .

(2)

Here time vector t = [0, 1, 2, . . . , H -2, H -1] T /H is defined on a discrete grid running from 0 to (H -1)/H, forecasting H steps ahead. Alternatively, the trend forecast in matrix form will then be:

y tr s, = Tθ f s, ,

where θ f s, are polynomial coefficients predicted by a FC network of layer of stack s described by equations (1); and T = [1, t, . . . , t p ] is the matrix of powers of t. If p is low, e.g. 2 or 3, it forces y tr s, to mimic trend.

Seasonality model. Typical characteristic of seasonality is that it is a regular, cyclical, recurring fluctuation. Therefore, to model seasonality, we propose to constrain g b s, and g f s, to belong to the class of periodic functions, i.e. y t = y t-∆ , where ∆ is a seasonality period. A natural choice for the basis to model periodic function is the Fourier series:

EQUATION

The seasonality forecast will then have the matrix form as follows:

y seas s, = Sθ f s, ,

where θ f s, are Fourier coefficients predicted by a FC network of layer of stack s described by equations (1); and S = [1, cos(2πt), . . . cos(2π H/2 -1 t)), sin(2πt), . . . , sin(2π H/2 -1 t))] is the matrix of sinusoidal waveforms. The forecast y seas s, is then a periodic function mimicking typical seasonal patterns.

The overall interpretable architecture consists of two stacks: the trend stack is followed by the seasonality stack. The doubly residual stacking combined with the forecast/backcast principle result in (i) the trend component being removed from the input window x before it is fed into the seasonality stack and (ii) the partial forecasts of trend and seasonality are available as separate interpretable outputs. Structurally, each of the stacks consists of several blocks connected with residual connections as depicted in Fig. 1 and each of them shares its respective, non-learnable g b s, and g f s, . The number of blocks is 3 for both trend and seasonality. We found that on top of sharing g b s, and g f s, , sharing all the weights across blocks in a stack resulted in better validation performance.

3.4. ENSEMBLING:

Ensembling is used by all the top entries in the M4-competition. We rely on ensembling as well to be comparable. We found that ensembling is a much more powerful regularization technique than the popular alternatives, e.g. dropout or L2-norm penalty. The addition of those methods improved individual models, but was hurting the performance of the ensemble. The core property of an ensemble is diversity. We build an ensemble using several sources of diversity. First, the ensemble models are fit on three different metrics: sMAPE, MASE and MAPE, a version of sMAPE that has only the ground truth value in the denominator. Second, for every horizon H, individual models are trained on input windows of different length: 2H, 3H, . . . , 7H, for a total of six window lengths. Thus the overall ensemble exhibits a multi-scale aspect. Finally, we perform a bagging procedure (Breiman, 1996) by including models trained with different random initializations. We use 180 total models to report results on the test set (please refer to Appendix B for the ablation of ensemble size). We use the median as ensemble aggregation function. 

4. RELATED WORK:

The approaches to TS forecasting can be split in a few distinct categories. The statistical modeling approaches based on exponential smoothing and its different flavors are well established and are often considered a default choice in the industry (Holt, 1957; 2004; Winters, 1960) . More advanced variations of exponential smoothing include the winner of M3 competition, the Theta method (Assimakopoulos & Nikolopoulos, 2000) that decomposes the forecast into several theta-lines and statistically combines them. The pinnacle of the statistical approach encapsulates ARIMA, auto-ARIMA and in general, the unified state-space modeling approach, that can be used to explain and analyze all of the approaches mentioned above (see Hyndman & Khandakar (2008) for an overview). More recently, ML/TS combination approaches started infiltrating the domain with great success, showing promising results by using the outputs of statistical engines as features. In fact, 2 out of top-5 entries in the M4 competition are approaches of this type, including the second entry (Montero-Manso et al., 2019) . The second entry computes the outputs of several statistical methods on the M4 dataset and combines them using gradient boosted tree (Chen & Guestrin, 2016) . Somewhat independently, the work in the modern deep learning TS forecasting developed based on variations of recurrent neural networks (Flunkert et al., 2017; Rangapuram et al., 2018b; Toubeau et al., 2019; Zia & Razzaq, 2018) being largely dominated by the electricity load forecasting in the multi-variate setup. A few earlier works explored the combinations of recurrent neural networks with dilation, residual connections and attention (Chang et al., 2017; Kim et al., 2017; Qin et al., 2017) . These served as a basis for the winner of the M4 competition (Smyl, 2020) . The winning entry combines a Holt-Winters style seasonality model with its parameters fitted to a given TS via gradient descent and a unique combination of dilation/residual/attention approaches for each forecast horizon. The resulting model is a hybrid model that architecturally heavily relies on a time-series engine. It is hand crafted to each specific horizon of M4, making this approach hard to generalize to other datasets.

5. EXPERIMENTAL RESULTS:

Our key empirical results based on aggregate performance metrics over several datasets-M4 ( M4 Team, 2018b; Makridakis et al., 2018b) , M3 (Makridakis & Hibon, 2000; Makridakis et al., 2018a) and TOURISM (Athanasopoulos et al., 2011) -appear in Table 1 . More detailed descriptions of the datasets are provided in Section 5.1 and Appendix A. For each dataset, we compare our results with best 5 entries for this dataset reported in the literature, according to the customary metrics specific to each dataset (M4: OWA and sMAPE, M3: sMAPE, TOURISM: MAPE). More granular dataset-specific results with data splits over forecast horizons and types of time series appear in respective appendices (M4: Appendix C.1; M3: Appendix C.2; TOURISM: Appendix C.3).

In Table 1 , we study the performance of two N-BEATS configurations: generic (N-BEATS-G) and interpretable (N-BEATS-I), as well as N-BEATS-I+G (ensemble of all models from N-BEATS-G and N-BEATS-I). On M4 dataset, we compare against 5 representatives from the M4 competition (Makri-dakis et al., 2018b) : each best in their respective model class. Pure ML is the submission by B. Trotta, the best entry among the 6 pure ML models. Statistical is the best pure statistical model by N.Z. Legaki and K. Koutsouri. ML/TS combination is the model by P. Montero-Manso, T. Talagala, R.J. Hyndman and G. Athanasopoulos, second best entry, gradient boosted tree over a few statistical time series models. ProLogistica is the third entry in M4 based on the weighted ensemble of statistical methods. Finally, DL/TS hybrid is the winner of M4 competition (Smyl, 2020) . On the M3 dataset, we compare against the Theta method ( Assimakopoulos & Nikolopoulos, 2000) , the winner of M3; DOTA, a dynamically optimized Theta model (Fiorucci et al., 2016) ; EXP, the most resent statistical approach and the previous state-of-the-art on M3 (Spiliotis et al., 2019) ; as well as ForecastPro, an off-the-shelf forecasting software that is based on model selection between exponential smoothing, ARIMA and moving average (Athanasopoulos et al., 2011; Assimakopoulos & Nikolopoulos, 2000) .

On the TOURISM dataset, we compare against 3 statistical benchmarks ( According to Table 1 , N-BEATS demonstrates state-of-the-art performance on three challenging non-overlapping datasets containing time series from very different domains, sampling frequencies and seasonalities. As an example, on M4 dataset, the OWA gap between N-BEATS and the M4 winner (0.821 -0.795 = 0.026) is greater than the gap between the M4 winner and the second entry (0.838 -0.821 = 0.017). Generic N-BEATS model uses as little prior knowledge as possible, with no feature engineering, no scaling and no internal architectural components that may be considered TS-specific. Thus the result in Table 1 leads us to the conclusion that DL does not need support from the statistical approaches or hand-crafted feature engineering and domain knowledge to perform extremely well on a wide array of TS forecasting tasks. On top of that, the proposed general architecture performs very well on three different datasets outperforming a wide variety of models, both generic and manually crafted to respective dataset, including the winner of M4, a model architecturally adjusted by hand to each forecast-horizon subset of the M4 data.

5.1 DATASETS Team, 2018b; Makridakis et al., 2018b) is the latest in an influential series of forecasting competitions organized by Spyros Makridakis since 1982 (Makridakis et al., 1982) . The 100k-series dataset is large and diverse, consisting of data frequently encountered in business, financial and economic forecasting, and sampling frequencies ranging from hourly to yearly. A table with summary statistics is presented in Appendix A.1, showing wide variability in TS characteristics.

M4 (M4

M3 (Makridakis & Hibon, 2000) is similar in its composition to M4, but has a smaller overall scale (3003 time series total vs. 100k in M4). A table with summary statistics is presented in Appendix A.2. Over the past 20 years, this dataset has supported significant efforts in the design of more optimal statistical models, e.g. Theta and its variants (Assimakopoulos & Nikolopoulos, 2000; Fiorucci et al., 2016; Spiliotis et al., 2019) . Furthermore, a recent publication (Makridakis et al., 2018a) based on a subset of M3 presented evidence that ML models are inferior to the classical statistical models.

TOURISM (Athanasopoulos et al., 2011) dataset was released as part of the respective Kaggle competition conducted by Athanasopoulos & Hyndman (2011) . The data include monthly, quarterly and yearly series supplied by both governmental tourism organizations (e.g. Tourism Australia, the Hong Kong Tourism Board and Tourism New Zealand) as well as various academics, who had used them in previous studies. A table with summary statistics is presented in Appendix A.3.

5.2. TRAINING METHODOLOGY:

We split each dataset into train, validation and test subsets. The test subset is the standard test set previously defined for each dataset (M4 Team, 2018a; Makridakis & Hibon, 2000; Athanasopoulos et al., 2011) . The validation and train subsets for each dataset are obtained by splitting their full train sets at the boundary of the last horizon of each time series. We use the train and validation subsets to tune hyperparameters. Once the hyperparameters are determined, we train the model on the full train set and report results on the test set. Please refer to Appendix D for detailed hyperparameter settings at the block level. N-BEATS is implemented and trained in Tensorflow (Abadi et al., 2015) . We share parameters of the network across horizons, therefore we train one model per horizon for each dataset. If every time series is interpreted as a separate task, this can be linked back to the multitask learning and furthermore to meta-learning (see discussion in Section 6), in which a neural network is regularized by learning on multiple tasks to improve generalization. We would like to stress that models for different horizons and datasets reuse the same architecture. Architectural hyperparameters (width, number of layers, number of stacks, etc.) are fixed to the same values across horizons and across datasets (see Appendix D). The fact that we can reuse architecture and even hyperparameters across horizons indicates that the proposed architecture design generalizes well across time series of different nature. The same architecture is successfully trained on the M4 Monthly subset with 48k time series and the M3 Others subset with 174 time series. This is a much stronger result than e.g. the result of S. Smyl (Makridakis et al., 2018b) who had to use very different architectures hand crafted for different horizons.

To update network parameters for one horizon, we sample train batches of fixed size 1024. We pick 1024 TS ids from this horizon, uniformly at random with replacement. For each selected TS id we pick a random forecast point from the historical range of length L H immediately preceding the last point in the train part of the TS. L H is a cross-validated hyperparameter. We observed that for subsets with large number of time series it tends to be smaller and for subsets with smaller number of time series it tends to be larger. For example, in massive Yearly, Monthly, Quarterly subsets of M4 L H is equal to 1.5; and in moderate to small Weekly, Daily, Hourly subsets of M4 L H is equal to 10. Given a sampled forecast point, we set one horizon worth of points following it to be the target forecast window y and we set the history of points of one of lengths 2H, 3H, . . . , 7H preceding it to be the input x to the network. We use the Adam optimizer with default settings and initial learning rate 0.001. While optimising the ensemble members relying on the minimization of sMAPE metric, we stop the gradient flows in the denominator to make training numerically stable. The neural network training is run with early stopping and the number of batches is determined on the validation set. The GPU based training of one ensemble member for entire M4 dataset takes between 30 min and 2 hours depending on neural network settings and hardware.

5.3. INTERPRETABILITY RESULTS:

Fig. 2 studies the outputs of the proposed model in the generic and the interpretable configurations.

As discussed in Section 3.3, to make the generic architecture presented in Fig. 1 interpretable, we constrain g θ in the first stack to have the form of polynomial (2) while the second one has the form of Fourier basis (3). Furthermore, we use the outputs of the generic configuration of N-BEATS as control group (the generic model of 30 residual blocks depicted in Fig. 1 is divided into two stacks) and we plot both generic (suffix "-G") and interpretable (suffix "-I") stack outputs side by side in Fig. 2 . The outputs of generic model are arbitrary and non-interpretable: either trend or seasonality or both of them are present at the output of both stacks. The magnitude of the output (peak-to-peak) is generally smaller at the output of the second stack. The outputs of the interpretable model exhibit distinct properties: the trend output is monotonic and slowly moving, the seasonality output is regular, cyclical and has recurring fluctuations. The peak-to-peak magnitude of the seasonality output is significantly larger than that of the trend, if significant seasonality is present in the time series.

Similarly, the peak-to-peak magnitude of trend output tends to be small when no obvious trend is present in the ground truth signal. Thus the proposed interpretable architecture decomposes its forecast into two distinct components. Our conclusion is that the outputs of the DL model can be made interpretable by encoding a sensible inductive bias in the architecture. Table 1 confirms that this does not result in performance drop.

6 DISCUSSION: CONNECTIONS TO META-LEARNING Meta-learning defines an inner learning procedure and an outer learning procedure. The inner learning procedure is parameterized, conditioned or otherwise influenced by the outer learning procedure (Bengio et al., 1991) . The prototypical inner vs. outer learning is individual learning in the lifetime of an animal vs. evolution of the inner learning procedure itself over many generations of individuals. To see the two levels, it often helps to refer to two sets of parameters, the inner parameters (e.g. synaptic weights) which are modified inside the inner learning procedure, and the outer parameters or meta-parameters (e.g. genes) which get modified only in the outer learning procedure.

N-BEATS can be cast as an instance of meta-learning by drawing the following parallels. The outer learning procedure is encapsulated in the parameters of the whole network, learned by gradient descent. The inner learning procedure is encapsulated in the set of basic building blocks and modifies the expansion coefficients θ f that basis g f takes as inputs. The inner learning proceeds through a sequence of stages, each corresponding to a block within the stack of the architecture. Each of the blocks can be thought of as performing the equivalent of an update step which gradually modifies the expansion coefficients θ f which eventually feed into g f in each block (which get added together to form the final prediction). The inner learning procedure takes a single history from a piece of a TS and sees that history as a training set. It produces forward expansion coefficients θ f (see Fig. 1 ), which parametrically map inputs to predictions. In addition, each preceding block modifies the input to the next block by producing backward expansion coefficients θ b , thus conditioning the learning and the output of the next block. In the case of the interpretable model, the meta-parameters are only in the FC layers because the g f 's are fixed. In the case of the generic model, the meta-parameters also include the V's which define the g f non-parametrically. This point of view is further reinforced by the results of the ablation study reported in Appendix B showing that increasing the number of blocks in the stack, as well as the number of stacks improves generalization performance, and can be interpreted as more iterations of the inner learning procedure.

7. CONCLUSIONS:

We proposed and empirically validated a novel architecture for univariate TS forecasting. We showed that the architecture is general, flexible and it performs well on a wide array of TS forecasting problems. We applied it to three non-overlapping challenging competition datasets: M4, M3 and TOURISM and demonstrated state-of-the-art performance in two configurations: generic and interpretable. This allowed us to validate two important hypotheses: (i) the generic DL approach performs exceptionally well on heterogeneous univariate TS forecasting problems using no TS domain knowledge, (ii) it is viable to additionally constrain a DL model to force it to decompose its forecast into distinct human interpretable outputs. We also demonstrated that the DL models can be trained on multiple time series in a multi-task fashion, successfully transferring and sharing individual learnings. We speculate that N-BEATS's performance can be attributed in part to it carrying out a form of meta-learning, a deeper investigation of which should be the subject of future work. A DATASET DETAILS

A.1 M4 DATASET DETAILS Table 2 outlines the composition of the M4 dataset across domains and forecast horizons by listing the number of time series based on their frequency and type (M4 Team, 2018b). The M4 dataset is large and diverse: all forecast horizons are composed of heterogeneous time series types (with exception of Hourly) frequently encountered in business, financial and economic forecasting. Summary statistics on series lengths are also listed, showing wide variability therein, as well as a characterization (smooth vs erratic) that follows Syntetos et al. (2005) , and is based on the squared coefficient of variation of the series. All series have positive observed values at all time-steps; as such, none can be considered intermittent or lumpy per Syntetos et al. (2005) .

A.2 M3 DATASET DETAILS:

Table 3 outlines the composition of the M3 dataset across domains and forecast horizons by listing the number of time series based on their frequency and type (Makridakis & Hibon, 2000) . The M3 is smaller than the M4, but it is still large and diverse: all forecast horizons are composed of heterogeneous time series types frequently encountered in business, financial and economic forecasting. Summary statistics on series lengths are also listed, showing wide variability in length, as well as a characterization (smooth vs erratic) that follows Syntetos et al. (2005) , and is based on the squared coefficient of variation of the series. All series have positive observed values at all time-steps; as such, none can be considered intermittent or lumpy per Syntetos et al. (2005) .

A.3 TOURISM DATASET DETAILS Table 4 outlines the composition of the TOURISM dataset across forecast horizons by listing the number of time series based on their frequency. Summary statistics on series lengths are listed, showing wide variability in length. All series have positive observed values at all time-steps. In contrast to M4 and M3 datasets, TOURISM includes a much higher fraction of erratic series. We performed an ablation study on the validation set, using sMAPE metric as performance criterion.

We addressed two specific questions with this study. First, Is stacking layers helpful? Second, Does the architecture based on the combination of layers with different basis functions results in better performance than the architecture using only one layer type?

Layer stacking. We start our study with the generic architecture that consists of stacks of one residual block of 5 FC layers each of the form Fig. 1 and we increase the number of stacks. Results presented in Table 5 confirm that increasing the number of stacks decreases error and at certain point the gain saturates. We would like to mention that the network having 30 stack of depth 5 is in fact a very deep network of total depth 150 layers.

Basis synergy. Stacking works well for the interpretable architecture as can be seen in Table 6 depicting the results of ablating the interpretable architecture configuration. Here we experiment with the architecture that is composed of 2 stacks, stack one is trend model and stack two is the seasonality model. Each stack has variable number of residual blocks and each residual block has 5 FC layers. We found that this architecture works best when all weights are shared within stack. We clearly see that increasing the number of layers improves performance. The largest network is 60 layers deep. On top of that, we observe that the architecture that consists of stacks based on different basis functions wins over the architecture based on the same stack. It looks like chaining stacks of different nature results in synergistic effects. This is logical as function classes that can be modelled by trend and seasonality stacks have small overlap.

B.2 ENSEMBLE SIZE:

Figure 3 demonstrates that increasing the ensemble size results in improved performance. Most importantly, according to Figure 3 , N-BEATS achieves state-of-the-art performance even if comparatively small ensemble size of 18 models is used. Therefore, computational efficiency of N-BEATS can be traded very effectively for performance and there is no over-reliance of the results on large ensemble size.

B.3 DOUBLY RESIDUAL STACKING:

In Section 3.2 we described the proposed doubly residual stacking (DRESS) principle, which is the topological foundation of N-BEATS. The topology is based on both (i) running a residual backcast connection and (ii) producing partial block-level forecasts that are further aggregated at stack and model levels to produce the final model-level forecast. In this section we conduct a study to confirm the accuracy effectiveness of this topology compared to several alternatives. The methodology underlying this study is that we remove either the backcast or partial forecast links or both and track how this affects the forecasting metrics. We keep the number of parameters in the network for each of the architectural alternatives fixed by using the same number of layers in the network (we used default hyperparameter settings reported in Table 18 ). The architectural alternatives are depicted in Figure 4 and described in detail below. N-BEATS-DRESS is depicted in Fig. 4a . This is the default configuration of N-BEATS using doubly residual stacking described in Section 3.2.

PARALLEL is depicted in Fig. 4b . This is the alternative where the backward residual connection is disabled and the overall model input is fed to every block. The blocks then forecast in parallel using the same input and their individual outputs are summed to make the final forecast.

NO-RESIDUAL is depicted in Fig. 4c . This is the alternative where the backward residual connection is disabled. Unlike PARALLEL, in this case the backcast forecast of the previous block is fed as input to the next block. Unlike the usual feed-forward network, in the NO-RESIDUAL architecture, each block makes a partial forecast and their individual outputs are summed to make the final forecast.

LAST-FORWARD is depicted in Fig. 4d . This is the alternative where the backward residual connection is active, however the model level forecast is derived only from the last block. So, the partial forward forecasts are disabled. This is the architecture that is closest to the classical residual network.

NO-RESIDUAL-LAST-FORWARD is depicted in Fig. 4f . This is the alternative where both backward residual and the partial forward connections are disabled. This is therefore a simple feed-forward network, but very deep.

The quantitative ablation study results on the M4 dataset are reported in Tables 7 8 9 10 . N-BEATS-DRESS model is essentially N-BEATS model in this study. For this study we used ensemble size of 18. Since the ensemble size is 18 for N-BEATS-DRESS, as opposed to 180 used for N-BEATS, the OWA metric reported in Table 9 for N-BEATS-DRESS is higher than the OWA reported for N-BEATS-G in Table 12 . Note that both results align well with OWA reported in Figure 3 for different ensemble sizes, as part of the ensemble size ablation conducted in Section B.2.

The results presented in Tables 7-10 demonstrate that the doubly residual stacking topology provides a clear overall advantage over the alternative architectures in which either backcast residual links or the partial forward forecast links are disabled. 

C DETAILED EMPIRICAL RESULTS:

C.1 DETAILED RESULTS: M4 DATASET Tables 11 and 12 present our key quantitative empirical results showing that the proposed model achieves the state of the art performance on the challenging M4 benchmark. We study the performance of two model configurations: generic (Ours-G) and interpretable (Ours-I), as well as Ours-I+G (ensemble of all models from Ours-G and Ours-I). We compare against 4 representatives from the M4 competition: each best in their respective model class. Best pure ML is the submission by B. Trotta, the best entry among the 6 pure ML models. Best statistical is the best pure statistical model by N.Z. Legaki and K. Koutsouri. Best ML/TS combination is the model by P. Montero-Manso, T. Talagala, R.J. Hyndman and G. Athanasopoulos, second best entry, gradient boosted tree over a few statistical time series models. Finally, DL/TS hybrid is the winner of M4 competition (Smyl, 2020) .

N-BEATS outperforms all other approaches on all the studied subsets of time series. The average OWA gap between our generic model and the M4 winner (0.821 -0.795 = 0.026) is greater than the gap between the M4 winner and the second entry (0.838 -0.821 = 0.017).

A more granular and detailed statistical analysis of our results on M4 is provided in Table 13 . This table first presents the sMAPE for N-BEATS, decomposed by M4 time series sub-type and sampling frequency (upper part). Then (lower part), it shows the average sMAPE difference between the N-BEATS results and the M4 winner (TS/DL hybrid by S. Smyl), adding the standard error of that difference (in parentheses); bold entries indicate statistical significance at the 99% level based on a two-sided paired t-test.

We note that each cross-section of the M4 dataset into horizon and type may be regarded as an independent mini-dataset. We observe that over those mini-datasets there is a preponderance of statistically significant differences between N-BEATS and Smyl (18 cases out of 31) to the advantage of N-BEATS. This provides evidence that (i) the improvement observed on average in Tables 11 and 12 is statistically significant and consistent over smaller subsets of M4 and (ii) N-BEATS generalizes well over time series of different types and sampling frequencies. 

For some of the methods, either average sMAPE was not reported or sMAPE for some of the splits was not reported in their respective publications. Below, we list those cases. BaggedETS.BC (Bergmeir et al., 2016) has not reported numbers on Others.

LGT (Smyl & Kuber, 2016) did not report results on Monthly and Quarterly data. According to the authors, the underlying RNN had problems dealing with raw seasonal data, the ETS based pre-processing was not effective and the LGT pre-processing was not computationally feasible given comparatively large number of time series and their comparatively large length (Smyl & Kuber, 2016) . Finally, EXP (Spiliotis et al., 2019) reported average performance computed using a different methodology than the default M3 and M4 methodology (source: personal communication with the authors). For the latter method we recomputed the Average sMAPE based on the previously reported Yearly, Quarterly and Monthly splits. To calculate it, we follow the M3, M4 and TOURISM competition methodology and compute the average metric as the average over all time series and over all forecast horizons. Given the performance metric values aggregated over Yearly, Quarterly and Monthly splits, the average can be computed straightforwardly as:

sMAPE Average = N Year N Tot sMAPE Year + N Quart N Tot sMAPE Quart + N Month N Tot sMAPE Month + N Others N Tot sMAPE Others .

(5)

Here N Tot = N Year + N Quart + N Month + N Others and N Year = 6 × 645, N Quart = 8 × 756, N Month = 18 × 1428, N Others = 8 × 174. It is clear that for each split, its N is the product of its respective number of time series and its largest forecast horizon. (Athanasopoulos & Hyndman, 2011) . Some of the participants chose to take part only in the second part. Therefore, In addition to entries present in Table 1 , we report competitors from (Athanasopoulos & Hyndman, 2011) that have missing results in Yearly competition. In particular, SaliMali team is the winner of the Quarterly/Monthly time series forecasting competition (Brierley, 2011) . Their approach is based on a weighted ensemble of statistical methods.

Teams Robert and Idalgo used unknown approaches. We can see from Table 15 that N-BEATS achieves state-of-the-art performance on all subsets of TOURISM dataset. On average, it is state of the art and it gains 4.2% over the best-known approach LeeCBaker, and 11.5% over auto-ARIMA.

The average metrics have not been reported in the original competition results (Athanasopoulos et al., 2011; Athanasopoulos & Hyndman, 2011) . Therefore, in Table 15 , we present the Average MAPE metric calculated by us based on the previously reported Yearly, Quarterly and Monthly splits. To calculate it, we follow the M4 competition methodology and compute the average metric as the average over all time series and over all forecast horizons. Given the performance metric values aggregated over Yearly, Quarterly and Monthly splits, the average can be computed straightforwardly as:

EQUATION

Here N Tot = N Year + N Quart + N Month and N Year = 4 × 518, N Quart = 8 × 427, N Month = 24 × 366. It is clear that for each split, its N is the product of its respective number of time series and its largest forecast horizon.

C.5 DETAILED RESULTS: COMPARE TO DEEPAR, DEEP STATE SPACE MODELS N-BEATS-I and N-BEATS-G correspond to the interpretable and generic model configurations defined in Section 3.3.

D.1 COMMON PARAMETERS:

L H is the coefficient defining the length of training history immediately preceding the last point in the train part of the TS that is used to generate training samples. For example, if for M4 Yearly the forecast horizon is 6 and L H is 1.5, then we consider 1.5 • 6 = 9 most recent points in the train dataset for each time series to generate training samples. A training sample from a given TS in M4 Yearly is then generated by choosing one of the most recent 9 points as an anchor. All the points preceding the anchor are used to create the input to N-BEATS, while the points following and including the anchor become training target. Target and history points that fall outside of the time series limits given the anchor position are filled with zeros and masked during the training. We observed that for subsets with large number of time series L H tends to be smaller and for subsets with smaller number of time series it tends to be larger. For example, in massive Yearly, Monthly, Quarterly subsets of M4 L H is equal to 1.5; and in moderate to small Weekly, Daily, Hourly subsets of M4 L H is equal to 10.

Iterations is the number of batches used to train N-BEATS.

Losses is the set of loss functions that is used to build ensemble. We observed on the respective validation sets that for M4 and M3 mixing models trained on a variety of metrics resulted in performance gain. In the case of TOURISM dataset training only on MAPE led to the best validation scores.

Sharing defines whether the coefficients in the fully-connected layers are shared. We observed that the interpretable model works best when weights are shared across stack, while generic model works best when none of the weights are shared.

Lookback period is the length of the history window forming the input to the model (please refer to Figure 1 ). This is the function of the forecast horizon length, H. In our experiments we mixed models with lookback periods 2H, 3H, 4H, 5H, 6H, 7H in one ensemble. As an example, for a forecast horizon length H = 8 and a lookback period 7H, the model's input will consist of the history window of 7 • 8 = 56 samples.

Batch is the batch size. We used batch size of 1024. We observed that the training was faster with larger batch sizes, however in our setup little gain was observed with batch sizes beyond 1024.

D.2 N-BEATS-I PARAMETERS:

S-width is the width of the fully connected layers in the blocks comprising the seasonality stack of the interpretable model (please refer to Figure 1 ).

S-blocks is the number of blocks comprising the seasonality stack of the interpretable model (please refer to Figure 1 ).

S-block-layers is the number of fully-connected layers comprising one block in the seasonality stack of the interpretable model (preceding the final fully-connected projection layers forming the backcast/forecast fork, please refer to Figure 1 ).

T-width is the width of the fully connected layers in the blocks comprising the trend stack of the interpretable model (please refer to Figure 1 ).

T-degree is the degree p of polynomial in the trend stack of the interpretable model (please refer to equation ( 2)).

T-blocks is the number of blocks comprising the trend stack of the interpretable model (please refer to Figure 1 ).

T-block-layers is the number of fully-connected layers comprising one block in the trend stack of the interpretable model (preceding the final fully-connected projection layers forming the backcast/forecast fork, please refer to Figure 1 ).

D.3 N-BEATS-G PARAMETERS:

Width is the width of the fully connected layers in the blocks comprising the stacks of the generic model (please refer to Figure 1 ).

Blocks is the number of blocks comprising the stack of the generic model (please refer to Figure 1 ).

Block-layers is the number of fully-connected layers comprising one block in the stack of the generic model (preceding the final fully-connected projection layers forming the backcast/forecast fork, please refer to Figure 1 ). First, we reproduce Fig. 5 for the convenience of the reader. Second, for each row in the figure, we produce a table showing the numeric values of each signal depicted in corresponding plots (please refer to Tables 19 20 21 22 23 24 ). We make sure that the names of signals in figure legends and in the table columns match, such that they can easily be cross-referenced. It can be clearly seen in Tables 19-24 that (i) traces STACK1-I and STACK2-I sum up to trace FORECAST-I, (ii) traces STACK1-G and STACK2-G sum up to trace FORECAST-G, (iii) traces FORECAST-I and FORECAST-G are overall very similar even though their components may significantly differ from each other.

:

With minor differences compared to the sMAPE definition used for M4. Please refer to Appendix A in(Makridakis & Hibon,

2000) for the mathematical definition.

https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014

https://archive.ics.uci.edu/ml/datasets/PEMS-SF

https://github.com/rofuyu/exp-trmf-nips16/blob/master/python/exp-scripts/datasets/download-data.sh

Figures and Tables:

y T +i | |y T +i | + | y T +i | ,

Figure 2: The outputs of generic and the interpretable configurations, M4 dataset. Each row is one time series example per data frequency, top to bottom (Yearly: id Y3974, Quarterly: id Q11588, Monthly: id M19006, Weekly: id W246, Daily: id D404, Hourly: id H344). The magnitudes in a row are normalized by the maximal value of the actual time series for convenience. Column (a) shows the actual values (ACTUAL), the generic model forecast (FORECAST-G) and the interpretable model forecast (FORECAST-I). Columns (b) and (c) show the outputs of stacks 1 and 2 of the generic model, respectively; FORECAST-G is their summation. Columns (d) and (e) show the output of the Trend and the Seasonality stacks of the interpretable model, respectively; FORECAST-I is their summation.

Figure 3: M4 test performance (OWA) as a function of ensemble size, based on N-BEATS-G. This figure shows that N-BEATS loses less than 0.5% in terms of OWA performance even if 10 times smaller ensemble size is used.

Figure 4: The architectural configurations used in the ablation study of the doubly residual stack. Symbol denotes unconnected output.

Figure 5: The outputs of generic and the interpretable configurations, M4 dataset. Each row is one time series example per data frequency, top to bottom (Yearly: id Y3974, Quarterly: id Q11588, Monthly: id M19006, Weekly: id W246, Daily: id D404, Hourly: id H344). The magnitudes in a row are normalized by the maximal value of the actual time series for convenience. Column (a) shows the actual values (ACTUAL), the generic model forecast (FORECAST-G) and the interpretable model forecast (FORECAST-I). Columns (b) and (c) show the outputs of stacks 1 and 2 of the generic model, respectively; FORECAST-G is their summation. Columns (d) and (e) show the output of the Trend and the Seasonality stacks of the interpretable model, respectively; FORECAST-I is their summation.

Performance on the M4, M3, TOURISM test sets, aggregated over each dataset. Evaluation metrics are specified for each dataset; lower values are better. The number of time series in each dataset is provided in brackets.

Table content: <table><tr><td colspan="2">M4 Average (100,000)</td><td colspan="2">M3 Average (3,003)</td><td colspan="2">TOURISM Average (1,311)</td></tr><tr><td/><td>sMAPE OWA</td><td/><td>sMAPE</td><td/><td>MAPE</td></tr><tr><td>Pure ML</td><td colspan="3">12.894 0.915 Comb S-H-D 13.52</td><td>ETS</td><td>20.88</td></tr><tr><td>Statistical</td><td colspan="2">11.986 0.861 ForecastPro</td><td>13.19</td><td>Theta</td><td>20.88</td></tr><tr><td>ProLogistica</td><td colspan="2">11.845 0.841 Theta</td><td>13.01</td><td>ForePro</td><td>19.84</td></tr><tr><td colspan="3">ML/TS combination 11.720 0.838 DOTM</td><td>12.90</td><td>Stratometrics</td><td>19.52</td></tr><tr><td>DL/TS hybrid</td><td colspan="2">11.374 0.821 EXP</td><td>12.71</td><td>LeeCBaker</td><td>19.35</td></tr><tr><td>N-BEATS-G</td><td>11.168 0.797</td><td/><td>12.47</td><td/><td>18.47</td></tr><tr><td>N-BEATS-I</td><td>11.174 0.798</td><td/><td>12.43</td><td/><td>18.97</td></tr><tr><td>N-BEATS-I+G</td><td>11.135 0.795</td><td/><td>12.37</td><td/><td>18.52</td></tr></table>

Athanasopoulos et al.,  2011): ETS, exponential smoothing with cross-validated additive/multiplicative model; Theta method; ForePro, same as ForecastPro in M3; as well as top 2 entries from the TOURISM Kaggle competition(Athanasopoulos & Hyndman, 2011): Stratometrics, an unknown technique; LeeCBaker(Baker & Howard, 2011), a weighted combination of Naïve, linear trend model, and exponentially weighted least squares regression trend.

Table content: <table/>

Composition of the M4 dataset: the number of time series based on their sampling frequency and type.

Table content: <table><tr><td/><td/><td/><td colspan="2">Frequency / Horizon</td><td/><td/><td/></tr><tr><td>Type</td><td colspan="6">Yearly/6 Qtly/8 Monthly/18 Wkly/13 Daily/14 Hrly/48</td><td>Total</td></tr><tr><td>Demographic</td><td>1,088</td><td>1,858</td><td>5,728</td><td>24</td><td>10</td><td>0</td><td>8,708</td></tr><tr><td>Finance</td><td>6,519</td><td>5,305</td><td>10,987</td><td>164</td><td>1,559</td><td>0</td><td>24,534</td></tr><tr><td>Industry</td><td>3,716</td><td>4,637</td><td>10,017</td><td>6</td><td>422</td><td>0</td><td>18,798</td></tr><tr><td>Macro</td><td>3,903</td><td>5,315</td><td>10,016</td><td>41</td><td>127</td><td>0</td><td>19,402</td></tr><tr><td>Micro</td><td>6,538</td><td>6,020</td><td>10,975</td><td>112</td><td>1,476</td><td>0</td><td>25,121</td></tr><tr><td>Other</td><td>1,236</td><td>865</td><td>277</td><td>12</td><td>633</td><td>414</td><td>3,437</td></tr><tr><td>Total</td><td colspan="2">23,000 24,000</td><td>48,000</td><td>359</td><td>4,227</td><td colspan="2">414 100,000</td></tr><tr><td>Min. Length</td><td>19</td><td>24</td><td>60</td><td>93</td><td>107</td><td>748</td><td/></tr><tr><td>Max. Length</td><td>841</td><td>874</td><td>2812</td><td>2610</td><td>9933</td><td>1008</td><td/></tr><tr><td>Mean Length</td><td>37.3</td><td>100.2</td><td>234.3</td><td>1035.0</td><td>2371.4</td><td>901.9</td><td/></tr><tr><td>SD Length</td><td>24.5</td><td>51.1</td><td>137.4</td><td>707.1</td><td>1756.6</td><td>127.9</td><td/></tr><tr><td>% Smooth</td><td>82%</td><td>89%</td><td>94%</td><td>84%</td><td>98%</td><td>83%</td><td/></tr><tr><td>% Erratic</td><td>18%</td><td>11%</td><td>6%</td><td>16%</td><td>2%</td><td>17%</td><td/></tr></table>

Composition of the M3 dataset: the number of time series based on their sampling frequency and type.

Table content: <table><tr><td/><td/><td colspan="2">Frequency / Horizon</td><td/><td/></tr><tr><td>Type</td><td colspan="3">Yearly/6 Quarterly/8 Monthly/18</td><td>Other/8</td><td>Total</td></tr><tr><td>Demographic</td><td>245</td><td>57</td><td>111</td><td>0</td><td>413</td></tr><tr><td>Finance</td><td>58</td><td>76</td><td>145</td><td>29</td><td>308</td></tr><tr><td>Industry</td><td>102</td><td>83</td><td>334</td><td>0</td><td>519</td></tr><tr><td>Macro</td><td>83</td><td>336</td><td>312</td><td>0</td><td>731</td></tr><tr><td>Micro</td><td>146</td><td>204</td><td>474</td><td>4</td><td>828</td></tr><tr><td>Other</td><td>11</td><td>0</td><td>52</td><td>141</td><td>204</td></tr><tr><td>Total</td><td>645</td><td>756</td><td>1,428</td><td>174</td><td>3,003</td></tr><tr><td>Min. Length</td><td>20</td><td>24</td><td>66</td><td>71</td><td/></tr><tr><td>Max. Length</td><td>47</td><td>72</td><td>144</td><td>104</td><td/></tr><tr><td>Mean Length</td><td>28.4</td><td>48.9</td><td>117.3</td><td>76.6</td><td/></tr><tr><td>SD Length</td><td>9.9</td><td>10.6</td><td>28.5</td><td>10.9</td><td/></tr><tr><td>% Smooth</td><td>90%</td><td>99%</td><td>98%</td><td>100%</td><td/></tr><tr><td>% Erratic</td><td>10%</td><td>1%</td><td>2%</td><td>0%</td><td/></tr></table>

Composition of the TOURISM dataset: the number of time series based on their sampling frequency.

Table content: <table><tr><td/><td/><td>Frequency / Horizon</td><td/></tr><tr><td/><td colspan="3">Yearly/4 Quarterly/8 Monthly/24</td><td>Total</td></tr><tr><td/><td>518</td><td>427</td><td>366</td><td>1,311</td></tr><tr><td>Min. Length</td><td>11</td><td>30</td><td>91</td></tr><tr><td>Max. Length</td><td>47</td><td>130</td><td>333</td></tr><tr><td>Mean Length</td><td>24.4</td><td>99.6</td><td>298</td></tr><tr><td>SD Length</td><td>5.5</td><td>20.3</td><td>55.7</td></tr><tr><td>% Smooth</td><td>77%</td><td>61%</td><td>49%</td></tr><tr><td>% Erratic</td><td>23%</td><td>39%</td><td>51%</td></tr></table>

Performance on the M4 test set, sMAPE. Lower values are better. The results are obtained on the ensemble of 18 generic models.

Table content: <table><tr><td/><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td></tr><tr><td>PARALLEL-G</td><td>13.279</td><td>9.558</td><td>12.510</td><td>3.691</td><td>11.538</td></tr><tr><td>NO-RESIDUAL-G</td><td>13.195</td><td>9.555</td><td>12.451</td><td>3.759</td><td>11.493</td></tr><tr><td>LAST-FORWARD-G</td><td>13.200</td><td>9.322</td><td>12.352</td><td>3.703</td><td>11.387</td></tr><tr><td colspan="2">NO-RESIDUAL-LAST-FORWARD-G 15.386</td><td>11.346</td><td>15.282</td><td>6.673</td><td>13.931</td></tr><tr><td>RESIDUAL-INPUT-G</td><td>13.264</td><td>9.545</td><td>12.316</td><td>3.692</td><td>11.438</td></tr><tr><td>N-BEATS-DRESS-G</td><td>13.211</td><td>9.217</td><td>12.122</td><td>3.636</td><td>11.251</td></tr></table>

Performance on the M4 test set, sMAPE. Lower values are better. The results are obtained on the ensemble of 18 interpretable models.

Table content: <table><tr><td/><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td></tr><tr><td>PARALLEL-I</td><td>13.207</td><td>9.530</td><td>12.500</td><td>3.710</td><td>11.510</td></tr><tr><td>NO-RESIDUAL-I</td><td>13.075</td><td>9.707</td><td>12.708</td><td>4.007</td><td>11.637</td></tr><tr><td>LAST-FORWARD-I</td><td>13.168</td><td>9.547</td><td>12.111</td><td>3.599</td><td>11.313</td></tr><tr><td colspan="2">NO-RESIDUAL-LAST-FORWARD-I 13.067</td><td>10.207</td><td>15.177</td><td>4.912</td><td>12.986</td></tr><tr><td>RESIDUAL-INPUT-I</td><td>13.104</td><td>9.716</td><td>12.814</td><td>4.005</td><td>11.697</td></tr><tr><td>N-BEATS-DRESS-I</td><td>13.155</td><td>9.286</td><td>12.009</td><td>3.642</td><td>11.201</td></tr></table>

Performance on the M4 test set, OWA. Lower values are better. The results are obtained on the ensemble of 18 generic models.

Table content: <table><tr><td/><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td></tr><tr><td>PARALLEL-G</td><td>0.780</td><td>0.832</td><td>0.852</td><td>0.844</td><td>0.822</td></tr><tr><td>NO-RESIDUAL-G</td><td>0.774</td><td>0.831</td><td>0.851</td><td>0.853</td><td>0.819</td></tr><tr><td>LAST-FORWARD-G</td><td>0.774</td><td>0.808</td><td>0.840</td><td>0.846</td><td>0.811</td></tr><tr><td colspan="2">NO-RESIDUAL-LAST-FORWARD-G 0.948</td><td>1.029</td><td>1.095</td><td>1.296</td><td>1.030</td></tr><tr><td>RESIDUAL-INPUT-G</td><td>0.779</td><td>0.831</td><td>0.840</td><td>0.844</td><td>0.817</td></tr><tr><td>N-BEATS-DRESS-G</td><td>0.776</td><td>0.800</td><td>0.823</td><td>0.835</td><td>0.803</td></tr></table>

Performance on the M4 test set, OWA. Lower values are better. The results are obtained on the ensemble of 18 interpretable models.

Table content: <table><tr><td/><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td></tr><tr><td>PARALLEL-I</td><td>0.776</td><td>0.831</td><td>0.857</td><td>0.845</td><td>0.821</td></tr><tr><td>NO-RESIDUAL-I</td><td>0.769</td><td>0.848</td><td>0.886</td><td>0.886</td><td>0.833</td></tr><tr><td>LAST-FORWARD-I</td><td>0.773</td><td>0.836</td><td>0.825</td><td>0.817</td><td>0.808</td></tr><tr><td colspan="2">NO-RESIDUAL-LAST-FORWARD-I 0.771</td><td>0.900</td><td>1.085</td><td>1.016</td><td>0.922</td></tr><tr><td>RESIDUAL-INPUT-I</td><td>0.771</td><td>0.848</td><td>0.892</td><td>0.887</td><td>0.836</td></tr><tr><td>N-BEATS-DRESS-I</td><td>0.771</td><td>0.805</td><td>0.819</td><td>0.836</td><td>0.800</td></tr></table>

Performance on the M4 test set, sMAPE. Lower values are better. Red -second best.

Table content: <table><tr><td/><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td></tr><tr><td>Best pure ML</td><td>14.397</td><td>11.031</td><td>13.973</td><td>4.566</td><td>12.894</td></tr><tr><td>Best statistical</td><td>13.366</td><td>10.155</td><td>13.002</td><td>4.682</td><td>11.986</td></tr><tr><td>Best ML/TS combination</td><td>13.528</td><td>9.733</td><td>12.639</td><td>4.118</td><td>11.720</td></tr><tr><td colspan="2">DL/TS hybrid, M4 winner 13.176</td><td>9.679</td><td>12.126</td><td>4.014</td><td>11.374</td></tr><tr><td>N-BEATS-G</td><td>13.023</td><td>9.212</td><td>12.048</td><td>3.574</td><td>11.168</td></tr><tr><td>N-BEATS-I</td><td>12.924</td><td>9.287</td><td>12.059</td><td>3.684</td><td>11.174</td></tr><tr><td>N-BEATS-I+G</td><td>12.913</td><td>9.213</td><td>12.024</td><td>3.643</td><td>11.135</td></tr></table>

Performance on the M4 test set, OWA and M4 rank. Lower values are better. Red -second best.

Table content: <table><tr><td/><td colspan="6">Yearly Quarterly Monthly Others Average Rank</td></tr><tr><td/><td>(23k)</td><td>(24k)</td><td>(48k)</td><td>(5k)</td><td>(100k)</td><td/></tr><tr><td>Best pure ML</td><td>0.859</td><td>0.939</td><td>0.941</td><td>0.991</td><td>0.915</td><td>23</td></tr><tr><td>Best statistical</td><td>0.788</td><td>0.898</td><td>0.905</td><td>0.989</td><td>0.861</td><td>8</td></tr><tr><td>Best ML/TS combination</td><td>0.799</td><td>0.847</td><td>0.858</td><td>0.914</td><td>0.838</td><td>2</td></tr><tr><td colspan="2">DL/TS hybrid, M4 winner 0.778</td><td>0.847</td><td>0.836</td><td>0.920</td><td>0.821</td><td>1</td></tr><tr><td>N-BEATS-G</td><td>0.765</td><td>0.800</td><td>0.820</td><td>0.822</td><td>0.797</td><td/></tr><tr><td>N-BEATS-I</td><td>0.758</td><td>0.807</td><td>0.824</td><td>0.849</td><td>0.798</td><td/></tr><tr><td>N-BEATS-I+G</td><td>0.758</td><td>0.800</td><td>0.819</td><td>0.840</td><td>0.795</td><td/></tr></table>

Performance decomposition on non-overlapping subsets of the M4 test set and comparison with the Smyl model results. Average sMAPE difference vs Smyl model, computed as N-BEATS -Smyl.

Table content: <table><tr><td/><td colspan="4">Demographic Finance Industry Macro</td><td>Micro</td><td>Other</td></tr><tr><td colspan="5">sMAPE per M4 series type and sampling frequency</td><td/><td/></tr><tr><td>Yearly</td><td>8.931</td><td>13.741</td><td colspan="4">16.317 13.327 10.489 13.320</td></tr><tr><td>Quarterly</td><td>9.219</td><td>10.787</td><td>8.628</td><td>8.576</td><td>9.264</td><td>6.250</td></tr><tr><td>Monthly</td><td>4.357</td><td>13.353</td><td colspan="4">12.657 12.571 13.627 11.595</td></tr><tr><td>Weekly</td><td>4.580</td><td>3.004</td><td>9.258</td><td colspan="2">7.220 10.425</td><td>6.183</td></tr><tr><td>Daily</td><td>6.351</td><td>3.467</td><td>3.835</td><td>2.525</td><td>2.299</td><td>2.885</td></tr><tr><td>Hourly</td><td/><td/><td/><td/><td/><td>8.197</td></tr><tr><td colspan="4">Standard error of the mean displayed in parenthesis.</td><td/><td/><td/></tr><tr><td colspan="5">Bold entries are significant at the 99% level (2-sided paired t-test).</td><td/><td/></tr><tr><td>Yearly</td><td colspan="6">-0.749 -0.337 -0.065 -0.386 -0.168 -0.157</td></tr><tr><td/><td>(0.119)</td><td>(0.065)</td><td>(0.087)</td><td>(0.085)</td><td>(0.056)</td><td>(0.140)</td></tr><tr><td>Quarterly</td><td colspan="6">-0.651 -0.281 -0.328 -0.712 -0.523 -0.029</td></tr><tr><td/><td>(0.085)</td><td>(0.047)</td><td>(0.043)</td><td>(0.060)</td><td>(0.051)</td><td>(0.083)</td></tr><tr><td>Monthly</td><td colspan="3">-0.185 -0.379 -0.419</td><td>0.089</td><td colspan="2">0.338 -0.279</td></tr><tr><td/><td>(0.023)</td><td>(0.034)</td><td>(0.036)</td><td>(0.039)</td><td>(0.034)</td><td>(0.162)</td></tr><tr><td>Weekly</td><td colspan="6">-0.336 -1.075 -0.937 -1.627 -3.029 -1.193</td></tr><tr><td/><td>(0.270)</td><td>(0.221)</td><td>(1.399)</td><td>(0.770)</td><td>(0.378)</td><td>(0.772)</td></tr><tr><td>Daily</td><td colspan="6">0.191 -0.098 -0.124 -0.026 -0.367 -0.037</td></tr><tr><td/><td>(0.231)</td><td>(0.018)</td><td>(0.025)</td><td>(0.057)</td><td>(0.013)</td><td>(0.015)</td></tr><tr><td>Hourly</td><td/><td/><td/><td/><td/><td>-1.132</td></tr><tr><td/><td/><td/><td/><td/><td/><td>(0.163)</td></tr></table>

Performance on the M3 test set, Average sMAPE, aggregate over all forecast horizons (Yearly: 1-6, Quarterly: 1-8, Monthly: 1-18, Other: 1-8, Average: 1-18). Lower values are better. Red -second best. † Numbers are computed by us.

Table content: <table><tr><td colspan="5">Yearly Quarterly Monthly Others Average</td></tr><tr><td>(645)</td><td>(756)</td><td>(1428)</td><td>(174)</td><td>(3003)</td></tr></table>

Performance on the TOURISM test set, Average MAPE, aggregate over all forecast horizons (Yearly: 1-4, Quarterly: 1-8, Monthly: 1-24, Average: 1-24). Lower values are better. Red -second best.

Table content: <table><tr><td colspan="4">Yearly Quarterly Monthly Average</td></tr><tr><td>(518)</td><td>(427)</td><td>(366)</td><td>(1311)</td></tr></table>

ND (7) performance of DeepAR, DeepState models published in(Rangapuram  et al., 2018a)  and N-BEATS.

Table content: <table><tr><td colspan="4">Table 17: ND Performance of DeepAR, Deep State Space, and N-BEATS models on M4-Hourly and</td></tr><tr><td>TOURISM datasets</td><td/><td/><td/></tr><tr><td/><td colspan="3">M4 (Hourly) TOURISM (Monthly) TOURISM (Quarterly)</td></tr><tr><td>DeepAR</td><td>0.09</td><td>0.107</td><td>0.11</td></tr><tr><td>DeepState</td><td>0.044</td><td>0.138</td><td>0.098</td></tr><tr><td>N-BEATS-G (ours)</td><td>0.023</td><td>0.097</td><td>0.080</td></tr><tr><td>N-BEATS-I (ours)</td><td>0.027</td><td>0.103</td><td>0.079</td></tr><tr><td>N-BEATS-I+G (ours)</td><td>0.025</td><td>0.099</td><td>0.077</td></tr></table>

Settings of hyperparameters across subsets of M4, M3, TOURISM datasets.

Table content: <table><tr><td/><td/><td/><td>M4</td><td/><td/><td/><td/><td/><td>M3</td><td/><td colspan="2">TOURISM</td><td/></tr><tr><td/><td colspan="13">Yly Qly Mly Wly Dly Hly Yly Qly Mly Other Yly Qly Mly</td></tr><tr><td>Parameter</td><td/><td/><td/><td/><td/><td colspan="2">N-BEATS-I</td><td/><td/><td/><td/><td/><td/></tr><tr><td>L H</td><td>1.5</td><td>1.5</td><td>1.5</td><td>10</td><td>10</td><td>10</td><td>20</td><td>5</td><td>5</td><td>20</td><td>20</td><td>10</td><td>20</td></tr><tr><td>Iterations</td><td colspan="3">15K 15K 15K</td><td colspan="3">5K 5K 5K</td><td colspan="2">50 6K</td><td>6K</td><td>250</td><td colspan="3">30 500 300</td></tr><tr><td>Losses</td><td/><td colspan="3">sMAPE/MAPE/MASE</td><td/><td/><td colspan="4">sMAPE/MAPE/MASE</td><td/><td>MAPE</td><td/></tr><tr><td>S-width</td><td/><td/><td/><td/><td/><td>2048</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>S-blocks</td><td/><td/><td/><td/><td/><td>3</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>S-block-layers</td><td/><td/><td/><td/><td/><td>4</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>T-width</td><td/><td/><td/><td/><td/><td>256</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>T-degree</td><td/><td/><td/><td/><td/><td>2</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>T-blocks</td><td/><td/><td/><td/><td/><td>3</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>T-block-layers</td><td/><td/><td/><td/><td/><td>4</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Sharing</td><td/><td/><td/><td/><td/><td colspan="2">STACK LEVEL</td><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">Lookback period</td><td/><td/><td colspan="5">2H, 3H, 4H, 5H, 6H, 7H</td><td/><td/><td/><td/><td/></tr><tr><td>Batch</td><td/><td/><td/><td/><td/><td>1024</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Parameter</td><td/><td/><td/><td/><td/><td colspan="2">N-BEATS-G</td><td/><td/><td/><td/><td/><td/></tr><tr><td>L H</td><td>1.5</td><td>1.5</td><td>1.5</td><td>10</td><td>10</td><td>10</td><td>20</td><td>20</td><td>20</td><td>10</td><td>5</td><td>10</td><td>20</td></tr><tr><td>Iterations</td><td colspan="3">15K 15K 15K</td><td colspan="3">5K 5K 5K</td><td colspan="3">20 250 10K</td><td>250</td><td colspan="3">30 100 100</td></tr><tr><td>Losses</td><td/><td colspan="3">sMAPE/MAPE/MASE</td><td/><td/><td colspan="4">sMAPE/MAPE/MASE</td><td/><td>MAPE</td><td/></tr><tr><td>Width</td><td/><td/><td/><td/><td/><td>512</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Blocks</td><td/><td/><td/><td/><td/><td>1</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Block-layers</td><td/><td/><td/><td/><td/><td>4</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Stacks</td><td/><td/><td/><td/><td/><td>30</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td>Sharing</td><td/><td/><td/><td/><td/><td>NO</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="2">Lookback period</td><td/><td/><td colspan="5">2H, 3H, 4H, 5H, 6H, 7H</td><td/><td/><td/><td/><td/></tr><tr><td>Batch</td><td/><td/><td/><td/><td/><td>1024</td><td/><td/><td/><td/><td/><td/><td/></tr><tr><td colspan="5">D HYPER-PARAMETER SETTINGS</td><td/><td/><td/><td/><td/><td/><td/><td/><td/></tr></table>

presents the hyperparameter settings used to train models on different subsets of M4, M3 and TOURISM datasets. A brief discussion of field names in the table is warranted.

Table content: <table/>

Peer review (iclr2020):
This goal of this paper is to present a strong empirical result showing that a *pure* machine learning based method can outperform all known methods on some of the most challenging time series forecasting benchmarks (TOURISM, M3 and especially M4). Since I am not from the field of forecasting, I can not be sure of this, but from my understanding these benchmark datasets are indeed challenging and the cited references back up the claims of the paper related to these datasets being important in the field. On the most challenging dataset (M4), the best known performing method combines RNNs with a traditional smoothing algorithm. The model proposed in this paper outperforms it without being combined with any classical approach, though it does utilize ensembling. The experimental setup is sound in my opinion, and the result appears to be of high potential significance. However, despite trying to go through Section 3 multiple times, the exact model architecture is not clear to me. Due to this reason, my current decision is a weak rejection since the model is a central contribution of the paper. I will be happy to increase my score if the authors can make the model description crystal clear. Even though my expertise is deep neural network architectures, I find it hard to follow the descriptions in Sec. 3. I faced the most difficulty understanding section 3.1, which obviously made the rest of subsections even harder to follow. Here are my main points of confusion: - One big issue is that the paper uses an illustration (Fig. 1) to explain the architecture instead of equations, but then uses symbols in the main text that do not appear on Fig. 1 at all such as g_theta. Is the *FC Stack (4 layers)* g_theta? - Where are (uppercase) phi functions? I could infer that these are the *FC* blocks but they should be labeled. - The Figure has the symbols g^b_theta and g^f_theta that do not appear in the text description. What exactly do they do? And is the theta that parameterizes each of them the same theta that parameterizes g_theta? How is this possible if g_theta is the *FC Stack (4 layers)*? - The description in second and third paragraph of Section 3.1 is very confusing and unclear. It should be replaced or augmented with equations using clearly defined symbols that match Figure 1. - More confusion stems from the use of the term *parameters* in (I believe) a different context than is used in neural networks, where *parameters* refers to connection weights. But here parameters are outputs of some functions, so either they are not connection weights or this is a fast-weight style architecture where outputs are weights [1], in which this should be made clear. - Design of the doubly residual architecture in Section 3.2 makes sense to me at a high level, but I feel it is still very hard to clearly understand and implement it. Again, use of equations to clearly define the computation would be very helpful. [1] Schmidhuber, Jürgen. *Learning to control fast-weight memories: An alternative to dynamic recurrent networks.* Neural Computation 4.1 (1992): 131-139. --- Update after rebuttal --- I am happy to see the paper greatly improved by the authors in their updates. My concerns related to the presentation of the model have been addressed, and I find the architecture much easier to understand. I also appreciate the detailed supplementary material that is likely to help readers interested in the area. Related to the areas I work in, I noticed the following missing references: Densenets: Lang, K.J. and Witbrock, M.J., 1988, June. Learning to tell two spirals apart. In Proceedings of the 1988 connectionist models summer school (No. 1989, pp. 52-59). Metalearning: Schmidhuber, J., 1987. Evolutionary principles in self-referential learning, or on learning how to learn: the meta-meta-... hook (Doctoral dissertation, Technische Universität München). While improving papers is generally the objective of the rebuttal phase, I suggest that the authors to not take this as an opportunity to submit unpolished papers in the first phase. That said, I have increased my rating to reflect my satisfaction with the current version of the paper.