NeuralProphet: Explainable Forecasting at Scale

We introduce NeuralProphet, a successor to Facebook Prophet, which set an industry standard for explainable, scalable, and user-friendly forecasting frameworks. With the proliferation of time series data, explainable forecasting remains a challenging task for business and operational decision making. Hybrid solutions are needed to bridge the gap between interpretable classical methods and scalable deep learning models. We view Prophet as a precursor to such a solution. However, Prophet lacks local context, which is essential for forecasting the near-term future and is challenging to extend due to its Stan backend. NeuralProphet is a hybrid forecasting framework based on PyTorch and trained with standard deep learning methods, making it easy for developers to extend the framework. Local context is introduced with auto-regression and covariate modules, which can be configured as classical linear regression or as Neural Networks. Otherwise, NeuralProphet retains the design philosophy of Prophet and provides the same basic model components. Our results demonstrate that NeuralProphet produces interpretable forecast components of equivalent or superior quality to Prophet on a set of generated time series. NeuralProphet outperforms Prophet on a diverse collection of real-world datasets. For short to medium-term forecasts, NeuralProphet improves forecast accuracy by 55 to 92 percent.


Background
Time series data is prominent in most industrial sectors. Though extensively studied in theory and applications such as economics, practical forecasting in industrial applications has not received widespread attention until recently.
Since many companies have significantly improved their data collection to the extent where data availability now exceeds their data analytic capabilities, how to process and forecast industrial time series is quickly becoming an important topic. Often the person tasked with solving a particular issue is not a time series expert but rather a practitioner forecaster.
Accurate time series forecasting is essential for decision-making processes, especially in infrastructure planning, budget allocation and supply chain management. Over-forecasts as well as under-forecasts can be both be costly. In applications, where forecasts inform business or operational decisions with potentially fatal consequences, it is often a requirement for the model to be explainable. A common approach is to break the forecast down into components which are individually interpretable. This allows for the forecast to be reviewed by a domain expert who may make adjustments when appropriate, a procedure known as human-in-the-loop.
Classical Time Series Methods. The field of forecasting was traditionally dominated by statistical techniques. This was most notable on many of the early forecasting competitions such as NN3, NN5 and M3 (Makridakis & Hibon, 2000;Crone, 2008). Classical models such as Auto-Regressive Integrated Moving Average (ARIMA) and Exponential Smoothing (ETS) have been well studied and provide interpretable components. However, their restrictive assumptions and parametric nature limit their performance in real world applications. A skillful forecasting expert can transform data and combine algorithms to satisfy specific conditions for better performance. This requires deep domain knowledge in the application itself and in classical time series modelling.
Machine Learning Methods. Machine Learning (ML) based models were constantly performing poorly at the early forecasting competitions. Neural Networks (NN) were once even labeled as not being competitive for forecasting (Hyndman, 2018). They were further criticized in literature for their black-box nature, as in Makridakis, Spiliotis & Assimakopoulos (2018). However, with the explosion in availability of large scale time series, NN-based data-driven techniques regained their popularity in forecasting. The amount of data available is no longer insufficient for ML and time series components introduced by the Prophet package with NN modules into a hybrid model enabling it to fit to non-linear dynamics. Two such NN components are the auto-regression and covariate modules. By using PyTorch as the backend, NeuralProphet can be updated with the latest innovation in deep learning. We provide a convenient tool for non-experts in forecasting by pre-selecting strong defaults and automating many of the modelling decisions. Advanced users may also find the package useful as they can incorporate domain knowledge with many customizaiton options. Overall, NeuralProphet abstracts a lot of forecasting domain knowledge and incorporates ML best practices so users can focus on the task at hand. in this paper. We introduce the NeuralProphet model and its components in detail and compare it to Prophet. The methods section describes the experimental setup and evaluation. The synthetic data results quantify the accuracy of the decomposed interpretable forecast components. The realworld benchmarks demonstrate the realistic out-of-the-box forecasting performance on a diverse set of univariate time series. Finally, we summarize the findings and provide suggestions on which model to use depending on the task at hand.

Model Components
A core concept of the NeuralProphet model is its modular composability. The model is composed of modules which each contribute an additive component to the forecast. Most components can also be configured to be scaled by the trend for a multiplicative effect. Each module has its individual inputs and modelling processes. However, all modules must produce h outputs, where h defines the number of steps to be forecasted into the future at once. These are added up as the predicted valuesŷ t , ...,ŷ t+h−1 for the time series future values y t , ..., y t+h−1 . If the model is only time-dependent, an arbitrary number of forecasts can be produced. In the following descriptions, that special case will be treated mathematically equivalent to a one-step ahead forecast with h = 1. In the following subsections we discuss each of the components in more detail.

