# Prediction Error and Forecasting Interval Analysis of Decision Trees with an Application in Renewable Energy Supply Forecasting

Authors: Xin Zhao; Xiaokai Nie
Journal: Complexity
 (2020)

---

## Abstract

Renewable energy has become popular compared with traditional energy like coal. The relative demand for renewable energy compared to traditional energy is an important index to determine the energy supply structure. Forecasting the relative demand index has become quite essential. Data mining methods like decision trees are quite effective in such time series forecasting, but theory behind them is rarely discussed in research. In this paper, some theories are explored about decision trees including the behavior of bias, variance, and squared prediction error using trees and the prediction interval analysis. After that, real UK grid data are used in interval forecasting application. In the renewable energy ratio forecasting application, the ratio of renewable energy supply over that of traditional energy can be dynamically forecasted with an interval coverage accuracy higher than 80% and a small width around 22, which is similar to its standard deviation.


---

## Body

## 1. Introduction

Renewable energy such as solar and wind has been playing an integral role in sustaining power supply and relieving the environment pollution and global warming crisis. With the increasing penetration of renewable energy, determining the amounts of renewable energy generation is critical to maintain the energy balance and the stability and reliability of power networks. Forecasting the mixing shares of the energy generation offers the guidance of setting up the power generation for each energy source and ensures the load demand of power networks to be satisfied [1, 2]. Data-based prediction methods, in particular machine learning methods, provide a promising solution to infer the required ratios of energy generation, among which decision tree is a well-recognized approach due to its satisfactory accuracy and interpretation [3–6].Although decision tree provides an effective method in forecasting, the theory explaining when and how it performs well is rarely discussed. The required ratios of renewable energy generation can be seen as a linear time series. In this context, we explore how the tree model performs in terms of the bias, variance, and prediction error. In addition, point prediction is not sufficient in time series prediction, so we also provide prediction interval choices like Gaussian and quantile intervals in theories with the application in renewable energy ratio forecasting.Decision tree [7] is a nonparametric supervised learning method used for discovery and prediction-oriented classification and regression. The goal is to create a model that predicts the value of a target variable by learning simple decision rules inferred from the data features. Decision tree, compared to other data mining methods, has its own advantages. (1) For casual relationship, it can deal with nonlinear models. In most cases, economics pay more attention to linear models, while if it is a nonlinear model, it will be transferred to be a linear model. In problems like consumer behavior analysis, the number of variables exceeds the normal extent to be tens or even hundreds, which will definitely lead to high correlation among variables. In that case, coefficients may have the wrong meaning in reality. Decision tree, however, provides variable importance ranking criteria, which helps a lot. (2) In terms of comprehensibility, it tends to be relatively better than “black-box” models like neural network, which means it can interpret data structure more clearly and help readers understand the information involved. These undoubtedly bring convenience to decision making in medical treatment [8–10], e-commerce, [11–13] and so on.We now explore the performance of trees when fitted to data generated from a linear model. The corresponding bias, variance, and prediction error between the fitted simplified tree and the true simple linear model will be calculated. Then, how those errors vary will be explored when the linear data distribution changes. The motivation is to explore how the trees perform under different distributions. Afterwards, prediction interval is proposed using Gaussian and quantile intervals, which explains why quantile interval is chosen in the study by Zhao et al. [14]. The simple linear model in use is(1)Y=α+βX+ɛ,where fX=α+βX is the true model. It is supposed that, throughout this paper, X∼Ua,b independently and ɛ∼N0,σ2. Uniform distribution guarantees that if the tree has k terminal nodes, the sample size in each node will be equal which is convenient in theory and simulation analysis. Decision tree analysis under the uniform distribution assumption includes the work by Hancock [15], Jackson and Servedio [16], and White and Liu [17]. Other distributions can also be considered but the analysis will be much more complex as the sample size of each terminal node depends on many parameters.The expected squared prediction error (SPE) is one of the important metrics to measure how well the trained model is applied to further unseen data. As shown in Hastie et al. [18], SPE of a regression fit f^X at an input point X=x0 is(2)SPEx0=EY−f^x02|X=x0=σ2+Ef^x0−fx02+Ef^x0−Ef^x02=σ2+Bias2f^x0+Varf^x0=irreducible error+bias2+variance.In (2), the first term is the variance of the target around its true mean fx0 and cannot be avoided no matter how well the fx0 is estimated, unless σ2=0. The second term is the squared bias, the amount by which the average of the estimate differs from the true mean; the last term is the variance, the expected squared deviation of f^x0 around its mean. Typically the more complex the model f^ is, the lower the (squared) bias but the higher the variance [18] will be.In Section2, the performance of regression trees is analyzed when fitted to data which simply follow a uniform distribution, with additive Gaussian noise. When we predict this time series using simplified trees, the prediction error is calculated and decomposed into variance and other errors. When Gaussian or uniform effect is strong, those errors have different kinds of behavior. Other exploration is conducted in Section 3 including the best tree depth with minimum prediction error and the performance of Gaussian and quantile prediction intervals under different conditions. A real interval forecasting application is conducted in Section 4. Conclusions are drawn in Section 5. All calculations were done using R [19]; ‘waveslim’ [20] was used for wavelet decomposition and ‘ctree’ [21] for CTree.


