# Pakistan CO2 Emission Modelling and Forecasting: A Linear and Nonlinear Time Series Approach

Authors: Kassim Tawiah; Muhammad Daniyal; Moiz Qureshi
Journal: Journal of Environmental and Public Health
 (2023)

---

## Abstract

Pakistan is considered among the top five countries with the highest CO2 emissions globally. This calls for pragmatic policy implementation by all stakeholders to bring finality to this alarming situation since it contributes greatly to global warming, thereby leading to climate change. This study is an attempt to make a comparative analysis of the linear time series models with nonlinear time series models to study CO2 emission data in Pakistan. These linear and nonlinear time series models were used to model and forecast future values of CO2 emissions for a short period. To assess and select the best model among these linear and nonlinear time series models, we used the root mean square error (RMSE) and the mean absolute error (MAE) as performance indicators. The outputs showed that the nonlinear machine learning models are the best among all other models, having the lowest RMSE and MAE values. Based on the forecasted value of the nonlinear machine learning neural network autoregressive model, Pakistan’s CO2 emissions will be 1.048 metric tons per capita by 2028. The increasing trend in emissions is a frightening and clear warning, suggesting that innovative policies must be initiated to reduce the trend. We encourage the Pakistan government to price CO2 emissions by companies and entities per ton, adapt electricity production from hydro, wind, and different sources with no emissions of CO2, initiate rigorous planting of more trees in the populated areas of Pakistan as forest covers, provide incentives to companies, organisations, institutions, and households to come out with clean technologies or use technologies with no CO2 emissions or those with lower ones, and fund more studies to develop clean and innovative technologies with less or no CO2 emissions.


---

## Body

## 1. Introduction

The collection of greenhouse gas (GHG) emissions in airspace is creating chaos across the globe and destroying life, economy, health, and food. Carbon (IV) oxide (CO2) is a major component of GHG [1]. The world is nowhere near securing a global temperature rise of less than two degrees Celsius as proposed in the Paris Agreement of the United Nations (UN) Framework Convention on Climate Change (UNFCCC). Using the 1990 criterion, most countries are still releasing more GHG into the atmosphere, with some at a standstill and others releasing less [1]. The UN Environmental Program (UNEP) Climate Action Note indicated that the world is in a state of climate emergency [1].Pakistan is labelled among the countries with the highest GHG emissions in the world by the UNEP [1]. It contributes to 1.02% of global GHG emissions. Pakistan emitted an estimated 504.59 million tons of GHG in 2018. It is one of the most impacted nations from unfriendly effects of environmental change as well as air contamination. GHG emissions are the leading cause of global warming, thereby becoming the maximum essential research issues within the fields of technological know-how and politics [1].In Southern Asia, Pakistan is ranked among the developing economies with swift economic positive growth. This is expected to continue for a long time. Agriculture is the principal dominant area in Pakistan’s economic growth irrespective of the growing industrial sector. However, this industrial drive, coupled with high population growth, has caused massive deforestation, resulting in less absorption of CO2 by plants, leading to a higher CO2 concentration in the atmosphere [2, 3].Pakistan is ranked first among countries in Asia facing mass deforestation as a result of bulk utilization of natural energy from the environment. This is causing serious degradation of the environment. This degradation of the environment, together with deforestation, has brought about global warming, leading to climate change. The critical factor of Pakistan’s commercial revolution is reworking its economy from organic economies based totally on human and animal energy to inorganic economies, which are completely dependent on fossil fuels. Interestingly, the by-product of fossil fuels is CO2, thereby contributing to the concentration of CO2 in the environment [2, 3].The quest to be an industrialized giant definitely comes at a price. Pakistan is geared towards a full-industrial revolution which encourages the use of fossil fuels to power the growing industry in the country. The growth in population means expanding cities by destroying green vegetation, which absorbs most of the CO2 in the atmosphere. The growing population also means increasing cars with carbon emissions being used in the country [2, 3].The Government of Pakistan has put in place a lot of programmes to reduce GHG emissions in the country. Great progress has been chalked up in climate action that ranges from policies, programmes, and nature-based solutions (NBS) to technology-based innovations [4]. Pakistan, recognizing the role of nature in climate modification and diminution, has put forward robust natural capital restoration plans including the Ten Billion Tree Tsunami Program as well as the Protected Drive. These programmes have contributed to enhancing livelihoods and giving opportunities to the vulnerable, including youth and women. In addition, Pakistan has introduced a number of policy actions focused on mitigating GHG emissions from high-emission sectors, such as energy and industry [4].Notwithstanding, this available information indicates a rise in CO2 emissions in Pakistan. This needs to be scientifically substantiated and verified to assist the government to restrategize in an attempt to achieve the targeted goal in the fight of reducing these emissions in order to achieve the Paris Agreement of the UNFCCC. This study provides a statistical and machine learning framework to delve deep into the discourse of CO2 emissions in Pakistan from 1960 to 2018. It provides CO2 emission forecasts for future and possible policy initiatives and directions to combat this issue based on scientific results to meet the UNFCCC target.


### 1.1. Literature Review

