The following post is from Hang Qian, Software Developer on the Econometrics Toolbox Team.
Global carbon emissions have increased dramatically since 1901. However, the annual growth rates were timevarying in the past century and slowed down in recent years, possibly due to international efforts on climate change under the Kyoto Protocol and Paris Agreement, and exogenous shocks by the pandemic.
This demo analyzes timevarying local trends of carbon emissions from fossil fuels and builds a multivariate dynamic factor model for carbon emissions from coal, gas, and oil.
Retrieve Data from Haver Analytics®
Connect to the Haver Analytics Environmental, Social, and Governance Indicators (ESG) database by specifying its database name file located on your local machine, for example, C:\DLX\DATA\ESG.dat.
h = haver(“C:\DLX\DATA\ESG.dat”);
h.DataReturnFormat = “timetable”;
The ESG database contains yearly measurements of carbon emissions from burning fossil fuels from 1901, among other series.
 C001BB — Total yearly carbon emissions from burning fossil fuels (billions of metric tons of CO2)
 C001BBTP — Carbonemission contribution from coal
 C001BBTG — Carbonemission contribution from gas
 C001BBTO — Carbonemission contribution from oil
CO2Emissions = fetch(h,[“C001BB” “C001BBTP” “C001BBTO” “C001BBTG”]);
CO2Emissions is a 4column cell array containing timetable objects with the carbonemission measurements for each retrieved series.
Visualize the Time Series
Specify variables for each time series.
dates = CO2Emissions{1}.Time;
fossil = CO2Emissions{1}.GlobalCarbonEmissionsFromFossilFuelCombustion_Industry_Bil_Metr;
coal = CO2Emissions{2}.GlobalCarbonEmissionsFromCoal_Bil_MetricTons_;
oil = CO2Emissions{3}.GlobalCarbonEmissionsFromOil_Bil_MetricTons_;
gas = CO2Emissions{4}.GlobalCarbonEmissionsFromGas_Bil_MetricTons_;
Visually inspect the variables to determine their dynamic properties.
ylabel(tl,“Carbon Emissions (bil. metric tons)”);
The time series plots indicate that carbon emissions have an upward trend, but a deterministic linear trend is inadequate for time series modeling.
The local function
paramMap is a parameter mapping function that defines the local trend model of a single variable in statespace form. Create a statespace model by passing the parameter mapping
paramMap to
ssm and estimate the unknown parameters by maximum likelihood using
estimate.
params0 = [0.1; 0.1; 0.5; 0.01];
EstMdl = estimate(Mdl,y,params0,Univariate=true,lb=zeros(4,1));
Method: Maximum likelihood (fmincon)
Sample size: 121
Logarithmic likelihood: 73.4142
Akaike info criterion: 138.828
Bayesian info criterion: 127.645
 Coeff Std Err t Stat Prob
—————————————————
c(1)  0.03895 0.00487 7.99895 0
c(2)  0.08024 0.00389 20.62245 0
c(3)  0.51501 0.98526 0.52271 0.60118
c(4)  0.02905 0.55409 0.05243 0.95818

 Final State Std Dev t Stat Prob
x(1)  10.00236 0.06387 156.60423 0
x(2)  9.96258 0.04528 220.01330 0
Based on the fitted model, the Kalman filter and smoother provide the state estimation. The first state variable corresponds to the local trend, and the second variable is the lagged variable.
States = smooth(EstMdl,y);
localTrend = States(:,1);
Forecast LocalTrend StateSpace Model
The estimated local trend facilitates the outofsample forecast. Without new observations, assume that the current state of the growth momentum passes on to the forecasting periods. The Kalman filter provides the forecasted carbon emission levels.
Forecast the local trend model into a horizon of length n. Compute approximate 95% forecast intervals. Plot the forecasts.
[yf,yfMSE] = forecast(EstMdl,n,y);
yfLB = yf – 2*sqrt(yfMSE);
yfUB = yf + 2*sqrt(yfMSE);
xf = datetime((year(dates(end))+1:year(dates(end))+n)’,1 ,1);
ylabel(“Carbon Emissions (bil. metric tons)”)
title(“InSample Fit and OutofSample Forecast”)
Dummy Variables for Intervention Effects
As countries are working towards carbon neutrality, it is of interest to consider the effects of the Kyoto Protocol and the Paris Agreement on carbon emissions. Also, exogenous shocks by the pandemic slowed down the world economy, which resulted in the reduction of greenhouse gas emissions. Create dummy variables that indicate interventions on climate.
Kyoto = year(dates) > 1997;
Paris = year(dates) > 2015;
Covid = year(dates) >= 2020;
Dummy = [Kyoto Paris Covid];
The local function paramMapDummy is a parameter mapping that includes the exogenous predictors in the observation equation of the localtrend statespace model. Create a statespace model for the localtrend model with intervention dummy variables, and then fit the model to the total carbon emissions variable fossil.
MdlDummy = ssm(@(params)paramMapDummy(params,fossil,Dummy));
params0Dummy = [params0; zeros(nvar,1)];
[~,estParams] = estimate(MdlDummy,fossil,params0Dummy, …
Univariate=true,lb=[zeros(4,1); Inf(nvar,1)]);
Method: Maximum likelihood (fmincon)
Sample size: 121
Logarithmic likelihood: 89.3099
Akaike info criterion: 164.62
Bayesian info criterion: 145.049
 Coeff Std Err t Stat Prob
