Point and interval prediction of crude oil futures prices based on chaos theory and multiobjective slime mold algorithm.

Crude oil is the most important energy source in the world, and fluctuations in oil prices can significantly influence investors, companies, and governments. However, crude oil prices have numerous characteristics, including randomness, sudden structural changes, intrinsic nonlinearity, volatility, and chaotic nature. This makes the accurate forecasting of crude oil prices a difficult and challenging task. In this paper, a hybrid prediction model for crude oil futures prices is proposed, the accuracy and robustness of which are demonstrated via controlled experiments and sensitivity analysis. This study uses a new data denoising method for data processing to improve the accuracy and stability of the predictions of crude oil prices. Furthermore, the chaotic time-series prediction method, shallow neural networks, linear model prediction methods, and deep learning methods are adopted as submodels. The results of interval forecasts with narrow widths and high prediction accuracies are derived by introducing a confidence interval adjustment coefficient. The results of the simulation experiments indicate that the proposed hybrid prediction model exhibits higher accuracy and efficiency, as well as better robustness of the forecasting than the control models. In summary, the proposed forecasting framework can derive accurate point and interval forecasts and provide a valuable reference for the price forecasting of crude oil futures.

Introduction

As the most important global energy source, crude oil has become a pivotal strategic resource. Owing to the COVID-19 pandemic (Li et al., 2021), the demand for crude oil saw an unprecedented decline of 9.3% (a decrease of 9.1 million barrels per day) in 2020. Nevertheless, according to the 2021 BP Statistical Review of World Energy (BP China, 2021), crude oil still accounts for the largest share (31.2%) of global primary energy consumption in 2020. The crude oil price is closely related to socioeconomic factors, international politics, and national security (Sun et al., 2020). For example, the three oil crises (in 1973, 1978, and 1990, respectively) caused by surges in crude oil prices led to a reduction in production and a significant slowdown in economic growth in industrialized countries worldwide (Qin, 2020). Therefore, the prediction of crude oil prices (COP) is crucial for investors, companies, and governments. The forecasting of the COP is also a challenging task because the COP can be influenced by multiple factors, including inflation, exchange rates, supply and demand, international politics, wars, and pandemics (Chai et al., 2018; You et al., 2021). COP time series have many characteristics, including intrinsic nonlinearity, randomness, sudden structural changes, volatility, and a chaotic nature (Cerqueti and Fanelli, 2021; Zhao et al., 2021). With the continuous development of machine learning (ML) models and optimization algorithms, the prediction of COP has become increasingly accurate. Because the international crude oil futures prices (COFP) are generally regarded as the reference prices for the international COP (Zhang et al., 2021), an increasing number of researchers have engaged themselves in the prediction of the COFP.

Based on a review of existing studies, forecasting methods applied to the COFP can be divided into two categories, i.e., single and hybrid models.

The single models used for COFP can be divided into ML and traditional time-series analysis models. Traditional time-series analysis methods such as the vector autoregression model (VAR) (Mirmirani and Li, 2004), generalized autoregressive conditional heteroskedasticity (GARCH) model (Agnolucci, 2009), and autoregressive integrated moving average model (ARIMA) (He et al., 2010) are widely used for COP forecasting. Although these models are effective in capturing the linear information in a time series, they fail to capture its nonlinear information to a certain extent (Abdollahi and Ebrahimi, 2020). Therefore, these models are not very effective when predicting the COFP owing to the nonlinear characteristics of COFP (Sun et al., 2018). Furthermore, these models have many strict assumptions (Li et al., 2019), but these assumptions are often difficult to practically satisfy. To compensate for the shortcomings of traditional time-series analysis models, an increasing number of ML methods have been used to forecast the COFP. ML models are not only free of strict assumptions, but they can also capture nonlinear information in a time series (Abedin, Chi, et al., 2021; Abedin, Moon, et al., 2021). Initially, a multilayer perceptron (MLP) (Guotai et al., 2017a), support vector machine (SVM) (Abedin et al., 2019; Abedin et al., 2019; Guo et al., 2012; Li, Chen, et al., 2020; Li, Wen, et al., 2020), back propagation neural network (BPNN) (Mingming and Jinliang, 2012), and extreme learning machine (ELM) (Wang, Athanasopoulos, et al., 2018) were used to improve the accuracy of COFP forecasts. In recent years, deep learning methods have attracted extensive attention owing to their excellent performance in terms of prediction accuracy and stability (Abedin, Chi, et al., 2021; Abedin, Moon, et al., 2021; Lv et al., 2022; Wang et al., 2022; Wang, Athanasopoulos, et al., 2018; Wang, Du, et al., 2018). Long short-term memory (LSTM) networks (Zhang et al., 2020), deep belief networks (DBN) (Zhang and Ci, 2020), and bidirectional long short-term memory networks (BiLSTM) (Abedin, Moon, et al., 2021) are gradually used in COFP predictions. However, there are still certain drawbacks in using ML models to predict the COFP, such as their poor generalization ability and local optimum problems (Wang et al., 2020).

To address these issues, many hybrid models have been proposed for COFP forecasting. A hybrid model is a combination of several single models through an optimization algorithm. It can inherit the advantages of every single model, thereby improving the prediction accuracy and stability (Guotai et al., 2017b; Hu et al., 2021; Jiang et al., 2020; Khalilpourazari and Doulabi, 2021). In existing studies, the most commonly used optimization algorithms are the genetic algorithm (GA) (Yang et al., 2019), particle swarm optimization (PSO) (Ribeiro et al., 2021), ant lion optimization algorithm (ALO) (P. et al., 2018), frog-leaping algorithm (FLA) (He et al., 2021), and whale optimization algorithm (WOA) (Lin and Zhang, 2021). Only one optimization objective can be set when these optimization algorithms are used to build hybrid models. To simultaneously consider more optimization objectives, researchers have proposed many multiobjective optimization algorithms, such as the multiobjective ant lion optimization algorithm (MOALO) (Wang, Du, et al., 2018) and multiobjective whale optimization algorithm (MOWOA) (Wang et al., 2017). Due to the fact hybrid models constructed using multiobjective optimization algorithms can improve forecasting stability and accuracy, an increasing number of researchers are applying such models to COFP forecasting (Abdollahi and Ebrahimi, 2020; Chai et al., 2018; Zhao et al., 2021). The COFP time series has chaotic characteristics (Wang et al., 2020). In the present study, the Lyapunov exponents of the Brent crude oil futures price (BCOFP) and the West Texas Intermediate crude oil futures price (WCOFP) are also calculated, the results of which (i.e., 0.0376 and 0.0036, respectively, both of which are greater than zero) indicate the chaotic characteristics of these time series. However, existing studies cannot sufficiently extract chaotic information from a time series to improve the COFP prediction accuracy. The characteristics and representative studies of the COFP forecasting models are listed in Table 1.

