Modeling the COVID-19 epidemic in Bolivia

1.   Introduction

The stock price is the primary factor affecting investors. Short-term stock price forecasts have important significance for reference to investors. The factors affecting stock price are very complex, which makes it difficult to analyze and predict stock price. The stock market can be seen as a nonlinear dynamic system under the action of multiple factors. Time series forecasting models are usually adopted for stock price forecasts. There are multiple regression models, exponential smoothing models, ARIMA models, neural network models, combined forecasting models, and grey forecasting models.

Many scholars have conducted a series of studies on stock price prediction. Shi ^[1] used ARMA and BP neural networks to obtain the linear and non-linear components of the stock price sequence and then used the Markov model to correct the error of the result. Vijh ^[2] used the random forest algorithm to integrate the indicators of multiple stocks, and then used BP neural network to output the predicted observations. Liu ^[3] proposes a new effective metric, Mean Profit Rate (MPR). The effectiveness of metrics is measured based on the correlation between the metric value and profit of the model. And, the relationship between the features of Bitcoin and the next day change in the price of Bitcoin using an Artificial Neural Network ensemble approach called Genetic Algorithm based Selective Neural Network Ensemble ^[4] is explored by Sin. Zhu ^[5] proposes an intelligent trading system using support vector regression optimized by genetic algorithms (SVR-GA) and multilayer perceptron optimized with GA (MLP-GA). And, an improved one-step ahead forecasting system and an improved multi-step ahead prediction ^[6] are proposed by Dong. Zhang ^[7] uses BP neural network to train stock prices in the form of spline division, and then output the predicted observations with data independent of the training samples. Teo ^[8] train the wavelet packet multi-layer perceptron neural network (WP-MLP) by backpropagation for time series prediction. Summarizing the above literature, using deep learning and machine learning methods to predict stock prices can achieve good prediction results. However, most methods do not consider the memory of the sequence when processing the original price sequence. Sang ^[9] studied the characteristics of long-term memory of fluctuations in the Chinese stock market, which means that the historical data of the stock market have an impact on the stock price at the current time. Since fractional calculus has power-law memory ^[10], the idea of fractional calculus can be used to obtain information about historical data. Based on the characteristics of fractional calculus, some models that combined the idea of fractional calculus with traditional time series forecasting methods have attracted the attention of many scholars. For example, the ARFIMA model, fractional grey model, etc. are the models that combine fractional ideas. Because the grey model has good performance in the uncertainty system, the prediction of the stock price is very consistent with the grey model.

The Grey model is an important part of the grey theory, and the most famous one is GM (1, 1) model. However, GM (1, 1) has many shortcomings. The background value is based on approximation, and it cannot handle high-dimensional data. And the ability of nonlinear fitting is insufficient. Because of the shortcomings, many scholars had improved the GM (1, 1) model ^{[11,12,13,14,15,16,17,18,19,20,21,22,23,24,25]}. Tien ^[13] proposed the GMC (1, m) model to solve the shortcomings of traditional GM (1, 1) that cannot handle high-dimensional data. Literature ^{[14,15,16,17,18,19]} proposed a variety of grey Bernoulli prediction models and grey models combined with kernel methods, which solved the shortcomings of insufficient nonlinear fitting capabilities. Literature ^[20,21] had solved the shortcoming that the background value is based on approximation. Literature ^{[22,23,24,25]} combined fractional calculus with the grey model and used historical data to solve many practical application problems.