## 2. Bias-Variance Exploration

### 2.1. Decomposition Background

For theith observation Yi, the (unconditional) expectation is(3)EYi=Eα+βXi+ɛi=α+βEXi+Eɛi=α+β⋅a+b2,and the variance is(4)VarYi=Varα+βXi+ɛi=β2VarXi+Varɛi=β2⋅b−a212+σ2.They both have no relationship toi. In that case, EY=EYi and VarY=VarYi. Accordingly, for N observations, the expectation and variance for the average Y¯ are shown in(5)EY¯=1NEY1+EY2+⋯+EYN,EY=α+β⋅a+b2,(6)VarY¯=1N2VarY1+VarY2+⋯+VarYN=1NVarY=1Nβ2⋅b−a212+σ2.


### 2.2. Decomposition in the Context of Decision Trees

In the context of decision trees, the fitted model isf^X in a simplified form is(7)f^X=Y¯i,i=1,2,…,k,where k is the number of terminal nodes in the tree f^X and Y¯i is the mean of y in terminal node i. In a tree with only the root node, k=1, and the fitted model is f^X=Y¯. Then, for point x0,(8)Bias2f^x0=Ef^x0−fx02=α+β⋅a+b2−α−β⋅x02=β2a+b2−x02,and the variance is(9)Varf^x0=Ef^x0−Ef^x02=EY¯−α+β⋅a+b22=EY¯2+α+β⋅a+b22−2EY¯α+β⋅a+b2=E2Y¯+VarY¯−α+β⋅a+b22=σ2N+β2b−a212N.Thus, theSPE at point x0 is(10)SPEf^x0=σ2+bias2f^x0+Varf^x0=σ2+β2a+b2−x02+σ2N+β2b−a212N=1+1Nσ2+β2a+b2−x02+β2b−a212N.Then, the mean squared prediction error (MSPE) is(11)MSPE=∫abSPEf^xPX=xdx=∫abSPEf^x1b−adx=1+1Nσ2+β2b−a212N+β2∫aba+b2−x21b−adx=1+1Nσ2+β2b−a212N+β2112b−a2=1+1Nσ2+β2b−a212.Comprising variance is(12)EVar=1Nσ2+β2b−a212,and squared bias is(13)EBias2=β2b−a212.Now the number of terminal nodes in the decision tree is extended fromk=1 to a general k; then, the MSPE, bias2, and variance for x∈a,b are equal to those for x∈a,a+b−a/k since the decision tree is assumed to make k equal terminal nodes with the same number of observations in each terminal node. In that case, for x∈a,b for a general k, the MSPE is(14)MSPE=1+kNσ2+β2b−a212k2,with variance as(15)EVar=kNσ2+β2b−a212k2and squared bias as(16)EBias2=β2b−a212k2.It is easy to see that with a lowerβ, b−a, and σ2 and higher N, variance, squared bias, and MSPE will all decrease.


### 2.3. Optimalk to Minimize MSPE

The ideal number of terminal nodes can be found by minimizing the MSPE with respect tok. Here k is a discrete integer, so the target k will be the nearest integer from the differentiate result. Calculating the first derivative of MSPE, we get(17)dMSPEkdk=σ2N−β2b−a26k3−β2b−a212Nk2,and the second derivative of MSPE is always positive. Therefore, we only need to solve(18)dMSPEkdk=0,k∈1,N.The real root of (18) is(19)kmin=β2b−a212σ23⋅N+N2−β2b−a2324σ23+N−N2−β2b−a2324σ23.Having(20)N≫βb−a18σ,we can approximate kmin by(21)kmin≈β2b−a2N6σ2.3In addition, the constraint for rootk is also kmin∈1,N. If kmin is not in 1,N, MSPE might always decrease.By substitutingkmin in (19) back into (16), we will get(22)EBias2=4812σ2βb−aN2/3,and it is easy to see that, with the increase of σ and βb−a, when N is fixed, EBias2 will increase. The others will be shown as figures.Accordingly, how will the ratiosEVar/MSPE, EBias2/MSPE, σ2/MSPE vary when parameters change? Since a, b, and β appear together, they are regarded as one parameter. For b and a, the thing that matters is their difference, so we use a=0 and only change b. Here, k is set to be kmin calculated using given parameters for (19), and if kmin does not exist, the results will not be shown. The results in Figure 1 (changing β2b−a2) and Figure 2 (changing σ2) show that, under both circumstances, MSPE, EVar, and EBias2 all increase.Figure 1
RatiosEVar/MSPE, EBias2/MSPE, σ2/MSPE with different β2b−a2. N=100 and α=0.
(a)(b)(c)Figure 2
RatiosEVar/MSPE, EBias2/MSPE, σ2/MSPE with different σ2. N=100 and α=0.
(a)(b)(c)In Figure1, when β2b−a2 gets bigger, X is more likely to be uniformly distributed and kmin increases as y is more accurately described with a uniform distribution; besides, the ratio of Var and Bias2 over MSPE gets larger while σ2 increases. In Figure 2, when σ2 gets bigger, the Gaussian distribution will play a bigger role in data generation and kmin decreases. That is why σ2/MSPE increases. EVar/MSPE and EBias2/MSPE generally decrease. The decrease speed slows with bigger b and β as expected.


