######Import of data
library(readr)
CHN_IND_CO2 <- read_csv("Desktop/CHN_IND_CO2.csv")
View(CHN_IND_CO2)
###### Useful R-function
library(forecast)
library(lmtest) 
library(tseries)
library(MLmetrics)
library(car)
################UNIVARIATE################
##########CHINA##########
#time-series
ts_CHN <- ts(CHN_IND_CO2$CHN, start= c(1961, 1),frequency = 1)
ts_CHN
#stationarity
adf.test(ts_CHN)
#autocorrelation functions
acf(ts_CHN)
#First order difference and tests
CO2_CHN <- diff(log(ts_CHN))
adf.test(CO2_CHN)
pp.test(CO2_CHN)

#outliers_CHN
boxplot(CO2_CHN, main = "Boxplot diff(log(CO2 emissions China))")
boxplot.stats(CO2_CHN)$out
print(CO2_CHN)

median_value_CHN <- median(CO2_CHN, na.rm = TRUE) # median value
CO2_CHINA <- CO2_CHN
CO2_CHINA[1] <- median_value_CHN
CO2_CHINA[6] <- median_value_CHN
CO2_CHINA[8] <- median_value_CHN
CO2_CHINA[9] <- median_value_CHN

#arima CHINA / tests / forecasting
acf(CO2_CHINA)
pacf(CO2_CHINA)
arima_CHN <- auto.arima(CO2_CHINA, d=0, seasonal =FALSE, stepwise = FALSE, approximation = FALSE)
summary(arima_CHN)
test_CHN <- arima(CO2_CHINA,order = c(1,0,2))
summary(test_CHN)
test_CHN_2 <- arima(CO2_CHINA,order = c(1,0,1))
summary(test_CHN_2)

#test for white noise
checkresiduals(arima_CHN) # no autocorr and normal distri + portmanteau test (ljung box) with non significant p value
#suggests for white noise --> portmanteau test: check theory
#autocorrelation
acf(arima_CHN$residuals) #within the threshold which suggests no auto correlation
#normality additional test
jarque.bera.test(arima_CHN$residuals) #not significant so we don't reject H0 that the residuals follow a normal distribution
#homoscedasticity test
library(lmtest)
residuals_AR_CHN <- lm(arima_CHN$residuals ~ 1)
gqtest(residuals_AR_CHN)#p-value not significant means we can't reject H0 of homoscedasticity
#structural breaks
plot(CO2_CHINA)
library(strucchange)
time <- CHN_IND_GDP$Year
time1962 <- time[-1]
sctest(CO2_CHINA ~ time1962, type="Chow", point= 14)
#WHITE NOISE

#forecast
F_AR_CHN <- forecast(arima_CHN, h = 3)
F_AR_CHN
#PLOT
# Real CO2 emissions data from 1961 to 2022
actual_CO2 <- c(570.63, 459.62, 456.78, 460.64, 500.29, 549.46, 460.23, 495.51, 607.68,
                807.95, 909.21, 968.65, 1008.29, 1028.10, 1183.21, 1226.42, 1340.83, 1492.78,
                1525.66, 1494.50, 1476.49, 1606.59, 1694.22, 1844.83, 1998.08, 2104.21, 2257.74,
                2425.89, 2463.65, 2484.85, 2606.10, 2730.79, 2921.65, 3103.74, 3361.64, 3508.82,
                3515.59, 3364.59, 3557.27, 3649.20, 3728.51, 4103.04, 4841.12, 5217.35, 5882.14,
                6494.34, 6983.58, 7501.50, 7891.09, 8620.63, 9532.41, 9779.35, 9956.38, 9998.67,
                9866.95, 9765.03, 10011.15, 10353.93, 10721.04, 10914.01, 11336.23, 11396.78)

# Forecasted logarithmic differences for 2023 to 2025
log_diff_forecast <- c(0.03341800,0.04400166,0.04798921)
# Convert log differences to actual CO2 forecasted values
forecasted_CO2 <- numeric(length(log_diff_forecast))
forecasted_CO2[1] <- actual_CO2[length(actual_CO2)] * exp(log_diff_forecast[1])
for (i in 2:length(log_diff_forecast)) {
  forecasted_CO2[i] <- forecasted_CO2[i-1] * exp(log_diff_forecast[i])
}
# Combining actual and forecasted CO2 emissions
years <- 1961:2025
CO2_emissions <- c(actual_CO2, forecasted_CO2)
# 80% confidence intervals for the forecasts
conf_intervals <- matrix(c(-0.015931945, 0.08276794, 
                           -0.008734738, 0.09673806,
                           -0.005210433, 0.10118885), 
                         ncol = 2, byrow = TRUE)
# Calculate lower and upper bounds of the confidence intervals
conf_lower <- numeric(length(forecasted_CO2))
conf_upper <- numeric(length(forecasted_CO2))
conf_lower[1] <- actual_CO2[length(actual_CO2)] * exp(conf_intervals[1, 1])
conf_upper[1] <- actual_CO2[length(actual_CO2)] * exp(conf_intervals[1, 2])
for (i in 2:length(log_diff_forecast)) {
  conf_lower[i] <- conf_lower[i-1] * exp(conf_intervals[i, 1])
  conf_upper[i] <- conf_upper[i-1] * exp(conf_intervals[i, 2])
}

# Combining actual and forecasted CO2 emissions with confidence intervals
lower_bound <- c(rep(NA, length(actual_CO2)), conf_lower)
upper_bound <- c(rep(NA, length(actual_CO2)), conf_upper)

# Plotting
library(ggplot2)

data <- data.frame(
  Year = years,
  CO2 = CO2_emissions,
  Lower = lower_bound,
  Upper = upper_bound
)

ggplot(data, aes(x = Year, y = CO2)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "purple", alpha = 0.5) +
  labs(title = "fossil CO2 Emissions of China (1961-2022) with forecast (2023-2025)", 
       x = "Year", 
       y = "CO2 Emissions (million tons)") +
  scale_x_continuous(breaks = seq(1961, 2025, by = 7), limits = c(1961, 2025)) +
  theme_minimal()

#Accuracy
accuracy(F_AR_CHN)
##########INDIA##########
#time-series
ts_IND <- ts(CHN_IND_CO2$IND, start= c(1961, 1),frequency = 1)
ts_IND

#stationarity
library(tseries)
adf.test(ts_IND)
#autocorrelation functions
acf(ts_IND)
#First order difference and tests
CO2_IND <- diff(log(ts_IND))
adf.test(CO2_IND)
pp.test(CO2_IND)

#structural breaks for India? + outliers removal
library(strucchange)
time <- CHN_IND_GDP$Year
time1962 <- time[-1]
sctest(CO2_IND ~ time1962, type="Chow", point= 58)

boxplot(CO2_IND, main = "Boxplot diff(log(CO2 emissions India))")
boxplot.stats(CO2_IND)$out
plot(CO2_IND)

median_value_IND <- median(CO2_IND, na.rm = TRUE) # median value
CO2_INDIA <- CO2_IND
CO2_INDIA[59] <- median_value_IND
adf.test(CO2_INDIA)
pp.test(CO2_INDIA)

#arima INDIA / tests / forecasting
library(forecast)
acf(CO2_INDIA)
pacf(CO2_INDIA)
arima_IND <- auto.arima(CO2_INDIA,d=0, seasonal =FALSE, stepwise = FALSE, approximation = FALSE)
summary(arima_IND)
#test for white noise
checkresiduals(arima_IND) # no autocorr and normal distri + portmanteau test (ljung box) with non significant p value
#suggests for white noise --> portmanteau test: check theory
#autocorrelation
acf(arima_IND$residuals) #within the threshold which suggests no auto correlation
#normality additional test
jarque.bera.test(arima_IND$residuals) #not significant so we don't reject H0 that the residuals follow a normal distribution
#homoscedasticity test
library(lmtest)
residuals_AR_IND <- lm(arima_IND$residuals ~ 1)
gqtest(residuals_AR_IND) #p-value not significant means we can't reject H0 of homoscedasticity
#WHITE NOISE

#forecast
F_AR_IND <- forecast(arima_IND, h = 3)
F_AR_IND

#PLOT
# Real CO2 emissions data from 1961 to 2022
actual_CO2 <- c(120.40, 132.58, 142.44, 139.49, 153.70, 159.37, 159.56, 174.07, 177.41, 181.72, 
                191.96, 203.04, 209.09, 215.85, 234.21, 244.75, 258.96, 263.15, 276.28, 291.71, 
                314.97, 325.38, 352.20, 361.56, 397.59, 426.31, 455.34, 491.69, 540.65, 578.00, 
                615.37, 655.45, 677.30, 714.06, 760.46, 823.62, 858.01, 875.77, 950.46, 977.53, 
                990.97, 1021.66, 1059.16, 1125.10, 1185.67, 1292.48, 1392.51, 1489.44, 1612.22, 
                1677.34, 1764.71, 1925.70, 1995.10, 2148.34, 2234.22, 2354.66, 2426.61, 2593.06, 
                2612.89, 2421.55, 2674.22, 2829.64)