Trend
A classic approach to modelling trend is to model it as the combination of an offset m and a growth rate k. The trend effect at a time t 1 is given by multiplying the growth rate by the difference in time since the starting point t 0 on top of the offset m.
NeuralProphet models the trend component with this classic approach, but allows the growth rate to change at a number of locations. Thus, the trend is modelled as a continuous piece-wise linear series. This results in an interpretable, yet non-linear form of trend modelling. It is simple to interpret, as in a segment between two points, the trend effect is given by the steady growth rate multiplied by the difference in time. We can generalize the trend by defining a time-dependent growth rate δ(t) and a time-dependent offset ρ(t).
The piece-wise linear trend only varies the growth rate at a finite number of changepoints. A set C of n c changepoints are defined at different times as C = (c 1 , c 2 , ..., c nc ). Between changepoints, the trend growth rate is kept constant. The first segment's growth rate and offset are given as δ 0 and ρ 0 respectively. Rate adjustments at each changepoint can be defined as a vector δ ∈ R n C , where δ j is the rate change at the j th change-point. The growth rate at time t is determined by adding the initial growth rate δ 0 with the summation of the rate adjustments at all the changepoints up to time step t. Each growth rate change δ j is a parameter to be fitted on the data. A corresponding vector of offset adjustments can be defined as ρ ∈ R nc . Similarly, the offset at time t is given by the initial offset ρ 0 and the sum of the offset adjustments at each change-point up to time t. However, the offset for a changepoint c j is not an independent parameter but defined as ρ j = −c j δ j . This particular offset definition makes the piece-wise linear series continuous. Further, we introduce a binary vector Γ(t) ∈ R nc representing whether the time t is past each changepoint.
Thus, we can define a vectorized equation for the trend T (t) at time t, which is given in equation 3: where, δ = (δ 1 , δ 2 , ..., δ n C ) ρ = (ρ 1 , ρ 2 , ..., ρ n C ) NeuralProphet provides a simple, semi-automatic mechanism for the selection of relevant changepoints. Given the number n C of desired change-points, n C equidistant points along the series are selected as initial changepoints. Optionally, their growth change rate parameters can be regularized during model training. This is is similar to a fully automatic changepoint selection, as only the most relevant changepoints will be selected, if any. The user can also opt to manually define the specific times of a custom number of changepoints. To avoid overfitting on a small number of final points, the final trend segment (after the last changepoint) is set to a larger set of observations (default: 15 % of training data). When making predictions into the unobserved future, the final growth rate is used to linearly extrapolate the trend.

