Neural forecasting: from lag MLPs to N-BEATS and N-HiTS

Statistical forecasting models such as exponential smoothing and ARIMA each write down a small, explicit story about how a series evolves: a level that adapts, a trend that persists, a seasonal pattern that repeats. Neural forecasting models take a different route. They turn the series into many input-output examples, then let a flexible function learn the mapping from recent history to the next stretch of the future. The story is no longer written by hand; it is learned from data, sometimes from thousands of related series at once.

This guide builds that idea from the ground up. We start by turning a series into supervised training windows, look at how multi-step outputs are produced, why every window is rescaled, and what a "global" model is. Then we climb through a plain multilayer perceptron on lags, recurrent encoders, probabilistic outputs in the style of DeepAR, and two architectures designed specifically for forecasting: N-BEATS and N-HiTS. The tutorial trains all of them with Nixtla's neuralforecast on the AirPassengers series and compares them honestly against a seasonal naive baseline. If forecast origins, backtests and prediction intervals are new to you, read the mathematical intuition guide first; the statistical models guide covers the baselines we compare against.

Forecasting as supervised learning

A neural network needs a fixed-size input and a fixed-size target. A time series is one long sequence, so we cut it into windows. Pick an input length $L$ (in neuralforecast this is input_size) and a horizon $h$. At every forecast origin $t$, the input is the last $L$ observed values and the target is the next $h$ values:

Take a tiny series of six monthly values, $112, 118, 132, 129, 121, 135$, with $L = 3$ and $h = 2$. Sliding the origin forward one step at a time gives two complete training pairs:

A real monthly series of 132 training months with $L = 36$ and $h = 12$ yields $132 - 36 - 12 + 1 = 85$ complete windows. That is not much data for a neural network, which is one reason these models shine when many series are trained together. Each window is one row of a supervised dataset; the network learns a function $f_\theta$ with parameters $\theta$ so that $f_\theta(\mathbf{x}_t)$ is close to $\mathbf{y}_t$ on average.

The information boundary from the intuition guide still applies. Every value in $\mathbf{x}_t$ must have been known at time $t$. Any statistic used to normalize the window, and any extra feature, must also respect that boundary. neuralforecast constructs these windows internally, which removes a common source of off-by-one leakage, but you still choose $L$, $h$ and the holdout.

flowchart LR
  S["Raw series y1 ... yT"] --> W["Slide a window of input_size L plus horizon h"]
  W --> N["Scale each window using only its L input values"]
  N --> F["Network f maps L inputs to h outputs"]
  F --> D["Undo the scaling on the h outputs"]
  D --> Lo["Loss against the true next h values"]
  Lo -.->|"gradient step"| F

Choosing $L$ is a modelling decision. It should cover at least one full seasonal cycle, so the network can see what last January looked like when it predicts next January. For monthly data with yearly seasonality, $L$ between 24 and 36 is a sensible start. Much longer windows give the network more context but fewer training examples from a short series.

Direct multi-output versus recursive forecasting

There are two main ways to produce $h$ future values from a model.

The direct multi-output strategy asks the network to emit all $h$ values at once. The output layer has $h$ units, and the training loss averages over every horizon step:

where $N$ is the number of training windows and $\ell$ is a pointwise loss such as absolute error. Each horizon step gets its own output unit, so the network can learn that twelve months ahead behaves differently from one month ahead. Predictions never feed back into the input, so errors do not compound.

The recursive (autoregressive) strategy trains a one-step model $g_\theta$ and feeds its own predictions back:

and so on until step $h$. With the tiny series above and $L = 3$, if the one-step model predicted $\hat{y} = 125$ for month 7, the input for month 8 would become $(121, 135, 125)$: a mix of real data and the model's own guess. If that guess is off by 10, the error enters the next input and can grow. The recursive strategy is natural for recurrent networks and for generating sample paths, and it can forecast any horizon without retraining, which is why DeepAR uses it.

