# Propensity with twang "bgm"
# First, cleaning the objects
rm(list=ls(all=TRUE))
ls()

# Fixing the directory of the workspace ...
setwd("C:/Users/juanq/OneDrive/Documentos/Paper Innovación/Biprobit")
getwd()

# Loading the database ...
innovation<-read.csv("Base_modelo_biprobit_full2.csv",  header = TRUE)

# Package installation
install.packages("twang", dependencies=T)
install.packages("survey", dependencies=T)
install.packages("doParallel", dependencies=T)
install.packages("tidyr", dependencies=T)
install.packages("mets", dependencies=T)
install.packages("GJRM", dependencies=T)

# Loading packages
library("twang")
library("survey")
library("doParallel")
library("tidyr")
library("mets")
library("GJRM")


# Parallelization
set.seed(123)

# Run model in parallel
cl <- makeCluster(detectCores())
registerDoParallel(cl)

# Running the model 
ps.Innova <- ps(dtrai ~ total_empleados + edad_empresa +
                  pge + dep_id + col_clien + col_univ + col_prov,
                data = innovation,
                n.trees=5000,
                interaction.depth=2,
                shrinkage=0.01,
                perm.test.iters=0,
                stop.method=c("es.mean","ks.max"),
                estimand = "ATE",
                sampw = innovation$fe_empresasx,
                verbose=FALSE)

##############☺
stopCluster(cl)

# Plotting
plot(ps.Innova)

# If we wish to focus on only one stopping rule, the plotting commands 
#also take a subset argument.

plot(ps.Innova, subset = 2)

####### Model summary 
summary(ps.Innova$gbm.obj,
        n.trees=ps.Innova$desc$ks.max.ATE$n.trees,
        plot=F) 

######### Balance 
lalonde.balance <- bal.table(ps.Innova)
lalonde.balance


summary(ps.Innova)

######### Propensity score creation

innovation$wm1 <- get.weights(ps.Innova, stop.method="es.mean")
design.ps <- svydesign(ids=~1, weights=~wm1, data=innovation)

### Database creation for the propensity score variable 
psc_modelo1<- innovation %>%
  dplyr::select(id_b, wm1)

write.csv(psc_modelo1, file = "C:/Users/juanq/OneDrive/Documentos/Paper Innovación/Biprobit/psc_modelo1.csv",
          row.names = FALSE)