# Forecasted logarithmic differences for 2023 to 2025
log_diff_forecast <- c(0.05809866, 0.07641965, 0.06424753)
# Convert log differences to actual CO2 forecasted values
forecasted_CO2 <- numeric(length(log_diff_forecast))
forecasted_CO2[1] <- actual_CO2[length(actual_CO2)] * exp(log_diff_forecast[1])
for (i in 2:length(log_diff_forecast)) {
  forecasted_CO2[i] <- forecasted_CO2[i-1] * exp(log_diff_forecast[i])
}
# Combining actual and forecasted CO2 emissions
years <- 1961:2025
CO2_emissions <- c(actual_CO2, forecasted_CO2)
# 80% confidence intervals for the forecasts
conf_intervals <- matrix(c(0.02561702, 0.09058030, 
                           0.04330396, 0.10953535, 
                           0.03066559, 0.09782946), 
                         ncol = 2, byrow = TRUE)
# Calculate lower and upper bounds of the confidence intervals
conf_lower <- numeric(length(forecasted_CO2))
conf_upper <- numeric(length(forecasted_CO2))
conf_lower[1] <- actual_CO2[length(actual_CO2)] * exp(conf_intervals[1, 1])
conf_upper[1] <- actual_CO2[length(actual_CO2)] * exp(conf_intervals[1, 2])
for (i in 2:length(log_diff_forecast)) {
  conf_lower[i] <- conf_lower[i-1] * exp(conf_intervals[i, 1])
  conf_upper[i] <- conf_upper[i-1] * exp(conf_intervals[i, 2])
}

# Combining actual and forecasted CO2 emissions with confidence intervals
lower_bound <- c(rep(NA, length(actual_CO2)), conf_lower)
upper_bound <- c(rep(NA, length(actual_CO2)), conf_upper)

# Plotting
library(ggplot2)

data <- data.frame(
  Year = years,
  CO2 = CO2_emissions,
  Lower = lower_bound,
  Upper = upper_bound
)

ggplot(data, aes(x = Year, y = CO2)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "purple", alpha = 0.5) +
  labs(title = "fossil CO2 Emissions of India (1961-2022) with forecast (2023-2025)", 
       x = "Year", 
       y = "CO2 Emissions (million tons)") +
  scale_x_continuous(breaks = seq(1961, 2025, by = 7), limits = c(1961, 2025)) +
  theme_minimal()
#Accuracy
accuracy(F_AR_IND)
################MULTIVARIATE################
############### CHINA ###############
# Data imports
library(readxl)
CHINA_MULT1 <- read_excel("Desktop/CHINA_MULT1.xlsx", 
                          range = "A39:AY101")
View(CHINA_MULT1)
#time series / stationarity / fit with CO2 WITHOUT OUTLIERS
#land sqm -> no: I(2)
ts_Land_sqm_CHN <- ts(CHINA_MULT1$`Agricultural land (sq. km)`, start= c(1961, 1),end=c(2021, 1),frequency = 1)
vector_Land_sqm_CHN <- as.numeric(ts_Land_sqm_CHN)
adf.test(vector_Land_sqm_CHN)
land_sqm_CHN <- diff(log(vector_Land_sqm_CHN))
adf.test(land_sqm_CHN)
land_sqm_CHN_2 <- diff(diff(log(vector_Land_sqm_CHN)))

CO2_CHN_MINUS1 <- head(CO2_CHN, -1)
model1_CHN <- lm(CO2_CHN_MINUS1 ~ land_sqm_CHN) 
summary(model1_CHN)
#agricult -> yes
ts_AGRI_CHN <- ts(CHINA_MULT1$`Agriculture, forestry, and fishing, value added (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_AGRI_CHN)
agri_CHN <- diff(log(ts_AGRI_CHN))
adf.test(agri_CHN)

model2_CHN <- lm(CO2_CHN ~ agri_CHN) 
summary(model2_CHN)
#fish -> no: I(2)
ts_FISH_CHN <- ts(CHINA_MULT1$`Capture fisheries production (metric tons)`, start= c(1961, 1),end=c(2021, 1),frequency = 1)
adf.test(ts_FISH_CHN)
vector_FISH_CHN <- as.numeric(ts_FISH_CHN)
fish_CHN <- diff(log(vector_FISH_CHN))
adf.test(fish_CHN)
fish_CHN_2 <- diff(diff(log(vector_FISH_CHN)))

CO2_CHN_MINUS1 <- head(CO2_CHN, -1)
model3_CHN <- lm(CO2_CHN_MINUS1 ~ fish_CHN) 
summary(model3_CHN)
#cereal -> yes
ts_Cereal_CHN <- ts(CHINA_MULT1$`Cereal production (metric tons)`, start= c(1961, 1),frequency = 1)
adf.test(ts_Cereal_CHN)
cereal_CHN <- diff(log(ts_Cereal_CHN))
adf.test(cereal_CHN)

model4_CHN <- lm(CO2_CHN ~ cereal_CHN) 
summary(model4_CHN)
#chg inventories -> yes
ts_chg_inventories_ADJ_CHN <- ts(CHINA_MULT1$`changes in inventories, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_chg_inventories_ADJ_CHN)
chg_inv_CHN <- diff(log(ts_chg_inventories_ADJ_CHN))
adf.test(chg_inv_CHN)

model5_CHN <- lm(CO2_CHN ~ chg_inv_CHN) 
summary(model5_CHN)
#crop prod -> yes
ts_crop_pod_CHN <- ts(CHINA_MULT1$`Crop production index (2014-2016 = 100)`, start= c(1961, 1),frequency = 1)
adf.test(ts_crop_pod_CHN)
crop_prod_CHN <- diff(log(ts_crop_pod_CHN))
adf.test(crop_prod_CHN)

model6_CHN <- lm(CO2_CHN ~ crop_prod_CHN) 
summary(model6_CHN)
#dec_chg -> no: I(2)
ts_DEC_CHN <- ts(CHINA_MULT1$`DEC alternative conversion factor (LCU per US$)`, start= c(1961, 1),frequency = 1)
adf.test(ts_DEC_CHN)
dec_chg_CHN <- diff(log(ts_DEC_CHN))
adf.test(dec_chg_CHN)
dec_chg_CHN_2 <- diff(diff(log(ts_DEC_CHN)))

model7_CHN <- lm(CO2_CHN ~ dec_chg_CHN) 
summary(model7_CHN)
#exp adj -> no: I(2)
ts_EXP_ADJ_CHN <- ts(CHINA_MULT1$`Exports, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_EXP_ADJ_CHN)
exp_ADJ_CHN <- diff(log(ts_EXP_ADJ_CHN))
adf.test(exp_ADJ_CHN)
exp_ADJ_CHN_2 <- diff(diff(log(ts_EXP_ADJ_CHN)))

model8_CHN <- lm(CO2_CHN ~ exp_ADJ_CHN) 
summary(model8_CHN)
#ext balance adj -> yes
ts_ext_balance_ADJ_CHN <- ts(CHINA_MULT1$`External balance, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_ext_balance_ADJ_CHN)
ext_balance_ADJ_CHN <- diff(ts_ext_balance_ADJ_CHN)
adf.test(ext_balance_ADJ_CHN)

model9_CHN <- lm(CO2_CHN ~ ext_balance_ADJ_CHN) 
summary(model9_CHN)
#fert -> yes
ts_FERT_CHN <- ts(CHINA_MULT1$`Fertilizer consumption (kilograms per hectare of arable land)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
adf.test(ts_FERT_CHN)
vector_FERT_CHN <- as.numeric(ts_FERT_CHN)
FERT_CHN <- diff(log(vector_FERT_CHN))
adf.test(FERT_CHN)

CO2_CHN_MINUS3 <- head(CO2_CHN, -1)
model10_CHN <- lm(CO2_CHN_MINUS3 ~ FERT_CHN) 
summary(model10_CHN)
#final cons exp adj-> no: I(2) 
ts_final_cons_exp_ADJ_CHN <- ts(CHINA_MULT1$`Final consumption expenditure, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_final_cons_exp_ADJ_CHN)
final_cons_exp_ADJ_CHN <- diff(log(ts_final_cons_exp_ADJ_CHN))
adf.test(final_cons_exp_ADJ_CHN)
final_cons_exp_ADJ_CHN_2 <- diff(diff(log(ts_final_cons_exp_ADJ_CHN)))

model11_CHN <- lm(CO2_CHN ~ final_cons_exp_ADJ_CHN) 
summary(model11_CHN)
#gdp -> yes
ts_GDP_CHN <- ts(CHINA_MULT1$`GDP (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_GDP_CHN)
GDP_CHN <- diff(log(ts_GDP_CHN))
adf.test(GDP_CHN)

model12_CHN <- lm(CO2_CHN ~ GDP_CHN) 
summary(model12_CHN)
#gni adj -> yes
ts_GNI_ADJ_CHN <- ts(CHINA_MULT1$`GNI, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_GNI_ADJ_CHN)
GNI_ADJ_CHN <- diff(log(ts_GNI_ADJ_CHN))
adf.test(GNI_ADJ_CHN)

