A hybridization of MODWT-SVR-DE model emphasizing on noise reduction and optimal parameter selection for prediction of CO₂ emission in Thailand

Abstract

In present, Thailand is a major manufacturing base in Southeast Asia, which tremendously produces many sources of greenhouse gases (GHGs). The emission CO₂ become a spotlight issue since Thailand legislated for controlling industrial carbon footprint, which resulted in carbon credit charged. In order to control and monitor those effects of GHGs, future accuracy of CO₂ emission and CO₂ emission equivalence play a significant role and can support critical making on suitable control and monitoring. In this research, a hybridization of MODWT-SVR-DE model is developed and proposed that employs maximal overlap discrete wavelet transform (MODWT) with first decomposition to reduce fluctuation of CO₂ emission and CO₂ emission equivalence. Afterward, the support vector regression is used to formulate complex model. Meanwhile, differential evolution is used to search given parameters of support vector regression. Moreover, the proposed model is compared to conventional forecasting models (i.e. ARIMA, Holt and simple exponential smoothing [SES]) and hybrid model of support vector regression and differential evolution. The empirical results indicated that the proposed model outperforms all candidate models and provides significant difference than candidate models at 0.05 significance levels. Consequently, the proposed model is able to be accurately applied in order to monitor the environment and support policy planners to take step with the helpful guideline.

Keywords:

• Noise reduction
• time series analysis
• optimal parameter selection
• environment monitoring

Reviewing Editor:

• Zhou Zude
• Wuhan University of Technology
• China

Subjects:

• Science
• Environmental Studies & Management
• Environment & Business
• Mathematics & Statistics
• Advanced Mathematics
• Applied Mathematics
• Statistics & Probability

1. Introduction

Nowadays, all people around the world are aware and concerned about climate change (Franco et al., 2022; Lu et al., 2022; Kim & Jang, 2022; Sovacool et al., 2021), which is caused from greenhouse gas (GHG) in forms of CO₂ emission and CO₂ emission equivalence (Raihan & Voumik, 2022; Batool et al., 2022; Raihan & Tuspekova, 2022; Nong et al., 2021; Yoro & Daramola, 2020). Consequently, all forms of CO₂ emission and CO₂ emission equivalence are limited and forced by law. For Thailand (Lu et al., 2022; Wongsapai & Daroon, 2021; Srikaummun et al., 2021; Jaiboon et al., 2021), many industrial estates have been rapidly growing, which tremendously produce a lot of GHG. However, amount of CO₂ emission and CO₂ emission equivalence in a form of GHG is difficult to timely measure. Moreover, the amount of GHG in a form of CO₂ emission and CO₂ emission equivalence is uncertainty in each period.

In order to monitor and control risks of greenhouse effect, future amount of GHG in a form of CO₂ emission and CO₂ emission equivalence is need to realize before making critical decision on reduce risks of greenhouse effect (Franco et al., 2022).

With the aim of extrapolating amount of CO₂ emission and CO₂ emission equivalence in a form of GHG in advance, many forecasting techniques have been developed and proposed to project future amount of CO₂ emission and CO₂ emission equivalence in recent years. One of the most popular techniques is univariate time series analysis, which is formulated from only previous observations in order to form the prediction function. Therefore, the developed technique is a convenient tool dealing with predicting CO₂ emission and CO₂ emission equivalence in situations of reality. The summary list of abbreviations and acronyms used in this research is presented in Table 1.

Table 1. List of abbreviations and acronyms is used in this research.

Acronyms	Forecasting method
ARIMA	Autoregressive integrated moving average
MAE	Mean absolute error
MAPE	Mean absolute percentage error
MODWT	Maximal overlap discrete wavelet transform
MODWT-SVR-DE	Hybrid model of maximal overlap discrete wavelet transform, support vector regression and differential evolution
RMSE	Root mean square error
RMSPE	Root mean square percentage error
SES	Simple exponential smoothing
sMAPE	Symmetric mean absolute error
SVR-DE	Hybrid model of support vector regression and differential evolution

In present, several conventional approaches of time series are still widely used to predict future amount of CO₂ emission and CO₂ emission equivalence. For instances, autoregressive integrated moving average (ARIMA) approaches (Cahyono et al., 2022; Javanmard & Ghaderi, 2022; Kumari & Singh, 2022; Xu et al., 2021; Ning et al., 2021; Leerbeck et al., 2020) provide well performance in linear forecasting problems while exponential smoothing without more complex approaches can provide robust of forecast in many literatures of science (Franco et al., 2022; Kumari & Singh, 2022). However, individual approach may not be sufficient to provide well performance in all situations of reality.