### 2.4. Simulation

In this simulation, a simplified tree model will be designed to confirm the theory results using simulated data. That is, when parameters of the simulated data change, the distribution ofX and y will also change. The question is, how will the statistics of Var, Bias2, MSE, and kmin change accordingly?In the simplified tree,X is evenly split into k intervals, i=1,2,…,k. For specific k, a, b, N, α, and β, we are going to calculate the statistics of MSPE, Var, and Bias2 for the ith interval in k from simulated data. Thus, for the ith interval, the x range is(23)Ri=a+i−1b−aK,a+ib−aK.The number of observations in intervali (i=1,2,…k−1) is ni:(24)ni=N−∑j=0i−1njk−i,defining n0=0 and nk=N−∑j=0k−1nj.(i)
Step 1: for the data (x, y) in Ri, we train a model from them as(25)f^ix=y¯,for simulatedy, and y¯ is the averaged value of y in Ri.(ii)
Step 2: repeat Step 1s times. Then, we have s trained models f^jx, j=1,2,…,s.(iii)
Step 3: simulate onex0 uniformly from the x range Ri. We are going to calculate the SPEx0, Varx0, and Bias2x0 for this specific x0.(iv)
Step 4: simulates values of yj using x0.(v)
Step 5: calculate the statistics ofSPEx0, Varx0, and Bias2x0 for this specific x0 as(26)SPEx0=1s∑j=1sf^ix0−yj2,Bias2x0=1s∑j=1sf^ix0−fx02,Varx0=variancef^ix0.(vi)
Step 6: repeat Step 3 to Step 5 forn.repeat times and calculate the mean of SPEx0, Varx0, and Bias2x0 as MSPEi, Biasi2, and Vari.Follow Step 1 to Step 6 for alli, i=1,2,…,k, and calculate the mean as MSPE, Bias2, and Var.The results of simulations with 200 trials are shown in Figures3 and 4. For Figure 3, kmin∈1,N, we have a minimum MSPE. However, when kmin∉1,N as in Figure 4, MSPE keeps decreasing.Figure 3
An example whenkmin exists. The x axis label is k: the number of trees splits. MSPE, Bias2, and Var are averaged calculations from 200 simulation trials. The black line is the MSPE, the blue line is the Bias2, the orange line is σ2, and the red line is the Var. The solid lines are from simulated data, and the dashed lines are theoretical calculations. The parameters values are given as follows: α=2, β=4, N=100, a=0, b=3, and σ=4.
(a)(b)(c)(d)Figure 4
An example whenkmin does not exist. The x axis label is k: the number of trees splits. MSPE, Bias2, and Var are averaged calculations from 200 simulation trials. The black line is the MSPE, the blue line is the Bias2, the orange line is σ2, and the red line is the Var. The solid lines are from simulated data, and the dashed lines are theoretical calculations. The parameters values are given as follows: α=2, β=12, N=100, a=0, b=4, and σ=0.2.
(a)(b)(c)(d)


## 2.1. Decomposition Background

For theith observation Yi, the (unconditional) expectation is(3)EYi=Eα+βXi+ɛi=α+βEXi+Eɛi=α+β⋅a+b2,and the variance is(4)VarYi=Varα+βXi+ɛi=β2VarXi+Varɛi=β2⋅b−a212+σ2.They both have no relationship toi. In that case, EY=EYi and VarY=VarYi. Accordingly, for N observations, the expectation and variance for the average Y¯ are shown in(5)EY¯=1NEY1+EY2+⋯+EYN,EY=α+β⋅a+b2,(6)VarY¯=1N2VarY1+VarY2+⋯+VarYN=1NVarY=1Nβ2⋅b−a212+σ2.


## 2.2. Decomposition in the Context of Decision Trees

