FORECASTING EUROPEAN UNION CO2 EMISSIONS USING AUTOREGRESSIVE INTEGRATED MOVING AVERAGE-AUTOREGRESSIVE CONDITIONAL HETEROSCEDASTICITY MODELS

In the past few decades, there are lot of discussions around global warming and climate change primarily due to the increased CO2 emissions generated by the consumption of fossil fuels, such as oil and natural gas. After an enormous effort, the EU-28 managed to reduce CO2 emissions in 2014 by 25.7% comparing to 1990 (Kyoto Protocol). This effort should continue in the future so that the EU-28 achieve a 40% reduction on CO2 emissions by 2030. The current paper aims at investigating the optimum model to forecast CO2 emissions in the EU-28. To achieve this aim an ARIMA(1,1,1)-ARCH(1) model was used, combined with the linear ARIMA model and the conditional variance of the ARCH model. The estimation of parameter optimisation of ARIMA(1,1,1)-ARCH(1) model was done with the Maximum Likelihood approach using the Marquardt (1963), and Berndt-Hall-Hall-Hausman (BHHH) algorithms and the three distributions (Normal, t-Student, Generalized error), whereas for the estimation of the covariance coefficient the reversed matrix by Hessian was used. Finally, in order to forecast the ARIMA(1,1,1)-ARCH(1) model, a dynamic as well as a static process was applied. The results of the forecasting revealed that the static procedure provides a better forecast comparing to the dynamic one.


INTRODUCTION
During the last decades, the increase of the greenhouse gas emissions is considered the biggest threat for global warming. The economic growth of developed countries pushes the intensive use of energy and the consumption of fossil fuels, which results to more residues and waste leading to environmental decomposition.
Data from the 1960s and 1970s, show that the concentration of CO 2 in the atmosphere is increasing significantly. Hence, scientists put a lot of pressure to governments to take action.
Unfortunately, it took the international community years to respond to this request. Back in 1988 the International Meteorology Organisation and the United Nations Environmental Program formed an intergovernmental committee for climate change.
The evaluation report published in 1990, mentioned that the problem of temperature rise is an existing one and owes to be dealt with promptly. This outcome pushed governments to establish the United Nations Framework Convention on Climate Change (UNFCCC), which was signed in Rio de Janeiro in 1992. UNFCCC as well as the Kyoto Protocol that followed, form the only international frameworks for tacking climate change The Kyoto Protocol is an international treaty (signed on 11 December 1997, but entered into force on 16 February 2005), which sets the principles of reducing greenhouse gas concentrations in the atmosphere which cause temperature rise in the planet. Based on Kyoto Protocol, countries which have signed the treaty, are bound to reduce greenhouse gas emissions on average by 5.2% between 2008 and 2012, comparing to 1990 levels. The protocol is based on the principle of common responsibilities in tackling climate change but acknowledges the different capacity to achieve this based on each country's economic growth.
Negotiations for resolving the global temperature increase problem were tough due to clashing interests. Consequently, opposing country-teams with diverging views were generated. For example, countries producing carbon, such as Japan, USA, Canada, Australia, New Zealand, but also members of OPEK with Russia and Norway, which support the development of oil and natural gas production, are affected by Kyoto Protocol because they need to reduce their production and instead switch to alternative energy sources. In addition, emerging makers such as China and India could not commit to reducing greenhouse gas. On the contrary, the European Union was the earnest supporter of Kyoto Protocol, which committed to reduce greenhouse gas emissions by 6% for 2012 and reduction of fossil fuel consumption by 16%.
Kyoto Protocol followed a second phase during the period 2013-2020 known as Doha Amendment 2012. Thirty-eight developed countries participate among them 28 member states of the European Union. In this second phase, participating states commit to reduce their emissions to a level 18% lower than that of 1990. The EU also has committed to reduce by 20% for that period.
The Paris agreement for climate change took place in December 2015. This conference set two requirements for its application. The first one concerned the confirmation of the UNFCCC, by at least 55 country-members and the second one concerned the minimum amount of greenhouse gas emissions, which each country need to follow.
The Kyoto Protocol applies to six greenhouse gases: Carbon dioxide (CO 2 ), methane (CH 4 ), nitrous oxide (N 2 O), hydrofluorocarbons (HFCs), perfluorocarbons (PFCs), and sulphur hexafluoride (SF 6 ). Carbon dioxide (CO 2 ) is a natural gas, which is defined by the photosynthesis into organic material. Most of the CO 2 emissions are due to fossil fuel consumption such as carbon, petroleum, natural gas and the burning of biomass. CO 2 emissions are also generated by the change in soil use, the car mobility and various others industrial activities. CO 2 is the main anthropogenic greenhouse gas that affects the radioactive balance of the earth. Moreover, CO 2 is the gas that form the basis in which other greenhouse gas are calculated, resulting in a dynamic overheating of the planet.
The majority of the countries, failed to achieve their commitments with regards to CO 2 emissions. The European Union is taking notable measures for reducing CO 2 emissions generated from the absorption of fossil fuels.
Given that there exists a two-way causal relationship between economic growth and CO 2 emissions, a reduction on the CO 2 emissions will have an unfavourable impact on the economic growth of European Countries Table A1. Table A2, presents the relationship between per capita CO 2 emissions and the per capita economic power of the 28 European Countries.
The amount of CO 2 emissions during the period 1990-2014 in EU, USA, China and the World is presented in the Figure 1.
From the Figure 1 we detect that since the Kyoto Protocol in 1990, the USA have reduced CO 2 emissions by 0.63 % pa and the EU by 1.19% pa. On the contrary, China has increased CO 2 emissions by 5.49% pa whereas the global CO 2 emissions are increasing by 0.73% pa.
The Table 1 presents the descriptive statistics for the per capita CO 2 emissions in the European Union, the USA, China and the World from 1990 (Kyoto Protocol) until 2014. The descriptive statistics mean, standard deviation (Std. Dev.) and coefficient of variation (CV) of these variables are recorded below in Table 1. Table 1 shows that the variability in the per capita CO 2 emissions is greater in China and smaller in in the case of the USA.
The Figure 2 presents the rate of CO 2 emissions in the EU-28 countries.
The Figure 2 shows an overall downward trend of the CO 2 emissions between 1990 and 1999, with the exception of a peak in 1996, when an exceptional of a cold winter which led to increased demand for heating. From 1999 to 2006, the CO 2 emission for the EU-28 was relative stable. The remainder of the paper is organized as follows: Section 2 provides a brief literature review. Section 3 presents the analysis of methodology. Section 4 summarizes the data. The empirical results are provided in Section 5 and Section 6 proposes the forecasting results. Finally, the last section offers the concluding remarks.