Table 1

Summary of available studies

Type	Subtype	Model	References	Advantages	Disadvantages
Single model	Linear prediction model	VAR	Mirmirani and Cheng Li (2004)	These models have low complexity, fast computational speed, and effective linear time series prediction ability	The nonlinear time series prediction ability of these models is poor
		ARIMA	He et al. (2010)
		GARCH	Agnolucci (2009)
	Shallow neural networks	MLP	Guotai et al. (2017a, 2017b)	These models have an excellent performance to fit simple functions and fast running speed	The ability of these models to fit complex functions is limited
		BPNN	Mingming and Jinliang (2012)
		SVM	Guo et al. (2012)
		ELM	Wang, Du, et al. (2018), Wang, Athanasopoulos, et al. (2018)
	Deep learning models	DBN	Zhang and Ci (2020)	Deep learning can approximate complex functions and is highly adaptive	These models require a lot of computation and a long-running time
		LSTM	Zhang et al. (2020)
		BiLSTM	Li et al. (2021)
Hybrid model	Single objective	GA	Yang et al. (2019)	Fast convergence and global optimization search	Only one objective can be achieved at a time
		PSO	Ribeiro et al. (2021)
		ALO	Reddy et al. (2018)
		FLA	He et al. (2021)
		WOA	Lin and Zhang (2021)
	Multiple objectives	MOALO	Wang, Du, et al. (2018)); Wang, Athanasopoulos, et al. (2018))	Multiple objectives can be achieved simultaneously	A good balance of convergence and diversity has not yet been achieved
		MOWOA	Wang et al. (2017)
		MOSMA	Li, Chen, et al. (2020)

In this paper, a hybrid prediction model applying time varying filtering for empirical mode decomposition (TVF_EMD) and a multiobjective slime mold algorithm (MOSMA), called TVF_EMD_MOSMA, is proposed for COFP forecasting. Based on this framework, the results of the point forecast (PF) and interval forecast (IF) are derived. First, we decompose the COFP time series using TVF_EMD and eliminate redundant noise series. Second, we obtain the PF results of COFP based on TVF_EMD_MOSMA. Finally, we apply the maximum likelihood estimate (MLE) to estimate the probability density function of the COFP residuals based on the fitting errors. Thereafter, we apply MOSMA to determine the confidence interval adjustment coefficient (CIAC). Subsequently, we obtain the IF results of COFP through the PF results, CIAC, and optimal distribution. The contributions of this study are as follows:

A new hybrid prediction model (TVF_EMD_MOSMA) for COFP is presented. The results of comparative experiments show that TVF_EMD_MOSMA exhibits high accuracy and stability of PF and IF for COFP. The predictive performance of TVF_EMD_MOSMA is efficient and robust.

A novel data denoising method (TVF_EMD) is used for COFP data processing. Based on the idea of “decomposition and combination,” the COFP time series is reconstructed after removing the high-frequency noise using TVF_EMD, thus improving the accuracy of the prediction model.

The proposed MOSMA is applied for COFP prediction for the first time. MOSMA applies an archive component to store all non-dominated Pareto solutions and implements multiobjective optimization based on non-dominated sorting and the crowding distance mechanism. Experimental results show that MOSMA can effectively enhance the prediction accuracy and stability of the hybrid model for COFP prediction.

To obtain an IF with a narrower width and higher prediction accuracy, the CIAC determined using MOSMA is added. The contradiction between the interval prediction accuracy and prediction interval width was balanced, significantly improving the IF performance.

The remainder of this paper is organized as follows. Section 2 introduces the proposed hybrid prediction model and its submodels. Section 3 describes the performance evaluation metrics and data used in the study. Section 4 presents the experimental results. Section 5 further discusses TVF_EMD_MOSMA and describes the results of the sensitivity analysis, the accuracy and stability improvement ratio, and forecasting effect analysis. Section 6 discusses the practical applications, limitations, and scope of future research work. Finally, Sect. 7 presents the conclusions of this study.

Research method

In this section, the basic models, interval estimation theory, and proposed hybrid prediction model are introduced.

Time varying filter empirical mode decomposition

Empirical mode decomposition (EMD) can decompose a COFP time series into a finite number of intrinsic mode function (IMF) signals based on the time-scale characteristics of the COFP time series (Wang, Athanasopoulos, et al., 2018; Wang, Du, et al., 2018). Since the decomposed time series facilitates the extraction of time-series features, EMD is extensively used in time-series forecasting (Wang and Wang, 2020). However, EMD is plagued by marginal effects and pattern confusion which leads to reduced accuracy in a time-series decomposition. Recursive empirical mode decomposition (REMD) has been proposed to alleviate the problems of marginal effects and pattern confusion in EMD (Wang and Wang, 2020). Recently, TVF_EMD was proposed to simultaneously solve the modal separation and intermittent operation problems. Simultaneously, the physical meaning of the model parameters in TVF_EMD is clear, which facilitates parameter selection (Wang, Niu, et al., 2021). The detailed calculation process of TVF_EMD is shown below.

Step 1: The frequency and instantaneous amplitude of the COFP time series are calculated using the Hilbert transform.

In Eqs. (1) and (2), and are obtained by interpolating and . In addition, and are the local minimum and maximum of respectively; and represent the instantaneous mean value and instantaneous envelope; and are obtained by interpolating and .

In Eq. (3), represents the local cut-off probability. Arrange to solve the intermission problem: Define a signal as and use the extreme point of as the node. By approximating the COFP time series through a B-spline interpolation, the approximate result is obtained.

