A fuel sales forecast method based on variational Bayesian structural time series

Abstract

Fuel prices, which are of broad concern to the general public, are always seen as a challenging research topic. This paper proposes a variational Bayesian structural time-series model (STM) to effectively process complex fuel sales data online and provide real-time forecasting of fuel sales. While a traditional STM normally uses a probability model and the Markov chain Monte Carlo (MCMC) method to process change points, using the MCMC method to train the online model can be difficult given a relatively heavy computing load and time consumption. We thus consider the variational Bayesian STM, which uses variational Bayesian inference to make a reliable judgment of the trend change points without relying on artificial prior information, for our prediction method. With the inferences being driven by the data, our model passes the quantitative uncertainties to the forecast stage of the time series, which improves the robustness and reliability of the model. After conducting several experiments by using a self-collected dataset, we show that compared with a traditional STM, the proposed model has significantly shorter computing times for approximate forecast precision. Moreover, our model improves the forecast efficiency for fuel sales and the synergy of the distributed forecast platform based on an architecture of network.

Keywords

Change point recognition hierarchical Bayesian model structural time-series model variational inference

1. Introduction

Due to the gradual opening of China’s fuel market, the sales model of the traditional regional market is no longer effective. Fuel retail terminals are facing increasingly fierce competition. To improve competitiveness in the fuel market, it is necessary for these terminals to promptly formulate effective sales strategies. Traditional sales strategies depend on the judgments of terminal sellers regarding the market situation to a large extent and are highly subjective and uncertain. Currently, there are many widely distributed fuel terminal stations, which leads to complex fuel sales data and makes it difficult to conduct a real-time analysis. To address these issues, we design an efficient distributed platform of network architecture, to support real-time data acquisition and analysis, those fuel terminal stations’ large-scale data intensive upload and distributed storage, distributed real-time model training and predict, the network architecture includes 5G mobile network, VPN network based on internet, local intranet. Therefore, the focus of this research is learning how to develop an efficient online forecast using complex fuel sales data and provide support for the timely adjustment of sales strategies, based on this platform. The fuel sales data used in this study are a typical economic time series. An important part of the analysis of the time series is the selection of a suitable probability model (or class of models) for the data.

In the past, traditional statistical models, such as autoregressive (AR) models, autoregressive integrated moving average (ARIMA) models and generalized autoregressive conditional heteroscedastic (GARCH) models, were generally adopted to forecast fuel market information. AR models are most widely studied because of their flexibility to model many stationary processes. However, AR coefficients change over time and when an input variable switches smoothly between two regimes. In practice, many parameters have to be optimized to build a consistent time-varying smooth transition autoregressive (STAR) model, and no efficient analytical optimization method is available. ARIMA models are used in many areas of science. Conejo et al. [5] applied this type of model for electricity price forecasting, with good results. ARIMA models are largely limited to capturing first-order nonstationarity and second-order (moment) nonstationarity, i.e., time-varying conditional variance or volatility. The subsequently developed GARCH model is used when variances in the time series differ. Garcia et al. [9] dealt with the heteroscedasticity problem with electricity price time-series data by using a GARCH model.

These linear models provide only a coarse approximation to real-world complex systems and generally fail to accurately predict the evolution of nonlinear and nonstationary processes.

Neural networks (NNs) and other nonlinear operator-based models have been used for nonlinear time-series forecasting in manufacturing systems and finance. For oil price prediction research, Kristjanpoller and Minutolo proposed a novel hybrid model, comprising an artificial neural network and a GARCH model, to analyze and predict oil price return volatility, where the artificial neural network better improves accuracy compared to several traditional oil price prediction models, such as the GARCH and ARFIMA models [19]. Mahdiani and Khamehchi [22] developed an oil price prediction method based on a neural network model using the genetic algorithm to optimize the model’s parameters, and the results showed that the model has significantly improved performance, especially when there is a small amount of training data or the variables show high fluctuations. Baruník and Malinská [1] used a generalized regression framework-based method in addressing oil price forecasting problems, while neural networks are also being used in such frameworks. By using this method with a real-time dataset, results show that the new method consistently has the smallest error [1]. Gupta and Nigam [12] proposed another novel method of predicting crude oil prices using the artificial neural network (ANN) that delivers significantly better results than other models by modifying time delay changes on the closed optimized results, which are then evaluated by the root mean square error. However, these models are relatively simple, consider fewer market factors and involve strict constraints. Therefore, it is impossible for them to accurately and effectively reflect the complexity of fuel sales data.

To address the abovementioned problem, a structural time-series model (STM)-based market forecast method was proposed by Harvey and Fernandes [14] and Harvey and Shephard [15] as continuous-time structural models for use in forecasting lead time demand for inventory control. The focus of STM is not to accurately represent the data-generation process but to decompose the series into different components and recreate the series. Compared with a traditional time-series model, an STM appears to be more accurate, more explicit and flexible in terms of describing economic variables. In particular, it is much easier for an STM to accurately replicate the characteristics of changes in the series. Therefore, using an STM, time series may be decomposed into components with multiple factors, and different methods can be used to process different factors [13]. The advantage of this model lies in the fact that it can process the data better than other models if there are missing values in the sample data or the data involve more complex changes such as trend changes and unexpected changes. Due to the flexibility of the STM, it has been used in research on sustainable development: Brodersen et al. [2] of Google captured the trend, seasonal and similar components of the target series by developing a Bayesian structural time-series (BSTS) model, and other studies [3,6,31,32,36] demonstrated how the BSTS model could be applied given different events. Recently, Taylor and Letham [34] of Facebook developed a state decomposition model with adjustable parameters. When forecasting fuel sales, data are often missing due to differences in how fuel enterprises develop their information systems. Therefore, strict requirements are proposed regarding the change-point processing ability of the forecast model in terms of obtaining a good forecast result with limited fuel sales data. A traditional STM uses the Markov chain Monte Carlo (MCMC) method to establish a change-point recognition model. The MCMC method is more accurate but time-consuming. It is impossible for such a method to allow the front-end and the back-end to collaborate for real-time forecasting of fuel sales.

