An Experimental Review on Deep Learning Architectures for Time Series Forecasting

In plain words

Deep learning is often used to predict the next values of a series of numbers, such as energy use, road traffic or website visits. There are many network designs for this, and each one has many settings to tune, so picking one is hard. This paper takes seven common designs, from a basic feed-forward network to recurrent and convolutional ones, trains a very large number of versions of each on twelve public collections of series, and records how accurate and how fast each one is.

The main lesson is practical. The recurrent LSTM was the most accurate when tuned well. The convolutional CNN came close and was much less sensitive to the settings, and it was also the quickest of the deep models. The basic feed-forward network was the least accurate and the fastest to train. The echo state network trained quickly but predicted slowly.

Read it as a map for choosing a starting point, not as a verdict for every problem. The authors conclude that the best choice depends on how accurate the forecasts must be, how fast they must come, and how much tuning you can afford.

TL;DR

A univariate forecasting benchmark of seven deep architectures on 12 public datasets with more than 50,000 series, each grid-searched and scored with WAPE §3.1 §3.2 tab. 8. Author-reported: LSTM has the best best-case accuracy on 4 of 12 datasets, CNN and LSTM rank best across configurations, and CNN gives the best accuracy-to-time balance §4.1 §4.3 tab. 11. The experimental grids and the moving-window protocol are reusable as a baseline. The headline "best" numbers are selected on the evaluation split, so treat them as optimistic (see Pitfalls).

Problem

The task is to forecast the next H values of one series from its last past values §2. The paper handles univariate series only §2. Existing reviews cover one architecture family or one data scenario, and new deep models are usually compared against classical methods rather than against other deep models §1. Classical ETS and ARIMA models need stationarity and a separate model per series, which makes them hard to scale to many series §2. Choosing a deep architecture and its many hyperparameters is described as a task that needs expertise and intuition §3.2.1.

Method

The paper has two parts: a literature review of MLP, recurrent and convolutional forecasters §2, tabs. 1–7, and an experimental comparison of seven architectures §3–§4. The implementable part is the experimental protocol below. The review is summarized in §2 and is not reproduced here.

Data and windows

Each dataset is a collection of series with a shared forecast horizon H. The study uses a global model per dataset that learns from all of its series §2, §3.1.

Each series is normalized with statistics from its own training part, using min-max scaling or z-score §3.3 eq. 2 eq. 3. The test values are scaled with the same training statistics preprocessing.py.

Training instances come from a moving window. The input is the past L = factor × H observations, the factor being 1.25, 2 or 3. The output is the next H observations. The window slides one step at a time. This is the multi-input multi-output (MIMO) strategy §3.3 fig. 2.

# Normalize each series with statistics of its own training part  [§3.3]
def normalize(series, method):
    if method == "minmax":
        # Scale to [0, 1] with the training minimum and maximum  [eq. 2]
        low = series.min()
        high = series.max()
        return (series - low) / (high - low)
    # Z-score with the standard deviation. The paper's eq. 3 divides by
    # max - min, but the code divides by the standard deviation  [eq. 3]
    # [code: experiments/preprocessing.py#L32-L36]
    return (series - series.mean()) / series.std()

# Turn one training series into MIMO (window, horizon) pairs  [§3.3] [fig. 2]
def make_windows(series, past_history, horizon):
    window_inputs = []
    window_targets = []
    # Slide one step at a time; keep only windows with a full horizon  [§3.3]
    # [code: experiments/preprocessing.py#L113-L119]
    for end in range(past_history, len(series) - horizon + 1):
        window_inputs.append(series[end - past_history : end])
        window_targets.append(series[end : end + horizon])
    return window_inputs, window_targets

For testing, each series contributes one forecast origin. The window is its last L observations, and the target is the held-out horizon of the same series §3.3 preprocessing.py.

Architectures

Every model maps a window of shape (L, 1) to H outputs through one dense layer with one neuron per horizon step §3.2.1. The paper names seven families and grids over their main parameters §3.2.1 tab. 9.

