improved CO2 target plots

renv.lock

 {
   "R": {
-    "Version": "4.0.3",
+    "Version": "4.3.2",
     "Repositories": [
       {
         "Name": "REPO_NAME",
@@ -11,7 +11,7 @@
   "Packages": {
     "MASS": {
       "Package": "MASS",
-      "Version": "7.3-53",
+      "Version": "7.3-60",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
@@ -22,15 +22,16 @@
         "stats",
         "utils"
       ],
-      "Hash": "d1bc1c8e9c0ace57ec9ffea01021d45f"
+      "Hash": "a56a6365b3fa73293ea8d084be0d9bb0"
     },
     "Matrix": {
       "Package": "Matrix",
-      "Version": "1.2-18",
+      "Version": "1.6-1.1",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
         "R",
+        "grDevices",
         "graphics",
         "grid",
         "lattice",
@@ -38,7 +39,7 @@
         "stats",
         "utils"
       ],
-      "Hash": "5eefb6f276ded77fad5fff409e4f0557"
+      "Hash": "4767d80e0dcdf2ee70bd5ed1e1d37c4a"
     },
     "R6": {
       "Package": "R6",
@@ -107,7 +108,7 @@
     },
     "boot": {
       "Package": "boot",
-      "Version": "1.3-25",
+      "Version": "1.3-28.1",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
@@ -115,7 +116,7 @@
         "graphics",
         "stats"
       ],
-      "Hash": "bd51734a754b6c2baf28b2d1ebc11e91"
+      "Hash": "9a052fbcbe97a98ceb18dbfd30ebd96e"
     },
     "cli": {
       "Package": "cli",
@@ -140,13 +141,13 @@
     },
     "codetools": {
       "Package": "codetools",
-      "Version": "0.2-16",
+      "Version": "0.2-19",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
         "R"
       ],
-      "Hash": "89cf4b8207269ccf82fbeb6473fd662b"
+      "Hash": "c089a619a7fae175d149d89164f8c7d8"
     },
     "colorspace": {
       "Package": "colorspace",
@@ -405,7 +406,7 @@
     },
     "lattice": {
       "Package": "lattice",
-      "Version": "0.20-41",
+      "Version": "0.21-9",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
@@ -416,7 +417,7 @@
         "stats",
         "utils"
       ],
-      "Hash": "fbd9285028b0263d76d18c95ae51a53d"
+      "Hash": "5558c61e0136e247252f5f952cdaad6a"
     },
     "lifecycle": {
       "Package": "lifecycle",
@@ -443,7 +444,7 @@
     },
     "mgcv": {
       "Package": "mgcv",
-      "Version": "1.8-33",
+      "Version": "1.9-0",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
@@ -456,7 +457,7 @@
         "stats",
         "utils"
       ],