To efficiently forecast fuel sales in real time, this paper proposes an STM based on variational Bayesian inference. The major contribution of this paper is that it establishes a model that identifies the change points in data by using variational Bayesian inference. Through experimental verification, this study shows that this method is more efficient than the traditional MCMC method. For approximate forecast precision, the modeling time of the former is reduced to approximately 1/200 of that of the latter. In theory, the MCMC method can be used to find the optimal solution with a sufficient number of iterations. However, in practice, the results obtained by using the MCMC method are problematic due to an excessively large variance, which greatly decreases the efficiency of convergence and prevents the immediate ending of MCMC convergence. By overcoming the shortcomings of the MCMC method, the variational method addresses the problem and can be used for optimization; this model uses the gradient descent algorithm for fast iterations, thus making real-time fast iterations possible for the time-series forecast model. Therefore, the variational method provides an efficient and feasible means for the extensive deployment and implementation of a collaborative front-end and back-end forecast of fuel sales data.

This paper contributes in the following two ways:

A variational BSTS model that uses variational Bayesian inference to make a reliable judgment on the trend change points without relying on artificial prior information is proposed for forecasting.

To improve the robustness and reliability of the model, data-driven inferences are implemented, and the quantitative uncertainties are passed to the forecast stage of the time series.

This study proves through experimental verification that compared with a traditional STM, the proposed model has a significantly shorter computing time for approximate forecast precision. The model can improve the efficiency of forecasting fuel sales and the synergy of the forecast platform.

2. Structural time-series forecast model for fuel sales

The STM was first proposed by Harvey and Fernandes [14]. An STM can approximate the state space and address unobserved aspects of the data, including trends, cycles, seasonality and irregularities. As these factors are entered into the model as unobservable variables, it is impossible to estimate the model parameters by using the traditional regression analysis method. Therefore, these factors are expressed in the form of a state space in statistical processing. An STM can better interpret the characteristics of economic variables and can deal with complex time series that occur in real business activities. This paper forecasts complex fuel sales information by using an STM. In consideration of the economic indicators of fuel sales, the established model includes three major components, namely, trends, seasonality and promotion shocks. The STM model can be decomposed as follows (1): $\begin{matrix} (1) & y (t) = g (t) + s (t) + h (t) + ε_{t} \end{matrix}$ where $y (t)$ represents the observed value of the time series at time t; $g (t)$ represents the trend function, which models the nonperiodic changes in the times series value; $s (t)$ represents periodic seasonal changes; $h (t)$ represents the pulse response at the time of promotion; and $ε_{t}$ represents random noise in the time series, which is set to show a normal distribution. Thus, the STM for fuel sales forecasting is composed of a trend growth model, seasonable cycle model, promotion time pulse model and random noise model.

2.1. Establishment of the trend growth model

The trend growth factor reflects real changes in the time series and is the most basic factor in the time series. The instability of the time series is caused by the instability of the trend growth factor. Taylor and Letham [34] claimed that the trend model is mainly a combination of two trend models, the saturated growth model and the piecewise linear model. The development of the market and human society follow unique laws. Due to the impacts of resources and the environment, it is impossible for the development the market and human society to continue indefinitely. When development reaches its limit for a certain stage, saturation will undoubtedly occur. Hutchinson [16] claimed that market coverage growth will reach saturation. As the sales volume of the fuel market is constrained by market demand and production capacity, the trend model must be a combination of two trend models, namely, a saturated growth model and piecewise linear model. Because the corresponding professional knowledge and experience are required to establish a saturation model for the fuel market, the saturated growth model is omitted, and only a piecewise linear model is incorporated here. The specific form of the piecewise linear model is expressed as follows (2): $\begin{matrix} (2) & g (t) = (k + a {(t)}^{T} δ) t + (m + a {(t)}^{T} γ) \end{matrix}$ where m is the time offset parameter, k is the base growth rate, and δ is a vector of the parameters of the amplitude of the adjustments in the growth rate. Here, $γ_{j} = - s_{j} δ_{j}$ , which represents the time when change points occur. The trend growth model is composed of a series of trend line segments described by k and a series of change points. These points are defined as an S-dimensional vector. S is the number of change points. $s_{j}$ represents the time when change points occur, where $j = 1, \dots, S$ . The variation in the growth rate is determined by $δ_{j}$ . The growth rate is adjusted at time t following each change-point time $s_{j}$ . The growth rate of the time is the base growth rate plus the variation in the growth rates for all previous change points and can be expressed as follows (3): $\begin{matrix} (3) & k + \sum_{j : t > s_{j}} δ_{j} \end{matrix}$

