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

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.

text: 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.

heading: 1 INTRODUCTION
text: 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?

heading: 1.1 SUMMARY OF CONTRIBUTIONS
text: 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).

heading: 2 PROBLEM STATEMENT
text: We consider the univariate point forecasting problem in discrete time. Given a length-H forecast horizon a length-T observed series history [y1, . . . ,yT ] ∈ RT , the task is to predict the vector of future values y ∈ RH = [yT+1,yT+2, . . . ,yT+H ]. For simplicity, we will later consider a lookback window of length t ≤ T ending with the last observed value yT to serve as model input, and denoted x ∈ Rt = [yT−t+1, . . . ,yT ]. 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):
sMAPE = 200 H
H
∑ i=1 |yT+i− ŷT+i| |yT+i|+ |ŷT+i| , MAPE = 100 H
H
∑ i=1 |yT+i− ŷT+i| |yT+i| ,
MASE = 1 H
H
∑ i=1
|yT+i− ŷT+i| 1
T+H−m ∑ T+H j=m+1 |y j− y j−m|
, OWA = 1 2
[ sMAPE
sMAPENaïve2 +
MASE
MASENaïve2
] .
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.

heading: 3 N-BEATS
text: 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.

heading: 3.1 BASIC BLOCK
text: 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
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 gb` 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`,1x`+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`,iv f i , x̂` =
dim(θ b` )
∑ i=1 θ b`,iv b i .
Here v fi and v b i are forecast and backcast basis vectors, θ f `,i is the i-th element of θ f ` . The function
of gb` and g f ` is to provide sufficiently rich sets {v f i } dim(θ f` ) i=1 and {vbi } dim(θ b` ) i=1 such that their respective outputs can be represented adequately via varying expansion coefficients θ f` and θ b ` . As shown below, gb` 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 gb` and g f ` are discussed in Section 3.3.

heading: 3.2 DOUBLY RESIDUAL STACKING
text: 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, x1 ≡ 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 gb` and g f ` for each layer, this makes the network more transparent to gradient flows. In a special situation of deliberate structure enforced in gb` and g f ` shared over a stack, explained next, this has the critical importance of enabling interpretability via the aggregation of meaningful partial forecasts.

heading: 3.3 INTERPRETABILITY
text: We propose two configurations of the architecture, based on the selection of gb` 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 gb` 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 gbs,` and g fs,` to be a polynomial of small degree p, a function slowly varying across forecast window:
ŷs,` = p
∑ i=0 θ fs,`,it 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̂trs,` = Tθ f s,`,
where θ fs,` are polynomial coefficients predicted by a FC network of layer ` of stack s described by equations (1); and T = [1, t, . . . , tp] is the matrix of powers of t. If p is low, e.g. 2 or 3, it forces ŷtrs,` 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 gbs,` and g f s,` to belong to the class of periodic functions, i.e. yt = yt−∆, where ∆ is a seasonality period. A natural choice for the basis to model periodic function is the Fourier series:
ŷs,` = bH/2−1c
∑ i=0 θ fs,`,i cos(2πit)+θ f s,`,i+bH/2c sin(2πit), (3)
The seasonality forecast will then have the matrix form as follows:
ŷseass,` = Sθ f s,`,
where θ fs,` 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πbH/2−1ct)),sin(2πt), . . . ,sin(2πbH/2−1ct))] is the matrix of sinusoidal waveforms. The forecast ŷseass,` 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 gbs,` and g f s,`. The number of blocks is 3 for both trend and seasonality. We found that on top of sharing gbs,` and g f s,`, sharing all the weights across blocks in a stack resulted in better validation performance.

heading: 3.4 ENSEMBLING
text: 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.

heading: 4 RELATED WORK
text: 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.

heading: 5 EXPERIMENTAL RESULTS
text: 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 (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.
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.

heading: 5.1 DATASETS
text: M4 (M4 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.
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.

heading: 5.2 TRAINING METHODOLOGY
text: 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 LH immediately preceding the last point in the train part of the TS. LH 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 LH is equal to 1.5; and in moderate to small Weekly, Daily, Hourly subsets of M4 LH 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.

heading: 5.3 INTERPRETABILITY RESULTS
text: 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.

heading: 6 DISCUSSION: CONNECTIONS TO META-LEARNING
text: 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.

heading: 7 CONCLUSIONS
text: 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.

heading: A DATASET DETAILS

heading: A.1 M4 DATASET DETAILS
text: 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).

heading: A.2 M3 DATASET DETAILS
text: 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).

heading: A.3 TOURISM DATASET DETAILS
text: 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.

heading: B ABLATION STUDIES

heading: B.1 LAYER STACKING AND BASIS SYNERGY
text: 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.

heading: B.2 ENSEMBLE SIZE
text: 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.

heading: B.3 DOUBLY RESIDUAL STACKING
text: 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–10. N-BEATSDRESS 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.

heading: C DETAILED EMPIRICAL RESULTS

heading: C.1 DETAILED RESULTS: M4 DATASET
text: 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.

heading: C.2 DETAILED RESULTS: M3 DATASET
text: Results for M3 dataset are provided in Table 14. The performance metric is calculated using the earlier version of sMAPE, defined specifically for the M3 competition:1
sMAPE = 200 H
H
∑ i=1 |yT+i− ŷT+i| yT+i + ŷT+i . (4)
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:
sMAPEAverage = NYear NTot sMAPEYear+ NQuart NTot sMAPEQuart+ NMonth NTot sMAPEMonth+ NOthers NTot sMAPEOthers .
(5)
Here NTot = NYear +NQuart +NMonth +NOthers and NYear = 6×645,NQuart = 8×756,NMonth = 18× 1428,NOthers = 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.
1With minor differences compared to the sMAPE definition used for M4. Please refer to Appendix A in (Makridakis & Hibon, 2000) for the mathematical definition.

heading: C.3 DETAILED RESULTS: TOURISM DATASET
text: Detailed results for the TOURISM competition dataset are provided in Table 15. The respective Kaggle competition was divided into two parts: (i) Yearly time series forecasting and (ii) Quarterly/Monthly time series forecasting (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:
MAPEAverage = NYear NTot MAPEYear+ NQuart NTot MAPEQuart+ NMonth NTot MAPEMonth . (6)
Here NTot = NYear +NQuart +NMonth and NYear = 4×518,NQuart = 8×427,NMonth = 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.

heading: C.4 DETAILED RESULTS: ELECTRICITY AND TRAFFIC DATASETS
text: In this experiment we are comparing the performances of MatFact (Yu et al., 2016), DeepAR (Flunkert et al., 2017) (Amazon Labs), Deep State (Rangapuram et al., 2018a) (Amazon Labs), Deep Factors (Wang et al., 2019) (Amazon Labs), and N-BEATS models on ELECTRICITY2 (Dua & Graff, 2017) and TRAFFIC3 (Dua & Graff, 2017) datasets. The results are presented in in Table 16.
Both datasets are aggregated to hourly data, but using different aggregation operations: sum for ELECTRICITY and mean for TRAFFIC. The hourly aggregation is done so that all the points available in (h− 1 : 00,h : 00] hours are aggregated to hour h, thus if original dataset starts on 2011-01-01 00:15 then the first time point after aggregation will be 2011-01-01 01:00. For the ELECTRICITY dataset we removed the first year from training set, to match the training set used in (Yu et al., 2016), based on the aggregated dataset downloaded from, presumable authors’, github repository4. We also made sure that data points for both ELECTRICITY and TRAFFIC datasets after aggregation match those used in (Yu et al., 2016). The authors of MatFact model were using the last 7 days of datasets as test set, but papers from Amazon are using different splits, where the split points are provided by a date. Changing split points without a well grounded reason adds uncertainties to the comparability of the models performances and creates challenges to the reproducibility of the results, thus we were trying to match all different splits in our experiments. It was especially challenging on TRAFFIC dataset, where we had to use some heuristics to find records dates; the dataset authors state: “ The measurements cover the period from Jan. 1st 2008 to Mar. 30th 2009” and “ We remove public holidays from the dataset, as well as two days with anomalies (March 8th 2009 and March 9th 2008) where all sensors were muted between 2:00 and 3:00 AM. ” , but we failed to match a part of the provided labels of week days to actual dates. Therefore, we had to assume that the actual list of gaps, which include holidays and anomalous days, is the following:
1. Jan. 1, 2008 (New Year’s Day)
2. Jan. 21, 2008 (Martin Luther King Jr. Day)
3. Feb. 18, 2008 (Washington’s Birthday)
4. Mar. 9, 2008 (Anomaly day)
5. May 26, 2008 (Memorial Day)
6. Jul. 4, 2008 (Independence Day)
7. Sep. 1, 2008 (Labor Day)
8. Oct. 13, 2008 (Columbus Day)
9. Nov. 11, 2008 (Veterans Day)
10. Nov. 27, 2008 (Thanksgiving)
11. Dec. 25, 2008 (Christmas Day)
12. Jan. 1, 2009 (New Year’s Day)
13. Jan. 19, 2009 (Martin Luther King Jr. Day)
14. Feb. 16, 2009 (Washington’s Birthday)
15. Mar. 8, 2009 (Anomaly day)
The first 6 gaps were confirmed by the gaps in labels, but the rest were more than 1 day apart from any public holiday of years 2008 and 2009 in San Francisco, California and US. More over the number of gaps we found in the labels provided by dataset authors is 10, while the number of days between Jan. 1st 2008 and Mar. 30th 2009 is 455, assuming that Jan. 1st 2008 was skipped from the values and labels we should end up with either 454−10 = 444 instead of 440 days or different end date. The metric is reported in Normalized deviation (ND) as in (Yu et al., 2016) which is equal to p50 loss used in DeepAR, Deep State, and Deep Factors papers.
2https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014 3https://archive.ics.uci.edu/ml/datasets/PEMS-SF 4https://github.com/rofuyu/exp-trmf-nips16/blob/master/python/exp-scripts/datasets/download-data.sh
ND = ∑i,t |Ŷit −Yit |
∑i,t |Yit | (7)
Contrary to Amazon models N-BEATS does not use any covariates, like day-of-week, hour-of-day, etc.
The N-BEATS architecture used in this experiment is exactly the same as used in M4, M3 and TOURISM datasets, the only difference is history size and the number of iterations. These parameters were chosen based on performance on validation set. Where the validation set consists of 7 consecutive days right before the test set. After the parameters are chosen the model is retrained on training set which includes the validation set, then tested on test set. The model is trained once and tested on test set using rolling window operation described in (Yu et al., 2016).

heading: C.5 DETAILED RESULTS: COMPARE TO DEEPAR, DEEP STATE SPACE MODELS
text: Table 17 compares ND (7) performance of DeepAR, DeepState models published in (Rangapuram et al., 2018a) and N-BEATS.

heading: D HYPER-PARAMETER SETTINGS
text: Table 18 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.
Subset names Yly, Qly, Mly, Wly, Dly, Hly, Other correspond to yearly, quarterly, monthly, weekly, daily, hourly and other frequency subsets defined in the original datasets.
N-BEATS-I and N-BEATS-G correspond to the interpretable and generic model configurations defined in Section 3.3.

heading: D.1 COMMON PARAMETERS
text: LH 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 LH 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 LH 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 LH is equal to 1.5; and in moderate to small Weekly, Daily, Hourly subsets of M4 LH 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.

heading: D.2 N-BEATS-I PARAMETERS
text: 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).

heading: D.3 N-BEATS-G PARAMETERS
text: 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).