Seasonality
Seasonality in NeuralProphet is handled by using Fourier terms (Harvey & Shephard, 1993) as was originally done in Prophet (Taylor & Letham, 2017 Fourier terms are a great tool for modelling seasonality as they produce smooth functions which are simple to interpret and stable to fit to data. However, Fourier terms only model deterministic seasonal shapes which are assumed to be fixed through time. A higher number of Fourier terms allows the model to fit a more complex seasonal pattern. Too much flexibility may lead to overfitting or to random patterns between observations. Thus, each Fourier term corresponds to a frequency proportional to j p , modelled by by a weighted combination of a sine and cosine transform. Every seasonality is associated with 2k number of coefficients. For time step t, the effect from all the seasonalities considered in the model can be indicated by S(t) in below Equation 5, where P refers to the set of all the periodicities.
Both additive and multiplicative seasonal patterns are supported. Each seasonal periodicity S p can individually be configured to be multiplicative, in which case the component is multiplied by The framework automatically activates daily, weekly and or yearly seasonality depending on data frequency and length. Each of these three types of seasonal periodicities is activated if the data frequency is higher resolution than the respective periodicity, and if at least two full periods of data are available. As an example, if the data is of daily frequency, the model will enable yearly seasonality if the data spans two years or more. Weekly seasonality will also be added if two or more weeks of data are available. Daily seasonality will not be activated, as the daily frequency is not of higher resolution to allow for intra-day seasonality. Unless otherwise specified, the default number of Fourier terms per seasonality are: k = 6 for p = 365.25 yearly, k = 3 for p = 7 weekly, and k = 6 for p = 1 daily seasonality.

Auto-Regression
Auto-regression (AR) refers to the process of regressing a variable's future value against its past values, a critical part of many forecasting applications as can be seen in Hyndman & Athanasopoulos (2014). The number of past values included is usually referred to as the order p of the AR(p) model. Hereby, a coefficient θ i is fitted for each past value. Each coefficient θ i controls the direction and magnitude of effect of a particular past value on the forecast. In a classic AR process, an intercept c and a white noise term t are included.
In many applications, we are interested to predict multiple steps into the future. We refer to the number of forecasts made at once as the forecast horizon h. A traditional AR model will only produce one prediction for a one-step-ahead forecast with h = 1. If a multi-step ahead forecast of length h > 1 is demanded, h models will have to be fitted, one for each forecast step. The AR module in NP is based on a modified version of AR-Net in Triebe, Laptev & Rajagopal (2019).
AR-Net is able to produce all h forecasts with one model, which can be of linear or non-linear nature. In any configuration, the p last observations of the target variable y t−1 , y t−2 , ..., y t−p , also referred to as lags, are the inputs to the module. The outputs are h values corresponding to the AR-effect for each forecast step A t (t), A t (t + 1), ..., A t (t + h − 1). Hereby, A t (t + 2) denotes the predicted AR-effect for forecastŷ t+2 at t + 2, three steps in the future, predicted at forecast origin t with observed data up to and including t − 1.
A t (t), A t (t + 1), ..., A t (t + h − 1) = AR-Net(y t−1 , y t−2 , ..., y t−p ) It is important to note, that each time we forecast at a specific origin, we obtain h predictions.
Thus, for a given point in time, we have up to h different predictions, each originating from a different forecast made in the past. They differ from each other based on the data available to the model at the time of the forecast. As an example, when we use an AR(5) model to forecast h = 3 steps into the future, at a given time t = 3, we will have 3 predicted values of the AR-effect for y t=3 . These three effects A t=1 (t = 3), A t=2 (t = 3), ..., A t=3 (t = 3) are each an estimation of the true AR-effect A(t = 3), with different 'age'. The oldest estimate A t=1 (t = 3) is 3 steps old, while the most recent estimate A t=3 (t = 3) is one step old. With Auto-Regression, any estimate will be at least one step old, as we assume to never know the true value of the series at the current time t but only the observed value at the last step t − 1 and older.
AR order. The most important parameter for this module is the number of past values to be regressed over, also known as the order p of the AR(p) model. This parameter should be chosen based on the approximate length of relevant context in past observations. In practice, it is hard to determine accurately and is commonly set to twice the innermost periodicity or twice the forecast horizon. Alternatively, a conservatively large order can be chosen when used in combination with regularization to obtain a sparse AR-model.
Linear AR. The default AR-Net configuration does not contain hidden layers, and is functionally identical to a classic AR model. In practice, it is a single layer NN with p inputs, h outputs, no biases, and no activation function. The single layer weights each regress a particular lag onto a particular forecast step. Thus, each weight can be matched to a corresponding coefficient of a collection of h classic AR(p) models, making it simple to interpret the model.
We can denote the model as a vector-matrix multiplication for the predicted AR effects y ∈ R h , with lagged observations as input x ∈ R p , and with weight matrix W ∈ R h×p , where W i,j denotes the coefficent defining the linear impact of lag y j on AR-effect A t (t = i) Deep AR. AR-Net based AR module can model non-linear dynamics when hidden layers are configured. In this case, the module trains a fully connected Neural Network (NN) with the specified number of hidden layers and dimensions. The addition of hidden layers may lead to a better forecasting accuracy, however it a partial trade-off in interpretability. Instead of being able to directly quantify the contribution of a particular past observation to a particular prediction, we can only observe the relative importance of a given past observation on all predictions. This can be approximated by comparing the sums of the absolute weights of the first layer for each input position.
The p last observations of the time series are the inputs to the first layer. After each hidden layer, the logits are passed through an activation function, in our case, a rectified linear unit (ReLU). The final layer outputs h logits, are not transformed by an activation function and have no bias. For l hidden layers with a hidden layer dimension of size d, we have: Hereby, the layer biases all have the same dimensions b ∈ R d , while the layer weights are W ∈ R d×d , except for the first W 1 ∈ R d×p and last W l+1 ∈ R h×d layer weights.
Sparse AR. AR-Net demonstrated that the correct order can be approximated by setting the order to a slightly larger than expected value when regularization is used to sparsify the model weights.
This allows for a more convenient model configuration while retaining interpretability.
The regularization function proposed by the authors of AR-Net is given in equation 8. Though we make this regularization function available, by default, our implementation uses a different regularization function, which we found to work better for a wider range of data. We use the regularization function introduced in equation 14 with parameters = 3 and α = 1, where θ is the vector of weights corresponding to the p AR model coefficients: where, As an example, when forecasting the next 24 steps of an hourly time series, the forecaster could choose to set the AR order to 100 with some regularization. After fitting, the AR module may then exhibit sparse coefficients with the most significant weights located on a few positions only. In this example, only positions 1, 2, 3, 24, 48, and 72 could have non-zero weights. These positions can be interpreted as the locations with an auto-correlation, and further, their weights can be examined to study how different combinations of past observation values affect the forecast.
For reference, the original regularization function proposed in AR-Net is given by: where the authors suggest setting c 1 ≈ 3 and c 2 ≈ 3.

Lagged Regressors
Lagged regressors are used to correlate other observed variables to our target time series. They are often referred to as covariates. Unlike future regressors, the future of lagged regressors is unknown to us. At the time t of forecasting, we only have access to their observed, past values up to and including t − 1.
Given a set of covariates X ∈ R T ×n l , we create a separate lagged regressor module for each of the m covariates x of lenght T . This allows to individually attribute the effect of each covariate on the predictions. Each lagged regressor module is functionally identical to the AR module, with the only difference being the inputs. Here, the p last observations of the covariate x are the inputs to the module (instead of the series y itself as in AR). The outputs are of identical form, each All considerations regarding the order, depth, sparsification, and interpretability of the AR-Net modules for each covariate are identical to the AR module, as discussed in section 2.1.3, if we substitute:

Future Regressors
To model future regressors, both past and future values of these regressors have to be known.
Given the set of future regressors as F ∈ R T ×n f , where n f is the number of regressors, the effect from all future regressors at time step t can be denoted by F (t) as in Equation 11, where d f stands for the coefficient of the model for future regressor f ∈ V. By default, future regressors have an additive effect, which can be configured to be multiplicative instead. where,

Events and Holidays
Effects from special events or holidays may occur sporadically. Such events are modelled analogous to future regressors, with each event e as a binary variable e ∈ [0, 1], signaling whether the event occurs on the particular day or not. For a set of events E ∈ R T ×ne with n e number of events and the length of the series T , the effect from all events at time step t can be denoted by E(t) in Equation 12, where z e denotes the coefficient of the model corresponding to the event e ∈ E. where, NeuralProphet allows the modelling of two types of events; 1) user defined 2) country specific holidays. Given a country name, its national holidays are automatically retrieved and added to the set of events E. Similar to seasonal effects, events can also be specified as either additive or multiplicative.
Additionally, for a given event at time t e , a window [t e − i, t e + j] of i + j days can be configured to be considered as special events of their own. For example, by setting a window of [−1, 0] for Christmas day, will allow the day before Christmas to have its own effect on the forecast. Hereby, a new variable is created for each day within the window around event and added to the set of events E.

