# Electric Load Forecast Using Combined Models with HP Filter-SARIMA and ARMAX Optimized by Regression Analysis Algorithm

Authors: Cui  Herui; Peng Xu; Mu Yupei
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Electric load in summer has a significant cyclical trend with temperature effects. In general, the parameters of the SARIMA and the SMA turn out to be nonsignificant in most cases. To address this issue, the hybrid time series model is utilized to extract the spectrum sequences with different frequencies. The original electric load series are first decomposed into the trend sequence “G” and the cycle sequence “C.” After that, a revised ARMAX model is proposed to deal with the two divided sequences. Finally, the combined models are tested by case study. The case study on electric load forecast in one city from China shows that the proposed model outperforms other four comparative models in terms of prediction accuracy. It proves that the combined model proposed by the authors is more accurate than those based on a single forecasting method.


---

## Body

## 1. Introduction

Load forecast is mainly used to predict power load in the next few days [1, 2], which plays an important role in the modern electricity Demand Side Management (DSM). Accurate forecast on short-term electric load is a crucial element to the dynamic operations of Advanced Electricity Demand Side Management (EDSM) and Advanced Power Information Systems (APIS) in improving the efficiency and the safety of power grid. To a certain extent, medium-term power load is affected by seasonal factors, summer temperature, and consumption peaks due to unexpected cases. In general, high temperature in summer may form a high “air conditioning load.” With the rise of temperature, air conditioning load increases as well. Temperature data fluctuate widely in summer. Electric load also changes a lot with temperature. Theoretically, it is difficult to forecast mutable data, such as summer air conditioning load. To address this problem, this paper proposes the combined model with HP Filter-SARIMA and ARMAX model optimized by regression analysis algorithm.Over the past decade, many forecasting methods have been put forward, such as time series, gray model, SVM, and artificial neural networks (ANN) [3]. Azadeh et al. explored seasonal fluctuation and nonlinearity in forecasting based on the fuzzy system and data mining techniques to analyze monthly electricity demand in Iran. This model is established on well developed statistical theories to show explicit relationships between input data and outputs. However, the SARIMA model does not perform well when electric power deviates greatly from the normal weekly pattern. It cannot react under abnormal load conditions before the flow deviation is detected [4]. Hamzacebi and Es predicted annual electricity consumption by the optimization GM(1,1) model, which is more effective in analyzing the long-term trend but has less effect in seasonal variation [5]. Zhu et al. used the HP filter to decompose the GDP sequence into the tendency item and the cyclical item [6]. He et al. also took advantage of the HP filter in energy price analysis to study the synchronization between the markets home and abroad [7]. Zhang et al. employed four improved adaptive coefficient approaches optimized by particle swarm optimization (PSO) to forecast daily mean wind speed, the simulated results of which showed that the PSO obtained an observable improvement in forecasting performance [8]. Guo et al. proposed a modified EMD-FNN model by combining empirical mode decomposition (EMD) with the ensemble learning paradigm of feedforward neural network (FNN), which had better accuracy than that based on the basic FNN and unmodified EMD-FNN [9].The multivariate ARIMAX model is hypothesized to improve the one-step-ahead forecasting accuracy of the univariate SARIMA model. Bierens and Broersma used the ARMAX model to study the relation between unemployment and interest rate. They found that the relationship is not confined only to the Netherlands; but it holds for USA, Canada, Japan, Germany, UK, and France [10]. Bordignon et al. analyzed combined versus individual forecasts for British electricity price prediction. It is found that combined forecasts are more accurate than or at least equivalent to individual ones [11]. Yan and Chowdhury adopted a hybrid midterm forecasting model based on the combination of both least squares support vector machine (LSSVM) and autoregressive moving average with external input (ARMAX) modules to forecast electricity market clearing price (MCP). It is shown that the proposed hybrid model can improve forecasting accuracy compared with the forecasting model using a single LSSVM [12]. Wang et al. proposed a two-stage model in estimating value-at-risk (VaR) based on ARMAX-GARCHSK and extreme value theory (EVT). It is shown in the empirical analysis that the ARMAX-GARCHSK-EVT model can rapidly reflect the most recent and relevant change of electricity prices, with accurate forecasts of VaR at all confidence levels, thereby presenting better dynamic characteristics [13]. Yang et al. proposed a new evolutionary programming (EP) approach to identify the autoregressive moving average with exogenous variable (ARMAX) model for hourly load demand forecasts from one day to one week ahead. The developed EP based load forecasting algorithm is verified by different types of data for Taiwan Power (Taipower) system and substation load as well as temperature values [14]. Huang et al. proposed a new particle swarm optimization (PSO) approach to identify the autoregressive moving average with exogenous variable (ARMAX) model for load forecasts. It is indicated by the testing results that the proposed PSO has the characteristic of high-quality solution, superior convergence, and shorter computation time [15].Wangdi et al. adapted ARIMAX model to determine predictors of malaria in the coming month. ARIMAX model is an extension of ARIMA modeling in an attempt to predict the malaria cases using the climatic factors and the number of cases in the previous month. The predictors in the model include the number of cases in the previous month, mean maximum and minimum temperature, relative humidity, and rainfall lagged in a month. It is shown by test results that prediction accuracy has been greatly improved [16].The above forecasting methods have obvious effects in dealing with the cyclical and trend data. However, the generalization capability is generally weak. Traditional time series forecasting methods can be used to predict short-term load data [17]. The forecasting accuracy of the method is not as good as the combined models [18]. Autoregressive integrated moving average with exogenous variables (ARIMAX) model is mainly aimed at forecasting the data under the influence of external factors. However, there are only a few literatures discussing this model for load forecast. With respect to the above active research, the combined models of HP Filter-SARIMA-Revised ARMAX optimized by regression analysis are proposed here to forecast short-term electric load. They have strong application value in summer load forecasting field.


## 2. Forecasting Methods

### 2.1. Basic Concepts

Some important concepts are displayed below.Stationarity. Set xt as time series, m∈Z∗,t1,t2,…,tm∈T,∀τ∈Z,  ∃Ft1,t2,…,tm(x1,x2,…,xm)=Ft1+τ,t2+τ,…,tm+τ(x1,x2,…,xm) and named xt for strictly stationary time series(1)Ext=μ,var⁡xt=σ2,cov⁡xt,xt-1=γl.White Noise. Time series xt meet the condition of ①∀t∈T,∃EXt=μ and ②∀t,s∈T; then(2)γt,s=σ2,t=s0,t≠snamed xt for white noise sequence or is displayed below:(3)Xt~WNμ,σ2.Autocorrelation Function. Consider(4)ρk=∑t=k+1nwt-w-wt-k-w-∑t=1nwt-w-2.Cointegration Theory. Cointegration theory was put forward by Engle and Granger [19]. Forecasting models can be developed without requiring that all sequences are stationary. Regression model is as follows:(5)yt=β0+∑i=1kβixit+εt.Cointegration relationship exists between the independent variable sequence{x1},…,{xk} and response variable sequence {yt}, with stable regression residuals sequence. There is a strong correlation between electric load and temperature in summer, which is shown in Figure 4. The cointegration relationship is evident, which needs to be tested at the second section of the combined forecasting models.


### 2.2. Forecasting Models

#### 2.2.1. ARMA

Wold’s decomposition theorem (1938) is the foundation of time series analysis, which states that any stationary time series can be decomposed into the purely deterministic component and the purely stochastic component. The latter one can be expressed as a moving average time series. The theorem is expressed as(6)Zt=Ztd+∑i=0∞ψiet-i,where Zt(d) is purely deterministic component; ∑i=0∞ψiet-i is a moving average time series. ARMA model (autoregressive moving average model) requires that sequence itself is stable. Thus the nonstationary sequence data is hard to be managed.ARMA model is defined as(7)zt=ϕ0+ϕ1zt-1+⋯+ϕpzt-p+εt-θ1εt-1-⋯-θqεt-q,ϕp≠0,θq≠0,Eεt=0,Varεt=σε2,Eεtεs=0,s≠t,Ezsεt=0,∀s<t.By introducing delay operatorB, it can also be presented as Φ(B)zt=Θ(B)εt.In terms of structure,ARIMAp,d,q models are the same as ARMAp,q models, where the time series has first been transformed by differencing, the order of which is specified by d. ARIMA model flowchart is as follows in Figure 1.Figure 1
ARIMA model flowchart.