Here, $a (t) \in 0, 1^{S}$ satisfies the following conditions, as shown in the following (4): $\begin{matrix} (4) & a_{j} (t) = \{\begin{matrix} 1 & t ⩾ s_{j}, \\ 0 & otherwise . \end{matrix} \end{matrix}$

A vector created by (3) and (4) is used to define the total variation of each change point. The growth rate at time t can be expressed as follows (5): $\begin{matrix} (5) & k + a {(t)}^{T} δ \end{matrix}$

Thus, all states at time t can be described by the following matrix (6): $\begin{array}{rcl} A & = & [\begin{matrix} t_{1} ⩾ s_{1} & t_{1} ⩾ s_{2} & \dots & t_{1} ⩾ s_{n} \\ t_{2} ⩾ s_{1} & t_{2} ⩾ s_{2} & \dots & t_{2} ⩾ s_{n} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ t_{k} ⩾ s_{1} & t_{k} ⩾ s_{2} & \dots & t_{k} ⩾ s_{n} \end{matrix}] \\ = & [\begin{matrix} 1 ⩾ 2 & 1 ⩾ 5 & \dots & 1 ⩾ n \\ 2 ⩾ 2 & 2 ⩾ 5 & \dots & 2 ⩾ n \\ ⋮ & ⋮ & ⋱ & ⋮ \\ k ⩾ 2 & k ⩾ 5 & \dots & k ⩾ n \end{matrix}] \\ (6) & = & [\begin{matrix} 0 & 0 & \dots & 0 \\ 1 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 1 & 1 & \dots & 1 \end{matrix}] \end{array}$

Equations (2), (4) and (5) can be combined to obtain the trend growth model for the STM for forecasting fuel sales according to (7): $\begin{matrix} (7) & \vec{g} = (k + A \vec{δ}) ⊙ \vec{t} + (m + A \vec{γ}) \end{matrix}$

2.2. Establishment of the seasonal cycle model

Due to the impacts of mankind’s commercial behaviors, fuel sales fluctuations tend to have a seasonal cycle, with multiple periods. To fit and forecast this seasonality, the cycle function of $s (t)$ is established as a seasonal cycle model. The Fourier series are employed in this paper for the cycle effect model. Here, P is the regular cycle for the expected time series (the annual data P is set as 365.25 and the weekly data as 7). The smooth seasonal effect can be approximated by the following (8): $\begin{matrix} (8) & s (t) = \sum_{n = 1}^{N} (a_{n} cos (\frac{2 π n t}{P}) + b_{n} sin (\frac{2 π n t}{P})) \end{matrix}$

When fitting the seasonal cycle factor using Fourier series, $2 N$ parameters $β = [a_{1}, b_{1}, \dots, a_{n}, b_{n}]$ must be estimated. The least squares method is used to solve for the parameters in (8) according to (9) and (10): $\begin{array}{l} (9) & a_{n} = \frac{2}{P} \sum_{i = 1}^{n} y_{i} cos \frac{2 π n t}{P}, n = 1, 2, \dots, m \\ (10) & b_{n} = \frac{2}{P} \sum_{i = 1}^{n} y_{i} sin \frac{2 π n t}{P}, n = 1, 2, \dots, m \end{array}$ where m is the largest integer and is not larger than $n / 2$ . The seasonal components with large impacts can be obtained by using (9) and (10), yielding (11): $\begin{matrix} (11) & G (t) = \sum_{n = 1}^{m} \sqrt{a_{n}^{2} + b_{n}^{2}} sin (\frac{2 n π t}{P} + φ) \end{matrix}$ where $tan (φ) = a_{n} / b_{n}$ . When there is an n that will allow for a large value of $\sqrt{a_{n}^{2} + b_{n}^{2}}$ , then a seasonal component exists. A seasonal vector regression matrix is constructed to represent each value of t in previous and future data. For annual seasonality, $N = 30$ . The seasonable vector regression matrix is expressed as follows (12): $\begin{matrix} (12) & R (t) = [cos (\frac{2 π (1) t}{365.25}), \dots, sin (\frac{2 π (30) t}{365.25})] \end{matrix}$

By using (12), the main component of the seasonal factor of the fuel sales forecast STM can be expressed as (13): $\begin{matrix} (13) & s (t) = R (t) β \end{matrix}$ where β obeys a normal distribution, and the value indicating the prior season’s smoothness can be calculated. The larger N is, the better the fit to the seasonal factor of the time series. However, there is an increased risk of model overfitting [7,24,30].

2.3. Impulse response function at the time of promotion

Due to the characteristics of fuel sales, fuel sales promotion and other behaviors will impact the sales time series to a large extent, and these shocks are periodic. Therefore, there is no way to model the impacts of these factors using a stale cycle function. Only the impacts of fuel promotion behaviors are considered in this paper. For each promotion period i, $D_{i}$ represents the set of dates for the promotion period. An indicator function is constructed to indicate whether time t occurs during a promotion, and one parameter $κ_{i}$ is allocated to each promotion period to forecast the corresponding changes. The regression matrix is established as follows (14): $\begin{matrix} (14) & W (t) = [1 (t \in D_{1}), \dots, 1 (t \in D_{i})] \end{matrix}$

To reflect the impact of the promotion on the model, a linear pulse response function is adopted. Assuming that the prior ρ obeys a normal distribution of $N (0, v^{2})$ , it is superimposed into (14), and the promotion pulse response function is obtained as follows (15): $\begin{matrix} (15) & h (t) = W (t) ρ \end{matrix}$

3. Variational Bayesian change-point testing model

To establish the fuel sales forecast STM, a fuel sales growth trend model is required to forecast fuel sales. The change point $s_{j}$ implied in (7) is used for online detection and represents the data point or time point at which the trend changes. Studies on change-point testing mainly follow either the probability or Bayesian school of thought. The former tends to test change points by using the methods of Choy inspection, Bayes method, Schwarz information criterion method, product partition model (PPM) [16], maximum likelihood, least squares or local comparisons. Due to the rapid development of big data forecasting in recent years, Bayesian theory is being widely applied. The main idea of the Bayesian inference method is to adopt parameters in the model with change points and regard the change points as random variables, then gradually derive the posterior distribution of the change points in the model by using prior distribution and sample distribution, and estimate other parameters according to the posterior distribution of the positions of change points. The advantages of the Bayesian inference method include the following: prior information is known in advance, the model inference is more accurate, and the posterior probability density function can be obtained for change-point positions, thus enabling a complete, comprehensive and profound understanding of the positions [8,33,37]. The key to Bayesian inference is to calculate the mean and variance of the posterior distribution, but it is difficult to directly obtain a sample of the posterior distribution through direct calculation when the posterior distribution is excessively large. Currently, the most widely used solution is the MCMC algorithm [20], which has low computing efficiency when using a large amount of data. Variational inference provides an approximate inference method that can be used for a large amount of data. Variational inference does not have a higher computing speed than the MCMC method, but it can obtain posterior inference results that are approximate to those obtained by using the MCMC method [4]. Additionally, as the inference results obtained by using this algorithm are based on the same distribution as those obtained by using the MCMC method, this algorithm can introduce uncertainty about the changes for parameter estimation. Additionally, in terms of efficiency, the algorithm can better facilitate the fast forecasting of massive stream data for business forecasts. The ability to introduce uncertainty about changes can improve the robustness of the model. The change points are identified through variational inference, which is completed in the following three steps:

The first step is establishing hierarchical prior information; the second step is inferring the posterior distribution using variational inference; and the third step is forecasting future data points using the posterior distribution.

3.1. Hierarchical Bayesian prior information modeling

In this paper, the change points can be automatically selected through the sparse prior δ. There are S change points among previous T data points. To make the change points sparse, the prior distribution of the growth rate and growth points is selected as $δ_{j} \sim Laplace (0, σ)$ , equivalent to the L1 regularization of the parameters. By parameterizing the scale, the change rate can be altered based on the flexibility of the control model. In fact, $δ_{j} \sim Laplace (0, σ)$ is a hierarchical prior where the percentage of change-point priors is the parameter. λ needs to be parameterized, and the posterior analytical solution needs to be solved using a hierarchical Bayesian model. It is assumed that the parameter obeys the distribution [21,28] given by the following (16): $\begin{matrix} (16) & p (θ, λ) = p (θ | λ) p (λ) \end{matrix}$

To prevent the parameter in (16) from being restricted by prior information, a vector containing a series of hyperparameters, $λ = [λ_{1}, λ_{2}, \dots, λ_{n}]$ , is introduced. This vector is used to minimize the adverse effects of subjective information. The marginal distribution of $p (θ)$ is thus obtained and expressed as follows (17): $\begin{matrix} (17) & p (θ) = \int p (θ, λ) d λ \end{matrix}$

To establish the hierarchical Bayesian model, it is necessary to parameterize each Bayesian layer. It is generally assumed that these parameters obey a multivariate normal distribution and are expressed as follows (18): $\begin{matrix} (18) & p (θ) = \prod_{i = 1}^{K} p (θ_{i} | μ_{θ_{i}}, τ_{θ_{i}}) = \prod_{i = 1}^{K} N (μ_{θ_{i}}, τ_{θ_{i}}^{- 1} I) \end{matrix}$

In the first layer, $μ_{θ} = 0$ and $τ_{θ}$ represents random variables that obey the gamma distribution. Thus, the hierarchical Bayesian model decomposes the prior distribution as follows (19): $\begin{matrix} (19) & \begin{matrix} p (θ | τ_{θ}, α_{τ}, β_{τ}) = p (θ | τ_{θ}) p (τ_{θ} | α_{τ}, β_{τ}) \\ α_{τ} > 0, β_{τ} > 0 \end{matrix} \end{matrix}$ where $α_{τ}$ and $β_{τ}$ represent the two parameters of the gamma distribution. From (16) and (17), the integral of $τ_{θ}$ in (19) can be removed, and a t-distribution is obtained [26]. The t-distribution allows a normal distribution to be approximated when the sample is small. The prior distribution of (19) can be changed to the following form (20): $\begin{array}{rcl} p (θ | α_{τ}, β_{τ}) & = & \int_{0}^{\infty} p (θ | τ_{θ}, α_{τ}, β_{τ}) d τ_{θ} \\ (20) & = & St (0, \frac{α_{τ}}{β_{τ}}, 2 α_{τ}) \end{array}$

The discussion above describes the method used to construct the hierarchical Bayesian model. Assume that the problem to be solved is a linear regression. Then, let $σ^{2}$ replace τ and β replace δ. The aim is to make the variables sparse by ensuring that the prior of the regression coefficient is a Laplace prior so the change points can be automatically selected. According to the Bayesian equation, the variation in the change points obeys the distribution [27,29] given by (21): $\begin{matrix} (21) & π (β | y) = π (y | β, σ^{2}) π (β | σ^{2}) π (σ^{2}) \end{matrix}$

Equation (22) is obtained by assuming the prior distribution of the regression coefficients obeys the Laplace prior: $\begin{matrix} (22) & π (β | σ^{2}) = \prod_{j = 1}^{p} \frac{λ}{2 \sqrt{σ^{2}}} exp {- \frac{λ | β_{j} |}{\sqrt{σ^{2}}}} \end{matrix}$

Expressing (22) as a mixed normal distribution with a mixed exponential distribution gives (23): $\begin{matrix} (23) & \frac{λ}{2 \sqrt{σ^{2}}} exp {- \frac{λ | β |}{\sqrt{σ^{2}}}} = \int_{0}^{\infty} \frac{1}{\sqrt{2 π σ^{2} τ^{2}}} exp {- \frac{β^{2}}{2 σ^{2} τ^{2}}} \frac{λ^{2}}{2} exp {- \frac{λ^{2}}{2} τ^{2}} d τ^{2} \end{matrix}$

The prior is set. For any given adjusting parameter, the posterior distribution $π (β, σ^{2} | y)$ is a single-peak posterior. The posterior distribution of the single peak can accelerate the convergence rate of sampling. The prior distribution of the parameter is a $σ^{2}$ exponential distribution or $π (σ^{2}) \sim exponential (λ)$ . The hierarchal Bayesian model is expressed from (24) to (28) [11]: $\begin{array}{l} (24) & y | X, β, σ^{2} \sim N (X β, σ^{2} I_{n}) \\ (25) & β | σ^{2}, τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2} \sim N (0_{p}, σ^{2} D_{τ}) \\ (26) & D_{τ} = diag (τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2}) \\ (27) & σ^{2}, τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2} \sim π (σ^{2}) \prod_{j = 1}^{p} \frac{λ^{2}}{2} exp {- \frac{λ^{2} τ_{j}^{2}}{2}} \\ (28) & σ^{2} \sim exponential (ρ) \end{array}$