Missing Data
Missing data is less of an issue when working with non-lagged input variables, as corresponding timesteps can simply be dropped. In doing so, one data sample will be lost per missing entry. With Auto-regression or lagged regressor modules in use, a missing data point would lead to h + p data samples being dropped due to h missing forecast targets and p missing lag values. For example, a single missing point leads to the loss of 13 samples for a model forecasting h = 3 steps ahead with AR(p = 10). Thus, we implemented a data imputation mechanism to avoid excessive data loss when working with incomplete data. The imputation mechanism follows the following heuristics: Data Imputation. If not specified or missing, events are assumed to not be happening. Missing Events are filled in with zeros, indicating their absence. All other real-valued regressor variables, including the time series itself, if autoregression is enabled, are imputed in a three step procedure.
First. The missing values are approximated by a bi-directional linear interpolation. Hereby, the last and the first known value before and after the missing values are used as anchor points for the interpolation. This is done for up to 5 missing values in each direction. If there are more than 10 missing values, they will remain NAN after this step. The amount of missing values to interpolate is user-settable.
Second. The remaining missing values are imputed with a centred rolling average. The rolling average is computed over a window of 30, and fills at most 20 consecutive missing values. The amount of missing values to fill with a rolling average is user-settable.
Third. If there are more than 30 consecutive missing values, the imputation algorithm aborts and instead drops all the missing datapoints.

Data Normalization
The type of normalization to apply to the time series can be set by the user. The available options are described in

Tabularization
We tabularize the time series data to create a pseudo independent and identically distributed dataset, as is required for SGD based training. Hereby, a data sample is created for each available time-stamp of the target time series 'y'. A sample includes a normalized time-stamp, and all inputs required by each configured module.
The modules Auto-regression and lagged covariates hereby extract the p values before a given time-stamp from their respective time series and store them in a vector of size p serving as input to the model. If forecasting multiple steps ahead with h > 1, the forecast targets also are stored in a h-sized vector for each time-stamp. This is not a memory-efficient approach to prepare the dataset to training. Hoever, we chose this procedure due to its simplicity and for its compute-efficiency when training.
Analogous, the Fourier terms are prepared by computing the 2k different sine and cosine transforms of the standardized time component, for each of the configured seasonal periodicities.

Training
One of the fundamental differences between Prophet and NeuralProphet is the fitting procedure.
Prophet utilizes L-BFGS (Nocedal & Wright, 2006), implemented in Stan (Carpenter et al., 2017) for fitting model parameters to the data. NeuralProphet retools Prophet, from the bottom up, replacing Stan with PyTorch, which is both flexible and easy to use.

NeuralProphet relies on a modern version of mini-batch stochastic gradient descent (SGD).
This fitting procedure is compatible with the vast majority of deep learning methods, can be scaled reliably to large datasets, and can accommodate more complex model components. Any model component which is trainable by SGD can be included as a module. This makes it easy for developers to extend the framework with new features, and to adopt new research.

Loss Function
The default loss function is the Huber loss, also known as smooth L1-loss, which can be seen in equation 13. Below a given threshhold β, it is equivalent to mean squared error (MSE). For values larger than β it is equivalent to the mean absolute error (MAE). Compared to a pure MSE loss, it is less sensitive to outliers and may help prevent exploding gradients Girshick (2015). We chose β = 1 as threshold. Users can however choose to alternatively opt for MSE or MAE loss, or any other available or self-implemented function matching the PyTorch loss function format as their loss function.

Regularization
We use a scaled and shifted log-transform of the absolute weight values as general regularization function. For a module with model weights θ it is parametrized as follows: where ∈ (0, ∞) sets the inverse of the offset and α ∈ (0, ∞) sets the scaling of the logtransform. Higher values of lead to a steeper curve for weights close to zero. The parameter could be described as a control of the eagerness to sparsify weights. Higher values of α lead to a flatter curve for weights of larger magnitude. The parameter α can be interpreted as a control of the acceptance of large weights. As default values, we suggest = 1 and α = 1 for most applications.
The regularization is applied in the per-module configured strength and added to the loss function, to be back-propagated. The regularization only commences after a specified percentage of training, by default, after 50%, and is thereafter linearly increased from zero to its full configured strength at the end of training.

