A robust maximum correntropy forecasting model for time series with outliers

The regression model

According to the linear theory, the auto-regression (AR) model of the series $\hat{Y}=\{{\hat{y}}_1,{\hat{y}}_2,\cdots {\hat{y}}_M\}$is shown in Eq. (1)

(1) $${\hat{y}}_t={\beta }_1{y}_{t-1}+{\beta }_2{y}_{t-2}+\cdots +{\beta }_N{y}_{t-N}+e_t$$where ${\hat{y}}_t,t=1,2,\cdots ,M$ is the forecasting load at the current moment t, M the number of group, Y_t−i the actual load at the past moment t − i, β_i the regression parameter, e_t the error, and N the regression order. Therefore, it can be expressed to Eq. (2)

(2) $$\left[\begin{array}{c}{\hat{y}}_1\\ {\hat{y}}_2\\ ⋮\\ {\hat{y}}_M\end{array}\right]={\beta }_1\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{M-1}\end{array}\right]+{\beta }_2\left[\begin{array}{c}{y}_{-1}\\ y_0\\ ⋮\\ y_{M-2}\end{array}\right]+\cdots +{\beta }_N\left[\begin{array}{c}{y}_{1-N}\\ y_{2-N}\\ ⋮\\ y_{M-N}\end{array}\right]+\left[\begin{array}{c}{e}_1\\ e_2\\ ⋮\\ e_M\end{array}\right]$$

The least square algorithm is adopted to minimize the sum of the squares of errors (Vovk et al., 2019; Midiliç, 2020). Let $\mathit{\Phi }=[{\beta }_1,{\beta }_2,\cdots ,{\beta }_N]$, $Z=\{z_1,z_2,\cdots ,z_M\}$ and $z_t=[y_{t-1},y_{t-2},\cdots y_{t-N}{]}^T$. Equation (3) is the objective function

(3) $$\underset{\mathit{\Phi }}{min}(\hat{Y}-\mathit{\Phi }Z{)}^2$$

Furthermore, Eq. (3) can also be written to Eq. (4)

(4) $$\underset{\mathit{\Phi }}{min}\sum _{t=1}^M{({\hat{y}}_t-\mathit{\Phi }{z}_t)}^2$$

Actually, some factors, such as the machine failure, human errors, can cause the outlier point in the recording process. Consequently, the accuracy can be reduced by the above method. In this work, we develop a robust regression model.

MCAR forecasting model

MCAR forecasting model

In the regression model, the mean square error of $\hat{Y}$ at the time t and the historical data $\mathit{\Phi }Z$ is the quadratic function of the convex curve along a straight line $\hat{Y}=\mathit{\Phi }Z$. However, the value away from $\hat{Y}=\mathit{\Phi }Z$ will increase mean error of samples, and make the regression parameters have larger error in MSE.

Traditional mean square error is measured in the way of global similarity. All samples in the joint space contribute significantly to the similarity. However, the correntropy is measured in a local way. Due to locality, the value of the similarity is determined by the kernel function along the line of $\hat{Y}=\mathit{\Phi }Z$ for random variables. Assume that $X=\hat{Y}-\mathit{\Phi }Z$, the Gaussian kernel function of X is shown in the Fig. 1.

From Fig. 1, we can see that if X = 0, Gaussian kernel $k_{\sigma }(\hat{Y}-\mathit{\Phi }Z)$ is maximized, and the residual at the origin is zero. Here, the maximum entropy is a kind of adaptive loss function, known as the maximal entropy criterion. It is suitable for the situation with non-Gaussian and large outlier value (Santamaria, Pokharel & Principe, 2006; Bessa, Miranda & Gama, 2009; Chen & Principe, 2012; He et al., 2011). In this work, we establish a robust multidimensional regression model based on the maximal entropy (MC).

Assume that the number of groups and order number of the model are M and N, respectively. For the convenience of calculation, the equality constraint ${\mathit{\Phi }}^T\mathit{\Phi }=1$ is introduced. Thereby, Eq. (9) is the constraint problem on MC

(9) $$\begin{array}{rl}& \underset{\mathrm{\Phi }}{max}\sum _{t=1}^M{k}_{\sigma }({\hat{y}}_t-\mathrm{\Phi }{z}_t)\\ & \mathrm{s}.\mathrm{t}.{\mathit{\Phi }}^T\mathit{\Phi }=1\end{array}$$

The optimization problem of Eq. (9) is nonlinear and non-convex, and cannot be solved directly. The conjugate convex function is introduced to solve the semi-quadratic form. Here, the auxiliary variables is introduced, and accordingly it can be simplified to Eq. (10)

(10) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}\sum _{t=1}^M{\omega }_t({\hat{y}}_t-\mathit{\Phi }{z}_t{)}^2\\ & \mathrm{s}.\mathrm{t}.{\mathit{\Phi }}^T\mathit{\Phi }=1\end{array}$$