Then, the posterior distribution of the unknown parameters is given by the following (29): $\begin{array}{l} π (β, σ^{2}, τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2} | y) \\ (29) & \propto π (y | β, σ^{2}) π (β | σ^{2}, τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2}) π (σ^{2}, τ_{1}^{2}, τ_{2}^{2}, \dots, τ_{p}^{2}) \end{array}$

The hierarchal Bayesian framework for the trend model capable of automatic change-point detection is shown in Fig. 1 below.

Fig. 1.

Hierarchal Bayesian framework for the trend model capable of automatic change-point detection. The value of sigma has a sparse prior distributed as a Laplace density with parameter tau.

3.2. Approximating the posterior variational inference model

After obtaining the hierarchical prior distribution of a probability model, either variational inference will be used to infer the posterior distribution, or the iterative prior is used to update the parameter. Under the Bayesian framework, the joint posterior distribution is expressed as follows [18]: $\begin{matrix} (30) & p (θ, λ | D) = \frac{p (D | θ, λ) p (θ | λ) p (λ)}{p (D)} \end{matrix}$

By removing the integral of (30), the conditional distribution for the marginal distribution can be further simplified as $p (θ | D)$ , and the normalization factor of the Bayesian equation, $p (D)$ , can be expressed as follows: $\begin{matrix} (31) & p (D) = \int_{θ} p (D | θ) p (θ) d θ \end{matrix}$

