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

# Load RColorBrewer
library(RColorBrewer)

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

# 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 data
profs <- as.data.frame(read_excel("siblings-extended.xlsx"))

sink("sink from demographic transition.txt")

#select quality observations

#we drop those with less than 2 sibships and year of death of father is unkwown
profsgood <- profs[profs$SibshipSize>2 | (profs$SibshipSize<3 & !is.na(profs$YearDeathFather)),]
print(paste("criterion 1 leads to drop",nrow(profs)-nrow(profsgood),"observations"))
tp<-nrow(profsgood)
#we drop those with one sibship and one descendent
profsgood <- profsgood[profsgood$SibshipSize>1 | profsgood$NumberDescendants>1,]
print(paste("criterion 2 leads to drop",tp-nrow(profsgood),"observations"))


#table of sibshipsizemale to calibrate phi
table(profsgood$SibshipSizeMale)
summary(profsgood$SibshipSizeMale)


# include a rank variable

profsgood$elder <- ifelse(!is.na(profsgood$RankMale) & profsgood$RankMale==1,1,0)

# include a social class variable

strtop <- gsub("[\r\n]", "","archbishop|knighted|marquess|peer|prince|regent|aristocrat|theologian|linguist|king|title|seigneur|professor|councillor|
bishop|mayor|doctor|rector|general|governor|lord|chief justice|colonel|captain|advocate|baron|physician|baronet|owner|lawyer|author|dr|prosecutor|
lecturer|deputy|magistrate|curator|senator|advisor|king's|medical|scholar|writer|count|advocat|attorney|magister|procurator|commissar|businessman|
physicist|marquis|engineer|manager|politician|landlord|artist|orientalist|entrepreneur|industrialist|duke|geographer|conferenceraad|biskop|anatomist|
consultant|conventcommissary|noble|lieutenant-colonel|confessor|reichsmarschall|comminister|barrister|commandant|councilor|co-director|commonwealth marsh|
prince-episcopal|humanist|magistrat|drtheology|philologist|pedagogue|notaire|banker|exchequer|hetman|geologist|adjutant-general|baroness|lagman|adlad|
eversti|oberst|prepositus|logician|chatelain|cardinal|nobleman|prelate|cartographer|councillor-secretary|statesman|counsellor|vice-admiral|assessor|
chancellor|attorney-postulant|investor|apothecary|judge|architect|antiquary|consul|sir|earl|knight|colonel|surgeon|laird|principal|historian|mathematician|
theologist|admiral|commander|jurist|scholast|notary|dean|marshall|chief marshal|provost|president|chamberlain")

strmid <- gsub("[\r\n]", "","elder|thesaurier|printer|clergyman|inventeur|poet|trumpeter|painter|violinist|oboist|shipmaster|musician|alderman|administrator|
monastic|auditor|reformer|preacher|priest|merchant|pastor|officer|trader|master|superintendent|vicar|secretary|clerk|reverend|minister|lieutenant|teacher|
bailiff|baillif|chaplain|pharmacist|inspector|quartermaster|deacon|manufacturer|skipper|schoolmaster|postmaster|foreman|tutor|councilman|chef|shipowner|bookseller|
treasurer|bourgeois|organist|civil servant|minister|lector|glazier|shipbuilder|surveyor|weighmaster|explorer|syndicus|commissioner|archivist|protocollist|churchman|
controller|turner|librarian|registrar|scribe|greffier|jeweler|wholesaler|burgher|adjutant|marksman|nurse|baccalaureus|papermaker|salesman|cavalry|real estate agent|min. of old kirk|
makler|clothier|kornett|adventurer|orfever|reverend|collector|brigadier|armourer|accounter|churchwarden|sheriff|farmer|hundertman|leaseholder|cornet|counsel")

strbot <- gsub("[\r\n]", "","minstrel|artisan|bierbrouwer|saddler|shoemaker|craftsman|policeman|soldier|cutler|factor|blacksmith|gunsmith|goldsmith|tobacconist|brewer|
cossack|piece caster|shear|cartridge|tailor|fisherman|innkeeper|gardener|baker|grocer|tanner|carpenter|assistant|engraver|maker|cutter|receiver|horseman|
employee|guard|tinsmith|comrade|weaver|coalminer|fishing|miller|slave|schmidt|coppersmith|janitor|beater|churchman|hunting|miner|shopkeeper|caretaker|
shipwright|vigneron|guards|non-commissioned|driver|peasantry|fabrics|burner|military|sweeper|forester|grocery|farm|tractor|lorch|crimond")

strunc <- gsub("[\r\n]", "","syndict|appellate|citizen|ministry|college|chancery|head|collegium|pharmacies|ordynek|consistory|seeker|writter|musices|council|
government|school|trade|administrative|confidant|traveler|buyer|communicator|registar|commercial|advice|secretariat|priech|contractual|senevier|
answerer|receiver|piece caster|shear|cartridge|grocery|farm|aids|crimond|hunting|non-commissioned|
tractor|fabrics|military|lorch")

profsgood$OccupationFather <- tolower(profsgood$OccupationFather)

#occupation dummies
nprof <- nrow(profsgood)
profsgood$soctop <- 0
profsgood$socmid <- 0
profsgood$socbot <- 0
profsgood$socna <- 0

for (i in (1:nprof))
{
if (is.na(profsgood[i,"OccupationFather"])) {profsgood[i,"socna"]<-1
} else {
if (grepl(strtop,profsgood[i,"OccupationFather"])) {profsgood[i,"soctop"]<-1
} else {
if (grepl(strmid,profsgood[i,"OccupationFather"])) {profsgood[i,"socmid"]<-1
} else { 
if (grepl(strbot,profsgood[i,"OccupationFather"])) {profsgood[i,"socbot"]<-1
} else { 
if (grepl(strunc,profsgood[i,"OccupationFather"])) {profsgood[i,"socna"]<-1}
}}}}

if (profsgood[i,"socna"]+profsgood[i,"soctop"]+profsgood[i,"socmid"]+profsgood[i,"socbot"]==0)
{print(profsgood[i,"OccupationFather"]);
profsgood[i,"socna"]<-1}
}

#check for missing occupation
print("check for missing occupation")
profsgood[profsgood$socna==1 & !is.na(profsgood$OccupationFather),c("PKey","OccupationFather")]



#save the results
write_xlsx(profsgood,"siblings-extended-good.xlsx")


profsgood2 <- profsgood[!is.na(profsgood$MaleDescendants) & !is.na(profsgood$NumberDescendants)  
                & profsgood$NumberDescendants>0,]


#LET us now do demographics

#parties from wrigley
countwrigley <- c(0.191693856,0.236878826,0.232893927,0.174008955,0.086280695,0.031133555,0.007892533)

# Simple Bar Plot
counts1b <- table(profs$SibshipSizeMale)/nrow(profs)
counts2b <- table(profsgood$SibshipSizeMale)/nrow(profsgood)

pdf("histogram-parities-all.pdf")
barplot(rbind(counts1b[1:7],t(countwrigley)),ylab="Frequency (%)",
  xlab="Number of male children",main="Distribution of Parities - males - full sample",beside=T,
legend=c("Profs","English pop (Wrigley)"))
dev.off()

pdf("histogram-parities-good.pdf")
barplot(rbind(counts2b[1:7],t(countwrigley)), ylab="Frequency (%)",
 xlab="Number of male children",main="Distribution of Parities - males - final sample",beside=T,
legend=c("Profs","English pop (Wrigley)"))
dev.off()


pdf("frequency-publications.pdf")
hist(profsgood$titles,breaks=c(0,.1,1,2,3,4,5,6,7),freq=T,main="",
ylab="Nb. Profs",xlab="Publications")
abline(v=median(profsgood$titles), lwd=3, lty='dashed')
dev.off()

# Create a migration dummy comparing country of birth and death

profsgood$migrant <- toupper(profsgood$Locations_1_Country2016) == toupper(profsgood$Locations_2_Country2016)

profsgood$sexratio <- profsgood$AmongWhichGirls/profsgood$SibshipSizeMale


# Build a color palette,
coul <- colorRampPalette(brewer.pal(31, "Set3"))(21)

# Pie Chart with Percentages
slices <- c(mean(profs$copenhagen),mean(profs$groningen),mean(profs$alund),mean(profs$ulund),mean(profs$uppsala),
mean(profs$acaddnk),mean(profs$rostock),
mean(profs$kiel),mean(profs$abo),mean(profs$dorpat),mean(profs$amsterdam),
mean(profs$Adanzig),mean(profs$edinburgh),mean(profs$konigsberg),
mean(profs$stockholm),mean(profs$glasgow),mean(profs$greifswald),
mean(profs$Gdanzig),mean(profs$petersburg),mean(profs$franeker),mean(profs$haarlem),
mean(profs$aedinburgh),mean(profs$auppsala),mean(profs$aberdeen),mean(profs$andrews),mean(profs$leiden),
mean(profs$hamburg),mean(profs$trondheim),mean(profs$utrecht),mean(profs$soro))
lbls <- c("U Copenhagen","U Groningen","A Lund","U Lund","U Uppsala","A Danish","U Rostock","U Kiel","U Abo","U Dorpat","C Amsterdam",
"A Danzig","U Edinburgh","U Konigsberg","A Stockholm","U Glasgow","U Greifswald","C Danzig",
"A St-Petersburg","U Franeker","A Haarlem","A Edinburgh","A Uppsala","U Aberdeen","U St Andrews","U Leiden","A Hesse",
"A Trondheim","U Utrecht","G Soro")
pct <- round(slices/sum(slices)*100)
lbls <- paste(lbls, pct) # add percents to labels
lbls <- paste(lbls,"%",sep="") # ad % to labels

pdf(file="share-institutions.pdf", width = 10, height = 7)
pie(slices,labels = lbls, col=coul,main=paste("Pie Chart of Institutions. Nobs=",nrow(profs)))
dev.off()


#pie with urban rural
slices <- c(nrow(profsgood[!is.na(profsgood$Urban) & profsgood$Urban==1,]),
nrow(profsgood[!is.na(profsgood$Urban) & profsgood$Urban==0,]),
nrow(profsgood[is.na(profsgood$Urban),]))
lbls <- c("Urban","Rural","NA")
pct <- round(slices/sum(slices)*100)
lbls <- paste(lbls, pct) # add percents to labels
lbls <- paste(lbls,"%",sep="") # ad % to labels
pdf(file="share-urban-rural.pdf", width = 10, height = 7)
pie(slices,labels = lbls, col=coul,main=paste("Pie Chart of Parental Background. Nobs=",nrow(profs)))
dev.off()


#time series plots

pts<-70 #number of points for which we compute the fertility
sha <- 1-pts/100  # range of quantiles at each point

profsgood2$FemaleDescendants <- profsgood2$NumberDescendants-profsgood2$MaleDescendants

#order by date
profsgood <- profsgood[order(profsgood$refdate),]
profs <- profs[order(profs$refdate),]

results <- data.frame(titlesmed=numeric(),titlesmean=numeric(),longevity=numeric(),fertility=numeric(),fertilitymean=numeric(),
                          fertilitymean2=numeric(),upperclass=numeric(),middleclass=numeric(),
                          lowerclass=numeric(),naclass=numeric(),shareurban=numeric(),
                          sexratio=numeric(),sexratio2=numeric(),sharereject=numeric(),stringsAsFactors = FALSE)

ns<-nrow(profsgood)
nsall<-nrow(profs)

#generate time series
for (i in (1:pts)){

t1 <- profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)),"refdate"]
t2 <- profsgood[as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"refdate"]

results[i,"sharereject"] <- 1-
nrow(profsgood[profsgood$refdate>=t1 & profsgood$refdate<t2,])/
nrow(profs[profs$refdate>=t1 & profs$refdate<t2,])

results[i,"shareurban"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"Urban"],na.rm=T)
results[i,"titlesmed"] <- median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"titles"])
results[i,"titlesmean"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"titles"])
results[i,"upperclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"soctop"],na.rm=T)
results[i,"middleclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"socmid"],na.rm=T)
results[i,"lowerclass"] <- mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"socbot"],na.rm=T)
results[i,"longevity"] <- median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"longevity"],na.rm=T)
results[i,"fertility"] <- 2*median(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"])
results[i,"fertilitymean"] <- 2*mean(
profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"])

results[i,"fertilitymean2"] <- 2*mean(
profsgood2[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"MaleDescendants"])

results[i,"sexratio"] <- sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"AmongWhichGirls"])/
               sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"SibshipSizeMale"]-1)

