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


# to import xl files
library(readxl)


# 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()))


profsgood <- as.data.frame(read_excel("siblings-extended-good.xlsx"))



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


seuil <- median(profsgood$titles)
seuil33 <- quantile(profsgood$titles,.33333)
seuil66 <- quantile(profsgood$titles,.66666)

profstop <- profsgood[profsgood$titles>seuil,]
profsbottom <- profsgood[profsgood$titles<=seuil,]


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

#plot of coefficient in recursive estimation

date <- (1:pts)
sign <-(1:pts)
t_test_result <- (1:pts)
p_values <- (1:pts)
t <- (1:pts)
maletop <- (1:pts)
malebottom <- (1:pts)


ns<-nrow(profsgood)

#  splitting at the median

for (i in (1:pts)){

smpl<-profsgood[as.integer(1+(i-1)*ns*(1-sha)/(pts-1)):as.integer((i-1)*ns*(1-sha)/(pts-1)+ns*sha),]

profstop <- smpl[smpl$titles>seuil,]
profsbottom <- smpl[smpl$titles<=seuil,]

 
  maletop[i]=mean(profstop[,"MaleDescendants"],na.rm=T)
  malebottom[i]=mean(profsbottom[,"MaleDescendants"],na.rm=T)
  

date[i]=round(mean(smpl$refdate),0)
t_test_result <- t.test(profstop[,"MaleDescendants"], profsbottom[,"MaleDescendants"])
p_values[i] <- t_test_result$p.value

}



profsgood[!is.na(profsgood$MaleDescendants) & !is.na(profsgood$NumberDescendants) & profsgood$MaleDescendants>profsgood$NumberDescendants,"PKey"]



pdf(file="background_children.pdf", width = 11, height = 7)

plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
     main="Male Children of scholars by level of publications",
     sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
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]))
lines(t,malebottom,lty=1,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
       col=c("coral2", "darkgreen"), lty=1,lwd=3,
       box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
  if (p_values[i] > 0.1) {
    rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightpink", alpha.f = 0.3), border = NA, add = TRUE)
  } else {
    rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgreen", alpha.f = 0.3), border = NA, add = TRUE)
  }
}

dev.off()


pdf(file="background_children.pdf", width = 11, height = 7)

plot(t,maletop,lwd=3,type="l",xlab="Year of birth", ylab="fertility (male)",frame.plot=TRUE,ylim=c(1,2.5),
     main="Male Children of scholars by level of publications",
     sub=paste("Final sample. Nobs=",nrow(profsgood),". ",pts," points, sample=",sha,sep=""),col="coral2",xaxt='n')
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]))
lines(t,malebottom,lty=4,col="darkgreen",lwd=3)
legend("topright", legend=c("Above median number of publications", "Below median number of publications"),
       col=c("coral2", "darkgreen"), lty=c(1,3),lwd=3,
       box.lty=0)
#polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
for (i in 1:length(sign)) {
  if (p_values[i] <= 0.1) {
    rect(xleft = t[i] - 0.5, xright = t[i] + 0.5, ybottom = .9, ytop = 2.6, col = adjustcolor("lightgray", alpha.f = 0.3), border = NA, add = TRUE)
  }
}

dev.off()