In the context of decision trees, the fitted model isf^X in a simplified form is(7)f^X=Y¯i,i=1,2,…,k,where k is the number of terminal nodes in the tree f^X and Y¯i is the mean of y in terminal node i. In a tree with only the root node, k=1, and the fitted model is f^X=Y¯. Then, for point x0,(8)Bias2f^x0=Ef^x0−fx02=α+β⋅a+b2−α−β⋅x02=β2a+b2−x02,and the variance is(9)Varf^x0=Ef^x0−Ef^x02=EY¯−α+β⋅a+b22=EY¯2+α+β⋅a+b22−2EY¯α+β⋅a+b2=E2Y¯+VarY¯−α+β⋅a+b22=σ2N+β2b−a212N.Thus, theSPE at point x0 is(10)SPEf^x0=σ2+bias2f^x0+Varf^x0=σ2+β2a+b2−x02+σ2N+β2b−a212N=1+1Nσ2+β2a+b2−x02+β2b−a212N.Then, the mean squared prediction error (MSPE) is(11)MSPE=∫abSPEf^xPX=xdx=∫abSPEf^x1b−adx=1+1Nσ2+β2b−a212N+β2∫aba+b2−x21b−adx=1+1Nσ2+β2b−a212N+β2112b−a2=1+1Nσ2+β2b−a212.Comprising variance is(12)EVar=1Nσ2+β2b−a212,and squared bias is(13)EBias2=β2b−a212.Now the number of terminal nodes in the decision tree is extended fromk=1 to a general k; then, the MSPE, bias2, and variance for x∈a,b are equal to those for x∈a,a+b−a/k since the decision tree is assumed to make k equal terminal nodes with the same number of observations in each terminal node. In that case, for x∈a,b for a general k, the MSPE is(14)MSPE=1+kNσ2+β2b−a212k2,with variance as(15)EVar=kNσ2+β2b−a212k2and squared bias as(16)EBias2=β2b−a212k2.It is easy to see that with a lowerβ, b−a, and σ2 and higher N, variance, squared bias, and MSPE will all decrease.


## 2.3. Optimalk to Minimize MSPE

The ideal number of terminal nodes can be found by minimizing the MSPE with respect tok. Here k is a discrete integer, so the target k will be the nearest integer from the differentiate result. Calculating the first derivative of MSPE, we get(17)dMSPEkdk=σ2N−β2b−a26k3−β2b−a212Nk2,and the second derivative of MSPE is always positive. Therefore, we only need to solve(18)dMSPEkdk=0,k∈1,N.The real root of (18) is(19)kmin=β2b−a212σ23⋅N+N2−β2b−a2324σ23+N−N2−β2b−a2324σ23.Having(20)N≫βb−a18σ,we can approximate kmin by(21)kmin≈β2b−a2N6σ2.3In addition, the constraint for rootk is also kmin∈1,N. If kmin is not in 1,N, MSPE might always decrease.By substitutingkmin in (19) back into (16), we will get(22)EBias2=4812σ2βb−aN2/3,and it is easy to see that, with the increase of σ and βb−a, when N is fixed, EBias2 will increase. The others will be shown as figures.Accordingly, how will the ratiosEVar/MSPE, EBias2/MSPE, σ2/MSPE vary when parameters change? Since a, b, and β appear together, they are regarded as one parameter. For b and a, the thing that matters is their difference, so we use a=0 and only change b. Here, k is set to be kmin calculated using given parameters for (19), and if kmin does not exist, the results will not be shown. The results in Figure 1 (changing β2b−a2) and Figure 2 (changing σ2) show that, under both circumstances, MSPE, EVar, and EBias2 all increase.Figure 1
RatiosEVar/MSPE, EBias2/MSPE, σ2/MSPE with different β2b−a2. N=100 and α=0.
(a)(b)(c)Figure 2
RatiosEVar/MSPE, EBias2/MSPE, σ2/MSPE with different σ2. N=100 and α=0.
(a)(b)(c)In Figure1, when β2b−a2 gets bigger, X is more likely to be uniformly distributed and kmin increases as y is more accurately described with a uniform distribution; besides, the ratio of Var and Bias2 over MSPE gets larger while σ2 increases. In Figure 2, when σ2 gets bigger, the Gaussian distribution will play a bigger role in data generation and kmin decreases. That is why σ2/MSPE increases. EVar/MSPE and EBias2/MSPE generally decrease. The decrease speed slows with bigger b and β as expected.


## 2.4. Simulation