In neuralforecast, MLP, N-BEATS and N-HiTS are direct multi-output models. The recurrent models (LSTM, GRU, RNN) default to recurrent=False, which uses the recurrent network as an encoder and a separate decoder that outputs all $h$ steps directly; setting recurrent=True switches to recursive prediction. DeepAR is recursive by design.

Scaling each window and training global models

Suppose one series sells around 100 units a month and another around 10,000. Feed raw values into a network trained on both and the large series dominates the gradient: an error of 50 on the large series swamps an error of 5 on the small one, even though the second is relatively ten times worse. Neural networks also train best when inputs are roughly centred and of unit scale, because weight initialization and learning-rate defaults assume it.

The fix is to normalize each input window by statistics computed from that window alone. With the robust scaler, which uses the median and the mean absolute deviation from the median,

the network predicts on the normalized scale, and the output is mapped back with $\hat{y} = \mathrm{med}(\mathbf{x}_t) + \mathrm{mad}(\mathbf{x}_t)\, \hat{\tilde{y}}$. Because the statistics come only from the $L$ input values, nothing about the target leaks into the scaling. In neuralforecast this is the per-model scaler_type argument. The installed version accepts identity (no scaling), standard (mean and standard deviation), robust (median and mean absolute deviation, less sensitive to outliers), minmax, minmax1, invariant and revin (reversible instance normalization with learnable affine parameters). The source code confirms that the target positions are masked out when these statistics are computed. A separate local_scaler_type on the NeuralForecast object instead scales each whole series once before windowing, which is useful but must be fitted on training data only.

For AirPassengers, which roughly quadruples in level over twelve years, per-window scaling means the network sees "the shape of the last three years relative to their own level". A pattern learned in the 1950s can then be reused at the higher level of 1960.

Scaling also enables global models: one set of weights $\theta$ trained on windows from many series $i = 1, \dots, M$,

A local model such as ETS estimates separate parameters for every series. A global model shares parameters, so patterns learned from series with long histories can help series with short ones, and the number of training windows grows with the number of series. This pooling is where neural forecasters usually earn their keep. With a single series, as in our tutorial, a neural model has to learn everything from fewer than a hundred windows, and a well-chosen statistical model is a strong competitor.

MLPs on lags and recurrent encoders

The simplest neural forecaster is a multilayer perceptron applied to the scaled lag vector. With one hidden layer,

where $W_1$ has shape hidden size by $L$ and $W_2$ has shape $h$ by hidden size. Think of each hidden unit as a learned detector: one may respond to "the series rose over the last three months", another to "this month resembles the same month last year". The output layer combines these detectors into $h$ forecasts. An MLP treats lag 1 and lag 12 as unrelated inputs; it has no built-in idea of order. That is a weakness for very long inputs but rarely a problem when $L$ is a few dozen.

A recurrent network reads the window in order and carries a hidden state $\mathbf{h}_s$ forward:

so the final state $\mathbf{h}_t$ is a learned summary of the whole input. A plain recurrent network struggles to carry information over many steps because gradients shrink or explode as they flow backward through repeated multiplications. The LSTM adds a memory cell $\mathbf{c}_s$ and gates that decide what to forget, write and expose:

where $\mathbf{f}_s$, $\mathbf{i}_s$ and $\mathbf{o}_s$ are sigmoid gates between 0 and 1, $\tilde{\mathbf{c}}_s$ is a candidate update, and $\odot$ is element-wise multiplication. The additive cell update gives gradients a path that is not repeatedly squashed. The backpropagation guide for CNNs and RNNs derives backpropagation through time and shows why the gates help. In neuralforecast's LSTM, the encoder produces hidden states and an MLP decoder maps them to the $h$ outputs.

Probabilistic outputs: likelihoods and quantiles

A point forecast cannot tell an inventory planner how much safety stock to hold. Neural models can output distributions in two common ways.

Parametric likelihoods, DeepAR style

DeepAR, introduced in the DeepAR paper, is an autoregressive recurrent network. At each step the LSTM state $\mathbf{h}_s$ is mapped to the parameters of a chosen distribution, and training maximizes the likelihood of the observed values. For a Gaussian output, two small linear heads produce a mean and a positive scale:

and the loss is the negative log-likelihood,

Read the two terms as a trade-off. The squared-error term punishes a mean far from the observation, but dividing by $\sigma_s^2$ lets the model claim uncertainty to soften the penalty. The $\log \sigma_s$ term charges a price for claiming uncertainty. Ignoring the constant, if $y_s = 110$ and $\mu_s = 100$, then $\sigma_s = 5$ gives $1.61 + 2.0 = 3.61$, while $\sigma_s = 10$ gives $2.30 + 0.5 = 2.80$, and $\sigma_s = 50$ gives $3.91 + 0.02 = 3.93$. The likelihood prefers a scale that matches the typical size of the errors.

Other distributions change only the head and the loss. A Student-t adds a degrees-of-freedom parameter $\nu$; small $\nu$ gives heavier tails, so occasional large shocks are less surprising. A negative binomial suits non-negative counts, such as units sold. Parameterized by mean $\mu$ and shape $\alpha$, its variance is

so uncertainty grows with the level, as it usually does for counts. At prediction time DeepAR draws a value from the predicted distribution, feeds it back as the next input, and repeats to step $h$. Many such sample paths give empirical quantiles. In neuralforecast you pick the family with DistributionLoss("Normal"), DistributionLoss("StudentT"), DistributionLoss("NegativeBinomial") and others, and the level argument chooses which intervals are reported.

Quantile (pinball) loss

The second approach skips distributional assumptions and trains the network to output chosen quantiles directly. For quantile level $\tau$ and residual $u = y - \hat{q}_\tau$, the pinball loss is

With $\tau = 0.9$ and $y = 100$: predicting $\hat{q} = 90$ under-forecasts by 10 and costs $0.9 \times 10 = 9$; predicting $\hat{q} = 110$ over-forecasts by 10 and costs $0.1 \times 10 = 1$. Under-forecasting is nine times more expensive, so the network learns to push the 0.9 quantile up until roughly 90% of observations fall below it. More formally, the expected loss $\mathbb{E}[\rho_\tau(y - q)]$ has derivative $F(q) - \tau$ with respect to $q$, where $F$ is the cumulative distribution of $y$; it is zero exactly when $F(q) = \tau$, which is the definition of the $\tau$ quantile. With $\tau = 0.5$ the loss is half the absolute error, so the median is the special case.

MQLoss in neuralforecast averages the pinball loss over several quantiles, and the network's output layer grows from $h$ to $h$ times the number of quantiles. MQLoss(level=[80, 90]) trains quantiles 0.05, 0.10, 0.50, 0.90 and 0.95. Because each quantile has its own output unit, nothing forces them to stay ordered; we will see a crossing in the tutorial.

N-BEATS: deep stacks of fully connected blocks

N-BEATS showed that a carefully structured stack of plain fully connected layers, with no recurrence or convolution, can be a strong univariate forecaster. The paper reports competitive results on the M3, M4 and tourism competitions using large ensembles of such models. Its two key ideas are the block and doubly residual stacking.

The block

Block $\ell$ receives an input vector $\mathbf{x}_\ell$ of length $L$. A stack of fully connected layers with ReLU activations (four in the paper) produces a hidden vector, which two linear heads turn into expansion coefficients:

Basis matrices then turn the coefficients into two outputs: a backcast $\hat{\mathbf{x}}_\ell = V^b_\ell \boldsymbol{\theta}^b_\ell$ of length $L$, which is the block's explanation of its own input, and a forecast $\hat{\mathbf{y}}_\ell = V^f_\ell \boldsymbol{\theta}^f_\ell$ of length $h$.

Doubly residual stacking

Blocks are chained so that each one sees only what the previous blocks failed to explain, and the forecasts are summed:

"Doubly residual" refers to these two paths: a residual connection on the input side, subtracting what has been explained, and additive aggregation on the output side. Imagine the input is a rising line plus a yearly wave. If the first block captures the rising line in its backcast, the second block receives a nearly flat wave and can focus on seasonality. The subtraction also gives gradients a short path to every block, which keeps deep stacks trainable. Blocks are grouped into stacks, and the final forecast is the sum over every block in every stack.

flowchart TD
  X["Input window x of length L"] --> B1
  subgraph S1["Stack 1: trend"]
    B1["Block 1<br/>4 FC layers<br/>theta b and theta f"]
  end
  subgraph S2["Stack 2: seasonality"]
    B2["Block 2<br/>4 FC layers<br/>theta b and theta f"]
  end
  B1 -->|"backcast 1"| M1["Subtract: x minus backcast 1"]
  X --> M1
  M1 --> B2
  B1 -->|"forecast 1"| Sum["Sum of block forecasts"]
  B2 -->|"forecast 2"| Sum
  Sum --> Y["Final forecast of length h"]

Generic and interpretable bases

In the generic configuration, $V^b_\ell$ and $V^f_\ell$ are ordinary learned matrices: maximum flexibility, no built-in meaning. The interpretable configuration fixes the basis so each stack can only produce a particular kind of shape.

Let $\mathbf{t} = (0, 1, \dots, h-1)/h$ be normalized time over the horizon. A trend block restricts its forecast to a low-degree polynomial,

with small $p$, often 2 or 3. With $p = 1$ the block can only output a straight line $\theta_0 + \theta_1 \mathbf{t}$, so it cannot accidentally absorb seasonal wiggles. A seasonality block uses a Fourier basis,

which can express any smooth periodic pattern over the horizon. Because the trend stack's output is literally a polynomial and the seasonality stack's output is literally a sum of sinusoids, you can plot each stack's contribution and read it as a decomposition. In neuralforecast, stack_types takes a list from identity (generic), trend and seasonality; n_polynomials sets the polynomial degree and n_harmonics controls how many Fourier terms are used. The tutorial uses one block of each type.

N-HiTS: multi-rate pooling and hierarchical interpolation

N-BEATS has a cost problem for long horizons. Every block's final layer outputs $h$ forecast values from a large hidden vector, and every block reads the full-resolution input. Forecasting 720 hourly steps from a long window means large output layers and a model prone to producing jagged, high-frequency forecasts. N-HiTS keeps the N-BEATS block-and-residual skeleton and adds two changes.

Multi-rate input pooling

Before its fully connected layers, block $\ell$ max-pools the input with kernel size $k_\ell$:

With $L = 36$ and $k_\ell = 4$, the block sees 9 values instead of 36, each the maximum of four consecutive months. Early blocks use large kernels and see a coarse, smoothed picture suited to slow movements. Later blocks use small kernels and see fine detail. Pooling shrinks the first layer and pushes each block toward a particular frequency band.

Hierarchical interpolation

Instead of predicting all $h$ values, block $\ell$ predicts only $\lceil h / r_\ell \rceil$ knots, where $r_\ell$ is its downsampling ratio (n_freq_downsample in the library). The full forecast comes from interpolating between knots. With linear interpolation, for a horizon position $\tau$ between neighbouring knots at times $t_1$ and $t_2$,

In the tutorial, $h = 12$ and the ratios are $6, 2, 1$: the first block outputs 2 knots, a coarse line across the year; the second outputs 6 knots, enough for a bimonthly shape; the last outputs all 12 for fine corrections. The residual structure still applies, so each block corrects what earlier blocks left behind, now organized from low to high frequency.

flowchart LR
  X["Input window L = 36"] --> P1["Block 1<br/>MaxPool kernel 4<br/>9 inputs"]
  P1 --> K1["MLP outputs 2 knots"]
  K1 --> I1["Interpolate to 12 steps"]
  X --> R1["Subtract backcast 1"]
  P1 -.-> R1
  R1 --> P2["Block 2<br/>MaxPool kernel 2<br/>18 inputs"]
  P2 --> K2["MLP outputs 6 knots"]
  K2 --> I2["Interpolate to 12 steps"]
  R1 --> R2["Subtract backcast 2"]
  P2 -.-> R2
  R2 --> P3["Block 3<br/>no pooling<br/>36 inputs"]
  P3 --> K3["MLP outputs 12 values"]
  I1 --> Sum["Sum of block forecasts"]
  I2 --> Sum
  K3 --> Sum

