Time series models play very important roles in many business decisions. In the current big data era, all walks of life are faced with the problem of modeling and time sequence prediction. For example, the e-commerce platform needs to predict the future sales of all commodities; in the pre-sales industry, both online and offline pre-sales require significant time series forecasting. These data are non-linearly correlated in time series, with most of them affected by product promotion, inventory situation and market competition among other, which may lead to outliers in time series. Meanwhile, sales during the holiday promotion period are relatively volatile and flat, resulting in an asymmetry of yield fluctuation and the rate of return usually does not follow the normal distribution, showing skewness and peak thick tail. Therefore, how could such data be modeled and predicted remains an open question.
Let’s start with classic time series models, such as the Autoregressive Integrated Moving Average model-Generalized Autoregressive Conditional Heteroskedast (ARIMA-GARCH) model  and normal Asymmetric Power Autoregressive Conditional Heteroskedasticity (APARCH) model  based on the Autoregressive Integrated Moving Average (ARIMA) model  . On the one hand, although the model can solve the heteroscedastic effect in the residual sequence, it provides no solution for singular information when data are applied to the Generalized Autoregressive Conditional Heteroskedast (GARCH) model  . On the other hand, the classical time series model in parameter estimation is usually based on the assumption of normal distribution, which does not fit well the distribution of volatility in practical applications. The explosive growth of new algorithm development makes this issue even more worthy of attention.
Therefore, in this study, an improved GARCH model, termed the Pro-GARCH model, was proposed to solve the problem of singular information in the data. The improved method is described below. First, in the GARCH (p, q) model, the rank of the Hessian matrix H is defined as . We performed QR decomposition on the Hessian matrix H, i.e. ; R (QR) represents the rank of the matrix after QR decomposition. Application of data to the model ( ) leads to the generation of singular information, and the rank of the matrix after QR decomposition is now increased by one, i.e. , so that , which allows to solve the problem of singular information. Since the rank of the matrix is changed after QR decomposition, the estimated value of a given parameter is not affected. Secondly, because the matrix is singular, the inverse matrix of the Hessian matrix was obtained by determining the generalized inverse of the matrix, yielding the Pro-GARCH model. Furthermore, a skewed-t APARCH model  based on the ARIMA model  was proposed. After assessing JD sales data, the results showed that the model was superior to the classical time series model in the accuracy of parameter estimation and prediction, and could more accurately describe the skewness problem in the sequence.
The remainder of the article is as follows. In the second part, the definitions of Pro-GARCH model and skewed-t APARCH model based on the ARIMA model will be provided. In the third part, the modeling process of skewed-t APARCH model based on the ARIMA model will be described. In the fourth part, we applied the novel and traditional classical time series models to JD sales data, respectively, and compared the results. Empirical analysis showed that the model is better than other models.
2.1. The Pro-GARCH Model
The Pro-GARCH (p, q) model we proposed is:
where, is the deterministic information fitting, is independent and follows the standard normal distribution, and , , and . The improved GARCH (p, q) model can also be rewritten as the ARMA (p, q) model for , i.e.
where,  .
With an outlier in the data, the actual sequence is not , but an observation sequence , defined as
where is the indicator function, and denote the magnitude and dynamic model of the outlier effect, respectively.
2.2. The Partial T-APARCH Model Based on the ARIMA Model
The skewed-t APARCH model  based on the ARIMA model  can be defined as:
where, , is an autoregressive coefficient polynomial for a stationary reversible ARIMA (p, q) model; is the moving smoothing coefficient polynomial for the stationary reversible ARIMA (p, q) model; is the conditional mean, represents a distribution with mean and variance of 0 and 1, respectively; ; is independent, identically distributed and follows the skewed distribution. The purpose of the power function in Equation (4) is to improve the transformation of model fitting.
1) We built the ARIMA model  as follows. a) Data preprocessing. First, the Pro-GARCH model and the improved outlier detection method  were used to detect IO type outliers of the data. Then, the iterative method  was used to process the outliers and generate new data; Secondly, the stability and pure randomness of the new time series data were tested. b) Model identification. After calculating the sample autocorrelation coefficient and partial correlation coefficient, the appropriate ARIMA model was selected to fit the observation sequence. c) Model prediction and diagnosis. The established model was used to predict future trends of time series values and evaluate the model by analyzing whether parameter estimates are significant and the residual is a white noise sequence.
2) Autocorrelation and heteroscedasticity test for the residual sequence of the ARIMA model were performed by a statistical method using the Portmanteau Q test and the LM test  .
3) Model identification. After constructing the ARIMA model, the residual sequence was modeled by APARCH (1, 1), selecting skewed-t distribution.
4) Model diagnosis. Whether the residual sequence was a white noise sequence and the parameter estimation significant was assessed.
5) Model prediction. The final fitted model was obtained and evaluated.
In this section, the construction process of the skewed-t Asymmetric Power Auto-regressive conditional heteroskedasticity (APARCH) model based on the ARIMA model is introduced in detail. Compared with the ARIMA-GARCH and normal APARCH models, respectively, based on the ARIMA model, the validity of the skewed-t APARCH model based on the ARIMA model was demonstrated. All simulations in this paper were performed in R.
4.1. Data Preprocessing
We analyzed the sales data of Jingdong. Figure 1 shows the presence of outliers in the data. Therefore, we first used the Pro-GARCH model and the improved outlier detection method to process outliers and generated new data. This result was satisfactory. Among them, the data obtained after processing the abnormal values are shown in Figure 2 and were recorded as .
4.2. Model Establishment
The pure randomness test results showed that the P value of the LB test statistic was very low under the first-order to sixth-order delay (see Table 1). Therefore, we determined that the sequence belonged to a non-white noise sequence and could model the data.
The stability of the sequence was verified by the timing diagram method. As shown in Figure 2, the sequence was non-stationary. After using the
Figure 1. Data with outliers.
Figure 2. Data after outlier processing.
Table 1. P value of LB test statistics.
ARIMA (p, d, q) model to fit the sequence , the first-order difference of the sequence is shown in Figure 3. As shown in Figure 3, the sequence after the first-order difference was stationary, so in the ARIMA (p, d, q) model, the order of the difference was 1, i.e. d = 1.
Figure 4 shows the autocorrelation (ACF) and partial correlation (PACF) plots of the sequence . The system’s automatic ordering was compared with the ACF and PACF graphs. With p = 1 and q = 1, the model fitting was most reasonable; therefore, the ARIMA (p, d, q) model most suitable for this sequence was the ARIMA (1, 1, 1) model. Regarding the heteroscedasticity test, since the P values of the LM and Portmanteau Q tests were low (see Table 2), the residuals were heteroscedastic. Due to space constraints, only the first six P values were reported in this paper.
Figure 3. Data after first-order difference.
Figure 4. Autocorrelation and partial correlation plots for .
Table 2. P values for LM and portmanteau Q tests.
Therefore, this study selected the ARIMA (1, 1, 1)-GARCH (1, 1) and APARCH (1, 1) models based on the ARIMA (1, 1, 1) model for modeling the sequence under normal and partial t distributions. Table 3 provides the parameter estimation results of the new model.
4.3. Evaluation Criteria
In order to evaluate the accuracy of the model, mean error (MSE), mean absolute error (MAE), root mean square error (RMSE) and mean absolute percentage error (MAPE) were used. The smaller the variance of each loss function, the
Table 3. Model parameter estimation results.
Note: “***”, “**”, “*”, “.” indicate that the parameters are significant.
smaller the prediction error, and the more accurate the prediction. The calculation formula was as follows:
where, is the realized volatility, estimated using high frequency data, and is the predicted volatility at time t; N is the number of the performance evaluation data.
4.4. Prediction Results
Tables 4-6 provide the standardized residual test of the ARIMA (1, 1, 1)-GARCH (1, 1) and APARCH (1, 1) models based on the ARIMA (1, 1, 1) model under the assumption of normal and partial t distributions, respectively. Table 7 shows the prediction effect of the 20 steps of the model. Figure 5 provides the standardized residual QQ diagrams of ARIMA (1, 1, 1)-GARCH (1, 1) and APARCH (1, 1) models based on the ARIMA (1, 1, 1) model in normal and partial t distributions, respectively.
Four evaluation indicators were assessed, including parameter estimation saliency, standardized residual test, evaluation criteria and standardized residual QQ map. 1) parameter estimation of the skewed-t APARCH (1, 1) model based on the ARIMA model was more significant. 2) Considering the standardized
Figure 5. (a) The standardized residual QQ plot of ARIMA (1, 1, 1)-GARCH (1, 1) model; (b) The standardized residual QQ diagrams of APARCH (1, 1) models based on the ARIMA (1, 1, 1) model in normal distributions; (c) the standardized residual QQ diagrams of APARCH (1, 1) models based on the ARIMA (1, 1, 1) model in partial t distributions.
Table 4. Standardized residual test of the ARIMA (1, 1, 1)-GARCH (1, 1) model.
Table 5. Standardized residual test of the normal APARCH model based on the ARIMA model.
Table 6. Standardized residual test of the skewed-t APARCH model based on the ARIMA model.
Table 7. Analysis of the prediction effect of the model.
residuals test, all three models completely eliminated the ARCH effect and the correlation between sequences, so they could not be rejected. 3) The skewed-t APARCH (1, 1) model based on the ARIMA model was slightly higher in accuracy compared with the other two models. 4) The skewed-t distribution had a better fitting effect in the standardized residual map, indicating that the influence of introducing bias on model fitting was significant. Therefore, the skewed-t APARCH (1, 1) model based on the ARIMA model had a better prediction ability. Therefore, the final prediction model was as follows:
where, the skewness is 0.925  , the model coefficient is greater than 0 and satisfies . Equations (9) and (10) are the mean and variance equations, respectively.
Further, the residual, predicted confidence interval and 95% confidence interval of the partial-t APARCH (1, 1) model based on the ARIMA (1, 1, 1) model are depicted in Figure 6. The residual was almost completely within the confidence interval, indicating that model prediction was more accurate; the volatility of the model is presented in Figure 7.
First, the Pro-GARCH model solves the singular information problem in data. Secondly, as shown in Figure 5, the skewed-t APARCH model based on the ARIMA model could better capture the peak thick tail, skewness and leverage effect in the sequence. Finally, Table 3 and Table 7 show that the model is superior to the ARIMA-GARCH and APARCH models based on the ARIMA model under the assumption of normal distribution in the significance of parameter estimation and accuracy of prediction, respectively.
Figure 6. Residual, predictive confidence interval and 95% confidence interval plots.
Figure 7. Model volatility.
In this study, the Pro-GARCH model and the improved outlier detection method were used, and the iterative method was used to process outliers in the time series to obtain a new time series. The ARIMA-GARCH and normal APARCH models based on the ARIMA model, and the skewed-t APARCH model based on the ARIMA model were compared. Some concluding observations can be summarized as follows:
1) Using the Pro-GARCH model and the improved outlier detection method to process data and selecting absolute deviation of the median (MAD) as a robust estimation of the standard deviation of the model, the location of outliers could be found most accurately;
2) Compared with the ARIMA-GARCH and normal APARCH models based on the ARIMA model, the skewed-t APARCH model based on the ARIMA model could better capture the spikes and thick tails, skewness and leverage effects, and the model had elevated prediction ability;
3) No prediction method could stand out in any time series. Although the skewed-t APARCH model based on the ARIMA model showed good predictive power, it did not achieve the expected results, and there were certain losses; this model is not flexible and cannot be applied to multiple products simultaneously. This is a huge challenge for time series modelers and requires further research.
This work is supported by the National Natural Science Foundation of China (11801294).