# Recurrent forecaster; cell is ERNN, LSTM, GRU or ESN  [§3.2.1]
def recurrent_forecaster(
    cell,
    layers,
    units,
    return_sequences,
    horizon,
    inputs,
):
    x = inputs  # (batch, past_history, 1)
    for index in range(layers):
        is_last = index == layers - 1
        # Only the last layer follows the return-sequence flag; earlier
        # layers return their full sequence  [code: experiments/models.py#L40-L49]
        returns_sequence = return_sequences if is_last else True
        x = cell(units, return_sequences=returns_sequence)(x)
    if return_sequences:
        # Flatten the (steps, units) output before the dense layer
        # [code: experiments/models.py#L54-L55]
        x = flatten(x)
    # One output neuron per horizon step  [§3.2.1]
    return dense(horizon)(x)

# 1D CNN forecaster; each block is (filters, kernel, pool)  [§3.2.1]
def cnn_forecaster(blocks, horizon, inputs):
    x = inputs  # (batch, past_history, 1)
    for filters, kernel, pool in blocks:
        # Same-padded 1D convolution with ReLU
        # [code: experiments/models.py#L200-L202]
        x = conv1d(filters, kernel, padding="same", activation="relu")(x)
        # Max-pool by `pool` only if pool is set and the length is still
        # above it; pool 0 means no pooling
        # [code: experiments/models.py#L204-L214]
        if pool and length(x) // pool > 1:
            x = max_pool1d(pool)(x)
    x = flatten(x)
    # One output neuron per horizon step  [§3.2.1]
    return dense(horizon)(x)

# Temporal convolutional forecaster  [§2.3.1] [§3.2.1]
# [code: experiments/models.py#L228-L266]
def tcn_forecaster(
    inputs,
    filters,
    kernel,
    stacks,
    dilations,
    return_sequences,
    horizon,
):
    # Dilated causal convolutions with residual blocks; `stacks` is the
    # "number of TCN layers" of the paper  [§3.2.1]
    x = tcn_block(inputs, filters, kernel, stacks, dilations)
    if return_sequences:
        x = flatten(x)
    # One output neuron per horizon step  [§3.2.1]
    return dense(horizon)(x)

Training

Each configuration is trained with Adam and a mean absolute error loss on the normalized targets §3.3 main.py. The paper does not name the training loss §3.3. The code runs for five epochs tab. 10 and caps each epoch at 10,000 steps main.py.

# Train one configuration on the MIMO windows  [§3.3] [tab. 10]
def train(model, train_inputs, train_targets, batch_size, learning_rate):
    # Adam with the learning rate from the grid  [§3.3]
    # [code: experiments/main.py#L171]
    optimizer = Adam(learning_rate=learning_rate)
    # Mean absolute error on normalized targets
    # [code: experiments/main.py#L177] [code: experiments/models.py#L279]
    loss = "mae"
    # One pass over the windows, capped at 10,000 steps per epoch
    # [code: experiments/main.py#L167-L169]
    steps = min(ceil(len(train_inputs) / batch_size), 10000)
    # Five epochs with shuffled batches  [tab. 10]
    # [code: experiments/main.py#L184-L192]
    model.fit(
        train_inputs,
        train_targets,
        epochs=5,
        steps_per_epoch=steps,
        shuffle=True,
    )

Evaluation and statistics

Accuracy is measured with WAPE on the denormalized held-out horizon §3.3.1 eq. 4. Efficiency is measured as training time and inference time per model §3.3.1.

# Weighted absolute percentage error for one series  [eq. 4]
def wape(actual, forecast):
    # Mean absolute error over the horizon  [eq. 4]
    absolute_error = mean(abs(actual - forecast))
    # Divide by the mean actual value, with a 1e-10 guard  [eq. 4]
    # [code: experiments/metrics.py#L58-L60]
    return absolute_error / (mean(actual) + 1e-10)

# Score one trained model on the test split  [§3.3]
# [code: experiments/main.py#L202-L215]
def evaluate(model, test_windows, test_targets, norm_params):
    # One forecast origin per series  [§3.3]
    forecasts = model.predict(test_windows)
    scores = []
    for index in range(len(test_windows)):
        # Undo the per-series normalization  [code: experiments/main.py#L207-L210]
        forecast = denormalize(forecasts[index], norm_params[index])
        actual = denormalize(test_targets[index], norm_params[index])
        # Skip series whose actual values are all (nearly) identical
        # [code: experiments/metrics.py#L163-L171]
        if not is_constant(actual):
            scores.append(wape(actual, forecast))
    # Average the per-series WAPE  [code: experiments/metrics.py#L185]
    return mean(scores)