When there is a large parameter space but limited samples, it is not easy to obtain (31). However, variational inference can be used for processing, and ξ is parameterized to approximate the posterior distribution [23,25] as follows (32): $\begin{matrix} (32) & p (θ | D) \approx q_{ξ} (θ) \in Ξ \end{matrix}$

Variational inference tries to make $q_{ξ} (θ)$ approximate $p (θ | D)$ as close as possible through continuous iterations. Usually, Kullback–Leibler (KL) divergence is used to measure the distance between the two distributions. Thus, variational inference uses an optimization algorithm for approximation. Here, the target function of the evidence lower bound (ELBO) is expressed as follows (33): $\begin{array}{rcl} p (θ | D) & ≃ & q_{ξ *} (θ) \\ = & arg min K L (q_{ξ} (θ) | p (θ | D)) \\ (33) & = & arg min \int q_{ξ} (θ) log \frac{q_{ξ} (θ)}{p (θ | D)} d θ \end{array}$

To minimize the KL divergence, it is converted into the ELBO that maximizes the following (34): $\begin{array}{rcl} q_{ξ *} (θ) & = & L (θ) \\ = & arg max E_{q_{ξ}} (θ) [log p (θ, D) - log q_{ξ} (θ)] \\ (34) & = & arg max E_{q_{ξ}} (θ) [log p (D | θ)] - K L (q_{ξ} (θ) | p (θ)) \end{array}$ where the first conditional likelihood is employed as a data item and usually enables the posterior distribution to interpret data by maximizing the expected log likelihood value, providing unbiased random estimates for log likelihood, and classifying a data set into multiple subsets to calculate the differentiable posterior distribution. In addition to accelerating the calculation process, this model updates at a higher frequency, thus making the convergence more robust. This model can effectively avoid the risk of being trapped in the local minimum, as shown below (35): $\begin{array}{rcl} log p (D | θ) & = & \sum_{i = 1}^{N} log p (y_{i} | x_{i}, θ) \\ (35) & \approx & \frac{N}{M} \sum_{i = 1}^{M} log p (y_{i} | x_{i}, θ) \end{array}$

The mean field variational inference (MFVI) method is used to control the computing complexity of the second item in (35). The variational distribution is expressed by using a factor distribution subject to the decomposition layer by layer, where each factor is determined by its own variation parameters, as given by (36): $\begin{matrix} (36) & q_{ξ_{i}} (θ_{i}) = \frac{exp (\int log p (D, θ) \prod_{j \neq i} q_{ξ_{i}} (θ_{j}) d θ_{j})}{\int exp (\int log p (D, θ) \prod_{j \neq i} q_{ξ_{i}} (θ_{j}) d θ_{j}) d θ_{i}} \end{matrix}$

Equations (35) and (36) are substituted into (34). Then, the ELBO target function is transformed into the following (37): $\begin{array}{rcl} L (θ, ξ) & = & \int \frac{N}{M} \sum_{i = 1}^{M} log p (y_{i} | x_{i}, θ) \prod_{j = 1}^{K} q_{ξ_{j}} (θ_{j}) d θ_{j} \\ (37) & - \sum_{j = 1}^{K} \int q_{ξ_{i}} (θ_{j}) log q_{ξ_{j}} (θ_{j}) d θ_{j} \end{array}$

When (37) satisfies the convergence conditions, the iterations are terminated, and the final optimization result is obtained. The variational algorithm can be used for online detection of the change points in the fuel sales growth trend model. Then, the fuel sales forecast STM is established for the efficient online detection of fuel sales [10,17,35,38].

4. Experimental verification

4.1. Fuel sales data acquisition for testing and descriptions

To verify the efficiency and effect of the fuel sales forecast using the variational BSTS, an experiment is conducted with previous fuel sales data. The experiment uses the granularity data of a fuel sales company, which is the branch of the China National Petroleum Corporation (CNPC), from January 2015 to December 2018. The data for the last month or December 2018 (30 days) are set as the testing data. The preliminary data are training data, which are used to forecast the fuel sales of the company. The data include information on the time, fuel sales and relevant promotion time. As some data are missing, data need to be interpolated based on the experiences of experts in relevant fields. Fuel sales data can be acquired and analyzed by using the online fuel sales data acquisition platform shown in Fig. 2. The function modules of the management platform mainly include comprehensive management platforms, collaborative offices, mobile applications and other main modules. Among them, the comprehensive management platform is capable of marketing data asset management and analysis.

Fig. 2.

Online fuel sales data acquisition distributed platform.

4.2. Experiment results