#### 2.2.2. HP Filter

Hodrick-Prescott Filter. Hodrick and Prescott first put forward Hodrick-Prescott filter (HP filter) method in an analytical paper about economic cycle in postwar America [20]. The sequence is divided into two groups: the sequence G with long-term trend, denoted as G={g1,g2,…,gn}, and the sequence C with short-term volatility, denoted as C={c1,c2,…,cn}. The relationship is counted as (8)yt=gt+ct.In this paper, HP filter is applied to the non-seasonally adjusted series, and the original sequence is divided into two sequences with the significant spectral frequency. The model is more accurate, which weakens the interaction among factors like seasonality, trend, cycle, and so forth [21].


#### 2.2.3. ARMAX

ARMAX Model. Supposing that the response variable yt and the input variable sequences x1,…,xk are all stationary, then the regression model is established in response to the input variable sequences and response sequences [22](9)yt=μ+∑i=1kΘiBΦiBBlixit+εt,εt=ΘBΦBat,abbreviated  as  ARMAX  model.The ARMAX model is an improvement of the ARIMA model with explanatory exogenous variables(X). The model is a combination of a regression model with an ARIMA model, which includes the advantages of both models. In the actual modeling process, a combined “HP Filter-SARIMA-Revised ARMAX” model is proposed to forecast the short-term electric load. The specific process is displayed below.


#### 2.2.4. Grey Prediction Model

The grey prediction processes are as follows [23].Step 1.
Consider first-order accumulated generating (1-AGO).
VariableX(0) is the original nonnegative sequence(10)X0=X01,X02,…,X0n.
First-order accumulated generating sequence ofX(0) is(11)X1=X11,X12,…,X1n,wherex1k=∑i=1kx0i,k=1,2,…,n.Step 2.
Consider operating quasi-smooth test and quasi-exponential law test toX(0):(12)ρk=x0kx1k-1,k=2,3,…,n.
Ifρ(k)<1 and ρ(k)∈[0,ε](ε<0.5), ρ(k) shows a decreasing trend. X(0) is quasi-smooth sequence, while X(1) has quasi-exponential law. Otherwise(13)x′0k=1n-k+1xk+xk+1+⋯+xn,k=1,2,…,n.Step 3.
Consider the establishment of first-order linear differential equations:(14)dx1dt+a^x1=u^.Step 4.
Consider building grey forecasting model of sequenceX(0):(15)x^0k+1=e-a^-1x0n-u^a^e-a^k,k=0,1,2,…,n.Step 5.
Test the model.


#### 2.2.5. MLP Neural Network

Multilayer perceptron neural network (MLP neural network) is one of the artificial neural networks which contains three processes, namely, training, testing, and validation [24]. Figure 2 is a multilayer neural network diagram(16)yk1=11+e-w1kTx-ak1,k=1,…,M1,y1=y11,…,yM11T,yk2=11+e-w2kTy1-ak2,k=1,…,M2,y2=y12,…,yM22T,⋮yout=Fx;W=wpTyp-1.Figure 2
MLP neural networks.Synaptic weight change rules for the neurons of the hidden layer are as follows:(17)wj1,it+1=wj1,it-c·∂Et∂wj1,i,∂Et∂wj1,i=2·Fxt,Wt-yt·∂∂wj1,i11+e-w1,iTxt-ai,∂∂wj1,i11+e-w1,iTxt-ai=e-w1,iTxt-ai1+e-w1,iTxt-ai2·∂∂wj1,i-w1,iTxt-ai,∂∂wj1,i-w1,iTxt-ai=-xjt,wj1,it+1=wj1,it-c·2·Fxt,Wt-yt·e-w1,iTxt-ai1+e-w1,iTxt-ai2·-xjt.Synaptic weight change rules for the output neuron are as follows:(18)wi2t+1=wi2t-c·∂Et∂wi2,∂Et∂wi2=2·Fxt,Wt-yt·11+e-w1,iTxt-ai.


#### 2.2.6. Regression Analysis

Regression analysis is used to forecast the value of one variable (dependent variable) based on other variables (independent variables).The simple linear regression model is displayed below:(19)y=β0+β1x1+β2x2+⋯+βkxk+ε,where dependent variable is denoted as y, independent variables are expressed as x1,x2,…,xk, and β0,β1,…,βk,ε are parameters.


## 2.1. Basic Concepts

Some important concepts are displayed below.Stationarity. Set xt as time series, m∈Z∗,t1,t2,…,tm∈T,∀τ∈Z,  ∃Ft1,t2,…,tm(x1,x2,…,xm)=Ft1+τ,t2+τ,…,tm+τ(x1,x2,…,xm) and named xt for strictly stationary time series(1)Ext=μ,var⁡xt=σ2,cov⁡xt,xt-1=γl.White Noise. Time series xt meet the condition of ①∀t∈T,∃EXt=μ and ②∀t,s∈T; then(2)γt,s=σ2,t=s0,t≠snamed xt for white noise sequence or is displayed below:(3)Xt~WNμ,σ2.Autocorrelation Function. Consider(4)ρk=∑t=k+1nwt-w-wt-k-w-∑t=1nwt-w-2.Cointegration Theory. Cointegration theory was put forward by Engle and Granger [19]. Forecasting models can be developed without requiring that all sequences are stationary. Regression model is as follows:(5)yt=β0+∑i=1kβixit+εt.Cointegration relationship exists between the independent variable sequence{x1},…,{xk} and response variable sequence {yt}, with stable regression residuals sequence. There is a strong correlation between electric load and temperature in summer, which is shown in Figure 4. The cointegration relationship is evident, which needs to be tested at the second section of the combined forecasting models.


## 2.2. Forecasting Models

### 2.2.1. ARMA

Wold’s decomposition theorem (1938) is the foundation of time series analysis, which states that any stationary time series can be decomposed into the purely deterministic component and the purely stochastic component. The latter one can be expressed as a moving average time series. The theorem is expressed as(6)Zt=Ztd+∑i=0∞ψiet-i,where Zt(d) is purely deterministic component; ∑i=0∞ψiet-i is a moving average time series. ARMA model (autoregressive moving average model) requires that sequence itself is stable. Thus the nonstationary sequence data is hard to be managed.ARMA model is defined as(7)zt=ϕ0+ϕ1zt-1+⋯+ϕpzt-p+εt-θ1εt-1-⋯-θqεt-q,ϕp≠0,θq≠0,Eεt=0,Varεt=σε2,Eεtεs=0,s≠t,Ezsεt=0,∀s<t.By introducing delay operatorB, it can also be presented as Φ(B)zt=Θ(B)εt.In terms of structure,ARIMAp,d,q models are the same as ARMAp,q models, where the time series has first been transformed by differencing, the order of which is specified by d. ARIMA model flowchart is as follows in Figure 1.Figure 1
ARIMA model flowchart.


### 2.2.2. HP Filter

Hodrick-Prescott Filter. Hodrick and Prescott first put forward Hodrick-Prescott filter (HP filter) method in an analytical paper about economic cycle in postwar America [20]. The sequence is divided into two groups: the sequence G with long-term trend, denoted as G={g1,g2,…,gn}, and the sequence C with short-term volatility, denoted as C={c1,c2,…,cn}. The relationship is counted as (8)yt=gt+ct.In this paper, HP filter is applied to the non-seasonally adjusted series, and the original sequence is divided into two sequences with the significant spectral frequency. The model is more accurate, which weakens the interaction among factors like seasonality, trend, cycle, and so forth [21].