LITERATURE REVIEW
The forecasting issues are important and are being applied in various scientific fields, such as economics, meteorology, medicine, mechanics, ecology and many more. The increasing trend of the CO 2 emissions in a global level due to human activity indicated the increased atmospheric concentration of CO 2 . Climate change, because of global warming, is one of the most prolific issues during the last years. Reddy et al. (1995) suggest in their research that the total average temperature increase will reach 3-4°C, doubling the CO 2 emissions concentration, whereas in 2007 the intergovernmental committee for climate change reported an increase of the temperature between 1.1 to 6.4°C in the next 100 years.
Developed countries have a higher share of global CO 2 emissions comparing to developing countries. Nakicenovic back in 1994, studies the prospect of greenhouse gas emissions in a rural context. His findings suggest that developing countries are responsible for less than the 16% of the CO 2 concentration due to their previous consumption from mineral sources of energy. Developed countries have a higher share of global emissions. Researchers so far have investigated the forecasting of CO 2 emissions in various countries. Lotfalipour et al. (2013) have examined the economic aspects of CO 2 emissions and their consequences for the case of Iran. In their study they apply Grey and autoregressive integrated moving average (ARIMA) models to forecast CO 2 emission in the period between 1990 and 2011. Their results suggest that Grey models provided better results in terms of forecasting CO 2 emissions. Based on their estimated results, the quantity of CO 2 emissions will reach 925.680.000 tons in 2020, which indicated an increase of 66% compared to 2010, which is highly significant. (2017), investigated the CO 2 emissions between 1972 and 2015 in the case of Bagladesh. The optimum prognostic model for the CO 2 emissions in the period under investigation was the ARIMA(0, 2, 1) model. The results suggested that the CO 2 emissions for years 2016, 2017 and 2018 will be 83.94657, 89.90464 and 96.28557 metric tons respectively.