Optimizer
Analogous to the loss function, any optimizer matching the signature of a PyTorch optimizer can be configured. As a reliable default, the AdamW optimizer (Loshchilov & Hutter, 2019) is used, unless otherwise specified. In our testing we observed AdamW to fit most reliably, with a slight tendency to overfit. We initialize AdamW with the configured learning rate, and set β = (0.9, 0.999), = 1e − 08, and a weight decay of 1e − 04.
As a fallback option, we also provide a classic SGD optimizer with momentum set to 0.9, and weight decay set to 1e − 4. In our testing, we found SGD to lead to better validation performance, but at the expense of occasional divergence. If a user is willing to fine-tune training related hyperparameters, SGD poses an attractive alternative to AdamW.

Learning Rate
As we do not want to require users to be experts in machine learning, we make the learning rate, an essential hyperparameter for any NN, optional. We achieve this by relying on a simple but reasonably effective approach to estimate a learning rate, introduced as a learning rate range test in Smith (2017).
A learning rate range test is executed for 100 + log 10 (10 + T ) * 50) iterations, starting at η = 1e − 7, ending at η = 1e + 2 with the configured batch-size. After each iteration, the learning rate is increased exponentially until the final learning rate is reached in the last iteration. Excluding the first 10 and last 5 iterations, the steepest learning rate is defined as the resulting learning rate.
The steepest learning rate is found by selecting the learning rate at the position which maximizes the negative gradient of the losses.

Batch Size
The batch size B is an optional parameter. If not user specified, the following heuristic will determine the batch size based on the length T of the dataset: B = 2 2+ log 10 (T ) B = min(T, max(16, min(256, B )))

Training Epochs
The number of training epochs N epoch is an optional parameter. If not user specified, the following heuristic will determine the number of epochs :

Scheduler
As all training related hyperparameters are automatically approximated, none of them are optimal, and may thus potentially lead to training issues. Nonetheless, a user expects the model to be reasonably fitted to the data without issues. In order to meet these expectations despite suboptimal hyperparameters, we rely on the training schedule called '1cycle' policy, which allows for 'superconvergence' of NN training, according to Smith & Topin (2018). Hereby, the initial learning rate η 100 is gradually increased up to the peak learning rate η, reached at 30% of training. Thereafter, the learning rate is annealed along a cosine curve down to η 5000 at the end of training.

Postprocessing
At the end of training or predicting, the normalized values on which the model operates are transformed back into the initial distribution. This, alongside with the many heuristically or procedurally set optional hyperparameters abstract a lot of the machine learning specific and time series specific decision making, making NeuralProphet approachable to forecasting practitioners.

Metrics
A few metrics are recorded by default: The configured loss function, RMSE, and MAE. Further, any user defined metrics will also be recorded over the training and validation runs.
where T = size of data

Forecast Presentation
When used to predict, the model returns a dataframe with a column for each forecast component. These are the additive (or multiplicative) contributions of each component. They are individually displayed for each forecasted step ahead. Each of the components is referring to its row's datestamp-located target and ordered by the age of the foregast. e.g. 'yhat3' refers to the predicted y-value for the current location, based on data available three steps ago. Similarly, 'ar2' refers to the AR-effect predicted for the current location, based on data available two steps ago.
There are different plotting utilities available for visualizing the forecast itself, the forecast decomposition, the model parameters, and for specific forecast horizons of interest.

Experimental Setup
The experiments in this work serve to contrast the capabilities of NeuralProphet with Prophet.
First, we compare the accuracy of the interpretable forecast components of Prophet and Neural-Prophet on a set of synthetic datasets. Second, we compare the performance of the two frameworks on a number of real-world datasets. In the following, these experiments are described in detail.

Decomposition Demonstration on Synthetic Data
These experiments contrast the ability of Prophet and NeuralProphet to decompose a time series into individual components. Since we are interested in the decomposition accuracy of the different modelling components, we must know their ground truth. Due to a lack of suitable realworld datasets with known underlying components, we create synthetic datasets as the sum of a selection of generated time series with added noise. This allows us to compare the estimated decomposed components of both models to the original components.

Synthetic Dataset Composition
We generate the individual component time series to be identical to the type of components that Prophet and NeuralProphet claim to model. Thus, ideally, both models should be able to perfectly decompose the time series into its components, with some error due to additive noise included in the synthetic datasets.
In Table 2 we display which components are included in the synthetic dataset of each experimental scenario. The comparison between the two models focuses on the accuracy of individual decomposed components. We also measure the overall accuracy and the compute time. Hereby, the emphasis is on the accuracy of the individual components as this quantifies the decomposition accuracy of the models.
As many real-world datasets exhibit auto-regressive properties or depend on lagged regressors, we also include three experiments with lagged components. However, only NeuralProphet claims to be able to model Auto-Regression and Lagged regressor components. In experiments including such lagged components, Prophet's prediction for these components will be marked as zero.
Experiment Trend Seasonality Events Future reg. AR Lagged reg.

