# Comparing Prophet and Deep Learning to ARIMA in Forecasting Wholesale Food
Prices

Lorenzo Menculini111Corresponding author<EMAIL_ADDRESS>Idea-re
S.r.l., Perugia, Italy  Andrea Marini  Idea-re S.r.l., Perugia, Italy
Massimiliano Proietti  Idea-re S.r.l., Perugia, Italy  Alberto Garinei
Department of Engineering Sciences, Guglielmo Marconi University, Rome, Italy
Idea-re S.r.l., Perugia, Italy  Alessio Bozza  Cancelloni Food Service S.p.A.,
Magione (PG), Italy  Cecilia Moretti  Marcello Marconi  Department of
Engineering Sciences, Guglielmo Marconi University, Rome, Italy  Idea-re
S.r.l., Perugia, Italy

###### Abstract

Setting sale prices correctly is of great importance for firms, and the study
and forecast of prices time series is therefore a relevant topic not only from
a data science perspective but also from an economic and applicative one. In
this paper we examine different techniques to forecast sale prices applied by
an Italian food wholesaler, as a step towards the automation of pricing tasks
usually taken care by human workforce. We consider ARIMA models and compare
them to Prophet, a scalable forecasting tool by Facebook based on a
generalized additive model, and to deep learning models exploiting Long
Short–Term Memory (LSTM) and Convolutional Neural Networks (CNNs). ARIMA
models are frequently used in econometric analyses, providing a good benchmark
for the problem under study. Our results indicate that ARIMA models and LSTM
neural networks perform similarly for the forecasting task under
consideration, while the combination of CNNs and LSTMs attains the best
overall accuracy, but requires more time to be tuned. On the contrary, Prophet
is quick and easy to use, but considerably less accurate.

##  1 Introduction

The main aim of firms is profit maximization. To achieve this goal, the
constant updating and forecasting of selling prices is of fundamental
importance for every company. Although the digital transformation is a
phenomenon that is involving all companies, from small to large, many of them
still update prices by hand through logics that are not always clear nor
objective and transparent, but rather based on the experience and expertise of
those in charge of updating the price list. On the other hand, the automation
of price prediction and update can provide a strong productivity boost by
freeing up human resources, which can thus be allocated to more creative and
less repetitive tasks. This also increases the morale and commitment of
employees; it also speeds up the achievement of goals, and improves accuracy
by minimizing human errors. Subjectivity is also reduced: once the operating
criteria have been established, forecast algorithms will keep behaving
consistently. This in turn means an improvement in compliance.

Besides the automation of price updates, the prediction of the sales prices
charged to customers in the short term also holds great value. In general,
organizations across all sectors of industry must undertake business capacity
planning to be efficient and competitive. Predicting the prices of products is
tightly connected to demand forecasting and therefore allows for a better
management of warehouse stocks. The current economic crisis caused by COVID-19
has highlighted the value of such management optimization, stressing the
importance of the companies’ ability to minimize inventory reserves and just-
in-time production models. Forecast models considered in this paper can
contribute to keeping the difference between wholesale purchase prices and
company’s sales prices under control, in view of maximizing the gross
operating income. They can therefore help companies avoid the risk of
incalculable losses and ultimately improve their contractual capacity.

The present work proposes to deal with these topics by investigating and
comparing different price forecasting models. The specific task we consider is
that of predicting the prices of three food products sold by a medium/large-
size local wholesaler based in central Italy. In such way, we investigate the
predictability of wholesale prices, comparing the performance of traditional
econometrics time series forecasting models with Facebook’s Prophet and
machine learning models. The main goal of this paper is therefore to develop a
forecasting model that could represent a first step towards the automation of
the price-setting process, thus effectively aiding the work of company
employees. In this way, it aims to be of practical use to companies for the
maintenance and management of price lists. Scalability and flexibility of the
models presented in this paper are also an important point: for the sake of
simplicity, we have applied the models to three different products, but we
underline that the same models and algorithms can be easily applied to any
product.

Time series forecasting has always been a major topic in data science with
plenty of applications. For a general review of some of the most used tools,
see for example [1]. Well-known traditional econometric methods are not always
appropriate to study and forecast big and noisy time series data. This has
generated particular interest in machine learning methods, bolstering data
driven approaches that include a wide range of methods that have the advantage
of not relying on prior assumptions and knowledge on data. See for example [2,
3, 4, 5] for reviews focusing on the application of deep learning [6] to time
series. Long short–term memory (LSTM) networks [7] and convolutional neural
networks (CNNs) [8] are almost ubiquitous in time series forecasting with
machine learning. CNNs are even more commonly used for image recognition and
feature extraction. However, the forecasting accuracy of standalone CNNs can
be relatively low [9].

The literature concerning economic time series prediction employing various
methos – for classical to artificial intelligence ones – is very rich.
Nevertheless, although we believe automatic updating mechanisms and
forecasting of sale prices are of uttermost relevance, the literature on these
topics is not as developed as one would expect. Most studies focus primarily
on the implementation of models for the analysis and forecasting of general
price levels (inflation) or commodity and stock market prices. The forecasting
of food prices in China was considered by the authors of [10, 11]. In
particular, Zou et al. [11] compared the performances of ARIMA, neural
networks (NNs) and a combination of the two to forecast wheat prices in the
Chinese market. Their findings showed that, overall, NNs perform best at the
task. Neural networks were also employed in [12] to forecast monthly wholesale
prices of two agricultural products. Ahumada and Cornejo [13] considered a
similar problem, also taking into account possible cross-dependencies of
different product prices. In [14] the author focused on sales forecasting
using machine learning models, a topic similar to the one considered in the
present paper. For more recent work on forecasting commodities prices see
[15], where the authors forecasted gold prices, and [16] where the Levenberg-
Marquardt Backpropagation (LM-BP) algorithm was applied to stock prices
prediction. Other authors used machine learning methods for inflation
forecasting [17, 18], also in comparison with more classical econometric
models [19]. Xue et al. [20] recently presented a high-precision short-term
forecasting model for financial market time series employing deep LSTM neural
networks, comparing them with other NN models. Their results showed that LSTM
deep neural networks have high forecasting accuracy for stock market time
series. In 2020, Kamalov [21] evaluated multilayer perceptrons, CNNs and LSTM
neural networks to forecast significant changes in stock prices for four major
US public companies, showing that these three methods yield better results
when compared to similar studies that forecast the direction of price change.
For models similar to the ones considered in this work and applied again to
stock indexes forecasting, see [22]. Hybrid ARIMA/neural network models where
instead studied by the authors of [23]. Stock prices have also been forecasted
using LSTMs in conjunction with the attention mechanism [24]. Machine learning
models using LSTMs and CNNs are of widespread use in time series forecasting,
well beyond the financial and economic realm. For recent work on time series
forecasting using machine learning outside the economic and financial area see
[25], an application to COVID-19 spreading forecasting, and [26] for an
application of deep learning to Influenza prevalence forecasting.