Since most time series including stock price series usually contain linear and non-linear factors, a single forecasting model cannot achieve good results. According to the high-dimensional, memorable, and non-linear characteristics of stock data, a Hausdorff fractional grey model with convolution, locally linear embedding, and BP neural network (HFGMC (1, m)-LLE-BP) model is proposed in this paper. HFGMC (1, m) model is improved from FGMC (1, m) in literature ^[25]. The FGMC (1, m) model has several disadvantages. There are factorial operators in the elements of the accumulation generation matrix in the FGMC (1, m) model. When the data length is long, it will have the disadvantages of complex calculation and information loss. The calculation of background value in the FGMC (1, m) model is approximately obtained by using the trapezoidal formula between the two values in the sequence, which will calculate background value ignore the local change trend of the sequence. These shortcomings degrade the prediction performance of FGMC (1, m). The operator derived from HFGMC (1, m) model combined with Hausdorff fractional derivative does not contain a factorial operator, which solves the shortcomings of complex calculation and information loss. Because of the disadvantages of the background value based on approximation, the HFGMC (1, m) model uses the Newton-Cotes discrete formula and the Newton mean difference interpolation polynomial to optimize. Because the Newton-Cotes discrete formula and the Newton mean difference interpolation polynomial consider the changing trend of local data, it can solve the shortcomings of approximate background values. HFGMC (1, m) is used to predict the linear component of stock price. The BP neural network is used to predict the nonlinear component of stock price. The prediction result is obtained by adding a linear component and a nonlinear component.

This paper uses the stock closing price as the forecast target to verify the performance of HGFMC (1, m)-LLE-BP. Compared with FGMC (1, m), the prediction accuracy of HFGMC (1, m) is much higher than that of FGMC (1, m). This is because FGMC (1, m) only retains local information, while HGFMC (1, m) retains global information. Therefore, the prediction accuracy of HGFMC (1, m) will be higher. Besides, thanks to the excellent nonlinear fitting ability of the BP neural network, it can fit the nonlinear residuals in the time series. The prediction results show that the BP neural network improves the performance of HGFMC (1, m). Finally, compare the HGFMC (1, m)-LLE-BP model with the FGMC (1, m)-LLE-BP model, BP neural network model, LSTM model and the ARIMA-BP hybrid model. The performance of the HGFMC (1, m)- LLE-BP model is far better than other models.

The rest of this paper is organized as follows. In Section 2, the specific framework of the proposed model is given. In Section 3, the FGMC (1, m) is described. In Section 4, the HFGMC (1, m) is given. In Section 5, the LLE algorithm and BP neural network are introduced. In Section 6, experiments have verified the effectiveness of the proposed model. In Section 7, the work of this paper is summarized.

2.   Specific frameworks of the proposed model

In this paper, the HFGMC (1, m) is proposed to predict the linear component of the stock price, and then the BP neural network is used to predict the nonlinear component of the stock price. Because stock price is affected by nonlinear high-dimensional stock data, this paper used the LLE algorithm to reduce the dimensionality of nonlinear high-dimensional stock data. The flow chart of the model proposed in this paper is shown in Figure 1.

The text description of the specific steps is as follows.

Step 1: Input parameters r, K. The parameter r is used to calculate the r-order accumulative generator matrix. The parameter K represents the number of neighboring points of each sample point, which is used for the calculation of locally linear embedding.

Step 2: Input the original data matrix. Suppose the original data matrix ${X^{\left( 0 \right)}}'$, $X' \in {R^{n \times m}}$. n is the length of the sequence, and m is the m stock indicators including the sequence to be predicted. Select the sequence to be predicted (usually the first column of the original data matrix), and the remaining sequences are used as related factors that affect the sequence to be predicted (usually other columns in the original data matrix except for the first column, and each column represents a related factor).

Step 3: Form the correlation factor matrix that affects the sequence to be predicted, and apply the locally linear embedding algorithm to reduce the dimension to d dimension. The sequence to be predicted and the correlation factor matrix after dimensionality reduction form the input matrix ${X^{\left( 0 \right)}}$. Then normalize the input matrix.

Step 4: Perform Hausdorff fractional accumulative generation operation (HFAGO) on the input matrix to obtain the Hausdorff fractional accumulative matrix.

Step 5: Combine the Newton-Cotes formula and use the Hausdorff fractional accumulation matrix to calculate the background value.

Step 6: After obtaining the optimized background value, combine the least square estimation to solve the parameters.

Step 7: Solve the whitening equation and get the time response function. Solve the in-sample time response sequence.

Step 8: Perform the inverse Hausdorff fractional accumulation generation operation (I-HFAGO) on the time response sequence within the sample to obtain the in-sample prediction sequence.