Rahman and Hasan
Pruethsan in 2007, analysed CO 2 emissions in Thailand using the VARIMAX approach during the period 2000-2015 and subsequently determined the VARIMAX(2, 1, 1) and VARIMAX(2, 1, 3) models as the optimal ones for the CO 2 forecast emissions in Thailand. The forecasting results (using the VARIMAX(2, 1, 1) model) show that CO 2 gas greenhouse   (2012) emissions will increase steadily and will reach 25.17% until 2025 comparing to 2016, whereas using the VARIMAX(2, 1, 3) model the CO 2 gas greenhouse emissions will increase steadily and will reach 41.51% until 2040 comparing to 2016.
Nyoni and Mutongi (2019) are using annual data to investigate CO 2 emissions in the case of China from 1960 to 2017 using the Box-Jenkins approach. The ARIMA(1, 2, 1) model was shown to be the most suitable one to forecast CO 2 emissions in the period under investigation. The study results reveal that CO 2 emissions in China will increase and reach approximately 10.000.000 kt in 2024. This result forms a warning of the Chinese government with regards to clinical change and the overheating that China causes in the world.
Finally, Nyoni and Bonga (2019) use 1960-2017 data and a Box-Jenkins approach to forecast CO 2 emissions in China. The study proposes the ARIMA(2, 2, 0) as the optimum one to forecast CO 2 emissions. It further suggests that by 2025, the annual CO 2 emissions in India will reach 3.890.000 kt. The results is critical to the Indian government with respect its short-term and long-term planning of climate change and global overheating.

ARIMA Models
An ARIMA model is a generalization of an autoregressive moving average (ARMA) used in econometrics. ARIMA is one of the type of models in the Box and Jenkins (1976) methodology for analysis and forecasting a time series.
The ARIMA(p,d,q,) can be expressed as:  Engle (1982) developed the ARCH model for testing the volatility of time series. The basic ARCH model consists of two equations, a conditional mean equation and a conditional variance equation.

ARCH(q) model
Both equations form a system that is estimated together with maximum likelihood (ML) method. So, ARCH model is an ARMA and can be written as follows: where u t is conditional mean of y t , and  t is the shock at time t.
The variance  t will be: where u t is a white noise with zero mean and variance of one where  t 2 is the conditional variance of y t , ω is a constant term, and q is the order of the ARCH terms,  > 0 ,  i ≥ 0 and i > 0 . Bollerslev (1986) extended the ARCH model in a new one that allows the errors of variance to depend on its own lags as well as lags of the squared errors. In other words, it allows the extension of conditional variance to follow an ARMA process.

GARCH(q,p) model
The GARCH model can be expressed as: where u t is conditional mean of y t , and  t is the shock at time t.
where u t is a white noise with zero mean and variance of one where  t 2 is the conditional variance of y t , ω is a constant term, q is the order of the ARCH terms, and p is the order of the GARCH terms. We assume that for every p ≥ 0 and q > 0, the parameters are unknown and since the variance is positive, then the following relations must be positive too  ≥ 0 , and  i ≥ 0 for every i=1,…, q and  j ≥ 0 for j = 1,…, p.

ARIMA-ARCH/GARCH Model
The ARIMA-ARCH/GARCH model is one model in which the variance of the error term of the ARIMA model follows an ARCH/ GARCH process. In other words, the ARIMA-ARCH/GARCH model is a non-linear time series model which combined the lineal model ARIMA with the conditional variance of the ARCH/ GARCH model.
For the ARIMA-ARCH/GARCH process to be suggested, the following two phases should be applied. The first one uses the best ARIMA model which fits the stationary and linear data of the time series, whereas the linear model residuals should contain the non-linear part of the data. The second phase uses the ARCH/ GARCH model in order to include the non-linear patters of the residuals. The model combined the ARIMA model with the ARCH/ GARCH which contains the non-linear patterns of the residuals (Dritsaki, 2018).
The process of parameter estimation of the ARIMA-ARCH/ GARCH model is achieved through the logarithmic function of ML through nonlinear least squares using Marquardt's algorithm (1963). The latter is presentd by the following function: where θ is the vector of the parameters that have to be estimated for the conditional mean, conditional variance and density function, z t denoting their density function, D z t ( ( ), ) θ υ , is the log-likelihood function of [ ( )] y t  , for a sample of T observation.

Conditional Distributions
Logarithmic function of ML used for parameters' estimation on volatility models for all theoretical distributions are the following: (Dritsaki, 2017;2019).
where θ is the vector of the parameters that have to be estimated for the conditional mean, conditional variance and density function, T is observations.

Diagnostic Checking of the Model ARIMA-ARCH/GARCH
Before we accept a fitted model and interpret its findings, it is essential to check whether the model is correctly specified, that is, whether the model assumptions are supported by the data. The diagnostic tests of ARIMA-ARCH/GARCH models are based on residuals. Residuals' normality test is employed with Jarque and Bera (1980) test. Ljung and Box (1978) (Q-statistics) statistic for all time lags of autocorrelation is used for the serial correlation test. Also, for the conditional heteroscedasticity test we use the squared residuals of autocorrelation function.