model13_CHN <- lm(CO2_CHN ~ GNI_ADJ_CHN) 
summary(model13_CHN)
#cap form adj-> yes
ts_Capital_Form_ADJ_CHN <- ts(CHINA_MULT1$`Gross capital formation, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_Capital_Form_ADJ_CHN)
Capital_Form_ADJ_CHN <- diff(log(ts_Capital_Form_ADJ_CHN))
adf.test(Capital_Form_ADJ_CHN)

model14_CHN <- lm(CO2_CHN ~ Capital_Form_ADJ_CHN) 
summary(model14_CHN)
#domestic savings adj -> yes
ts_domestic_savings_ADJ_CHN <- ts(CHINA_MULT1$`Gross domestic savi ngs, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_domestic_savings_ADJ_CHN)
domestic_savings_ADJ_CHN <- diff(log(ts_domestic_savings_ADJ_CHN))
adf.test(domestic_savings_ADJ_CHN)

model15_CHN <- lm(CO2_CHN ~ domestic_savings_ADJ_CHN) 
summary(model15_CHN)
#gross cap form adj -> yes
ts_GrossF_cap_form_ADJ_CHN <- ts(CHINA_MULT1$`Gross fixed capital formation, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_GrossF_cap_form_ADJ_CHN)
GrossF_cap_form_ADJ_CHN <- diff(log(ts_GrossF_cap_form_ADJ_CHN))
adf.test(GrossF_cap_form_ADJ_CHN)

model16_CHN <- lm(CO2_CHN ~ GrossF_cap_form_ADJ_CHN) 
summary(model16_CHN)
#gne adj -> yes
ts_GNE_ADJ_CHN <- ts(CHINA_MULT1$`Gross national expenditure, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_GNE_ADJ_CHN)
GNE_ADJ_CHN <- diff(log(ts_GNE_ADJ_CHN))
adf.test(GNE_ADJ_CHN)

model17_CHN <- lm(CO2_CHN ~ GNE_ADJ_CHN) 
summary(model17_CHN)
#house npish adj -> no: I(2)
ts_house_NPISH_EXP_ADJ_CHN <- ts(CHINA_MULT1$`Houselholds and Npish consumption expenditure, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_house_NPISH_EXP_ADJ_CHN)
house_NPISH_EXP_ADJ_CHN <- diff(log(ts_house_NPISH_EXP_ADJ_CHN))
adf.test(house_NPISH_EXP_ADJ_CHN)
house_NPISH_EXP_ADJ_CHN_2 <- diff(diff(log(ts_house_NPISH_EXP_ADJ_CHN)))

model18_CHN <- lm(CO2_CHN ~ house_NPISH_EXP_ADJ_CHN) 
summary(model18_CHN)
#imports adj -> yes
ts_IMPORTS_ADJ_CHN <- ts(CHINA_MULT1$`Imports of goods, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_IMPORTS_ADJ_CHN)
IMPORTS_ADJ_CHN <- diff(log(ts_IMPORTS_ADJ_CHN))
adf.test(IMPORTS_ADJ_CHN)

model19_CHN <- lm(CO2_CHN ~ IMPORTS_ADJ_CHN) 
summary(model19_CHN)
#industry adj -> yes
ts_industry_CHN <- ts(CHINA_MULT1$`Industry (including construction), value added (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_industry_CHN)
industry_CHN <- diff(log(ts_industry_CHN))
adf.test(industry_CHN)

model20_CHN <- lm(CO2_CHN ~ industry_CHN) 
summary(model20_CHN)
#inflation -> yes
ts_infl_CHN <- ts(CHINA_MULT1$`Inflation, GDP deflator (annual %)`, start= c(1961, 1),frequency = 1)
adf.test(ts_infl_CHN)
infl_CHN <- diff(ts_infl_CHN)
adf.test(infl_CHN)

model21_CHN <- lm(CO2_CHN ~ infl_CHN) 
summary(model21_CHN)
#land cereal-> yes
ts_land_cereal_CHN <- ts(CHINA_MULT1$`Land under cereal production (hectares)`, start= c(1961, 1),frequency = 1)
adf.test(ts_land_cereal_CHN)
land_cereal_CHN <- diff(log(ts_land_cereal_CHN))
adf.test(land_cereal_CHN)

model22_CHN <- lm(CO2_CHN ~ land_cereal_CHN) 
summary(model22_CHN)
#merch exports adj -> yes
ts_merch_exp_CHN <- ts(CHINA_MULT1$`Exports of marchandise, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_merch_exp_CHN)
merch_exp_CHN <- diff(log(ts_merch_exp_CHN))
adf.test(merch_exp_CHN)

