Traditional forecasting models like ARIMA and ETS treat each time series separately, creating individual models. These are local models. They may not suit large-scale forecasting because maintaining individual models is impractical. Also, using exogenous covariates in traditional models often requires careful crafting and tuning.
Neural networks (NNs) can address these issues. They learn a compressed representation of input data, allowing a single model for many time series. This enables learning complex interactions between various data sources. These are global models. They simplify forecasting by removing the need for separate models and handcrafted features. Neural networks have a long history dating back to the 1960s but have become popular for forecasting only recently in specialised labs. Examples include Amazon DeepAR (Salinas et al., 2019), Google TFT (Lim et al., 2021), and Nixtla NHITS (Challu et al., 2023).
This chapter starts with a brief introduction to forecasting with neural networks in Section 14.1. Then, it introduces a core neural network architecture, the Multilayer Perceptron (MLP), in Section 14.2. It then covers more complex neural network architectures in Section 14.3. The later sections (14.4-14.7) explain how to forecast with neural networks using the NeuralForecast library.
14.1 Neural forecasting
When using neural networks to forecast, we aim to learn a (compressed) representation of our input data (Bengio et al., 2013). Mathematically, a neural network is a function f_\theta : X \mapsto Y, with X the input/feature space and Y the dependent variable space. We consider the setting with X = {y_{[0:t]}, x_{[0:t+h]}} and Y = {y_{[t+1:t+h]}}, where h is the forecast horizon, y is the target time series, and x are exogenous covariates. The forecasting task is to estimate the following conditional distribution:
P(y_{[t+1:t+h]} | y_{[0:t]}, x_{[0:t+h]}) = f_\theta (y_{[0:t]}, x_{[0:t+h]}) \tag{14.1}.
We can learn the function f_\theta by minimising an optimisation objective, commonly referred to as loss function. We can formalise the optimisation objective:
L = f(Y, \hat{Y}_\theta) \tag{14.2} ,
in which \hat{Y}_\theta denotes the forecasts produced by our neural network f_\theta. A frequently used optimisation objective when forecasting with neural networks is the Mean Absolute Error (MAE) (we will detail a few common optimisation objectives in Section 14.5):
\text{MAE}=\frac{1}{N} \sum_{i=0}^N |Y_i - \hat{Y}_{i, \theta}| \tag{14.3} .
We obtain the best forecast when the optimisation objective is minimised, which in the case of the MAE means the forecast \hat{Y}_{i, \theta} equals the actual observed value Y_i. The most used procedure to minimise the optimisation objective is called Gradient Descent (GD). In GD, we iteratively update the learnable parameters of the neural network by using the derivative of the optimisation objective with respect to the learnable parameters, i.e.
\theta_{i+1} \leftarrow \theta_i - \alpha \frac{\delta L}{\delta \hat{Y}_\theta} \tag{14.4},
where \alpha denotes the learning rate, a hyperparameter that allows us to adjust the rate of how fast the network learns the representation of the input data. We will discuss the optimisation of these hyperparameters in Section 14.7. The key idea is that if we perform a sufficient amount of iterations of GD, we keep updating our parameters \theta of the neural network in a direction that minimises the optimisation objective.
To summarise,
When forecasting with neural networks we aim to learn a (compressed) representation of our input data;
We can learn this representation by minimising an optimisation objective, or loss function;
The optimisation objective is commonly iteratively minimised using Gradient Descent (GD);
The optimisation procedure involves hyperparameters, which are tuning knobs that control the learning procedure.
14.2 Multilayer perceptron
A neural network can be thought of as a network of “neurons” which are organised in layers. The predictors (or inputs) form the bottom layer, and the forecasts (or outputs) form the top layer. There may also be intermediate layers containing “hidden neurons”.
The simplest networks contain no hidden layers and are equivalent to linear regressions. Figure 14.1 shows the neural network version of a linear regression with four predictors. We refer to the coefficients of these predictors as “weights” or “parameters”. The forecasts are obtained by a linear combination of the inputs. The weights are selected in the neural network framework using Gradient Descent that minimises a loss function such as the Mean Absolute Error (MAE) or the Mean Squared Error (MSE). Of course, in this simple example, we can use linear regression which minimises MSE and is a much more efficient method of training the model.
Figure 14.1: A simple neural network equivalent to a linear regression.
Once we add an intermediate layer with hidden neurons, the neural network becomes non-linear. A simple example is shown in Figure 14.2.
Figure 14.2: A neural network with four inputs and one hidden layer with three hidden neurons.
This is known as a Multilayer Perceptron (MLP), where each layer of nodes receives inputs from the previous layers. The outputs of the nodes in one layer are inputs to the next layer. The inputs to each node are combined using a weighted linear combination. The result is then modified by a nonlinear function before being output. For example, the inputs into each hidden neuron in Figure 14.2 are combined linearly to give
z_j = b_j + \sum_{i=1}^4 w_{i,j} x_i.
In the hidden layer, this is then modified using a nonlinear function such as a sigmoid,
s(z) = \frac{1}{1+e^{-z}},
to give the input for the next layer. This tends to reduce the effect of extreme input values, thus making the network somewhat robust to outliers.
The parameters b_1,b_2,b_3 and w_{1,1},\dots,w_{4,3} are “learned” (or estimated) from the data. The weights are initialised with random values, and are subsequently iteratively updated using Gradient Descent as described in Section 14.1. Consequently, there is an element of randomness in the predictions produced by a neural network.
The number of hidden layers and the number of nodes in each hidden layer are hyperparameters that must be specified in advance when trying to make predictions.
Example: Forecasting air passengers with an MLP
The AirPassengers dataset contains data on the number of US airline passengers from 1949 to 1960 for two different airlines. We split the data into a training and test set, where we use the last 12 months for testing. We use an MLP that takes in the last 24 months of observations (input_size =24), and the task is to forecast the next 12 months (horizon h=12).
# We use the last 12 months (1960) for testingtest_mask = AirPassengersPanel["ds"] >="1960"Y_train_df = AirPassengersPanel[~test_mask]Y_test_df = AirPassengersPanel[test_mask].reset_index(drop=True)# Instantiate the MLPmodel = MLP(h=12, input_size=24, scaler_type="robust")# Instantiate the forecasternf = NeuralForecast(models=[model], freq="M")# Fit the model to the training setnf.fit(df=Y_train_df)# Generate predictionsforecasts = nf.predict()
We can verify in Figure 14.3 that the MLP is able to predict the test set for both Airlines fairly accurately.
# Plot the predictions and the training datafig, axes = plt.subplots(2, sharex=True, figsize=(8, 7))plot_series(AirPassengersPanel, forecasts, xlabel="", ylabel="", palette="black_and_blue", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.3: AirPassengers MLP example.
14.3 Modern neural network architectures
The MLP set out previously is one of the simplest neural network architectures. Nowadays, there is a vast amount of research into the best neural network architectures for forecasting. NeuralForecast provides access to many of these modern architectures, listed in Table 14.1.
Table 14.1: Table of neural forecasting methods available in NeuralForecast.
Model
Family
Univariate / Multivariate
Forecast Type
Exogenous
0
Autoformer
Transformer
Univariate
Direct
F
1
BiTCN
CNN
Univariate
Direct
F/H/S
2
DeepAR
RNN
Univariate
Recursive
F/S
3
DeepNPTS
MLP
Univariate
Direct
F/H/S
4
DilatedRNN
RNN
Univariate
Recursive
F/H/S
5
FEDformer
Transformer
Univariate
Direct
F
6
GRU
RNN
Univariate
Recursive
F/H/S
7
HINT
Any
Both
Both
F/H/S
8
Informer
Transformer
Multivariate
Direct
F
9
iTransformer
Transformer
Multivariate
Direct
-
10
KAN
KAN
Univariate
Direct
-
11
LSTM
RNN
Univariate
Recursive
F/H/S
12
MLP
MLP
Univariate
Direct
F/H/S
13
MLPMultivariate
MLP
Multivariate
Direct
F/H/S
14
NBEATS
MLP
Univariate
Direct
-
15
NBEATSx
MLP
Univariate
Direct
F/H/S
16
NHITS
MLP
Univariate
Direct
F/H/S
17
NLinear
MLP
Univariate
Direct
-
18
PatchTST
Transformer
Univariate
Direct
-
19
RNN
RNN
Univariate
Recursive
F/H/S
20
SOFTS
MLP
Multivariate
Direct
-
21
StemGNN
GNN
Multivariate
Direct
-
22
TCN
CNN
Univariate
Recursive
F/H/S
23
TFT
Transformer
Univariate
Direct
F/H/S
24
TiDE
MLP
Univariate
Direct
F/H/S
25
TimeMixer
MLP
Multivariate
Direct
-
26
TimeLLM
LLM
Univariate
Direct
-
27
TimesNet
CNN
Univariate
Direct
F
28
TSMixer
MLP
Multivariate
Direct
-
29
TSMixerx
MLP
Multivariate
Direct
F/H/S
30
VanillaTransformer
Transformer
Univariate
Direct
F
In Table 14.1, we distinguish the following columns:
Model: The model name.
Family: The main neural network architecture underpinning the model. As you can see, there are many contemporary models that employ a variation of the MLP described above. Other common architecture families are Convolutional Neural Network (CNN), Transformers, and Recurrent Neural Network (RNN). It is beyond the scope of this book to explain all these building blocks in detail, and the interested reader is referred to the documentation of NeuralForecast for further explanation of each architecture.
Univariate / Multivariate: A multivariate model explicitly models the interactions between multiple time series in a dataset and will provide predictions for multiple time series concurrently. In contrast, a univariate model trained on multiple time series implicitly models interactions between multiple time series and provides predictions for single time series concurrently. Multivariate models are typically computationally expensive and empirically do not necessarily offer better forecasting performance compared to using a univariate model.
Forecast Type: Direct forecast models are models that produce all steps in the forecast horizon at once. In contrast, recursive forecast models predict one-step ahead, and subsequently use the prediction to compute the next step in the forecast horizon, and so forth. Direct forecast models typically suffer less from bias and variance propagation as compared to recursive forecast models, whereas recursive models can be computationally less expensive.
Exogenous: Whether the model accepts exogenous variables. This can be exogenous variables that contain information about the past and future (F), about the past only (historical, H), or that contain static information (static, S).
We turn to the AirPassengers dataset once again. Now, we train it using two more modern architectures, for example NHITS() and RNN(). We find that the results are comparable to the results we obtained using the relatively simple MLP architecture, which is due to the relative simplicity of the AirPassengers dataset (i.e., it doesn’t require a very sophisticated model to get near perfect forecasting performance on this dataset).
fig, axes = plt.subplots(2, sharex=True, figsize=(8, 5.5))plot_series(AirPassengersPanel, forecasts, xlabel="", ylabel="", palette="black_and_2color", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.4: AirPassengers NHITS and RNN example.
14.4 Scaling the data
In the AirPassengers example, we included scaler_type='robust'. This setting controls how data gets scaled before entering the neural network. Scaling becomes necessary because input ranges can vary widely, which may cause performance issues during training with gradient descent. Common scaling options include:
identity: no scaling
standard: scaled by mean and standard deviation
robust: a variation of standard scaling that handles outliers better. It removes the median and scales the data according to interquartile range (IQR). The IQR is the range between the 1st quartile (25th quantile) and the 3rd quartile (75th quantile).
Each technique scales time series and exogenous covariates based on the t values observed in the input window.
Let’s look at our MLP example when using identity scaling (i.e., no scaling at all). We can visually verify that the forecast has worsened, and it seems there is a bias present in the forecast.
model = MLP(h=12, input_size=24, scaler_type="identity")nf = NeuralForecast(models=[model], freq="M")nf.fit(df=Y_train_df)forecasts = nf.predict()
fig, axes = plt.subplots(2, sharex=True, figsize=(8, 5.5))plot_series(AirPassengersPanel, forecasts, xlabel="", ylabel="", palette="black_and_blue", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.5: AirPassengers MLP example using identity scaling.
14.5 Optimisation objectives
When forecasting with neural networks, we aim to learn a function f_\theta : X \mapsto Y, using a ground truth Y, such that we can generate forecasts \hat{Y}, i.e. f(X) = \hat{Y}. Ideally, there is no difference between the ground truth and the forecast: a perfect forecast without forecast error. Of course, with real data, this is impossible.
We measure the difference between the ground truth and the forecast using a loss function, which is our optimisation objective. For example, we can define the Mean Absolute Error:
\text{MAE}=\frac{1}{N} \sum_{i=0}^N |Y_i - \hat{Y}_i|
A perfect forecast has a zero MAE, and a very bad forecast a high MAE. Consequently, minimising the MAE by adjusting the weights of our neural network will improve the forecasting performance.
There are many different optimisation objectives. By default, NeuralForecast uses MAE as loss function. However, we can also use the Mean Squared Error (MSE):
\text{MSE}=\frac{1}{N} \sum_{i=0}^N (Y_i - \hat{Y}_i)^2
The MSE is more susceptible to outliers as a result of the quadratic term and should commonly be used with scaled data. Let’s see how our AirPassengers example looks like when using MSE as loss function.
Note that minimising MSE leads to forecasts equal to the mean of the forecast distribution, while minimising MAE leads to forecasts equal to the median of the forecast distribution.
fig, axes = plt.subplots(2, sharex=True, figsize=(8, 7))plot_series(AirPassengersPanel, forecasts, xlabel="", ylabel="", palette="black_and_blue", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.6: AirPassengers NHITS trained with MSE example.
Probabilistic optimisation objectives
In the previous example, we showed how to generate point predictions. By changing the optimisation objective, however, we can also generate probabilistic forecasts, in order to model the full conditional distribution of our forecasting task (rather than only the median (using MAE as loss function) or the mean (using MSE as loss function)).
For example, we can assume that our forecast distribution follows a Gaussian (‘Normal’) distribution by specifying a DistributionLoss. As can be seen, this brings the benefit that we now have a notion of uncertainty around the forecast, indicated by the 80- and 90- percentile shaded regions in the plot.
fig, axes = plt.subplots(2, sharex=True, figsize=(8, 7))# DistributionLoss accepts a level arg; the default is level=[80, 90]plot_series(AirPassengersPanel, forecasts, models=['NHITS'], level=[80, 90], xlabel="", ylabel="", palette="black_and_blue", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.7: AirPassengers NHITS trained with Gaussian DistributionLoss.
The prediction intervals shown in Figure 14.7 are generated by NHITS, which directly predicts future values and their variances in a single pass since we are using DistributionLoss(distribution="Normal"). However, this approach can lead to underestimated prediction intervals when the model overfits the training data. Since the model learns distribution parameters based on the training data, narrow intervals may result from the model capturing patterns too closely without properly accounting for future uncertainty.
Unlike some traditional statistical models, where uncertainty increases over time through recursive forecasting, neural models like NHITS rely on what they learn about the data’s distribution during training. When the training data has relatively low variability, this can further contribute to narrow prediction intervals.
14.6 Exogenous variables
A key benefit of neural networks is that we can add all sorts of different inputs to the forecasting model. For example, rather than only the values of the time series themselves, we can also provide the model with other sources of information that may be beneficial for our forecast. These exogenous covariates can be typically categorised into three distinct types:
Static covariates: covariates that are static throughout the forecast horizon. For example, the identifiers Airline1 and Airline2 in our AirPassengers dataset.
Historic covariates: covariates for which only the historic data is available. For example, the outside temperature can be a historic covariate for a task of forecasting ice cream demand: we don’t know the future outside temperature.
Future covariates: covariates for which we know their future values. For example, the day of the week can be used as a future covariate when trying to predict retail demand: we know the day of the week in advance when making forecasts for a particular day.
In NeuralForecast, each of these three types can be separately specified.
# Read the datadf = pd.read_csv("data/EPF_FR_BE.csv", parse_dates=["ds"])static_df = pd.read_csv("data/EPF_FR_BE_static.csv")futr_df = pd.read_csv("data/EPF_FR_BE_futr.csv", parse_dates=["ds"])# Specify the modelhorizon =24# day-ahead daily forecastmodel = BiTCN( h=horizon, input_size=5* horizon, futr_exog_list=["gen_forecast", "week_day"], # Future exogenous variables hist_exog_list=["system_load"], # Historical exogenous variables stat_exog_list=["market_0", "market_1"], # Static exogenous variables scaler_type="robust", max_steps=300, # Limit max_steps to control compute time)# Instantiate the forecasternf = NeuralForecast(models=[model], freq="H")# Fit the model to the training setnf.fit(df=df, static_df=static_df)# Generate predictionsforecasts = nf.predict(futr_df=futr_df)
Figure 14.8: BiTCN trained with exogenous variables.
14.7 Hyperparameter optimisation
Forecasting with neural networks often requires a practitioner to choose certain settings for training the forecasting model. These settings are commonly referred to as hyperparameters. We’ve seen a few of them already:
the number of hidden layers and number of nodes in an MLP;
the optimisation objective;
the input size (i.e., the amount of historic timesteps from a time series that the neural network considers for creating the forecast);
the scaler;
and so forth.
Ideally, a practitioner doesn’t have to choose these settings herself, but rather have an optimisation algorithm pick them. To this end, we can perform hyperparameter optimisation. NeuralForecast can perform the hyperparameter optimisation automatically based on a default configuration. To do so, it makes use of Auto models, which are versions of the models from Table 14.1 that allow one to automatically optimise hyperparameters of the neural network. For a specified number of iterations, an Auto model will intelligently select a set of hyperparameters from the eligible range and evaluate these on a validation set. At the end of the training process, the best performing hyperparameters are returned and can be used for the downstream forecasting application.
test_mask = AirPassengersPanel["ds"] >="1960"Y_train_df = AirPassengersPanel[~test_mask]Y_test_df = AirPassengersPanel[test_mask].reset_index(drop=True)# Get the default hyperparameter range for NHITS.# Update random_seed and n_pool_kernel_size with custom values.# Also limit max_steps to control compute time.nhits_config = {**AutoNHITS.get_default_config(h=12, backend="ray"),"random_seed": tune.randint(1, 10),"n_pool_kernel_size": tune.choice([[2, 2, 2], [16, 8, 1]]),"max_steps": tune.choice([100]),}# Instantiate the AutoNHITS model with the chosen configuration.# We optimise for 12 steps# i.e. we evaluate 10 different hyperparameter configurations.model = AutoNHITS(h=12, num_samples=10, config=nhits_config)# Instantiate the forecasternf = NeuralForecast(models=[model], freq="M")# Fit the model to the training setnf.fit(df=Y_train_df)# Generate predictionsforecasts = nf.predict()
fig, axes = plt.subplots(2, sharex=True, figsize=(8, 7))plot_series(AirPassengersPanel, forecasts, xlabel="", ylabel="", palette="black_and_blue", rm_legend=False, ax=axes, legend_loc="outside lower center")axes[-1].set(xlabel="Month [1M]")fig.suptitle("Number of passengers for different airlines", x=.54, fontsize=12)fig.supylabel("Passengers")fig
Figure 14.9: AirPassengers NHITS with automatic hyperparameter optimisation.
14.8 Exercises
Use the LSTM model to forecast 24-steps ahead for all the series in the pedestrian data set.
Consider the weekly data on US finished motor gasoline products supplied (millions of barrels per day) (series us_gasoline).
Use the MLP model to forecast the next 13 weeks of data.
Repeat the previous exercise using the AutoMLP and the AutoNHITS models from NeuralForecast.
Compare the accuracy of the models you generated.
Compare the accuracy of the models against the models of Chapter 12, Exercise 2.
14.9 Further reading
Bengio et al. (2013) provides an overview of representation learning.
Mariet & Kuznetsov (2019) for a more formal introduction to sequence-to-sequence modelling for time series.
Hewamalage et al. (2021) for a review of recurrent neural networks.
14.10 Used modules and classes
NeuralForecast
NeuralForecastclass - Core forecasting engine for neural network models
AutoNHITS()model - Automatic hyperparameter optimisation for NHITS
UtilsForecast
plot_series()function - For visualising time series data and forecasts
DatasetsForecast
AirPassengersPaneldataset - Classic airline passengers dataset for examples
Bengio, Y, Courville, A, & Vincent, P. (2013). Representation learning: A review and new perspectives. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8), 1798–1828.
Challu, C, Olivares, KG, Oreshkin, BN, Ramirez, FG, Mergenthaler Canseco, M, & Dubrawski, A. (2023). Nhits: Neural hierarchical interpolation for time series forecasting. Proceedings of the AAAI Conference on Artificial Intelligence, 37, 6989–6997.
Hewamalage, H, Bergmeir, C, & Bandara, K. (2021). Recurrent neural networks for time series forecasting: Current status and future directions. International Journal of Forecasting, 37(1), 388–427.
Lim, B, Arık, SÖ, Loeff, N, & Pfister, T. (2021). Temporal fusion transformers for interpretable multi-horizon time series forecasting. International Journal of Forecasting, 37(4), 1748–1764.
Mariet, Z, & Kuznetsov, V. (2019). Foundations of Sequence-to-Sequence Modeling for Time Series. The 22nd International Conference on Artificial Intelligence and Statistics, 408–417.
Salinas, D, Flunkert, V, Gasthaus, J, & Januschowski, T. (2019). DeepAR: Probabilistic forecasting with autoregressive recurrent networks. International Journal of Forecasting.