The world is being confronted with significant difficulties connected with an unnatural weather change and discharge of ozone-depleting substances in large amounts. In 2017, energy businesses represented 46% of all CO2 discharges worldwide [5]. Bokde et al. [5] proposed two-time series disintegration techniques for transient determination of CO2 outflow with a bunch of cutting-edge models for a 48-hour skyline power market. By estimating benchmarks for France, they showed that their new technique has a mean outright rate blunder, that is, 25% lower than that of existing models. Furthermore, they illustrated the use of the conjecture for booking adaptable power utilization in Germany, Norway, Denmark, and Poland. Their results showed that booking an adaptable block of 4 hour of power utilization at a 24-hour stretch can on average diminish subsequent CO2 discharge in all countries [5].The carbon exchanging market has turned into a powerful weapon in easing fossil fuel by-products in China, and carbon cost is at the centre of its activity. Thus, the carbon exchanging market fills in as an imperative part in gauging carbon cost precisely ahead of time [6]. Sun and Li [6] imaginatively investigated a group of driven long short-term memory network (LSTM) models in light of reciprocal outfit experimental mode decay (ROEMD) for carbon cost gauging, applying it to each of the eight carbon exchanging piloting centres in China. ROEMD was first executed for mode change to disintegrate the first convoluted mode into a bunch of straightforward modes. Accordingly, LSTM was utilized to demonstrate the planning between time-slacking factors and every model’s objective quality, thereby constructing numerous LSTM models for comparison. At last, converse ROEMD computation was acquainted with the coordinates’ expected consequences of the multimode on eventual outcomes. Its viable application at the same time embraced every one of the eight carbon piloting centres in China, covering their compared carbon cost information over a significantly extended period. The acquired outcomes showed that the proposed model had adequate precision in carbon cost anticipation in China in contrast to the single LSTM model as well as other ordinary machine learning models. Moreover, as indicated by the extent of its application, the creative model displayed solid dependability and comprehensiveness [6].Wu and Liu [7] proposed a multiresolution combo-predicting model that uses ensemble empirical mode decomposition (EEMD-ADD) to improve the predicting accuracy of the cost of carbon. They achieved it by decomposing carbon price sequencing using EEMD into numerous intrinsic framework functions (IFFs) by dividing these IFFs into several high-frequency compartments, low-frequency compartments, and drift compartments after which the least square vector support machine (LSSVM), particle swarm optimization LSSVM (PSO-LSSVM), and bat algorithm-LSSVM (BA-LSSVM) were utilized to forecast three compartments, respectively. Similarly, Yun et al. [8] built a robust carbon price-predicting model of complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN)-long short-term memory (LSTM) by blending the merits of CEEMDAN in decaying multiresolution time-frequency carbon price waves and the LSTM model by installing monetary waves.A report by Eckstein et al. [9] positioned Pakistan as the fifth generally impacted country because of environmental changes throughout the course of recent years, while the author in [10] positioned Pakistan’s air as the second generally contaminated in 2020. To diminish these effects, the nation needs to take broad variation and alleviation measures while changing the energy area towards dirtying and carbon impartial choices [11].Existing studies imply that greater than two-thirds of GHG emissions in Pakistan come from fossil power-related CO2 emissions. Consequently, forecasting CO2 emissions is crucial for adjusting regulations to mitigate weather changes. The forecasting of CO2 emissions has been mentioned dramatically in [12].Pao and Tsai [13] implemented the grey model to obtain CO2 emissions in Brazil between 1980 and 2007. A combination of models based on the optimization algorithms of quantum harmony and discounted mean square has been established and employed in CO2 emission forecasting for the world’s top five countries with high CO2 emissions [14]. Auffhammer and Carson [15] came up with China’s CO2 emission path using a set of panel data in a provincial-level framework [16–19]. The application of the traditional grey model and the upgraded grey model has been implemented to forecast CO2 emissions [20, 21].From qualitative and quantitative perspectives, numerous discussions have been conducted on the relationship between consumption of energy, economic growth, and CO2 emissions [22]. Al-mulali and Sab [23] investigated the influence of the consumption of basic energy and emissions of CO2 on sixteen developing countries, focusing on the impact on the development of their economies. Wang et al. [24] examined the significant factors contributing to CO2 emissions in China’s road freight transport. Zhang and Nian [25] computed the emissions of CO2 and examined the factors contributing to emissions in China’s transport sector. In China, Lin et al. [26] developed a new bootstrap ARDL bounding test for CO2 emissions and industrial development. Pao and Tsai [27] conducted an examination of the dynamism of the causal relationship, linking the emission of pollutants, consumption of energy, and their output in Britain, Russia, India, China, and other countries over a definite period [28–34].Little work has been conducted on providing mathematical and statistical modelling techniques for modelling and predicting CO2 emissions in Pakistan and where the country stands in the future concerning this important issue. It is against this backdrop that this paper provides a modelling and forecasting technique based on linear and nonlinear time series models. This study predicts the amount of CO2 emissions in Pakistan for almost a decade to assist the Pakistan government and other stakeholders to plan ahead to devise ways of mitigating the effects of CO2 emissions or possibly curtailing them.Section2 of the paper presents the source of data, data characteristics, and the modelling and forecasting techniques applied. Section 3 presents the discussion and the results, while the last section presents the conclusions and recommendations of the study.


## 1.1. Literature Review

The world is being confronted with significant difficulties connected with an unnatural weather change and discharge of ozone-depleting substances in large amounts. In 2017, energy businesses represented 46% of all CO2 discharges worldwide [5]. Bokde et al. [5] proposed two-time series disintegration techniques for transient determination of CO2 outflow with a bunch of cutting-edge models for a 48-hour skyline power market. By estimating benchmarks for France, they showed that their new technique has a mean outright rate blunder, that is, 25% lower than that of existing models. Furthermore, they illustrated the use of the conjecture for booking adaptable power utilization in Germany, Norway, Denmark, and Poland. Their results showed that booking an adaptable block of 4 hour of power utilization at a 24-hour stretch can on average diminish subsequent CO2 discharge in all countries [5].The carbon exchanging market has turned into a powerful weapon in easing fossil fuel by-products in China, and carbon cost is at the centre of its activity. Thus, the carbon exchanging market fills in as an imperative part in gauging carbon cost precisely ahead of time [6]. Sun and Li [6] imaginatively investigated a group of driven long short-term memory network (LSTM) models in light of reciprocal outfit experimental mode decay (ROEMD) for carbon cost gauging, applying it to each of the eight carbon exchanging piloting centres in China. ROEMD was first executed for mode change to disintegrate the first convoluted mode into a bunch of straightforward modes. Accordingly, LSTM was utilized to demonstrate the planning between time-slacking factors and every model’s objective quality, thereby constructing numerous LSTM models for comparison. At last, converse ROEMD computation was acquainted with the coordinates’ expected consequences of the multimode on eventual outcomes. Its viable application at the same time embraced every one of the eight carbon piloting centres in China, covering their compared carbon cost information over a significantly extended period. The acquired outcomes showed that the proposed model had adequate precision in carbon cost anticipation in China in contrast to the single LSTM model as well as other ordinary machine learning models. Moreover, as indicated by the extent of its application, the creative model displayed solid dependability and comprehensiveness [6].Wu and Liu [7] proposed a multiresolution combo-predicting model that uses ensemble empirical mode decomposition (EEMD-ADD) to improve the predicting accuracy of the cost of carbon. They achieved it by decomposing carbon price sequencing using EEMD into numerous intrinsic framework functions (IFFs) by dividing these IFFs into several high-frequency compartments, low-frequency compartments, and drift compartments after which the least square vector support machine (LSSVM), particle swarm optimization LSSVM (PSO-LSSVM), and bat algorithm-LSSVM (BA-LSSVM) were utilized to forecast three compartments, respectively. Similarly, Yun et al. [8] built a robust carbon price-predicting model of complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN)-long short-term memory (LSTM) by blending the merits of CEEMDAN in decaying multiresolution time-frequency carbon price waves and the LSTM model by installing monetary waves.A report by Eckstein et al. [9] positioned Pakistan as the fifth generally impacted country because of environmental changes throughout the course of recent years, while the author in [10] positioned Pakistan’s air as the second generally contaminated in 2020. To diminish these effects, the nation needs to take broad variation and alleviation measures while changing the energy area towards dirtying and carbon impartial choices [11].Existing studies imply that greater than two-thirds of GHG emissions in Pakistan come from fossil power-related CO2 emissions. Consequently, forecasting CO2 emissions is crucial for adjusting regulations to mitigate weather changes. The forecasting of CO2 emissions has been mentioned dramatically in [12].Pao and Tsai [13] implemented the grey model to obtain CO2 emissions in Brazil between 1980 and 2007. A combination of models based on the optimization algorithms of quantum harmony and discounted mean square has been established and employed in CO2 emission forecasting for the world’s top five countries with high CO2 emissions [14]. Auffhammer and Carson [15] came up with China’s CO2 emission path using a set of panel data in a provincial-level framework [16–19]. The application of the traditional grey model and the upgraded grey model has been implemented to forecast CO2 emissions [20, 21].From qualitative and quantitative perspectives, numerous discussions have been conducted on the relationship between consumption of energy, economic growth, and CO2 emissions [22]. Al-mulali and Sab [23] investigated the influence of the consumption of basic energy and emissions of CO2 on sixteen developing countries, focusing on the impact on the development of their economies. Wang et al. [24] examined the significant factors contributing to CO2 emissions in China’s road freight transport. Zhang and Nian [25] computed the emissions of CO2 and examined the factors contributing to emissions in China’s transport sector. In China, Lin et al. [26] developed a new bootstrap ARDL bounding test for CO2 emissions and industrial development. Pao and Tsai [27] conducted an examination of the dynamism of the causal relationship, linking the emission of pollutants, consumption of energy, and their output in Britain, Russia, India, China, and other countries over a definite period [28–34].Little work has been conducted on providing mathematical and statistical modelling techniques for modelling and predicting CO2 emissions in Pakistan and where the country stands in the future concerning this important issue. It is against this backdrop that this paper provides a modelling and forecasting technique based on linear and nonlinear time series models. This study predicts the amount of CO2 emissions in Pakistan for almost a decade to assist the Pakistan government and other stakeholders to plan ahead to devise ways of mitigating the effects of CO2 emissions or possibly curtailing them.Section2 of the paper presents the source of data, data characteristics, and the modelling and forecasting techniques applied. Section 3 presents the discussion and the results, while the last section presents the conclusions and recommendations of the study.


