Skip to content
Snippets Groups Projects
Commit e12d6977 authored by marius-lucas braun's avatar marius-lucas braun
Browse files

refactored preprocessing of welfare effects data, added comments

parent b5bb5739
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 152 additions and 191 deletions
Plots_Workspace.RData deleted 100644 → 0
+ 0
− 0
View file @ b5bb5739
File deleted
prepared/welf_esubd.rds 0 → 100644
+ 0
− 0
View file @ e12d6977
File added
prepared/welf_esubd_10p.rds deleted 100644 → 0
+ 0
− 0
View file @ b5bb5739
File deleted
prepared/welf_esubm.rds 0 → 100644
+ 0
− 0
View file @ e12d6977
File added
prepared/welf_esubm_10p.rds deleted 100644 → 0
+ 0
− 0
View file @ b5bb5739
File deleted
prepared/welf_esubva.rds 0 → 100644
+ 0
− 0
View file @ e12d6977
File added
prepared/welf_esubva_10p.rds deleted 100644 → 0
+ 0
− 0
View file @ b5bb5739
File deleted
scripts/histograms.R
+ 16
− 181
View file @ e12d6977
# This R code generates histograms of welfare effects based on sensitivity
# analyses of the Huebler et al. (2022) consumer split CGE model.
# WelWelfare effects are evaluated for 1000 random draws from +- 10 % intervals
# araround sector-specific estimates of the following parameters:
# - esubm(i): Armington elasticities
# - esubd(i): domestic import elasticities
# - esubva(j): input elasticities between production factors
# This R code generates histograms showing the distribution of welfare
# effects resulting from the distributional sensitivity analysis of the
# main results of Hübler et al. (2022). Separate histograms are generated
# for each policy scenario and income group. Bootstrap 95% confidence
# intervals are computed using the percentile method.
#
# Marius Braun, September 2022
install.packages("extrafont")
install.packages("Rcpp")
install.packages("openxlsx")
install.packages("moments")
install.packages("confintr")
# Marius Braun, April 2023
# load packages
library(ggplot2)
library(readr)
library(extrafont)
......@@ -22,69 +15,20 @@ library(Rcpp)
library(tictoc)
library(moments)
library(confintr)
library(stringr)
library(dplyr)
font_import()
loadfonts(device = "win")
# set random seed (Mersenne prime number because why not)
set.seed(127)
# clear workspace
rm(list = ls())
# esubm(i): Armington elasticities
# load parameter space
paramspace_esubm_10p = read.csv("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubm(i)\\parameter_space_esubm(i)_interval_10p.csv")
welf_esubm_10p = as.data.frame(matrix(, nrow = 1000, ncol = 6))
colnames(welf_esubm_10p) = c("policy_lo", "policy_mi", "policy_hi", "cbam_lo", "cbam_mi", "cbam_hi")
tic()
for(i in 1:nrow(paramspace_esubm_10p)){
# load output file for specific set of parameter values
# workaround to deal with trailing zeros and create correct filepaths
mAGR = toString(ifelse(paramspace_esubm_10p[i, 1] %% 1 == 0, format(round(paramspace_esubm_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 1] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 1] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 1], 2), nsmall = 2), paramspace_esubm_10p[i, 1]))))
mCOA = toString(ifelse(paramspace_esubm_10p[i, 2] %% 1 == 0, format(round(paramspace_esubm_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 2] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 2] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 2], 2), nsmall = 2), paramspace_esubm_10p[i, 2]))))
mCRU = toString(ifelse(paramspace_esubm_10p[i, 3] %% 1 == 0, format(round(paramspace_esubm_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 3] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 3] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 3], 2), nsmall = 2), paramspace_esubm_10p[i, 3]))))
mNGA = toString(ifelse(paramspace_esubm_10p[i, 4] %% 1 == 0, format(round(paramspace_esubm_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 4] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 4] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 4], 2), nsmall = 2), paramspace_esubm_10p[i, 4]))))
mPET = toString(ifelse(paramspace_esubm_10p[i, 5] %% 1 == 0, format(round(paramspace_esubm_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 5] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 5] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 5], 2), nsmall = 2), paramspace_esubm_10p[i, 5]))))
mFOO = toString(ifelse(paramspace_esubm_10p[i, 6] %% 1 == 0, format(round(paramspace_esubm_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 6] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 6] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 6], 2), nsmall = 2), paramspace_esubm_10p[i, 6]))))
mMIN = toString(ifelse(paramspace_esubm_10p[i, 7] %% 1 == 0, format(round(paramspace_esubm_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 7] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 7] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 7], 2), nsmall = 2), paramspace_esubm_10p[i, 7]))))
mPAP = toString(ifelse(paramspace_esubm_10p[i, 8] %% 1 == 0, format(round(paramspace_esubm_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 8] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 8] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 8], 2), nsmall = 2), paramspace_esubm_10p[i, 8]))))
mCHE = toString(ifelse(paramspace_esubm_10p[i, 9] %% 1 == 0, format(round(paramspace_esubm_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 9] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 9] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 9], 2), nsmall = 2), paramspace_esubm_10p[i, 9]))))
mNMM = toString(ifelse(paramspace_esubm_10p[i, 10] %% 1 == 0, format(round(paramspace_esubm_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 10] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 10] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 10], 2), nsmall = 2), paramspace_esubm_10p[i, 10]))))
mIRS = toString(ifelse(paramspace_esubm_10p[i, 11] %% 1 == 0, format(round(paramspace_esubm_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 11] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 11] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 11], 2), nsmall = 2), paramspace_esubm_10p[i, 11]))))
mNFM = toString(ifelse(paramspace_esubm_10p[i, 12] %% 1 == 0, format(round(paramspace_esubm_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 12] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 12] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 12], 2), nsmall = 2), paramspace_esubm_10p[i, 12]))))
mMAN = toString(ifelse(paramspace_esubm_10p[i, 13] %% 1 == 0, format(round(paramspace_esubm_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 13] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 13] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 13], 2), nsmall = 2), paramspace_esubm_10p[i, 13]))))
mELE = toString(ifelse(paramspace_esubm_10p[i, 14] %% 1 == 0, format(round(paramspace_esubm_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 14] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 14] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 14], 2), nsmall = 2), paramspace_esubm_10p[i, 14]))))
mTRN = toString(ifelse(paramspace_esubm_10p[i, 15] %% 1 == 0, format(round(paramspace_esubm_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 15] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 15] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 15], 2), nsmall = 2), paramspace_esubm_10p[i, 15]))))
mCON = toString(ifelse(paramspace_esubm_10p[i, 16] %% 1 == 0, format(round(paramspace_esubm_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 16] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 16] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 16], 2), nsmall = 2), paramspace_esubm_10p[i, 16]))))
mSER = toString(ifelse(paramspace_esubm_10p[i, 17] %% 1 == 0, format(round(paramspace_esubm_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 17] %% 0.1 == 0, format(round(paramspace_esubm_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubm_10p[i, 17] %% 0.01 == 0, format(round(paramspace_esubm_10p[i, 17], 2), nsmall = 2), paramspace_esubm_10p[i, 17]))))
f = paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubm(i)\\mAGR=", mAGR, "\\", "mCOA=", mCOA, "\\", "mCRU=", mCRU, "\\", "mNGA=", mNGA, "\\", "mPET=", mPET, "\\", "mFOO=", mFOO, "\\", "mMIN=", mMIN, "\\", "mPAP=", mPAP, "\\", "mCHE=", mCHE, "\\", "mNMM=", mNMM, "\\", "mIRS=", mIRS, "\\", "mNFM=", mNFM, "\\", "mMAN=", mMAN, "\\", "mELE=", mELE, "\\", "mTRN=", mTRN, "\\", "mCON=", mCON, "\\", "mSER=", mSER, "\\output.xlsx")
if(file.exists(f)){
welfare = read.xlsx(paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubm(i)\\mAGR=", mAGR, "\\", "mCOA=", mCOA, "\\", "mCRU=", mCRU, "\\", "mNGA=", mNGA, "\\", "mPET=", mPET, "\\", "mFOO=", mFOO, "\\", "mMIN=", mMIN, "\\", "mPAP=", mPAP, "\\", "mCHE=", mCHE, "\\", "mNMM=", mNMM, "\\", "mIRS=", mIRS, "\\", "mNFM=", mNFM, "\\", "mMAN=", mMAN, "\\", "mELE=", mELE, "\\", "mTRN=", mTRN, "\\", "mCON=", mCON, "\\", "mSER=", mSER, "\\output.xlsx"), 135)
# get welfare effects
welf_esubm_10p[i, 1] = welfare[10, 4]
welf_esubm_10p[i, 2] = welfare[19, 4]
welf_esubm_10p[i, 3] = welfare[28, 4]
welf_esubm_10p[i, 4] = welfare[11, 4]
welf_esubm_10p[i, 5] = welfare[20, 4]
welf_esubm_10p[i, 6] = welfare[29, 4]
}
# remove variables to keep workspace clean
remove(mAGR, mCOA, mCRU, mNGA, mPET, mFOO, mMIN, mPAP, mCHE, mNMM, mIRS, mNFM, mMAN, mELE, mTRN, mCON, mSER, welfare, f)
}
remove(i)
toc()
# welfare effects for random draws from +- 10 percent interval
# welf_esubm_10p = read.csv("C:\\Users\\Marius Braun\\Sensitivity analyses\\Plots\\welf_esubm(i)_10p.csv")
welf_esubm_10p$policy_lo = as.numeric(welf_esubm_10p$policy_lo) * 100
welf_esubm_10p$policy_mi = as.numeric(welf_esubm_10p$policy_mi) * 100
welf_esubm_10p$policy_hi = as.numeric(welf_esubm_10p$policy_hi) * 100
welf_esubm_10p$cbam_lo = as.numeric(welf_esubm_10p$cbam_lo) * 100
welf_esubm_10p$cbam_mi = as.numeric(welf_esubm_10p$cbam_mi) * 100
welf_esubm_10p$cbam_hi = as.numeric(welf_esubm_10p$cbam_hi) * 100
# set random seed
set.seed(127)
#### 1. esubm(i) ####
# consumer split: plot welfare effects of EU climate policy and CBAM separately for income groups
# EU climate policy
......@@ -233,63 +177,7 @@ remove(ci_esubm_cbam_hi, cbam_hi_mean_esubm_minus_error_margin, cbam_hi_mean_esu
# esubd(i): domestic import elasticities
# load parameter space
paramspace_esubd_10p = read.csv("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubd(i)\\parameter_space_esubd(i)_interval_10p.csv")
welf_esubd_10p = as.data.frame(matrix(, nrow = 1000, ncol = 6))
colnames(welf_esubd_10p) = c("policy_lo", "policy_mi", "policy_hi", "cbam_lo", "cbam_mi", "cbam_hi")
tic()
for(i in 1:nrow(paramspace_esubd_10p)){
# load output file for specific set of parameter values
# workaround to deal with trailing zeros and create correct filepaths
dAGR = toString(ifelse(paramspace_esubd_10p[i, 1] %% 1 == 0, format(round(paramspace_esubd_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 1] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 1] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 1], 2), nsmall = 2), paramspace_esubd_10p[i, 1]))))
dCOA = toString(ifelse(paramspace_esubd_10p[i, 2] %% 1 == 0, format(round(paramspace_esubd_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 2] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 2] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 2], 2), nsmall = 2), paramspace_esubd_10p[i, 2]))))
dCRU = toString(ifelse(paramspace_esubd_10p[i, 3] %% 1 == 0, format(round(paramspace_esubd_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 3] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 3] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 3], 2), nsmall = 2), paramspace_esubd_10p[i, 3]))))
dNGA = toString(ifelse(paramspace_esubd_10p[i, 4] %% 1 == 0, format(round(paramspace_esubd_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 4] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 4] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 4], 2), nsmall = 2), paramspace_esubd_10p[i, 4]))))
dPET = toString(ifelse(paramspace_esubd_10p[i, 5] %% 1 == 0, format(round(paramspace_esubd_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 5] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 5] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 5], 2), nsmall = 2), paramspace_esubd_10p[i, 5]))))
dFOO = toString(ifelse(paramspace_esubd_10p[i, 6] %% 1 == 0, format(round(paramspace_esubd_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 6] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 6] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 6], 2), nsmall = 2), paramspace_esubd_10p[i, 6]))))
dMIN = toString(ifelse(paramspace_esubd_10p[i, 7] %% 1 == 0, format(round(paramspace_esubd_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 7] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 7] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 7], 2), nsmall = 2), paramspace_esubd_10p[i, 7]))))
dPAP = toString(ifelse(paramspace_esubd_10p[i, 8] %% 1 == 0, format(round(paramspace_esubd_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 8] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 8] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 8], 2), nsmall = 2), paramspace_esubd_10p[i, 8]))))
dCHE = toString(ifelse(paramspace_esubd_10p[i, 9] %% 1 == 0, format(round(paramspace_esubd_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 9] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 9] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 9], 2), nsmall = 2), paramspace_esubd_10p[i, 9]))))
dNMM = toString(ifelse(paramspace_esubd_10p[i, 10] %% 1 == 0, format(round(paramspace_esubd_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 10] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 10] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 10], 2), nsmall = 2), paramspace_esubd_10p[i, 10]))))
dIRS = toString(ifelse(paramspace_esubd_10p[i, 11] %% 1 == 0, format(round(paramspace_esubd_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 11] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 11] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 11], 2), nsmall = 2), paramspace_esubd_10p[i, 11]))))
dNFM = toString(ifelse(paramspace_esubd_10p[i, 12] %% 1 == 0, format(round(paramspace_esubd_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 12] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 12] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 12], 2), nsmall = 2), paramspace_esubd_10p[i, 12]))))
dMAN = toString(ifelse(paramspace_esubd_10p[i, 13] %% 1 == 0, format(round(paramspace_esubd_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 13] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 13] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 13], 2), nsmall = 2), paramspace_esubd_10p[i, 13]))))
dELE = toString(ifelse(paramspace_esubd_10p[i, 14] %% 1 == 0, format(round(paramspace_esubd_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 14] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 14] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 14], 2), nsmall = 2), paramspace_esubd_10p[i, 14]))))
dTRN = toString(ifelse(paramspace_esubd_10p[i, 15] %% 1 == 0, format(round(paramspace_esubd_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 15] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 15] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 15], 2), nsmall = 2), paramspace_esubd_10p[i, 15]))))
dCON = toString(ifelse(paramspace_esubd_10p[i, 16] %% 1 == 0, format(round(paramspace_esubd_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 16] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 16] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 16], 2), nsmall = 2), paramspace_esubd_10p[i, 16]))))
dSER = toString(ifelse(paramspace_esubd_10p[i, 17] %% 1 == 0, format(round(paramspace_esubd_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 17] %% 0.1 == 0, format(round(paramspace_esubd_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubd_10p[i, 17] %% 0.01 == 0, format(round(paramspace_esubd_10p[i, 17], 2), nsmall = 2), paramspace_esubd_10p[i, 17]))))
f = paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubd(i)\\dAGR=", dAGR, "\\", "dCOA=", dCOA, "\\", "dCRU=", dCRU, "\\", "dNGA=", dNGA, "\\", "dPET=", dPET, "\\", "dFOO=", dFOO, "\\", "dMIN=", dMIN, "\\", "dPAP=", dPAP, "\\", "dCHE=", dCHE, "\\", "dNMM=", dNMM, "\\", "dIRS=", dIRS, "\\", "dNFM=", dNFM, "\\", "dMAN=", dMAN, "\\", "dELE=", dELE, "\\", "dTRN=", dTRN, "\\", "dCON=", dCON, "\\", "dSER=", dSER, "\\output.xlsx")
if(file.exists(f)){
welfare = read.xlsx(paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubd(i)\\dAGR=", dAGR, "\\", "dCOA=", dCOA, "\\", "dCRU=", dCRU, "\\", "dNGA=", dNGA, "\\", "dPET=", dPET, "\\", "dFOO=", dFOO, "\\", "dMIN=", dMIN, "\\", "dPAP=", dPAP, "\\", "dCHE=", dCHE, "\\", "dNMM=", dNMM, "\\", "dIRS=", dIRS, "\\", "dNFM=", dNFM, "\\", "dMAN=", dMAN, "\\", "dELE=", dELE, "\\", "dTRN=", dTRN, "\\", "dCON=", dCON, "\\", "dSER=", dSER, "\\output.xlsx"), 135)
# get welfare effects
welf_esubd_10p[i, 1] = welfare[10, 4]
welf_esubd_10p[i, 2] = welfare[19, 4]
welf_esubd_10p[i, 3] = welfare[28, 4]
welf_esubd_10p[i, 4] = welfare[11, 4]
welf_esubd_10p[i, 5] = welfare[20, 4]
welf_esubd_10p[i, 6] = welfare[29, 4]
}
# remove variables to keep workspace clean
remove(dAGR, dCOA, dCRU, dNGA, dPET, dFOO, dMIN, dPAP, dCHE, dNMM, dIRS, dNFM, dMAN, dELE, dTRN, dCON, dSER, welfare, f)
}
remove(i)
toc()
# welfare effects for random draws from +- 10 percent interval
# welf_esubd_10p = read.csv("C:\\Users\\Marius Braun\\Sensitivity analyses\\Plots\\welf_esubd(i)_10p.csv")
welf_esubd_10p$policy_lo = as.numeric(welf_esubd_10p$policy_lo) * 100
welf_esubd_10p$policy_mi = as.numeric(welf_esubd_10p$policy_mi) * 100
welf_esubd_10p$policy_hi = as.numeric(welf_esubd_10p$policy_hi) * 100
welf_esubd_10p$cbam_lo = as.numeric(welf_esubd_10p$cbam_lo) * 100
welf_esubd_10p$cbam_mi = as.numeric(welf_esubd_10p$cbam_mi) * 100
welf_esubd_10p$cbam_hi = as.numeric(welf_esubd_10p$cbam_hi) * 100
#### 2. esubd(i) ####
# consumer split: plot welfare effects of EU climate policy and CBAM separately for income groups
# EU climate policy
......@@ -439,60 +327,7 @@ remove(ci_esubd_cbam_hi, cbam_hi_mean_esubd_minus_error_margin, cbam_hi_mean_esu
# esubva(j): input elasticities between production factors
# load parameter space
paramspace_esubva_10p = read.csv("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubva(j)\\parameter_space_esubva(j)_interval_10p.csv")
welf_esubva_10p = as.data.frame(matrix(, nrow = 100, ncol = 6))
colnames(welf_esubva_10p) = c("policy_lo", "policy_mi", "policy_hi", "cbam_lo", "cbam_mi", "cbam_hi")
tic()
for(i in 1:nrow(paramspace_esubva_10p)){
# load output file for specific set of parameter values
# workaround to deal with trailing zeros and create correct filepaths
vAGR = toString(ifelse(paramspace_esubva_10p[i, 1] %% 1 == 0, format(round(paramspace_esubva_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 1] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 1], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 1] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 1], 2), nsmall = 2), paramspace_esubva_10p[i, 1]))))
vCOA = toString(ifelse(paramspace_esubva_10p[i, 2] %% 1 == 0, format(round(paramspace_esubva_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 2] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 2], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 2] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 2], 2), nsmall = 2), paramspace_esubva_10p[i, 2]))))
vCRU = toString(ifelse(paramspace_esubva_10p[i, 3] %% 1 == 0, format(round(paramspace_esubva_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 3] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 3], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 3] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 3], 2), nsmall = 2), paramspace_esubva_10p[i, 3]))))
vNGA = toString(ifelse(paramspace_esubva_10p[i, 4] %% 1 == 0, format(round(paramspace_esubva_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 4] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 4], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 4] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 4], 2), nsmall = 2), paramspace_esubva_10p[i, 4]))))
vPET = toString(ifelse(paramspace_esubva_10p[i, 5] %% 1 == 0, format(round(paramspace_esubva_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 5] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 5], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 5] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 5], 2), nsmall = 2), paramspace_esubva_10p[i, 5]))))
vFOO = toString(ifelse(paramspace_esubva_10p[i, 6] %% 1 == 0, format(round(paramspace_esubva_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 6] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 6], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 6] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 6], 2), nsmall = 2), paramspace_esubva_10p[i, 6]))))
vMIN = toString(ifelse(paramspace_esubva_10p[i, 7] %% 1 == 0, format(round(paramspace_esubva_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 7] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 7], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 7] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 7], 2), nsmall = 2), paramspace_esubva_10p[i, 7]))))
vPAP = toString(ifelse(paramspace_esubva_10p[i, 8] %% 1 == 0, format(round(paramspace_esubva_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 8] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 8], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 8] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 8], 2), nsmall = 2), paramspace_esubva_10p[i, 8]))))
vCHE = toString(ifelse(paramspace_esubva_10p[i, 9] %% 1 == 0, format(round(paramspace_esubva_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 9] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 9], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 9] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 9], 2), nsmall = 2), paramspace_esubva_10p[i, 9]))))
vNMM = toString(ifelse(paramspace_esubva_10p[i, 10] %% 1 == 0, format(round(paramspace_esubva_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 10] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 10], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 10] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 10], 2), nsmall = 2), paramspace_esubva_10p[i, 10]))))
vIRS = toString(ifelse(paramspace_esubva_10p[i, 11] %% 1 == 0, format(round(paramspace_esubva_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 11] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 11], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 11] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 11], 2), nsmall = 2), paramspace_esubva_10p[i, 11]))))
vNFM = toString(ifelse(paramspace_esubva_10p[i, 12] %% 1 == 0, format(round(paramspace_esubva_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 12] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 12], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 12] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 12], 2), nsmall = 2), paramspace_esubva_10p[i, 12]))))
vMAN = toString(ifelse(paramspace_esubva_10p[i, 13] %% 1 == 0, format(round(paramspace_esubva_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 13] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 13], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 13] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 13], 2), nsmall = 2), paramspace_esubva_10p[i, 13]))))
vELE = toString(ifelse(paramspace_esubva_10p[i, 14] %% 1 == 0, format(round(paramspace_esubva_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 14] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 14], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 14] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 14], 2), nsmall = 2), paramspace_esubva_10p[i, 14]))))
vTRN = toString(ifelse(paramspace_esubva_10p[i, 15] %% 1 == 0, format(round(paramspace_esubva_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 15] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 15], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 15] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 15], 2), nsmall = 2), paramspace_esubva_10p[i, 15]))))
vCON = toString(ifelse(paramspace_esubva_10p[i, 16] %% 1 == 0, format(round(paramspace_esubva_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 16] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 16], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 16] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 16], 2), nsmall = 2), paramspace_esubva_10p[i, 16]))))
vSER = toString(ifelse(paramspace_esubva_10p[i, 17] %% 1 == 0, format(round(paramspace_esubva_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 17] %% 0.1 == 0, format(round(paramspace_esubva_10p[i, 17], 1), nsmall = 1), ifelse(paramspace_esubva_10p[i, 17] %% 0.01 == 0, format(round(paramspace_esubva_10p[i, 17], 2), nsmall = 2), paramspace_esubva_10p[i, 17]))))
f = paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubva(j)\\vAGR=", vAGR, "\\", "vCOA=", vCOA, "\\", "vCRU=", vCRU, "\\", "vNGA=", vNGA, "\\", "vPET=", vPET, "\\", "vFOO=", vFOO, "\\", "vMIN=", vMIN, "\\", "vPAP=", vPAP, "\\", "vCHE=", vCHE, "\\", "vNMM=", vNMM, "\\", "vIRS=", vIRS, "\\", "vNFM=", vNFM, "\\", "vMAN=", vMAN, "\\", "vELE=", vELE, "\\", "vTRN=", vTRN, "\\", "vCON=", vCON, "\\", "vSER=", vSER, "\\output.xlsx")
if(file.exists(f)){
welfare = read.xlsx(paste0("C:\\Users\\Marius Braun\\Sensitivity analyses\\esubva(j)\\vAGR=", vAGR, "\\", "vCOA=", vCOA, "\\", "vCRU=", vCRU, "\\", "vNGA=", vNGA, "\\", "vPET=", vPET, "\\", "vFOO=", vFOO, "\\", "vMIN=", vMIN, "\\", "vPAP=", vPAP, "\\", "vCHE=", vCHE, "\\", "vNMM=", vNMM, "\\", "vIRS=", vIRS, "\\", "vNFM=", vNFM, "\\", "vMAN=", vMAN, "\\", "vELE=", vELE, "\\", "vTRN=", vTRN, "\\", "vCON=", vCON, "\\", "vSER=", vSER, "\\output.xlsx"), 135)
# get welfare effects
welf_esubva_10p[i, 1] = welfare[10, 4]
welf_esubva_10p[i, 2] = welfare[19, 4]
welf_esubva_10p[i, 3] = welfare[28, 4]
welf_esubva_10p[i, 4] = welfare[11, 4]
welf_esubva_10p[i, 5] = welfare[20, 4]
welf_esubva_10p[i, 6] = welfare[28, 4]
}
# remove variables to keep workspace clean
remove(vAGR, vCOA, vCRU, vNGA, vPET, vFOO, vMIN, vPAP, vCHE, vNMM, vIRS, vNFM, vMAN, vELE, vTRN, vCON, vSER, welfare, f)
}
remove(i)
toc()
# welfare effects for random draws from +- 10 percent interval
welf_esubva_10p$policy_lo = as.numeric(welf_esubva_10p$policy_lo) * 100
welf_esubva_10p$policy_mi = as.numeric(welf_esubva_10p$policy_mi) * 100
welf_esubva_10p$policy_hi = as.numeric(welf_esubva_10p$policy_hi) * 100
welf_esubva_10p$cbam_lo = as.numeric(welf_esubva_10p$cbam_lo) * 100
welf_esubva_10p$cbam_mi = as.numeric(welf_esubva_10p$cbam_mi) * 100
welf_esubva_10p$cbam_hi = as.numeric(welf_esubva_10p$cbam_hi) * 100
#### 3. esubva(j) ####
# consumer split: plot welfare effects of EU climate policy and CBAM separately for income groups
# EU climate policy
......
scripts/paramspaces.R
+ 28
− 9
View file @ e12d6977
# This R code generates random draws from a +-10% interval around the
# elasticities of substitution esubd(i) for each sector in the Pothen and
# Huebler model (2018). It then writes these draws to a CSV file to be
# used as the parameter space for sensitivity analysis in Snakemake.
# This R code generates parameter spaces for the distributional
# sensivitity analyses of the main results of Hübler et al. (2022).
# More specifically, 1000 random draws from a +-10% interval around
# each of the sector-specific elasiticity estimates are generated.
# This results in 1000 sets of sectoral parameter values to be used
# in a parameter sweep in Snakemake.
#
# A separate parameter space is generated for each of the three sets
# of sector-specific elasticity parameters present in the model:
# - esubd(i): elasticities between domestically produced versus
# imported goods
# - esubm(i): Armington elasticities
# - esubva(j): elasticities between production factors
#
# Marius Braun, May 2022
# Marius Braun, April 2023
# load packages
library(ggplot2)
library(readr)
library(extrafont)
library(openxlsx)
library(Rcpp)
library(tictoc)
library(moments)
library(confintr)
library(stringr)
library(dplyr)
# clear workspace
rm(list = ls())
num_draws = 1000 # number of draws from interval
set.seed(127) # set random seed
variation = 0.1
num_decimals = 3
variation = 0.1 # size of variation
num_decimals = 3 # max. number of decimals (important due to file path length restrictions)
# create vector of elasticity names
elasticities = c("esubd", "esubm", "esubva")
for(i in 1:length(elasticities)) {
# load CSV file with baseline parameter values
filename = paste("intervals_10p", elasticities[i], sep = "_")
......@@ -58,5 +78,4 @@ for(i in 1:length(elasticities)) {
row.names = FALSE
)
}
rm(i, j, filename, file, paramspace)
\ No newline at end of file
rm(i, j, filename, file, paramspace)
\ No newline at end of file
scripts/welf_data.R 0 → 100644
+ 107
− 0
View file @ e12d6977
# This R code extracts welfare effects from the output files of the
# distributional sensitivity analysis of the main results of Hübler
# et al. (2022). Welfare effects are then stored in a separate data
# frame for esubd(i), esubm(i) and esubva(j) and saved in RDS files
# to be used in subsequent analyses.
#
# Marius Braun, April 2023
# load packages
library(ggplot2)
library(readr)
library(extrafont)
library(openxlsx)
library(Rcpp)
library(tictoc)
library(moments)
library(confintr)
library(stringr)
library(dplyr)
# clear workspace
rm(list = ls())
# this is where you put the output files from the sensitivity analysis
# (use separate subdirectories for esubd(i), esubm(i), and esubva(j))
dir = "C:/Users/Marius Braun/output_sensitivity"
# string vectors for elasticity, policy and income group names
elasticities = c("esubd", "esubm", "esubva")
policies = c("policy", "cbam")
inc_groups = c("lo", "mi", "hi")
for(i in 1:length(elasticities)) {
# load file paths of output files
filenames = as.data.frame(
list.files(
paste(
dir,
elasticities[i],
sep = "/"
),
pattern = "output.xlsx",
full.names = T,
recursive = T)
)
colnames(filenames) = "files"
# create empty data frame to store welfare effects
welf_name = paste(
"welf",
elasticities[i],
sep = "_"
)
welf = as.data.frame(
matrix(,
nrow = nrow(filenames),
ncol = length(policies) * length(inc_groups)
)
)
# name columns of data frame according to policy and income group
for(j in 1:length(policies)) {
for(k in 1:length(inc_groups)) {
col_num = k +
(as.numeric(policies[j] == "cbam") * length(inc_groups))
colnames(welf)[col_num] = paste(
policies[j],
inc_groups[k],
sep = "_"
)
}
}
# get welfare effects from output files
for(l in 1:nrow(filenames)) {
f = filenames$files[l]
if(file.exists(f)) {
welfare = read.xlsx(f, sheet = "welfp") %>%
filter(TOC == "DEU" &
str_detect(X3,
paste(
paste0("\\b",
policies),
collapse = "|")
)
)
welfare = welfare %>% arrange(desc(X3))
}
welf[l, ] = t(welfare$X4)
}
# name welfare effects data frame
assign(
x = welf_name,
value = welf
)
# save welfare effects data frame as RDS file
saveRDS(get(welf_name),
paste0(
"prepared/",
welf_name,
".rds")
)
}
rm(i, j, k, l, f, col_num, welf_name, welf, filenames,
welfare, dir, elasticities, inc_groups, policies)
\ No newline at end of file
setup.R
+ 1
− 1
View file @ e12d6977
......@@ -20,7 +20,7 @@ renv::init(bare = TRUE)
# Install the packages
install.packages(c(
"ggplot2", "readr", "extrafont", "openxlsx", "Rcpp", "tictoc", "moments",
"confintr", "dplyr"
"confintr", "dplyr", "stringr"
))
# Take a snapshot of the renv
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment