
## remove (almost) everything in the working environment.
rm(list = ls())

# to import xl files
library(readxl)
#to write excel files
require(writexl)

#to do clustered SE
if(!require(estimatr)){
    install.packages("estimatr")
    require(estimatr)}

#alternative to stargazer
if(!require(texreg)){
    install.packages("texreg")
    require(texreg)}

# Set working directory based on the this.path package, which reliably detects the path of the running script. and projects.
if (!require(this.path)) {
  install.packages("this.path")
  require(this.path)
}
setwd(dirname(this.path()))


# load persons-activities data
base <- as.data.frame(read_excel("Profs.xlsx"))
#remove 9 occupations
base<-base[!(base$Field=="9"),]

#load siblings table
profs <- as.data.frame(read_excel("Siblings.xlsx"))

#load rural/urban dimension (based on Bairoch)
rural<- as.data.frame(read_excel("rural-urban.xlsx"))

#Opening the dictionary containing notoriety metrics for each individual
basew <- as.data.frame(read_excel("Dictionary-full.xlsx"))

#Opening the files with sons of professors (updated from de la Croix, Goni, JOEG, 2024)
son1 <- as.data.frame(read_excel("Families.xlsx"))
son2 <- as.data.frame(read_excel("Families-bis.xlsx"))
list_sons <- unique(c(son1$Son,son2$Son))

#list of institutions
listi <- c("AcadHessehamburg-1775","Ucopenhagen-1475","Ugroningen-1612","Ulund-1666","Uuppsala-1477","AcadCopenhagen-1742","Urostock-1419",
"Ukiel-1652","Uabo-1640","Udorpat-1632","Aamsterdam-1632","AcadDantzig-1743","Uedinburgh-1582","Ukonigsberg-1544","AcadStockholm-1739",
"Uglasgow-1451","Ugreifswald-1456","Gdanzig-1558","AcadStpetersburg-1724","AcadLund-1778","Ufraneker-1585","AcadHaarlem-1752",
"AcadEdinburgh-1783","AcadUppsala-1728","Uaberdeennew-1593","Uaberdeenold-1495","Ustandrews-1411",
"Uleiden-1575","AcadTrondheim-1766","Uutrecht-1636","GymSoro-1586")
#select institutions
base<-base[base$IKey %in% listi,]


nprofs<- nrow(base)
ninstit <- length(listi)