In this paper we compare the performance of standard Autoregressive Integrated
Moving Average (ARIMA) models [27], which we take as a benchmark, to Prophet –
a forecasting tool developed by Facebook and based on a Generative Additive
Model (GAM) [28] – and machine learning models exploiting LSTMs, both on their
own and in combination with CNNs. ARIMA univariate models are considered a
standard reference model in econometrics. The compared models are rather
different, as are the datasets that they accept in input, making the
comparison interesting. On one hand, Prophet’s driving principles are
simplicity and scalability; it is specifically tailored for business
forecasting problems and handles missing data very well by construction. On
the other, the NN models we construct allow for a multivariate regression,
fully exploiting all the collected data, but also require some data pre-
processing, as does ARIMA. Prophet has been compared to ARIMA models for the
prediction of stock prices [29] and bitcoin [30].

Our results indicate that the combination of CNNs and LSTMs yields the most
accurate results for all the three products, but require the longest and
computationally more expensive tuning. On the contrary, Prophet performances
were not brilliant, but model tuning and data preparation were particularly
quick. ARIMA and LSTM-only neural networks showed good performances both in
terms of accuracy and time required for model selection and training.

The rest of the paper proceeds as follows. Section 2 introduces the dataset
features, discussing its properties and some pre-processing steps that were
taken on it; it also briefly presents the three models under consideration,
their set-up and tuning. In Section 3 we give the results obtained with the
three approaches and compare them. We conclude in Section 4 by discussing the
results of this study and providing an outlook on future perspectives in the
light of the paper findings.

##  2 Materials and Methods

###  2.1 Dataset description and preparation

For this study we had access to a dataset comprising a total of approximately
260,000 food order records, reporting the following information: _date of
order_ , _order number_ , _unit price_ , _article code_ , _sold quantity_ ,
_customer code_ , _offer_ (if present) and _offer type_ , _unitary cost_. The
records were collected by the wholesaler in a period ranging from year 2013 to
2021. For the study conducted in this paper, we decided to focus on the three
products with the most records, namely _Carnaroli rice 1kg $\times$ 10_
(henceforth product 1), _Gorgonzola cheese 1/8 of wheel 1.5 kg_ (product 2)
and _Cured aged ham 6.5 kg_ (product 3). The forecasting task considered in
this work was to predict the average selling price for the following week, for
each of the selected products.

As a first thing, we chose to leave out all data following the outbreak of the
COVID-19 pandemic. This was motivated by the huge impact that the lockdowns
and restrictions imposed by the authorities had on the food and catering
sector, introducing a major shock in sales trends at all scales. Therefore, we
excluded all records dated later than March 9, 2020 (last day before the first
national lockdown in Italy).

A preliminary data analysis revealed that the dataset contained a good number
of outliers: for some of them, it appeared evident that this was due to
incorrect typing of the product sale price. To improve the quality of the
dataset, we calculated the _z-score_ of each record based on its price as
$z=(p-\bar{p}^{(w)})/\sigma^{(w)}\,,$ where $p$ is the unit sale price and
$\bar{p}^{(w)}$ and $\sigma^{(w)}$ are the mean and standard deviation for the
selected product, weighted by the quantity sold in each order. Then, we
filtered out all records with $|z|>4$. Figure 1 shows the price distribution
for product 2, after the filtering.

Figure 1: Price distribution (after _z-score_ filtering) for product 2: to the
left, boxplot, and to the right, violin plot.

In view of the subsequent time series forecasting and as a further step in
dealing with data inaccuracies, we decided to resample the dataset with weekly
frequency. This was done by cumulating the number of orders in each window and
calculating the average sale price for each week. For later use in neural
network models, when resampling data we also kept track of the following
fields in the dataset: _number served customers_ , _number of orders_ ,
_number of orders on sale_ , (weighted) _average product cost_ , and
(weighted) _price standard deviation_. Table 1 summarizes the main features of
the resampled price time series for each of the products. In Figures 2, 3 and
4 we display the time series of sale prices and sold quantities after
resampling. All prices, here and everywhere in the paper, are intended in
euros (€).

Table 1: Mean and standard deviation of the weekly prices time series. Product | 1 | 2 | 3
---|---|---|---
Mean (€) | 1.99 | 7.64 | 7.27
Std (€) | 0.36 | 0.40 | 0.26

Figure 2: (a) Unit price and (b) sold quantity time series for product 1 after
resampling with weekly frequency.

Figure 3: (a) Unit price and (b) sold quantity time series for product 2 after
resampling with weekly frequency.

Figure 4: (a) Unit price and (b) sold quantity time series for product 3 after
resampling with weekly frequency.

We then split the weekly dataset in the following way:

  * •

_Training dataset_ : years 2013–2017,

  * •

_Validation dataset_ : year 2018,

  * •

_Test dataset_ : years 2019–March 2020

This choice was made in order to make the training set as large as possible,
while also having validation and test sets that cover at least one year of
data.

As can be seen in Figures 2, 3 and 4, even after resampling the price time
series have missing periods. Therefore, they are not suited for being given as
input to ARIMA and neural network forecasting models, as they require evenly
spaced data. To overcome this problem, we adopt the following strategy for all
products:

  1. 1.

