##################################################
################ PREPARATION #####################
##################################################
# Clear environment
rm(list=setdiff(ls(), "path_to"))

# Load data
gcd = readRDS(path_to("gcd"))


##################################################
################ STANDARDIZE #####################
##################################################
# Prepare
vars = c(
  "coop_belief", "risktaking", "patience", "altruism", "norm"
)
weight_var = "wgt_expctrl"

# Standardize!
for (v in vars) {
  mu = Hmisc::wtd.mean(gcd[, v], gcd[, weight_var])
  omega = sqrt(Hmisc::wtd.var(gcd[, v], gcd[, weight_var]))
  gcd[, paste0(v, "_std")] = (gcd[, v, drop=T] - mu)/omega
  var_label(gcd[, paste0(v, "_std")]) = paste0("Standarized version of variable:", v, "; using weights: ", weight_var)
}

# Income percentile variables
dist = Hmisc::wtd.quantile(
  gcd$income_pc_log[gcd$treatment == 0],
  weights=gcd$wgt_expctrl_pop15[gcd$treatment == 0],
  probs=seq(0,1,0.01)
)
labels = as.numeric(gsub("[^0-9]", "", names(dist[!duplicated(dist)]))) + 1
labels = labels[-length(labels)]
dist = dist[!duplicated(dist)]
gcd$income_pctile_pop15 = as.numeric(as.character(cut(gcd$income_pc_log, breaks=dist, include.lowest=T, labels = labels)))
gcd$income_pctile_10 = gcd$income_pctile_pop15/10
var_label(gcd$income_pctile_pop15) = "Rank (0-100) in global income distribution (distribution of ctrl group, weighted with wgt_expctrl_pop15)"
var_label(gcd$income_pctile_10) = "Rank (0-100) in global income distribution divided by 10 for better interpretability (distribution of ctrl group, weighted with wgt_expctrl_pop15)"

# Export
saveRDS(gcd, path_to("gcd"))