Why is this cheaper? Take $h = 720$ with a ratio of 24: the output head shrinks from 720 units to 30, and with a 512-unit hidden layer that is roughly 15,000 output weights instead of 369,000 for that block. Pooling shrinks the input layers in the same way. The paper reports better accuracy than N-BEATS on several long-horizon benchmarks at a fraction of its compute; those are the authors' reported numbers on their benchmarks. On a 12-step monthly horizon the savings are small, and you should not expect N-HiTS to dominate N-BEATS there.

Implement it: setup and data

The code below was run in a fresh virtual environment on a CPU-only Mac with Python 3.14, neuralforecast 3.3.0, statsforecast 2.1.1, PyTorch 2.15.0 and PyTorch Lightning 2.6.6. Newer or older versions may change defaults or outputs.

python3 -m venv .venv
source .venv/bin/activate
pip install neuralforecast statsforecast

AirPassengersDF stores timestamps as month ends (1949-01-31, 1949-02-28, and so on), so the frequency string must be "ME". Passing "MS" (month start) produces forecasts dated on the first of each month, which then silently fail to join with the test rows. The series also arrives as float32 with a numeric id; we convert both for cleaner output.

import logging, warnings
import numpy as np
import pandas as pd
from neuralforecast import NeuralForecast
from neuralforecast.models import MLP, LSTM, DeepAR, NBEATS, NHITS
from neuralforecast.losses.pytorch import MAE, MQLoss, DistributionLoss
from neuralforecast.utils import AirPassengersDF
from statsforecast import StatsForecast
from statsforecast.models import SeasonalNaive

warnings.filterwarnings("ignore")
logging.getLogger("lightning.pytorch").setLevel(logging.ERROR)

df = AirPassengersDF.copy()
df["unique_id"] = "AirPassengers"
df["y"] = df["y"].astype("float64")

h, m = 12, 12                    # horizon and season length
train = df.iloc[:-h]             # 1949-01 .. 1959-12 (132 months)
test = df.iloc[-h:]              # 1960-01 .. 1960-12 (12 months)

# Shared settings: 36 lags, robust per-window scaling, short CPU training.
quiet = dict(accelerator="cpu", enable_progress_bar=False,
             enable_model_summary=False, logger=False)
common = dict(h=h, input_size=36, max_steps=300, random_seed=1,
              scaler_type="robust", **quiet)

The quiet dictionary is passed through to the PyTorch Lightning trainer. random_seed fixes weight initialization and window sampling, but exact reproducibility across machines and library versions is not guaranteed.

Implement it: train five models and compare

All five models go into one NeuralForecast object. Each is trained separately; the object just handles windowing, batching and the output frame. Hidden sizes are reduced from the defaults to keep training to well under a minute on a laptop CPU.

models = [
    MLP(num_layers=2, hidden_size=256, **common),
    LSTM(encoder_hidden_size=64, decoder_hidden_size=64, **common),
    DeepAR(loss=DistributionLoss("StudentT", level=[80, 90]),
           valid_loss=MAE(), lstm_hidden_size=64, **common),
    NBEATS(stack_types=["identity", "trend", "seasonality"],
           mlp_units=[[256, 256]] * 3, **common),
    NHITS(n_pool_kernel_size=[4, 2, 1], n_freq_downsample=[6, 2, 1],
          mlp_units=[[256, 256]] * 3, **common),
]
nf = NeuralForecast(models=models, freq="ME")
nf.fit(df=train)
fc = nf.predict()                # one row per future month, one column per model

sf = StatsForecast(models=[SeasonalNaive(season_length=m)], freq="ME")
sn = sf.forecast(df=train, h=h)
fc = fc.merge(sn, on=["unique_id", "ds"]).merge(test, on=["unique_id", "ds"])