Step 9: Subtract the in-sample prediction sequence and the in-sample true value sequence to obtain the in-sample residual sequence. The mean of the sum of the residual sequence elements is the mean absolute error (MAE). Update r and K, and then judge whether the parameters r and K exceed the search range. If it exceeds the search range, select r, K that minimize MAE to solve the in-sample prediction sequence, and then go to step 10. If it does not exceed the search range, repeat Step 2 to Step 9.

Step 10: Input the out-of-sample correlation factor matrix R, and then use the time response function to solve the out-of-sample time response sequence. Then the out-of-sample time response function obtains the out-of-sample prediction sequence through I-HFAGO as a linear component.

Step 11: The in-sample residual sequence uses BP neural network to do nonlinear fitting, and then the out-of-sample residual sequence is output as the nonlinear component.

Step 12: Add the linear component and the non-linear component to get the combined prediction sequence.

The above steps can be divided into 6 modules. Step 1 to Step 3 is the data preprocessing part. Step 4 to Step 7 is the modeling process of the HFGMC (1, m). Step 8 to Step 9 is the residual calculation part. Step 10 is to calculate the linear component through the HFGMC (1, m). Step 11 is to calculate the non-linear part through BP neural network. Step 12 gives the prediction result.

3.   The FGMC (1, m) model

Literature ^[25] gave the specific derivation process of FGMC (1, m). Let the original data matrix be ${X^{\left( 0 \right)}}$, $X \in {R^{n \times m}}$. n is the length of the sequence, and m is the indicators including the sequence to be predicted. Then perform the r-order accumulative generation operation (r-AGO) on each column of the original data matrix to obtain the r-order accumulative generation sequence. The mathematical form is defined as

where $\left( {\begin{array}{*{20}{c}} {t - i + r - 1} \\ {t - i} \end{array}} \right){\rm{ = }}\frac{{\left( {t - i + r - 1} \right)\left( {t - i + r - 1} \right)...\left( {r + 1} \right)r}}{{\left( {t - i} \right)!}}{\rm{ }}$, and r is called the fractional order. r should be taken in the interval (0, 1). However, r can be an arbitrary real number. The matrix form of the above equations can be written as

where ${A^r}_{n \times n}$ is called the r-order accumulative matrix, and its form is

Then establish the whitening equation. The specific mathematical form is

where $\tau$ is the delay period and u is the control coefficient. When $\tau$ and u are equal to 0, the traditional GM (1, m) model can be obtained. But the GM (1, m) had been proved inaccurate in ^[13]. When $\tau$ is equal to 0 and m is equal to 1, it can degenerate into a GM (1, 1) model.

By integrating the Eq (5) and using the trapezoidal formula in the interval [k-1, k], the discrete form of the existing fractional grey model can be obtained as

where $Z_j^{\left( r \right)}$ is called the background values, which are defined as

The parameter ${b_1}, {b_2}, \ldots , {b_m}, u$ of Eq (5) can be obtained by least square estimation. Then substitute it into the whitening equation of Eq (4), and solve the differential equation to obtain the time response function. The form of the time response function is

The convolution operation of the second term on the right side of the Eq (9) can be discretized using the Gaussian formula to obtain the following discrete equation. The specific mathematical form is

where $f\left( \upsilon \right) = \sum\limits_{j = 2}^m {{b_j}x_j^{\left( r \right)}\left( \upsilon \right)} + u$. Finally, perform the inverse r-order accumulative generation operation (r-IAGO) to obtain the predicted sequence

4.   The HFGMC (1, m)

Accumulate generation operation (AGO) is the core step of the grey model, which can eliminate the instability caused by fluctuations in the original data. From the r-order accumulation matrix in Eq (4), it can be seen that the original sequence needs to be multiplied by a factorial multiplier when the accumulation operation is performed. When the length of the sequence is long, the calculation in the r-order accumulation matrix will be very complicated. Besides, it can be seen from Eqs (7) and (8) that the calculation of the background value is only obtained by using a simple trapezoidal calculation formula. However, most systems do not grow linearly, which leads to calculation errors when using FGMC (1, m) models. The error is shown in Figure 2.

4.1.   Derivation of Hausdorff Fractional Accumulate Generation Operation (HFAGO)