In this simulation, a simplified tree model will be designed to confirm the theory results using simulated data. That is, when parameters of the simulated data change, the distribution ofX and y will also change. The question is, how will the statistics of Var, Bias2, MSE, and kmin change accordingly?In the simplified tree,X is evenly split into k intervals, i=1,2,…,k. For specific k, a, b, N, α, and β, we are going to calculate the statistics of MSPE, Var, and Bias2 for the ith interval in k from simulated data. Thus, for the ith interval, the x range is(23)Ri=a+i−1b−aK,a+ib−aK.The number of observations in intervali (i=1,2,…k−1) is ni:(24)ni=N−∑j=0i−1njk−i,defining n0=0 and nk=N−∑j=0k−1nj.(i)
Step 1: for the data (x, y) in Ri, we train a model from them as(25)f^ix=y¯,for simulatedy, and y¯ is the averaged value of y in Ri.(ii)
Step 2: repeat Step 1s times. Then, we have s trained models f^jx, j=1,2,…,s.(iii)
Step 3: simulate onex0 uniformly from the x range Ri. We are going to calculate the SPEx0, Varx0, and Bias2x0 for this specific x0.(iv)
Step 4: simulates values of yj using x0.(v)
Step 5: calculate the statistics ofSPEx0, Varx0, and Bias2x0 for this specific x0 as(26)SPEx0=1s∑j=1sf^ix0−yj2,Bias2x0=1s∑j=1sf^ix0−fx02,Varx0=variancef^ix0.(vi)
Step 6: repeat Step 3 to Step 5 forn.repeat times and calculate the mean of SPEx0, Varx0, and Bias2x0 as MSPEi, Biasi2, and Vari.Follow Step 1 to Step 6 for alli, i=1,2,…,k, and calculate the mean as MSPE, Bias2, and Var.The results of simulations with 200 trials are shown in Figures3 and 4. For Figure 3, kmin∈1,N, we have a minimum MSPE. However, when kmin∉1,N as in Figure 4, MSPE keeps decreasing.Figure 3
An example whenkmin exists. The x axis label is k: the number of trees splits. MSPE, Bias2, and Var are averaged calculations from 200 simulation trials. The black line is the MSPE, the blue line is the Bias2, the orange line is σ2, and the red line is the Var. The solid lines are from simulated data, and the dashed lines are theoretical calculations. The parameters values are given as follows: α=2, β=4, N=100, a=0, b=3, and σ=4.
(a)(b)(c)(d)Figure 4
An example whenkmin does not exist. The x axis label is k: the number of trees splits. MSPE, Bias2, and Var are averaged calculations from 200 simulation trials. The black line is the MSPE, the blue line is the Bias2, the orange line is σ2, and the red line is the Var. The solid lines are from simulated data, and the dashed lines are theoretical calculations. The parameters values are given as follows: α=2, β=12, N=100, a=0, b=4, and σ=0.2.
(a)(b)(c)(d)


## 3. Prediction Interval

Instead of point prediction, a prediction interval is also desirable especially for time series with high variance. If both the point prediction and the prediction interval can be provided, we will be more confident for the prediction. This study also helps us decide the proper prediction interval method for decision-tree-based regression problems. Gaussian-based prediction interval and quantile interval are compared under different parameters distributions.


### 3.1. Probability Function ofY

Since our linear model,(27)Y=α+βX+ɛ,is the sum of uniform and Gaussian distributions, the probability function for Y is(28)PYy=∫abPY=y|X=xPXxdx=∫ab12πσexp−y−α+βx22σ21b−adx.By lettingt=α+βx−y/σ, we obtain(29)PYy=1βb−a∫α+βa−y/σα+βb−y/σ12πexp−t22dt=1βb−aΦα+βb−yσ−Φα+βa−yσ.Now we get the probability ofY as (29). However, PYy is in a complex form meaning that the parameters are not easily solvable in theory by a given value for PYy.


### 3.2. Prediction Interval as a Gaussian Distribution

If we want to get the prediction interval, sayy1,y2 for Y at 1−p level, the theoretical way is to obtain y1 and y2 from the equations(30)∫−∞y1PYydy=p2,∫y2∞PYydy=1−p2.However, the integral ofΦ is not analytically solvable without approximating Φ with other suitable expressions. The results will also be quite complex. If we know the parameters values, then y1 and y2 can easily be found numerically.From Figure5, if the uniform (Gaussian) distribution plays a main role, then Y can be approximately described by a uniform (Gaussian) distribution. Under the conditions that β is not too large, σ is not too small, and k is 1 (with only one interval), we will approximate the distribution of Y as a Gaussian distribution NμY,σY2:(31)Y∼Nα+β⋅a+b2,β2⋅b−a212+σ2.Figure 5
Probability density function ofY with different parameter values. The red line is the theoretical probability, and the black line is the simulated probability.
(a)(b)Then, the prediction interval under95% criteria for this Gaussian distribution is around(32)Y¯−1.961+1Nβ2⋅b−a212+σ2,Y¯+1.961+1Nβ2⋅b−a212+σ2.Then, for a generalk, the prediction interval becomes(33)f^−1.96×1+kNβ2⋅b−a212k2+σ2,f^+1.96×1+kNβ2⋅b−a212k2+σ2,which is the form(34)f^−1.96×RMSPE,f^+1.96×RMSPE,a typical Gaussian prediction interval.