model23_CHN <- lm(CO2_CHN ~ merch_exp_CHN) 
summary(model23_CHN)
#merch imports adj -> yes
ts_merch_imp_CHN <- ts(CHINA_MULT1$`Imports of marchandise, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_merch_imp_CHN)
merch_imp_CHN <- diff(log(ts_merch_imp_CHN))
adf.test(merch_imp_CHN)

model24_CHN <- lm(CO2_CHN ~ merch_imp_CHN) 
summary(model24_CHN)
#exch rate -> no: I(2)
ts_exch_rate_CHN <- ts(CHINA_MULT1$`Official exchange rate (LCU per US$, period average)`, start= c(1961, 1),frequency = 1)
adf.test(ts_exch_rate_CHN)
exch_rate_CHN <- diff(log(ts_exch_rate_CHN))
adf.test(exch_rate_CHN)
exch_rate_CHN_2 <- diff(diff(log(ts_exch_rate_CHN)))

model25_CHN <- lm(CO2_CHN ~ exch_rate_CHN) 
summary(model25_CHN)
#cropland -> yes
ts_cropland_CHN <- ts(CHINA_MULT1$`Permanent cropland (% of land area)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_cropland_CHN <- as.numeric(ts_cropland_CHN)
adf.test(vector_cropland_CHN)
cropland_CHN <- diff(log(vector_cropland_CHN))
adf.test(cropland_CHN)

CO2_CHN_MINUS3 <- head(CO2_CHN, -1)
model26_CHN <- lm(CO2_CHN_MINUS3 ~ cropland_CHN) 
summary(model26_CHN)
#fresh water -> yes
ts_fresh_water_CHN <- ts(CHINA_MULT1$`Renewable internal freshwater resources per capita (cubic meters)`, start= c(1961, 1),end=c(2020,1),frequency = 1)
vector_fresh_water_CHN <-as.numeric(ts_fresh_water_CHN)
adf.test(vector_fresh_water_CHN)
fresh_water_CHN <- diff(log(vector_fresh_water_CHN))
adf.test(fresh_water_CHN)

CO2_CHN_MINUS4 <- head(CO2_CHN, -2)
model27_CHN <- lm(CO2_CHN_MINUS4 ~ fresh_water_CHN) 
summary(model27_CHN)
#rural pop -> yes
ts_rural_pop_CHN <- ts(CHINA_MULT1$`Rural population (% of total population)`, start= c(1961, 1),end=c(2022,1), frequency = 1)
adf.test(ts_rural_pop_CHN)
rural_pop_CHN <- diff(log(ts_rural_pop_CHN))

model28_CHN <- lm(CO2_CHN ~ rural_pop_CHN) 
summary(model28_CHN)
#services adj -> no: I(2)
ts_services_ADJ_CHN <- ts(CHINA_MULT1$`Services, value added, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_services_ADJ_CHN)
services_ADJ_CHN <- diff(log(ts_services_ADJ_CHN))
adf.test(services_ADJ_CHN)
services_ADJ_CHN_2 <- diff(diff(log(ts_services_ADJ_CHN)))

model29_CHN <- lm(CO2_CHN ~ services_ADJ_CHN) 
summary(model29_CHN)
#fish prod -> no: I(2)
ts_fish_prod_CHN <- ts(CHINA_MULT1$`Total fisheries production (metric tons)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_fish_prod_CHN <-as.numeric(ts_fish_prod_CHN)
adf.test(vector_fish_prod_CHN)
fish_prod_CHN <- diff(log(vector_fish_prod_CHN))
adf.test(fish_prod_CHN)
fish_prod_CHN_2 <- diff(diff(log(vector_fish_prod_CHN)))

CO2_CHN_MINUS1 <- head(CO2_CHN, -1)
model30_CHN <- lm(CO2_CHN_MINUS1 ~ fish_prod_CHN) 
summary(model30_CHN)
#trade -> yes
ts_trade_CHN <- ts(CHINA_MULT1$`Trade (% of GDP)`, start= c(1961, 1),frequency = 1)
adf.test(ts_trade_CHN)
trade_CHN <- diff(log(ts_trade_CHN))
adf.test(trade_CHN)

model31_CHN <- lm(CO2_CHN ~ trade_CHN) 
summary(model31_CHN)
#aqua prod -> yes
ts_aqua_prod_CHN <- ts(CHINA_MULT1$`Aquaculture production (metric tons)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_aqua_prod_CHN <-as.numeric(ts_aqua_prod_CHN)
adf.test(vector_aqua_prod_CHN)
aqua_prod_CHN <- diff(log(vector_aqua_prod_CHN))
adf.test(aqua_prod_CHN)

CO2_CHN_MINUS6 <- head(CO2_CHN, -1)
model32_CHN <- lm(CO2_CHN_MINUS6 ~ aqua_prod_CHN) 
summary(model32_CHN)
#arable land -> no: I(2)
ts_arable_land_CHN <- ts(CHINA_MULT1$`Arable land (hectares)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_arable_land_CHN <-as.numeric(ts_arable_land_CHN)
adf.test(vector_arable_land_CHN)
arable_land_CHN <- diff(log(vector_arable_land_CHN))
adf.test(arable_land_CHN)
arable_land_CHN_2 <- diff(diff(log(vector_arable_land_CHN)))

CO2_CHN_MINUS1 <- head(CO2_CHN, -1)
model33_CHN <- lm(CO2_CHN_MINUS1 ~ arable_land_CHN) 
summary(model33_CHN)

#VARIABLE(S) SELECTION
#setting up data: 1961 to 2020 data series
l_CO2_CHN <- head(CO2_CHN, -2)
l_agri_CHN <- head(agri_CHN, -2)
l_cereal_CHN <- head(cereal_CHN,-2)
l_chg_inv_CHN <- head(chg_inv_CHN, -2)
l_crop_prod_CHN <- head(crop_prod_CHN, -2)
l_fert_CHN <- head(FERT_CHN, -1)
l_gdp_CHN <- head(GDP_CHN, -2)
l_gni_CHN <- head(GNI_ADJ_CHN, -2)
l_cap_form_CHN <- head(Capital_Form_ADJ_CHN, -2)
l_dom_sav_CHN <- head(domestic_savings_ADJ_CHN, -2)
l_GrossF_cap_CHN <- head(GrossF_cap_form_ADJ_CHN, -2)
l_gne_CHN <- head(GNE_ADJ_CHN, -2)
l_imports_CHN <- head(IMPORTS_ADJ_CHN, -2)
l_industry_CHN <- head(industry_CHN, -2)
l_infl_CHN <- head(infl_CHN, -2)
l_merch_exp_CHN <- head(merch_exp_CHN, -2)
l_merch_imp_CHN <- head(merch_imp_CHN, -2)
l_cropland_CHN <- head(cropland_CHN, -1)
l_fresh_water_CHN <- fresh_water_CHN
l_rural_pop_CHN <- head(rural_pop_CHN, -2)
l_aqua_prod_CHN <- head(aqua_prod_CHN, -1)
l_land_cereal_CHN <- head(land_cereal_CHN, -2)
l_trade_CHN <- head(trade_CHN,-2)
l_ext_balance_CHN <- head(ext_balance_ADJ_CHN,-2)

#VARIABLE(S) SELECTION with Granger
grangertest(CO2_CHN~agri_CHN)
grangertest(CO2_CHN~cereal_CHN)
grangertest(CO2_CHN~chg_inv_CHN)
grangertest(CO2_CHN~crop_prod_CHN)
grangertest(l_CO2_CHN~l_fert_CHN)
grangertest(CO2_CHN~GDP_CHN)
grangertest(CO2_CHN~GNI_ADJ_CHN)
grangertest(CO2_CHN~Capital_Form_ADJ_CHN)
grangertest(CO2_CHN~domestic_savings_ADJ_CHN)
grangertest(CO2_CHN~GrossF_cap_form_ADJ_CHN)
grangertest(CO2_CHN~GNE_ADJ_CHN)
grangertest(CO2_CHN~IMPORTS_ADJ_CHN)
grangertest(CO2_CHN~industry_CHN)
grangertest(CO2_CHN~infl_CHN)
grangertest(CO2_CHN~merch_exp_CHN)
grangertest(CO2_CHN~merch_imp_CHN)
grangertest(l_CO2_CHN~l_cropland_CHN)
grangertest(l_CO2_CHN~l_fresh_water_CHN)
grangertest(CO2_CHN~rural_pop_CHN)
grangertest(l_CO2_CHN~l_aqua_prod_CHN)
grangertest(CO2_CHN~land_cereal_CHN)
grangertest(CO2_CHN~trade_CHN)
grangertest(CO2_CHN~ext_balance_ADJ_CHN)

# LASSO selection
#data set creation
df_CHN <- data.frame(co2 = l_CO2_CHN, b=l_agri_CHN,
  c=l_cereal_CHN,d=l_chg_inv_CHN,e=l_crop_prod_CHN,f=l_fert_CHN,g=l_gdp_CHN,h=l_gni_CHN,i=l_cap_form_CHN, j=l_dom_sav_CHN,k=l_GrossF_cap_CHN,
  l=l_gne_CHN,m=l_imports_CHN,n=l_industry_CHN,o=l_infl_CHN, p=l_merch_exp_CHN,q=l_merch_imp_CHN,r=l_cropland_CHN,s=l_fresh_water_CHN,t=l_rural_pop_CHN,
  u=l_aqua_prod_CHN,v=l_land_cereal_CHN,w=l_trade_CHN,x=l_ext_balance_CHN)
x_CHN <- as.matrix(df_CHN[, c("b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w","x")])  # Explanatory variables
y_CHN <- df_CHN$co2 
x_CHN <- as.matrix(x_CHN)
y_CHN <- as.numeric(y_CHN)
#LASSO model and optimal lambda
library(glmnet)
lasso_CHN <- cv.glmnet(x_CHN, y_CHN, alpha = 1)
lambda_optimal_CHN <- lasso_CHN$lambda.min
lambda_optimal_CHN
# Fit of final LASSO model using the optimal lambda
final_lasso_CHN <- glmnet(x_CHN, y_CHN, alpha = 1, lambda = lambda_optimal_CHN)
coef(final_lasso_CHN)

#variance explanation for lasso 
lasso1_CHN <- lm(l_CO2_CHN ~ l_agri_CHN+ l_industry_CHN+l_cap_form_CHN) 
summary(lasso1_CHN)

lasso2_CHN <- lm(l_CO2_CHN ~ l_agri_CHN+ l_industry_CHN+l_cap_form_CHN +l_merch_exp_CHN) 
summary(lasso2_CHN)

lasso3_CHN <- lm(l_CO2_CHN ~ l_agri_CHN+l_industry_CHN+ l_cap_form_CHN+l_merch_exp_CHN+l_gne_CHN) 
summary(lasso3_CHN)
#multicollinearity
library(car)
vif(lm(l_CO2_CHN ~ l_agri_CHN+l_industry_CHN+ l_cap_form_CHN +l_merch_exp_CHN))

#removal of outliers for granger and LASSO selected variables
outlier_gne_chn <- boxplot.stats(GNE_ADJ_CHN)$out 
outlier_gne_chn
median_gne_CHN <-median(GNE_ADJ_CHN, na.rm = TRUE)
GNE_ADJ_clean_CHN <- GNE_ADJ_CHN
GNE_ADJ_clean_CHN[GNE_ADJ_CHN %in% outlier_gne_chn] <- median_gne_CHN

outliers1_CHN <- boxplot.stats(agri_CHN)$out 
outliers1_CHN
outliers2_CHN <- boxplot.stats(industry_CHN)$out
outliers2_CHN
outliers3_CHN <- boxplot.stats(Capital_Form_ADJ_CHN)$out
outliers3_CHN
outliers4_CHN <- boxplot.stats(merch_exp_CHN)$out
outliers4_CHN

median_value1_CHN <- median(agri_CHN, na.rm = TRUE)
agri_clean_CHN <- agri_CHN
agri_clean_CHN[agri_CHN %in% outliers1_CHN] <- median_value1_CHN

median_value2_CHN <- median(industry_CHN, na.rm = TRUE)
industry_clean_CHN <- industry_CHN
industry_clean_CHN[industry_CHN %in% outliers2_CHN] <- median_value2_CHN

median_value3_CHN <- median(Capital_Form_ADJ_CHN, na.rm = TRUE)
Capital_Form_ADJ_clean_CHN <- Capital_Form_ADJ_CHN
Capital_Form_ADJ_clean_CHN[Capital_Form_ADJ_CHN %in% outliers3_CHN] <- median_value4_CHN

median_value4_CHN <- median(merch_exp_CHN, na.rm = TRUE)
merch_exp_clean_CHN <- merch_exp_CHN
merch_exp_clean_CHN[merch_exp_CHN %in% outliers4_CHN] <- median_value5_CHN

CO2_CHINA <- CO2_CHN
CO2_CHINA[1] <- median_value_CHN
CO2_CHINA[6] <- median_value_CHN
CO2_CHINA[8] <- median_value_CHN
CO2_CHINA[9] <- median_value_CHN

#VECM FOR GRANGER-CAUSAL variable
#dataset
dataset_CHN0 <- data.frame(Value1 = coredata(CO2_CHINA), Value2 = coredata(GNE_ADJ_clean_CHN))
names(dataset_CHN0) <- c("CO2_CHINA","GNE")

#LAG LENGTH determination for VECM
library(vars)
var0_CHN <- VAR(dataset_CHN0, type="const",ic="AIC")
var0_CHN$p

#number of cointegreation relationships -> 0: VECM not possible
vecm.test_CHN0 <- ca.jo(dataset_CHN0,ecdet="trend", spec= "transitory") 
summary(vecm.test_CHN0)

#VECM FOR LASSO variables
#dataset 
dataset_CHN <- data.frame(Value1 = coredata(CO2_CHINA),Value2 = coredata(agri_clean_CHN),
                                     Value3 = coredata(industry_clean_CHN),Value4 = coredata(Capital_Form_ADJ_clean_CHN),  Value5 = coredata(merch_exp_clean_CHN))
names(dataset_CHN) <- c("CO2_CHINA", "AGRI", "INDUSTRY","CAP_FORM", "MERCH_EXP")

#LAG LENGTH determination
library(vars)
var1_CHN <- VAR(dataset_CHN, type="const",ic="AIC")
var1_CHN$p
#number of cointegreation relationships
vecm.test_CHN <- ca.jo(dataset_CHN,ecdet="trend", spec= "transitory") 
summary(vecm.test_CHN)

#model estimation
library(tsDyn)
vecm_CHN <- VECM(dataset_CHN, lag = 1, r = 1, estim = "ML")
summary(vecm_CHN) #vecm model
t(vecm_CHN$model.specific$beta)

#misspecification tests - white noise test
resi_vecm_chn <- vecm_CHN$residuals[,1]
#ljung-box
checkresiduals(resi_vecm_chn) # no autocorr and normal distri + portmanteau test (ljung box) with non significant p value
#suggests for white noise --> portmanteau test
#autocorrelation
acf(resi_vecm_chn) #within the threshold which suggests no auto correlation
#normality additional test
jarque.bera.test(resi_vecm_chn) #not significant so we don't reject H0 that the residuals follow a normal distribution
#homoscedasticity test
library(lmtest)
residuals_tvecm_CHN <- lm(resi_vecm_chn ~ 1)
gqtest(residuals_tvecm_CHN)#p-value not significant means we can't reject H0 of homoscedasticity
#WHITE NOISE

#forecast
fore_CHN = predict(vecm_CHN,n.ahead = 8)
fore_CHN

#PLOT
library(ggplot2)
library(dplyr)
# Historical CO2 emissions data (1961-2022)
CO2_emissions <- c(570.63, 459.62, 456.78, 460.64, 500.29, 549.46, 460.23, 495.51, 607.68,
                   807.95, 909.21, 968.65, 1008.29, 1028.10, 1183.21, 1226.42, 1340.83, 1492.78,
                   1525.66, 1494.50, 1476.49, 1606.59, 1694.22, 1844.83, 1998.08, 2104.21, 2257.74,
                   2425.89, 2463.65, 2484.85, 2606.10, 2730.79, 2921.65, 3103.74, 3361.64, 3508.82,
                   3515.59, 3364.59, 3557.27, 3649.20, 3728.51, 4103.04, 4841.12, 5217.35, 5882.14,
                   6494.34, 6983.58, 7501.50, 7891.09, 8620.63, 9532.41, 9779.35, 9956.38, 9998.67,
                   9866.95, 9765.03, 10011.15, 10353.93, 10721.04, 10914.01, 11336.23, 11396.78)
# Forecasted differences of logarithms (2023-2030)
forecast_diffs <- c(0.013767009, 0.011900019, 0.006198723, 0.011857598, 0.010098715, 0.011424546, 0.010363445, 0.011169904)

# Convertion: differences of logarithms to levels
last_CO2 <- log(CO2_emissions[length(CO2_emissions)])
forecast_levels <- exp(cumsum(forecast_diffs) + last_CO2)
# Combination historical and forecasted data
years <- 1961:2030
CO2_all <- c(CO2_emissions, forecast_levels)
# Extractaction residuals from the VECM model and 95% confidence intervals
residuals_vecm <- residuals(vecm_CHN)
error_margin <- qt(0.8,276) * sd(residuals_vecm[, "CO2_CHINA"]) #80% IC, df=5*61 observations -29 parameters

lower_bound <- exp(cumsum(forecast_diffs - error_margin) + last_CO2)
upper_bound <- exp(cumsum(forecast_diffs + error_margin) + last_CO2)
lower_bound <- c(rep(NA, length(CO2_emissions)), lower_bound)
upper_bound <- c(rep(NA, length(CO2_emissions)), upper_bound)

data <- data.frame(
  Year = years,
  CO2 = CO2_all,
  Lower = lower_bound,
  Upper = upper_bound
)

# Plot the data
ggplot(data, aes(x = Year, y = CO2)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "orange", alpha = 0.5) +
  labs(title = "CO2 Emissions of China (1961-2022) with forecast (2023-2030)", 
       x = "Year", 
       y = "CO2 Emissions (million tons)") +
  scale_x_continuous(breaks = seq(1961, 2030, by = 6), limits = c(1961, 2030)) +
  theme_minimal()

#accuracy
fitted_values_CHN <- fitted(vecm_CHN)
fitted_values_CHN_CO2 <- fitted_values_CHN[,1]
actual_values_CHN <- c(CO2_CHINA[-c(1, 2)]) 

mae_CHN <- mean(abs(actual_values_CHN - fitted_values_CHN_CO2))
rmse_CHN <- sqrt(mean((actual_values_CHN - fitted_values_CHN_CO2)^2))
mape_CHN <- mean(abs((actual_values_CHN - fitted_values_CHN_CO2) / actual_values_CHN)) * 100

print(paste("MAE: ", mae_CHN))
print(paste("RMSE: ", rmse_CHN))
print(paste("MAPE: ", mape_CHN, "%"))

############### INDIA ###############
#import
library(readxl)
IND_MULT1 <- read_excel("Desktop/IND_MULT1.xlsx", 
                        range = "A40:AO102")
View(IND_MULT1)

#time series / stationarity / fit with CO2 WITHOUT OUTLIERS
#land sqm -> yes
ts_Land_sqm_IND <- ts(IND_MULT1$`Agricultural land (sq. km)`, start= c(1961, 1),end=c(2021, 1),frequency = 1)
vector_Land_sqm_IND <- as.numeric(ts_Land_sqm_IND)
adf.test(vector_Land_sqm_IND)
land_sqm_IND <- diff(log(vector_Land_sqm_IND))
adf.test(land_sqm_IND)

CO2_IND_MINUS1 <- head(CO2_IND, -1)

model1_IND <- lm(CO2_IND_MINUS1 ~ land_sqm_IND) 
summary(model1_IND)
#agricult -> yes
ts_AGRI_IND <- ts(IND_MULT1$`Agriculture, forestry, and fishing, value added (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_AGRI_IND)
agri_IND <- diff(log(ts_AGRI_IND))
adf.test(agri_IND)

model2_IND <- lm(CO2_IND ~ agri_IND) 
summary(model2_IND)
#fish -> yes
ts_FISH_IND <- ts(IND_MULT1$`Capture fisheries production (metric tons)`, start= c(1961, 1),end=c(2021, 1),frequency = 1)
adf.test(ts_FISH_IND)
vector_FISH_IND <- as.numeric(ts_FISH_IND)
fish_IND <- diff(log(vector_FISH_IND))
adf.test(fish_IND)

model3_IND <- lm(CO2_IND_MINUS1 ~ fish_IND) 
summary(model3_IND)
#cereal -> yes
ts_Cereal_IND <- ts(IND_MULT1$`Cereal production (metric tons)`, start= c(1961, 1),frequency = 1)
adf.test(ts_Cereal_IND)
cereal_IND <- diff(log(ts_Cereal_IND))
adf.test(cereal_IND)

model4_IND <- lm(CO2_IND ~ cereal_IND) 
summary(model4_IND)
#chg inventories -> yes
ts_chg_inventories_IND <- ts(IND_MULT1$`Changes in inventories (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_chg_inventories_IND)
chg_inv_IND <- diff(ts_chg_inventories_IND)
adf.test(chg_inv_IND)

model5_IND <- lm(CO2_IND ~ chg_inv_IND) 
summary(model5_IND)
#crop prod -> yes
ts_crop_pod_IND <- ts(IND_MULT1$`Crop production index (2014-2016 = 100)`, start= c(1961, 1),frequency = 1)
adf.test(ts_crop_pod_IND)
crop_prod_IND <- diff(log(ts_crop_pod_IND))
adf.test(crop_prod_IND)

model6_IND <- lm(CO2_IND ~ crop_prod_IND) 
summary(model6_IND)
#dec_chg -> yes
ts_DEC_IND <- ts(IND_MULT1$`DEC alternative conversion factor (LCU per US$)`, start= c(1961, 1),frequency = 1)
ts_DEC_IND
adf.test(ts_DEC_IND)
dec_chg_IND <- diff(log(ts_DEC_IND))
adf.test(dec_chg_IND)

model7_IND <- lm(CO2_IND ~ dec_chg_IND) 
summary(model7_IND)
#exp adj -> no: I(2)
ts_EXP_ADJ_IND <- ts(IND_MULT1$`Exports of goods and services (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_EXP_ADJ_IND)
exp_ADJ_IND <- diff(log(ts_EXP_ADJ_IND))
adf.test(exp_ADJ_IND)
exp_ADJ_IND_2 <- diff(diff(log(ts_EXP_ADJ_IND)))

model8_IND <- lm(CO2_IND ~ exp_ADJ_IND) 
summary(model8_IND)
#ext balance adj -> yes
ts_ext_balance_ADJ_IND <- ts(IND_MULT1$`External balance, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_ext_balance_ADJ_IND)
ext_balance_ADJ_IND <- diff(ts_ext_balance_ADJ_IND)
adf.test(ext_balance_ADJ_IND)

model9_IND <- lm(CO2_IND ~ ext_balance_ADJ_IND) 
summary(model9_IND)
#fert -> yes
ts_FERT_IND <- ts(IND_MULT1$`Fertilizer consumption (kilograms per hectare of arable land)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
adf.test(ts_FERT_IND)
vector_FERT_IND <- as.numeric(ts_FERT_IND)
FERT_IND <- diff(log(vector_FERT_IND))
adf.test(FERT_IND)

CO2_IND_MINUS3 <- head(CO2_IND, -1)
model10_IND <- lm(CO2_IND_MINUS3 ~ FERT_IND) 
summary(model10_IND)
#final cons exp adj-> yes
ts_final_cons_exp_ADJ_IND <- ts(IND_MULT1$`Final consumption expenditure (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_final_cons_exp_ADJ_IND)
final_cons_exp_ADJ_IND <- diff(log(ts_final_cons_exp_ADJ_IND))
adf.test(final_cons_exp_ADJ_IND)

model11_IND <- lm(CO2_IND ~ final_cons_exp_ADJ_IND) 
summary(model11_IND)
#gdp -> yes
ts_GDP_IND <- ts(IND_MULT1$`GDP (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_GDP_IND)
GDP_IND <- diff(log(ts_GDP_IND))
adf.test(GDP_IND)

model12_IND <- lm(CO2_IND ~ GDP_IND) 
summary(model12_IND)
#gni adj -> yes
ts_GNI_ADJ_IND <- ts(IND_MULT1$`GNI (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_GNI_ADJ_IND)
GNI_ADJ_IND <- diff(log(ts_GNI_ADJ_IND))
adf.test(GNI_ADJ_IND)

model13_IND <- lm(CO2_IND ~ GNI_ADJ_IND) 
summary(model13_IND)
#cap form adj-> yes
ts_Capital_Form_ADJ_IND <- ts(IND_MULT1$`Gross capital formation (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_Capital_Form_ADJ_IND)
Capital_Form_ADJ_IND <- diff(log(ts_Capital_Form_ADJ_IND))
adf.test(Capital_Form_ADJ_IND)

model14_IND <- lm(CO2_IND ~ Capital_Form_ADJ_IND) 
summary(model14_IND)
#domestic savings adj -> yes
ts_domestic_savings_ADJ_IND <- ts(IND_MULT1$`Gross domestic savings (% of GDP)`, start= c(1961, 1),frequency = 1)
adf.test(ts_domestic_savings_ADJ_IND)
domestic_savings_ADJ_IND <- diff(log(ts_domestic_savings_ADJ_IND))
adf.test(domestic_savings_ADJ_IND)

model15_IND <- lm(CO2_IND ~ domestic_savings_ADJ_IND) 
summary(model15_IND)
#gross fixed cap form adj -> yes
ts_GrossF_cap_form_ADJ_IND <- ts(IND_MULT1$`Gross fixed capital formation (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_GrossF_cap_form_ADJ_IND)
GrossF_cap_form_ADJ_IND <- diff(log(ts_GrossF_cap_form_ADJ_IND))
adf.test(GrossF_cap_form_ADJ_IND)

model16_IND <- lm(CO2_IND ~ GrossF_cap_form_ADJ_IND) 
summary(model16_IND)
#gne adj -> yes
ts_GNE_ADJ_IND <- ts(IND_MULT1$`Gross national expenditure (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_GNE_ADJ_IND)
GNE_ADJ_IND <- diff(log(ts_GNE_ADJ_IND))
adf.test(GNE_ADJ_IND)

model17_IND <- lm(CO2_IND ~ GNE_ADJ_IND) 
summary(model17_IND)
#house npish adj -> yes
ts_house_NPISH_EXP_ADJ_IND <- ts(IND_MULT1$`Households and NPISHs Final consumption expenditure (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_house_NPISH_EXP_ADJ_IND)
house_NPISH_EXP_ADJ_IND <- diff(log(ts_house_NPISH_EXP_ADJ_IND))
adf.test(house_NPISH_EXP_ADJ_IND)

model18_IND<- lm(CO2_IND ~ house_NPISH_EXP_ADJ_IND) 
summary(model18_IND)
#imports adj -> yes
ts_IMPORTS_ADJ_IND <- ts(IND_MULT1$`Imports of goods and services (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_IMPORTS_ADJ_IND)
IMPORTS_ADJ_IND <- diff(log(ts_IMPORTS_ADJ_IND))
adf.test(IMPORTS_ADJ_IND)

model19_IND <- lm(CO2_IND ~ IMPORTS_ADJ_IND) 
summary(model19_IND)
#industry adj -> yes
ts_industry_IND <- ts(IND_MULT1$`Industry (including construction), value added (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_industry_IND)
industry_IND <- diff(log(ts_industry_IND))
adf.test(industry_IND)

model20_IND <- lm(CO2_IND ~ industry_IND) 
summary(model20_IND)
#inflation -> yes
ts_infl_IND <- ts(IND_MULT1$`Inflation, GDP deflator (annual %)`, start= c(1961, 1),frequency = 1)
adf.test(ts_infl_IND)
infl_IND <- diff(ts_infl_IND)
adf.test(infl_IND)

model21_IND <- lm(CO2_IND ~ infl_IND) 
summary(model21_IND)
#land cereal-> yes
ts_land_cereal_IND <- ts(IND_MULT1$`Land under cereal production (hectares)`, start= c(1961, 1),frequency = 1)
adf.test(ts_land_cereal_IND)
land_cereal_IND <- diff(log(ts_land_cereal_IND))
adf.test(land_cereal_IND)

model22_IND <- lm(CO2_IND ~ land_cereal_IND) 
summary(model22_IND)
#merch exports adj -> no: I(2)
ts_merch_exp_IND <- ts(IND_MULT1$`Exports of marchandise, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_merch_exp_IND)
merch_exp_IND <- diff(log(ts_merch_exp_IND))
adf.test(merch_exp_IND)
merch_exp_IND_2 <- diff(diff(log(ts_merch_exp_IND)))

model23_IND <- lm(CO2_IND ~ merch_exp_IND) 
summary(model23_IND)
#merch imports adj-> no: I(2)
ts_merch_imp_IND <- ts(IND_MULT1$`Imports of marchandise, ADJ`, start= c(1961, 1),frequency = 1)
adf.test(ts_merch_imp_IND)
merch_imp_IND <- diff(log(ts_merch_imp_IND))
adf.test(merch_imp_IND)
merch_imp_IND_2 <- diff(diff(log(ts_merch_imp_IND)))

model24_IND <- lm(CO2_IND ~ merch_imp_IND) 
summary(model24_IND)
#exch rate -> yes
ts_exch_rate_IND <- ts(IND_MULT1$`Official exchange rate (LCU per US$, period average)`, start= c(1961, 1),frequency = 1)
adf.test(ts_exch_rate_IND)
exch_rate_IND <- diff(log(ts_exch_rate_IND))
adf.test(exch_rate_IND)

model25_IND <- lm(CO2_IND ~ exch_rate_IND) 
summary(model25_IND)
#cropland -> no: I(2)
ts_cropland_IND <- ts(IND_MULT1$`Permanent cropland (% of land area)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_cropland_IND <- as.numeric(ts_cropland_IND)
adf.test(vector_cropland_IND)
cropland_IND <- diff(log(vector_cropland_IND))
adf.test(cropland_IND)
cropland_IND_2 <- diff(diff(log(vector_cropland_IND)))

CO2_IND_MINUS3 <- head(CO2_IND, -1)
model26_IND <- lm(CO2_IND_MINUS3 ~ cropland_IND) 
summary(model26_IND)
#fresh water -> no: I(2)
ts_fresh_water_IND <- ts(IND_MULT1$`Renewable internal freshwater resources per capita (cubic meters)`, start= c(1961, 1),end=c(2020,1),frequency = 1)
vector_fresh_water_IND <-as.numeric(ts_fresh_water_IND)
adf.test(vector_fresh_water_IND)
fresh_water_IND <- diff(log(vector_fresh_water_IND))
adf.test(fresh_water_IND)
fresh_water_IND_2 <- diff(diff(log(vector_fresh_water_IND)))

CO2_IND_MINUS4 <- head(CO2_IND, -2)
model27_IND <- lm(CO2_IND_MINUS4 ~ fresh_water_IND) 
summary(model27_IND)
#rural pop -> no: I(3)
ts_rural_pop_IND <- ts(IND_MULT1$`Rural population (% of total population)`, start= c(1961, 1), frequency = 1)
adf.test(ts_rural_pop_IND)
rural_pop_IND <- diff(log(ts_rural_pop_IND))
adf.test(rural_pop_IND)
rural_pop_IND_2 <- diff(diff(diff(log(ts_rural_pop_IND))))

model28_IND <- lm(CO2_IND ~ rural_pop_IND) 
summary(model28_IND)
#services adj -> no: I(2)
ts_services_IND <- ts(IND_MULT1$`Services, value added (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_services_IND)
services_IND <- diff(log(ts_services_IND))
adf.test(services_IND)
services_IND_2 <- diff(diff(log(ts_services_IND)))

model29_IND <- lm(CO2_IND ~ services_IND) 
summary(model29_IND)
#fish prod -> yes
ts_fish_prod_IND <- ts(IND_MULT1$`Total fisheries production (metric tons)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_fish_prod_IND <-as.numeric(ts_fish_prod_IND)
adf.test(vector_fish_prod_IND)
fish_prod_IND <- diff(log(vector_fish_prod_IND))
adf.test(fish_prod_IND)

CO2_IND_MINUS1 <- head(CO2_IND, -1)
model30_IND <- lm(CO2_IND_MINUS1 ~ fish_prod_IND) 
summary(model30_IND)
#trade -> no: I(2)
ts_trade_IND <- ts(IND_MULT1$`Trade (% of GDP)`, start= c(1961, 1),frequency = 1)
adf.test(ts_trade_IND)
trade_IND <- diff(log(ts_trade_IND))
adf.test(trade_IND)
trade_IND_2 <- diff(diff(log(ts_trade_IND)))

model31_IND <- lm(CO2_IND ~ trade_IND) 
summary(model31_IND)
#aqua prod -> no: I(2)
ts_aqua_prod_IND <- ts(IND_MULT1$`Aquaculture production (metric tons)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_aqua_prod_IND <-as.numeric(ts_aqua_prod_IND)
adf.test(vector_aqua_prod_IND)
aqua_prod_IND <- diff(log(vector_aqua_prod_IND))
adf.test(aqua_prod_IND)
aqua_prod_IND_2 <- diff(diff(log(vector_aqua_prod_IND)))

CO2_IND_MINUS6 <- head(CO2_IND, -1)
model32_IND <- lm(CO2_IND_MINUS6 ~ aqua_prod_IND) 
summary(model32_IND)
#arable land -> yes
ts_arable_land_IND <- ts(IND_MULT1$`Arable land (hectares)`, start= c(1961, 1),end=c(2021,1),frequency = 1)
vector_arable_land_IND <-as.numeric(ts_arable_land_IND)
adf.test(vector_arable_land_IND)
arable_land_IND <- diff(log(vector_arable_land_IND))
adf.test(arable_land_IND)

CO2_IND_MINUS6 <- head(CO2_IND, -1)
model33_IND <- lm(CO2_IND_MINUS6 ~ arable_land_IND) 
summary(model33_IND)
#exports as capacity to import -> yes
ts_exp_imp_IND <- ts(IND_MULT1$`Exports as a capacity to import (constant LCU)`, start= c(1961, 1),frequency = 1)
adf.test(ts_exp_imp_IND)
exp_imp_IND <- diff(log(ts_exp_imp_IND))
adf.test(exp_imp_IND)

model34_IND <- lm(CO2_IND ~ exp_imp_IND) 
summary(model34_IND)
#Infl, consumer prices -> yes
ts_cons_prices_IND <- ts(IND_MULT1$`Inflation, consumer prices (annual %)`, start= c(1961, 1),frequency = 1)
adf.test(ts_cons_prices_IND)
cons_prices_IND <- diff(ts_cons_prices_IND)
adf.test(cons_prices_IND)

model35_IND <- lm(CO2_IND ~ cons_prices_IND) 
summary(model35_IND)

#VARIABLE(S) SELECTION
#setting up data: 1961 to 2020 data series
l_CO2_IND <- head(CO2_IND, -2)
l_land_sqm_IND <- head(land_sqm_IND, -1)
l_agri_IND <- head(agri_IND, -2)
l_fish_IND <- head(fish_IND, -1)
l_cereal_IND <- head(cereal_IND,-2)
l_chg_inv_IND <- head(chg_inv_IND, -2)
l_crop_prod_IND <- head(crop_prod_IND, -2)
l_dec_chg_IND <- head(dec_chg_IND, -2)
l_ext_balance_ADJ_IND <-head(ext_balance_ADJ_IND,-2)
l_fert_IND <- head(FERT_IND, -1)
l_final_cons_exp_ADJ_IND <- head(final_cons_exp_ADJ_IND, -2)
l_gdp_IND <- head(GDP_IND, -2)
l_gni_IND <- head(GNI_ADJ_IND, -2)
l_cap_form_IND <- head(Capital_Form_ADJ_IND, -2)
l_dom_sav_IND <- head(domestic_savings_ADJ_IND, -2)
l_GrossF_cap_IND <- head(GrossF_cap_form_ADJ_IND, -2)
l_gne_IND <- head(GNE_ADJ_IND, -2)
l_house_NPISH_EXP_ADJ_IND <- head(house_NPISH_EXP_ADJ_IND, -2)
l_imports_IND <- head(IMPORTS_ADJ_IND, -2)
l_industry_IND <- head(industry_IND, -2)
l_infl_IND <- head(infl_IND, -2)
l_land_cereal_IND <- head(land_cereal_IND, -2)
l_exch_rate_IND <- head(exch_rate_IND, -2)
l_fish_prod_IND <- head(fish_prod_IND, -1)
l_arable_land_IND <- head(arable_land_IND, -1)
l_exp_imp_IND <- head(exp_imp_IND, -2)
l_cons_prices_IND <- head(cons_prices_IND,-2)

#VARIABLE(S) SELECTION with Granger
grangertest(l_CO2_IND~l_land_sqm_IND)
grangertest(CO2_IND~agri_IND)
grangertest(l_CO2_IND~l_fish_IND)
grangertest(CO2_IND~cereal_IND)
grangertest(CO2_IND~chg_inv_IND)
grangertest(CO2_IND~crop_prod_IND)
grangertest(CO2_IND~dec_chg_IND)
grangertest(CO2_IND~ext_balance_ADJ_IND)
grangertest(l_CO2_IND~l_fert_IND)
grangertest(l_CO2_IND~l_final_cons_exp_ADJ_IND)
grangertest(CO2_IND~GDP_IND)
grangertest(CO2_IND~GNI_ADJ_IND)
grangertest(CO2_IND~Capital_Form_ADJ_IND)
grangertest(CO2_IND~domestic_savings_ADJ_IND)
grangertest(CO2_IND~GrossF_cap_form_ADJ_IND)
grangertest(CO2_IND~GNE_ADJ_IND)
grangertest(CO2_IND~house_NPISH_EXP_ADJ_IND)
grangertest(CO2_IND~IMPORTS_ADJ_IND)
grangertest(CO2_IND~industry_IND)
grangertest(l_CO2_IND~l_infl_IND)
grangertest(CO2_IND~land_cereal_IND)
grangertest(l_CO2_IND~l_cropland_IND)
grangertest(l_CO2_IND~l_fish_prod_IND)
grangertest(l_CO2_IND~l_arable_land_IND)
grangertest(CO2_IND~exp_imp_IND)                          
grangertest(CO2_IND~cons_prices_IND)

# -> not possible to do VECM with granger methodology, we try with LASSO
#data set creation
df_IND <- data.frame(co2 = l_CO2_IND, b = l_agri_IND,
  c=l_fish_IND,d=l_cereal_IND,e=l_chg_inv_IND,f=l_crop_prod_IND,g=l_dec_chg_IND,h=l_ext_balance_ADJ_IND,i=l_fert_IND,j=l_final_cons_exp_ADJ_IND,k=l_gdp_IND,l=l_gni_IND,m=l_cap_form_IND,n=l_dom_sav_IND,o=l_GrossF_cap_IND,
  p=l_gne_IND,q=l_house_NPISH_EXP_ADJ_IND,r=l_imports_IND,s=l_industry_IND,t=l_infl_IND,u=l_land_cereal_IND,v=l_fish_prod_IND,w=l_arable_land_IND,
  x=l_exp_imp_IND,y=l_cons_prices_IND)
x_IND <- as.matrix(df_IND[, c("b","c","d","e","f","g","h","i","j","k","l","m","n","o","p","q","r","s","t","u","v","w","x","y")])  # Explanatory variables
y_IND <- df_IND$co2 
x_IND <- as.matrix(x_IND)
y_IND <- as.numeric(y_IND)
#LASSO model and optimal lambda
library(glmnet)
lasso_IND <- cv.glmnet(x_IND, y_IND, alpha = 1)
lambda_optimal_IND <- lasso_IND$lambda.min
lambda_optimal_IND

# Fit of final LASSO model using the optimal lambda
final_lasso_IND <- glmnet(x_IND, y_IND, alpha = 1, lambda = lambda_optimal_IND)
coef(final_lasso_IND)

#variance explanation for lasso 
lasso1_IND <- lm(l_CO2_IND ~ l_agri_IND+ l_ext_balance_ADJ_IND+l_industry_IND+l_dom_sav_IND) 
summary(lasso1_IND)

lasso2_IND <- lm(l_CO2_IND ~ l_agri_IND+ l_final_cons_exp_ADJ_IND+l_industry_IND+l_dom_sav_IND+l_dec_chg_IND) 
summary(lasso2_IND)

lasso3_IND <- lm(l_CO2_IND ~ l_ext_balance_ADJ_IND+l_industry_IND+l_dom_sav_IND) 
summary(lasso3_IND)
#multicollinearity
library(car)
vif(lm(l_CO2_IND ~ l_agri_IND+ l_final_cons_exp_ADJ_IND+l_industry_IND+l_dom_sav_IND+l_dec_chg_IND))

#removal of outliers for LASSO selected variables
outliers1_IND <- boxplot.stats(agri_IND)$out 
outliers1_IND
outliers2_IND <- boxplot.stats(final_cons_exp_ADJ_IND)$out
outliers2_IND
outliers3_IND <- boxplot.stats(industry_IND)$out
outliers3_IND
outliers4_IND <- boxplot.stats(domestic_savings_ADJ_IND)$out
outliers4_IND
outliers5_IND <- boxplot.stats(dec_chg_IND)$out
outliers5_IND

median_value1_IND <- median(agri_IND, na.rm = TRUE)
agri_clean_IND <- agri_IND
agri_clean_IND[agri_IND %in% outliers1_IND] <- median_value1_IND

median_value2_IND <- median(final_cons_exp_ADJ_IND, na.rm = TRUE)
final_cons_exp_ADJ_clean_IND <- final_cons_exp_ADJ_IND
final_cons_exp_ADJ_clean_IND[final_cons_exp_ADJ_IND %in% outliers2_IND] <- median_value2_IND

median_value3_IND <- median(industry_IND, na.rm = TRUE)
industry_clean_IND <- industry_IND
industry_clean_IND[industry_IND %in% outliers3_IND] <- median_value3_IND

median_value4_IND <- median(domestic_savings_ADJ_IND, na.rm = TRUE)
domestic_savings_ADJ_clean_IND <- domestic_savings_ADJ_IND
domestic_savings_ADJ_clean_IND[domestic_savings_ADJ_IND %in% outliers4_IND] <- median_value4_IND

median_value5_IND <- median(dec_chg_IND, na.rm = TRUE)
dec_chg_clean_IND <- dec_chg_IND
dec_chg_clean_IND[dec_chg_IND %in% outliers5_IND] <- median_value5_IND

median_value_IND <- median(CO2_IND, na.rm = TRUE) # median value
CO2_INDIA <- CO2_IND
CO2_INDIA[3] <- median_value_IND
CO2_INDIA[59] <- median_value_IND

#VECM MODEL: dataset creation
dataset_IND <- data.frame(Value1 = coredata(CO2_INDIA),Value2 = coredata(agri_clean_IND),
                          Value3 = coredata(final_cons_exp_ADJ_clean_IND), Value4 = coredata(industry_clean_IND), Value5 = coredata(domestic_savings_ADJ_clean_IND), Value6 = coredata(dec_chg_clean_IND))
names(dataset_IND) <- c("CO2_INDIA", "AGRI","FINAL_CONS","INDUSTRY","DOM_SAVINGS", "DEC_CHG")

#LAG LENGTH determination
library(vars)
var1_IND <- VAR(dataset_IND, type="const",ic="AIC")
var1_IND$p
#number of cointegreation relationships
vecm.test_IND <- ca.jo(dataset_IND,ecdet="trend", spec= "transitory") 
summary(vecm.test_IND)

#model estimation
library(tsDyn)
vecm_IND <- VECM(dataset_IND, lag = 1, r = 3, estim = "ML")
summary(vecm_IND) #vecm model
t(vecm_IND$model.specific$beta)

#misspecification tests - white noise test
resi_vecm_ind <- vecm_IND$residuals[,1]
#ljung-box
checkresiduals(resi_vecm_ind) # no autocorr and normal distri + portmanteau test (ljung box) with non significant p value
#suggests for white noise --> portmanteau test
#autocorrelation
acf(resi_vecm_ind) #within the threshold which suggests no auto correlation
#normality additional test
jarque.bera.test(resi_vecm_ind) #not significant so we don't reject H0 that the residuals follow a normal distribution
#homoscedasticity test
library(lmtest)
residuals_tvecm_IND <- lm(resi_vecm_ind ~ 1)
gqtest(residuals_tvecm_IND)#p-value not significant means we can't reject H0 of homoscedasticity
#WHITE NOISE

#forecast
fore_IND = predict(vecm_IND,n.ahead = 8)
fore_IND

#PLOT
library(ggplot2)
library(dplyr)
# Historical CO2 emissions data (1961-2022)
CO2_emissions <- c(120.40, 132.58, 142.44, 139.49, 153.70, 159.37, 159.56, 174.07, 
                   177.41, 181.72, 191.96, 203.04, 209.09, 215.85, 234.21, 244.75, 
                   258.96, 263.15, 276.28, 291.71, 314.97, 325.38, 352.20, 361.56, 
                   397.59, 426.31, 455.34, 491.69, 540.65, 578.00, 615.37, 655.45, 
                   677.30, 714.06, 760.46, 823.62, 858.01, 875.77, 950.46, 977.53, 
                   990.97, 1021.66, 1059.16, 1125.10, 1185.67, 1292.48, 1392.51, 
                   1489.44, 1612.22, 1677.34, 1764.71, 1925.70, 1995.10, 2148.34, 
                   2234.22, 2354.66, 2426.61, 2593.06, 2612.89, 2421.55, 2674.22, 
                   2829.64)
# Forecasted differences of logarithms (2023-2030)
forecast_diffs <- c(0.07509972, 0.07684168, 0.06284641, 0.07513904, 0.07532490, 
                    0.07081696, 0.07306676, 0.07421015)
# Convertion: differences of logarithms to levels
last_CO2 <- log(CO2_emissions[length(CO2_emissions)])
forecast_levels <- exp(cumsum(forecast_diffs) + last_CO2)
# Combination historical and forecasted data
years <- 1961:2030
CO2_all <- c(CO2_emissions, forecast_levels)
# Extract residuals from the VECM model and 95% confidence intervals
residuals_vecm <- residuals(vecm_IND)
error_margin <- qt(0.8,315) * sd(residuals_vecm[, "CO2_INDIA"]) #80% IC, df=6*61 observations -51 parameters

lower_bound <- exp(cumsum(forecast_diffs - error_margin) + last_CO2)
upper_bound <- exp(cumsum(forecast_diffs + error_margin) + last_CO2)
lower_bound <- c(rep(NA, length(CO2_emissions)), lower_bound)
upper_bound <- c(rep(NA, length(CO2_emissions)), upper_bound)

data <- data.frame(
  Year = years,
  CO2 = CO2_all,
  Lower = lower_bound,
  Upper = upper_bound
)

# Plot the data
ggplot(data, aes(x = Year, y = CO2)) +
  geom_line(color = "blue") +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "orange", alpha = 0.5) +
  labs(title = "CO2 Emissions of India (1961-2022) with forecast (2023-2030)", 
       x = "Year", 
       y = "CO2 Emissions (million tons)") +
  scale_x_continuous(breaks = seq(1961, 2030, by = 6), limits = c(1961, 2030)) +
  theme_minimal()

#accuracy
fitted_values_IND <- fitted(vecm_IND)
fitted_values_IND_CO2 <- fitted_values_IND[,1]
actual_values_IND <- c(CO2_INDIA[-c(1, 2)]) 

mae_IND <- mean(abs(actual_values_IND - fitted_values_IND_CO2))
rmse_IND <- sqrt(mean((actual_values_IND - fitted_values_IND_CO2)^2))
mape_IND <- mean(abs((actual_values_IND - fitted_values_IND_CO2) / actual_values_IND)) * 100

print(paste("MAE: ", mae_IND))
print(paste("RMSE: ", rmse_IND))
print(paste("MAPE: ", mape_IND, "%"))