The paper ranks architectures with a Friedman test and a Holm–Bonferroni post-hoc procedure at α = 0.05. It builds eleven rankings: accuracy at each model's best, its average accuracy, the spread of accuracy, its worst accuracy, the same three views for training and inference time, and an overall average of those rankings §3.3.2 §4.3 fig. 6. For configuration choices it uses a paired Wilcoxon signed-rank test at α = 0.05 §3.3.2 §4.5 tab. 16 tab. 17.

Hyperparameters

Grids as searched by the paper tab. 9 tab. 10. Code values come from the released experiments/parameters.json and the model files.

Best settings per architecture, as selected by the Wilcoxon tests. ** means significantly best at p ≤ 0.05, * means a tendency at p ≤ 0.2, and = means no difference tab. 17.

Best training values per architecture tab. 18. ** and * as above.

Results

Author-reported. Setup: seven architectures, 12 univariate datasets with more than 50,000 series in total §3.1 tab. 8, grids of tab. 9 and tab. 10, WAPE on the held-out horizon eq. 4. The experiments took about four months, on five NVIDIA TITAN Xp GPUs, an Intel i7-8700 CPU and an additional cloud GPU §4 intro. The only baselines are the other six architectures; no classical statistical models are run §3.2.

Best WAPE per dataset across each architecture's configurations tab. 11. The bold cells are the best per dataset. Dataset numbers follow tab. 8: 1 CIF2016o12, 2 CIF2016o6, 3 ExchangeRate, 4 M3, 5 M4, 6 NN5, 7 SolarEnergy, 8 Tourism, 9 Traffic, 10 Traffic-metr-la, 11 Traffic-perms-bay, 12 WikiWebTraffic.

Headline findings, all author-reported:

• Best-cell accuracy: LSTM is best on 4 of 12 datasets. ERNN, GRU and TCN are best on 2 each, and ESN and CNN on 1 each §4.1 tab. 11. The paper reports that the best-cell differences are not statistically significant, except for MLP, which is worst §4.3 fig. 6.
• Mean WAPE across configurations: CNN ranks first and LSTM second, and CNN and LSTM also have the smallest spread §4.1 §4.3 fig. 6 fig. 3.
• Speed: MLP trains fastest. CNN is the most efficient deep model in the rankings, and ESN has slow inference §4.2 fig. 4 fig. 5.
• Trade-off: LSTM gives the best accuracy with adequate time, and CNN gives the best time-accuracy balance §4.2 fig. 5. The paper concludes that CNN and LSTM are the most suitable choices overall §4.3.
• Within recurrent models, fewer stacked layers did better, except for ESN (4 layers, full sequences) §4.4.2 tab. 17 fig. 7.
• CNN: four layers without pooling was best; single-layer models ranked last §4.4.3 tab. 14.
• TCN: a single TCN layer with kernel 6 and 64 filters, returning the last output only, ranked first §4.4.3 tab. 15.

Paper vs code

Code is the authors' repository at commit 2fc5468bc7717e42feb686f85fac2fd7ec84c457 (MIT) main.py. I read the code. I did not reproduce the numbers.