### 3.3. Prediction Simulation Using Gaussian Prediction Interval and Quantile Interval

In this simulation, we explore the performance of Gaussian prediction intervals and quantile intervals under different parameter combinations. The parameters includeσ, b−a, β, and k. When the other parameters are fixed, a higher σ means a stronger Gaussian distribution effect, in which case, Gaussian prediction interval may work well. When β2b−a2 is large, the uniform distribution plays a bigger role. Then, Gaussian prediction interval may not work so well. Both Gaussian prediction interval and quantile interval are influenced by the observation size of the terminal node. When the sample size is large, they can have stable performance, but when sample size is small, performance differs.The Gaussian prediction interval in use is(35)f^−cRMSPE,f^+cRMSPE,where c is 1.96 and RMSPE is the root mean squared error estimated from the training data in each terminal node.The quantile intervalL,U comes from the 0.025 and 0.975 quantiles of each terminal node from the training data.(i)
Step 1: training data generation.Using given parametersα, β, a, b, σ, N, data are generated according to the model(36)Y=α+βX+ɛ.Therefore, we get the true fitted values forY.(ii)
Step 2: RMSPE and quantiles from training data.From this training data, the trained model, RMSPE, and quantiles are calculated as in the following steps.(i)
Step 2.1: model training.For training dataA (the rest of data B is the test data), we sort the data x in an ascending order, so y will also be rearranged following x, and then Atraining is divided into k roughly successive equal folds, making a total of N observations. The number of observations in fold i (i=1,2,…k) is ni:(37)ni=N−∑j=0i−1njk−i,definingn0=0.For theith fold in A, giving xi and yi, the predicted value will be(38)y^=y¯i,$x∈minxi,maxxi$,in the trees context. The predicted value of a tree model is the averaged response values of each terminal node. Samples being split into those terminal nodes will have the corresponding averaged value as the predicted value.(ii)
Step 2.2: RMSPE and quantile calculation.When the model for eachi is trained as modelAi, the predicted values for y in A will be y^. Then, the RMSPE for the training data is(39)RMSPE=∑r=1Ny−y^2N.The quantile intervalsL and U are the 0.025 and 0.975 quantiles of the ith training data y.(iii)
Step 3: test data generation and model testing.Using the same parametersα, β, a, b, and σN as in Step 1, data are generated according to(40)Y=α+βX+ɛ.Then, the test dataB are put into modelA and the coverage is computed as(41)1N∑r=1NIy^−cRMSPE<y<y^+cRMSPE.(iv)
Step 4: repeating Steps 1 to 3.Repeat Steps 1 to 3s times to get an averaged coverage.Using parametersa=0, α=2, and s=200, the results are shown in Figures 6.Figure 6
Prediction coverage using RMSPE. The black dash line is 0.95. The dash lines are coverage. The solid lines are width. The red lines represent Gaussian prediction interval, and the blue lines represent quantile prediction interval. The right-handy axis (purple) is the axis for width. N=10000.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)The results show that quantile interval coverages are closer to the 0.95 reference line for fixedσ, b, and β. Gaussian prediction interval is only closer to the 0.95 coverage when σ is large; otherwise, it is larger than 0.95 at the cost of wider width. When k is chosen as the best kmin, the coverages get closer to the 0.95 reference line as σ increases for both quantile and Gaussian prediction intervals. However, when the uniform distribution effect gets stronger, the coverages all go far away from 0.95. Accordingly, when the number of observations for each terminal node is large and the data distribution is not obviously Gaussian, quantile intervals are suggested. When the data follows obvious Gaussian distribution, Gaussian prediction intervals are recommended.


## 3.1. Probability Function ofY

Since our linear model,(27)Y=α+βX+ɛ,is the sum of uniform and Gaussian distributions, the probability function for Y is(28)PYy=∫abPY=y|X=xPXxdx=∫ab12πσexp−y−α+βx22σ21b−adx.By lettingt=α+βx−y/σ, we obtain(29)PYy=1βb−a∫α+βa−y/σα+βb−y/σ12πexp−t22dt=1βb−aΦα+βb−yσ−Φα+βa−yσ.Now we get the probability ofY as (29). However, PYy is in a complex form meaning that the parameters are not easily solvable in theory by a given value for PYy.


## 3.2. Prediction Interval as a Gaussian Distribution