The Hausdorff fractional derivative is defined as ^[26].

Since the first difference can usually be approximately defined by the first derivative as

Therefore, the difference of Hausdorff fractional derivative can be approximately defined by Hausdorff fractional derivative as

Let $f\left( t \right) = {x^{\left( 0 \right)}}\left( 1 \right), f\left( {t - 1} \right) = {x^{\left( 0 \right)}}\left( 2 \right), \ldots , f\left( 1 \right) = {x^{\left( 0 \right)}}\left( t \right)$ (This is the reverse mapping of the sequence when doing AGO operations), then 1-order reverse accumulate operation (1-RAGO) can be defined as

1-RAGO is essentially the same as 1-AGO because it just does a simple inverse mapping of the sequence.

Eq (15) combined with Eq (13) to get

In the background of the grey prediction model, the difference and accumulate operation are reversible ^[24]. Eq (16) can be also written as

Without loss of generality, HFAGO is expressed as ${\nabla ^r}$ and satisfies the following relationship

Combine Eq (17) with Eq (14), then move the term to the right to get

Combining Eq (17), Eq (15) is simplified to

Equation (20) is HFAGO. Compared with Eq (2), HFAGO does not contain a factorial multiplier. To understand the difference between r-AGO and HFAGO more intuitively, the next step is to visually compare r-AGO and HFAGO. Assuming

where $t$ is equal to $1, 2, \ldots , n$. The $\Gamma \left( \bullet \right)$ in Eq (21) is the gamma function. Using Eqs (21) and (22), Eq (4) can be rewritten as

Assuming that r = 0.5, visualized Eqs (21) and (22). The visual comparison is shown in Figure 3.

In Figure 3, the solid line represents the change of the multiplier overtime when r-AGO is performed. The dotted line represents the change of the multiplier over time when performing HFAGO. The analysis shows that the r-AGO is 0 when t > 172. This is because the factorial number is too large, so calculate the value of the denominator in matlab2020a will get INF. It will result in the loss of subsequent information and make the calculation result inaccurate. The visualization curve of HFAGO is very similar to the traditional r-AGO. On this basis, HFAGO can retain global information and reduce the amount of calculation.

4.2.   Optimization of background value based on Newton-Cotes formula

Integrate both sides of Eq (5) on the interval [k-1, k] at the same time.

To optimize the error caused by the trapezoidal formula, it is necessary to use the Newton-Cotes formula to calculate the optimized background value for the integral term in Eq (25). Assume the integration area [a, b] is divided into l equal parts. Select node x_i = a+ih to construct a discrete integral formula, where h is the step size. Then, the constructed discrete integral formula is called the Newton-Cotes formula. The mathematical form is defined as

The C^(l) in Eq (26) is called the Cotes coefficient. When l = 1, it is a trapezoidal formula. when l = 4, the corresponding quadrature formula is

This paper divided the integration area into 4 equal parts. Then, the integral term in Eq (25) is constructed using Eq (27) to obtain the optimized background value. The optimized background value is

It can be noticed that ${x^{\left( r \right)}}\left( {k + \frac{1}{4}} \right)$, ${x^{\left( r \right)}}\left( {k + \frac{1}{2}} \right)$, ${x^{\left( r \right)}}\left( {k + \frac{3}{4}} \right)$ do not exist for discrete data, and the Newton means difference interpolation polynomial can be used to solve this problem. The form of the Newton means difference interpolation polynomial is

where $f\left[ {{x_0}, {x_1}, \ldots , {x_n}} \right]$ is the (n-1) order mean difference of the function $f\left( x \right)$, and n represents the data points involved in the construction of the interpolation polynomial. However, the more points involved in the construction of the polynomial are not the better. When n > 8, the Runge effect usually appears ^[27]. This paper will use every 4 data points to construct a Newton mean difference interpolation polynomial.

5.   Local linear embedding algorithm and BP neural network

In this section, the LLE algorithm ^[28] and BP neural network ^[29] are introduced.

5.1.   Local linear embedding (LLE) dimensionality reduction algorithm