The experiment conducted in this paper considers only the positions of the change points. When assuming the number of change points, one change point is designated at a certain time interval in advance, and then, the change rate is realized through the sparse prior based on the variation in the change points and through hierarchical Bayesian to reduce the variation to 0 and automatically select the change points. Here, the number of change points is set to 15, and the variation in the sparse change points is shown in Fig. 3. Through sparse variation, the change points automatically selected are obtained, as shown in Fig. 4 below. Figures 3 and 4 show that among the 15 change points that are designated on the time axis and uniformly distributed, eight are automatically selected through the regularization of the sparse prior. Figure 4 also shows that there are change points automatically selected by the algorithm near each inflection point, and the long-term trend change points are around the change points with large variations.

Fig. 3.

Sparse change-point variation. Fuel sales have suffered a persistent depression since the fourth quarter of 2018.

Fig. 4.

Locations of the sparse change points. The first plot shows the result of the whole trend function as a reasonable continuous trend. The middle plot shows the growth rate at time t. The last plot shows the offset at time t. These can adjust the growth rate’s height.

To verify the estimated efficiency of the fuel forecast method used with the variational BSTS for the variation, the MCMC method and variational inference are used to estimate the parameter posterior of the model.

Figure 5 shows the results when simulating 1,000 samples (2,000 iterations) by using the MCMC method for 39 minutes. Figure 6 shows the results when simulating 10,000 samples (2,000 iterations) by using the MCMC method for 480 minutes (8 h).

Fig. 5.

Simulation of 1,000 samples by using the Markov chain Monte Carlo (MCMC) method. (Left: kernel estimate posterior density plots; right: sampling trace plots). The chain exhibits very high autocorrelation according to the trace plot. An issue of high autocorrelation between posterior samples can be raised. Postprocessing algorithms require samples to be independent.

A comparison of Figs 5 and 6 shows that the distribution function curves of the latter appear much smoother. Additionally, the distribution of the sample values indicate that the latter has a random state and no significant correlation and trend. Therefore, the convergence effect is much better when simulating 10,000 samples using the MCMC method. However, using the MCMC method for 10,000 samples took 480 minutes. Despite its better convergence effect, its efficiency fails to satisfy the requirements of the existing massive stream data for the model to forecast online training in a real-time manner.

Fig. 6.

Simulation of 10,000 samples by using the Markov chain Monte Carlo (MCMC) method. (Left: kernel estimate posterior density plots, right: sampling trace plots). The chain is exploring the sample space well, and autocorrelation between posterior samples is no longer significant.

Variational inference is used in the experiment for parameter inference, and 10,000 samples are simulated, as shown in Fig. 7. In only 2 minutes and 40 seconds, results similar to those obtained using the MCMC method with 10,000 samples are obtained. The distribution curves in the figure show that the distribution shape obtained through variational inference is close to that obtained with its counterpart using the MCMC method.

Fig. 7.

Simulation of 10,000 samples by using variational inference. (Left: kernel estimate posterior density plots; right: sampling trace plots). As shown in Fig. 6, autocorrelation between the samples is not significant, and this figure shows similar parameter estimations.

Figure 8 shows that the variables estimated by using the variational method converge earlier, and the estimated variable mean is generally close to that obtained by using the MCMC method. Figures 9 and 10 show the sample distribution of the estimated variables obtained by using both methods. The covariances of the estimated variables of the two methods are basically similar. Obviously, the simulated sample distribution obtained by using the variational inference method is more concentrated, which increases the speed of attaining convergence. The results of the experiment show that the variational inference after processing the Bayesian posterior distribution not only obtains excellent results but also has obvious computing advantages when using large amounts of sample data.

Fig. 8.

Estimated variable means obtained by using the Markov chain Monte Carlo (MCMC) method (left) and variational inference (right). Y-axis: number of samples; X-axis: iterations.

Fig. 9.

Distribution of the estimated variable samples by using the Markov chain Monte Carlo (MCMC) method. Five parameters are shown (sigma, tau, m, k, and delta). The diagonal cells are the marginal distributions of the individual parameters; the other cells are the joint distribution of the two parameters. There are negative correlations between m and k and between m and delta. There is a positive correlation between k and delta.

Fig. 10.

Distribution of the estimated variable samples by using variational inference. Five parameters are shown (sigma, tau, m, k, and delta). Similar to Fig. 9, the joint distribution of m, k and delta exhibits a correlation, but its shape is narrower than the shape in Fig. 9, and this explains the quick convergence of variational inference.

The variational BSTS model and the traditional STM are now used to forecast the experimental data. The forecast results for three gas stations are shown in Figs 11 and 12. The results of the experiment are shown in Table 1. In terms of forecast precision, the performance of the variational inference structural time-series model (VI-STM) model is similar to that of the MCMC-STM model. By contrast, the computing time of the VI-STM model is significantly shorter than that of MCMC-STM and can ensure accurate online forecasting of complex fuel sales data.

Fig. 11.

Forecast by using the Markov chain Monte Carlo structural time-series model (MCMC-STM). The red line represents the actual values, and the blue line represents the forecasted values. The gray band around the blue line is called the “uncertainty interval”, and similar to the confidence interval, it indicates the range of forecasted values produced by the model.

Fig. 12.

Forecast by using the variational Bayesian structural time-series (BSTS) model. The red line represents the actual values, and the blue line represents the forecasted values. The gray band around the blue line is called the “uncertainty interval”, and similar to the confidence interval, it indicates the range of forecasted values produced by the model.

Table 1

Forecast precision of the different time-series models

Model	Gas Station	MAPE	MAE	RMSE	$R^{2}$
MCMC-STM	Station1	7.4%	578.43	743.22	0.56
	Station2	8.8%	1150.95	2552.43	0.55
	Station3	10.8%	633.93	2091.28	0.72