# MASE scale: in-sample MAE of the seasonal naive forecast on the training data
y_tr = train["y"].to_numpy()
scale = np.mean(np.abs(y_tr[m:] - y_tr[:-m]))

rows = []
for c in ["SeasonalNaive", "MLP", "LSTM", "DeepAR", "NBEATS", "NHITS"]:
    e = fc["y"] - fc[c]
    rows.append((c, e.abs().mean(), np.sqrt((e ** 2).mean()), e.abs().mean() / scale))
print(pd.DataFrame(rows, columns=["model", "MAE", "RMSE", "MASE"]).round(2))

Fitting all five models took about 45 seconds on CPU. The results on the 1960 holdout:

These numbers come from a single run with random_seed=1 and 300 training steps, using the versions listed above. Change the seed, the step count or the library version and individual entries can move by several units; treat the ordering among MLP, LSTM, N-BEATS and N-HiTS as unsettled. MASE below 1 means a model beat the in-sample seasonal naive error; the seasonal naive itself scores 1.57 on the holdout because the series grew strongly during 1960, and repeating 1959 under-forecasts every month.

Implement it: intervals, cross-validation and automatic search

Reading DeepAR's interval

DistributionLoss adds columns for the median and the requested intervals, such as DeepAR-lo-80 and DeepAR-hi-80. In this run, the 80% interval contained the actual value in 6 of 12 test months. In the summer peak, the July 1960 actual of 622 sat above the interval of 505.6 to 535.3. That is a strongly overconfident interval, together with a weak point forecast: with only 300 steps on 85 windows, the recursive likelihood model is clearly undertrained here. Do not read this as a verdict on DeepAR in general; it was designed for large collections of related series.

Quantile N-HiTS with MQLoss

qmodel = NHITS(h=h, input_size=36, loss=MQLoss(level=[80, 90]),
               n_pool_kernel_size=[4, 2, 1], n_freq_downsample=[6, 2, 1],
               mlp_units=[[256, 256]] * 3, scaler_type="robust",
               max_steps=300, random_seed=1, **quiet)
nfq = NeuralForecast(models=[qmodel], freq="ME")
nfq.fit(df=train)
q = nfq.predict().merge(test, on=["unique_id", "ds"])
print(q[["ds", "y", "NHITS-lo-80", "NHITS-median", "NHITS-hi-80"]].round(1))
for lv in [80, 90]:
    inside = (q.y >= q[f"NHITS-lo-{lv}"]) & (q.y <= q[f"NHITS-hi-{lv}"])
    print(lv, "coverage", inside.mean())

The median forecast from this model had MAE 12.29, better than any point model in the table. But both the 80% and 90% intervals covered only 6 of 12 months, with average widths of about 22 and 25 passengers (in thousands). The intervals are too narrow. Look also at January 1960: the "upper" 80% bound was 416.2 while the median was 420.9. The quantiles crossed, because nothing in MQLoss forces the output units to stay ordered. Sorting the quantile columns row by row is a common practical fix; a better fix is more data and longer training. Twelve correlated test months are also far too few to estimate coverage reliably, which is one reason to backtest over several origins.

Rolling-origin cross-validation

cv_model = NHITS(h=h, input_size=36, n_pool_kernel_size=[4, 2, 1],
                 n_freq_downsample=[6, 2, 1], mlp_units=[[256, 256]] * 3,
                 scaler_type="robust", max_steps=300, random_seed=1, **quiet)
cv = NeuralForecast(models=[cv_model], freq="ME").cross_validation(
    df=df, n_windows=3, step_size=12, refit=True)
cvs = StatsForecast(models=[SeasonalNaive(season_length=12)], freq="ME").cross_validation(
    df=df, h=12, n_windows=3, step_size=12)
cv = cv.merge(cvs.drop(columns=["y"]), on=["unique_id", "ds", "cutoff"])
print(cv.groupby("cutoff").apply(lambda g: pd.Series({
    "NHITS_MAE": (g.y - g.NHITS).abs().mean(),
    "SNaive_MAE": (g.y - g.SeasonalNaive).abs().mean()})).round(2))