The LLE algorithm was proposed by S.T. Roweis and L.K. Saul in 2000. The LLE algorithm is an important dimensionality reduction algorithm in flow learning. The principle is to assume that the data in a high-dimensional space is linear in a region, and certain data can be expressed linearly using data samples in its neighborhood. Then map the data to a low-dimensional space, and try to keep the weight coefficients of the linear relationship before and after the projection the same.

The algorithm steps are as follows.

Step 1: Find the K neighboring points of each sample point ${X_i}$ in the high-dimensional space, and then calculate the Euclidean distance between them.

Step 2: Solve the covariance matrix between neighboring points.

Step 3: Solve the weight coefficient vector of each sample point.

Step 4: The weight coefficient vector forms the weight coefficient matrix W, and then calculates the matrix M that can minimize the reconstruction error function.

Step 5: Solve d+1 smallest eigenvalues in matrix M. The eigenvector corresponding to the eigenvalues from 2 to d+1 are the original data matrix X corresponding to the mapping space with the smallest reconstruction error.

5.2.   BP neural network

BP neural network was proposed by Hinton and McClelland in 1986. The core idea is based on the multi-layer forward neural network combined with the error loss gradient descent direction. The error loss gradient descent direction is combined with the error backpropagation algorithm to update the weights of the multilayer neural network. After enough iterations, the error loss converges to the minimum and finally can approximate any nonlinear function.

6.   Experiment analysis

Through the above method and combined with the process in Figure 1, experiments are carried out to verify the performance of the model proposed in this paper. The first step is to select experimental examples. The second step is to select parameters. The third step is to compare the predicted value with the real value. The last step is to evaluate the performance of the model proposed in this paper through error indicators.

6.1.   Experimental examples

This paper selected the stock closing price as the forecast target and forecasted the closing price for the next 10 days to verify the performance of the model proposed in this paper. This paper divides the predicted stocks into two categories. The first category is stocks with large fluctuations in the closing price of the stock. For example, within the range of a given trading day (excluding public holidays), The stocks with the highest closing price minus the lowest closing price greater than 5 are the first type of stocks. The second category is stocks with small fluctuations in the closing price of the stock. For example, within the range of a given trading day, the stocks with the highest closing price minus the lowest closing price less than 2 are the second type of stocks. In order to verify the effectiveness of this model for the above two types of stocks, this paper selects 6 stocks from the Shanghai Stock Exchange (code prefix is SH) and Shenzhen Stock Exchange (code prefix is SZ) as prediction objects. The 6 stocks are, SZ002607, SZ002567, SZ000796, SH600586, SZ002770 and SZ000573. SZ002607, SZ002567 and SZ000796 belong to the first category of stocks. SH600586, SZ002770 and SZ000573 belong to the second category of stocks. Stock data were crawled on NetEase Finance and Economics website using a web crawler program.

Figure 4 shows the stock closing price trend chart for a total of 256 trading days from 2019-1-16 to 2020-2-10. The closing price of SH600586, SZ002770 and SZ000573 has a small variation range, while the closing price of SZ002607, SZ002567 and SZ000796 has a large variation range.

Table 1 shows the data used for modeling in SH600586, and the modeling data of the remaining 5 stocks are shown in the Supplementary. In this paper, the highest price, lowest price, opening price, previous closing price, floor change, floor, exchange, trading volume, trading amount, total market value and circulating market value are selected as indicators that affect the closing price of stocks. Indicators such as the highest price, lowest price, opening price, previous closing, floor change, floor and the exchange can reflect the changing trend of the closing price of stocks in the recent period. Indicators such as trading volume, trading amount, total market value and circulating market value can reflect the financial status of the stock company corresponding to the stock. Imitating the method of literature ^[30], the above indicators are reduced by dimensionality reduction algorithm, and then input into the HGFMC (1, m) model.

6.2.   Discussion of hyperparameters