for (i in 1:nprofs)
{
id<-base[i,"PKey"]

base[i,"fatherprof"]<- ifelse(id %in% list_sons,1,0)

weight=1/nrow(base[base$PKey==id,])  # nb of universities where id was appointed
base[i,"titles"] <-  asinh(as.numeric(basew[basew$PKey==id,"titles"]) )
base[i,"longevity"] <-  base[i,"YearDied"]-  base[i,"YearBorn"]
base[i,"ageentry"] <- ifelse(!is.na(pmin(base[i,]$Beginning,base[i,]$DateApprox,na.rm=T)),pmin(base[i,]$Beginning,base[i,]$DateApprox,na.rm=T),
ifelse(!is.na(base[i,]$YearBorn),base[i,]$YearBorn+30,
ifelse(!is.na(base[i,]$End),base[i,]$End,base[i,]$YearDied)))-  base[i,"YearBorn"]
base[i,"leiden"]<- ifelse("Uleiden-1575" %in% base[i,"IKey"],1,0)
base[i,"ulund"]<- ifelse("Ulund-1666" %in% base[i,"IKey"],1,0)
base[i,"copenhagen"]<- ifelse("Ucopenhagen-1475" %in% base[i,"IKey"],1,0)
base[i,"acaddnk"]<- ifelse("AcadCopenhagen-1742" %in% base[i,"IKey"],1,0)
base[i,"groningen"]<- ifelse("Ugroningen-1612" %in% base[i,"IKey"],1,0)
base[i,"uppsala"]<- ifelse("Uuppsala-1477" %in% base[i,"IKey"],1,0)
base[i,"rostock"]<- ifelse("Urostock-1419" %in% base[i,"IKey"],1,0)
base[i,"kiel"]<- ifelse("Ukiel-1652" %in% base[i,"IKey"],1,0)
base[i,"abo"]<- ifelse("Uabo-1640" %in% base[i,"IKey"],1,0)
base[i,"dorpat"]<- ifelse("Udorpat-1632" %in% base[i,"IKey"],1,0)
base[i,"amsterdam"]<- ifelse("Aamsterdam-1632" %in% base[i,"IKey"],1,0)
base[i,"Adanzig"]<- ifelse("AcadDantzig-1743" %in% base[i,"IKey"],1,0)
base[i,"edinburgh"]<- ifelse("Uedinburgh-1582" %in% base[i,"IKey"],1,0)
base[i,"konigsberg"]<- ifelse("Ukonigsberg-1544" %in% base[i,"IKey"],1,0)
base[i,"stockholm"]<- ifelse("AcadStockholm-1739" %in% base[i,"IKey"],1,0)
base[i,"glasgow"]<- ifelse("Uglasgow-1451"  %in% base[i,"IKey"],1,0)
base[i,"greifswald"]<- ifelse("Ugreifswald-1456"  %in% base[i,"IKey"],1,0)
base[i,"Gdanzig"]<- ifelse("Gdanzig-1558"  %in% base[i,"IKey"],1,0)
base[i,"petersburg"]<- ifelse("AcadStpetersburg-1724"  %in% base[i,"IKey"],1,0)
base[i,"alund"]<- ifelse("AcadLund-1778"  %in% base[i,"IKey"],1,0)
base[i,"franeker"]<- ifelse("Ufraneker-1585"  %in% base[i,"IKey"],1,0)
base[i,"haarlem"]<- ifelse("AcadHaarlem-1752"  %in% base[i,"IKey"],1,0)
base[i,"aedinburgh"]<- ifelse("AcadEdinburgh-1783"  %in% base[i,"IKey"],1,0)
base[i,"auppsala"]<- ifelse("AcadUppsala-1728" %in% base[i,"IKey"],1,0)
base[i,"aberdeen"]<- ifelse("Uaberdeennew-1593" %in% base[i,"IKey"] | "Uaberdeenold-1495" %in% base[i,"IKey"],1,0)
base[i,"andrews"]<- ifelse("Ustandrews-1411" %in% base[i,"IKey"],1,0)
base[i,"hamburg"]<- ifelse("AcadHessehamburg-1775" %in% base[i,"IKey"],1,0)
base[i,"trondheim"]<- ifelse("AcadTrondheim-1766" %in% base[i,"IKey"],1,0)
base[i,"utrecht"]<- ifelse("Uutrecht-1636" %in% base[i,"IKey"],1,0)
base[i,"soro"]<- ifelse("GymSoro-1586" %in% base[i,"IKey"],1,0)


base[i,"refdate"]<- base[i,"YearBorn"]

base[i,"geni"] <- ifelse(grepl("geni.com",base[i,"Genealogy"]),1,0)
base[i,"geneanet"] <- ifelse(grepl("geneanet",base[i,"Genealogy"]),1,0)
}

#define fields dummies
base$theo <- ifelse(base$Field==1 | (!is.na(base$Field2) & base$Field2==1) | (!is.na(base$Field3) & base$Field3==1),1,0) 
base$law <- ifelse(base$Field==2 | (!is.na(base$Field2) & base$Field2==2) | (!is.na(base$Field3) & base$Field3==2),1,0) 
base$med <- ifelse(base$Field==4 | (!is.na(base$Field2) & base$Field2==4) | (!is.na(base$Field3) & base$Field3==4),1,0) 
base$sci <- ifelse(base$Field==5 | (!is.na(base$Field2) & base$Field2==5) | (!is.na(base$Field3) & base$Field3==5),1,0) 

#normalizing sum of dummies to one
base$nfield <-base$theo +base$law + base$med + base$sci
base$theo <- base$theo/ifelse(base$nfield>1,base$nfield,1)
base$law<- base$law/ifelse(base$nfield>1,base$nfield,1)
base$med <- base$med/ifelse(base$nfield>1,base$nfield,1)
base$sci <- base$sci/ifelse(base$nfield>1,base$nfield,1)

#######################GENEALOGIES

sink("sink from build data base.txt")

#list missing genealogies
print("missing genealogies:")
base[is.na(base$Genealogy),"PKey"]


#people with genealogy but no siblings
list1=unique(base[!is.na(base$Genealogy) & nchar(base$Genealogy)>13,"PKey"])
list2=unique(profs$PKey)

#people with a genealogy but no entry in siblings
setdiff(list1,list2)

#we remove those with no sibshipsize
print(paste("We remove",nrow(profs[is.na(profs$SibshipSize) | is.na(profs$AmongWhichGirls),]),"genealogies with no sibshipsize or no gender"))
profs <- profs[!is.na(profs$SibshipSize) & !is.na(profs$AmongWhichGirls),]


#select profs with a genealogy
listprofs<-profs$PKey

print(paste("we have",length(listprofs),"genealogies with some information on sibshipsize"))

base$exist <- ifelse(base$PKey %in% listprofs,1,0)

print(paste("we have a total of",nrow(base),"occupations"))
print(paste(sum(base$exist),"of them can be linked to a genealogy"))