Step 2: Determine the cut-off condition, . If , is considered an IMF. Otherwise, set and repeat Steps 1 and 2.

Step 3: Using these above steps, decompose the COFP time series into multiple IMFs.

Volterra adaptive filter based on phase space reconstruction

To extract the chaotic information in the COFP time series, the Volterra adaptive filter is used. The following introduces the principle of the phase space reconstruction, the calculation of the optimal embedding dimension and delay time, and the principle of the Volterra adaptive filter.

Phase space reconstruction

Takens theorem states that a single variable chaotic time series can be reconstructed into a multidimensional phase space, and this space can contain the chaotic features of the original time series. Thus, the laws and properties of a chaotic time series can be accurately captured (Lin and Zhang, 2021).

Suppose that the COFP time series is and that an m-dimensional vector is formed through the delay time:, where is the optimal embedding dimension, is the predicted delay time, is the phase point in m-dimensional phase space, and is the number of phase points,. In addition, describe the evolutionary trajectory of a dynamical system in the phase space, and thus, the chaotic behavior of a dynamical system can be studied in reconstructed m-dimensional phase space.

Determination of delay time

In this study, the mutual information (MI) method is used to determine the delay time of the COFP time series. The probability of occurrence of in the time series is defined as ; the probability of the occurrence of in time series is defined as ; and the joint probability of occurrence of and in the two series is defined as , where and can be solved based on the probability of occurrence in their respective time series. The joint probability can then be obtained by counting the lattice on the plane . The MI function is as follows.

We apply the when the MI function takes the first minimal value point as the delay time.

Determination of the optimal embedding dimension

In this study, the false nearest neighbor method (FNN) is used to determine the optimal embedding dimension. The calculation procedure for the FNN is as follows.

Step 1: In the embedding space with embedding dimension m, find the Euclidean distance nearest neighbor of all points. The Euclidean distance between and is calculated using Eq. (5).

Step 2: When any pair of nearest neighbors satisfies the following criterion, it is an FNN point.

In Eq. (6), and represent the squared distance between any pair of nearest neighbors at the optimal embedding dimensions of m + 1 and m, respectively, and represents the threshold value.

Step 3: When , the ratio of the FNN points to the total number of phase points is calculated, and m is gradually increased until the ratio is less than 5%. The chaotic attractor geometry is considered to be completely opened and is the optimal embedding dimension at this time.

Volterra adaptive filter prediction model

The Volterra adaptive filter prediction model can predict a chaotic time series using only a small sample of data, and it can automatically track the chaotic motion trajectory. This model has high prediction accuracy for a chaotic time series (Qiao et al., 2020).

Assuming that the input variable is and the output variable is , the Volterra adaptive second-order filtering model is as follows:

The coefficient and input vectors are expressed through Eqs. (8) and (9), respectively.

Based on Eqs. (8) and (9), Eq. (7) can be expressed as

ARIMA model

The ARIMA model is a traditional time-series analysis model that is widely used for time-series forecasting (Ribeiro et al. 2021). The model provides a better prediction for a linearly smoothed time series.

Here, is expressed through Eq. (11).

In Eq. (11), ; indicates the autoregressive coefficient polynomial in ARIMA, and indicates the moving smoothing coefficient polynomial.

ELM

ELM is a feedforward neural network with a single hidden layer. This model structure comprises an input layer, implicit layer, and output layer, similar to an ANN (Lin and Zhang, 2021). The layers are connected to each other using a characteristic mapping function. Information from the input layer is processed by the implicit layer and passed to the output layer, which then derives the calculated value according to the mapping function. Although the random initialization of the parameters improves the generalization of the ELM, it also requires the ELM to add a large number of nodes to achieve accurate training. For large samples, several nodes consume an excessive number of computational resources and may cause overfitting.

Bidirectional long-short term memory model

The bidirectional long-short term memory model (BiLSTM) is divided into two independent LSTMs, and the input sequences are input into the two LSTMs in forward and inverse order for feature extraction. LSTM is a recurrent neural network, and it is proposed to solve the gradient disappearance and explosion problems. The design concept of BiLSTM is to simultaneously obtain the characteristics of the data with information between the past and future (Wang et al., 2020). BiLSTM outperforms a single LSTM approach in terms of efficiency and performance.

MOSMA

The slime mold algorithm (SMA) is a population-based metaheuristic algorithm proposed by Li, Chen, et al. (2020). The slime mold can establish the best pathway for connecting food in a relatively superior manner through a combination of positive and negative feedback. Therefore, SMA adjusts the search path and obtains the optimal result based on a positive and negative feedback system. SMA simulates three different morphotypes in the hunting process of slime mold: finding food, wrapping food, and approaching the food morphotype. The mathematical model of SMA is as follows:where

In Eqs. (12) and (13), indicates the current position of the slime mold, indicates the number of current iterations, and represent the lower and upper bounds of the search range, is the vibration parameter, represents the weight of the slime mold, denotes the optimal fitness, denotes the worst fitness, and denotes the optimal fitness in all iterations. denotes the sequence of sorted fitness values (ascending in the minimum value problem), and and denote a random value within the range of [0,1].

SMA can effectively solve many practical problems (Li, Chen, et al., 2020; Li, Wen, et al., 2020). MOSMA is proposed to effectively achieve multiple goals. It is a multiobjective improvement algorithm based on the SMA algorithm, using an elite non-dominated sorting method to estimate the Pareto optimal solutions. In addition, to ensure the diversity of Pareto optimal solutions, MOSMA added a crowding distance mechanism to increase the coverage of all objectives (Premkumar et al., 2021). The steps for an elitist non-dominated sorting approach are as follows:

Step 1: Calculate the non-dominated results of the objective function.

Step 2: Sort the non-dominated results using non-dominated sorting.

Step 3: Find the non-dominated ranking of all non-dominated results to determine the optimal solution.

The crowding distance () is calculated through the following formula:

In Eq. (14), and denote the maximum and minimum values of the objective function, respectively. The non-dominated ranking and crowding distance () are used to determine the optimal solution.

To optimize the prediction accuracy and robustness of the hybrid prediction model, the multiobjective functions are defined as follows:

In Eq. (15), and denote the actual and predicted COPs, respectively. In addition, denotes the mean of .

Interval forecasting method

Although the PF provides an explicit value, the reliability of providing this value is not given by the PF. Therefore, the IF is proposed to fill in this gap. The IF can provide considerable information to users regarding forecast results (Sun et al., 2020). The common distributions in the field of energy price forecasting are the Gumbel, generalized extreme value (GEV), and gamma distribution (Jiang et al., 2021; Wang, Niu, et al., 2021). In this study, the MLE is used to fit the optimal distribution of the prediction error series. However, the interval width and prediction accuracy of the IF are irreconcilable. To obtain an IF with a narrow width and high prediction accuracy, the CIACs ( and ) are added (Jiang et al., 2021), and are determined using MOSMA. The calculation formulas of the day’s IF are as shown below:where and are the lower and upper bounds of the confidence interval at the confidence level (), denotes the day’s PF result, and and are the quantile and quantile of the optimal distribution, respectively. In addition, is the prediction error sequence of TVF_EMD_MOSMA, and and are the CIACs determined through MOSMA. The objective functions for the IF in MOSMA can be defined as follows:where and denote the actual and predicted COFP, respectively; and and are the lower and upper bounds of the confidence interval at the confidence level (). When the objective functions take the minimum value, the position of the slime mold is the values of and .

Framework of proposed TVF_EMD_MOSMA

The proposed TVF_EMD_MOSMA is introduced in this section. The overall framework is shown in Fig. 1. In this study, the point and interval prediction results of WCOFP and BCOFP are derived using TVF_EMD_MOSMA. The detailed steps are as follows:

Fig. 1

Framework of the proposed TVF_EMD_MOSMA

Using the TVF_EMD method, the high-frequency noise in the COFP time series is filtered out. Consequently, a smooth COFP time series is obtained. Thereafter, the processed COFP time series is used in the TVF_EMD_MOSMA prediction model.

Step 1: Data processing

The processed COFP time series is input into the four submodels (Volterra adaptive filter, ARMIA, ELM, and BiLSTM) as an input variable, and the output is . Here, denotes the predicted result from the submodel. In addition, is combined as the final predicted value by a set of weights determined using MOSMA.

Step 2: Point forecast

The optimal distribution of the prediction error series is determined using the MLE. The CIACs ( and ) are determined using MOSMA. Finally, the confidence intervals with confidence levels of 90%, 95%, and 99% are calculated by Eqs. (15) and (16).

Step 3: Interval forecast

Studied data and performance evaluation metrics

In this section, the studied data (original, training set, and test set data) are presented in detail. The performance evaluation metrics for the PF and IF are also presented.

Studied data