VI-STM	Station1	9.5%	688.70	1526.56	0.55
	Station2	12.2%	1452.23	2565.66	0.54
	Station3	14.16%	822.08	2092.72	0.70
ARIMA	Station1	24.3%	1897.35	2982.47	0.49
	Station2	28.8%	3796.28	4769.71	0.46
	Station3	30.0%	1786.53	3140.76	0.62
GARCH	Station1	20.2%	1733.74	3687.91	0.52
	Station2	29.6%	3901.97	5072.43	0.47
	Station3	26.9%	1601.61	2811.23	0.63
LSTM	Station1	16.6%	1296.13	2687.70	0.53
	Station2	18.6%	2451.59	3723.44	0.52
	Station3	18.4%	1095.23	2592.69	0.67

MAPE-Mean absolute percentage error, MAE-Mean absolute error, RMSE-Root mean squared error, MCMC-STM-Markov chain Monte Carlo structural time-series model, VI-STM- Variational inference structural time-series model, ARIMA-Autoregressive integrated moving average, GARCH-Generalized autoregressive conditional heteroscedastic, LSTM-Long and short-term memory.

To comprehensively evaluate the forecasting precision of the variational STM in forecasting fuel sales, the experimental data were assessed with the ARIMA time-series model, the GARCH time-series model and an LSTM model. In this section, we also set the last month of 2018 as the forecasting period, and the data from other periods were used for training.

Different $ARIMA (p, d, q)$ values may have an influence in forecasting for the same batch of data; by comparing the mean absolute percentage error (MAPE) of the model, we chose $ARIMA (1, 1, 1)$ , which yielded the lowest MAPE. Additionally, $GARCH (1, 1)$ was selected based on the same principle. LSTM models are deep learning methods employed for time-series prediction. Through experiments, we found that the accuracy and generalization ability of LSTM can be balanced by selecting the appropriate number of epochs. When we trained the LSTM model with 1,000 epochs, the forecasting performance of the model was optimal. The $R^{2}$ fitness values for the training data are shown in Table 1. Although the LSTM fitness is close to that of VI-STM, the MAPE of LSTM is obviously much worse than that of VI-STM. Moreover, the ARIMA and GARCH models displayed the worst performance in fitting and predicting complex fuel sales data. The results indicate that STM models are obviously better than other traditional statistical models and LSTM methods in predicting complex fuel sales data. The forecast precision of the STM model based on variational Bayesian inference is slightly lower than that of the STM model based on the MCMC method which is shown in Table 1. However, considering that the variational method greatly improves the computing efficiency, the gap between the two models in terms of forecast precision can be accepted.

5. Conclusions

To efficiently forecast fuel sales data online, a variational Bayesian structural time-series fuel forecast model is established considering the characteristics of the data, variational Bayesian inference theory and a structural time-series. The model is also used for the real-time forecasting of fuel sales time-series data. The model uses variational inference to process the change points. Compared with the traditional STM, this model greatly reduces the computing time. Furthermore, we design an efficient distributed platform of network architecture to support time-series fuel forecast model, more effectively forecasts fuel sales, has greater synergy with the forecast platform and satisfies the current demand for large-scale computing and online model training. Finally, it needs to be mentioned that variational Bayesian just get approximate results, through the experiment we found that VI-STM’s forecasting performance is unstable for different fuel terminal stations, particularly when time series includes much missing data and much outliers, it’s performance of prediction is poor, our next work will be focused on solving these problems.

References

Baruník and

Malinská, Forecasting the term structure of crude oil futures prices with neural networks, Applied Energy 164 (2016), 366–379. doi:10.1016/j.apenergy.2015.11.051.

K.H.

Brodersen,

Gallusser,

Koehler,

Remy and

S.L.

Scott, Inferring causal impact using Bayesian structural time-series models, Annals of Applied Statistics 9 (2015), 247–274. doi:10.1214/14-AOAS788.

Chai,

Lu,

Hu,

Wang,

K.K.

Lai and

Liu, Analysis and Bayes statistical probability inference of crude oil price change point, Technological Forecasting and Social Change 126 (2018), 271–283. doi:10.1016/j.techfore.2017.09.007.

Cheng,

S.T.S.

Bukkapatnam,

L.M.

Raff,

Hagan and

Komanduri, Monte Carlo simulation of carbon nanotube nucleation and growth using nonlinear dynamic predictions, Chemical Physics Letters 530 (2012), 81–85. doi:10.1016/j.cplett.2012.01.067.

A.J.

Conejo,

M.A.

Plazas,

Espinola and

A.B.

Molina, Day-ahead electricity price forecasting using the wavelet transform and ARIMA models, IEEE Transactions on Power Systems 20 (2005), 1035–1042. doi:10.1109/TPWRS.2005.846054.

J.G.

De Gooijer and

R.J.

Hyndman, 25 years of time series forecasting, International Journal of Forecasting 22 (2006), 443–473. doi:10.1016/j.ijforecast.2006.01.001.

D.F.

Findley,

K.C.

Wills and

B.C.

Monsell, Seasonal adjustment perspectives on “Damping seasonal factors: Shrinkage estimators for the X-12-ARIMA program”, International Journal of Forecasting 20 (2004), 551–556. doi:10.1016/j.ijforecast.2004.03.008.

Fryzlewicz and

Subba Rao, Multiple-change-point detection for auto-regressive conditional heteroscedastic processes, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 76 (2014), 903–924. doi:10.1111/rssb.12054.

R.C.

Garcia,

Contreras,

vanAkkeren and

J.B.C.

Garcia, A GARCH forecasting model to predict day-ahead electricity prices, IEEE Transactions on Power Systems 20 (2005), 867–874. doi:10.1109/TPWRS.2005.846044.

10.

E.S.

Gardner, Exponential smoothing: The state of the art, Journal of Forecasting 4 (1985), 1–28. doi:10.1002/for.3980040103.

11.

S.J.

Gershman and

D.M.