Since the time series for all products have a long window with no data in the
second half of year 2013, we do not consider this empty period and start right
after it;

  2. 2.

When occasional weeks with no data occur, we take the average of the preceding
and following week prices and interpolate.

In this way we were able to fill in all empty weeks – in fact, after
resampling the missing records were very sparse. Note that the above procedure
is only necessary for preparing the dataset for ARIMA and NN models, as
Prophet has no problems in handling missing datapoints. The size of the
datasets for each product, both before and after removal of empty periods, is
summarized in Table 2.

Table 2: Dataset size (# of datapoints) for each product and forecasting model. Product | Model | Train | Valid | Test
---|---|---|---|---
1 | Prophet | 240 | 52 | 62
ARIMA & NN | 211
2 | Prophet | 241 | 52 | 62
ARIMA & NN | 211
3 | Prophet | 242 | 52 | 62
ARIMA & NN | 212

###  2.2 ARIMA models

ARIMA models [27] are among the most simple and used econometric approaches to
univariate time series modeling. In this work, we implemented non-seasonal
ARIMA models, neglecting the modulation effects of holidays and using
therefore pure trend lines.

In econometrics, it is quite customary when dealing with price variables to
transform prices through a logarithmic map, since this generally leads to
better results. We decided to follow this approach when using the ARIMA
modelling, thus working with $\log(p_{t})$ in the model fitting. As a first
step we checked the stationarity properties of the time series. We performed
the _Augmented Dikey-Fuller_ unit root test using the built-in method in the
_statsmodels_ Python package. The results we obtained are qualitatively
similar for all the three products we considered: for the $\log(p_{t})$ time
series one cannot reject the null hypothesis of the presence of a unit root,
signalling the non-stationarity of the series. First differencing the series,
i.e. considering $\Delta\log(p_{t})=\log(p_{t})-\log(p_{t-1})$, makes it
stationary. Thus the $\log(p_{t})$ series are integrated of order one, and
accordingly the models we considered are ARIMA$(p,1,q)$.

In order to have a rough indication on the AR orders, $p$’s, and on the MA
orders, $q$’s, we computed the sample autocorrelation function (ACF) and the
partial autocorrelation function (PACF) for $\Delta\log(p_{t})$. Recall that

  * •

for an exact MA($q$), ACF is zero for lags larger than $q$;

  * •

for an exact AR($p$), PACF is zero for lags larger than $p$.

As an example we show the plots of these functions for product 2 in Figure 5.

Figure 5: Autocorrelation function and partial autocorrelation function for
$\Delta\log(p_{t})$ of product 2.

In the ARIMA model selection and fitting procedures, we used a different
dataset splitting scheme with respect to the one described above, in that the
_training_ set comprised years 2013–2018 (i.e. the union of the former
training and validation sets). This is since we decided not to use the
validation set to selected the hyperparameters $p$ and $q$, instead exploiting
the _Bayesian Information Criterion_ (BIC) as a metric for model comparison
[31]. Hence we took into account different combinations of $p$ and $q$ around
the values suggested by the ACF and PACF plots, and eventually we selected the
model with least BIC.

###  2.3 Prophet

Prophet is an open-source tool provided by Facebook Inc., available both in
Python and R. For the current analysis, the Python package (with Python 3.9)
was used. As explained by the authors [32], the idea leading to Prophet was to
develop a flexible forecasting tool which is easy to both use and tune. The
underlying model features a decomposable time series with three components:
growth (or trend) $g(t)$, seasonality $s(t)$ and holidays $h(t)$ (if present).
In the present case, there are no obvious holidays to consider, as the
wholesaler’s customers are mainly restaurants and hotels, which tend to stay
open during holidays. The time series is therefore decomoposed as

$y(t)=g(t)+s(t)+\epsilon_{t}\,,$ (1)

where $\epsilon_{t}$ encodes variations that are not taken into account by the
model, and which are assumed to be normally distributed [32]. The Prophet
model can be seen as a GAM [28]. In this framework, forecasting is phrased as
a curve-fitting task, with time as the only regressor, so the model is
univariate.

The trend function adopted for the problem under study is a piecewise linear
function written as

$g(t)=\big{(}k+\sum_{i:t>s_{i}}\delta_{i}\big{)}t+\big{(}m+\sum_{j:t>s_{j}}\gamma_{j}\big{)}\,,$
(2)

where $k$ is a scalar coefficient, $s_{i}$ are the _trend changepoints_ –
_i.e._ $S$ times $s_{1}$, $s_{2}$, … , $s_{S}$ at which the angular
coefficient of the trend is allowed to change – $\delta_{i}$ are the rate
adjustements, and $\gamma_{j}=-s_{j}\delta_{j}$ are parameters used to make
the function continuous. The algorithm starts with a number $S=25$ of
potential changepoints, placed in the first $80\%$ of the time series in order
to avoid responding to fluctuations in the last part of the series. Then, the
actual changepoints are selected by putting a sparse prior of the kind
$\delta_{j}\sim\text{Laplace}(0,\tau)$, with $\tau$ (tunable hyperparameter)
regulating222Using a Laplace prior is equivalent to L1-regularization. the
magnitudes rate adjustments. A larger $\tau$ means the model has more power to
fit trend changes.

As for seasonality, Prophet accounts for it using Fourier series, namely

$s(t)=\sum_{n=1}^{N}\left(a_{n}\cos\left(\frac{2\pi
nt}{P}\right)+b_{n}\sin\left(\frac{2\pi nt}{P}\right)\right)\,.$ (3)

Since we considered weekly data with no other obvious expected seasonality
effects but yearly ones, we had $P=365.25\text{d}$. For weekly seasonality,
the truncation parameter is set to $N=10$ by the authors of [32] when
modelling yearly seasonality, and we follow this specification. When
performing the fit, a smoothing prior $\beta\sim N(0,\sigma^{2})$ is imposed
on the $2N$ components $\beta=(a_{1},\ldots,a_{n},b_{1},\ldots,b_{n})^{T}$,
with $\sigma$ a second hyperparameter (essentially acting as an
L2-regularization parameter).

Prophet fits its GAM using the L-BFGS quasi-Newton optimization method of [33]
in a Bayesian setting, finding a maximum _a posteriori_ estimate.

###  2.4 Neural Networks

The advent of artificial intelligence, in particular machine learning, has led
to the development of a set of techniques that have proved to be very useful
in many different areas. One breakthrough has certainly been deep learning
[6], which has revolutionized our way of handling and exploiting information
contained in data. Deep learning can effectively detect and model hidden
complexity in data, automatically extracting features that should otherwise be
extracted manually by dataset inspection.

A standard choice when facing problems involving time series is that of using
LSTM neural networks, a kind of recurrent neural networks (RNNs) devised by
Hochreiter and Schmidhuber in 1997 [7]. Like all RNNs, they can by
construction handle data endowed with temporal structure, while also providing
a way to deal with the vanishing gradient problem [34]. Here we will describe
the application of LSTM NNs to the problem under study, both on their own and
in combination with CNNs. Indeed, standard LSTM NNs for time series
forecasting can be enriched with one-dimensional convolutional layers that
sequentially apply a unidimensional filter to the time series. Convolutions
can be seen as non-linear transformations on the time series data. This
enhances the model’s capability to learn discriminative features which are
useful for the forecasting and that can be fed to the LSTM layers that follow.
The models developed in this work were trained and tested using Python 3.9 and
TensorFlow 2.5.

Unlike in the ARIMA and Prophet case, with NNs we can exploit a larger
fraction of the information available in the dataset by setting up a
multivariate regression. However, since dates cannot be used as input
variables to a NN, we made an addition to the fields listed in Section 2.1,
performing a time embedding to provide information about seasonality. We did
this by adding the columns

week_cos $\displaystyle=\cos(2\pi w/52.1429)\,,$ (4a) week_sin
$\displaystyle=\sin(2\pi w/52.1429)\,,$ (4b)

where $w$ is the week number ($w=0,1,\ldots,52$). Therefore, we had a total of
9 input columns that were passed to the NN models. A sample of the input
dataset for one product is shown in Table 3.

Table 3: A slice of the dataset used to generate input data for the NN models. quantity | customers | orders | on sale | cost | week_cos | week_sin | p_std | avg_price
---|---|---|---|---|---|---|---|---
10 | 1 | 1 | 0 | 1.40 | 0.990 | -0.141 | 0 | 1.41
0 | 0 | 0 | 0 | 1.40 | 1.000 | -0.0214 | 0 | 1.55
70 | 6 | 6 | 0 | 1.40 | 0.993 | 0.120 | 0.0690 | 1.69
220 | 17 | 18 | 0 | 1.40 | 0.971 | 0.239 | 0.0580 | 1.75
230 | 14 | 15 | 0 | 1.39 | 0.935 | 0.353 | 0.0685 | 1.86

To construct the actual training dataset (made of a tensor $\bf{x}_{t}$ and a
scalar $y_{t}$, for each time $t$) that could be fed to LSTM neural networks,
we then performed the following steps:

  * •

reshape data so that at each time $t$, $\bf{x}_{t}$ is a $n\times 9$ tensor
containing the $n$ last values of each time series in Table 3;

  * •

set $y_{t}=\Delta_{t+1}=p_{t+1}-p_{t}$ as the variable to be used in the cost
function.

In this way, the model learns to predict $y$ (price variation) at each time
based on information about the last $n$ timesteps. Predicting the increment of
the quantity of interest instead of the quantity itself is a well-known way to
improve performance when training multivariate machine learning models.
Moreover, we checked through the _Augmented Dikey-Fuller_ test that the
$\Delta_{t}$ time series was stationary. The number $n$ of timesteps used
depends on the model and will be specified later.

The NN models tried in this paper for predicting product prices fall in two
classes: those using only LSTM layers and those with CNN layers before the
LSTM. We denoted these classes A and B, respectively.

##  3 Results

In this section we report results obtained with the three different approaches
to forecasting – namely ARIMA, Prophet and deep learning – studied in this
work.

###  3.1 ARIMA results

We considered first ARIMA models, to provide a standard performance benchmark
with which to compare the other models developed in the rest of this work.

The best ARIMA models for the three products are given in Table 4. As outlined
in Section 2.2, they were selected by considering the series of the price
logarithms $\log(p_{t})$ and using a least BIC criterion [31]. For the sake of
comparison with the other models, we transformed back the $\log(p_{t})$ series
to the $p_{t}$ series to compute the root mean squared error (RMSE) between
the predicted and observed increments

$\hat{\Delta}(t)=\hat{p}_{t}-p_{t-1}\,,\qquad\Delta(t)=p_{t}-p_{t-1}\,,$ (5)

$\hat{p}_{t}$ being the predicted price at time $t$.

Table 4: Selected ARIMA models for the three products and associated RMSE on the entire training+validation set. Product | selected model | (tr+v)RMSE
---|---|---
1 | ARIMA(2,1,0) | 0.097
2 | ARIMA(0,1,2) | 0.232
3 | ARIMA(3,1,1) | 0.211

In all cases, we checked also that the Ljung-Box statistics [35, 36] for 1-,6-
and 12-lag residual autocorrelations do not reject the null hypothesis, so the
residuals can be considered approximately white noise.

The results obtained on the test set by the selected ARIMA models are
summarized in Table 5. We also report values for the MAE (mean absolute error)
and MAPE (mean absolute percent error).

Table 5: Performance of the best ARIMA models on the test set. The MAPE is
computed for price time series $p_{t}$, not for the $\Delta(t)$ time series.
ARIMA results
---
Product | 1 | 2 | 3
RMSE | 0.0758 | 0.173 | 0.215
MAE | 0.0581 | 0.132 | 0.159
MAPE | 0.0348 | 0.0178 | 0.0222

###  3.2 Prophet results

####  3.2.1 Prophet grid search

As suggested in the Prophet documentation and reviewed in Section 2.3, one can
tune the $\tau$ (trend changepoints prior scale) and $\sigma$ (seasonality
prior scale) hyperparameters so that the model fits data as well as possible.
We did so by performing a grid search over $\tau$ and $\sigma$ in the
following way: for each combination of $\tau\in\\{0.005,0.01,0.05,0.1,0.5\\}$
and $\sigma\in\\{0.01,0.05,0.1,0.5,1,2\\}$, we started by fitting the model
over the training dataset, and predicted the price for the following week
(first datapoint in the validation set). We calculated the squared error
between predicted and observed price. Then, we moved on to the second
datapoint in the validation set, performed a new fit using also the first
validation set datapoint, and predicted the price for the following week. The
whole process was repeated until the validation dataset was exhausted. For
each product, the configuration with least RMSE was selected, yielding the
results shown in Table 6.

Table 6: Prophet grid search results. vRMSE indicates the RMSE on the validation set, calculated as described in the main text. Product | $\tau$ | $\sigma$ | vRMSE
---|---|---|---
1 | 0.5 | 0.01 | 0.0831
2 | 0.1 | 0.01 | 0.293
3 | 0.5 | 1.0 | 0.215

####  3.2.2 Prophet forecasting

After selecting the best values of the parameters for each product, we
employed them to specify the Prophet model in the test phase. Testing took
place in the following way: we started by fitting the model over the entire
training plus validation dataset, predicting the first data entry in the test
dataset and calculating the squared error. We repeated the procedure for all
entries in the test dataset, each time using all previous history, and
calculated the achieved RMSE at the end of the process. Note that by employing
this procedure, also the test dataset is progressively used to fit the model.
We plot in Figure 6 the result of the fit over the entire dataset for product
1, i.e. the function that would be used to predict the unit price for the
first week following the end of the test dataset.

Figure 6: Prophet fit over the entire product 1 dataset. The actual fit is
shown in light blue, while trend changepoints are shown in red.

Figure 7 shows instead the plot of the predicted and observed increments –
$\hat{\Delta}(t)$ and $\Delta(t)$ as defined in Eq. 5 – in the case of product
2. The performance of Prophet in forecasting the weekly price time series is
summarized in Table 7. As done for ARIMA models, besides the RMSE parameter
used in the fine tuning, we also report values for the MAE and MAPE.

Figure 7: Comparison of Prophet forecasts and observed data for price
variations in the test set, for product 2. The corresponding RMSE is 0.220.
Table 7: Performance of the tuned Prophet models on the test dataset. Prophet
results
---
Product | 1 | 2 | 3
RMSE | 0.0812 | 0.220 | 0.350
MAE | 0.0694 | 0.165 | 0.301
MAPE | 0.0414 | 0.0224 | 0.0424

###  3.3 Neural Networks results

####  3.3.1 NN grid search

As a third forecasting tool, we studied deep neural networks. Two different
classes of models – class A and B as introduced in Section 2.4 – were
analyzed. We performed two different grid searches to select the best model in
each class, as we now briefly describe. The common features between the
classes were: the usage of an MSE cost function and of the _Adam_ optimization
algorithm [37], with learning rate (or _stepsize_)
$\alpha\in\\{0.0005,0.001\\}$; the adoption of an early stopping procedure,
monitoring the cost function on the validation set with a patience of 5
epochs, while also setting an upper bound of 150 training epochs. Finally,
data were standardized using a _z-score_ normalization333Specifically, the
_MinMaxScaler_ function of the _scikit-learn_ Python package was employed.,
and we used a batch size of 32.

For class A, we trained NN models with the following architecture and
hyperparameters: a number $l\in\\{1,2,3\\}$ of LSTM layers with
$n_{u}\in\\{32,64,96\\}$ neurons each and _normal Glorot_ weight
initialization [38]. Each LSTM layer was followed by a dropout layer, with
dropout rate $r\in\\{0.1,0.2,0.3\\}$. The output layer consisted of a single
neuron with linear activation function, again with normal Glorot
initialization. For this class of models, we used a number of timesteps $n=4$.
The grid search over the hyperparameters $l$, $n_{u}$, $r$, $\alpha$ was
performed by initializing and training each model ten times, monitoring the
cost function on the validation set and recording the best result obtained for
every configuration. The best performing models for each product are reported
in Table 8.

Table 8: Results of the grid search on class A (LSTM only). trRMSE and vRMSE are the RMSEs computed on the training and validation sets, respectively. Product | $l$ | $n_{u}$ | $r$ | $\alpha$ | trRMSE | vRMSE
---|---|---|---|---|---|---
1 | 3 | 32 | 0.1 | 0.001 | 0.0964 | 0.0612
2 | 3 | 64 | 0.3 | 0.001 | 0.163 | 0.173
3 | 1 | 96 | 0.1 | 0.0005 | 0.156 | 0.142

The second class of models, class B, consisted in a combination of CNN layers
and LSTM layers, as done for example in [15] and [22]. We added two one-
dimensional convolutional layers, with a pooling layer in between. Each of the
convolutional layers had $f\in\\{10,20,30\\}$ _output filters_ , _kernel size_
$k_{s}\in\\{2,4\\}$ and _relu_ (rectified linear unit) activation function.
Moreover, we tried to use _same_ , _causal_ or no padding in each of the
conv1D layers: we dub the corresponding hyperparameters pad1, pad2. The 1D
_average pooling_ layer had _pool size_ equal to 2 and no padding. This block
was then followed by the same LSTM layers as for the first class of models.
The reason behind adding CNN layers is that they can help make better use of
the data history, improving the ability of LSTM networks to learn from long
series of past data. The grid search was indeed performed with varying numbers
of timesteps $n\in{4,8,12}$. As for the previous class, each model was
initialized and trained twice. Results for the grid search over the
hyperparameters $l$, $n_{u}$, $r$, $\alpha,f,k_{s},n$ are shown in Table 9.
Figure 8 shows the trajectory of training and validation cost functions for
the best-performing model in the case of product 1.

Table 9: Results of the grid search on class B (CNN + LSTM). trRMSE and vRMSE are the RMSEs computed on the training and validation sets, respectively. Product | $l$ | $n_{u}$ | $r$ | $\alpha$ | $f$ | $k_{s}$ | pad1 | pad2 | $n$ | trRMSE | vRMSE
---|---|---|---|---|---|---|---|---|---|---|---
1 | 1 | 64 | 0.3 | 0.0005 | 20 | 2 | causal | causal | 8 | 0.0770 | 0.0553
2 | 3 | 64 | 0.1 | 0.0005 | 20 | 2 | no | same | 12 | 0.175 | 0.165
3 | 1 | 32 | 0.2 | 0.001 | 20 | 2 | same | same | 8 | 0.148 | 0.132
Figure 8: Trajectory of training and validation losses, as a function of the
training epoch, for product 1. Note that the cost functions are calculated on
rescaled data, therefore they cannot be directly compared to the values
appearing on the first line of Table 9.

###  3.4 NN forecasting

We then used the models specified by the hyperparameter choices in Tables 8
and 9 to forecast price variations on the test set. We proceeded by training
the each selected model on the entire training _and_ validation set, for a
number of epochs corresponding to the epoch at which the training had stopped
during the grid search. Table 10 reports the test set performances of NN
models obtained in this way. Although they played no role in model training
and validation, here we also report the values of the MAE and MAPE metrics
obtained on the test set. We already observe that class B models always
outperform class A models in forecasting prices for all the three products.
Figure 9 displays the forecasts made by the best NN model and compares it to
the actual price variations, in the case of product 3.

Table 10: Performance of the fine-tuned NN models on the test dataset. NN
results
---
Product | 1 | 2 | 3
class | A | B | A | B | A | B
RMSE | 0.0617 | 0.0613 | 0.179 | 0.162 | 0.219 | 0.200
MAE | 0.0498 | 0.0511 | 0.135 | 0.126 | 0.157 | 0.150
MAPE | 0.0299 | 0.0305 | 0.0181 | 0.0168 | 0.0221 | 0.0212
Figure 9: Comparison of NN forecasts and observed data for price variations in
the test set, for product 3. The corresponding RMSE is 0.200.

###  3.5 Result comparison

We can now compare the results obtained with the various models that have been
tried. The comparison can be made by looking at Tables 5, 7 and 10

We begin by noting that Prophet performances are considerably poorer than the
ARIMA benchmark and neural networks. On one hand, this could be expected as
Prophet is based on a model with relatively low complexity that prioritizes
ease of use and of tuning. On the other, as mentioned in Section 3.2.2,
Prophet was progressively fitted also on the test set, effectively using more
data than the other models for the fit, so its performances have arguably
benefited from that. Nevertheless, one point that needs to be stressed is that
we were able to provide Prophet with the entire weekly price time series
comprising all of 2013 data, with no need for particular data pre-processing.
Therefore, although Prophet’s performances were not brilliant for the problem
studied in this paper, it could still be useful in certain contexts where
quick, preliminary forecasts are needed.

Turning to the other models, we observe that those featuring both CNN and LSTM
layers (class B) yielded the most accurate forecasts for all the three
products. At the same time, tuning them required by far the longest time
(approximately 20 hrs for each product using a NVIDIA RTX 3060 graphic card)
with moderate improvement on purely LSTM (class A) models: considering the
RMSE metric, the improvement amounts to about 9% for both product 2 and 3,
while for product 1 the two classes yield similar results (class B has lower
RMSE but higher MAE and MAPE). On the other hand, class A models required a
considerably lower computational effort (of the order of 30 minutes to select
the model hyperparameters for each product). Notice that once the grid search
was concluded, training the best model on the entire training plus validation
dataset required only a few minutes, both for class A and B models.

ARIMA models performed well if we take into account both the achieved values
of the metrics and time necessary for tuning. They were less accurate than
class B models, yielding RMSEs that were 23% higher RMSE for product 1, 7%
higher for product 2, and 8% higher for product 3444Similar considerations
apply if we instead look at MAE and MAPE metrics.. However, they required
about the same tuning time as class A models, and performed better for product
2 and 3: the RMSE obtained by ARIMA models was 23% higher for product 1, but
2% lower for products 2 and 3. We remark that ARIMA is univariate, while a
multivariate input dataset was used to train and test deep learning models:
this highlights the effectiveness of the ARIMA approach for the problem under
study, while at the same time suggesting that the additional data fields used
to perform the multivariate analysis were not so informative for price
prediction. We add that machine learning models were tuned over a larger
parameter space than the others: the search grids were made of 54 and
8978555Factoring out the three different choices of the timesteps number $n$,
which are not really part of the NN models but refer to the way data is
prepared, we are left with 2916 configurations for class B models.
hyperparameter configurations for class A and class B NN models respectively,
versus a maximum of 5 ARIMA configurations and the 25 of Prophet.

Another important aspect to consider in comparing the models is that dataset
size strongly affects the performance of machine learning models. To this
regard, we note that the sizes of the datasets were not very large (just over
200 for each product, as seen in Table 2), hence one could expect especially
the NN models performance to further improve when more historical data is made
available.

##  4 Discussion

In this paper, we have discussed the application of different methods to the
forecast of wholesale prices. We put a standard econometric model (ARIMA) side
by side with two different approaches to time series forecasting. These were
rather diverse both in kind and complexity, going from a simple additive model
using a piecewise linear trend and Fourier series seasonality (Prophet) to
deep learning models featuring both convolutional and LSTM layers. The
findings showed that while Prophet was quick to set-up and tune, requiring no
data pre-processing, it was not able to come close to the performance of the
other, well-established time-series forecasting models. Instead, we found that
the best deep learning models performed better than ARIMA, but also required
much longer times for the hyperparameter tuning.

The work done in this paper can be extended in many directions. First, it
would be interesting to carry out a similar analysis also for the forecasting
of sales, and to consider data with daily frequency instead of weekly. Sales
forecasting with higher frequency can indeed be relevant for wholesalers and
retailers. Second, a more refined version of the study would distinguish
between different customers, as the selling strategies adopted by the
wholesaler do certainly differ when dealing with customers of various kinds
and sizes. Customer profiling is an extremely interesting and active avenue of
applied research, which can clearly enhance companies’ returns. Therefore, in
the near future we plan to carry on the analysis by combining customer
classification algorithms with time-series forecasting: understanding how and
when each customer buys, how specific products are treated differently by
different customers (identifying _complementary_ and _substitute_ goods for
each of them), as well as relating the price elasticity of demand for
different product/customers, are all aspects that could lead to economic
benefits for the wholesaler. To this end, one would want to compare further
machine learning models to dynamic panel data where prices for each product
and customer and price demand elasticity are considered. The approach could
lead to an essential gain in the accuracy of the forecast, and would be an
important contribution to sales prediction analysis, which is becoming an
increasingly important part of modern business intelligence [39, 40]. Another
aspect of sure relevance would be to evaluate algorithms and models on data
including the effects of the COVID-19 outbreak, to both learn how the market
has transformed and help modifying selling strategies in adapting to the
current rapidly changing situation.

The limits encountered in the application of the forecasting tools examined in
this work encourage the evaluation of further models that could bring together
the advantages of each approach. Finally, it would be relevant to apply one of
the most exciting recent advances in machine learning, namely the attention
mechanism [41], to the forecasting task we have considered. Work on this
particular topic is definitely attractive and recently made its appearance in
this field [42].

##### Note:

The paper is an output of the project “Cancelloni Big Data Analytics (BDA)” at
Cancelloni Food Service S.p.A. Data used in this work are property of
Cancelloni Food Service S.p.A. and cannot be disclosed.

## References

  * [1] Ratnadip Adhikari and R.. Agrawal  “An Introductory Study on Time Series Modeling and Forecasting”  In _CoRR_ abs/1302.6613, 2013  arXiv: http://arxiv.org/abs/1302.6613
  * [2] Hassan Ismail Fawaz et al.  “Deep learning for time series classification: a review”  In _Data mining and knowledge discovery_ 33.4 Springer, 2019, pp. 917–963  DOI: https://doi.org/10.1007/s10618-019-00619-1
  * [3] Omer Berat Sezer, Mehmet Ugur Gudelek and Ahmet Murat Ozbayoglu  “Financial time series forecasting with deep learning : A systematic literature review: 2005–2019”  In _Applied Soft Computing_ 90, 2020, pp. 106181  DOI: https://doi.org/10.1016/j.asoc.2020.106181
  * [4] Bryan Lim and Stefan Zohren  “Time-series forecasting with deep learning: a survey”  In _Philosophical Transactions of the Royal Society A_ 379.2194 The Royal Society Publishing, 2021, pp. 20200209  DOI: http://doi.org/10.1098/rsta.2020.0209
  * [5] Pedro Lara-Benítez, Manuel Carranza-García and José C Riquelme  “An Experimental Review on Deep Learning Architectures for Time Series Forecasting”  In _International Journal of Neural Systems_ 31.03 World Scientific, 2021, pp. 2130001  DOI: https://doi.org/10.1142/S0129065721300011
  * [6] Yann LeCun, Yoshua Bengio and Geoffrey Hinton  “Deep learning”  In _Nature_ 521.7553 Nature Publishing Group, 2015, pp. 436–444  DOI: https://doi.org/10.1038/nature14539
  * [7] Sepp Hochreiter and Jürgen Schmidhuber  “Long Short-Term Memory”  In _Neural Computation_ 9.8, 1997, pp. 1735–1780  DOI: 10.1162/neco.1997.9.8.1735
  * [8] Yann LeCun and Yoshua Bengio  “Convolutional networks for images, speech, and time series”  In _The handbook of brain theory and neural networks_ 3361.10, 1995
  * [9] Emir Alibašić, Boris Fažo and Ivica Petrović  “A new approach to calculating electrical energy losses on power lines with a new improved three-mode method”  In _Tehnički vjesnik_ 26.2 Strojarski fakultet u Slavonskom Brodu; Fakultet elektrotehnike, računarstva …, 2019, pp. 405–411
  * [10] Zou Haofei, Xia Guoping, Yang Fangting and Yang Han  “A neural network model based on the multi-stage optimization approach for short-term food price forecasting in China”  In _Expert Systems with Applications_ 33.2 Elsevier, 2007, pp. 347–356  DOI: https://doi.org/10.1016/j.eswa.2006.05.021
  * [11] HF Zou, GP Xia, FT Yang and HY Wang  “An investigation and comparison of artificial neural network and time series models for Chinese food grain price forecasting”  In _Neurocomputing_ 70.16-18 Elsevier, 2007, pp. 2913–2923  DOI: https://doi.org/10.1016/j.neucom.2007.01.009
  * [12] Girish K Jha and Kanchan Sinha  “Agricultural price forecasting using neural network model: An innovative information delivery system”  In _Agricultural Economics Research Review_ 26.347-2016-17087, 2013, pp. 229–239  DOI: 10.22004/ag.econ.162150
  * [13] Hildegart Ahumada and Magdalena Cornejo  “Forecasting food prices: The case of corn, soybeans and wheat”  In _International Journal of Forecasting_ 32.3 Elsevier, 2016, pp. 838–848  DOI: https://doi.org/10.1016/j.ijforecast.2016.01.002
  * [14] Bohdan M. Pavlyshenko  “Machine-Learning Models for Sales Time Series Forecasting”  In _Data_ 4.1, 2019  DOI: 10.3390/data4010015
  * [15] Ioannis E Livieris, Emmanuel Pintelas and Panagiotis Pintelas  “A CNN–LSTM model for gold price time-series forecasting”  In _Neural computing and applications_ 32.23 Springer, 2020, pp. 17351–17360  DOI: https://doi.org/10.1007/s00521-020-04867-x
  * [16] Li Zhang et al.  “Prediction of stock prices based on LM-BP neural network and the estimation of overfitting point by RDCI”  In _Neural Computing and Applications_ 30.5 Springer, 2018, pp. 1425–1444  DOI: 10.1007/s00521-017-3296-x
  * [17] Gour Sundar Mitra Thakur, Rupak Bhattacharyya and Seema Sarkar Mondal  “Artificial Neural Network Based Model for Forecasting of Inflation in India”  In _Fuzzy Information and Engineering_ 8.1, 2016, pp. 87–100  DOI: https://doi.org/10.1016/j.fiae.2016.03.005
  * [18] Livia Paranhos  “Predicting Inflation with Neural Networks”, 2021  arXiv:2104.03757 [econ.EM]
  * [19] Gustavo Silva Araujo and Wagner Piazza Gaglianone  “Machine learning methods for inflation forecasting in Brazil: new contenders versus classical models”, 2020
  * [20] Xue Yan, Wang Weihan and Miao Chang  “Research on financial assets transaction prediction model based on LSTM neural network”  In _Neural Computing and Applications_ 33.1 Springer, 2021, pp. 257–270  DOI: 10.1007/s00521-020-04992-7
  * [21] Firuz Kamalov  “Forecasting significant stock price changes using neural networks”  In _Neural Computing and Applications_ 32.23 Springer, 2020, pp. 17655–17667  DOI: 10.1007/s00521-020-04942-3
  * [22] Yaping Hao and Qiang Gao  “Predicting the trend of stock market index using the hybrid neural network based on multiple time scale feature learning”  In _Applied Sciences_ 10.11 Multidisciplinary Digital Publishing Institute, 2020, pp. 3961  DOI: https://doi.org/10.3390/app10113961
  * [23] Yi Xiao, Jin Xiao and Shouyang Wang  “A hybrid model for time series forecasting”  In _Human Systems Management_ 31.2 IOS Press, 2012, pp. 133–143
  * [24] Jiayu Qiu, Bin Wang and Changjun Zhou  “Forecasting stock prices with long-short term memory neural network based on attention mechanism”  In _PLOS ONE_ 15.1 Public Library of Science, 2020, pp. 1–15  DOI: 10.1371/journal.pone.0227222
  * [25] Vinay Kumar Reddy Chimmula and Lei Zhang  “Time series forecasting of COVID-19 transmission in Canada using LSTM networks”  In _Chaos, Solitons & Fractals_ 135, 2020, pp. 109864  DOI: https://doi.org/10.1016/j.chaos.2020.109864
  * [26] Neo Wu, Bradley Green, Xue Ben and Shawn O’Banion  “Deep Transformer Models for Time Series Forecasting: The Influenza Prevalence Case”  In _CoRR_ abs/2001.08317, 2020  URL: https://arxiv.org/abs/2001.08317
  * [27] George EP Box, Gwilym M Jenkins, Gregory C Reinsel and Greta M Ljung  “Time series analysis: forecasting and control”  John Wiley & Sons, 2015
  * [28] Trevor Hastie and Robert Tibshirani  “Generalized Additive Models”  In _Statistical Science_ 1.3 Institute of Mathematical Statistics, 1986, pp. 297–310  DOI: 10.1214/ss/1177013604
  * [29] Wint Nyein Chan  “Time Series Data Mining: Comparative Study of ARIMA and Prophet Methods for Forecasting Closing Prices of Myanmar Stock Exchange”  In _J. Comput. Appl. Res._ 1, 2020, pp. 75–80
  * [30] Işil Yenidoğan et al.  “Bitcoin forecasting using ARIMA and prophet”  In _2018 3rd International Conference on Computer Science and Engineering (UBMK)_ , 2018, pp. 621–624  IEEE  DOI: 10.1109/UBMK.2018.8566476
  * [31] Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani  “An introduction to statistical learning”  Springer, 2013  DOI: DOI https://doi.org/10.1007/978-1-4614-7138-7
  * [32] Sean J. Taylor and Benjamin Letham  “Forecasting at Scale”  In _The American Statistician_ 72.1 Taylor & Francis, 2018, pp. 37–45  DOI: 10.1080/00031305.2017.1380080
  * [33] Richard H. Byrd, Peihuang Lu, Jorge Nocedal and Ciyou Zhu  “A Limited Memory Algorithm for Bound Constrained Optimization”  In _SIAM Journal on Scientific Computing_ 16.5, 1995, pp. 1190–1208  DOI: 10.1137/0916069
  * [34] Razvan Pascanu, Tomas Mikolov and Yoshua Bengio  “On the difficulty of training recurrent neural networks”  In _Proceedings of the 30th International Conference on Machine Learning_ 28.3, Proceedings of Machine Learning Research  Atlanta, Georgia, USA: PMLR, 2013, pp. 1310–1318  arXiv: http://proceedings.mlr.press/v28/pascanu13.html
  * [35] George EP Box and David A Pierce  “Distribution of residual autocorrelations in autoregressive-integrated moving average time series models”  In _Journal of the American statistical Association_ 65.332 Taylor & Francis, 1970, pp. 1509–1526  DOI: 10.2307/2284333
  * [36] Greta M Ljung and George EP Box  “On a measure of lack of fit in time series models”  In _Biometrika_ 65.2 Oxford University Press, 1978, pp. 297–303  URL: http://www.jstor.org/stable/2335207
  * [37] Diederik P Kingma and Jimmy Ba  “Adam: A method for stochastic optimization”  In _3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings_ , 2014  URL: http://arxiv.org/abs/1412.6980
  * [38] Xavier Glorot and Yoshua Bengio  “Understanding the difficulty of training deep feedforward neural networks”  In _Proceedings of the thirteenth international conference on artificial intelligence and statistics_ 9, Proceedings of Machine Learning Research, 2010, pp. 249–256  JMLR WorkshopConference Proceedings  URL: http://proceedings.mlr.press/v9/glorot10a.html
  * [39] John T Mentzer and Mark A Moon  “Sales forecasting management: a demand management approach”  Sage Publications, 2004
  * [40] G Peter Zhang  “Neural networks in business forecasting”  IGI global, 2004
  * [41] Ashish Vaswani et al.  “Attention is all you need”  In _Advances in neural information processing systems_ 30, 2017, pp. 5998–6008  arXiv: https://proceedings.neurips.cc/paper/2017/file/3f5ee243547dee91fbd053c1c4a845aa-Paper.pdf
  * [42] Vijay Ekambaram et al.  “Attention based Multi-Modal New Product Sales Time-series Forecasting”  In _Proceedings of the 26th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_ Association for Computing Machinery, 2020, pp. 3110–3118  DOI: https://doi.org/10.1145/3394486.3403362