print(paste("this involves",length(unique(base$PKey)),"unique persons"))


#select only those with a genealogy
baseg<-base[base$PKey %in% listprofs,]
nprofs<- nrow(baseg)
nprofs

base$published <- ifelse(base$titles>0,1,0)
baseg$published <- ifelse(baseg$titles>0,1,0)

print(paste("Number of institutions",length(unique(base$IKey))))

reg1 <- lm_robust(titles~refdate, data=base,clusters=PKey)
reg2 <- lm_robust(titles~refdate+longevity+ageentry+fatherprof, data=base,clusters=PKey)
reg3 <- lm_robust(titles~refdate+longevity+ageentry+fatherprof+theo+law+med+sci, data=base,clusters=PKey)
reg4 <- lm_robust(titles~refdate+longevity+ageentry+fatherprof+theo+law+med+sci+WeakLink+exist, data=base,clusters=PKey)
reg5 <- lm_robust(titles~refdate+longevity+ageentry+fatherprof+theo+law+med+sci+WeakLink+exist+copenhagen+groningen+alund+ulund+uppsala+acaddnk+rostock+
                                                kiel+abo+dorpat+amsterdam+Adanzig+edinburgh+konigsberg+
                                                stockholm+glasgow+greifswald+Gdanzig+petersburg+franeker+haarlem+
                                                aedinburgh+auppsala+aberdeen+andrews+leiden+hamburg+trondheim+
                                                utrecht+soro, data=base,clusters=PKey)

texreg(list(reg1,reg2,reg5),include.ci = F,digits=4,
stars = c(0.01, 0.05, 0.1),booktabs=T)

write_xlsx(base,"sample-for-balance.xlsx")

#read info from siblings
for (i in 1:nprofs)
{
id<-baseg[i,"PKey"]
baseg[i,"agedeathfather"] <- -( baseg[i,"YearBorn"]-as.numeric(profs[profs$PKey==id,"YearDeathFather"]))
if (!is.na(baseg[i,"agedeathfather"]) & baseg[i,"agedeathfather"] < -1) {print(id)}
baseg[i,"SibshipSize"] <-  profs[profs$PKey==id,"SibshipSize"]
baseg[i,"AmongWhichGirls"] <-  profs[profs$PKey==id,"AmongWhichGirls"]
baseg[i,"HalfSiblings"] <-  profs[profs$PKey==id,"HalfSiblings"]
baseg[i,"NumberDescendants"] <-  profs[profs$PKey==id,"NumberDescendants"]
baseg[i,"YearDeathFather"] <-  profs[profs$PKey==id,"YearDeathFather"]
baseg[i,"OccupationFather"] <-  profs[profs$PKey==id,"OccupationFather"]
baseg[i,"RankMale"] <-  profs[profs$PKey==id,"RankMale"]
baseg[i,"Rank"] <-  profs[profs$PKey==id,"Rank"]
baseg[i,"MaleDescendants"] <-  profs[profs$PKey==id,"MaleDescendants"]
}

baseg$SibshipSizeMale <- baseg$SibshipSize - baseg$AmongWhichGirls

if (i<nprofs){print(paste("------------Error with PKey",id))}  # if the loop stops before being completed


#build table of results

result <- data.frame(Name=character(), Nscholars=numeric(),N2scholars=numeric(),
              share=numeric(),stringsAsFactors = FALSE)

for (i in 1:ninstit)
{
result[i,]$Name <- listi[i]
result[i,]$Nscholars <- nrow(base[base$IKey==listi[i],])
result[i,]$N2scholars <- nrow(baseg[baseg$IKey==listi[i],])
}
result$share <- round(result$N2scholar/result$Nscholar,2)

#print statistics by institution
print(result[,c("Name","Nscholars","N2scholars","share")], row.names = FALSE)
print(c("TOTAL",sum(result$Nscholars),sum(result$N2scholars),round(sum(result$N2scholars)/sum(result$Nscholars),2)), row.names = FALSE)


rm(profs)
profs<-baseg

print("Profs with no refdate")
profs[is.na(profs$refdate),] 
profs<-profs[!is.na(profs$refdate),]  #we keep only those with a date (some guys with field=9 do not have one)

print("Profs with empty sibshipsize")
profs[profs$SibshipSizeMale==0,]  # this is a mistake, to be checked
profs<- profs[profs$SibshipSizeMale!=0,]

profs$SibshipSizeMale<-profs$SibshipSize-profs$AmongWhichGirls

#including rural/urban dimension
profs <- merge(profs,rural,by = "BirthPlace", all.x = TRUE)

write_xlsx(profs,"siblings-extended.xlsx")

sink()