If we want to get the prediction interval, sayy1,y2 for Y at 1−p level, the theoretical way is to obtain y1 and y2 from the equations(30)∫−∞y1PYydy=p2,∫y2∞PYydy=1−p2.However, the integral ofΦ is not analytically solvable without approximating Φ with other suitable expressions. The results will also be quite complex. If we know the parameters values, then y1 and y2 can easily be found numerically.From Figure5, if the uniform (Gaussian) distribution plays a main role, then Y can be approximately described by a uniform (Gaussian) distribution. Under the conditions that β is not too large, σ is not too small, and k is 1 (with only one interval), we will approximate the distribution of Y as a Gaussian distribution NμY,σY2:(31)Y∼Nα+β⋅a+b2,β2⋅b−a212+σ2.Figure 5
Probability density function ofY with different parameter values. The red line is the theoretical probability, and the black line is the simulated probability.
(a)(b)Then, the prediction interval under95% criteria for this Gaussian distribution is around(32)Y¯−1.961+1Nβ2⋅b−a212+σ2,Y¯+1.961+1Nβ2⋅b−a212+σ2.Then, for a generalk, the prediction interval becomes(33)f^−1.96×1+kNβ2⋅b−a212k2+σ2,f^+1.96×1+kNβ2⋅b−a212k2+σ2,which is the form(34)f^−1.96×RMSPE,f^+1.96×RMSPE,a typical Gaussian prediction interval.


## 3.3. Prediction Simulation Using Gaussian Prediction Interval and Quantile Interval

In this simulation, we explore the performance of Gaussian prediction intervals and quantile intervals under different parameter combinations. The parameters includeσ, b−a, β, and k. When the other parameters are fixed, a higher σ means a stronger Gaussian distribution effect, in which case, Gaussian prediction interval may work well. When β2b−a2 is large, the uniform distribution plays a bigger role. Then, Gaussian prediction interval may not work so well. Both Gaussian prediction interval and quantile interval are influenced by the observation size of the terminal node. When the sample size is large, they can have stable performance, but when sample size is small, performance differs.The Gaussian prediction interval in use is(35)f^−cRMSPE,f^+cRMSPE,where c is 1.96 and RMSPE is the root mean squared error estimated from the training data in each terminal node.The quantile intervalL,U comes from the 0.025 and 0.975 quantiles of each terminal node from the training data.(i)
Step 1: training data generation.Using given parametersα, β, a, b, σ, N, data are generated according to the model(36)Y=α+βX+ɛ.Therefore, we get the true fitted values forY.(ii)
Step 2: RMSPE and quantiles from training data.From this training data, the trained model, RMSPE, and quantiles are calculated as in the following steps.(i)
Step 2.1: model training.For training dataA (the rest of data B is the test data), we sort the data x in an ascending order, so y will also be rearranged following x, and then Atraining is divided into k roughly successive equal folds, making a total of N observations. The number of observations in fold i (i=1,2,…k) is ni:(37)ni=N−∑j=0i−1njk−i,definingn0=0.For theith fold in A, giving xi and yi, the predicted value will be(38)y^=y¯i,$x∈minxi,maxxi$,in the trees context. The predicted value of a tree model is the averaged response values of each terminal node. Samples being split into those terminal nodes will have the corresponding averaged value as the predicted value.(ii)
Step 2.2: RMSPE and quantile calculation.When the model for eachi is trained as modelAi, the predicted values for y in A will be y^. Then, the RMSPE for the training data is(39)RMSPE=∑r=1Ny−y^2N.The quantile intervalsL and U are the 0.025 and 0.975 quantiles of the ith training data y.(iii)
Step 3: test data generation and model testing.Using the same parametersα, β, a, b, and σN as in Step 1, data are generated according to(40)Y=α+βX+ɛ.Then, the test dataB are put into modelA and the coverage is computed as(41)1N∑r=1NIy^−cRMSPE<y<y^+cRMSPE.(iv)
Step 4: repeating Steps 1 to 3.Repeat Steps 1 to 3s times to get an averaged coverage.Using parametersa=0, α=2, and s=200, the results are shown in Figures 6.Figure 6
Prediction coverage using RMSPE. The black dash line is 0.95. The dash lines are coverage. The solid lines are width. The red lines represent Gaussian prediction interval, and the blue lines represent quantile prediction interval. The right-handy axis (purple) is the axis for width. N=10000.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)The results show that quantile interval coverages are closer to the 0.95 reference line for fixedσ, b, and β. Gaussian prediction interval is only closer to the 0.95 coverage when σ is large; otherwise, it is larger than 0.95 at the cost of wider width. When k is chosen as the best kmin, the coverages get closer to the 0.95 reference line as σ increases for both quantile and Gaussian prediction intervals. However, when the uniform distribution effect gets stronger, the coverages all go far away from 0.95. Accordingly, when the number of observations for each terminal node is large and the data distribution is not obviously Gaussian, quantile intervals are suggested. When the data follows obvious Gaussian distribution, Gaussian prediction intervals are recommended.