### 2.2.3. ARMAX

ARMAX Model. Supposing that the response variable yt and the input variable sequences x1,…,xk are all stationary, then the regression model is established in response to the input variable sequences and response sequences [22](9)yt=μ+∑i=1kΘiBΦiBBlixit+εt,εt=ΘBΦBat,abbreviated  as  ARMAX  model.The ARMAX model is an improvement of the ARIMA model with explanatory exogenous variables(X). The model is a combination of a regression model with an ARIMA model, which includes the advantages of both models. In the actual modeling process, a combined “HP Filter-SARIMA-Revised ARMAX” model is proposed to forecast the short-term electric load. The specific process is displayed below.


### 2.2.4. Grey Prediction Model

The grey prediction processes are as follows [23].Step 1.
Consider first-order accumulated generating (1-AGO).
VariableX(0) is the original nonnegative sequence(10)X0=X01,X02,…,X0n.
First-order accumulated generating sequence ofX(0) is(11)X1=X11,X12,…,X1n,wherex1k=∑i=1kx0i,k=1,2,…,n.Step 2.
Consider operating quasi-smooth test and quasi-exponential law test toX(0):(12)ρk=x0kx1k-1,k=2,3,…,n.
Ifρ(k)<1 and ρ(k)∈[0,ε](ε<0.5), ρ(k) shows a decreasing trend. X(0) is quasi-smooth sequence, while X(1) has quasi-exponential law. Otherwise(13)x′0k=1n-k+1xk+xk+1+⋯+xn,k=1,2,…,n.Step 3.
Consider the establishment of first-order linear differential equations:(14)dx1dt+a^x1=u^.Step 4.
Consider building grey forecasting model of sequenceX(0):(15)x^0k+1=e-a^-1x0n-u^a^e-a^k,k=0,1,2,…,n.Step 5.
Test the model.


### 2.2.5. MLP Neural Network

Multilayer perceptron neural network (MLP neural network) is one of the artificial neural networks which contains three processes, namely, training, testing, and validation [24]. Figure 2 is a multilayer neural network diagram(16)yk1=11+e-w1kTx-ak1,k=1,…,M1,y1=y11,…,yM11T,yk2=11+e-w2kTy1-ak2,k=1,…,M2,y2=y12,…,yM22T,⋮yout=Fx;W=wpTyp-1.Figure 2
MLP neural networks.Synaptic weight change rules for the neurons of the hidden layer are as follows:(17)wj1,it+1=wj1,it-c·∂Et∂wj1,i,∂Et∂wj1,i=2·Fxt,Wt-yt·∂∂wj1,i11+e-w1,iTxt-ai,∂∂wj1,i11+e-w1,iTxt-ai=e-w1,iTxt-ai1+e-w1,iTxt-ai2·∂∂wj1,i-w1,iTxt-ai,∂∂wj1,i-w1,iTxt-ai=-xjt,wj1,it+1=wj1,it-c·2·Fxt,Wt-yt·e-w1,iTxt-ai1+e-w1,iTxt-ai2·-xjt.Synaptic weight change rules for the output neuron are as follows:(18)wi2t+1=wi2t-c·∂Et∂wi2,∂Et∂wi2=2·Fxt,Wt-yt·11+e-w1,iTxt-ai.


### 2.2.6. Regression Analysis

Regression analysis is used to forecast the value of one variable (dependent variable) based on other variables (independent variables).The simple linear regression model is displayed below:(19)y=β0+β1x1+β2x2+⋯+βkxk+ε,where dependent variable is denoted as y, independent variables are expressed as x1,x2,…,xk, and β0,β1,…,βk,ε are parameters.


## 2.2.1. ARMA

Wold’s decomposition theorem (1938) is the foundation of time series analysis, which states that any stationary time series can be decomposed into the purely deterministic component and the purely stochastic component. The latter one can be expressed as a moving average time series. The theorem is expressed as(6)Zt=Ztd+∑i=0∞ψiet-i,where Zt(d) is purely deterministic component; ∑i=0∞ψiet-i is a moving average time series. ARMA model (autoregressive moving average model) requires that sequence itself is stable. Thus the nonstationary sequence data is hard to be managed.ARMA model is defined as(7)zt=ϕ0+ϕ1zt-1+⋯+ϕpzt-p+εt-θ1εt-1-⋯-θqεt-q,ϕp≠0,θq≠0,Eεt=0,Varεt=σε2,Eεtεs=0,s≠t,Ezsεt=0,∀s<t.By introducing delay operatorB, it can also be presented as Φ(B)zt=Θ(B)εt.In terms of structure,ARIMAp,d,q models are the same as ARMAp,q models, where the time series has first been transformed by differencing, the order of which is specified by d. ARIMA model flowchart is as follows in Figure 1.Figure 1
ARIMA model flowchart.


## 2.2.2. HP Filter

Hodrick-Prescott Filter. Hodrick and Prescott first put forward Hodrick-Prescott filter (HP filter) method in an analytical paper about economic cycle in postwar America [20]. The sequence is divided into two groups: the sequence G with long-term trend, denoted as G={g1,g2,…,gn}, and the sequence C with short-term volatility, denoted as C={c1,c2,…,cn}. The relationship is counted as (8)yt=gt+ct.In this paper, HP filter is applied to the non-seasonally adjusted series, and the original sequence is divided into two sequences with the significant spectral frequency. The model is more accurate, which weakens the interaction among factors like seasonality, trend, cycle, and so forth [21].


## 2.2.3. ARMAX

ARMAX Model. Supposing that the response variable yt and the input variable sequences x1,…,xk are all stationary, then the regression model is established in response to the input variable sequences and response sequences [22](9)yt=μ+∑i=1kΘiBΦiBBlixit+εt,εt=ΘBΦBat,abbreviated  as  ARMAX  model.The ARMAX model is an improvement of the ARIMA model with explanatory exogenous variables(X). The model is a combination of a regression model with an ARIMA model, which includes the advantages of both models. In the actual modeling process, a combined “HP Filter-SARIMA-Revised ARMAX” model is proposed to forecast the short-term electric load. The specific process is displayed below.


## 2.2.4. Grey Prediction Model

The grey prediction processes are as follows [23].Step 1.
Consider first-order accumulated generating (1-AGO).
VariableX(0) is the original nonnegative sequence(10)X0=X01,X02,…,X0n.
First-order accumulated generating sequence ofX(0) is(11)X1=X11,X12,…,X1n,wherex1k=∑i=1kx0i,k=1,2,…,n.Step 2.
Consider operating quasi-smooth test and quasi-exponential law test toX(0):(12)ρk=x0kx1k-1,k=2,3,…,n.
Ifρ(k)<1 and ρ(k)∈[0,ε](ε<0.5), ρ(k) shows a decreasing trend. X(0) is quasi-smooth sequence, while X(1) has quasi-exponential law. Otherwise(13)x′0k=1n-k+1xk+xk+1+⋯+xn,k=1,2,…,n.Step 3.
Consider the establishment of first-order linear differential equations:(14)dx1dt+a^x1=u^.Step 4.
Consider building grey forecasting model of sequenceX(0):(15)x^0k+1=e-a^-1x0n-u^a^e-a^k,k=0,1,2,…,n.Step 5.
Test the model.


## 2.2.5. MLP Neural Network