First, the delay period parameter $\tau$ is discussed. It can be noted that there are 10 trading days' indicator data after the serial number 256. The indicator data of these 10 trading days is used to predict the closing price of the next 10 days. Since the HGFMC (1, m) model is a time series forecasting model that processes high-dimensional data, the model needs a priori data of the next 10 days to predict the stock closing price in the next 10 days. Considering that investors are not sensitive to changes in stock data ^[31] (except for market crises), this makes investors have a certain delay in responding to stock data. In other words, the stock data at the moment will have a time delay to affect investor sentiment, and then investor sentiment will affect the stock price. It can be seen from Eq (10) that HGFMC (1, m) integrates the delay function and shows the relationship between the stock price and various indicators. Therefore, when used Eq (10) to solve the time response function, the delay period is set to 10. The setting of the delay period $\tau$ is equal to the amount of prior data required.

Then the dimension d of the low-dimensional mapping space of the LLE algorithm is discussed. Literature ^[30,32,33] shows that setting the dimension d of the low-dimensional space to 2, 3, and 4 can well express the feature space of the data. Therefore, this paper uses dimension 3 as the dimension of the low-dimensional mapping space of the LLE algorithm.

6.3.   Discussion on the stability of stock prices

The stationary test of stock prices is discussed in this section. This paper uses the ADF test to test the stability of stock prices.

The ADF test can check whether there is a unit root in the sequence. If there is a unit root, the sequence is not stationary. First, ADF test is performed on the original price series.

As shown in Table 2, the original sequence of 6 stocks cannot reject the existence of unit roots at the significance levels of 1, 5, and 10%. Therefore, the original sequence of 6 stocks can be considered unstable. Since the stock price series is not stationary, it is necessary to stabilize the series before forecasting.

From Eqs (1) and (20), it can be known that after the original sequence uses HFAGO, the random walk of the price appears as the increment between data points on the accumulative generation sequence. Therefore, the original sequence can be considered stationary after using HFAGO.

Figure 5 shows the stabilization process of SZ002770's stock price. The left image of Figure 5 shows the original sequence. The right figure of Figure 5 shows the accumulative generation sequence. The stationary test of the accumulative generation sequence of 6 stocks also uses the ADF test for further verification.

It can be seen from Table 3 that SZ000573, SZ002567, SZ002607, SZ002770 and SH600586 reject the hypothesis at the significance levels of 1, 5, and 10%, while SZ00796 rejects the hypothesis at the significance levels of 5% and 10%. Therefore, it can be considered that the accumulative generation sequence obtained by the original sequence after using HFAGO is stationary.

6.4.   Select parameters

It can be seen from Figure 1 that the establishment of the model requires the selection of appropriate parameters r and K. This paper established the contour diagram between the parameters r, K and the mean absolute error (MAE), and then found the parameter that could minimize the MAE.

Figure 6 is a contour diagram between the parameters r, K and MAE. In the figures, the horizontal axis represents the order r, and the vertical axis represents the number of K neighboring points of each sample point. The bluer the color, the smaller the MAE. The yellower the color, the larger the MAE. This paper enumerated the order r in the range of 0 to 1 in steps of 0.1 and enumerated the number of neighboring points K of the sample points from 0 to 100 in steps of 5. Then, select the search range that minimizes the MAE. Finally, use the exhaustive method within the range to get the relative optimal parameters.

Table 4 shows the optimal parameter range of each stock. The search step for parameter r is 0.01, and the search step for parameter K is 5. For all stocks, the optimal search range of parameter K is 45 to 100. For SZ000573, SZ002770 and SH60586, the search range of parameter r is 0.25 to 0.6. For SZ000796, the search range of parameter r is 0.1 to 0.45. For SZ002567 and SZ002607, the search range of parameter r is 0.35 to 0.7.

Table 5 shows the result of parameters. For SZ000573, when the order r is 0.30 and K is 60, the MAE is 0.6120. For SZ000796, when the order r is 0.14 and K is 50, the MAE is 1.2227. For SZ002567, when the order r is 0.35 and K is 95, the MAE is 2.4767. For SZ002607, when the order r is 0.35 and K is 65, the MAE is 1.1929. For SZ002770, when the order r is 0.44 and K is 65, the MAE is 0.8013. For SH600586, when the order r is 0.41 and K is 55, the MAE is 0.6336.

6.5.   Discussion between HFGMC (1, m) and FGMC (1, m)