WCOFP and BCOFP have a large international influence and are often regarded as the reference prices for the international oil spot market. Therefore, the study sample includes the daily settlement prices of the WCOFP and BCOFP from 4 January 2010 to 30 September 2021. These data are available on Investing.com (https://cn.investing.com/com-modities/crude-oil-historical-data). The period from 4 January 2010 to 30 September 2020 is selected as the training set for training the prediction model. The remaining data (1 October 2020 to 30 September 2021) are used as the test set to evaluate the model prediction performance. More details of the dataset are presented in Table 2.

Table 2

Detailed description of the studied data (Data from https://cn.investing.com/com-modities/crude-oil-historical-data)

Type	Dataset	Date	Number	Statistics ($)
				Maximum	Minimum	Mean	Std	Kurtosis	Skewness
WTI	All samples	2010/01/04–2021/09/30	3053	113.9300	− 37.6300	68.8658	22.4751	1.9931	0.1061
	Training samples	2010/01/04–2020/09/30	2790	113.9300	− 37.6300	69.7469	23.0555	1.9098	0.0329
	Testing samples	2020/10/01–2021/09/30	263	76.2500	35.7900	59.5181	11.3902	2.0058	− 0.5233
Brent	All samples	2010/01/04–2021/09/30	3034	126.6500	19.3300	76.1399	26.1419	1.7798	0.2194
	Training samples	2010/01/04–2020/09/30	2776	126.6500	19.3300	77.4311	26.7414	1.6886	0.1232
	Testing samples	2020/10/01–2021/09/30	258	79.5300	37.4600	62.2462	11.4875	2.0842	− 0.6089

Performance evaluation metrics

In this section, a system of metrics for evaluating the performances of PF and IF is presented (Wang, Niu, et al., 2021; Wang Wang et al., 2021; Wang et al., 2021; Wang et al., 2021). A detailed description of the metrics is presented in Table 3. The abbreviations used in this study are described in "Appendix 1".

Table 3

Description of the metrics used by the model

Type	Metric	Full name	Calculation formula
PF	MAPE	Mean absolute percentage error	MAPE=1n∑i=1nOP_i-FP_i/OP_i×100%
	RMSE	Root mean square error	RMSE=1n∑i=1nOP_i-FP_i^2
	MAE	Mean absolute error	MAE=1n∑i=1nOP_i-FP_i
	MdAPE	Median absolute percentage error	MdAPE=medianOP_i-FP_i/OP_i×100%
IF	FICP	Prediction interval coverage probability	PICP=1n∑i=1nB_i,B_i=1,ifOP_i∈Lb_i,Ub_i0,else
	AIS	Average interval score	AIS=∑i=1nAi(1-α),Ai(1-α)=-2α·IWi(1-α)-4Lbi1-α-FP_i,ifFP_i<Lbi1-α-2α·IWi(1-α),ifFP_i∈Lbi1-α,Ubi1-α-2α·IWi(1-α)-4FP_i-Ubi1-α,ifFP_i>Ubi1-α
	FINAW	Prediction interval normalized average width	PINAW=1n∑u=1nUbi(1-α)-Lbi1-αmaxOP_u,u=1,2,...,n-minOP_u,u=1,2,...,n

and denote the actual and predicted COFPs, respectively; ; and and are the lower and upper bounds of the confidence interval at the confidence level (). In addition, MAPE and MdAPE are used to reflect the predictive accuracy of the model, and the smaller the value, the higher the accuracy; RMSE and MAE are used to reflect the predictive robustness of the model, and the smaller the value, the higher the robustness; FICP is used to reflect the accuracy of the IF, and the larger the value, the higher the accuracy; FINAW is used to measure the width of the prediction interval, and the smaller the value, the higher the accuracy; and AIS is a composite indicator that reflects the accuracy and stability of the interval forecast, and the larger the value, the higher the accuracy

Experiments and analysis

In this section, two control experiments to test the performance of the proposed TVF_EMD_MOSMA prediction model are introduced in detail. Experiment I demonstrated the excellent PF performance of the proposed hybrid prediction model by setting up multiple sets of controlled trials. Experiment II determined the optimal distribution of the prediction error series using MLE. Confidence intervals with confidence levels of 90%, 95%, and 99% were calculated using MOSMA. The IF performance of the proposed hybrid prediction model was tested.

Experiment configuration

Experimental environment

All experiments in this study were conducted using MATLAB R2020a. The experimental platform was a laptop computer with a 64-bit 1.80 GHz AMD Ryzen7 4800U CPU with a Radeon graphics card and 16 GB of RAM, running on Windows 10.

Model parameter settings

The proposed hybrid prediction model consists of four submodels (Volterra adapter filter, ARIMA, ELM, and BiLSTM) and MOSMA optimization algorithms. Three types of control models have been developed: (1) Benchmark models: BPNN, LSTM, Volterra adapter filter, ARIMA, ELM, BiLSTM, TVF_EMD_BPNN, TVF_EMD_Voterra, TVF_EMD_ARIMA, TVF_EMD_ELM, and TVF_EMD_BiLSTM. (2) The models using different data denoising methods: EMD_MOSMA and REMD_MOSMA. (3) The models using different optimization algorithms: TVF_EMD_SFL, TVF_EMD_SMA, TVF_EMD_MOALO, and TVF_EMD_MOWOA. The detailed parameter settings for the models are listed in Table 4.

Table 4

Parameter settings of the model and optimization algorithm

Method	Parameter	Symbol	Value	Reason
TVF_EMD	Bandwidth threshold	w_b	0.08	Common value
	Max number of IMFs	N_I	50	Preset
	B-spline order	O_b	26	Common value
BPNN	Hidden layer nodes number	N_h	10,15,10	Trial–error manner
	Iteration number	N_i	10,000	Trial–error manner
	Learning rate	L_r	0.1	Trial–error manner
	Training requirements precision	P_tr	0.0001	Preset
LSTM	Hidden layer nodes number	N_h	50	Trial–error manner
	Iteration number	N_i	1000	Trial–error manner
	Learning rate	L_r	0.001	Trial–error manner
Volterra	Time delay	T_d	21(WTI), 14(Brent)	The MI method
	Embedding dimension	D_e	10(WTI), 11(Brent)	The FNN method
	Order	O_n	3	Common value
ELM	Iteration number	N_i	10,000	Preset
	Hidden layer nodes number	N_h	5	Trial–error manner
BiLSTM	Variable learning rate	L_r	0.001(0.5)	Trial–error manner
	Iteration number	N_i	3000	Trial–error manner
	Hidden layer nodes number	N_h	100	Trial–error manner
FLA, SMA	Iteration number	N_i	100	Trial–error manner
	Population size	S_p	200	Trial–error manner
MOALO, MOWOA, MOSMA	Population size	S_p	200	Trial–error manner
	Archive size	S_a	100	Trial–error manner
	Iteration number	N_i	100	Trial–error manner

Experiment I: PF result analysis

Three types of control trials were conducted to demonstrate the superiority of the PF performance of the proposed hybrid prediction model. The results of Experiment I are shown in Fig. 2.

Fig. 2

Results of Experiment I

Comparison with benchmark models

Multiple benchmark models were used as control models to demonstrate that the TVF_EMD_MOSMA outperformed the single model. The PF performances of TVF_EMD_MOSMA and multiple benchmark models are presented in Table 5.

Table 5

PF performance of TVF_EMD_MOSMA and benchmark models

Model	MAPE		MdAPE		MAE		RMSE
	WCOFP	BCOFP	WCOFP	BCOFP	WCOFP	BCOFP	WCOFP	BCOFP
Volterra	1.9240	2.0442	1.7847	1.7114	1.3266	1.4964	1.5559	1.8056
BPNN	1.6181	1.5095	1.3078	1.1053	0.9470	0.9187	1.2434	1.2259
LSTM	2.2360	2.6364	1.6559	1.8858	1.2266	1.4416	1.6239	2.1953
BiLSTM	1.9095	1.5275	1.6279	1.1342	1.1372	0.9336	1.4324	1.2803
ARIMA	2.1412	1.5358	1.6513	1.2029	1.2295	0.9259	1.5927	1.2467
ELM	1.6058	1.5823	1.3558	1.2901	0.9392	0.9756	1.2365	1.2750
TVF_EMD_Volterra	1.2246	1.2122	0.8981	0.8980	0.8654	0.8565	1.1264	1.1004
TVF_EMD_BiLSTM	1.2033	1.1705	0.9863	0.8619	0.7063	0.7225	0.9160	0.9825
TVF_EMD_ARIMA	1.0259	0.9916	0.7460	0.7342	0.5951	0.6039	0.8046	0.8117
TVF_EMD_ELM	0.9317	1.0481	0.7231	0.9022	0.5431	0.6413	0.7114	0.8209
TVF_EMD_MOSMA	0.7538	0.8281	0.5122	0.7109	0.5306	0.5862	0.7283	0.7488

Bold is used in Table to highlight the results of the “TVF_EMD_MOSMA” model proposed in this paper

By comparing the model using the TVF_EMD denoising method and the model without denoising, the PF’s accuracy and stability of the model using TVF_EMD are better than those of the model without denoising. For instance, the evaluation metrics of the ELM using TVF_EMD for WCOFP are , , , and . However, the evaluation metrics of the ELM without denoising for WCOFP are ,,, and . The MAPE, MdAPE, MAE, and RMSE of the TVF_EMD_ELM are much smaller than those of the ELM, indicating that the PF performances of ELM are improved by using the TVF_EMD denoising method.

The benchmark models show different predictive performances for the different datasets. For WCOFP, TVF_EMD_ELM achieved the best PF performance. However, TVF_EMD_ARIMA was the best-performing benchmark model for BCOFP. Therefore, to obtain the optimal prediction performance for different datasets, the proposed hybrid prediction model combines submodels using the optimal weights determined by MOSMA.

For the BCOFP, all evaluation metrics of TVF_EMD_MOSMA are smaller than those of all the benchmark models (,,, and ). For WCOFP, the model comparison results are similar to those of BCOPF, demonstrating that the TVF_EMD_MOSMA prediction model achieves a better PF performance than all benchmark models.

Comparison of different denoising methods

To compare the different denoising methods (TVF_EMD, EMD, and REMD), EMD_MOSMA and REMD_MOSMA were set as the control models. The PF performances of TVF_EMD_MOSMA and the control models with different denoising methods are listed in Table 6.

Table 6

The PF performance of TVF_EMD_MOSMA and the control models

Model	MAPE		MdAPE		MAE		RMSE
	WCOFP	BCOFP	WCOFP	BCOFP	WCOFP	BCOFP	WCOFP	BCOFP
EMD_MOSMA	1.0895	1.4658	0.7828	1.1880	0.7730	1.0298	0.9996	1.2920
REMD_MOSMA	1.1222	1.1112	0.8565	0.9461	0.7701	0.7985	0.9789	1.0035
TVF_EMD_FLA	1.0441	1.4731	0.9211	1.0938	0.7351	1.0483	0.9102	1.5689
TVF_EMD_SMA	0.8556	1.2238	0.7289	0.9758	0.6005	0.8636	0.8229	1.1538
TVF_EMD_MOALO	0.8476	0.8401	0.7517	0.6798	0.5952	0.5958	0.7435	0.7526
TVF_EMD_MOWOA	1.0890	1.4657	0.7823	1.1819	0.7728	1.0295	0.9990	1.2919
TVF_EMD_MOSMA	0.7538	0.8281	0.5122	0.7109	0.5306	0.5862	0.7283	0.7488

Bold is used in Table to highlight the results of the “TVF_EMD_MOSMA” model proposed in this paper

By comparing the evaluation indicators of TVF_EMD_MOSMA, EMD_MOSMA, and REMD_MOSMA, the evaluation indicators (MAPE, MAE, MdAPE, and RMSE) of TVF_EMD_MOSMA were smaller than those of EMD_MOSMA and REMD_MOSMA. The results indicate that the TVF_EMD method shows excellent performance. In other words, TVF_EMD is more suitable for data processing of COFP prediction than the other denoising methods (EMD and REMD).

Comparison of different optimization algorithms

In this section, different optimization algorithms are compared. To demonstrate the superiority of MOSMA, multiple control models (TVF_EMD_SFL, TVF_EMD_SMA, TVF_EMD_MOALO, and TVF_EMD_MOWOA) were constructed by keeping the submodels and denoising methods unchanged and solely changing the optimization algorithm. The parameter settings for these models are listed in Table 4. The PF performances of TVF_EMD_MOSMA and the control models with different optimization algorithms are listed in Table 6.

Comparing the evaluation metrics of TVF_EMD_MOSMA and the control models using different optimization algorithms, only the MdAPE of TVF_EMD_MOSMA for BCOFP is larger than that of the optimal model (TVF_EMD_MOALO) by 0.04. For other cases, TVF_EMD_MOSMA has the smallest evaluation metric values (MAPE, MAE, MdAPE, and RMSE) for the WCOFP and BCOFP. In other words, MOSMA achieves higher stability and accuracy than SFL, SMA, MOALO, and MOWOA. Experiment I helps conclude that the proposed hybrid prediction model is more effective for the PF of COFP.

Experiment II: IF result analysis

In this section, the optimal distribution of the prediction error series is determined using MLE. The confidence intervals with confidence levels of 90%, 95%, and 99% and the IF performance of the proposed hybrid prediction model are presented. The results of Experiment II are shown in Fig. 3.

Fig. 3

Results of experiment II

Selection of the optimal distribution function for COFP prediction error series

The IF can provide more information to users regarding the forecast results. The common distributions in the field of energy price forecasting are the Gumbel, GEV, and gamma distributions (Jiang et al., 2021; Wang, Niu, et al., 2021). Therefore, we chose these distribution functions to fit the distribution characteristics of the COFP prediction error series using the MLE. The and RMSE are used as evaluation metrics to evaluate the fitting effect of the distribution function. And is the correlation coefficient of the fitted distribution (Gumbel, GEV, and gamma) with empirical distribution (observations). The results of fitting the distribution functions of the WCOFP prediction error series (WCOFP_E) and BCOFP prediction error series (BCOFP_E) are listed in Table 7.

Table 7

Results of fitting the distribution functions of WCOFP_E and BCOFP_E

Data	Metric	Gumbel	GEV	Gamma
WCOFP_E	RMSE	0.2910	0.7220	0.2510
	R	0.9740	0.9250	0.9770
BCOFP_E	RMSE	0.3140	1.1060	0.1240
	R	0.9630	0.8430	0.9870

For WCOFP_E, the of the Gumbel and gamma distributions are close, and the RMSE of the gamma distribution is the smallest. Therefore, the gamma distribution is chosen as the distribution function of WCOFP_E. For the distribution function of WCOFP_E, the shape parameter is 5.2 and the inverse scale parameter is 0.4, as estimated using the MLE.

For BCOFP_E, the of the gamma distribution was the largest, and the RMSE of the gamma distribution was the smallest. Therefore, the gamma distribution is considered to be the optimal distribution of BCOFP_E. For the distribution function of BCOFP_E, the shape parameter is 16.63 and the inverse scale parameter is 0.19, as estimated using the MLE.

IF performance of TVF_EMD_MOSMA

After determining the optimal distribution function of COFP_E using the MLE, the CIACs were added to the hybrid prediction model to enable the IF with a narrow width and high prediction accuracy. The confidence intervals with confidence levels of 90%, 95%, and 99% were derived using the PF results, optimal distribution, and CIACs. The IF performance of the proposed TVF_EMD_MOSMA is presented in Table 8.

Table 8

IF performance of the proposed TVF_EMD_MOSMA

Data	Confidence level 1-α%	AIS	FICP	FINAW
BCOFP	90	− 0.5943	0.9008	0.1423
	95	− 0.3441	0.9542	0.1676
	99	− 0.0944	0.9924	0.2706
WCOFP	90	− 0.6363	0.9043	0.1638
	95	− 0.4012	0.9468	0.1718
	99	− 0.0868	0.9894	0.3112

The AIS is used to measure the accuracy of the IF, and a larger AIS indicates a higher prediction accuracy. The FICP is also used to measure the accuracy of IF, and it is the frequency at which observations fall into the prediction intervals. FINAW is used to measure the width of the confidence interval, with smaller values indicating a better IF performance. As presented in Table 8, the value of the FICP is extremely close to the confidence level, indicating that the IF results in this study are reasonable. The values of FINAW and AIS are quite small when the given confidence level is achieved, which indicates that the IF results show high accuracy and narrow interval width. In other words, the IF performance of the proposed hybrid prediction model is excellent.

Discussion

TVF_EMD_MOSMA is further discussed in this section, including the sensitivity analysis, the accuracy and stability improvement ratio, and the forecasting effect analysis.

Sensitivity analysis

Sensitivity analysis was employed to measure the robustness of the model’s predictive performance; in particular, sensitivity analysis was used to measure the effect of a change in the model parameters on the prediction results. This study considered the effect of the variation of two parameters, population size and iteration number, on the model’s prediction performance. The sensitivity analysis indicators are calculated as follows (Wu et al., 2022).where is the evaluation metric (MAPE, MdAPE, MAE, and RMSE) of experiment, is the mean of , and is the number of experiments.

In this study, the population size takes individual values from (100, 150, 200*, and 250), and the number of iterations takes individual values from (50, 100*, 150, and 200). The * indicates the optimal parameter of MOSMA. The pattern indicates that the population size changes and the number of iterations remains unchanged at 100. The pattern indicates that the number of iterations changes and the population size remains unchanged at 200. The results of the sensitivity analysis are listed in Table 9.

Table 9

Results of the sensitivity analysis

Data	Pattern	AsMAPE	AsMdAPE	AsMAE	AsRMSE
WCOFP	S_pN_i¯¯	0.0004	0.0025	0.0002	0.0001
	S_p¯¯N_i	0.0004	0.0059	0.0002	0.0001
BCOFP	S_pN_i¯¯	0.0002	0.0022	0.0001	0.0004
	S_p¯¯N_i	0.0001	0.0008	0.0001	0.0000

The values of all sensitivity analysis indicators are extremely small. Take as an example, is 0.0004. However, the smallest MAPE for WCOFP is 0.8281, which is more than two thousand times as large as . Therefore, the effect of the parameter changes on the prediction performance is small. In other words, the prediction performance of the proposed hybrid prediction model is robust.

Accuracy and stability improvement ratio

The improvement of the proposed TVF_EMD_MOSMA was measured in comparison with the control model. and were proposed to measure the improvements in the prediction accuracy and stability, respectively. and were calculated as follows:where and denote the MAPE and MAE of the control model, respectively, and and denote the MAPE and MAE of the proposed TVF_EMD_MOSMA. The results of and are presented in Table 10.

Table 10

Results of and

Model	WCOFP		BCOFP
	IR_a(%)	IR_s(%)	IR_a(%)	IR_s(%)
Volterra	60.8226	60.0070	59.4903	60.8264
BPNN	53.4161	43.9726	45.1408	36.1938
LSTM	66.2892	56.7463	68.5897	59.3378
BiLSTM	60.5251	53.3454	45.7872	37.2133
ARIMA	64.7967	56.8469	46.0802	36.6911
ELM	53.0586	43.5107	47.6648	39.9139
TVF_EMD_Volterra	38.4473	38.6899	31.6862	31.5614
TVF_EMD_BiLSTM	37.3578	24.8855	29.2525	18.8610
TVF_EMD_ARIMA	26.5256	10.8402	16.4859	2.9233
TVF_EMD_ELM	19.0994	2.3140	20.9895	8.5942
EMD_MOSMA	30.8147	31.3629	43.5053	43.0763
REMD_MOSMA	32.8307	31.1044	25.4770	26.5874
TVF_EMD_SFL	27.8064	27.8241	43.7852	44.0809
TVF_EMD_SMA	11.9011	11.6462	32.3337	32.1214
TVF_EMD_MOALO	11.0696	10.8594	1.4284	1.6113
TVF_EMD_MOWOA	30.7829	31.3451	43.5014	43.0597
Average	39.0965	33.4563	37.5749	32.6658

Compared with the control model, the proposed TVF_EMD_MOSMA showed a significant improvement in the prediction accuracy and stability. For the WCOFP, TVF_EMD_MOSMA improved the prediction accuracy and stability by an average of 39.0965% and 33.4563%, respectively. For the BCOFP, TVF_EMD_MOSMA improved the prediction accuracy and stability by an average of 37.5749% and 32.6658%, respectively. Compared with LSTM, TVF_EMD_MOSMA achieved the largest improvement in the prediction accuracy ( and ). Compared with Volterra, TVF_EMD_MOSMA has the largest improvement in prediction stability ( and ). Compared with the benchmark model, the model using different denoising methods, and the model using different optimization algorithms, the proposed hybrid prediction model significantly improved prediction accuracy and stability.

Forecasting effect analysis

The first- and second- order effectiveness ( and , respectively) were introduced to measure the forecasting effect of the model (Wang et al., 2021b, c, d; Wang, Niu, et al., 2021).where and denote the actual and predicted COFP, respectively, and is the number of predicted values. In addition, and are larger, indicating a higher predictive efficiency of the model. The evaluation results regarding the forecasting efficiency of the multiple models are listed in Table 11.

Table 11

Evaluation results of multiple models’ forecasting efficiency

Model	WCOFP		BCOFP
	TE_1	TE_2	TE_1	TE_2
TVF_EMD_MOSMA	0.9925	0.9853	0.9921	0.9856
TVF_EMD_MOALO	0.9915	0.9850	0.9916	0.9851
TVF_EMD_SMA	0.9914	0.9832	0.9878	0.9768
TVF_EMD_ELM	0.9907	0.9825	0.9908	0.9827
TVF_EMD_ARIMA	0.9898	0.9804	0.9901	0.9812
TVF_EMD_SFL	0.9896	0.9820	0.9853	0.9696
TVF_EMD_MOWOA	0.9891	0.9803	0.9853	0.9741
EMD_MOSMA	0.9891	0.9802	0.9853	0.9741
REMD_MOSMA	0.9888	0.9799	0.9889	0.9805
TVF_EMD_BiLSTM	0.9880	0.9779	0.9883	0.9776
TVF_EMD_Volterra	0.9878	0.9776	0.9879	0.9782
ELM	0.9839	0.9704	0.9842	0.9710
BPNN	0.9838	0.9704	0.9849	0.9718
BiLSTM	0.9809	0.9667	0.9847	0.9704
Volterra	0.9808	0.9688	0.9796	0.9661
ARIMA	0.9786	0.9610	0.9846	0.9710
LSTM	0.9776	0.9559	0.9736	0.9367

The results in Table 11 demonstrate that the first- and second- order effectiveness ( and , respectively) of the proposed TVF_EMD_MOSMA are the highest. In other words, the proposed hybrid prediction model achieved the highest prediction efficiency in comparison with the other control models. By comprehensively considering the sensitivity analysis, accuracy and stability improvement ratio, and forecasting effect analysis, the proposed TVF_EMD_MOSMA can achieve a higher forecasting accuracy, better forecasting stability, and higher forecasting efficiency compared with other control models. Therefore, the proposed hybrid prediction model is reliable, valid, and significant.

Practical applications and limitations of the model

In this section, the practical applications and limitations of the proposed TVF_EMD_MOSMA, as well as future research in this area, are presented.

Practical applications

COPs are closely related to socioeconomic, international, and national security. However, international COPs are volatile and uncertain, owing to various factors. The proposed TVF_EMD_MOSMA for the prediction of COFP is reliable, valid, and significant. It can provide valuable reference information for investors (Medina–Olivares et al., 2021), companies, and governments.

By considering a combination of the spot prices and the forecast prices of COFP, investors can decide on an investment strategy to achieve their profit goals.

Companies that use crude oil as a feedstock or produce oil for sale can apply hedging operations based on the forecast results of COFP to control their production costs and sales risks.

Governments can decide on the import, export, storage, and use of crude oil according to the forecast results of the COFP to ensure the stability of domestic oil prices. Stable domestic oil prices are crucial for the smooth development of economies and society.

Limitations and future research

The proposed TVF_EMD_MOSMA uses only the COFP time series to predict future prices. COFP can be influenced by multiple factors, including inflation, exchange rates, supply and demand, international politics, wars, and epidemics. Hence, these factors should be considered in future studies. In the future, the proposed TVF_EMD_MOSMA should be extended to the forecasting of other energy prices, such as coal and natural gas prices.

Conclusion

Crude oil is the most important energy source in the world, and fluctuations in crude oil prices can have a significant impact on investors, companies, and governments. Therefore, accurate prediction of COFP is crucial. In this paper, the TVF_EMD_MOSMA prediction model is proposed to improve the accuracy and robustness of the prediction. A new data denoising method, TVF_EMD, is used for COFP data processing. The chaotic time-series prediction method, shallow neural networks, linear model prediction methods, and deep learning methods are adopted as submodels for COFP prediction. The predicted values of submodels are combined with the optimal weight that is determined using MOSMA. The results of IF with a narrower width and higher prediction accuracy were derived by introducing CIACs determined using MOSMA. The conclusions of this study are as follows.

The new data denoising method, TVF_EMD, can significantly improve the prediction accuracy of COFP. Comparison experiments helped determine that the prediction accuracy of the model using TVF_EMD was significantly higher than that of other denoising methods (EMD and REMD).

The chaotic time-series prediction method, shallow neural network, linear prediction model, and deep learning method were adopted as submodels. Combining their prediction results with MOSMA can obtain accurate and stable prediction results. The MAPE and MAE of WCOFP and BCOFP were ,,, and , respectively.

The PF performance of the proposed TVF_EMD_MOSMA is highly robust. The sensitivity analysis demonstrated that the variation in the model parameters had a slight effect on the prediction performance. The maximum value of the sensitivity analysis indicator is 0.0059, which is extremely small.

The IF performance of the proposed TVF_EMD_MOSMA is excellent. By introducing the CIAC determined using MOSMA, the contradiction between the prediction accuracy and interval width is balanced.

Table 12

List of abbreviations

Abbreviation	The full name
AIS	Average interval score
ALO	Ant lion optimization algorithm
ARIMA	Autoregressive integrated moving average model
BiLSTM	Bidirectional long short-term memory networks
BCOFP	Brent crude oil futures price
BPNN	Back propagation neural network
COFP	Crude oil futures price
COFP_E	COFP prediction error series
COP	Crude oil price
CIAC	Confidence interval adjustment coefficient
DBN	Deep belief network
EMD	Empirical mode decomposition
ELM	Extreme learning machine
FICP	Prediction interval coverage probability
FLA	Frog leaping algorithm
FINAW	Prediction interval normalized average width
FNN	The false nearest neighbor method
GARCH	Generalized autoregressive conditional heteroskedasticity model
GEV	Generalized extreme value
IF	Interval forecasts
LSTM	Long short-term memory network
MdAPE	Median absolute percentage error
MAE	Mean absolute error
MLE	Maximum likelihood estimate
MAPE	Mean absolute percentage error
MI	The mutual information method
ML	Machine learning model
MLP	Multi-Layer Perceptron
MOALO	Multiobjective ant lion optimization algorithm
MOSMA	Multiobjective slime mold algorithm
MOWOA	Multiobjective whale optimization algorithm
PF	Point forecasts
PSO	Particle swarm optimization
RMSE	Root mean square error
REMD	Recursive empirical mode decomposition
GA	Genetic algorithm
SMA	Slime mold algorithm
SVM	Support vector machine
TVF_EMD	Time varying filtering for empirical mode decomposition
VAR	Vector autoregression model
WCOFP	West texas intermediate crude oil futures price
WOA	Whale optimization algorithm