Multilayer perceptron neural network (MLP neural network) is one of the artificial neural networks which contains three processes, namely, training, testing, and validation [24]. Figure 2 is a multilayer neural network diagram(16)yk1=11+e-w1kTx-ak1,k=1,…,M1,y1=y11,…,yM11T,yk2=11+e-w2kTy1-ak2,k=1,…,M2,y2=y12,…,yM22T,⋮yout=Fx;W=wpTyp-1.Figure 2
MLP neural networks.Synaptic weight change rules for the neurons of the hidden layer are as follows:(17)wj1,it+1=wj1,it-c·∂Et∂wj1,i,∂Et∂wj1,i=2·Fxt,Wt-yt·∂∂wj1,i11+e-w1,iTxt-ai,∂∂wj1,i11+e-w1,iTxt-ai=e-w1,iTxt-ai1+e-w1,iTxt-ai2·∂∂wj1,i-w1,iTxt-ai,∂∂wj1,i-w1,iTxt-ai=-xjt,wj1,it+1=wj1,it-c·2·Fxt,Wt-yt·e-w1,iTxt-ai1+e-w1,iTxt-ai2·-xjt.Synaptic weight change rules for the output neuron are as follows:(18)wi2t+1=wi2t-c·∂Et∂wi2,∂Et∂wi2=2·Fxt,Wt-yt·11+e-w1,iTxt-ai.


## 2.2.6. Regression Analysis

Regression analysis is used to forecast the value of one variable (dependent variable) based on other variables (independent variables).The simple linear regression model is displayed below:(19)y=β0+β1x1+β2x2+⋯+βkxk+ε,where dependent variable is denoted as y, independent variables are expressed as x1,x2,…,xk, and β0,β1,…,βk,ε are parameters.


## 3. Establishing Process of the Combined Models

### 3.1. Modeling Steps

Step 1.
The original sequenceyt or xt, defined as the superposition of the waves with different frequencies, can be divided into the sequence G and the sequence C [25], the separation process of which must satisfy the minimum loss function principle:(20)min⁡∑t=1nYtGC-YtG2+λ∑Yt+1G-YtG-YtG-Yt-1G2.Step 2.
The sequenceG is the function of time. The error term is denoted as ε1,ε2,…,εn. The method of ordinary least squares (OLS) is used in making polynomial curve fitting to minimize the sum of squared errors. Polynomial fitting can be represented as(21)yG=a1tn+a2tn-1+⋯+an-1t+an.
The sequences’ predictions can be got and defined asg^t.Step 3.
The stationary property of the sequenceC is tested. If it is stationary, the correlation analysis can be used on the sequence; If it is stationary, the correlation analysis can be used on the sequence; if otherwise, the difference transformation is conducted on the sequence until it is stationary [26, 27]. The stationary sequence is denoted as wc.Step 4.
Through the correlation analysis, the autocorrelation and moving average items can be acquired. The ARIMA model is presented as(22)yc=φp-1LΦP-1L1-L-d1-Ls-DθqLΘQLεc.Step 5.
Consider testing of the rationality of the model by the residual sequence.Step 6.
By ARIMA model, the predictions of the sequences can be denoted asc^t.Step 7.
The final prediction result can be got according to the principle of HP filter:(23)y^t=g^t+c^t.Step 8.
Consider exploring the correlation coefficient between the stationary sequences “y^t” and “X^t” to determine the structure of the improved ARMAX model. This step is an improved version for traditional ARMAX model. The revised ARMAX model can be calculated as follows: (24)yt∗=μ+∑i=1kΘiBΦiBBlixit∗+εt.Step 9.
Consider fitting residual sequenceεt:(25)εt=ΘBΦBat,where at is a zero mean white noise sequence.Based on above steps, the combined “HP Filter-SARIMA-Revised ARMAX” model can be applied in load forecasting process.


### 3.2. Modeling Flowchart

The forecasting processes of the combined models are shown in Figure3.Figure 3
The combined models with HP Filter-SARIMA and the revised ARMAX.Figure 4
Electric load and temperature figure.


## 3.1. Modeling Steps

Step 1.
The original sequenceyt or xt, defined as the superposition of the waves with different frequencies, can be divided into the sequence G and the sequence C [25], the separation process of which must satisfy the minimum loss function principle:(20)min⁡∑t=1nYtGC-YtG2+λ∑Yt+1G-YtG-YtG-Yt-1G2.Step 2.
The sequenceG is the function of time. The error term is denoted as ε1,ε2,…,εn. The method of ordinary least squares (OLS) is used in making polynomial curve fitting to minimize the sum of squared errors. Polynomial fitting can be represented as(21)yG=a1tn+a2tn-1+⋯+an-1t+an.
The sequences’ predictions can be got and defined asg^t.Step 3.
The stationary property of the sequenceC is tested. If it is stationary, the correlation analysis can be used on the sequence; If it is stationary, the correlation analysis can be used on the sequence; if otherwise, the difference transformation is conducted on the sequence until it is stationary [26, 27]. The stationary sequence is denoted as wc.Step 4.
Through the correlation analysis, the autocorrelation and moving average items can be acquired. The ARIMA model is presented as(22)yc=φp-1LΦP-1L1-L-d1-Ls-DθqLΘQLεc.Step 5.
Consider testing of the rationality of the model by the residual sequence.Step 6.
By ARIMA model, the predictions of the sequences can be denoted asc^t.Step 7.
The final prediction result can be got according to the principle of HP filter:(23)y^t=g^t+c^t.Step 8.
Consider exploring the correlation coefficient between the stationary sequences “y^t” and “X^t” to determine the structure of the improved ARMAX model. This step is an improved version for traditional ARMAX model. The revised ARMAX model can be calculated as follows: (24)yt∗=μ+∑i=1kΘiBΦiBBlixit∗+εt.Step 9.
Consider fitting residual sequenceεt:(25)εt=ΘBΦBat,where at is a zero mean white noise sequence.Based on above steps, the combined “HP Filter-SARIMA-Revised ARMAX” model can be applied in load forecasting process.


## 3.2. Modeling Flowchart

The forecasting processes of the combined models are shown in Figure3.Figure 3
The combined models with HP Filter-SARIMA and the revised ARMAX.Figure 4
Electric load and temperature figure.


## 4. Load Forecasting with ARMA Model

### 4.1. Data Source

Figure4 shows the daily maximum power load and the maximum temperature in a city in China from May 1st to July 15th. The specific information about the city is not allowed to be shared here. In this paper, the classical time series models and the combined models are applied to forecast load. The prediction results of different models are compared in Section 5.It is shown in Figure4 that load data has clear cyclical fluctuations by observing the sequences. The sequences are obviously nonstationary.


### 4.2. Establishing Seasonal-ARIMA Model

The ADF unit root test demonstrates that the original sequence is nonstationary, while the first-order difference of the original sequence is stationary under the 5% significant level [28]. Therefore, the additive model is used to adjust the seasonal sequence. The analysis on the partial autocorrelation and the autocorrelation is shown in Figure 5.Figure 5
Autocorrelation figure.It is shown in Figure5 that the autocorrelation of daily peak power load data is in stationary and periodic series. Partial autocorrelation shows that only the first-order partial autocorrelation coefficient is significantly greater than two times standard deviation [29]. The rest partial autocorrelation coefficient rapidly declines to zero, making random fluctuations within two times standard deviation ranges. Thus it may be regarded as the first-order truncation [30].After determining the model order by parameter significance testing, the time series model obtained isARIMA(2,1,11)(1,1,0)12 [31] in Tables 1 and 2. In this model, the nonseasonal autoregressive items are AR(1) and AR(2), the nonseasonal moving average items are MA(1) and MA(11), and the seasonal autoregressive items are SAR(1).Table 1
The parameter estimating of the ARIMA model.
ARIMA model
AR(1)
AR(2)
AR(11)
SAR(1)
ARIMA(11,1,13)(1,1,1)12
Coefficient(t-statistic)
−0.742(−3.265)
−0.225(−1.699)
0.345(4.954)
0.543(4.956)
ARIMA(2,1,11)(1,1,0)12
Coefficient(t-statistic)
−0.368(−2.353)
−0.277(−2.488)
−0.660(−3.991)
MA(1)
MA(11)
MA(13)
SMA(1)
ARIMA(11,1,13)(1,1,1)12
Coefficient(t-statistic)
−1.148(−17.508)
−0.421(−6.468)
0.536(7.110)
0.401(1.638)
ARIMA(2,1,11)(1,1,0)12
Coefficient(t-statistic)
0.614(5.127)
0.341(3.251)Table 2
Comparison of AIC and SBC.
R
2
AIC
SBC
ARIMA(11,1,13)(1,1,1)12
0.677
12.991
13.206
ARIMA(2,1,11)(1,1,0)12
0.397
13.460
13.588In the process of model order selection,MA(11) is greatly influential in modeling and predicting. This phenomenon shows that the errors caused by long-term observations still affect the current monthly electricity consumption to some extent [32]. It is inferred that monthly electricity consumption may have the characteristic of long-term memory. This conjecture is confirmed by long-term memory test using the Rescaled Range Analysis (R/S) method, which shows the Hurst exponent: (26)H=log⁡R/S-log⁡Klogn,H=0.837.According to the criteria,0.5<0.837<1, it is known that the monthly electricity consumption has the long-term memory characteristic. Therefore, the current forecast is influenced by distant observations [33, 34].