## 4. Real Application

We have explored the performance of decision trees under different circumstances. A real application is conducted in this section. The data come from UK Gridwatch (http://www.gridwatch.templar.co.uk/), which are the demand data of grid and the supply data of each energy source. The time series begin from year 2011 to year 2020, making a total of 953824 observations with a record every 5 minutes. The details are shown in Figure 7.Figure 7
The supply and demand of energy in UK from year 2011 to year 2020. The time gap is 5 minutes.From the figure, we can see that the demand of grid changes in period as expected since there are peak and valley values daily and seasonally. The general trend of grid demand changes a little. Some kinds of energy like wind and biomass increase a lot in supply these years; they will be more frequently used in the future than traditional energy like coal as they are more environmentally friendly. We construct a metricratio to measure the ratio of other energy supply over that of coal. By deleting observations which have none or zero values of coal, we have 847922 observations left, as shown in Figure 8.Figure 8
The ratio of other energy supply over that of coal.We average the time seriesratio from a frequency of 5 minutes to a daily basis, ending with 2954 observations left. A forecasting method is conducted on ratio to help us know how much renewable energy is needed in the near future. The interval forecasting method we use is from our designed method, Zhao et al. [14], called ctreeone, which uses the tree method ctree in a dynamic interval forecasting context. The different parameter we choose is 7 for time gap (weekly dynamic forecasting), leaving the other parameters unchanged.Interval forecasting provides not only the point forecasting results but also the prediction interval that the predicted point belongs to. Small change of theratio often happens, which influences a little the the energy supply and demand system, so no action is needed in this circumstance. When the predicted ratio changes a lot, out of a preset limit, an alarm may be raised to help the system accommodate to the new circumstance, for example, by producing more renewable energies in advance to meet the instant demand. The interval forecasting model provides such an alerting system to adjust the energy production.The results are shown in Table1 and Figure 9. The coverage and width make a good balance; that is, a higher coverage costs a relatively higher width. We end with a coverage of 80.31% and a suitable width of 22.95 which is similar to the standard deviation of ratio of 19.78.Table 1
Results of interval forecasting on the time seriesratio.
Coverage (%)Averaged widthComputational time (s)80.3122.9511.86The experiment is run on a 3.2 GHz 8-Core Intel Xeon W processor.Figure 9
The dynamic interval forecasting result. The black lines are the raw daily ratio time series; the red lines are the predicted values; the blue line is the upper forecasted interval; and the yellow line is the lower forecasted interval. The purple line describes whether the data are coved by the interval or not. The circle means the model is retrained on that time point. The upper and lower forecasted intervalU,L give an interval that in most cases the real future value will be in, which is similar to the 95% confidence interval that the fitted value is most likely to be covered in.


## 5. Conclusion

In this paper, the data are constructed using a simple model that includes both Gaussian and uniform distributions. We explore the squared prediction error in the context of trees and decompose that error into bias, variance, and irreducible error. The bias decreases when the tree gets bigger. However, for squared prediction error and variance, the relationship is not monotonic. We also calculate the best tree depth with a minimum mean squared prediction error. When Gaussian effect dominates, the best tree depth density decreases. However, when uniform effect dominates, the best tree depth increases. Under both circumstances, mean squared error, variance, and bias all increase.After that, two options are given for the prediction interval using Gaussian prediction interval or quantile interval. When Gaussian distribution is obviously dominant, Gaussian prediction intervals are suggested. Otherwise, quantile intervals are suggested, which is also why quantile intervals are chosen as the prediction intervals in our regression application, although they both perform poorly when the uniform distribution is quite strong. When the number of observations is small in the terminal node, both interval constructions perform poorly in terms of coverage.In the real data application, we applied our method to the UK grid energy supply and demand data to forecast the ratio of renewable energy supply over that of coal. We have good forecasting results as 80.31% in interval coverage and 22.95 in interval width. The method can be extended to other models as well besides decision trees.We use the model decision tree for interval forecasting. In practice, other models can also be considered. For example, Hall et al. [22] used multiple nonlinear regression to forecast and analyze the changes of climate and weather dynamics and proposed a simple model averaging approach to reduce model and prediction uncertainty. Besides decision tree, other dynamic regression models can also be considered, for example, Gu et al. [23] used dynamic regression model to predict the dynamics of a specific space weather index and proposed a new approach for prediction uncertainty analysis using point-cloud model parameters. Dynamic regression model was also applied to social dynamic behavior modeling and analysis [24].In the future research, the model can be applied to more kinds of datasets to test its generation ability. In the simulation, instead of linear model, nonlinear model can also be considered to test the tree performance.



---

*Source: 3567894-2020-10-26.xml*