Actually, the weighted function can reduce the large error term and the adverse effect of the outlier on the optimization result. We define the matrix $R=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(w)$, where $w=[{\omega }_1,{\omega }_2,\cdots ,{\omega }_M]$, and hence Eq. (10) is equivalent to Eq. (11)

(11) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}R{\Vert \hat{Y}-\mathit{\Phi }Z\Vert }^2\\ & \mathrm{s}.\mathrm{t}.{\mathit{\Phi }}^T\mathit{\Phi }=1\end{array}$$

The solution of Eq. (11) is a non-convex quadratic programming problem, which is difficult to solve. When the initial condition of $R(0)=diag(1)$, $\mathit{\Phi }$ becomes the optimal solution, and the above Eq. (11) is transformed into a homogeneous constraint programming of Eq. (12)

(12) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}{\Vert \gamma \hat{Y}-\mathit{\Phi }Z\Vert }^2\\ & \mathrm{s}.\mathrm{t}.{\gamma }^2=1,{\mathit{\Phi }}^T\mathit{\Phi }=1\end{array}$$

It is equivalent to Eq. (13)

(13) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}[\begin{array}{cc}{\mathit{\Phi }}^T& \gamma ]\end{array}\left[\begin{array}{cc}{Z}^{T}Z& -Z^T\hat{Y}\\ -{\hat{Y}}^{T}Z& {\Vert Z\Vert }^2\end{array}\right]\left[\begin{array}{c}{\mathit{\Phi }}^T\\ \gamma \end{array}\right]\\ & \mathrm{s}.\mathrm{t}.{\gamma }^2=1,{\mathit{\Phi }}^T\mathit{\Phi }=1\end{array}$$

Let $\xi =[{\mathit{\Phi }}^T\gamma {]}^T,B=\left[\begin{array}{cc}{I}_{N\times N}& 0\\ 0& 0\end{array}\right]$, and $C=\left[\begin{array}{cc}{Z}^{T}Z& -Z^T\hat{Y}\\ -{{\hat{Y}}_{}}^{T}Z& {\Vert Z\Vert }^2\end{array}\right]$. Equation (13) can simplified to Eq. (14)

(14) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}{\xi }^{T}C\xi \\ & \mathrm{s}.\mathrm{t}.{\xi }^{T}B\xi =1\end{array}$$

Using the semi-definite relaxation (SDR) method, the objective function and constraint conditions in Eq. (14) are equivalent to Eq. (15)

(15) $$\begin{array}{rl}& {\xi }^{T}C\xi =\mathrm{T}\mathrm{r}\{{\xi }^{T}C\xi \}=\mathrm{T}\mathrm{r}\{C\xi {\xi }^T\}\\ & {\xi }^{T}B\xi =\mathrm{T}\mathrm{r}\{{\xi }^{T}B\xi \}=\mathrm{T}\mathrm{r}\{B\xi {\xi }^T\}\end{array}$$where Tr{} is the trace of the matrix. In Eq. (16), we define the matrix

(16) $$\mathit{\omega }=\xi {\xi }^T$$where $\mathit{\omega }$ is a symmetric positive semi-definite (PSD) matrix with a rank of 1. The resulting semi-definite relaxation optimization constraint is shown in Eq. (17)

(17) $$\begin{array}{rl}& \underset{\mathit{\Phi }}{min}\mathrm{T}\mathrm{r}\{C\mathit{\omega }\}\\ & \mathrm{s}.\mathrm{t}.\mathrm{T}\mathrm{r}\{B\mathit{\omega }\}=1,\mathit{\omega }\ge 0\end{array}$$

Furthermore, let Eq. (18) is the eigen-decomposition of the matrix $\mathit{\omega }$

(18) $$\mathit{\Phi }=VR{V}^T$$where $V=[v_1,v_2,\dots v_M]$ is the eigenvector of $\mathit{\omega }$, and $R=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(r_1,r_2\cdots ,r_M)$ is the corresponding eigenvalue. $\mathit{\omega }(1)=r_1{v}_1{v_1}^T$ is closest to $\mathit{\omega }$ when the rank is 1, and so $\xi$ is estimated in Eq. (19)

(19) $$\xi =\sqrt{r_1{v}_1}$$

Finally, the parameter $\mathit{\Phi }$ is the first N value of $\xi$.

In summary, the steps of the robust MCAR forecasting method are as follows:

Step (1): Set the maximum number of iteration, and initialize $R(0)=diag(1)$;

Step (2): Solve the optimization constraint problem in Eq. (17) and $\xi$ by Eq. (19);

Step (3): According to the Silverman specification: $\sigma =1.06\times min({\sigma }_e,R/1.34)\times I^{-0.2}$ ( ${\sigma }_e$ is the standard deviation of $\hat{Y}-\mathit{\Phi }Z$, and R is the quartile difference) get solution of ${\omega }_t$;