For being in need of more accuracy in present, many combined approaches are developed and proposed to achieve those situations. The combined approaches take advantage of a forecasting technique to reduce disadvantage of one another. A combined approach emphasizing on optimal parameter selection is a well-established approach and is easy to find in literatures surrounding many problems of time series forecasting. One of the powerful hybridizations is a combined approach of support vector regression and differential evolution, namely SVR-DE model, which utilizes support vector regression to build complex and adjustable functions to describe time series pattern whether linear or nonlinear function (Li et al., 2022; Ehteram et al., 2021; Ahmadi et al., 2019). The advantages of nonlinearity and complexity can provide better performance than that of conventional approaches of statistical technique (e.g. exponential smoothing and ARIMA) in real-world complex problems. Since statistical techniques are based on linear formulation and many prior assumptions, then these approaches are not adjustable to manipulate model and cannot also explain nonlinear relation of time series. Furthermore, support vector regression is formulated from principle of structural risk minimization, which provides global optimal solution due to convex problem. Although support vector regression performs well in forecasting, its performance depends heavily on appropriate parameter selection (Javanmard & Ghaderi, 2022; Adnan et al., 2022; Ehteram et al., 2021; Ghazvini et al., 2020). Consequently, the differential evolution (Hou et al., 2018; Hu & Chen, 2018; Zhang et al., 2016), a recent evolution algorithm as effective global optimization algorithm that has many advantages with simplicity, robustness, and reliability of implementation, is used to search for the most proper parameters of support vector regression within given search space. In addition, the combined approach of support vector regression and differential evolution can achieve more accuracy than traditional forecasting approaches. However, the combined approach may not be efficient and successful in all situations of actuality with regard to many noise problems.

In this regard, a combined approach including data pre-processing is introduced to address those problems, which emphasizes on preliminary process on datasets by decomposing time series into more stationary and regular subseries. The subseries are generally clearer to analyze with removing the irrelevant and redundant components of time series. Although data pre-processing techniques are a useful tool, model manipulations in each subseries are quite difficult as well. Therefore, this combined approach uses benefit of data pre-processing techniques emphasizing only on first level of decomposition to overcome the limitation of support vector regression concerning noise problems in time series analysis. One of the several data pre-processing techniques is maximal overlap discrete wavelet transform (MODWT), which is a well-known technique of modified version of discrete wavelet transform due to capability to address the circular shift effect (Yang et al., 2022; Shams et al., 2021; Seo et al., 2017).

In this research, a hybridization of MODWT-SVR-DE model is developed and proposed to forecast amount of CO₂ emission source of several sectors. With the intention of covering both perspectives, data pre-processing and optimal parameter selection techniques are developed and proposed to achieve more accuracy and precision of forecast. The proposed model is evaluated its performance of forecasting based on five accuracy measures and is compared to conventional models that are ARIMA, Holt’s model and simple exponential smoothing (SES). Moreover, the proposed model is compares to hybridization of SVR-DE model based on five accuracy measures as well.

For accuracy measures, both scale-dependent and scale-independent metrics are mean absolute error (MAE), root mean square error (RMSE), mean absolute percentage error (MAPE), symmetric mean absolute error (sMAPE) and root mean square percentage error (RMSPE). Those measures are still preferred to use in many literatures. Furthermore, scale-independent metrics (e.g. MAPE, sMAPE and RMSPE) are adopted with Friedman test and post hoc test to identify significant difference between forecasting performances.

2. Material and methodologies

2.1. Datasets of CO₂ emission and CO₂ emission equivalence

The datasets of CO₂ emission and CO₂ emission equivalence are categorized into several sectors and are secondary data of open online database from Climate Watch (2022). Washington, D.C.: World Resource Institute based on data source of Climate Analysis Indicators Tool (CAIT), which are time series datasets from 1990 to 2019. The summary of datasets is demonstrated in Figure 1.

Figure 1. Summary of CO₂ emission and CO₂ emission equivalence datasets.

2.2. Simple exponential smoothing

The simplest form of exponential smoothing is SES, which is a weighted average of previous observations with exponentially decreasing weights over time from present observation. The SES is suitable for time series dataset with no clear trend or seasonal pattern. The mathematical formulation is demonstrated as Equations (1)(1) $$l_{t-1}=\alpha y_{t-1}+\left(1-\alpha \right)l_{t-2}$$(1) and (2)(2) $$y_t=l_{t-1}+{\epsilon }_t$$(2) . (1) $$l_{t-1}=\alpha y_{t-1}+\left(1-\alpha \right)l_{t-2}$$(1) (2) $$y_t=l_{t-1}+{\epsilon }_t$$(2) where $y_t$ and ${\epsilon }_t$ are actual observation and error at time t, respectively. $l_t$ is the level or smoothed value of time series at time t. $\alpha$ is smoothing parameter between 0 and 1.