—————————————————
c(1)  0.05716 0.00594 9.62038 0
c(2)  0.05491 0.00369 14.87828 0
c(3)  0.52095 0.75995 0.68550 0.49303
c(4)  0.02591 0.84557 0.03064 0.97555
c(5)  0.17017 0.32139 0.52948 0.59648
c(6)  0.13397 0.15493 0.86474 0.38719
c(7)  0.73923 0.10720 6.89565 0

 Final State Std Dev t Stat Prob
x(1)  11.09709 0.04838 229.38155 0
x(2)  10.74625 0.03542 303.37964 0
The last three parameters correspond to the intervention effects. Empirical results confirm that carbon emissions have been reduced under the Kyoto Protocol and Paris Agreement. The pandemic also has a significant negative effect on greenhouse gas emissions.
Common Local Trend for Three Series
Assume that a common stochastic trend drives global carbon emissions from coal, gas, and oil. This is a dynamic factor model, in which the common local trend serves as an unobserved factor underlying three observed time series. To prepare for the multivariate extension of the local trend model, reparametrize the local trend as
MdlFactor = ssm(@(params)paramMapFactor(params,Y));
params0Factor = [0.1*ones(3,1); 0.1*ones(3,1); 0.5; 0.01];
estimate(MdlFactor,Y,params0Factor,Univariate=true, …
lb=zeros(size(params0Factor)));
Method: Maximum likelihood (fmincon)
Sample size: 121
Logarithmic likelihood: 271.547
Akaike info criterion: 527.094
Bayesian info criterion: 504.728
 Coeff Std Err t Stat Prob
—————————————————–
c(1)  0.00499 0.00076 6.53820 0
c(2)  0.00687 0.00103 6.66628 0
c(3)  0.00335 0.00050 6.68231 0
c(4)  0.16843 0.01224 13.75993 0
c(5)  0.45911 0.03079 14.90903 0
c(6)  0.01370 0.00061 22.48582 0
c(7)  0.73514 0.06504 11.30249 0
c(8)  0.00205 0.00172 1.19137 0.23351

 Final State Std Dev t Stat Prob
x(1)  641.47037 2.89157 221.84121 0
x(2)  627.97142 2.14629 292.58522 0
Conclusion
A feature of the statespace model is to incorporate the unobserved components like local levels and trends in the state variables. As a nonstationary model, the randomwalk local trend captures the timevarying trends that provide a good fit for the global carbon emission data.
Local Functions
function [A,B,C,D,mean0,Cov0] = paramMap(params)
mean0 = [intercept;interceptslope];
function [A,B,C,D,mean0,Cov0,stateType,yDeflate] = paramMapDummy(params,y,Dummy)
dummyCoeff = params(5:end);
mean0 = [intercept;interceptslope];
yDeflate = y – Dummy*dummyCoeff;
function [A,B,C,D,mean0,Cov0,stateType,YDeflate] = paramMapFactor(params,Y)
C = [sigmaU, zeros(3,1)];
linear = [intercept;0;0] + [slope;0;0] .* (1:size(Y,1));
Comments
To leave a comment, please click here to sign in to your MathWorks Account or create a new one.