### 4.3. Integrative Models Using HP Filter

HP filter is applied in the analysis of original sequence. The decomposition results are shown in Figure6, in which the blue curve is the original sequence, while the red curve is the long-term trend. It can be found that the growth rate of the power consumption is mainly constant from May 1st to July 15th. The green curve changes cyclically and irregularly. As time goes on, the fluctuation range is larger.Figure 6
Separated electric load sequence by HP filter.According to HP filter, the original electric load sequence can be separated into the sequenceG with long-term trend and the sequence C with other fluctuant properties [35].Data separation results using HP filter are shown in Figures7 and 8. It can be seen that the sequence G approximates a smooth curve, while the curve of the sequence C fluctuates up and down around zero.Figure 7
Trend sequence of electric load.Figure 8
Cycle sequence of electric load.In the perspective of statistics, the sequenceG can be transformed into a time-related sequence. The independent variable t represents time, while the sequence G shows the dependent variable of the system.The fitting polynomial of electric load data is(27)G=-0.0001779687t2+0.0290533421t+1.7012003425.Fitting figure of load is shown in Figure9.Figure 9
Fitting polynomial of electric load.The sequenceG of temperature from May 1st to July 15th can be forecasted based on fitting polynomial [36], while the sequence C is predicted based on ARMA(1,2).The same procedure can be applied to temperature data.Based on HP filter, the original temperature sequence can be separated into the sequenceG′ with long-term trend and the sequence C′ with other fluctuant properties in Figure 10.Figure 10
Separated temperature sequence by HP filter.Data separation results using HP filter are shown in Figures11 and 12. It is shown that the sequence G′ looks like a smooth curve, while the curve of the sequence C′ fluctuates up and down around zero.Figure 11
Trend sequence of temperature.Figure 12
Cycle sequence of temperature.The fitting polynomial of temperature data is(28)G=0.0000000040t4-0.0000006727t3+0.0000198680t2+0.2283366535t+22.505274.The sequenceG′ of the temperature from May 1st to July 15th can be forecasted based on the fitting polynomial, while the sequence C′ is predicted based on ARMA(1,1). Fitting figure of temperature is shown in Figure 13.Figure 13
Fitting polynomial of temperature.Finally, the input value of revised ARMAX model is obtained based on the HP filter principley^t=g^t+c^t, x^t=g^t+c^t.


### 4.4. Load Forecasting with Combined Models

#### 4.4.1.y^t and x^t Sequences Model

Firstly, they^t model is established. It is shown in the test that y^t is a stationary white noise sequence; thus the fitting model is(29)y^=εt.It is indicated in Table3 that x^ is a stationary white noise sequence.Table 3
Autocorrelation check for white noise.
To lag
Pr > ChiSq
6
0.1156
12
0.2341The best order forARMA(p,q) model is [AR(0),MA(1)], which is shown in Table 4.Table 4
Minimum information criterion.
MA(0)
MA(1)
MA(2)
MA(3)
MA(4)
MA(5)
AR(0)
−4.32533
−4.42021
−4.30792
−4.23921
−4.17551
−4.07372
AR(1)
−4.36563
−4.29661
−4.18864
−4.11551
−4.0525
−3.95355
AR(2)
−4.32658
−4.20127
−4.07613
−3.99036
−3.92849
−3.82853
AR(3)
−4.23955
−4.11437
−3.99
−3.8682
−3.87325
−3.82193
AR(4)
−4.17601
−4.05104
−3.92805
−3.8028
−3.76851
−3.75611
AR(5)
−4.07686
−3.95249
−3.82807
−3.83598
−3.71155
−3.79342
Minimum table value: BIC(0,1) = −4.42021Secondly, thex^ model is set up. The test results obtained by SAS show that x^ is a stationary white noise sequence [37]. Therefore, the fitting model is ARMA(0,1) or MA(1) model. Fitting parameter of x^ is shown in Table 5.Table 5
Fitting parameter ofx^.
Parameter
Estimate
t-value
Approx. Pr >t
Lag
MA(1,1)
0.39758
2.06
0.0495
1The final fitting model is shown below:(30)x^t=1-0.39758Bεt.


#### 4.4.2. Computing Load Data with Revised ARMAX Model