• Z-score: the paper's eq. 3 divides by max − min eq. 3. The code divides by the standard deviation preprocessing.py.
• CNN depth: the paper grids over 1, 2 and 4 layers tab. 9 and describes the 4-layer best model §4.4.3. The shipped grid stops at 3 layers, and its 3-layer block uses kernels 7, 5, 2 parameters.json. The paper gives no 3-layer kernels.
• Training grid: the paper grids over two learning rates, three past-history factors, two normalizations and two batch sizes tab. 10. The shipped parameters.json uses one learning rate, one past-history factor, one normalization and batch sizes 32 and 64 parameters.json.
• Training loss: the paper does not name the loss §3.3. The code uses MAE main.py models.py.
• Step cap: the paper says nothing about a cap. The code caps steps per epoch at 10,000 main.py parameters.json.
• Validation: the paper has no validation split. The code passes the test windows as validation_data main.py.
• Activations: the paper does not state them. The code uses no activation in the MLP's dense layers models.py, ReLU in the CNN models.py models.py, and a linear activation in the TCN models.py.
• ESN head: the paper describes one dense output layer §3.2.1. The code adds a Dense(32) layer before the output for ESN only models.py.
• Pooling: the paper says a pooling factor of 2 is considered §3.2.1, while the grid includes 0 tab. 9. The code skips pooling when the sequence is too short models.py.
• WAPE aggregation: the paper defines WAPE per series eq. 4. The code averages the per-series values and skips constant series metrics.py.
• Test split: the paper says the test set is the last horizon of each series §3.3. The code loads a separate test.csv per dataset from URLs in data/datasets.json generate_data.py.
• Not checked: results/results.csv and results/appendix.pdf ship in the repository. I did not reconcile them with tab. 11.

Pitfalls

• Selection on the evaluation split. The paper reports, for each configuration, only the best WAPE over the training grid §4.4. The code passes the test windows as validation data main.py. The best-cell numbers in tab. 11 are therefore best-of-grid on the split that is reported. Compare architectures by mean and spread §4.1 fig. 3 rather than by single best cells.
• Grid count. The paper says 6,432 configurations per dataset §1. Its experimental setup section gives the same 6,432 figure across all 12 datasets and reports more than 38,000 models §3.2. These statements disagree. The grids in tab. 9 and tab. 10 give 134 architectures × 24 training settings = 3,216 per dataset (my count). 3,216 × 12 = 38,592, which matches "more than 38,000". The paper does not explain the factor of two.
• Repetitions and seeds. The paper says each configuration was "trained several times" §4.4. The code runs each configuration once main.py. No random seed is set in experiments/ (grep), so exact numbers are not reproducible.
• Linear MLP as implemented. The code's MLP stacks dense layers with no activation models.py. As written, that model is linear in its input. The paper's remarks about MLP error are not therefore evidence about non-linear MLPs.
• One origin per series. Each series contributes one test window preprocessing.py. The ranking rests on one forecast per series. Small collections such as CIF2016 and ExchangeRate have few series tab. 8, so their numbers are noisy.
• Epoch cap. Five epochs with the 10,000-step cap tab. 10 main.py can mean each epoch sees only part of the windows on large datasets. This follows from reading the code. The paper does not describe it.
• Zero and near-zero means. WAPE divides by the mean of the actual values plus 1e-10 metrics.py. Only constant series are skipped metrics.py. Non-constant series with a near-zero mean can produce very large per-series values that dominate a mean.
• Normalization guards. The min-max guard for a zero range does take effect preprocessing.py. The z-score guard does not: it assigns 1e-10 to a local variable and then divides by the unguarded value, so a constant series divides by zero preprocessing.py. Test values are not clipped, so they can fall outside the training range preprocessing.py.
• Text contradicts table. The paper states that returning the full sequence before the dense layer hurts TCN §4.5 §4.4.3. Its TCN selection table lists return-sequence as "True*" tab. 17. The best TCN in tab. 15 has return-sequence False.
• Statistical analysis not in the repo. The paper's Friedman, Holm–Bonferroni and Wilcoxon analyses §3.3.2 have no code in the repository (grep finds no friedman, wilcoxon or holm in the Python files).
• Timing is confounded. The paper names mixed hardware §4 intro. The code selects a GPU by index when more than one is present main.py. The paper does not say which device ran which model, so time comparisons between architectures may mix hardware.
• Shipped data script is partial. generate_data.py skips every dataset except SolarEnergy generate_data.py. Remove that line to regenerate all datasets.
• TCN parameter naming. The paper calls the TCN depth "number of layers" §3.2.1. The code calls it nb_stacks parameters.json. Map one to the other when you implement the grid.

Lineage

• Builds on: the paper's own review of MLP, Elman RNN, LSTM, echo state, GRU, CNN and TCN forecasters §2. The cited originals are not on the wiki.
• Compared against: only the other six architectures in the same grid §3.2. No classical statistical forecasters are run in the experiments.
• Superseded by: not known from content already on the wiki.