To make the distribution of the data input to FGMC (1, m) the same as HGFMC (1, m), the parameters of the data input to FGMC (1, m) also need to be set according to Table 5. The comparison of HFGMC (1, m) and FGMC (1, m) uses the mean error between the in-sample predicted sequence and the real sequence as an indicator. Table 6 shows the MAE comparison between HFGMC (1, m) and FGMC (1, m).

It can be seen from Table 6 that the improvement of FGMC (1, m) by HGFMC (1, m) is very obvious, and the average MAE is reduced by 7.3778.

6.6.   Discussion of BP Neural Network

In this section, the improvement of the model by the BP neural network is discussed. It can be known from Eq (10) that the out-sample forecast data can be regarded as a linear combination of various indicators, so a single HGFMC (1, m) can only predict the linear part of the future data. Due to the strong nonlinear fitting ability of the BP neural network, it can fit the residual function in the in-sample prediction sequence. Therefore, using BP neural network can get the nonlinear part of future data. Next, compare the model performance of HGFMC (1, m) and HGFMC (1, m) combined with the BP neural network.

As shown in Table 7, except for the slight decrease and no change in the forecast accuracy of SZ000796 and SH600586, the forecast accuracy of other stocks has been significantly improved.

6.7.   Forecast result

This paper chose FGMC (1, m)-LLE-BP, long short-term memory neural network (LSTM), BP neural network (BP), ARIMA-BP ^[34] as the comparison model. BPNN and LSTMNN used 10 lagging data of the stock closing price to input the neural network in the form of a vector, and output the current value. For example, the input and output samples of SH600586 are shown in Table 8.

The vector x in Table 8 represents the input vector, and y represents the output. This paper used a sliding window method to obtain all training samples. The forecast results of each model are shown in Figure 7.

From the analysis of the prediction results in Figure 7, it is concluded that both SZ000573, SZ002770 and SH600586 with a small change in stock trend or SZ000796, SZ002567 and SZ002607 with a large change in stock trend predicted by the proposed model are very close to the true value. However, there are some errors between the predicted results of the comparison models and the true values. The specific model performance analysis is given in the error analysis section.

6.8.   Error analysis

The forecast error analysis in this paper used mean absolute error (MAE), root mean square error (RMSE) and mean absolute percentage error (MAPE) as evaluation indicators.

where $\hat x\left( i \right)$ is the predicted value, and $x\left( i \right)$ is the true value.

From the analysis of Tables 9–14, it can be seen that when the model proposed in this paper predicted SZ000573, SZ002770 and SH600586 with small changes in trend, the error between the predicted value and the true value is on average below 0.21. When predicted SZ000796, SZ002567 and SZ002607 with large changes in trend, the average error between the predicted value and the true value is less than 0.78. The performance is stable.

When FGMC (1, m)-LLE-BP predicted SZ000573, SZ002770, and SH600586 with small changes in trend, the error between the predicted value and the true value is on average below 0.66. When predicted SZ000796, SZ002567 and SZ002607 with large changes in trend, the average errors are 1.01, 0.91 and 3.36 respectively. The performance is not stable.

When LSTM neural network model predicted the SZ000573, SZ002770, and SH600586 with small trend change, the average errors are 0.37, 0.16 and 0.49 respectively. When predicted the SZ000796, SZ002567 and SZ002607 with large changes in trend, the average errors are 1.29, 0.51 and 3.23 respectively. This shows that LSTMNN has volatility in stocks with large changes in predicted trends and stocks with small changes in predicted trends.

When BP neural network model predicted the SZ000573, SZ002770, and SH600586 with small trend change, the average errors are 0.31, 0.19 and 0. 90 respectively. When predicted the SZ000796, SZ002567 and SZ002607 with large changes in trend, the average errors are 1.53, 0.65 and 3.79 respectively.

When ARIMA-BP predicted the SZ000573, SZ002770, and SH600586 with small trend change, the average errors are 0.21, 0.14 and 0. 48 respectively. When predicted the SZ000796, SZ002567 and SZ002607 with large changes in trend, the average errors are 1.72, 0.91 and 1. 90 respectively.