The above model is used to filter input variable sequencex^1,…,x^k and the response variable sequence y^t, which is followed by the calculation of mutual relationship number between the independent variables and the response variable after filtration.The regression analysis in Table6 shows that the final regression coefficient is 0.62871. The statistics test is conducted with residual sequence, showing that the residual sequence is stationary white noise sequence (Pr>0.05) [38, 39]. The fitted model for residual sequence is εt=at, where at is zero mean white noise sequence.Table 6
Parameter estimates of REG procedure.
Variable
Dependent variable
Parameter estimates
t-value
Pr >t
y
^
x
^
0.62871
6.75
<0.0001The number of 0-order delay mutual relationship is significantly nonzero, which means that there is no hysteretic effect between response sequence and input sequence. Thus the model should be treated in the same period [40].The statistics test is operated with residual sequence, showing that the residual sequence is stationary white noise sequence. The fitted model for residual sequence isεt=at, where at is zero mean white noise sequence [41–43]. It is known that there is significant correlation in the zero order between the two sequences. The same period model between y^ and x^ is established based on the results in Table 6:(31)Δy^t=0.62871Δx^t+εt,εt=at.The load from 16th to 31st is forecasted according to the combined models, HP Filter-SARIMA-Revised ARMAX. The second column in Table7 is the actual load data (104 kw·h). The third column is the actual temperature (°C). The fourth is the prediction value of SARIMA model. The fifth is prediction value of ARIMA model. The sixth is prediction value of grey system theory. The seventh is prediction value of MLP neural network.Table 7
Predictions of five models.
Obs.
Actual load
Actual temperature
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
16th Jul.
4585
35.0
4570
−0.3272%
4748
3.5551%
4287
−6.4995%
4587
0.0436%
4428
−3.4242%
17th Jul.
4733
39.4
4699
−0.7184%
4990
5.4300%
4210
−11.0501%
4457
−5.8314%
4558
−3.6974%
18th Jul.
4824
39.9
4799
−0.5182%
5158
6.9237%
4320
−10.4478%
4657
−3.4619%
4721
−2.1352%
19th Jul.
5012
41.9
4989
−0.4589%
5520
10.1357%
4580
−8.6193%
4851
−3.2123%
4518
−9.8563%
20th Jul.
4293
31.9
4308
0.3494%
4824
12.3690%
4127
−3.8668%
4427
3.1214%
4212
−1.8868%
21st Jul.
4470
39.0
4510
0.8949%
4750
6.2640%
4215
−5.7047%
4639
3.7808%
4498
0.6264%
22nd Jul.
4589
40.8
4588
−0.0218%
4650
1.3293%
4218
−8.0846%
4418
−3.7263%
4357
−5.0556%
23rd Jul.
4630
33.3
4599
−0.6695%
4598
−0.6911%
4324
−6.6091%
4417
−4.6004%
4500
−2.8078%
24th Jul.
4497
35.9
4485
−0.2668%
4325
−3.8248%
4129
−8.1832%
4213
−6.3153%
4367
−2.8908%
25th Jul.
4116
29.3
4120
0.0972%
3998
−2.8669%
3950
−4.0330%
3854
−6.3654%
4007
−2.6482%
26th Jul.
3835
26.4
3827
−0.2086%
3560
−7.1708%
3648
−4.8761%
3712
−3.2073%
3714
−3.1551%
27th Jul.
3726
27.1
3697
−0.7783%
3589
−3.6769%
3457
−7.2195%
3320
−10.8964%
3421
−8.1857%
28th Jul.
3797
33.3
3785
−0.3160%
3498
−7.8746%
3321
−12.5362%
3341
−12.0095%
3516
−7.4006%
29th Jul.
3845
35.5
3821
−0.6242%
3642
−5.2796%
3214
−16.4109%
3517
−8.5306%
3578
−6.9441%
30th Jul.
4144
35.8
4129
−0.3620%
3859
−6.8774%
3357
−18.9913%
3524
−14.9614%
3874
−6.5154%
31st Jul.
4061
36.0
4054
−0.1724%
3745
−7.7813%
3458
−14.8486%
3642
−10.3177%
3821
−5.9099%In our work, the prediction data points and the errors for five models from July 16th to July 31st (training set) have been conducted to assess the models’ fit performance. Table8 lists the comparison of simulation performance among five models for training set, and Figures 14 and 15 show the corresponding error radar chart and histograms for direct observation.Table 8
The comparison of error among five models for training set.
Model
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
RMSE (root mean square error)
21.47818
282.50409
433.36301
301.64787
224.33025
MAE (mean absolute error)
18.4375
248.1875
396.375
261.9375
195.1875
MAPE (mean absolute percentage error)
0.424
5.7531
9.2488
6.2738
4.5712Figure 14
The error radar chart for combined model and other comparative models.Figure 15
The error histogram for combined model and other comparative models.The error criterion indicators are expressed as follows:(32)ERROR=y^t-ytyt×100,t=1,2,…,n,MAPE=1n∑t=1ny^t-ytyt×100,MAE=∑t=1Ty^t-ytT,RMSE=1n∑t=1ny^t-yt2.Forecast graphic is displayed in Figure16.Figure 16
Forecasting graphics of different methods.It is shown in residual stationarity and white noise test that the residual is stationary white noise sequence,εt=at. There is a second-order delay correlation between y^ and x^. The final fitting model is(33)y^t=4.51539+0.97436x^t-2+εt,εt=at.


## 4.1. Data Source

Figure4 shows the daily maximum power load and the maximum temperature in a city in China from May 1st to July 15th. The specific information about the city is not allowed to be shared here. In this paper, the classical time series models and the combined models are applied to forecast load. The prediction results of different models are compared in Section 5.It is shown in Figure4 that load data has clear cyclical fluctuations by observing the sequences. The sequences are obviously nonstationary.


## 4.2. Establishing Seasonal-ARIMA Model

The ADF unit root test demonstrates that the original sequence is nonstationary, while the first-order difference of the original sequence is stationary under the 5% significant level [28]. Therefore, the additive model is used to adjust the seasonal sequence. The analysis on the partial autocorrelation and the autocorrelation is shown in Figure 5.Figure 5
Autocorrelation figure.It is shown in Figure5 that the autocorrelation of daily peak power load data is in stationary and periodic series. Partial autocorrelation shows that only the first-order partial autocorrelation coefficient is significantly greater than two times standard deviation [29]. The rest partial autocorrelation coefficient rapidly declines to zero, making random fluctuations within two times standard deviation ranges. Thus it may be regarded as the first-order truncation [30].After determining the model order by parameter significance testing, the time series model obtained isARIMA(2,1,11)(1,1,0)12 [31] in Tables 1 and 2. In this model, the nonseasonal autoregressive items are AR(1) and AR(2), the nonseasonal moving average items are MA(1) and MA(11), and the seasonal autoregressive items are SAR(1).Table 1
The parameter estimating of the ARIMA model.
ARIMA model
AR(1)
AR(2)
AR(11)
SAR(1)
ARIMA(11,1,13)(1,1,1)12
Coefficient(t-statistic)
−0.742(−3.265)
−0.225(−1.699)
0.345(4.954)
0.543(4.956)
ARIMA(2,1,11)(1,1,0)12
Coefficient(t-statistic)
−0.368(−2.353)
−0.277(−2.488)
−0.660(−3.991)
MA(1)
MA(11)
MA(13)
SMA(1)
ARIMA(11,1,13)(1,1,1)12
Coefficient(t-statistic)
−1.148(−17.508)
−0.421(−6.468)
0.536(7.110)
0.401(1.638)
ARIMA(2,1,11)(1,1,0)12
Coefficient(t-statistic)
0.614(5.127)
0.341(3.251)Table 2
Comparison of AIC and SBC.
R
2
AIC
SBC
ARIMA(11,1,13)(1,1,1)12
0.677
12.991
13.206
ARIMA(2,1,11)(1,1,0)12
0.397
13.460
13.588In the process of model order selection,MA(11) is greatly influential in modeling and predicting. This phenomenon shows that the errors caused by long-term observations still affect the current monthly electricity consumption to some extent [32]. It is inferred that monthly electricity consumption may have the characteristic of long-term memory. This conjecture is confirmed by long-term memory test using the Rescaled Range Analysis (R/S) method, which shows the Hurst exponent: (26)H=log⁡R/S-log⁡Klogn,H=0.837.According to the criteria,0.5<0.837<1, it is known that the monthly electricity consumption has the long-term memory characteristic. Therefore, the current forecast is influenced by distant observations [33, 34].


## 4.3. Integrative Models Using HP Filter

HP filter is applied in the analysis of original sequence. The decomposition results are shown in Figure6, in which the blue curve is the original sequence, while the red curve is the long-term trend. It can be found that the growth rate of the power consumption is mainly constant from May 1st to July 15th. The green curve changes cyclically and irregularly. As time goes on, the fluctuation range is larger.Figure 6
Separated electric load sequence by HP filter.According to HP filter, the original electric load sequence can be separated into the sequenceG with long-term trend and the sequence C with other fluctuant properties [35].Data separation results using HP filter are shown in Figures7 and 8. It can be seen that the sequence G approximates a smooth curve, while the curve of the sequence C fluctuates up and down around zero.Figure 7
Trend sequence of electric load.Figure 8
Cycle sequence of electric load.In the perspective of statistics, the sequenceG can be transformed into a time-related sequence. The independent variable t represents time, while the sequence G shows the dependent variable of the system.The fitting polynomial of electric load data is(27)G=-0.0001779687t2+0.0290533421t+1.7012003425.Fitting figure of load is shown in Figure9.Figure 9
Fitting polynomial of electric load.The sequenceG of temperature from May 1st to July 15th can be forecasted based on fitting polynomial [36], while the sequence C is predicted based on ARMA(1,2).The same procedure can be applied to temperature data.Based on HP filter, the original temperature sequence can be separated into the sequenceG′ with long-term trend and the sequence C′ with other fluctuant properties in Figure 10.Figure 10
Separated temperature sequence by HP filter.Data separation results using HP filter are shown in Figures11 and 12. It is shown that the sequence G′ looks like a smooth curve, while the curve of the sequence C′ fluctuates up and down around zero.Figure 11
Trend sequence of temperature.Figure 12
Cycle sequence of temperature.The fitting polynomial of temperature data is(28)G=0.0000000040t4-0.0000006727t3+0.0000198680t2+0.2283366535t+22.505274.The sequenceG′ of the temperature from May 1st to July 15th can be forecasted based on the fitting polynomial, while the sequence C′ is predicted based on ARMA(1,1). Fitting figure of temperature is shown in Figure 13.Figure 13
Fitting polynomial of temperature.Finally, the input value of revised ARMAX model is obtained based on the HP filter principley^t=g^t+c^t, x^t=g^t+c^t.


## 4.4. Load Forecasting with Combined Models

### 4.4.1.y^t and x^t Sequences Model

Firstly, they^t model is established. It is shown in the test that y^t is a stationary white noise sequence; thus the fitting model is(29)y^=εt.It is indicated in Table3 that x^ is a stationary white noise sequence.Table 3
Autocorrelation check for white noise.
To lag
Pr > ChiSq
6
0.1156
12
0.2341The best order forARMA(p,q) model is [AR(0),MA(1)], which is shown in Table 4.Table 4
Minimum information criterion.
MA(0)
MA(1)
MA(2)
MA(3)
MA(4)
MA(5)
AR(0)
−4.32533
−4.42021
−4.30792
−4.23921
−4.17551
−4.07372
AR(1)
−4.36563
−4.29661
−4.18864
−4.11551
−4.0525
−3.95355
AR(2)
−4.32658
−4.20127
−4.07613
−3.99036
−3.92849
−3.82853
AR(3)
−4.23955
−4.11437
−3.99
−3.8682
−3.87325
−3.82193
AR(4)
−4.17601
−4.05104
−3.92805
−3.8028
−3.76851
−3.75611
AR(5)
−4.07686
−3.95249
−3.82807
−3.83598
−3.71155
−3.79342
Minimum table value: BIC(0,1) = −4.42021Secondly, thex^ model is set up. The test results obtained by SAS show that x^ is a stationary white noise sequence [37]. Therefore, the fitting model is ARMA(0,1) or MA(1) model. Fitting parameter of x^ is shown in Table 5.Table 5
Fitting parameter ofx^.
Parameter
Estimate
t-value
Approx. Pr >t
Lag
MA(1,1)
0.39758
2.06
0.0495
1The final fitting model is shown below:(30)x^t=1-0.39758Bεt.


### 4.4.2. Computing Load Data with Revised ARMAX Model

The above model is used to filter input variable sequencex^1,…,x^k and the response variable sequence y^t, which is followed by the calculation of mutual relationship number between the independent variables and the response variable after filtration.The regression analysis in Table6 shows that the final regression coefficient is 0.62871. The statistics test is conducted with residual sequence, showing that the residual sequence is stationary white noise sequence (Pr>0.05) [38, 39]. The fitted model for residual sequence is εt=at, where at is zero mean white noise sequence.Table 6
Parameter estimates of REG procedure.
Variable
Dependent variable
Parameter estimates
t-value
Pr >t
y
^
x
^
0.62871
6.75
<0.0001The number of 0-order delay mutual relationship is significantly nonzero, which means that there is no hysteretic effect between response sequence and input sequence. Thus the model should be treated in the same period [40].The statistics test is operated with residual sequence, showing that the residual sequence is stationary white noise sequence. The fitted model for residual sequence isεt=at, where at is zero mean white noise sequence [41–43]. It is known that there is significant correlation in the zero order between the two sequences. The same period model between y^ and x^ is established based on the results in Table 6:(31)Δy^t=0.62871Δx^t+εt,εt=at.The load from 16th to 31st is forecasted according to the combined models, HP Filter-SARIMA-Revised ARMAX. The second column in Table7 is the actual load data (104 kw·h). The third column is the actual temperature (°C). The fourth is the prediction value of SARIMA model. The fifth is prediction value of ARIMA model. The sixth is prediction value of grey system theory. The seventh is prediction value of MLP neural network.Table 7
Predictions of five models.
Obs.
Actual load
Actual temperature
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
16th Jul.
4585
35.0
4570
−0.3272%
4748
3.5551%
4287
−6.4995%
4587
0.0436%
4428
−3.4242%
17th Jul.
4733
39.4
4699
−0.7184%
4990
5.4300%
4210
−11.0501%
4457
−5.8314%
4558
−3.6974%
18th Jul.
4824
39.9
4799
−0.5182%
5158
6.9237%
4320
−10.4478%
4657
−3.4619%
4721
−2.1352%
19th Jul.
5012
41.9
4989
−0.4589%
5520
10.1357%
4580
−8.6193%
4851
−3.2123%
4518
−9.8563%
20th Jul.
4293
31.9
4308
0.3494%
4824
12.3690%
4127
−3.8668%
4427
3.1214%
4212
−1.8868%
21st Jul.
4470
39.0
4510
0.8949%
4750
6.2640%
4215
−5.7047%
4639
3.7808%
4498
0.6264%
22nd Jul.
4589
40.8
4588
−0.0218%
4650
1.3293%
4218
−8.0846%
4418
−3.7263%
4357
−5.0556%
23rd Jul.
4630
33.3
4599
−0.6695%
4598
−0.6911%
4324
−6.6091%
4417
−4.6004%
4500
−2.8078%
24th Jul.
4497
35.9
4485
−0.2668%
4325
−3.8248%
4129
−8.1832%
4213
−6.3153%
4367
−2.8908%
25th Jul.
4116
29.3
4120
0.0972%
3998
−2.8669%
3950
−4.0330%
3854
−6.3654%
4007
−2.6482%
26th Jul.
3835
26.4
3827
−0.2086%
3560
−7.1708%
3648
−4.8761%
3712
−3.2073%
3714
−3.1551%
27th Jul.
3726
27.1
3697
−0.7783%
3589
−3.6769%
3457
−7.2195%
3320
−10.8964%
3421
−8.1857%
28th Jul.
3797
33.3
3785
−0.3160%
3498
−7.8746%
3321
−12.5362%
3341
−12.0095%
3516
−7.4006%
29th Jul.
3845
35.5
3821
−0.6242%
3642
−5.2796%
3214
−16.4109%
3517
−8.5306%
3578
−6.9441%
30th Jul.
4144
35.8
4129
−0.3620%
3859
−6.8774%
3357
−18.9913%
3524
−14.9614%
3874
−6.5154%
31st Jul.
4061
36.0
4054
−0.1724%
3745
−7.7813%
3458
−14.8486%
3642
−10.3177%
3821
−5.9099%In our work, the prediction data points and the errors for five models from July 16th to July 31st (training set) have been conducted to assess the models’ fit performance. Table8 lists the comparison of simulation performance among five models for training set, and Figures 14 and 15 show the corresponding error radar chart and histograms for direct observation.Table 8
The comparison of error among five models for training set.
Model
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
RMSE (root mean square error)
21.47818
282.50409
433.36301
301.64787
224.33025
MAE (mean absolute error)
18.4375
248.1875
396.375
261.9375
195.1875
MAPE (mean absolute percentage error)
0.424
5.7531
9.2488
6.2738
4.5712Figure 14
The error radar chart for combined model and other comparative models.Figure 15
The error histogram for combined model and other comparative models.The error criterion indicators are expressed as follows:(32)ERROR=y^t-ytyt×100,t=1,2,…,n,MAPE=1n∑t=1ny^t-ytyt×100,MAE=∑t=1Ty^t-ytT,RMSE=1n∑t=1ny^t-yt2.Forecast graphic is displayed in Figure16.Figure 16
Forecasting graphics of different methods.It is shown in residual stationarity and white noise test that the residual is stationary white noise sequence,εt=at. There is a second-order delay correlation between y^ and x^. The final fitting model is(33)y^t=4.51539+0.97436x^t-2+εt,εt=at.


## 4.4.1.y^t and x^t Sequences Model