S-AL
S-TSAL Due to the freedom to choose offsets in their additive components, the models may find interpretable components which match the dynamics of the underlying components, but shifted by an offset. Thus, before computing any component-wise metrics, we zero-center all components. from a random uniform distribution such that a j ∈ [0, 1) and b j ∈ [0, 1) and k = 5. For the experiments we create two independent seasonal components: monthly S 30 (t) and yearly S 365 (t).

Generated Component Time Series
Both are added to the aggregate series and both are individually evaluated in their decomposition accuracy.

E: Events. The event component is a binary series with 25 occurrences of an event at random
locations. Analogous to equation 12, we have E = (e 1 ) ∈ R T ×1 , with e 1 = (e 1 1 , ..., e 1 T ) ∈ R T where e 1 i = 1 if an event of type e 1 occurs at timestamp t i and 0 otherwise. The initial effect strength z e = 1 will effectively become z e ≥ 1 Ncomponents after scaling the main series.
F: Future Regressor. The future regressor component is generated as an AR(3) process with coefficients φ 1 = 0.2, φ 2 = 0.3, φ 3 = −0.5 and additive white noise. After scaling the series to a range of [0, 1], the regressor is added to the aggregate series with a weight of d f = 1, which will be Ncomponents after scaling the main series.
A: Auto-Regression. The Auto-regressive component is drawn from an AR(2) process with coefficients φ 1 = 0.3, φ 2 = 0.3, and additive white noise. This series is generated independent of the the other components. This is not identical to how NeuralProphet models the auto-regressive behavior.
NeuralProphet does not decompose the inputs to auto-regression module, but rather uses the raw time series values as inputs.
L: Lagged Regressor. First, an independent time series x is generated with an AR(2) process, analogous to future regressors. Next, we create the lagged regressor effect series L(t) ∈ R T , where each effect depends on a weighted combination of the last three observations from series x : where c 1 , c 2 , ..., c lag are weights from a random uniform distribution on (0, 1]. Finally, we utilize the series L(t) as the lagged regressor component. m: multiplicative mode. In multiplicative mode, we multiply each component element-wise by the trend. For example, in S-mTSEF:

Main Benchmark on Real-World Datasets
We evaluate the forecasting accuracy of Prophet and NeuralProphet on a diverse collection of datasets covering many univariate forecasting tasks. The datasets span data lengths from one hundred to half a million samples, while data record intervals span minutely to monthly frequencies.
This wide spectrum of univariate time series is in stark contrast to most benchmarking efforts focused on NN-based methods. The results represent realistic out-of-the-box performance across applications such as energy demand, tourism, environmental factors, and retail sales among others.
In this section, we present the selected datasets and the evaluation procedure in more detail.

Selected Datasets
The datasets are selected to represent a wide range of univariate datasets. Figure 1a shows that the lengths of datasets range from one hundred to one million samples, with most datasets having around 10,000 samples. The data observation frequencies are distributed from monthly to minutely sampling frequencies. Figure 1b bins minutely to half-hourly frequencies together and displays the count of each type of frequency, with hourly frequencies being the most common.
Our collection of testing datasets covers a wide variety of subjects, including environmental, energy-related, social, and economic data, with environmental and energy-related datasets being the most common (Figure 1c).
Quantitative criteria for inclusion in our collection of testing datasets were: Univariate time series of minutely to daily sampling frequency, length of 1000 samples or more, and permissive license. We also include the four time series which Prophet uses as demonstration time series. Two of those datasets do not meet the inclusion criteria, as they are of monthly frequency and much less than 1000 samples long. They are included nonetheless, to make sure the datasets are not biased towards the proposed model NeuralProphet, but rather towards the comparison model Prophet.
Further details on each dataset can be found in Table 8 (Appendix A).
(a) Boxplot of length for data frequencies.
Preprocessing. Most of the selected datasets do not contain missing values, but for those that do, we impute the values using linear interpolation unless otherwise described. However, for datasets which are mostly zero or near-zero-valued, such as data covering solar, parking, wind, and air quality subjects, we fill missing values with zeros. Missing time-stamps are also regarded as missing data.
Further, we manually preprocess each dataset and fix minor misalignments, such as time-stamps of 01:59 instead of 02:00 on data of hourly sampling frequency.

Benchmark Models
In this work, we compare the performance of time series forecasting packages which support full model automation, offer interpretability and scale to large datasets, without requiring domain knowledge in machine learning or time series. Prophet was the first such package, which we benchmark against its only fully featured successor, NeuralProphet.
Hyperparameters. Given that we aim to demonstrate the real world performance a non-expert forecasting practitioner can reasonably expect, we do not manually tune any hyperparameters. We only manually set hyperparameters in the synthetic benchmarking section if no automatic setting is available, or if needed to make models comparable.
We mention any manually set parameters explicitly when discussing results. Occurrences of manual model configuration, but not hyperparameter tuning, include: When studying the effect of adding auto-regression to the model, we set different numbers of lags p, which are all shown in the results. For multi-step ahead forecasts, we further need to set the forecast horizon h. In a reverse ablation study we further examine the impact of adding hidden layers to the model, hereby we also set the learning rate to the same value of 0.001 for all of the model configurations and all datasets.
These ablation studies of number of lags, number of forecasts and neural network configurations will also serve to quantify the sensitivity of NeuralProphet to parameter choices.

