# 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/Bases Revisar")
getwd()

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

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

# Dropping the observations corresponding to firms either being a university or without employees 
innovation <- innovation %>%
  filter(!(total_empleados == 0 | id_b == 25734))

table(innovation$dfinan_publ)

innovation$dfinan_publ <- as.numeric(
  innovation$dfinan_publ == "Recibe transf directa"
)

# Parallelization 
set.seed(123)

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

# Running the model 
ps.Innova <- ps(dfinan_publ ~ 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 of 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/Bases Revisar/psc_modelo_test1.csv", 
          row.names = FALSE)