results[i,"sexratio2"] <- sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"FemaleDescendants"])/
               sum(profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),"MaleDescendants"])
}

results$naclass=1-results$upperclass-results$middleclass-results$lowerclass



date <- (1:pts)
for (i in (1:pts)){
date[i]=round(quantile(profsgood$refdate,i/100+.15)[[1]],0)
}

#export longevity for later treatment
outt <- data.frame(
  date = date,
  longevity = unlist(results$longevity)
)

write.csv(outt, "UTHC.csv", row.names = FALSE)




pdf(file="socialclass.pdf", width = 10, height = 7)
plot(results[,"upperclass"],type="l",xlab="Year of Birth", ylab="Share",lwd=2,
     frame.plot=TRUE,ylim=c(0,.6),col="#000000", lty=1,
     main="Proportions of social classes over time",sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(0,.6, by=.1),labels=c(0,.1,.2,.3,.4,.5,.6))
lines(results[,"middleclass"],col="#0072B2",lwd=2,lty=2)
lines(results[,"lowerclass"],col="#D55E00",lwd=3,lty=3)
lines(results[,"naclass"],col="#009E73",lwd=2,lty=4)
legend(10, .3, legend=c("Upper class", "Middle class","Lower class","NA"),
       col=c("#000000","#0072B2","#D55E00","#009E73"), lwd=2, lty=c(1,2,3,4), box.lty=0)
dev.off()






pdf(file="ntitles.pdf", width = 10, height = 7)
plot(results[,"titlesmed"],type="l",xlab="Year of Birth", ylab="Share",lwd=2,col=coul[2],
frame.plot=TRUE,main="Number of titles (asinh)",ylim=c(2.6,4),
     sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(2.6,4, by=.2),labels=c(2.6,2.8,3,3.2,3.4,3.6,3.8,4))
lines(results[,"titlesmean"],col=coul[4],lwd=2)
legend(10,4 , legend=c("Median", "Mean"),
       col=c(coul[2],coul[4]),  lwd=2,box.lty=0)
dev.off()

pdf(file="shareurban.pdf", width = 10, height = 7)
plot(results[,"shareurban"],type="l",xlab="Year of Birth", ylab="Share",
frame.plot=TRUE,main="Share of scholars born in cities",ylim=c(.3,.6),
     sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(.3,.6, by=.1),labels=c(0.3,0.4,0.5,0.6))
dev.off()


pdf(file="fertility.pdf", width = 10, height = 7)
plot(results[,"fertility"],type="l",xlab="Year of Birth", ylab="Median sibshipsize",ylim=c(4,8),
frame.plot=TRUE,main="Median size of families",
     sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(5,7, by=1),labels=c(5,6,7))