For establishing appropriate model of SES, the fitted model is evaluated by using SES function of R programing language (Hyndman et al., 2022) based on criterion of the lowest mean square error to provide proper smoothing parameter and initial estimated observation of first observation of time series.

2.3. Holt’s model

With regard to trend pattern using exponential smoothing, Holt’s model is a form of exponential smoothing with attracting equation of trend to the SES. Hence, the mathematical form is able to describe dataset with trend and is demonstrated as Equations (3)–(5). (3) $$l_{t-1}=\alpha y_{t-1}+\left(1-\alpha \right)\left(l_{t-2}+b_{t-2}\right)$$(3) (4) $$b_{t-1}=\beta \left(l_{t-1}-l_{t-2}\right)+\left(1-\beta \right)b_{t-2}$$(4) (5) $$y_t= l_{t-1}+b_{t-1}+{\epsilon }_t$$(5) where $y_t$ and ${\epsilon }_t$ are actual observation and error at time t. The $l_t$ and $b_t$ are the estimate of level and trend of the time series at time t, respectively. $\alpha$ and $\beta$ are smoothing parameter for level and trend between 0 and 1, respectively.

A holt function of R programming (Hyndman et al., 2022) is used to formulate the fitted model based on the lowest mean square error and provides suitable smoothing parameters and also estimated values of level and trend.

2.4. Autoregressive integrated moving average

ARIMA is well-known approach for linear forecasting problems, which is the general form of autoregressive moving average (ARMA) with differencing term to transform non-stationary time series to stationary series in order to estimate parameters of ARIMA model. The mathematical approach with mean $\mu$ is described as Equation (6)(6) $$\left(1-\sum _{i=1}^p{\varphi }_i{B}^i\right){\left(1-B\right)}^d\left(y_t-\mu \right)=\left(1-\sum _{j=1}^q{\theta }_j{B}^j\right){\epsilon }_t$$(6) . (6) $$\left(1-\sum _{i=1}^p{\varphi }_i{B}^i\right){\left(1-B\right)}^d\left(y_t-\mu \right)=\left(1-\sum _{j=1}^q{\theta }_j{B}^j\right){\epsilon }_t$$(6) where $y_t$ and ${\epsilon }_t$ is the actual observation and error at time t, respectively. $p$ and $q$ are the order of autoregressive part and moving average part, respectively. $d$ is the order of differencing part, which is integer as well.

Since several ARIMA models can mimic a time series dataset, then the most suitable ARIMA model is chosen from the lowest Akaike information criterion for small sample size. The fitted model is received from auto.arima function of R programming language (Hyndman & Khandakar, 2008).

2.5. Support vector regression

The support vector regression is an extended form of support vector machine to deal with regression problems, which mathematical expression is presented as Equation (7)(7) $$f\left(x\right)=\sum _{i=1}^T\left({\alpha }_i-{\alpha }_i^{\ast }\right)K\left(x,x_i\right)+b$$(7) . (7) $$f\left(x\right)=\sum _{i=1}^T\left({\alpha }_i-{\alpha }_i^{\ast }\right)K\left(x,x_i\right)+b$$(7) where ${\alpha }_i$ and ${\alpha }_i^{\ast }$ are the so-called Lagrange multipliers. $(\cdot ,\cdot )$ denotes vector inner product. $\left(\cdot ,\cdot \right)$ is a scalar threshold; $K(x,x_i)$ is kernel function.

In this research, the types of kernel function satisfying Mercer’s condition are defined as Equations (8)(8) $$\text{Linear}:\hspace{1em}K\left(x,x_i\right)=x^T{x}_i$$(8) and (9)(9) $$\text{Radial basis}:\hspace{1em}K\left(x,x_i\right)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(-\gamma {\Vert x-x_i\Vert }^2\right)$$(9) . (8) $$\text{Linear}:\hspace{1em}K\left(x,x_i\right)=x^T{x}_i$$(8) (9) $$\text{Radial basis}:\hspace{1em}K\left(x,x_i\right)=\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\left(-\gamma {\Vert x-x_i\Vert }^2\right)$$(9)

2.6. Combined approach of support vector regression and differential evolution

The combined approach is developed with selecting proper parameters of support vector regression, which are input lag and also hyperplane based on advantage of differential evolution to investigate the appropriate parameters within given search space In other word, the combined approach can reduce a risk of using improper parameters of support vector regression and can reduce computational time consumption to model prediction function. The algorithm of the SVR-DE model is demonstrated as follows:

2.6.1. The differential evolution generates initial parameters of support vector regression dealing with input lag and hyperplane as described in Equation (10)(10) $${\theta }_{i,G}=\left[{\theta }_{1,i,G},{\theta }_{2,i,G},{\theta }_{3,i,G},\dots ,{\theta }_{d,i,G}\right] i=1,2,3,\dots ,N$$(10) . (10) $${\theta }_{i,G}=\left[{\theta }_{1,i,G},{\theta }_{2,i,G},{\theta }_{3,i,G},\dots ,{\theta }_{d,i,G}\right] i=1,2,3,\dots ,N$$(10) where $\theta$ is a vector of parameters of support vector regression, $d$ is dimensional parameters of support vector regression, $N$ is the size of population and $G$ is the number of generation.

2.6.2. The support vector regression utilizes the given input lag, which is gained from differential evolution to build prediction function by rearranging historical observations into m columns of the historical observations as following. $$\left[\begin{array}{ccccc}{y}_1& y_2& y_3& \cdots & y_m\\ y_2& y_3& y_4& \cdots & y_{m+1}\\ y_3& y_4& y_5& \cdots & y_{m+2}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ y_{t-m+1}& \cdots & y_{t-2}& y_{t-1}& y_t\end{array}\right]$$

The first m − 1 columns of the historical observations are used as input lag, while the last column of the historical observations is uses as target data.

The prediction function of support vector regression is formulated by using a set of hyperplane in the given vector of parameters as Equation (11)(11) $${\widehat{y}}_t=f\left(y_{t-1},y_{t-2},y_{t-3},\dots ,y_{t-m+1}| \theta \right)$$(11) (11) $${\widehat{y}}_t=f\left(y_{t-1},y_{t-2},y_{t-3},\dots ,y_{t-m+1}| \theta \right)$$(11) $f$ is prediction function of support vector regression, and ${\widehat{y}}_t$ is predicted value at time t.

2.6.3. Forecast one step ahead to be future amount of CO₂ emission or CO₂ emission equivalence as described in Equation (12)(12) $$y_{t+1}={\widehat{y}}_{t+1}+{\epsilon }_{t+1}$$(12) . (12) $$y_{t+1}={\widehat{y}}_{t+1}+{\epsilon }_{t+1}$$(12) where $y_{t+1},$ ${\widehat{y}}_{t+1}$ and ${\epsilon }_{t+1}$ are the actual observation, forecast value and error at time t + 1, respectively.

2.6.4. Compare forecast value to the actual value in test data in order to compute MAPE.

2.6.5. An initial mutant vector of parameters $v_{i,G}$ is produced by choosing randomly three members of the population, ${\theta }_{r_0,G},$ ${\theta }_{r_1,G},$ and ${\theta }_{r_3,G}$ as Equation (13)(13) $$v_{i,G}={\theta }_{i,G}+\left({\theta }_{\mathit{\text{best}},G}-{\theta }_{i,G}\right)+{\theta }_{r_0,G}+0.8\times \left({\theta }_{r_1,G}-{\theta }_{r_2,G}\right)$$(13) . (13) $$v_{i,G}={\theta }_{i,G}+\left({\theta }_{\mathit{\text{best}},G}-{\theta }_{i,G}\right)+{\theta }_{r_0,G}+0.8\times \left({\theta }_{r_1,G}-{\theta }_{r_2,G}\right)$$(13) where ${\theta }_{i,G}$ and ${\theta }_{\mathit{\text{best}},G}$ are the ith member of the population at current generation and the best member with the best fitness, respectively. $r_0,$ $r_1$ and $r_2$ are randomly selected numbers within the population size. G is the number of generations; and $i=1,2,3,\dots ,N.$