Forecast Evaluation
In the following, we describe the forecasting task and how the forecast performance is evaluated.

Forecast Horizon
In this study, we are interested in univariate time series forecasting covering typical forecasting tasks. Given a desired number of steps to forecast into the future, we predict the next few consecutive values of the series, starting from the present timestamp. The number of next consecutive values to be forecasted from the last observed value is also called forecast horizon. If auto-regression is configured, the model will regress over the last few observations of the series itself. We treat each given forecasting horizon as a separate forecasting task.
A forecasting task with a horizon of ten targets involves predicting the next ten targets starting at the current timestamp t, i.e. obtaining a series {ŷ t ,ŷ t+1 , ...,ŷ t+9 } containing the predictions over the forecast horizon. In our empirical evaluation over real-world datasets, we evaluate the forecast horizons of [∞, 1, 3, 15, 60]. Hereby, a horizon of ∞ refers to forecasting the entire test set at once, as is possible with model based on time-features only, such as Prophet and default NeuralProphet.

Expanding Origin Backtest
Evaluation procedures in time series forecasting are diverse in nature. Naming conventions have only been establish for evaluations of classical time series models (Tashman (2000)). Conventions for evaluating machine learning based time series forecasting models are currently being established.
We utilize a method based on a fusion of cross-validation and backtesting, which originate from machine learning and trading algorithm evaluation. This method, called expanding origin backtest, has recently been adopting by researchers and companies trying to evaluate machine learning based models on time series data (Uber (2018(Uber ( , 2020).
The procedure consists of an outer and an inner loop. The term expanding origin backtest mostly refers to the outer loop, but is often used to describe the whole procedure. The inner loop is usually based on a variation of the classical evaluation procedures. In our case, the inner loop is a rolling origin forecast evaluation with non-constant holdout size ( (Svetunkov, 2021)) with updating but without refitting of the model. In classical time series literature, refitting is referred to as reparametrization, while updating the model refers to updating the inputs to the model to produce the next forecast without refitting of model parameters. This means that a model will be re-fitted only k times per evaluation procedure.
A more verbose name for our forecast evaluation procedure could be: K-fold expanding origin backtest with an inner rolling origin multi-step forecast evaluation with updating but without refitting of the model.  Figure 3: Expanding Origin Backtest with 5 folds. Test sets contain 10% of the time series each.
Outer Loop: K-Fold Expanding Origin Backtest With Constant Holdout Size. Analogous to k-fold cross-validation, the overall dataset is replicated into k folds, each containing training and testing data, covering different parts of the overall dataset. The model to be evaluated is fitted separately on each fold's training data and evaluated on the same folds testing data. This signifies that the model Unlike k-fold cross-validation, the folds are not randomized. Each fold has a training set with values strictly occurring before the same fold's test set, a visualization can be seen in Fig. 3.
We perform a 5-fold expanding origin backtest with test sets of 10 percent and a forward rolling of the origin on 5 percent per fold. This means that each test fold will share 50 percent of their data with the previous fold, and that each train fold will grow by 5 percent, while the size of the test fold remains constant. Hereby, the first fold will use the first 70 percent of data for training and the next 10 percent for test. The fifth fold will utilize 90 percent of data for training and the last 10 percent of data for testing. Thus, a model is tested five times on 10 percent of the data, covering the last 30 percent of data over all tests.
Inner Loop: Rolling Origin With Updating But Without Refitting. The inner loop is repeated for each train-test fold. First, the model is fitted over the train set. The test set evaluation is identical to a rolling origin forecast evaluation with a non-constant holdout size ( (Svetunkov, 2021)) with updating but without refitting of the model.
Hereby, the forecast origin rolls forward by 1 step at a time, without re-training the model with the new data. This repeats through a single test fold shown in Figure 3 until the end of the test fold is reached. Auto-regressive models relying on lagged observations as inputs are updated with the new data point to produce the next forecast. If a model produces multiple forecasts per forecast origin, the average accuracy of all forecast steps is recorded. For models with unlimited forecast steps which do not depend on being updated with new observations, we forecast over the entire test set at once.
The single-forecast version of this procedure is also known as time series cross-validation in Hyndman & Athanasopoulos (2014), which we consider unfortunate naming, as it is not to be confused with the outer loop (Expanding Origin Backtest) which is more closely related to k-fold cross-validation.

Metrics
We calculate metrics on each i-step-ahead prediction versus corresponding actual values, and then average this error over all steps. For evaluation we use two metrics: Mean Absolute Scaled Error (MASE) and Root Mean Squared Scaled Error (RMSSE), as they are defined in Hyndman & Koehler (2006). MSSE is given as the RMSE of the evalutated method divided by RMSE of the Naïve prediction. This is the squared analogue of mean absolute scaled error (MASE). Hereby, Naïve refers to predicting the next observation to be identical to the previous value observed.
MASE and RMSSE metrics allow to evaluate the performance of a method compared to Naïve forecasting directly, as a value smaller one signifies an improvement over Naïve. Further, as the metrics are scaled quantities, different models can be compared across different datasets.
where T = size of train data J = size of test data

Computational Resources
Unless otherwise mentioned, all models were trained and evaluated on a regular laptop with an 8-core CPU (2020 Macbook Pro with M1 chip).

Demonstration of Interpretable Component Decomposition
The results are displayed separately for experiments S-TS, S-EF, S-TSEF, and S-mTSEF without lagged components and for experiments S-AL, S-TSAL, and S-TSEFAL including lagged components. In the former, we expect both models to perform identically, and in the latter, we expect NeuralProphet to perform better.
On data without lagged components, NeuralProphet performs identical or slightly better than Prophet, as can be seen in figure 4 and figure 5a. The overall goodness of fit on the 'y' component (the sum of the components) is comparable or slightly better. The biggest difference is seen on S-mTSEF with multiplicative components, where NeuralProphet has a significantly lower error and standard deviation of error.
On experiments involving lagged components, NeuralProphet performs significantly better than Prophet, as seen in figure 4 and figure 5a. Most of the difference comes from a better modelling of lagged regressors.  The fit on the series itself is not as relevant as the component-wise decomposition accuracy, as the overall goodness of fit does not protect from overfitting. The component-wise accuracy, is an adequate measure, as a model that can capture the true underlying dynamics will generalize well.

Component-Wise Decomposition Accuracy
Without Lagged Components. When compared on the component-wise accuracy NeuralProphet performs identical or better than Prophet, as can be seen in figure 5a. Both models fit the monthly seasonality, events and future regressors near perfectly. Both models mostly struggle to find a good fit for the yearly seasonality. This is likely explained by the limited number of full annual periods in the data to fit the model on. A significant difference is observed in the goodness of fit for the trend component, where Neuralprophet exhibits a fraction of the error of Prophet.
Including Lagged Components. When modelling lagged components, NeuralProphet offers a significant improvement decomposing lagged regressors. However, the AR component is not significantly more accurate than Prophet's zero-prediction. This is due to the fact, the NeuralProphet does not decompose the inputs to the AR module, but rather uses the raw time series values, while in this experiment, we evaluate the performance based on an independently added AR process component.
On real applications, such a distinction largely does not matter.
Zooming in on the performance for individual experimental setups in Fig. 7 (Appendix B), we   Table 3 shows that training for NeuralProphet is on average 4.0 times slower compared to Prophet. However, both exhibit a large standard deviation, with Prophet having a long tail upwards and NeuralProphet downwards. Prediction time exhibits the inverse relationship, with NeuralProphet computing 13.5 times faster compared to Prophet. NeuralProphet has a tight spread of the prediction times with no upward outliers, while Prophet has a large standard deviation with a long tail upwards, as seen in figure 6.  These observations mean that NeuralProphet may take some extra resources to be fit, but will need significantly less resources in deployment. Further, as NeuralProphet is faster in computing predictions by a magnitude, it makes it possible to deploy it in time-sensitive applications where the next prediction needs to be reliably computed in a fraction of a second. Thanks to being fully implemented in PyTorch, the NeuralProphet model may be parallelized in the future, allowing it to be deployed on GPUs. We would expect this to close the gap in training time, and the gap in prediction time to widened.

Main Benchmark Results
The overall results are displayed in table 4. We would like to highlight three observations: First, Prophet and NeuralProphet perform near identical in their default mode. However, both perform significantly worse than Naive one step ahead. Second, NeuralProphet forecasts improve substantially when configuring any amount of lags, consistently performing better than Naive forecasting one step ahead. More lags generally lead to a better performance. Third, when forecasting multiple steps ahead, even when we increase the horizon to 60 steps, NeuralProphet with any number of lags still performs better than Prophet or NeuralProphet in default configuration.
We also observe that the specific number of lags chosen has a small impact on the forecast accuracy, suggesting that the model is not sensitive to optimal parameter choices.
All model parameters which were set for any of the models are explicitly displayed in    subject type, where only a NN-based model outperforms Prophet and default NeuralProphet.
Measurements of human activity tend to have strong seasonality and have a quickly fading impact of local context. When forecasting into the distant future, statistical patterns tend to dominate human activity data.

Task-dependent Strength of Prophet and NeuralProphet
In Table 7 we present our subjective advice on which model to use depending on the task at hand. These are qualitative suggestions based on the empirical evidence we found when evaluating the two models against each other.  we suggest to use for the given task.

Conclusion
NeuralProphet is the first hybrid forecasting framework that meets the industry standards for explainability and simplicity-of-use set by Facebook Prophet. As NeuralProphet is based on PyTorch and trained with standard deep learning methods, any algorithm trainable by mini-batch Stochastic Gradient Descent can be included as a module, which makes it easy for developers to extend the framework with new features, and to adopt new research.