This forecasts 1958, 1959 and 1960 from origins at the end of 1957, 1958 and 1959. With refit=True the network is retrained at each origin; the default refit=False trains once on data before the first origin and reuses those weights, which is cheaper but penalizes later folds. Our run gave:

The first fold is the important one. In 1958 air travel grew slowly, so "same as last year" was excellent, and N-HiTS, extrapolating the earlier growth, was almost three times worse. A single holdout year would have hidden this. Notice also that the 1959-12 fold uses the same training data as the holdout experiment yet reports 13.70 rather than 18.09. We did not track down which internal difference in the cross-validation path causes this, but it is a concrete reminder that a small neural model's error on one year can move by several units for incidental reasons.

Automatic hyperparameter search

Every architecture above has many knobs: input_size, learning rate, scaler, hidden sizes, pooling kernels. The Auto* classes (AutoMLP, AutoLSTM, AutoNBEATS, AutoNHITS and others) wrap a model in a search. Each trial samples a configuration, trains on the training data minus a validation tail of val_size steps, and is scored on that validation tail; the best configuration is refitted and used for prediction. The backend is Ray Tune by default or Optuna, and each class ships a default search space. We ran a deliberately tiny Optuna search, four trials of 200 steps each, so the result says more about the mechanics than about the best possible model:

from neuralforecast.auto import AutoNHITS

def config(trial):
    return {
        "input_size": trial.suggest_categorical("input_size", [24, 36]),
        "learning_rate": trial.suggest_float("learning_rate", 1e-4, 1e-2, log=True),
        "scaler_type": trial.suggest_categorical("scaler_type", ["robust", "standard"]),
        "max_steps": 200,
        "random_seed": 1,
        **quiet,
    }

auto = AutoNHITS(h=h, config=config, backend="optuna", num_samples=4, loss=MAE())
nfa = NeuralForecast(models=[auto], freq="ME")
nfa.fit(df=train, val_size=12)          # the last 12 training months score each trial
print(nfa.models[0].results.best_trial.params)
pa = nfa.predict().merge(test, on=["unique_id", "ds"])
print((pa.y - pa.AutoNHITS).abs().mean())

The best trial chose input_size=36, a learning rate of about 0.00023 and the standard scaler, with validation MAE 11.30 on 1959. Its 1960 test MAE was 14.44. With so few trials and a one-year validation window, that improvement could easily be luck. Real searches use dozens of trials, and the validation tail should be long enough to be representative.

How to read these results

Start with the baseline. Every neural point model except DeepAR beat seasonal naive on the 1960 holdout by a wide margin. That says these networks learned the growth that seasonal naive ignores. It does not say they beat a good statistical model: the statistical models guide fits ETS and ARIMA variants on exactly this split, and you should put the two tables side by side before drawing conclusions.

Look at the spread among the four neural point models. N-BEATS, N-HiTS and LSTM landed within about 1.3 MAE of each other. On one series and twelve test points, that spread is well within what a different seed or training length can produce. A credible ranking needs several origins, several seeds, and ideally many series.

Read the cross-validation table as a lesson about regimes. The neural model shines when the future continues the recent growth and suffers when growth pauses. Seasonal naive does the opposite. An ensemble or a model choice made per origin over many folds is more trustworthy than a winner chosen on one year.

Read intervals with suspicion. Both probabilistic models here produced intervals too narrow to reach their nominal coverage, and one showed quantile crossing. With small data and short training, neural uncertainty estimates are often overconfident. Check coverage across many backtest windows before handing an interval to a planner.

Pitfalls

Frequency mismatches. A wrong freq string shifts future timestamps. Here "MS" versus "ME" produced forecasts that matched none of the test dates, and the metrics silently became missing values. Always check the forecast dates and the row count after merging.

Too few windows. With 132 training months and $L = 36$ there are only 85 complete windows, and a validation tail removes more. Large networks will memorize them. Shrink the model, shorten the input, or pool many series.