heading: E DETAILED SIGNAL TRACES OF INTERPRETABLE INPUTS PRESENTED IN
text: FIGURE 2
The goal of this section is to show the detailed traces (numeric values) of signals visualized in Fig. 2. This is to demonstrate that even though the StackT-I (Fig. 2 (d)) and StackS-I (Fig. 2 (e)) provide response lines different from the counterparts in Stack1-G (Fig. 2 (b)) and Stack2-G (Fig. 2 (c)), the summations in the combined line (Fig. 2 (a)) can still be very similar.
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– 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.

Review summaries:
The paper proposes a DL architecture that achieves better performance on time series prediction. The proposed architecture is relatively straightforward and composes residual blocks. While the paper does achieve superior results, a lot of the text is devoted to comparing to prior work and arguing that DL approaches can do better than hand-crafted approaches, instead of focussing on the importance of specific technical contributions made in the paper. 

The paper investigates a pure deep learning architecture for Univariate time series analysis by simply ensembling feed-forward networks, along with the residual stacking mechanism for fluid learning. Each of the generic block consists of 4 FC layers followed by the use of forward and backward predictor to have forecast and backcast output of the original input. These blocks forms the stacks, where each stack provides the residuals and the forecast responses further to the next stacks, which ultimately provide the global forecast. To make the internal stack outputs interpretable, assumptions are imposed on the trend model, which follows a polynomial function of time vector, and seasonality model which follows periodic Fourier series. Further, the ensembling of models based on different metrics and input windows is used for better accuracy.