Step (4): Perform Steps (2) and Step (3) until the termination condition is reached, and find the regression parameters $\mathit{\Phi }$;

Step (5): Forecast the data at the next moment based on the parameter $\mathit{\Phi }$ in Step (4).

It is worth pointing out that, unlike previous methods (Bashir & El-Hawary, 2009), the proposed MCAR model does not need to judge whether there is outliers in the data set, but directly constructs and trains the model based on the data set, since the MCAR model can automatically reduce the impact of outliers by the maximal entropy.

Comparison of MCAR and regression models

In this section, the performance of the proposed MCAR model is compared with the differential autoregressive (MDAR) model and the AR model. Here, to guarantee a fair comparison, parameters of three methods are set to the same value as shown in Table 1. The training data sets are the same, which are taken from the January data set. Firstly, the performance of three methods on normal data (no outlier) is verified. Results of the relative errors (Re) of the 50 different forecasting points are illustrated in Fig. 3. It can be seen that, the relative error of the AR and the MDAR model have a large error at the forecasting point 7, 29 and 34, due to the sudden increase or decrease of the actual data.

Table 2 shows the corresponding performance indexes. Here, the MAPE of MCAR is 4.74%, while that of the AR and MDAR models are greater than 7%. Meanwhile, the other performance indexes of MCAR are relatively small, indicating that the proposed MCAR model is superior to the other regression models.

Next, the robustness of the MCAR model on the outlier is further verified. Fig. 4 shows the comparison of forecasting results from January 6 (day 1) to January 20 (day 15). It can be observed that, after the training sample with the outlier, the forecasting values of multiple points of AR and MDAR deviate greatly. Moreover, the results of MCAR are still good compared with the actual valuables, which mean that the proposed MCAR model in this work has higher robustness.

The comparisons of the relative errors (Re) between the forecasting result and the actual one are shown in Fig. 5. It can be seen that, AR and MDAR models have large forecasting errors. The relative error of the MCAR model remains relatively low, indicating that it is less affected by the outlier. From the performance indexes in Table 3, it can be seen that the MAPE of MCAR model is 4.88%, which is significantly smaller than that of the AR and MDAR models. In addition, compared with results with no outlier in Table 2, performance indexes did not increase significantly after the outlier data was added. This shows that, the MCAR model proposed in this work is robust to the outlier and can improve the forecasting accuracy.

Parameters selection of MCAR

The influences of parameter kernel width σ, model order N and group number M on the accuracy are discussed.

The influences of parameters of the Gaussian kernel on the forecasting performance are analyzed. Actually, for different data sets, the width of Gaussian kernel is different. Here, the Silverman rule is used to determinate the kernel width. The performance indexes are calculated using the forecasting valuables obtained from each kernel width, compared with the Silverman rule. It can be seen from Table 4 that, RMSE, MAPE and MAE of Silverman rule are smaller than the corresponding valuables with other kernel width. Here, the model considers the standard deviation and quartile potential difference of error sequence synthetically, realizing the restriction of increasing the correntropy to the larger error estimate value, and improves the robustness of the correntropy measurement method to the outlier value.

Then, the AICc criterion is applied to optimize the model order in this work (Chang et al., 2018). Equation (23) is the AICc criterion

(23) $$AICc=T\ln RSSN+2N+{\displaystyle \frac{2N(N+1)}{T-N-1}}$$where T is the number of samples, and N the model order. $RSSN={\Vert Y-\mathit{\Phi }Z\Vert }^2$ is the forecasting error. Here, when T is small, the constraint ability of the AICc criterion on the number of parameters is strengthened, and so it is applicable to the case that the number of samples is small in this work.

Table 5 shows the performance indexes of the normal data and the data with the outlier with different model orders. Here, the number of groups is set to 4. It can be seen that, the performance indexes are smallest when the order N = 4, which indicates that the forecasting results are dependent on the model order. As shown in Table 6, the minimum value of the model order is 34.07 of AICc model. Meanwhile, the optimal model order is N = 4. This is consistent with the experimental result in Table 5.

Lastly, the number of groups in the MCAR model is also a variable parameter which can affect on the forecasting performance. Table 7 shows the performance indexes of MCAR model with different groups. Here, kernel width is optimized by the Silverman rule and the model order is fixed to 4. It can be seen that, whether or not the load data with the outlier, the forecasting error is minimum when the number of groups is 4 in this experiment. Here, when the number of groups is small, the model cannot find the similarity of load series. However, with the further increase of the number of groups, the similarity of the load series decreases, and the distribution difference of the time series becomes larger. It should be noted that, the optimal number of groups is different for different data sets. Thereby, in the actual modeling process, the number of groups needs to be selected according to the characteristics of the dataset.