## 2. Materials and Methods

### 2.1. Data

The data are made up of CO2 emissions in Pakistan from 1960 to 2018. The amount of CO2 was measured in metric tons per capita according to international standards. The data were obtained from the World Bank data indicators of Pakistan (https://databank.worldbank.org/). Figure 1 shows CO2 emissions (metrics tons per capita) from 1960 to 2018 as obtained from the World Bank. It can be observed from the figure that there is an increasing trend, which is an indication of a seasonal trend in the series. Seasonality must therefore be removed from the series before applying our modelling and forecasting techniques. The descriptive statistics in Table 1 show that the average CO2 emissions from 1960 to 2018 were 0.572 metric tons per capita, with a minimum CO2 emission of 0.314 metric tons per capita and a maximum of 0.981 metric tons per capita.Figure 1
Time plot of CO2 in Pakistan from 1960 to 2018 (https://databank.worldbank.org/).Table 1
Descriptive statistics of CO2 in Pakistan from 1960 to 2018.
StatisticEstimateMean0.572Median0.567Maximum0.982Minimum0.308Standard deviation0.191First quartile0.390Third quartile0.735Skewness0.259Kurtosis1.845


### 2.2. Methods

We present a brief outline of all linear and nonlinear time series models implemented in this study. In what follows, we illustrate the linear time series followed by nonlinear machine learning models.


#### 2.2.1. Linear Time Series Models

In this section, we give a brief outline of linear time series models, such as the Box–Jenkins ARIMA, trigonometric seasonality, Box–Cox transformation, ARMA errors, trend and seasonal model, naïve model, and exponential smoothing model.(1) Box–Jenkins ARIMA Modelling Approach. The ARIMA (p,d,q) [35, 36] can be represented by(1)X^t=μ+θ1x^t−1+…+θpx^t−p+…+α1et−1+…+αqet−q+et.With lagged valuesθ1y^t−1,…,θpy^t−p and lagged errors α1et−1,…,αqet−q. The p,d, and q of the series x^t are the autoregressive term order, the differencing degree of the series, and the moving average order term, respectively. The white noise et has mean 0 and variance σ2. The differencing of the series x^t can be conducted once or more times.The seasonal multiplicative ARIMA model has the form ARIMAp,d,q×P,D,Q written as(2)θpKδPEf1−Ed1−EfDX^t=∆+ϑqEμQEfεt,with θpE=1−θ1E−…θpEp; δPEf=1−δ1Ef−…−δPEPf;ϑqK=1−ϑ1E−…ϑqEq;andμQEf=1−μ1Kf−…−μQEQf.In the formula,f and K are the seasonality frequency and the shift balanced operator, respectively. The seasonal difference and ordinary differencing degrees are represented by D and d, respectively. δPEfandθpE are the seasonal autoregressive polynomial of order P and regular autoregressive polynomial of order p, respectively while μQEfandϑqE are the seasonal moving average of order Q and polynomials of the regular moving average of order q, respectively. Similarly, ∆=ρ1−θ1−…θp1−δ1−…−δP, where the mean of the process 1−Ed1−EfDX^t is ρ. The white noise εt has mean 0 and variance σ2.Figure2 shows the diagrammatical display of the steps involved in this methodology. In Figure 2, it can be noticed that there are four steps involved in ARIMA modelling. We discuss the main steps involved in this method briefly.Figure 2
Steps involved in the Box–Jenkins methodology(source: Anvari et al. [37]).Step 1.
Identification of the ARIMA Model
This step requires the series to be stationary and parameters independent of time. In time series analysis, it is often found that the series is not white noised in advance. To make the series white noised, differencing is required. To test whether the series is stationary or not, we use the augmented Dickey–Fuller test. Once the series becomes stationary, graphical tools such as the autocorrelation function (ACF) and partial autocorrelation function (PACF) are implemented to identify the order of the candidate model.Step 2.
Model Estimation
With the help of the visual display of the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the series, the appropriate candidate model for the series is estimated. Different combinations of candidate models are applied, and the final model is selected based on the criteria of accuracy of parameters.Step 3.
Diagnostic Checks
Here, the selected candidate model passes through various tests, which include the ACF plot, QQ-Norm plot, and Shapiro–Wilk to verify the normality of residuals. Once the model passes all these tests, the model is considered one of the best models for the next step.Step 4.
Forecasting
In this step, the candidate model fulfilling the three steps discussed earlier is used for predicting the future values of the series. The efficiency of the forecasted model is measured by the root mean square error (RMSE) and the mean absolute error (MAE). The mathematical expressions for RMSE and MAE are given by(3)RMSE=1t−n∑i=n+1txi−x^i2,MAE=1t−n∑i=N+1txi−x^i,
wherex1,…,xn and xn+1,…,xt are the partitions of the data for training and forecasting, respectively. In our case, we used the complete data for both training and forecasting. The RMSE and MAE assess the quality as well as summarize the model for superior forecasting. The model with the smallest RMSE and MAE values is chosen for forecasting [5]. Other indicators of measuring efficiency can also be implemented.
(2) Trigonometric seasonality, Box-Cox transformation, ARMA errors, Trend, and Seasonal (TBATS) Modelling Approach
The trigonometric seasonality, Box–Cox transformation, ARMA errors, trend, and seasonal component (TBATS) model is mainly used to predict the time series with complex seasonal patterns by exponential smoothing. The TBATS model utilizes the representation of trigonometric seasonal components based on the Fourier series. The TBATS model [38, 39] is given by(4)xtα=τt−1+βkt−1+∑i=1Tαt−nii+ϑt,τt=τt−1+βkt−1+γϑt,kt=βkt−1+ωϑt,ϑt=∑i=1pδiϑt−1+∑i=1qθiεt−1+εt,wherextα is the Box–Cox-transformed time series at moment t and αt−nii is the ith seasonal component. τt is the local level, while kt is the trend with damping. ϑt is ARMA p,q and εt is the white noise from a Gaussian process. T is the seasonality amount, whereas ni is the length of the ith seasonal period. The Box–Cox transformation is α, while β is trending damping. Also, ω, γ represent the smoothing and δi,θi are the ARMA p,q coefficients. εt is the white noise. The seasonal component of the TBATS model is given by(5)zti=∑i=1lizj,ti,zj,ti=zj,t−1icosξi+zj,t−1∗isinξi+ϕ1iϑt,zj,t∗i=−zj,t−1isinξi+zj,t−1∗icosξi+ϕ2iϑt,whereξi=2πj/ni and ϕ1i, ϕ2i are the seasonal smoothing.
(3) Naïve Modelling Approach. The naïve modelling approach utilizes the last year’s genuine data as a figure for the current year. The following year’s estimate t+1 will be equivalent to the current year’s information [33]. The supposition that its interest lies later in the period will be equivalent to the interest in the last period. The expression of the naïve method [37, 40, 41] is(6)X^t+k=X^t, The seasonal naïve model for the period t+k is of the form(7)X^t+k=X^t+k−ng+1,where n is the period of seasonality and g is the integer part of k−1/n, which is the number of years prior to the forecast year t+k.
(4) Exponential Smoothing (ETS) Modelling Approach. The exponential smoothing (ETS) technique is identically connected to one or more stochastic models. It likewise follows the property of robustness and treated as a model of critical extrapolation. The thought behind the ETS is that the estimates from this approach are the weighted normal average of past data, and weights exponentially decay with time. That is, the current data have bigger loads in contrast to the previous ones. In this way, ETS delivers savvy forecasting in contrast to typical techniques [37, 42, 43].
The simplest form of ETS [37, 42, 43] is given by(8)x^t+1t=βxt+β1−βxt−1+β1−β2xt−2+…,where β,0≤β≤1, is the parameter for smoothing and t+1 is the time for one-step-forward average forecast for the series x1,…,xt. β controls the rate of decrease in weights.
For ETS (A, N, and N), the additive model [37, 42, 43] is represented mathematically by(9)x^t=Ft−1+εt,(10)Ft=Ft−1+βεt,where (9) is the equation for the data and (10) is for the transition.
The ETS (A, N, and N) multiplicative model [37, 42, 43] is of the form(11)x^t=Ft−11+εt,Ft=Ft−11+βεt,where Ft is the new level and Ft−1 is the previews level. εt is white noise with mean 0 and variance, σ2.
We, therefore, present the nonlinear machine learning time series model after thoroughly discussing linear time series models.


#### 2.2.2. Nonlinear Machine Learning Time Series Models

Here, we briefly discuss neural network autoregressive and multilayer perceptron nonlinear machine learning time series models.(1) Neural network autoregressive (NNAR) modelling approach. A neural network autoregressive (NNAR) model is obtained when the lagged values of the time series are input to a neural network. NNAR (p, r) denotes p delayed inputs and r nodes. NNAR (p, 0) equals ARIMA (p, 0,0) without stationarity parameterization. The mathematical expression for NNAR (p,r) [36, 44–46] is of the form(12)hy=β0+∑r=1Rβrszk0+∑j=1pzrjyj,Model (12) can be constructed in two different stages, with the first being the activationR. Also, the hidden layer Br is computed as(13)Br=br=szr0+∑j=1pzrjyj,where r=1,…,Randyj=yt−1,…,yt−p.s is a previously defined nonlinear activation function. Every Br acts as a unique transformation characteristic bry. Mathematically, the hidden layer,(14)hy=β0+∑r=1RβrBr,receives R instigations.The NNAR model assumes a sigmoid function,(15)fω=11+exp⁡−ω,used in the linear function translation of 0 probability to a probability of 1.(2) Multilayer Perceptron (MLP) Model. The structure of the multilayer perceptron (MLP) model [47–50] is composed of the input layer, hidden layer, and output layer with limited nonlinear functions that are differentiable. In the MLP, information in the form of a nonlinear differentiable function is processed by artificial neurons from the input layer through the hidden layer to the output layer, yielding a response in the form of a disjointed feed-forward network algorithm. The mathematical expression of the MLP is given by(16)tx=ks∑p=0Kg1p0k∑n=0Ngpnicn+dn,where cn,dn,k,ks, tx, gpni, and g1p0 are the network inputs, bias of the network, activation function of intermediate layers, output layer activation function, output signal, weights of the intermediate layer, and connections of output neurons, respectively. Figure 3 shows the MLP model with a single hidden layer.Figure 3
Structure of multilayer perceptron (MLP) modelling with a single hidden layer.The RMSE and MAE will be used for all model comparisons.R packages, forecastHybrid, forecast, Tplyr, caTools, T-series, nnet, TTR, and ForecastTB [51–57] were used for all modelling and forecasting.


## 2.1. Data

The data are made up of CO2 emissions in Pakistan from 1960 to 2018. The amount of CO2 was measured in metric tons per capita according to international standards. The data were obtained from the World Bank data indicators of Pakistan (https://databank.worldbank.org/). Figure 1 shows CO2 emissions (metrics tons per capita) from 1960 to 2018 as obtained from the World Bank. It can be observed from the figure that there is an increasing trend, which is an indication of a seasonal trend in the series. Seasonality must therefore be removed from the series before applying our modelling and forecasting techniques. The descriptive statistics in Table 1 show that the average CO2 emissions from 1960 to 2018 were 0.572 metric tons per capita, with a minimum CO2 emission of 0.314 metric tons per capita and a maximum of 0.981 metric tons per capita.Figure 1
Time plot of CO2 in Pakistan from 1960 to 2018 (https://databank.worldbank.org/).Table 1
Descriptive statistics of CO2 in Pakistan from 1960 to 2018.
StatisticEstimateMean0.572Median0.567Maximum0.982Minimum0.308Standard deviation0.191First quartile0.390Third quartile0.735Skewness0.259Kurtosis1.845


## 2.2. Methods

We present a brief outline of all linear and nonlinear time series models implemented in this study. In what follows, we illustrate the linear time series followed by nonlinear machine learning models.


### 2.2.1. Linear Time Series Models

In this section, we give a brief outline of linear time series models, such as the Box–Jenkins ARIMA, trigonometric seasonality, Box–Cox transformation, ARMA errors, trend and seasonal model, naïve model, and exponential smoothing model.(1) Box–Jenkins ARIMA Modelling Approach. The ARIMA (p,d,q) [35, 36] can be represented by(1)X^t=μ+θ1x^t−1+…+θpx^t−p+…+α1et−1+…+αqet−q+et.With lagged valuesθ1y^t−1,…,θpy^t−p and lagged errors α1et−1,…,αqet−q. The p,d, and q of the series x^t are the autoregressive term order, the differencing degree of the series, and the moving average order term, respectively. The white noise et has mean 0 and variance σ2. The differencing of the series x^t can be conducted once or more times.The seasonal multiplicative ARIMA model has the form ARIMAp,d,q×P,D,Q written as(2)θpKδPEf1−Ed1−EfDX^t=∆+ϑqEμQEfεt,with θpE=1−θ1E−…θpEp; δPEf=1−δ1Ef−…−δPEPf;ϑqK=1−ϑ1E−…ϑqEq;andμQEf=1−μ1Kf−…−μQEQf.In the formula,f and K are the seasonality frequency and the shift balanced operator, respectively. The seasonal difference and ordinary differencing degrees are represented by D and d, respectively. δPEfandθpE are the seasonal autoregressive polynomial of order P and regular autoregressive polynomial of order p, respectively while μQEfandϑqE are the seasonal moving average of order Q and polynomials of the regular moving average of order q, respectively. Similarly, ∆=ρ1−θ1−…θp1−δ1−…−δP, where the mean of the process 1−Ed1−EfDX^t is ρ. The white noise εt has mean 0 and variance σ2.Figure2 shows the diagrammatical display of the steps involved in this methodology. In Figure 2, it can be noticed that there are four steps involved in ARIMA modelling. We discuss the main steps involved in this method briefly.Figure 2
Steps involved in the Box–Jenkins methodology(source: Anvari et al. [37]).Step 1.
Identification of the ARIMA Model
This step requires the series to be stationary and parameters independent of time. In time series analysis, it is often found that the series is not white noised in advance. To make the series white noised, differencing is required. To test whether the series is stationary or not, we use the augmented Dickey–Fuller test. Once the series becomes stationary, graphical tools such as the autocorrelation function (ACF) and partial autocorrelation function (PACF) are implemented to identify the order of the candidate model.Step 2.
Model Estimation
With the help of the visual display of the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the series, the appropriate candidate model for the series is estimated. Different combinations of candidate models are applied, and the final model is selected based on the criteria of accuracy of parameters.Step 3.
Diagnostic Checks
Here, the selected candidate model passes through various tests, which include the ACF plot, QQ-Norm plot, and Shapiro–Wilk to verify the normality of residuals. Once the model passes all these tests, the model is considered one of the best models for the next step.Step 4.
Forecasting
In this step, the candidate model fulfilling the three steps discussed earlier is used for predicting the future values of the series. The efficiency of the forecasted model is measured by the root mean square error (RMSE) and the mean absolute error (MAE). The mathematical expressions for RMSE and MAE are given by(3)RMSE=1t−n∑i=n+1txi−x^i2,MAE=1t−n∑i=N+1txi−x^i,
wherex1,…,xn and xn+1,…,xt are the partitions of the data for training and forecasting, respectively. In our case, we used the complete data for both training and forecasting. The RMSE and MAE assess the quality as well as summarize the model for superior forecasting. The model with the smallest RMSE and MAE values is chosen for forecasting [5]. Other indicators of measuring efficiency can also be implemented.
(2) Trigonometric seasonality, Box-Cox transformation, ARMA errors, Trend, and Seasonal (TBATS) Modelling Approach
The trigonometric seasonality, Box–Cox transformation, ARMA errors, trend, and seasonal component (TBATS) model is mainly used to predict the time series with complex seasonal patterns by exponential smoothing. The TBATS model utilizes the representation of trigonometric seasonal components based on the Fourier series. The TBATS model [38, 39] is given by(4)xtα=τt−1+βkt−1+∑i=1Tαt−nii+ϑt,τt=τt−1+βkt−1+γϑt,kt=βkt−1+ωϑt,ϑt=∑i=1pδiϑt−1+∑i=1qθiεt−1+εt,wherextα is the Box–Cox-transformed time series at moment t and αt−nii is the ith seasonal component. τt is the local level, while kt is the trend with damping. ϑt is ARMA p,q and εt is the white noise from a Gaussian process. T is the seasonality amount, whereas ni is the length of the ith seasonal period. The Box–Cox transformation is α, while β is trending damping. Also, ω, γ represent the smoothing and δi,θi are the ARMA p,q coefficients. εt is the white noise. The seasonal component of the TBATS model is given by(5)zti=∑i=1lizj,ti,zj,ti=zj,t−1icosξi+zj,t−1∗isinξi+ϕ1iϑt,zj,t∗i=−zj,t−1isinξi+zj,t−1∗icosξi+ϕ2iϑt,whereξi=2πj/ni and ϕ1i, ϕ2i are the seasonal smoothing.
(3) Naïve Modelling Approach. The naïve modelling approach utilizes the last year’s genuine data as a figure for the current year. The following year’s estimate t+1 will be equivalent to the current year’s information [33]. The supposition that its interest lies later in the period will be equivalent to the interest in the last period. The expression of the naïve method [37, 40, 41] is(6)X^t+k=X^t, The seasonal naïve model for the period t+k is of the form(7)X^t+k=X^t+k−ng+1,where n is the period of seasonality and g is the integer part of k−1/n, which is the number of years prior to the forecast year t+k.
(4) Exponential Smoothing (ETS) Modelling Approach. The exponential smoothing (ETS) technique is identically connected to one or more stochastic models. It likewise follows the property of robustness and treated as a model of critical extrapolation. The thought behind the ETS is that the estimates from this approach are the weighted normal average of past data, and weights exponentially decay with time. That is, the current data have bigger loads in contrast to the previous ones. In this way, ETS delivers savvy forecasting in contrast to typical techniques [37, 42, 43].
The simplest form of ETS [37, 42, 43] is given by(8)x^t+1t=βxt+β1−βxt−1+β1−β2xt−2+…,where β,0≤β≤1, is the parameter for smoothing and t+1 is the time for one-step-forward average forecast for the series x1,…,xt. β controls the rate of decrease in weights.
For ETS (A, N, and N), the additive model [37, 42, 43] is represented mathematically by(9)x^t=Ft−1+εt,(10)Ft=Ft−1+βεt,where (9) is the equation for the data and (10) is for the transition.
The ETS (A, N, and N) multiplicative model [37, 42, 43] is of the form(11)x^t=Ft−11+εt,Ft=Ft−11+βεt,where Ft is the new level and Ft−1 is the previews level. εt is white noise with mean 0 and variance, σ2.
We, therefore, present the nonlinear machine learning time series model after thoroughly discussing linear time series models.


### 2.2.2. Nonlinear Machine Learning Time Series Models

Here, we briefly discuss neural network autoregressive and multilayer perceptron nonlinear machine learning time series models.(1) Neural network autoregressive (NNAR) modelling approach. A neural network autoregressive (NNAR) model is obtained when the lagged values of the time series are input to a neural network. NNAR (p, r) denotes p delayed inputs and r nodes. NNAR (p, 0) equals ARIMA (p, 0,0) without stationarity parameterization. The mathematical expression for NNAR (p,r) [36, 44–46] is of the form(12)hy=β0+∑r=1Rβrszk0+∑j=1pzrjyj,Model (12) can be constructed in two different stages, with the first being the activationR. Also, the hidden layer Br is computed as(13)Br=br=szr0+∑j=1pzrjyj,where r=1,…,Randyj=yt−1,…,yt−p.s is a previously defined nonlinear activation function. Every Br acts as a unique transformation characteristic bry. Mathematically, the hidden layer,(14)hy=β0+∑r=1RβrBr,receives R instigations.The NNAR model assumes a sigmoid function,(15)fω=11+exp⁡−ω,used in the linear function translation of 0 probability to a probability of 1.(2) Multilayer Perceptron (MLP) Model. The structure of the multilayer perceptron (MLP) model [47–50] is composed of the input layer, hidden layer, and output layer with limited nonlinear functions that are differentiable. In the MLP, information in the form of a nonlinear differentiable function is processed by artificial neurons from the input layer through the hidden layer to the output layer, yielding a response in the form of a disjointed feed-forward network algorithm. The mathematical expression of the MLP is given by(16)tx=ks∑p=0Kg1p0k∑n=0Ngpnicn+dn,where cn,dn,k,ks, tx, gpni, and g1p0 are the network inputs, bias of the network, activation function of intermediate layers, output layer activation function, output signal, weights of the intermediate layer, and connections of output neurons, respectively. Figure 3 shows the MLP model with a single hidden layer.Figure 3
Structure of multilayer perceptron (MLP) modelling with a single hidden layer.The RMSE and MAE will be used for all model comparisons.R packages, forecastHybrid, forecast, Tplyr, caTools, T-series, nnet, TTR, and ForecastTB [51–57] were used for all modelling and forecasting.


## 2.2.1. Linear Time Series Models

In this section, we give a brief outline of linear time series models, such as the Box–Jenkins ARIMA, trigonometric seasonality, Box–Cox transformation, ARMA errors, trend and seasonal model, naïve model, and exponential smoothing model.(1) Box–Jenkins ARIMA Modelling Approach. The ARIMA (p,d,q) [35, 36] can be represented by(1)X^t=μ+θ1x^t−1+…+θpx^t−p+…+α1et−1+…+αqet−q+et.With lagged valuesθ1y^t−1,…,θpy^t−p and lagged errors α1et−1,…,αqet−q. The p,d, and q of the series x^t are the autoregressive term order, the differencing degree of the series, and the moving average order term, respectively. The white noise et has mean 0 and variance σ2. The differencing of the series x^t can be conducted once or more times.The seasonal multiplicative ARIMA model has the form ARIMAp,d,q×P,D,Q written as(2)θpKδPEf1−Ed1−EfDX^t=∆+ϑqEμQEfεt,with θpE=1−θ1E−…θpEp; δPEf=1−δ1Ef−…−δPEPf;ϑqK=1−ϑ1E−…ϑqEq;andμQEf=1−μ1Kf−…−μQEQf.In the formula,f and K are the seasonality frequency and the shift balanced operator, respectively. The seasonal difference and ordinary differencing degrees are represented by D and d, respectively. δPEfandθpE are the seasonal autoregressive polynomial of order P and regular autoregressive polynomial of order p, respectively while μQEfandϑqE are the seasonal moving average of order Q and polynomials of the regular moving average of order q, respectively. Similarly, ∆=ρ1−θ1−…θp1−δ1−…−δP, where the mean of the process 1−Ed1−EfDX^t is ρ. The white noise εt has mean 0 and variance σ2.Figure2 shows the diagrammatical display of the steps involved in this methodology. In Figure 2, it can be noticed that there are four steps involved in ARIMA modelling. We discuss the main steps involved in this method briefly.Figure 2
Steps involved in the Box–Jenkins methodology(source: Anvari et al. [37]).Step 1.
Identification of the ARIMA Model
This step requires the series to be stationary and parameters independent of time. In time series analysis, it is often found that the series is not white noised in advance. To make the series white noised, differencing is required. To test whether the series is stationary or not, we use the augmented Dickey–Fuller test. Once the series becomes stationary, graphical tools such as the autocorrelation function (ACF) and partial autocorrelation function (PACF) are implemented to identify the order of the candidate model.Step 2.
Model Estimation
With the help of the visual display of the autocorrelation function (ACF) and partial autocorrelation function (PACF) of the series, the appropriate candidate model for the series is estimated. Different combinations of candidate models are applied, and the final model is selected based on the criteria of accuracy of parameters.Step 3.
Diagnostic Checks
Here, the selected candidate model passes through various tests, which include the ACF plot, QQ-Norm plot, and Shapiro–Wilk to verify the normality of residuals. Once the model passes all these tests, the model is considered one of the best models for the next step.Step 4.
Forecasting
In this step, the candidate model fulfilling the three steps discussed earlier is used for predicting the future values of the series. The efficiency of the forecasted model is measured by the root mean square error (RMSE) and the mean absolute error (MAE). The mathematical expressions for RMSE and MAE are given by(3)RMSE=1t−n∑i=n+1txi−x^i2,MAE=1t−n∑i=N+1txi−x^i,
wherex1,…,xn and xn+1,…,xt are the partitions of the data for training and forecasting, respectively. In our case, we used the complete data for both training and forecasting. The RMSE and MAE assess the quality as well as summarize the model for superior forecasting. The model with the smallest RMSE and MAE values is chosen for forecasting [5]. Other indicators of measuring efficiency can also be implemented.
(2) Trigonometric seasonality, Box-Cox transformation, ARMA errors, Trend, and Seasonal (TBATS) Modelling Approach
The trigonometric seasonality, Box–Cox transformation, ARMA errors, trend, and seasonal component (TBATS) model is mainly used to predict the time series with complex seasonal patterns by exponential smoothing. The TBATS model utilizes the representation of trigonometric seasonal components based on the Fourier series. The TBATS model [38, 39] is given by(4)xtα=τt−1+βkt−1+∑i=1Tαt−nii+ϑt,τt=τt−1+βkt−1+γϑt,kt=βkt−1+ωϑt,ϑt=∑i=1pδiϑt−1+∑i=1qθiεt−1+εt,wherextα is the Box–Cox-transformed time series at moment t and αt−nii is the ith seasonal component. τt is the local level, while kt is the trend with damping. ϑt is ARMA p,q and εt is the white noise from a Gaussian process. T is the seasonality amount, whereas ni is the length of the ith seasonal period. The Box–Cox transformation is α, while β is trending damping. Also, ω, γ represent the smoothing and δi,θi are the ARMA p,q coefficients. εt is the white noise. The seasonal component of the TBATS model is given by(5)zti=∑i=1lizj,ti,zj,ti=zj,t−1icosξi+zj,t−1∗isinξi+ϕ1iϑt,zj,t∗i=−zj,t−1isinξi+zj,t−1∗icosξi+ϕ2iϑt,whereξi=2πj/ni and ϕ1i, ϕ2i are the seasonal smoothing.
(3) Naïve Modelling Approach. The naïve modelling approach utilizes the last year’s genuine data as a figure for the current year. The following year’s estimate t+1 will be equivalent to the current year’s information [33]. The supposition that its interest lies later in the period will be equivalent to the interest in the last period. The expression of the naïve method [37, 40, 41] is(6)X^t+k=X^t, The seasonal naïve model for the period t+k is of the form(7)X^t+k=X^t+k−ng+1,where n is the period of seasonality and g is the integer part of k−1/n, which is the number of years prior to the forecast year t+k.
(4) Exponential Smoothing (ETS) Modelling Approach. The exponential smoothing (ETS) technique is identically connected to one or more stochastic models. It likewise follows the property of robustness and treated as a model of critical extrapolation. The thought behind the ETS is that the estimates from this approach are the weighted normal average of past data, and weights exponentially decay with time. That is, the current data have bigger loads in contrast to the previous ones. In this way, ETS delivers savvy forecasting in contrast to typical techniques [37, 42, 43].
The simplest form of ETS [37, 42, 43] is given by(8)x^t+1t=βxt+β1−βxt−1+β1−β2xt−2+…,where β,0≤β≤1, is the parameter for smoothing and t+1 is the time for one-step-forward average forecast for the series x1,…,xt. β controls the rate of decrease in weights.
For ETS (A, N, and N), the additive model [37, 42, 43] is represented mathematically by(9)x^t=Ft−1+εt,(10)Ft=Ft−1+βεt,where (9) is the equation for the data and (10) is for the transition.
The ETS (A, N, and N) multiplicative model [37, 42, 43] is of the form(11)x^t=Ft−11+εt,Ft=Ft−11+βεt,where Ft is the new level and Ft−1 is the previews level. εt is white noise with mean 0 and variance, σ2.
We, therefore, present the nonlinear machine learning time series model after thoroughly discussing linear time series models.


## 2.2.2. Nonlinear Machine Learning Time Series Models

Here, we briefly discuss neural network autoregressive and multilayer perceptron nonlinear machine learning time series models.(1) Neural network autoregressive (NNAR) modelling approach. A neural network autoregressive (NNAR) model is obtained when the lagged values of the time series are input to a neural network. NNAR (p, r) denotes p delayed inputs and r nodes. NNAR (p, 0) equals ARIMA (p, 0,0) without stationarity parameterization. The mathematical expression for NNAR (p,r) [36, 44–46] is of the form(12)hy=β0+∑r=1Rβrszk0+∑j=1pzrjyj,Model (12) can be constructed in two different stages, with the first being the activationR. Also, the hidden layer Br is computed as(13)Br=br=szr0+∑j=1pzrjyj,where r=1,…,Randyj=yt−1,…,yt−p.s is a previously defined nonlinear activation function. Every Br acts as a unique transformation characteristic bry. Mathematically, the hidden layer,(14)hy=β0+∑r=1RβrBr,receives R instigations.The NNAR model assumes a sigmoid function,(15)fω=11+exp⁡−ω,used in the linear function translation of 0 probability to a probability of 1.(2) Multilayer Perceptron (MLP) Model. The structure of the multilayer perceptron (MLP) model [47–50] is composed of the input layer, hidden layer, and output layer with limited nonlinear functions that are differentiable. In the MLP, information in the form of a nonlinear differentiable function is processed by artificial neurons from the input layer through the hidden layer to the output layer, yielding a response in the form of a disjointed feed-forward network algorithm. The mathematical expression of the MLP is given by(16)tx=ks∑p=0Kg1p0k∑n=0Ngpnicn+dn,where cn,dn,k,ks, tx, gpni, and g1p0 are the network inputs, bias of the network, activation function of intermediate layers, output layer activation function, output signal, weights of the intermediate layer, and connections of output neurons, respectively. Figure 3 shows the MLP model with a single hidden layer.Figure 3
Structure of multilayer perceptron (MLP) modelling with a single hidden layer.The RMSE and MAE will be used for all model comparisons.R packages, forecastHybrid, forecast, Tplyr, caTools, T-series, nnet, TTR, and ForecastTB [51–57] were used for all modelling and forecasting.


## 3. Results and Discussion

From the time series plotted in Figure3, the series is not stationary. We verified this by applying stationarity at the level (0) and found that the series is not stationary at the level (0). This means that the mean, variance, and covariance are not stationary over some time. This violation of nonstationarity was removed by taking the first difference of the series. We implemented the Dickey–Fuller (DF) test without drift, the Phillips–Perron (PP) test without drift, to verify whether our series is stationary or not, and the structural break (Chow test) [57] test to check whether there is a structural change in the series. The Chow test indicated that there is no structural change in the series as the p value is greater than 0.05. As a result, the series becomes stationary at the level (0). The output of these tests is given in Table 2 for the level (0).Table 2
The augmented Dickey–Fuller, Phillips–Perron, and structural break (Chow test) tests for CO2 emission data at the level (0).
ItemEstimateDickey–Fuller−3.0238Lag order3p value0.161Phillips–Perron unit root test−9.6469Lag order3p value0.5404Structural break test(Chow test)p value1.00The Phillips–Perron unit root test: 9.6469.From Table2, the series is not stationary at the level (0), so we implemented a single differencing technique to make the series stationary. After making the series stationary, we searched for the best candidate model for the CO2 emission data by producing the correlogram of the differenced series. The visual autocorrelation (ACF) and partial autocorrelation (PACF) plots are presented in Figure 4.Figure 4
Correlogram of the first difference of the series of CO2 emissions.In Figure4, the ACF and PACF suggest that the estimated ARIMA model for modelling CO2 emissions includes two terms from the autoregressive (AR) part and one term from the moving average (MA) part. In the correlogram, the suitable ARIMA model for the series is ARIMA (4,1,7) since the lags of ACF and PACF are out of the boundaries on 5 and 7.Using the complete dataset, we presented the visual single horizon forecast of the ARIMA, ETS, TBATS, naïve, NNAR, and MLP models in Figure5, Figure 6, Figure 7, Figure 8, Figure 9, and Figure 10, respectively, to assess their visual performance. We traced the model performance visually to show the fitted versus the observed values. Figure 5, Figure 6, Figure 7, and Figure 8 show the visuals of linear time series modelling and forecasting, while Figure 9 and 10 show the visuals of the modelling and forecasting of nonlinear machine learning time series models.Figure 5
The ARIMA model for modelling and forecasting CO2 emissions.Figure 6
The ETS model for modelling and forecasting CO2 emissions.Figure 7
The TBATS model for modelling and forecasting CO2 emissions.Figure 8
The naïve model for modelling and forecasting CO2 emissions.Figure 9
The NNAR model for modelling and forecasting CO2 emissions.Figure 10
The MLP model for modelling and forecasting CO2 emissions.We used thennfor and neuralnet packages in R to perform the analyses for the nonlinear machine learning time series models. The packages adjust itself for the optimum number of parameters according to the behaviour of the series.Although all models showed a similar trend as the original series (observed), the nonlinear machine learning time series models had a more closer trend compared to the linear time series models. This is an indication that the nonlinear machine learning models have superior forecasting qualities than the linear time series models.Except for ETS and naïve models (Figures6 and 8) which forecasted stable CO2 emissions in subsequent years, all other models forecasted increasing CO2 emissions in Pakistan in subsequent years. This trend is worrying and devastating and similar to that observed by Butt et al. [11].We estimated the RMSE and MAE for both linear and nonlinear time series models to illustrate the numerical performance of our models. Table3 shows the RMSE and MAE estimates of all models using the complete dataset. From the table, we observe that RMSE values decrease in the order TBATS, naïve, ETS, ARIMA (4, 1, 7), MLP, and NNAR (1, 1), respectively, with the TBATS having the largest value and NNAR (1, 1) having the least value. For MAE values, the decreasing order is of the form TBATS, naïve, ETS, ARIMA (4, 1, 7), MLP, and NNAR (1, 1), respectively, with the TBATS having the largest value and NNAR (1, 1) having the least value. We noticed that both nonlinear machine learning time series models had smaller RSME and MAE values than all linear time series models, indicating their superiority to linear time series models. The NNAR (1, 1) model’s least RSME and MAE values of 0.0007 and 0.0005, respectively, is a clear evidence that the machine learning nonlinear NNAR (1, 1) model possesses a superior forecasting quality and data assessment than the rest of the time series models. As a result, we utilized the NNAR (1, 1) model for forecasting CO2 emissions in Pakistan from 2019 to 2028. The forecast values as illustrated in Figure 9 show that CO2 emissions in Pakistan will be 1.048 metric tons per capita by 2028. This depicts a rapid increase in CO2 emissions from previous years. These results are in line with those of Aftab et al. [2] and Faruque et al. [3].Table 3
RMSE and MAE estimates for linear and nonlinear time series models for the complete data.
Candidate modelsRMSEMAENaïve0.03040.0236ETS0.03020.0232ARIMA (4, 1, 7)0.02390.0189TBATS0.03080.0240NNAR (1, 1)0.00070.0005MLP0.00200.0017


## 4. Conclusion

The present study related to forecasting and modelling CO2 emissions in Pakistan clearly warns that the current level of efforts made by the Government of Pakistan is not enough to capture the changes in climate change by reducing CO2 emissions. The increasing trend in emissions is a frightening and clear warning to policymakers. The increasing CO2 emissions is a serious threat to healthy living and fuels global warming vis-à-vis climate change. We call for pragmatic policy implementation by all stakeholders to bring finality to this alarming situation. To get a clearer picture of its increasing trend and to look for an accurate prediction of future CO2 emissions in Pakistan, this study employed a linear and nonlinear time series modelling approach to study CO2 emission data in Pakistan from 1960 to 2018, namely, ARIMA, naïve, TBATS, ETS, NNAR, and MLP. These linear and nonlinear time series models were used to forecast CO2 emissions for almost a decade. We utilized the complete data to present visual single horizon forecast for all models for easy visual comparison. Apart from the ETS and naïve models, which forecasted stable CO2 emissions in subsequent years, all other models forecasted increasing CO2 emissions in Pakistan in subsequent years. Our results reflected the findings of Ullah et al. [58]. Studies in China and India showed similar CO2 emission trends [15, 59, 60]. The RMSE and the MAE were used to check and assess the quality of the forecasting of models using the complete dataset. We observed that both nonlinear machine learning time series models had smaller RMSE and MAE values than all linear time series models, indicating their superiority to them. The NNAR (1, 1) model had the least RMSE and MAE values for the complete dataset. We, therefore, chose it as our best model for forecasting. Based on the forecasted value of the NNAR (1, 1) model, Pakistan’s CO2 emissions will be 1.048 metric tons per capita by 2028. These results are in line with those of Aftab et al. [2] and Faruque et al. [3]. To achieve the targeted goal in the fight of reducing CO2 emissions, in order to accomplish the Paris Agreement of the UNFCCC, we propose that the Pakistan government should adapt the following:(i)
We should charge US$ 1000 per ton of CO2 emitted by companies and entities in every quarter of the year [61, 62]. This amount, which is set out by international authorities to curtail CO2 emissions, should be increased with time to discourage the use of high emitting fuels and encourage the use of low emitting ones. When this is initiated, it can reduce CO2 emissions by over fifty percent [61, 62].(ii)
We should produce electricity from hydro, wind, and different sources with no emissions of CO2. This can reduce CO2 emissions in excess of seventy percent as well as make electricity cheaper [61].(iii)
We should initiate rigorous planting of more trees, at least five million per year, in the populated areas of Pakistan as forest covers to absorb a greater chunk of emitted CO2 in the atmosphere. Studies reveal that five million trees in a particular area can absorb seventy percent of CO2 concentrates in the atmosphere in the area [61]. Therefore, this could reduce the concentration of CO2 in the atmosphere by over seventy percent of the recent value [61].(iv)
We should provide, as a matter of urgency, incentives to companies, organisations, institutions, and households to come out with clean technologies or use technologies with no CO2 emissions or those with lower ones [61].(v)
We should fund more studies to develop clean and innovative technologies with less or no emissions [61]. This will go a long way to reduce emissions by over forty percent.Our study can be used as a suggestion for forecasting CO2 emissions depending on different and dynamic settings. New policies for future purposes can be formulated and devised in light of this study. Our study can be extended to comparative analysis with other machine learning models using different activation functions.



---

*Source: 5903362-2023-01-31.xml*