Forecasting Performance Measures
In order to compare the forecasting performance of ARIMA-ARCH/GARCH models we use the following statistics: • Mean absolute error (MAE) • MSE has the following formula MSE gives an overall idea of the error occurred during forecasting.
• The mean absolute percentage error (MAPE) This measure represents the percentage of average absolute error occurred. It is independent of the scale of measurement, but affected by data transformation.

DATA
Annual data for CO 2 emissions (CO 2 ) in the E.U (metric tons per capita), are downloaded from the World Bank's development indicators. The data is for the period from 1960 to 2014. Figure 3 presents the course of the per capita CO 2 emissions in the European Union between 1960 and 2014 at the level. Figure 3 shows that the CO 2 emissions at the EU present a random walk. Therefore, we will check for the stationarity of the series and its Figure 4 autocorrelation.
From Figure 4 shows that the auto-regression coefficients decline with rapid pace, which implied that the series is non-stationary.
We then apply the aforementioned tests afresh, in order to investigate the presence of stationarity in the first differences of the series. Figure 5, shows the first differences of the CO 2 emissions.
From Figure 5, we observe that the CO 2 emissions present intense fluctuations in their first differences, which is a possible indication of stationarity. We then test the stationarity with the auto-correlation Figure 6.
Figures 5 and 6 both show that the series is likely to be stationary in its first differences.
The confirmation of the series stationarity is achieved by applying the unit root tests Dickey and Fuller (1979;1981) and Phillips and Perron (1998). The results of Table 2, confirm that the series is stationary in the first differences. The next step is to determine the ARIMA(p, q) model, based on the results from Figure 6. The parameters p και q of the ARIMA model could be determined rom the partial autocorrelation and auto-correlation coefficients comparing them respectively with the critical value ± =± =± 2 2 55 0 262 n . .
Moreover, to test for autocorrelation we use the Ljung and Box (1978) test determined by: Looking at the values of the partial autocorrelation and autocorrelation coefficients ( Figure 6) the value for p is between 0 < p < 1 and respectively, the value for q will be between 0 < q < 2. Using the values above, we chose the best ARIMA(p, d, q) model from the lowest values of AIC, SC, and HQ criteria. Table 3 shows the values for p and q.   C,T C C,T CO 2 EE -1.700(0) -1.787(0) -1.884(4) -1.786(1) D CO 2 EE -5.284(0)* -6.854(0)* -5.502(4)* -6.870(2)* * , ** and *** show significant at 1%, 5% and 10% levels respectively. (2) The numbers within parentheses followed by ADF statistics represent the lag length of the dependent variable used to obtain white noise residuals.
(3) The lag lengths for ADF equation were selected using Schwarz information criterion. (4) Mackinnon (1996) critical value for rejection of hypothesis of unit root applied. (5) The numbers within brackets followed by PP statistics represent the bandwidth selected based on Newey and West (1994) method using Bartlett Kernel. (6) C=Constant, T=Trend. (7) Δ=First differences Results from Table 3, reveal that based on Akaike (AIC), Schwartz (SIC) and Hannan-Quinn (HQ) criteria, ARΙMA(1,1,1) model is the most suitable one.
The estimation of the ARIMA(1,1,1) model is achieved by the ML method, whereas the optimisation of the model will be achieved using the Berndt-Hall-Hall-Hausman (BHHH), algorithm. The covariance coefficient will be estimated with the inverse Hessian matrix. Table 4, shows the results of estimating the ARIMA(1,1,1) model. Table 4, show that there aren't any problems with the significance of the coefficients. Also, the coefficient for the error variation estimation, shown as SIGMASQ, is ρ s = 0.058 and is also statistical significant. The inverse roots of the model are AR = 0.94 and MA = 0.76 and are presented in the following Figure 7. Figure 7 shows that the inverse roots of the model (inverted AR, MA roots), are within the inverted unit cycle, which confirms that the series under review is stationary.

Results from
We then test for heteroscedasticity (ARCH(q)) from the squared residuals of the model above. Table 5 shows the following results. Table 5 reveals that both auto-correlation coefficients and partial autocorrelation coefficients after the first order, are not statistical significant. Therefore, the ARCH or GARCH procedures should be considered.

