# Hybrid Load Forecasting Using Gaussian Process Regression and Novel Residual Prediction
## 4 Energy Load Forecasting Method  
For the proposed energy load prediction, several forecasting models were tested. The best match for the above-presented historical data was an AR model combined with GPR and an error-adjusting function specific to the datasets. First, a simple AR model is presented, beginning with the mathematical representation [32]:  
\[y_{m}=\ ct+\sum\
olimits_{i=1}^{n}c_{i}y_{m-1}\ +\ e_{m} \tag{6}\]  
where \(c_{i}\) is the method coefficient for the current iteration and represents how the previous value influences the current value, \(e_{m}\) represents the error for the \(m\)th value, and \(ct\) is the first coefficient of the function with a constant value. To compute the order of the AR model, a partial autocorrelation function is used with Yule-Walker equations:  
\[p_{i}=\sum\limits_{j=1}^{k}\Phi_{pj}*\sum\limits_{t=1}^{M-i}(y_{m}-\overline{y })(y_{m+1}-\overline{y})\frac{1}{ct*M} \tag{7}\]  
where \(k\) is the value at which the correlation function \(\Phi pk\) cuts off after exceeding the model order, \(M\) represents the number of observations, and \(ct\) is the mean value iteration of the finite time series. The next step for the method after determining the order of the model is to estimate the model coefficients. To accomplish this, the first equation is expressed in vector form to express the interconnection between the input and the output:  
\[\underbrace{\left[ \begin{array}{c}y_{m}\\ y_{m-1}\\ \vdots\\ y_{m-\omega+1}\end{array} \right]}_{y}=\underbrace{\left[ \begin{array}{cccc} 1&y_{m-1}&y_{m-2}&\ldots&y_{m-n}\\ 1&y_{m-2}&y_{m-3}&\ldots&y_{m-n+1}\\ \vdots&\vdots&\vdots&\vdots&\vdots\\ 1&y_{m-\omega}&y_{m-\omega+1}&\ldots&y_{m-n+\omega}\\ \end{array} \right]}_{x}\left[ \begin{array}{c}c_{t}\\ c_{1}\\ c_{n}\end{array} \right]+\left[ \begin{array}{c}e_{1}\\ e_{2}\\ \vdots\\ e_{n+1}\end{array} \right] \tag{8}\]  
Figure 6: Comparison between reconstructed empirical wavelet transform (EWT) signal and the original signal. where \(\omega\) represents the number of input/output pairs that are constructed using the original time series. Having written the equation in vector form, we can determine the coefficients of the model using the least squares method:  
\[C=\left(X^{T}X\right)^{-1}X^{T}y \tag{9}\]  
Next, the proposed GPR model is presented [33]. It can be written as \(f(\mathbf{x})\sim\mathcal{GP}(\mathrm{m}(\mathbf{x}),\mathrm{k}(\mathbf{x}, \mathbf{x}^{\prime}))\), an expression of its mean and covariance functions. The expression of the output regarding the \(n\)-order input vector is computed as follows:  
\[y=f(x)+e \tag{10}\]  
where \(e\) represents the noise, and may be approximated with a Gaussian distribution of zero-mean value, a variance of \(\sigma_{n}^{2}\), \(x\) is the input vector, and \(y\) is the scalar output:  
\[x_{i}=\left[y_{m-i},y_{m-i+1},\to y_{m-i+n-1}\right]^{T} \tag{11}\]  
The model, together with the noise assumption, determines the likelihood that represents the probability density based on the parameters. The objective is to construct a training set \(D=\left[(x_{i},y_{i})\middle|i=1,2\ldots,\omega\right]\) using the input vector and scalar output described above:  
\[p(y\middle|\boldsymbol{f})=\mathcal{N}\left(y\middle|\boldsymbol{f},\sigma_{n }^{2}I\right) \tag{12}\]  
The likelihood presented in the above expression is given by the Gaussian distribution with \(y\middle|\boldsymbol{f}\) mean and \(\sigma_{n}^{2}I\) variance, \(\boldsymbol{y}=\left[y_{1},y_{2}\ldots y_{\omega}\right]^{T}\),\(f=\left[\mathrm{f}(x_{1}),\mathrm{f}(x_{2})\ldots\mathrm{f}(x_{\omega}) \right]^{T}\), and \(I\) being the unit matrix of \(\omega\) order.