No scaling on trending data. With scaler_type="identity", which is the default for several models, a series that quadruples in level presents the network with values it has never seen at forecast time. Use a temporal scaler, or a local_scaler_type fitted on training data only.

Judging on one seed and one holdout. Neural training is stochastic. Report the mean and spread over several seeds and several backtest origins.

Recursive models and long horizons. Recursive models such as DeepAR feed their own predictions back. Errors and overly narrow sampling distributions compound with the horizon.

Validation data leaking into model choice. If you tune hyperparameters by looking at the test year, the test score is no longer an honest estimate. Use val_size or cross-validation folds that end before the test period, and look at the test once.

Trusting the default refit=False in backtests. It is faster, but the model at the last origin has not seen the most recent years. Know which one you used before comparing against a baseline that is effectively refitted at every origin.

Quantile crossing and miscalibration. Sort quantiles if you must use them, measure coverage by horizon, and remember that pinball-trained quantiles are only as good as the amount of data behind them.

When to choose neural forecasters

flowchart TD
  A{"Many related series<br/>or long histories?"} -->|"no, one short series"| B["Start with statistical models<br/>ETS, ARIMA, Theta"]
  A -->|"yes"| C{"Need probabilistic forecasts?"}
  C -->|"yes, counts or heavy tails"| D["DeepAR style likelihood<br/>NegativeBinomial or StudentT"]
  C -->|"yes, specific quantiles"| E["N-HiTS or N-BEATS with MQLoss"]
  C -->|"point forecasts"| F{"Long horizon?"}
  F -->|"yes"| G["N-HiTS"]
  F -->|"no"| H["N-BEATS, MLP or LSTM<br/>tune with Auto classes"]
  B -.->|"baseline for every branch"| E

Neural forecasters are most worthwhile when you have many related series to train one global model, when histories are long, when the horizon is long enough that N-HiTS's efficiency matters, or when you need exogenous inputs and shared patterns that local models cannot express. They need more care than a statistical model: scaling, seeds, training length, and validation design all matter, and training costs compute.

For a single short series like AirPassengers, a well-specified ETS or ARIMA model is usually just as accurate, cheaper, and easier to explain. Use the neural models here as a learning exercise and as candidates, not as presumed winners. Choose N-BEATS with interpretable stacks when stakeholders want to see trend and seasonal components. Choose N-HiTS for long horizons or long inputs. Choose a likelihood model when the decision depends on a distribution and the data justify a parametric family. When attention and covariates matter, the transformer models guide is the next step; when you have little data of your own, the foundation models guide covers pretrained zero-shot forecasters.

Practice exercises

1. Rerun the five-model comparison with seeds 1 through 5. Report the mean and standard deviation of MAE per model. Does the ordering among N-BEATS, N-HiTS and LSTM survive?
2. Set scaler_type="identity" for the MLP and compare. Explain the change using the level shift between the 1950s and 1960.
3. Change input_size to 12, 24 and 48 for N-HiTS. Which values cannot see a full seasonal cycle, and how does that show up in the forecasts?
4. Replace N-BEATS's stacks with ["identity", "identity", "identity"] (generic). Compare accuracy, then fit the interpretable version again and plot each stack's contribution on the test year.
5. Train DeepAR with DistributionLoss("Normal") and with 1,000 steps. Measure 80% coverage over the three cross-validation folds rather than one year.
6. Sort the MQLoss quantile columns row by row and recompute coverage. Does sorting fix calibration, or only the crossing?
7. Run cross_validation with refit=False and refit=True and explain the differences fold by fold.

References

• Oreshkin et al., N-BEATS: Neural basis expansion analysis for interpretable time series forecasting
• Challu et al., N-HiTS: Neural Hierarchical Interpolation for Time Series Forecasting
• Salinas et al., DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks
• NeuralForecast documentation
• StatsForecast documentation
• Hyndman and Athanasopoulos, Forecasting: Principles and Practice, 3rd edition