# Script encoding: UTF-8 # R version: 3.6.3 # Author: Emmanuel Aramendia # This script provides the code for building the figures as shown in the paper: # Title: Moving from final to useful stage in energy-economy analysis: a critical assessment # doi: 10.1016/j.apenergy.2020.116194 # It assumes that the granger_causality_tests.R script has been run before and that a results data frame is available in the working environment. # The script is decomposed as follows: # (1) Importing libraries # (2) Setting up colours # (3) Tweaking the results data frame # (4) Calculating total number of valid tests using each energy metric, each gdp metric, each number of lags, each procedure... # (5) (5) General results graph # (6) Graph - Influence of the econometric procedure # (7) Graph - Influence of the energy metric # (8) Graph - Influence of the gdp variable # (9) Graph - Influence of the number of lags for useful exergy and all gdp metrics #### (1) Importing libraries #### library(viridis) # script working with version 0.5.1 library(dplyr) # script working with version 1.0.2 library(glue) # script working with version 1.4.2 #### (2) Setting up colours #### my_colours <- c("#FED976", "#FEB24C", "#FD8D3C", "#FC4E2A") #### (3) Tweaking the results data frame #### results$lags <- as.numeric(results$lags) # Setting up causality indicators # - bic stands for bilateral causality # - ey stands for causality from energy consumption growth to economic growth # - ye stands for causality from economic growth to energy consumption growth # - nc stands for no causality detected results$bic[(results$p_e_to_gdp <= 0.05) & (results$p_gdp_to_e <= 0.05)] <- 1 results$bic[!((results$p_e_to_gdp <= 0.05) & (results$p_gdp_to_e <= 0.05))] <- 0 results$ey[(results$p_e_to_gdp <= 0.05) & (results$p_gdp_to_e >= 0.05)] <- 1 results$ey[!((results$p_e_to_gdp <= 0.05) & (results$p_gdp_to_e >= 0.05))] <- 0 results$ye[(results$p_e_to_gdp >= 0.05) & (results$p_gdp_to_e <= 0.05)] <- 1 results$ye[!((results$p_e_to_gdp >= 0.05) & (results$p_gdp_to_e <= 0.05))] <- 0 results$nc[(results$p_e_to_gdp >= 0.05) && (results$p_gdp_to_e >= 0.05)] <- 1 results$nc[!((results$p_e_to_gdp >= 0.05) & (results$p_gdp_to_e >= 0.05))] <- 0 #### (4) Calculating total number of valid tests using each energy metric, each gdp metric, each number of lags, each procedure... #### # Needed for calculating percentage # Number of tests conducted n <- dplyr::filter(results, validity == 1) %>% dplyr::count() %>% dplyr::pull() # Number of tests conducted with TY procedure n_ty <- dplyr::filter(results, (validity == 1) & (granger_test == "toda-yamamoto")) %>% dplyr::count() %>% dplyr::pull() # Number of tests conducted with grangertest procedure n_granger <- dplyr::filter(results, (validity == 1) & (granger_test == "grangertest")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using final energy n_fe <- dplyr::filter(results, (validity == 1) & (energy == "g_FEG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using useful energy n_ue <- dplyr::filter(results, (validity == 1) & (energy == "g_UEG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using useful exergy n_ux <- dplyr::filter(results, (validity == 1) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using rgdpna n_rgdpna <- dplyr::filter(results, (validity == 1) & (gdp == "g_rgdpna")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using rgdpe n_rgdpe <- dplyr::filter(results, (validity == 1) & (gdp == "g_rgdpe")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using rgdpo n_rgdpo <- dplyr::filter(results, (validity == 1) & (gdp == "g_rgdpo")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 1 lag n_1 <- dplyr::filter(results, (validity == 1) & (lags == 1) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 2 lags n_2 <- dplyr::filter(results, (validity == 1) & (lags == 2) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 3 lags n_3 <- dplyr::filter(results, (validity == 1) & (lags == 3) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 4 lags n_4 <- dplyr::filter(results, (validity == 1) & (lags == 4) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 5 lags n_5 <- dplyr::filter(results, (validity == 1) & (lags == 5) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 6 lags n_6 <- dplyr::filter(results, (validity == 1) & (lags == 6) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 7 lags n_7 <- dplyr::filter(results, (validity == 1) & (lags == 7) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() # Number of tests using 8 lags n_8 <- dplyr::filter(results, (validity == 1) & (lags == 8) & (energy == "g_UXG")) %>% dplyr::count() %>% dplyr::pull() #### (5) General results graph #### # Calculations number_nc <- dplyr::filter(results, (validity == 1) & (nc == 1)) %>% dplyr::count() %>% dplyr::pull() number_ye <- dplyr::filter(results, (validity == 1) & (ye == 1)) %>% dplyr::count() %>% dplyr::pull() number_ey <- dplyr::filter(results, (validity == 1) & (ey == 1)) %>% dplyr::count() %>% dplyr::pull() number_bic <- dplyr::filter(results, (validity == 1) & (bic == 1)) %>% dplyr::count() %>% dplyr::pull() first_plot <- c(number_nc / n * 100, number_ye / n * 100, number_ey / n * 100, number_bic / n * 100) # Actual graph x11() par( mar = c(2.5, 4.1, 1.5, 2.1), cex = 0.95) x <- barplot(first_plot, names.arg = c("No causality", "Conservation", "Growth", "Feedback"), ylim = c(0,100), ylab = "Tests supporting each hypothesis (%)", col = my_colours) text(x, y = first_plot + 3, labels = c(glue("n = {number_nc}"), glue("n = {number_ye}"), glue("n = {number_ey}"), glue("n = {number_bic}"))) dev.off() #### (6) Graph - Influence of the econometric procedure #### # Calculations ty <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (granger_test == "toda-yamamoto")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (granger_test == "toda-yamamoto")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (granger_test == "toda-yamamoto")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (granger_test == "toda-yamamoto")) %>% dplyr::count() ) causality_proc <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (granger_test == "causality")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (granger_test == "causality")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (granger_test == "causality")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (granger_test == "causality")) %>% dplyr::count() ) granger_proc <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (granger_test == "grangertest")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (granger_test == "grangertest")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (granger_test == "grangertest")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (granger_test == "grangertest")) %>% dplyr::count() ) procedure <- rbind(ty / n_ty * 100, causality_proc / n_granger * 100, granger_proc / n_granger * 100) # Actual graph x11() par( mar = c(2.5, 4.1, 1.5, 2.1), cex = 1.05) x <- barplot(t(as.matrix(procedure)), col = my_colours, beside = TRUE, names.arg = c("Toda-Yamamoto", "Causality", "Granger test"), ylim = c(0,100), ylab = "Tests supporting each hypothesis (%)") # t is for transposing matrix text(x, y = t(as.matrix(procedure)) + 3, labels = c(glue("{ty[[1,1]]}"), glue("{ty[[1,2]]}"), glue("{ty[[1,3]]}"), glue("{ty[[1,4]]}"), glue("{causality_proc[1,1]}"), glue("{causality_proc[1,2]}"), glue("{causality_proc[1,3]}"), glue("{causality_proc[1,4]}"), glue("{granger_proc[1,1]}"), glue("{granger_proc[1,2]}"), glue("{granger_proc[1,3]}"), glue("{granger_proc[1,4]}")), cex = 0.85) legend("topleft", legend = c("No causality", "Conservation", "Growth", "Feedback"), fill = my_colours, box.lty = 0, cex = 1) dev.off() #### (7) Graph - Influence of the energy metric #### # Calculations final_energy <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (energy == "g_FEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (energy == "g_FEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (energy == "g_FEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (energy == "g_FEG")) %>% dplyr::count() ) useful_energy <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (energy == "g_UEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (energy == "g_UEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (energy == "g_UEG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (energy == "g_UEG")) %>% dplyr::count() ) useful_exergy <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (energy == "g_UXG")) %>% dplyr::count() ) energy_var <- rbind(final_energy / n_fe * 100, useful_energy / n_ue * 100, useful_exergy / n_ux * 100) # Actual graph x11() par( mar = c(2.5, 4.1, 1.5, 2.1), cex = 1.05) x <- barplot(t(as.matrix(energy_var)), col = my_colours, beside = TRUE, names.arg = c("Final energy", "Useful energy", "Useful exergy"), ylim = c(0,100), ylab = "Tests supporting each hypothesis (%)") # t is for transposing matrix text(x, y = t(as.matrix(energy_var)) + 3, labels = c(glue("{final_energy[1,1]}"), glue("{final_energy[1,2]}"), glue("{final_energy[1,3]}"), glue("{final_energy[1,4]}"), glue("{useful_energy[1,1]}"), glue("{useful_energy[1,2]}"), glue("{useful_energy[1,3]}"), glue("{useful_energy[1,4]}"), glue("{useful_exergy[1,1]}"), glue("{useful_exergy[1,2]}"), glue("{useful_exergy[1,3]}"), glue("{useful_exergy[1,4]}")), cex = 0.85) legend("topright", legend = c("No causality", "Conservation", "Growth", "Feedback"), fill = my_colours, box.lty = 0, cex = 1) dev.off() #### (8) Graph - Influence of the gdp variable #### # Calculations rgdpna <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (gdp == "g_rgdpna")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (gdp == "g_rgdpna")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (gdp == "g_rgdpna")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (gdp == "g_rgdpna")) %>% dplyr::count() ) rgdpe <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (gdp == "g_rgdpe")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (gdp == "g_rgdpe")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (gdp == "g_rgdpe")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (gdp == "g_rgdpe")) %>% dplyr::count() ) rgdpo <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (gdp == "g_rgdpo")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (gdp == "g_rgdpo")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (gdp == "g_rgdpo")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (gdp == "g_rgdpo")) %>% dplyr::count() ) gdp_var <- rbind(rgdpna / n_rgdpna * 100, rgdpo / n_rgdpo * 100, rgdpe / n_rgdpe * 100) # Actual graph x11() par( mar = c(2.5, 4.1, 1.5, 2.1), cex = 1.05) x <- barplot(t(as.matrix(gdp_var)), col = my_colours, beside = TRUE, names.arg = c("rgdpna", "rgdpo", "rgdpe"), ylim = c(0,100), ylab = "Tests supporting each hypothesis (%)") # t is for transposing matrix text(x, y = t(as.matrix(gdp_var)) + 3, labels = c(glue("{rgdpna[1,1]}"), glue("{rgdpna[1,2]}"), glue("{rgdpna[1,3]}"), glue("{rgdpna[1,4]}"), glue("{rgdpo[1,1]}"), glue("{rgdpo[1,2]}"), glue("{rgdpo[1,3]}"), glue("{rgdpo[1,4]}"), glue("{rgdpe[1,1]}"), glue("{rgdpe[1,2]}"), glue("{rgdpe[1,3]}"), glue("{rgdpe[1,4]}")), cex = 0.85) legend("topright", legend = c("No causality", "Conservation", "Growth", "Feedback"), fill = my_colours, box.lty = 0, cex = 1) dev.off() #### (9) Graph - Influence of the number of lags for useful exergy and all gdp metrics #### # Calculations lag_1 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 1) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 1) & (energy == "g_UXG")) %>% dplyr::count() ) lag_2 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 2) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 2) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 2) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 2) & (energy == "g_UXG")) %>% dplyr::count() ) lag_3 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 3) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 3) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 3) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 3) & (energy == "g_UXG")) %>% dplyr::count() ) lag_4 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 4) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 4) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 4) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 4) & (energy == "g_UXG")) %>% dplyr::count() ) lag_5 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 5) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 5) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 5) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 5) & (energy == "g_UXG")) %>% dplyr::count() ) lag_6 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 6) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 6) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 6) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 6) & (energy == "g_UXG")) %>% dplyr::count() ) lag_7 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 7) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 7) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 7) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 7) & (energy == "g_UXG")) %>% dplyr::count() ) lag_8 <- data.frame( results %>% dplyr::filter((validity ==1) & (nc ==1) & (lags == 8) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ye ==1) & (lags == 8) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (ey ==1) & (lags == 8) & (energy == "g_UXG")) %>% dplyr::count(), results %>% dplyr::filter((validity ==1) & (bic ==1) & (lags == 8) & (energy == "g_UXG")) %>% dplyr::count() ) number_lags <- rbind(lag_1 / n_1 * 100, lag_2 / n_2 * 100, lag_3 / n_3 * 100, lag_4 / n_4 * 100, lag_5 / n_5 * 100, lag_6 / n_6 * 100, lag_7 / n_7 * 100, lag_8 / n_8 * 100) # Actual graph x11() par( mar = c(2.5, 4.1, 1.5, 2.1), cex = 1.15) x <- barplot(t(as.matrix(number_lags)), col = my_colours, beside = TRUE, names.arg = c("1 lag", "2 lags", "3 lags", "4 lags", "5 lags", "6 lags", "7 lags", "8 lags"), ylim = c(0,105), ylab = "Tests supporting each hypothesis (%)") # t is for transposing matrix text(x, y = t(as.matrix(number_lags)) + 3, labels = c(glue("{lag_1[1,1]}"), glue("{lag_1[1,2]}"), glue("{lag_1[1,3]}"), glue("{lag_1[1,4]}"), glue("{lag_2[1,1]}"), glue("{lag_2[1,2]}"), glue("{lag_2[1,3]}"), glue("{lag_2[1,4]}"), glue("{lag_3[1,1]}"), glue("{lag_3[1,2]}"), glue("{lag_3[1,3]}"), glue("{lag_3[1,4]}"), glue("{lag_4[1,1]}"), glue("{lag_4[1,2]}"), glue("{lag_4[1,3]}"), glue("{lag_4[1,4]}"), glue("{lag_5[1,1]}"), glue("{lag_5[1,2]}"), glue("{lag_5[1,3]}"), glue("{lag_5[1,4]}"), glue("{lag_6[1,1]}"), glue("{lag_6[1,2]}"), glue("{lag_6[1,3]}"), glue("{lag_6[1,4]}"), glue("{lag_7[1,1]}"), glue("{lag_7[1,2]}"), glue("{lag_7[1,3]}"), glue("{lag_7[1,4]}"), glue("{lag_8[1,1]}"), glue("{lag_8[1,2]}"), glue("{lag_8[1,3]}"), glue("{lag_8[1,4]}")), cex = 0.85) legend("topright", legend = c("No causality", "Conservation", "Growth", "Feedback"), fill = my_colours, box.lty = 0, cex = 1) dev.off()