Blei, A tutorial on Bayesian nonparametric models, Journal of Mathematical Psychology 56 (2012), 1–12. doi:10.1016/j.jmp.2011.08.004.

12.

Gupta and

Nigam, Crude oil price prediction using artificial neural network, Procedia Computer Science 170 (2020), 642–647. doi:10.1016/j.procs.2020.03.136.

13.

Harvey, Forecasting with unobserved components time series models, in: Handbook of Economic Forecasting,

Elliott,

C.W.J.

Granger and

Timmermann, eds, Elsevier, New York, NY, 2006, pp. 327–412. doi:10.1016/S1574-0706(05)01007-4.

14.

A.C.

Harvey and

Fernandes, Time series models for count or qualitative observations, Journal of Business & Economic Statistics 7 (1989), 407–417.

15.

A.C.

Harvey and

Shephard, 10 structural time series models, Handbook of Statistics 11 (1993), 261–302. doi:10.1016/S0169-7161(05)80045-8.

16.

G.E.

Hutchinson, An Introduction to Population Ecology, Yale University Press, New Haven, CT, 1978.

17.

Ko and

Fox, GP-BayesFilters: Bayesian filtering using Gaussian process prediction and observation models, Autonomous Robots 27 (2009), 75–90. doi:10.1007/s10514-009-9119-x.

18.

Koutsoyiannis and

Montanari, Statistical analysis of hydroclimatic time series: Uncertainty and insights, Water Resources Research 43 (2007), W05429.

19.

Kristjanpoller and

M.C.

Minutolo, Forecasting volatility of oil price using an artificial neural network-GARCH model, Expert Systems with Applications 65 (2016), 233–241. doi:10.1016/j.eswa.2016.08.045.

20.

Lavielle and

Lebarbier, An application of MCMC methods for the multiple change-points problem, Signal Processing 81 (2001), 39–53. doi:10.1016/S0165-1684(00)00189-4.

21.

C.B.

Lee, Bayesian analysis of a change-point in exponential families with applications, Computational Statistics & Data Analysis 27 (1998), 195–208. doi:10.1016/S0167-9473(98)00009-7.

22.

M.R.

Mahdiani and

Khamehchi, A modified neural network model for predicting the crude oil price, Intellectual Economics 10 (2016), 71–77. doi:10.1016/j.intele.2017.02.001.

23.

R.M.

Neal, Bayesian Learning for Neural Networks, Lecture Notes in Statistics, Springer, New York, NY, 1996.

24.

Paap,

P.H.

Franses and

Hoek, Mean shifts, unit roots and forecasting seasonal time series, International Journal of Forecasting 13 (1997), 357–368. doi:10.1016/S0169-2070(97)00023-X.

25.

Park and

Casella, The Bayesian lasso, Journal of the American Statistical Association 103 (2012), 681–686. doi:10.1198/016214508000000337.

26.

Perreault,

Bernier,

Bobée and

Parent, Bayesian change-point analysis in hydrometeorological time series. Part 1. The normal model revisited, Journal of Hydrology 235 (2000), 221–241. doi:10.1016/S0022-1694(00)00270-5.

27.

Perreault,

Bernier,

Bobée and

Parent, Bayesian change-point analysis in hydrometeorological time series. Part 2. Comparison of change-point models and forecasting, Journal of Hydrology 235 (2000), 242–263. doi:10.1016/S0022-1694(00)00271-7.

28.

Perreault,

Haché,

Slivitzky and

Bobée, Detection of changes in precipitation and runoff over eastern Canada and U.S. using a Bayesian approach, Stochastic Environmental Research and Risk Assessment (SERRA) 13 (1999), 201–216. doi:10.1007/s004770050039.

29.

Perreault,

É.

Parent,

Bernier,

Bobée and

Slivitzky, Retrospective multivariate Bayesian change-point analysis: A simultaneous single change in the mean of several hydrological sequences, Stochastic Environmental Research and Risk Assessment 14 (2000), 243–261. doi:10.1007/s004770000051.

30.

Pfeffermann,

Morry and

Wong, Estimation of the variances of X-11 ARIMA seasonally adjusted estimators for a multiplicative decomposition and heteroscedastic variances, International Journal of Forecasting 11 (1995), 271–283. doi:10.1016/0169-2070(94)00573-U.

31.

S.L.

Scott and

H.R.

Varian, Predicting the present with Bayesian structural time series, International Journal of Mathematical Modelling and Numerical Optimisation 5 (2014), 4–23. doi:10.1504/IJMMNO.2014.059942.

32.

Sorjamaa,

Hao,

Reyhani,

Ji and

Lendasse, Methodology for long-term prediction of time series, Neurocomputing 70 (2007), 2861–2869. doi:10.1016/j.neucom.2006.06.015.

33.

R.R.

Sundararajan and

Pourahmadi, Nonparametric change point detection in multivariate piecewise stationary time series, Journal of Nonparametric Statistics 30 (2018), 926–956. doi:10.1080/10485252.2018.1504943.

34.

S.J.

Taylor and

Letham, Forecasting at scale, The American Statistician 72 (2018), 37–45. doi:10.1080/00031305.2017.1380080.

35.

Truong,

Oudre and

Vayatis, Selective review of offline change point detection methods, Signal Processing 167 (2020), 107299. doi:10.1016/j.sigpro.2019.107299.

36.

C.L.

Wu and

K.W.

Chau, Prediction of rainfall time series using modular soft computingmethods, Engineering Applications of Artificial Intelligence 26 (2013), 997–1007. doi:10.1016/j.engappai.2012.05.023.

37.

Xing,

Yuan and

Zhou, A mixtured localized likelihood method for GARCH models with multiple change-points, Review of Economics and Finance 8 (2017), 44–60.

38.

Zhang,

Geng and

Lai, Multiple change-points estimation in linear regression models via sparse group lasso, IEEE Transactions on Signal Processing 63 (2015), 2209–2224. doi:10.1109/TSP.2015.2411220.