dev.off()



pdf(file="fertilitymean.pdf", width = 10, height = 7)
plot(results[,"fertilitymean"],type="l",xlab="Year of Birth", ylab="Mean sibshipsize",ylim=c(4,8),
frame.plot=TRUE,main="Mean sibshipsize",
     sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(5,7, by=1),labels=c(5,6,7))
dev.off()

pdf(file="sexratio.pdf", width = 10, height = 7)
plot(results[,"sexratio"],type="l",xlab="Year of Birth", ylab="Sex ratio",
     frame.plot=TRUE,ylim=c(.7,1),main="Sex Ratio (F/M)",
     sub=paste("Final sample. Nobs=",nrow(profsgood)),axes=FALSE)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=c(as.integer(date[1]),date[10],date[20],
                   date[30],date[40],date[50],date[60],date[pts]))
axis(2,at=seq(0.6,1, by=.1),labels=c(.6,.7,.8,.9,1))
dev.off()


proba_boy <- 1.05 / 2.05

parities <- 1:5

tsex_ratio <- sapply(parities, function(i) {
  (1 + proba_boy * i) / ((1 - proba_boy) * i)
})

esex_ratio <- sapply(parities, function(i) {
  sum(profsgood[profsgood$SibshipSize == i + 1, "SibshipSizeMale"]) /
    sum(profsgood[profsgood$SibshipSize == i + 1, "AmongWhichGirls"])
})

tsex_ratio_limit <- proba_boy / (1 - proba_boy)

sex_ratio_table <- data.frame(
  Parity = c(as.character(parities + 1), "Infinity"),
  Empirical_sex_ratio = c(esex_ratio, NA),
  Theoretical_sex_ratio = c(tsex_ratio, tsex_ratio_limit)
)

print("sex ratio analysis")
print(sex_ratio_table)

#rank analysis
print("rank analysis")
table(profsgood$RankMale)

nrow(profsgood[!is.na(profsgood$RankMale) & profsgood$RankMale==1,])

nrow(profsgood[!is.na(profsgood$RankMale),])

nrow(profsgood[!is.na(profsgood$RankMale) & profsgood$RankMale==1,])/
nrow(profsgood[!is.na(profsgood$RankMale),])

sink()