-      "Hash": "4794fe2a7f2a8190c52e87658b8a3d79"
+      "Hash": "f55d28806bb2238d6aaceb58d82edaaa"
     },
     "moments": {
       "Package": "moments",
@@ -478,7 +479,7 @@
     },
     "nlme": {
       "Package": "nlme",
-      "Version": "3.1-149",
+      "Version": "3.1-163",
       "Source": "Repository",
       "Repository": "CRAN",
       "Requirements": [
@@ -488,7 +489,7 @@
         "stats",
         "utils"
       ],
-      "Hash": "7c24ab3a1e3afe50388eb2d893aab255"
+      "Hash": "8d1938040a05566f4f7a14af4feadd6b"
     },
     "openxlsx": {
       "Package": "openxlsx",

scripts/histograms.R

@@ -15,6 +15,7 @@ library(moments)
 library(confintr)
 library(stringr)
 library(dplyr)
+library(tibble)
 font_import()
 loadfonts(device = "win")
 
@@ -25,99 +26,139 @@ rm(list = ls())
 set.seed(127)
 
 # create elasticity, policy and income group names
-elasticities = c("esubd", "esubm", "esubva")
+params = c(
+  #"esubd",
+  #"esubm",
+  #"esubva",
+  "CO2factor"
+)
 policies = c("policy", "cbam")
 inc_groups = c("lo", "mi", "hi")
 
-for(i in 1:length(elasticities)) {
+for(i in 1:length(params)) {
   # load welfare effects data frame
-  welf_name = paste("welf", elasticities[i], sep = "_")
-  
+  welf_name = paste("welf", params[i], sep = "_")
   assign(
     x = welf_name,
     value = readRDS(
-         paste0("prepared/", welf_name, ".rds")
+      paste0("prepared/", welf_name, ".rds")
     )
   )
+
+  # create subdirectory for figures
+  dir.create(
+    file.path("figures", params[i]),
+    showWarnings = F
+  )
   
-  # compute bootstrap 95% confidence intervals
-  # and create histograms
+  # create plots
   for(j in 1:length(policies)) {
     for(k in 1:length(inc_groups)) {
+      # get welfare effect from dataframe
       welf_effect_name = paste(
         policies[j], inc_groups[k], sep = "_"
       )
       welf_effect = (get(welf_name))[welf_effect_name]
       
-      # compute bootstrap 95% confidence intervals
-      ci_value = ci_mean(
-        welf_effect[, ],
-        type = "bootstrap",
-        boot_type = "perc",
-        R = 1000
-      )
+      # if parameter is CO2factor: line plot, else histogram
+      plot_type = ifelse(params[i] == "CO2factor", "plot", "hist")
       
-      # create histogram with mean and 95% confidence intervals
-      plot_name = paste("hist",
-                        elasticities[i],
-                        welf_effect_name,
-                        sep = "_"
-                        )
-      plot_value = ggplot(
-        data = get(welf_name),
-        aes(x = welf_effect[, ])) +
-        # histogram
-        geom_histogram(
-          aes(y = after_stat(density)),
-          color = "black",
-          fill = "#AFC4DE",
-          bins = 15
+      # create plot name
+      plot_name = paste(
+        plot_type,
+        params[i],
+        welf_effect_name,
+        sep = "_"
+      )
+
+      # if parameter is CO2factor: create line plot
+      if(params[i] == "CO2factor") {
+        plot_value = ggplot(
+          data = get(welf_name)
+        ) +
+          # line plot
+          geom_line(
+            aes(
+              x = CO2factor,
+              y = welf_effect[, ],
+              group = 1
+            ),
+            color = "#36638B",
+            linewidth = 1
           ) +
-        # mean
-        geom_vline(
-          xintercept = mean(welf_effect[, ]),
-          color = "#36638B",
-          linetype = "dashed",
-          size = 1
+          scale_x_reverse(limits = c(NA, 0.8)) +
+          # label axes
+          labs(
+            x = "CO2 target",
+            y = "Welfare change / %"
           ) +
-        # lower bound of 95% confidence intervals
-        geom_vline(
-          xintercept = ci_value$interval[1],
-          color = "black",
-          linetype = "dashed",
-          size = 1
+          theme_minimal(base_size = 18.5)
+      } else { # else: create histogram with mean and 95% confidence intervals
+        # compute bootstrap 95% confidence intervals
+        ci_value = ci_mean(
+          welf_effect[, ],
+          type = "bootstrap",
+          boot_type = "perc",
+          R = 1000
+        )
+
+        plot_value = ggplot(
+          data = get(welf_name),
+          aes(x = welf_effect[, ])
+        ) +
+          # histogram
+          geom_histogram(
+            aes(y = after_stat(density)),
+            color = "black",
+            fill = "#AFC4DE",
+            bins = 15
+            ) +
+          # mean
+          geom_vline(
+            xintercept = mean(welf_effect[, ]),
+            color = "#36638B",
+            linetype = "dashed",
+            linewidth = 1
           ) +
-        # lower bound of 95% confidence intervals
-        geom_vline(
-          xintercept = ci_value$interval[2],
-          color = "black",
-          linetype = "dashed",
-          size = 1
+          # lower bound of 95% confidence intervals
+          geom_vline(
+            xintercept = ci_value$interval[1],
+            color = "black",
+            linetype = "dashed",
+            linewidth = 1
           ) +
-        # label axes
-        labs(
-          x = "Welfare change / %",
-          y = "Density"
+          # lower bound of 95% confidence intervals
+          geom_vline(
+            xintercept = ci_value$interval[2],
+            color = "black",
+            linetype = "dashed",
+            linewidth = 1
           ) +
-        # Kernel density estimation
-        geom_density(
-          alpha = 0.2,
-          color = "#36638B",
-          size = 1
+          # label axes
+          labs(
+            x = "Welfare change / %",
+            y = "Density"
           ) +
-        theme_minimal(base_size = 18.5)
-      
+          # Kernel density estimation
+          geom_density(
+            alpha = 0.2,
+            color = "#36638B",
+            linewidth = 1
+          ) +
+          theme_minimal(base_size = 18.5)
+      }
+
       # save histogram as PDF
       ggsave(
-        paste(
+        file.path(
           "figures",
-          elasticities[i],
-          paste0(plot_name, ".pdf"),
-          sep = "/")
+          params[i],
+          paste0(plot_name, ".pdf")
+        )
       )
     }
   }
 }
 rm(i, j, k, welf_effect_name, welf_effect, welf_name,
-   ci_value, plot_name, plot_value, elasticities,
+   ci_value, plot_name, plot_value, params,
    inc_groups, policies)
\ No newline at end of file

scripts/welf_data.R

@@ -125,6 +125,10 @@ welf_output = foreach(
       "/",
       simplify = T
     )[, 1]
+    welf = welf %>%
+      tibble::rownames_to_column(var = "CO2factor") %>%
+      arrange(desc(CO2factor))
+    welf$CO2factor = as.numeric(welf$CO2factor)
   }
   
   # name welfare effects data frame

setup.R

@@ -20,7 +20,7 @@ renv::init(bare = TRUE)
 # Install the packages
 install.packages(c(
   "ggplot2", "readr", "extrafont", "openxlsx", "Rcpp", "tictoc", "moments",
-  "confintr", "dplyr", "stringr", "foreach", "doSNOW", "parallel"
+  "confintr", "dplyr", "stringr", "foreach", "doSNOW", "parallel", "tibble"
   ))
 
 # Take a snapshot of the renv