Combining the analysis of three error indicators, the model proposed in this paper performs best, followed by ARIMA-BP and FGMC (1, m)-LLE-BP. However, neither the LSTM nor the BP neural network model of a single prediction model can achieve the desired effect. This shows that when predicting time series with linear and non-linear factors, the combined model can achieve good results.

From Figure 8, it can be seen that the model proposed in this paper has stable and excellent performance for 6 stocks with different trends, while the performance of the other 4 comparison models is only good in some stocks. Combining the analysis of three evaluation indicators, the model proposed in this paper performed best on the MAE, RMSE, and MAPE indicators when predicted SZ002607, SZ002567, SZ000796, SH600586, SZ002770 and SZ000573. Besides, the model proposed in this paper performs better than FGMC (1, m)-LLE-BP, which shows that the HFGMC (1, m) has a certain improvement in the accuracy of prediction.

6.9.   Limitations and future works

In the discussion in Section 6.2, it is pointed out that a model that deals with high-dimensional data requires certain prior data. However, it is very difficult to obtain real-time a priori data for forecasting objects with very high real-time requirements such as stock data. Generally, the high-dimensional grey model obtains the prior data by staggering the time axis of the prediction object (closing price) and prior data (indicators) (as shown in Figure 9). Therefore, the practice of staggering the time axis determines that the model cannot be used for long-term forecasting. When the time axis is staggered longer, the error may increase significantly. In addition to staggering the time axis, the use of sliding time windows to obtain prior data is also widely used in deep learning time series forecasting ^[2,3,30] (as shown in Table 8). In addition, the price of individual stocks (the experimental case of this paper) is different from stock indexes (Dow Jones, Nasdaq, Shanghai Composite Index, etc.). The price of individual stocks is usually predicted for about 10 days or 10 units of time ^[4,30,35], and the stock index reflects the general laws of the market, so long-term forecasts are feasible ^{[36,37,38,39]}. However, whether it is individual stock prices or stock indexes, the problem of prior data cannot be solved when using high-dimensional models. This is a limitation of high-dimensional models.

Since directly staggering the time axis may lose the information of various indicators in the time axis, the next step can consider using the maximum entropy model to integrate all the information in the time axis as the prior information of the prediction object. Through the above work, the accuracy of the model may be more accurate.

7.   Conclusions

Because the stock price is affected by many factors, the stock market can be regarded as a non-linear dynamic and grey system under the action of many factors. It is very difficult to predict stock price. There are many time series forecasting methods, but many methods do not consider that the stock market has long-term memory.

According to the characteristics of the stock market, the HFGMC (1, m)-LLE-BP is proposed. The HFGMC (1, m) is based on the FGMC (1, m) model, and combined the Hausdorff fractional derivative and Newton-Cotes formula. The HFGMC (1, m) model is used to solve the linear component in the prediction sequence. Then, the BP neural network is used to extract the non-linear component in the prediction sequence. Besides, in order to compress the redundant information in high-dimensional data, the HFGMC (1, m) model combined the locally linear embedding algorithm to reduce the dimensionality of factors that affect the stock price. Finally, the linear component and the non-linear component are added to obtain the combined prediction sequence. In the experimental part, the model proposed in this paper predicts the closing price of 6 stocks with different trends in the next 10 days. Based on the above experimental results and discussion, the following conclusions are drawn:

(i) FGMC (1, m) model has some shortcomings. Hausdorff fractional derivative and Newton-Cotes formula can improve FGMC (1, m) model. The HFGMC (1, m) model improves the prediction accuracy of time series.

(ii) Because the stock price is affected by a variety of linear and nonlinear factors, a single time prediction model cannot accurately predict the stock price. A combined forecasting model has advantages in forecasting time series affected by complex factors.

(iii) The model proposed in this paper is very stable and accurate in predicting the stock price with different historical trends. Results show that the model proposed in this paper can be widely used in the prediction of simple or complex time series.

Acknowledgments

This work is supported by National Nature Science Foundation under Grant 61862062, 61104035.

Conflict of interest

The authors declared that they have no conflict of interest.