2.6.6. The crossover operation produces a trial vector $u_{i,G}$ as Equation (14)(14) $$u_{i,G}=\{\begin{array}{l}{v}_{j,i,G} \mathit{\text{if ran}}{d}_{ j,i} \le 0.5\text{or}j= I_{\mathit{\text{rand}}} \\ {\theta }_{j,i,G} \mathit{\text{otherwise}}\end{array}$$(14) . (14) $$u_{i,G}=\{\begin{array}{l}{v}_{j,i,G} \mathit{\text{if ran}}{d}_{ j,i} \le 0.5\text{or}j= I_{\mathit{\text{rand}}} \\ {\theta }_{j,i,G} \mathit{\text{otherwise}}\end{array}$$(14) where G is the number of generations; $i=1,2,3,\dots ,N; j=1,2,3,\dots ,d,$ $\mathit{\text{ran}}{d}_{j,i}\sim U\left[0,1\right],$ $I_{\mathit{\text{rand}}}$ is a random integer from $\left[1,2,\dots ,d\right].$

2.6.7. The differential evolution employs a greedy selection operator as Equation (15)(15) $${\theta }_{i,G+1}=\{\begin{array}{l}{u}_{i,G} \mathit{\text{if}}f\left(u_{i,G}\right)\le f\left({\theta }_{i,G}\right) \\ {\theta }_{i,G} \mathit{\text{otherwise}}\end{array}$$(15) . (15) $${\theta }_{i,G+1}=\{\begin{array}{l}{u}_{i,G} \mathit{\text{if}}f\left(u_{i,G}\right)\le f\left({\theta }_{i,G}\right) \\ {\theta }_{i,G} \mathit{\text{otherwise}}\end{array}$$(15) where $f({\theta }_{i,G})$ is equal to MAPE of the target vector. $f(u_{i,G})$ is the MAPE of the trial vector. G is the number of generations and $i=1,2,3,\dots ,N$

2.6.8. Repeat 2.6.2–2.6.7 until the criterion is met.

Finally, the combined approach of support vector regression and differential evolution returns the optimal parameters of support vector regression.

2.7. Combined approach of proposed model

The proposed model is a complex combined approach emphasizing on both noise reduction and optimal parameter selection, which aims at utilizing MODWT to reduce noise from time series based on first decomposition into scaling and wavelet coefficient series. Afterward, the support vector regression formulates complex function to describe patterns concerning scaling coefficients and wavelet coefficients. In the meantime, differential evolution algorithm is utilized to search appropriate parameters of support vector regression. The algorithm of the proposed model is presented as follows:

2.7.1. The differential evolution creates initial parameters of support vector regression, which are input lag of scaling coefficients, input lag of wavelet coefficients, hyperplane of scaling coefficients, and hyperplane of wavelet coefficients as demonstrated in Equation (16)(16) $${\theta }_{i,G}=\left[{\theta }_{1,i,G},{\theta }_{2,i,G},{\theta }_{3,i,G},\dots ,{\theta }_{d,i,G}\right]i=1,2,3,\dots ,N$$(16) . (16) $${\theta }_{i,G}=\left[{\theta }_{1,i,G},{\theta }_{2,i,G},{\theta }_{3,i,G},\dots ,{\theta }_{d,i,G}\right]i=1,2,3,\dots ,N$$(16) where $\theta$ is a vector of parameters of support vector regression, $d$ is dimensional parameters of support vector regression, $N$ is the size of population and $G$ is the number of generation.

2.7.2. The Original time series is decomposed into scaling and wavelet coefficient series by using Equations (17)–(19). (17) $${\tilde{V}}_{j,t}=\sum _{l=0}^{L-1}{\tilde{g}}_l{\tilde{V}}_{j-1,t-2^{j-1}l\mathrm{\text{mod}}N}, t=0,1,2,3,\dots ,N-1$$(17) (18) $${\tilde{W}}_{j,t}=\sum _{l=0}^{L-1}{\tilde{h}}_l{\tilde{V}}_{j-1,t-2^{j-1}l\mathrm{\text{mod}}N}, t=0,1,2,3,\dots ,N-1$$(18) (19) $$y_t={\tilde{V}}_t+{\tilde{W}}_t$$(19) where ${\tilde{V}}_{j,t}$ and ${\tilde{W}}_{j,t}$ are the $j^{\mathit{\text{th}}}$ level scaling coefficients and the $j^{\mathit{\text{th}}}$ level wavelet coefficients, respectively. The ${\tilde{h}}_l$ and ${\tilde{g}}_l$ are wavelet and scaling filters of MODWT. Letting ${\tilde{V}}_{0,t}=X_t$ and variable L is the length of either wavelet filter $({\tilde{h}}_l)$ or scaling filter $({\tilde{g}}_l).$ $y_t$ is actual observation at time t.

2.7.3. The support vector regression uses the given input lag of scaling coefficients to rearrange the series into m columns. $$\left[\begin{array}{ccccc}{\tilde{V}}_1& {\tilde{V}}_2& {\tilde{V}}_3& \cdots & {\tilde{V}}_m\\ {\tilde{V}}_2& {\tilde{V}}_3& {\tilde{V}}_4& \cdots & {\tilde{V}}_{m+1}\\ {\tilde{V}}_3& {\tilde{V}}_4& {\tilde{V}}_5& \cdots & {\tilde{V}}_{m+2}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\tilde{V}}_{t-m+1}& \cdots & {\tilde{V}}_{t-2}& {\tilde{V}}_{t-1}& {\tilde{V}}_t\end{array}\right]$$

2.7.4. The support vector regression adopts the given input lag of wavelet coefficients in order to rearrange into n columns. $$\left[\begin{array}{ccccc}{\tilde{W}}_1& {\tilde{W}}_2& {\tilde{W}}_3& \cdots & {\tilde{W}}_n\\ {\tilde{W}}_2& {\tilde{W}}_3& {\tilde{W}}_4& \cdots & {\tilde{W}}_{n+1}\\ {\tilde{W}}_3& {\tilde{W}}_4& {\tilde{W}}_5& \cdots & {\tilde{W}}_{n+2}\\ ⋮& ⋮& ⋮& ⋮& ⋮\\ {\tilde{W}}_{t-n+1}& \cdots & {\tilde{W}}_{t-2}& {\tilde{W}}_{t-1}& {\tilde{W}}_t\end{array}\right]$$

2.7.5. The support vector regression exploits the given parameters of hyperplane to model complex function for scaling coefficient series as presented in Equation (20)(20) $${\widehat{\tilde{V}}}_t=f\left({\tilde{V}}_{t-1},{\tilde{V}}_{t-2},{\tilde{V}}_{t-3},\dots ,{\tilde{V}}_{t-m+1}|\theta \right)$$(20) . (20) $${\widehat{\tilde{V}}}_t=f\left({\tilde{V}}_{t-1},{\tilde{V}}_{t-2},{\tilde{V}}_{t-3},\dots ,{\tilde{V}}_{t-m+1}|\theta \right)$$(20) where $f$ is prediction function of support vector regression for scaling coefficient series.

2.7.6. The given parameters of hyperplane are used to formulate complex function of support vector regression, which is described as Equation (21)(21) $${\widehat{\tilde{W}}}_t=f\left({\tilde{W}}_{t-1},{\tilde{W}}_{t-2},{\tilde{W}}_{t-3},\dots ,{\tilde{W}}_{t-n+1}|\theta \right)$$(21) . (21) $${\widehat{\tilde{W}}}_t=f\left({\tilde{W}}_{t-1},{\tilde{W}}_{t-2},{\tilde{W}}_{t-3},\dots ,{\tilde{W}}_{t-n+1}|\theta \right)$$(21) where $f$ is prediction function of support vector regression for wavelet coefficient series.

2.7.7. Aggregate predicted scaling coefficient and predicted wavelet coefficient to be a future amount of CO₂ emission or CO₂ emission equivalence as presented in Equation (22)(22) $$y_{t+1}={\widehat{\tilde{V}}}_{t+1}+{\widehat{\tilde{W}}}_{t+1}+{\epsilon }_{t+1}$$(22) . (22) $$y_{t+1}={\widehat{\tilde{V}}}_{t+1}+{\widehat{\tilde{W}}}_{t+1}+{\epsilon }_{t+1}$$(22) where $y_{t+1},$ ${\widehat{\tilde{V}}}_{t+1},$ ${\widehat{\tilde{W}}}_{t+1}$ and ${\epsilon }_{t+1}$ are the actual observation, predicted scaling coefficient, predicted wavelet coefficient and error at time t + 1, respectively.

2.7.8. An initial mutant vector of parameters $v_{i,G}$ is generated by selecting randomly three members of the population, ${\theta }_{r_0,G},$ ${\theta }_{r_1,G},$ and ${\theta }_{r_3,G}$ as Equation (23)(23) $$v_{i,G}={\theta }_{i,G}+\left({\theta }_{\mathit{\text{best}},G}-{\theta }_{i,G}\right)+{\theta }_{r_0,G}+0.8\times \left({\theta }_{r_1,G}-{\theta }_{r_2,G}\right)$$(23) . (23) $$v_{i,G}={\theta }_{i,G}+\left({\theta }_{\mathit{\text{best}},G}-{\theta }_{i,G}\right)+{\theta }_{r_0,G}+0.8\times \left({\theta }_{r_1,G}-{\theta }_{r_2,G}\right)$$(23) where ${\theta }_{i,G}$ and ${\theta }_{\mathit{\text{best}},G}$ are the ith member of the population at current generation and the best member with the best fitness, respectively. $r_0,$ $r_1$ and $r_2$ are randomly selected numbers within the population size. G is the number of generations; and $i=1,2,3,\dots ,N.$

2.7.9. The crossover operation creates a trial vector $u_{i,G}$ as Equation (24)(24) $$u_{i,G}=\{\begin{array}{l}{v}_{j,i,G} \mathit{\text{if ran}}{d}_{j,i}\hspace{0.17em}\le \hspace{0.17em}0.5\text{or}j= I_{\mathit{\text{rand}}}\\ {\theta }_{j,i,G} \mathit{\text{otherwise}}\end{array}$$(24) . (24) $$u_{i,G}=\{\begin{array}{l}{v}_{j,i,G} \mathit{\text{if ran}}{d}_{j,i}\hspace{0.17em}\le \hspace{0.17em}0.5\text{or}j= I_{\mathit{\text{rand}}}\\ {\theta }_{j,i,G} \mathit{\text{otherwise}}\end{array}$$(24) where G is the number of generations; $i=1,2,3,\dots ,N; j=1,2,3,\dots ,d,$ $\mathit{\text{ran}}{d}_{j,i}\sim U\left[0,1\right],$ $I_{\mathit{\text{rand}}}$ is a random integer from $\left[1,2,\dots ,d\right].$

2.7.10. A greedy selection operator is utilized as Equation (25)(25) $${\theta }_{i,G+1}=\{\begin{array}{l}{u}_{i,G} \mathit{\text{if}}f\left(u_{i,G}\right)\le f\left({\theta }_{i,G}\right) \\ {\theta }_{i,G} \mathit{\text{otherwise}}\end{array}$$(25) . (25) $${\theta }_{i,G+1}=\{\begin{array}{l}{u}_{i,G} \mathit{\text{if}}f\left(u_{i,G}\right)\le f\left({\theta }_{i,G}\right) \\ {\theta }_{i,G} \mathit{\text{otherwise}}\end{array}$$(25) where $f({\theta }_{i,G})$ is equal to MAPE of the target vector. $f(u_{i,G})$ is the MAPE of the trial vector. G is the number of generations and $i=1,2,3,\dots ,N$

2.7.11. Repeat 2.7.2–2.7.10 until the criterion is met.

Finally, the proposed model returns optimal parameters of support vector regression dealing with both scaling and wavelet coefficients.

3. Cross validation

In order to evaluate forecasting performance of those forecasting models, datasets of CO₂ emission or CO₂ emission equivalence are divided into two datasets, which are 70% and 30% of total data for training and test datasets, respectively. For modeling, the training data are utilized to formulate prediction function and forecast one-step ahead value to compare to the first observation of test data with the purpose of calculating error. After the first observation of test data is known, then it is combined with the current training data for updating training data for modeling and predicting one-step ahead value to compare next observation of test data. This process is continuously conducted to the last observation of test data. The mathematical expressions of five accuracy measures are described as Equations (26)–(30). (26) $$\text{MAE}=\frac{\sum _{t=1}^n\left|y_t-{\widehat{y}}_t\right|}{n}$$(26) (27) $$\text{RMSE}=\sqrt{\frac{\sum _{t=1}^n{\left(y_t-{\widehat{y}}_t\right)}^2}{n}}$$(27) (28) $$\text{RMSPE}=\sqrt{\frac{\sum _{t=1}^n{\left(y_t-{\widehat{y}}_t/y_t\times 100\right)}^2}{n}}$$(28) (29) $$\text{MAPE}=\frac{\sum _{t=1}^n\left|y_t-{\widehat{y}}_t\right|/y_t}{n}\times 100$$(29) (30) $$\text{sMAPE}=\frac{\sum _{t=1}^{n}2\times \left|y_t-{\widehat{y}}_t\right|/\left(y_t+{\widehat{y}}_t\right)}{n}\times 100$$(30) where $y_t$ and ${\widehat{y}}_t$ are actual observation and forecast value at time t, respectively.

4. Results and discussion

For experiments of this research by using actual data of CO₂ emission or CO₂ emission equivalence, all empirical results of the five accuracy measures are exploited to indicate better performance of forecasting models, which are summarized and presented in Table 2.

Table 2. The summary of five accuracy measures based on CO₂ emission or CO₂ emission equivalence by sector.

Sector	Gas	Model	MAE	MAPE	RMSE	sMAPE	RMSPE
Transportation	CO2	ARIMA	2.132	3.226	2.534	3.248	3.777
		Holt	2.164	3.274	2.561	3.293	3.820
		SES	2.613	3.917	3.195	4.019	4.686
		SVR-DE	2.028	3.071	3.132	3.184	4.690
		Proposed model	1.701	2.549	2.266	2.588	3.405
Energy	CO2	ARIMA	6.866	2.866	7.856	2.837	3.287
		Holt	6.818	2.846	7.829	2.818	3.275
		SES	6.080	2.498	7.776	2.543	3.208
		SVR-DE	4.925	2.025	6.700	2.057	2.751
		Proposed model	2.433	0.995	3.457	0.999	1.413
Electricity	CO2	ARIMA	4.465	4.144	5.207	4.081	4.853
		Holt	4.928	4.559	5.831	4.479	5.409
		SES	3.724	3.426	4.900	3.474	4.497
		SVR-DE	3.384	3.121	4.739	3.166	4.337
		Proposed model	2.732	2.519	4.019	2.484	3.756
Manufactory	CO2	ARIMA	4.098	8.206	4.774	7.931	9.724
		Holt	3.812	7.639	4.515	7.367	9.220
		SES	3.907	7.726	4.496	7.727	8.888
		SVR-DE	2.812	5.501	3.424	5.589	6.562
		Proposed model	1.563	3.185	2.139	3.180	4.377
Bunker fuels	CO2	ARIMA	0.815	5.241	0.906	5.280	5.847
		Holt	0.690	4.384	0.810	4.379	5.146
		SES	0.787	4.928	0.909	5.032	5.531
		SVR-DE	0.471	2.826	0.755	2.916	4.290
		Proposed model	0.299	1.876	0.476	1.874	3.073
Total excluding LUCF	CO2	ARIMA	7.044	2.769	7.908	2.743	3.113
		Holt	6.745	2.652	7.530	2.631	2.968
		SES	6.219	2.408	8.008	2.453	3.121
		SVR-DE	5.226	2.027	6.664	2.049	2.596
		Proposed model	2.589	1.018	3.106	1.023	1.238
Energy	All GHG	ARIMA	7.206	2.857	8.285	2.829	3.293
		Holt	7.213	2.857	8.209	2.833	3.259
		SES	6.358	2.482	8.191	2.527	3.210
		SVR-DE	5.038	1.968	6.509	1.985	2.554
		Proposed model	4.249	1.647	7.044	1.671	2.727
Electricity	All GHG	ARIMA	4.450	4.106	5.228	4.045	4.844
		Holt	4.887	4.495	5.731	4.421	5.281
		SES	3.805	3.478	4.952	3.527	4.519
		SVR-DE	3.375	3.095	4.982	3.163	4.530
		Proposed model	2.753	2.468	4.416	2.532	3.873

The bold values are meant to be the least error.

According to results in Table 2, each conventional model can perform well in each situation compared to one another. Furthermore, the MAPE is lower than 10%, which indicates more accuracy and can support reliable forecast. However, it is difficult to indicate that one forecasting model outperforms one another due to its advantage and disadvantage of each model. Furthermore, the combined approach of support vector regression and differential evolution outperforms almost accuracy measures compared to conventional models. Although the SVR-DE model performs well than conventional models, the proposed model outperforms all forecasting models in almost accuracy measures.

With regard to support significant difference of forecasting performances, Friedman test, which is nonparametric analysis of variance, is conducted to identify that at least one is significantly different from one another. Afterward, the post hoc pairwise multiple comparison of Conover test is conducted to identify significant difference between forecasting performances. The summary of Friedman test is presented in Figures 2–4.

Figure 2. The summary of Friedman test based on MAPE.

Figure 3. The summary of Friedman test based on sMAPE.

Figure 4. The summary of Friedman test based on RMSPE.

Referencing p value of Friedman test based on MAPE is less than 0.05, it is sufficient to conclude that at least model is significantly different from one another at 0.05 significance levels. Moreover, the Conover’s test reveals that the proposed model is significantly different from all forecasting models at 0.05 significance levels. Furthermore, the SVR-DE model is significantly different from all conventional models at 0.05 significance levels, although it cannot outperform the proposed model.

As the results of Friedman test depended on sMAPE, the p value is also less than 0.05 and demonstrates that at least performance is not equal to one another at 0.05 significance levels. In addition, the proposed model is still significantly different from all forecasting models at 0.05 significance levels. In the meantime, the SVR-DE model is still significantly different from all traditional models at 0.05 significance levels.

From results of Friedman test in Figure 4, it indicates that at least performance is not equal to one another at 0.05 significance levels because p value is less than 0.05. Based on Conover’s test, empirical results reveal that proposed model is still significantly better than all forecasting models at 0.05 significance levels. Conversely, the SVR-DE can significantly overcome only all traditional models at 0.05 significance levels.

5. Conclusion

According to all results of evaluations, the proposed model emphasizing on both noise reduction and optimal parameter selection techniques is significantly superior to all forecasting models. This evidence can support to conclude that noise reduction technique can remove the irrelevant and redundant components of time series, which can provide more explicit and stationary time series than original time series. Moreover, the differential evolution algorithm can search the most proper parameters of support vector regression within given search space. Meanwhile, the SVR-DE model significantly outperforms all conventional models. In other words, these combined approaches can improve more accuracy and provide more reliable forecast to support decision making on those problems of CO₂ emission or CO₂ emission equivalence.

Consequently, the proposed model can be a useful tool for extrapolating future amount of CO₂ emission or CO₂ emission equivalence to support policy planner to take step with helpful guideline. However, the proposed model can provide one year ahead forecast, which may not be suitable for middle or long-term planning.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This research was supported by Research and Graduate Studies, Khon Kaen University.