#####################################################################################################
##
##                                    Causal Machine Learning
##
## E. Flachaire, 16/12/2021
##
#####################################################################################################


###### Application 1: Do poor countries catch up rich countries?

library(hdm)
data("GrowthData") # the 2nd column is a vector of one #
y=as.matrix(GrowthData)[,1, drop=F]
d=as.matrix(GrowthData)[,3, drop=F]
X=as.matrix(GrowthData)[,-c(1,2,3), drop=F]

LS.fit = lm(y~d+X)
PO.fit=rlassoEffect(X, y, d, method="partialling out")
DS.fit=rlassoEffect(X, y, d, method="double selection")

LS=summary(LS.fit)$coefficients[2,]
PO=summary(PO.fit)$coefficients[1,]
DS=summary(DS.fit)$coefficients[1,]
rbind(ols=LS,double.selection=DS, partialling.out=PO)

rlasso(y~d+X, post=TRUE)$coefficients[2]


###### Application 2: The effect of gender on wage    @@@@@@@@

library(hdm)
data(cps2012)
y <- cps2012$lnw
X <- model.matrix(~-1+female+female:(widowed+divorced+separated+nevermarried+hsd08+hsd911+hsg+cg+ad+mw+so+we+exp1+exp2+exp3)+(widowed+divorced+separated+nevermarried+hsd08+hsd911+hsg+cg+ad+mw+so+we+exp1+exp2+exp3)^2,data = cps2012)
X<-X[,which(apply(X,2,var)!=0)] #exclude  constant variables 
index.gender <- grep("female", colnames(X)) 
effects.female<-rlassoEffects(x=X,y=y,index=index.gender) 
summary(effects.female)

# Generic approach to generate all covariates
Xcol=colnames(cps2012)[4:18]
dcol=colnames(cps2012)[3]
Xvar=paste(Xcol, collapse = "+")
Xint=paste("(",paste(Xcol,collapse="+"),")^2",sep="")
fmla=paste("~-1+",dcol,"+",dcol,":(",Xvar,")+",Xint,sep="")
X<-model.matrix(as.formula(fmla), data=cps2012)


##### Application 3: Effect of unemployment insurance bonus on unemployment duration (DML)

library(DoubleML)
library(mlr3)
# Initialization of the Data-Backend
data=fetch_bonus(return_type="data.table")
y="inuidur1"
d="tg"
x=c("female","black","othrace","dep1","dep2","q2","q3","q4","q5","q6","agelt35","agegt54","durable","lusd","husd")
dml_data=DoubleMLData$new(data,y_col=y,d_cols=d,x_cols=x)
##  Initialization of the Causal PLR Model
set.seed(31415) #required to replicate sample split
learner_g=lrn("regr.ranger",num.trees=500,min.node.size=2,max.depth=5)
learner_m=lrn("regr.ranger",num.trees=500,min.node.size=2,max.depth=5)
dml_plr=DoubleMLPLR$new(dml_data,
                        ml_m = learner_m,
                        ml_g = learner_g,
                        score = "partialling out",
                        n_folds = 5, n_rep = 1)
dml_plr
dml_plr$fit()
dml_plr$summary()

## Initialization of the Causal IRM Model
# Classifier for propensity score
learner_classif_m = lrn("classif.ranger", num.trees = 500,min.node.size = 2, max.depth = 5)
dml_irm=DoubleMLIRM$new(dml_data,
                        ml_m = learner_classif_m,
                        ml_g = learner_g,
                        score = "ATE", # or "ATTE"
                        n_folds = 5, n_rep = 1)
dml_irm
dml_irm$fit()
dml_irm$summary()


#####