# Hybrid Load Forecasting Using Gaussian Process Regression and Novel Residual Prediction
## 4 Energy Load Forecasting Method
Marginal distribution p(\(\boldsymbol{f}\)) is defined as having mean zero and a Gram K matrix:  
\[\boldsymbol{p}(\boldsymbol{f})=\mathcal{N}\left(\boldsymbol{f}|\boldsymbol{0},K\right),\ \ \ K_{i,j}=k\left(x_{j},x_{j}\right) \tag{13}\]  
With the above knowledge, the marginal distribution of \(y\) is computed:  
\[\boldsymbol{p}(\boldsymbol{y})=\int\ p(\boldsymbol{y}\middle|\boldsymbol{f})p( \boldsymbol{f})df=\mathcal{N}\left(\boldsymbol{f}|\boldsymbol{0},K+\sigma_{n} ^{2}I\right) \tag{14}\]  
To ground the predictions in theory, for the predicted output \(y_{*}\) with the input vector \(x_{*}\), all the parameters' values are averaged and weighted according to their posterior probability:  
\[\left[ \begin{array}{c}\boldsymbol{y}\\ y_{*}\end{array} \right]\sim\mathcal{N}\left(\boldsymbol{0},\left[ \begin{array}[ ]{cc}K_{y}&\boldsymbol{k_{*}}\\ \boldsymbol{k_{*}^{T}}&k_{**+}\sigma_{n}^{2}\end{array} \right]\right) \tag{15}\]  
If there are \(\omega\) training data points, then \(\boldsymbol{k_{*}}\) is the \(\omega\) order matrix of the covariances evaluated with all the data points for training and testing. Having considered the constraints of the conditioning Gaussians [33], the prediction can be determined:  
\[p\left(y_{*}|\boldsymbol{y}\right)\sim\mathcal{N}\left(k_{*}^{T}\mathrm{K}_{y }^{-1}\boldsymbol{y}\middle|\boldsymbol{k_{*}}-k_{*}^{T}\mathrm{K}_{y}^{-1} \boldsymbol{k_{*}}+\sigma_{n}^{2}\right) \tag{16}\]  
where \(k_{**}=k(x_{*},x_{*})\) and \(K_{y}\) is the covariance matrix. The kernel, or covariance, function is that which defines the distance and similarity between the training data points. After testing various functions, the best performing kernel function was the RQ:  
\[k(x,x^{\prime})=\sigma_{n}^{2}\left(1+\frac{(x-x^{\prime})^{2}}{2\alpha l^{2} }\right)^{-\alpha} \tag{17}\]  
This can also be represented as an infinite sum of squared exponential covariance functions with different length scale characteristics (l) and scale mixture (\(\alpha\)). RQ kernel was determined as the best choice to minimize the distance between two neighbors through experimental testing. It was decided that for this specific data string from all tested functions the best match was RQ. After integrating the kernel function, the prediction mean and covariance may be calculated:  
\[\left\{ \begin{array}{c}m(\textbf{x}_{*})=\textbf{c}^{T}\textbf{x}_{*}+\textbf{ k}_{*}^{T}K_{y}^{-1}(\textbf{y}-\textbf{m})\\ \sigma^{2}(\textbf{x}_{*})=k_{**}-\textbf{k}_{*}^{T}K_{y}^{-1}\textbf{k}_{*}+ \sigma_{n}^{2}\end{array} \right. \tag{18}\]  
This paper's forecasting method is based on a hybrid of the autoregressive model and a double GPR method, which takes a predicted error into consideration:  
\[y=h(x)+\textbf{c}^{T}\textbf{x}+\epsilon(x) \tag{19}\]  
where, \(h(x)\sim\mathcal{GP}(0,\)k(x,x^{\prime}))\), \(\epsilon\) represents the method coefficients matrix and \(\epsilon(x)\) represents the predicted residuals. This formulation indicates that the output is the summation of a zero-mean Gaussian process, a set of fixed basis functions, and the zero-mean process of the residuals.