Firstly, they^t model is established. It is shown in the test that y^t is a stationary white noise sequence; thus the fitting model is(29)y^=εt.It is indicated in Table3 that x^ is a stationary white noise sequence.Table 3
Autocorrelation check for white noise.
To lag
Pr > ChiSq
6
0.1156
12
0.2341The best order forARMA(p,q) model is [AR(0),MA(1)], which is shown in Table 4.Table 4
Minimum information criterion.
MA(0)
MA(1)
MA(2)
MA(3)
MA(4)
MA(5)
AR(0)
−4.32533
−4.42021
−4.30792
−4.23921
−4.17551
−4.07372
AR(1)
−4.36563
−4.29661
−4.18864
−4.11551
−4.0525
−3.95355
AR(2)
−4.32658
−4.20127
−4.07613
−3.99036
−3.92849
−3.82853
AR(3)
−4.23955
−4.11437
−3.99
−3.8682
−3.87325
−3.82193
AR(4)
−4.17601
−4.05104
−3.92805
−3.8028
−3.76851
−3.75611
AR(5)
−4.07686
−3.95249
−3.82807
−3.83598
−3.71155
−3.79342
Minimum table value: BIC(0,1) = −4.42021Secondly, thex^ model is set up. The test results obtained by SAS show that x^ is a stationary white noise sequence [37]. Therefore, the fitting model is ARMA(0,1) or MA(1) model. Fitting parameter of x^ is shown in Table 5.Table 5
Fitting parameter ofx^.
Parameter
Estimate
t-value
Approx. Pr >t
Lag
MA(1,1)
0.39758
2.06
0.0495
1The final fitting model is shown below:(30)x^t=1-0.39758Bεt.


## 4.4.2. Computing Load Data with Revised ARMAX Model

The above model is used to filter input variable sequencex^1,…,x^k and the response variable sequence y^t, which is followed by the calculation of mutual relationship number between the independent variables and the response variable after filtration.The regression analysis in Table6 shows that the final regression coefficient is 0.62871. The statistics test is conducted with residual sequence, showing that the residual sequence is stationary white noise sequence (Pr>0.05) [38, 39]. The fitted model for residual sequence is εt=at, where at is zero mean white noise sequence.Table 6
Parameter estimates of REG procedure.
Variable
Dependent variable
Parameter estimates
t-value
Pr >t
y
^
x
^
0.62871
6.75
<0.0001The number of 0-order delay mutual relationship is significantly nonzero, which means that there is no hysteretic effect between response sequence and input sequence. Thus the model should be treated in the same period [40].The statistics test is operated with residual sequence, showing that the residual sequence is stationary white noise sequence. The fitted model for residual sequence isεt=at, where at is zero mean white noise sequence [41–43]. It is known that there is significant correlation in the zero order between the two sequences. The same period model between y^ and x^ is established based on the results in Table 6:(31)Δy^t=0.62871Δx^t+εt,εt=at.The load from 16th to 31st is forecasted according to the combined models, HP Filter-SARIMA-Revised ARMAX. The second column in Table7 is the actual load data (104 kw·h). The third column is the actual temperature (°C). The fourth is the prediction value of SARIMA model. The fifth is prediction value of ARIMA model. The sixth is prediction value of grey system theory. The seventh is prediction value of MLP neural network.Table 7
Predictions of five models.
Obs.
Actual load
Actual temperature
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
Forecast
Error (%)
16th Jul.
4585
35.0
4570
−0.3272%
4748
3.5551%
4287
−6.4995%
4587
0.0436%
4428
−3.4242%
17th Jul.
4733
39.4
4699
−0.7184%
4990
5.4300%
4210
−11.0501%
4457
−5.8314%
4558
−3.6974%
18th Jul.
4824
39.9
4799
−0.5182%
5158
6.9237%
4320
−10.4478%
4657
−3.4619%
4721
−2.1352%
19th Jul.
5012
41.9
4989
−0.4589%
5520
10.1357%
4580
−8.6193%
4851
−3.2123%
4518
−9.8563%
20th Jul.
4293
31.9
4308
0.3494%
4824
12.3690%
4127
−3.8668%
4427
3.1214%
4212
−1.8868%
21st Jul.
4470
39.0
4510
0.8949%
4750
6.2640%
4215
−5.7047%
4639
3.7808%
4498
0.6264%
22nd Jul.
4589
40.8
4588
−0.0218%
4650
1.3293%
4218
−8.0846%
4418
−3.7263%
4357
−5.0556%
23rd Jul.
4630
33.3
4599
−0.6695%
4598
−0.6911%
4324
−6.6091%
4417
−4.6004%
4500
−2.8078%
24th Jul.
4497
35.9
4485
−0.2668%
4325
−3.8248%
4129
−8.1832%
4213
−6.3153%
4367
−2.8908%
25th Jul.
4116
29.3
4120
0.0972%
3998
−2.8669%
3950
−4.0330%
3854
−6.3654%
4007
−2.6482%
26th Jul.
3835
26.4
3827
−0.2086%
3560
−7.1708%
3648
−4.8761%
3712
−3.2073%
3714
−3.1551%
27th Jul.
3726
27.1
3697
−0.7783%
3589
−3.6769%
3457
−7.2195%
3320
−10.8964%
3421
−8.1857%
28th Jul.
3797
33.3
3785
−0.3160%
3498
−7.8746%
3321
−12.5362%
3341
−12.0095%
3516
−7.4006%
29th Jul.
3845
35.5
3821
−0.6242%
3642
−5.2796%
3214
−16.4109%
3517
−8.5306%
3578
−6.9441%
30th Jul.
4144
35.8
4129
−0.3620%
3859
−6.8774%
3357
−18.9913%
3524
−14.9614%
3874
−6.5154%
31st Jul.
4061
36.0
4054
−0.1724%
3745
−7.7813%
3458
−14.8486%
3642
−10.3177%
3821
−5.9099%In our work, the prediction data points and the errors for five models from July 16th to July 31st (training set) have been conducted to assess the models’ fit performance. Table8 lists the comparison of simulation performance among five models for training set, and Figures 14 and 15 show the corresponding error radar chart and histograms for direct observation.Table 8
The comparison of error among five models for training set.
Model
HP Filter-SARIMA-Revised ARMAX
SARIMA
ARIMA
Grey system
MLP neural network
RMSE (root mean square error)
21.47818
282.50409
433.36301
301.64787
224.33025
MAE (mean absolute error)
18.4375
248.1875
396.375
261.9375
195.1875
MAPE (mean absolute percentage error)
0.424
5.7531
9.2488
6.2738
4.5712Figure 14
The error radar chart for combined model and other comparative models.Figure 15
The error histogram for combined model and other comparative models.The error criterion indicators are expressed as follows:(32)ERROR=y^t-ytyt×100,t=1,2,…,n,MAPE=1n∑t=1ny^t-ytyt×100,MAE=∑t=1Ty^t-ytT,RMSE=1n∑t=1ny^t-yt2.Forecast graphic is displayed in Figure16.Figure 16
Forecasting graphics of different methods.It is shown in residual stationarity and white noise test that the residual is stationary white noise sequence,εt=at. There is a second-order delay correlation between y^ and x^. The final fitting model is(33)y^t=4.51539+0.97436x^t-2+εt,εt=at.


## 5. Conclusion

Based on the above analysis, it is shown that the combined HP Filter-SARIMA-Revised ARMAX models can effectively forecast electric load with external variables. The structure of the model is determined by the statistical regression analysis and the least squares method. The process strictly follows the rules oft-test, AIC, and SBC. The combined models are more accurate than the single forecasting method for short-term electric load forecasting.A total of four traditional forecasting models are applied to forecast electric loads in this paper. It has been proved in empirical analysis that the combined models have small relative error compared with the traditional methods. The prediction accuracy of the combined models is greatly improved. The significance of parameters has been greatly enhanced. With the existing literatures and the analysis in this paper, researchers may find out that combined models have better prediction performance in dealing with special data compared with a single model.



---

*Source: 386925-2015-12-16.xml*