EMPIRICAL RESULTS
Since we found that there exists a first order ARCH-GARCH procedure (Table 5), we could move into the model specification followed by its estimation of the conditional mean and the  conditional variance. Estimation of the ARIMA(1,1,1)-ARCH(1) or ARIMA(1,1,1)-GARCH(1,1) model is achieved with the ML method using the BHHH algorithm, the steps of the Marquardt method (1963), the three distributions (Normal, student's, generalized error), while for the co-variance coefficient the inverse matix by Hessian is applied. Coefficient estimation as well as residual tests with respect to normality, auto-regression and conditional heteroscedasticity, are presented in Table 6.
From Table 6 shows that the ARIMA(1,1,1)-ARCH(1) model with the GED distribution is the most appropriate one (as it has the highest LogL value). All model coefficients are statistical significant and do not present any issues at the diagnostic issues.
Therefore, we could use this model for forecasting purposes. Table 7 presents the results of the estimation of this model.
The estimate of the regression in the ARIMA(1,1,1)-ARCH(1) model could be presented as: D CO 2 EE t = 0.933*D CO 2 EE t-1 -0.846e t-1 (conditional mean equation)  t 2 =0.018+0.956 e 2 t-1 (conditional variance equation) Figure 8 shows the actual and fitted values of the series, as well as the residuals of the fitted model at the 95% confidence interval.   (24) is the Q-statistic of correlogram of squared residuals at twenty-four lags. (5) ARCH-X 2 (1) for autoregressive conditional heteroskedasticity, (6) the persistence is calculated as (α 1 +β 1 ) for the GARCH model Looking at Figure 8, the residuals show that there is an ARCH procedure in the data.

FORECASTING
In order to forecast the ARIMA(1,1,1)-ARCH(1,1) models we use both the dynamic (n-step ahead forecasts) and static (one step-ahead forecast) procedure. The dynamic procedure computes forecasting for periods after the first sample period, using the former fitted values from the lags of dependent variable and ARMA terms. The static procedure uses actual values of the dependent variable. In the following diagram, we present the criteria for the evaluation of forecasting using the dynamic and static forecast respectively (Dritsaki, 2019).  The Figure 9 indicate that the static procedure gives better results rather than the dynamic (MSE and MAE are lower in the static rather than the dynamic process). Since ARIMA(1,1,1)-ARCH(1,1) model is fit to the CO 2 data, therefore we can use to forecast values for the next 6 years out-of sample (from 2015 to 2020). The forecasted values of CO 2 are given in Table 8. The forecasted values indicate that the CO 2 present fluctuations until the end of 2020. The great drop as shown by the results of the current study for 2020, is consistent with the commitment which the EU promised by the Kyoto protocol, as well as the amendment of Doha in 2012.

SUMMARY AND CONCLUSION
The degradation of environment becomes the recrudescence of the environment via the exhaustion of sources such as air, water and soil. This degradation is a consequence of a combination of an already big and constantly growing population, the continuing economic growth, the technological exhaustion of natural resources and the pollutant technology. The results of these damages to the human lifestyle and prosperity have caused a great amount of concern.
Many studies showed that there is an influence of economic growth to environmental degradation. The correlation between  per capita GDP and CO 2 emissions is positive, implying that the increasing per capita GDP leads to increase of CO 2 emissions. No turning point is found at which emissions start to decrease. Market economy mechanisms according to studies' results are not sufficient in the decline of CO 2 emissions. For this reason, legal regulations are required to avoid further environmental degradation.
The purpose of this paper is to model and forecast CO 2 emissions of 28 member countries of EU based on annual data (from 1960 until 2014). Using ARIMA(1,1,1)-ARCH(1) model, ML methodology, Marquardt's algorithm methodology (1963) and BHHH, we forecasted CO 2 emissions for the next 6 years (2015-2020). The results of forecasting showed that CO 2 emissions will display fluctuations until the end of 2020. The year 2020 will present a considerable decrease of CO 2 emissions reaching 33.8% less than the year 1990 (Kyoto Protocol) and will cover by far the commitment of EU countries on the above Treaty. Moreover, after the biblical disasters worldwide (the burning of Amazon forest which covers 60% of the total rainforest, the lack of drinkable water) more countries such as USA, China, India and OPEC countries can adapt with last years' phenomena and decrease CO 2 emissions, avoiding planet's major disaster. According to our examined model, European Union will manage and reduce more CO 2 emissions by 40% in relation to 1990 and once more will be consistent with Kyoto protocol.