## 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
library(estimatr)
library(dplyr)

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


 profs <- as.data.frame(read_excel("siblings-extended.xlsx"))
 profsgood <- as.data.frame(read_excel("siblings-extended-good.xlsx"))
wc <- as.data.frame(read_excel("Dictionary-worldcat.xlsx"))

profsgood$nworks <- asinh(wc$nworks[match(profsgood$PKey,wc$PKey)])
sum(is.na(profsgood$nworks))

summary(profsgood$titles)
summary(profsgood$nworks)




#here we run the regression on full sample to see if everything is right

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




#################################
profsgood  <- profsgood[!is.na(profsgood$longevity),]
profs  <- profs[!is.na(profs$longevity),]


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



# clusters
profsgood$cluster<-profsgood$PKey
profs$cluster<-profs$PKey


################################################################################
# MAIN REGRESSION regression with quality observations
##################################################################################

smpl0 <- profsgood[!is.na(profsgood$nworks),]
ns<-nrow(smpl0)

eqvec <- list()
#initialization
bet <- (1:pts)
se <- (1:pts)
sbet <- (1:pts)
sse <- (1:pts)
rrr <- (1:pts)
t<- bet
df <- (1:pts)
high10 <- (1:pts)
low10 <- (1:pts)
high05 <- (1:pts)
low05 <- (1:pts)

for (i in (1:pts)){

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

  eq1 <- lm_robust(nworks ~ SibshipSizeMale+
                                longevity+ageentry +law+med+theo+sci+ 
                                WeakLink+ HalfSiblings+ Urban+fatherprof+geni+ soctop+socmid+ 
                                copenhagen+groningen+franeker+alund+ulund+uppsala+acaddnk+rostock+kiel+
                                abo+dorpat+amsterdam+Adanzig+Gdanzig+edinburgh+konigsberg+utrecht+
                                stockholm+glasgow+greifswald+petersburg+haarlem+aedinburgh+auppsala+aberdeen+andrews+leiden+hamburg+
                                trondheim+soro,
                                data=smpl,clusters=smpl$cluster)
  bet[i] <- coef(eq1)[[2]]
  se[i] <- sqrt(diag(vcov(eq1)))[[2]]
  sbet[i] <- coef(eq1)[[3]]
  sse[i] <- sqrt(diag(vcov(eq1)))[[3]]
  rrr[i] <- summary(eq1)$r.squared

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.95,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.975,df[i])*se[i]

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

pdf(file="recursive-preferred-good-worldcat.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.18,.13),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",nrow(smpl0)),
     ylab="Reg. coef. of Sibshipsize (male) on publications",frame.plot=TRUE,cex=1.3,
     main="(H) Rolling Regression - Profs with genealogy of good quality - Worldcat" ,axes=FALSE,lwd=3)
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=c(-.15,-.1,-.05,0,.05,.1),labels=c(-.15,-.1,-.05,0,.05,.1))
polygon(x=c(t,rev(t)),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA)
polygon(x=c(t,rev(t)),y=c(high10,rev(low10)), col =  adjustcolor("darkgray", alpha.f = 0.20), border = NA)
lines(t,rep(0,pts),lty=3,lwd=2)
lines(t,rep(0.05,pts),lty=3,lwd=1)
lines(t,rep(0.1,pts),lty=3,lwd=1)
lines(t,rep(-0.05,pts),lty=3,lwd=1)
lines(t,rep(-0.1,pts),lty=3,lwd=1)
lines(t,rep(-0.15,pts),lty=3,lwd=1)
dev.off()

