#############################################################################

### CLARE: a Causal Machine Learning Approach to Resilience Estimation ###    

#############################################################################     

#################################################################################################################################################################
#################################################################################################################################################################

##
##  This script: Causal Machine Learning Analysis - Estimation of causal relationships and variable importance weights ##   
##                                                            
##  Date: October 2025 ##
##                                                            
##  Object: Main analysis for the forecasting exercise (using EDDI drought data) ##
##

################################################################################################################################################################
################################################################################################################################################################

#============================================================================================================================================================#
# This R script carries out the estimation of shock effects and variable weights using the causal forest algorithm (Wager & Athey, 2018; Athey et al., 2019) #
#============================================================================================================================================================#
#========================================================================================================================================================================#
# Final outputs: 1) "results_forecasting_main.csv" (variable importance weights to be used in Stata for CLARE estimation); 2) Table B.1 & Figures B.1 and B.2 in the WP  #
#========================================================================================================================================================================#
# install.packages("randomForest")
# install.packages("grf")
# install.packages("sufrep")
library(sufrep)
library(grf)
library(randomForest)
library(dplyr)

options(grf.legacy.seed = TRUE)
set.seed(1234)
Sys.setenv(OMP_NUM_THREADS = 1, MKL_NUM_THREADS = 1, OPENBLAS_NUM_THREADS = 1)

# Causal forest analysis (Wager & Athey, 2018; Athey et al., 2019) of the impact of droughts on food insecurity: Forecasting Analysis.
dataset <- read.csv("../Data/dataset_start.csv") # please change working directory to load this file
names(dataset)
set.seed(1234) 
dataset$hhid <- as.factor(dataset$hhid)
dataset$wave <- as.factor(dataset$wave)
encoder <- make_encoder(dataset[, c(6:23)], dataset$hhid, method="means") 
X3 <- encoder(dataset[, c(6:23)], dataset$hhid)
names(X3)
encoder2 <- make_encoder(dataset[, c(6:23)], dataset$wave, method="means") 
X4 <- encoder2(dataset[, c(6:23)], dataset$wave)
X3 <- cbind(X3,X4[19:36])
Y.forest = regression_forest(X3, dataset$food_insecurity, clusters = dataset$hhid, tune.parameters = "all", num.threads = 12)
Y.hat = predict(Y.forest)$predictions
W.forest = regression_forest(X3, dataset$shock1, clusters = dataset$hhid, tune.parameters = "all", num.threads = 12)
W.hat = predict(W.forest)$predictions
# Causal forest estimation
cf = causal_forest(dataset[, c(6:23)], dataset$food_insecurity, dataset$shock1,
                   Y.hat = Y.hat,
                   W.hat = W.hat,
                   clusters = dataset$hhid,
                   tune.parameters = "all",
                   num.threads = 12)
# Prediction
dataset2 <- read.csv("../Data/dataset_all.csv") 
dataset2$hhid <- as.factor(dataset2$hhid)
dataset2$wave <- as.factor(dataset2$wave)
encoderT <- make_encoder(dataset2[, c(6:23)], dataset2$hhid, method="means") 
X6 <- encoderT(dataset2[, c(6:23)], dataset2$hhid)
encoderT2 <- make_encoder(dataset2[, c(6:23)], dataset2$wave, method="means") 
X7 <- encoderT2(dataset2[, c(6:23)], dataset2$wave)
X6 <- cbind(X6,X7[19:36])
ymatrix <- X6
xmatrix <- as.matrix(dataset2[, c(6:23)])
Y.hat2 = predict(Y.forest, newdata = ymatrix)$predictions
W.hat2 = predict(W.forest, newdata = ymatrix)$predictions
tau.hat = predict(cf, newdata = xmatrix)
mu.hat.0 <- Y.hat2 - W.hat2 * tau.hat
mu.hat.1 <- Y.hat2 + (1 - W.hat2) * tau.hat
tau.hat = tau.hat*100
hist(tau.hat$predictions, ylim= c(0,10000), breaks = c(-15, -10, -5, 0,  5,  10,  15,  20,  25), freq=TRUE) # This is Figure B.1 in the paper
test_calibration(cf) # This is Table B.1 in the WP

data <- cbind(dataset2,tau.hat, mu.hat.0, mu.hat.1,W.hat2, Y.hat2)
write.csv(data,"../Outputs/dataset_allexp_forecasting_main.csv", na = "", row.names = FALSE)
# get ATE
ATE <- average_treatment_effect(cf, target.sample = "overlap")
ATE # ATE signals a large positive and significant effects of drought on the probability of being food insecure
# getting variable importance
varimp = variable_importance(cf)
varimp
vi <- cf %>% # vi contains the variable importance measures (summing up to 1) which we will use as data-driven weights for computing CLARE
  variable_importance() %>%
  as.data.frame() %>%
  mutate(variable = colnames(cf$X.orig)) %>%
  arrange(desc(V1)) # This is Figure B.2 in the paper
# export the results
write.csv(vi,"../Outputs/results_forecasting_main.csv", na = "", row.names = FALSE)
# now we move to Stata for CLARE aggregation and evaluation