## 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(texreg) #alternative to stargazer
library(boot)  #to compute SE of correlation
library(dplyr)

# Load RColorBrewer
library(RColorBrewer)
# Build a color palette,
coul <- colorRampPalette(brewer.pal(26, "Set3"))(21)

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



# Defining the extensive margin of publications

profsgood$extensive <- ifelse(profsgood$titles==0, 0, 1)

#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

#test that institutions dummies are ok
profsgood$test <- profsgood$copenhagen+profsgood$groningen+profsgood$franeker+profsgood$alund+profsgood$ulund+
             profsgood$uppsala+profsgood$acaddnk+profsgood$rostock+profsgood$kiel+profsgood$haarlem+
             profsgood$abo+profsgood$dorpat+profsgood$amsterdam+profsgood$Adanzig+profsgood$Gdanzig+profsgood$edinburgh+
             profsgood$konigsberg+profsgood$stockholm+profsgood$glasgow+profsgood$greifswald+profsgood$petersburg+
             profsgood$aedinburgh+profsgood$auppsala+profsgood$aberdeen+profsgood$andrews+profsgood$leiden+
             profsgood$hamburg+profsgood$trondheim+profsgood$utrecht+profsgood$soro

mean(profsgood$test)


profsgood[is.na(profsgood$NumberDescendants),"NumberDescendants"] <-0 # to be changed when NA clarified (needs access to myheritage)

profsgood$Urban <- ifelse(is.na(profsgood$Urban) | profsgood$Urban==0,0,1)
profsgood$Girl <-ifelse(profsgood$AmongWhichGirls>0,1,0)
profsgood$GM <- profsgood$AmongWhichGirls/profsgood$SibshipSizeMale

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


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

sink("statdes.txt")
require(stargazer)
stargazer(profsgood)
sink()

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


#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)
datemin <- (1:pts)
datemed <- (1:pts)
datemax <- (1:pts)
nobser <- (1:pts)
mediandate <- (1:pts)
meanfertility <- (1:pts)
sdfertility <- (1:pts)
meanlongevity <- (1:pts)
meanageentry <- (1:pts)
meanlaw <- (1:pts)
meanmed <- (1:pts)
meantheo <- (1:pts)
meansci <- (1:pts)
meanWeakLink <- (1:pts)
meanHalfSiblings <- (1:pts)
meanUrban <- (1:pts)
meanfatherprof <- (1:pts)
meangeni <- (1:pts)
meansoctop <- (1:pts)
meansocmid <- (1:pts)
meancopenhagen <- (1:pts)
meangroningen <- (1:pts)
meanfraneker <- (1:pts)
meanalund <- (1:pts)
meanulund <- (1:pts)
meanuppsala <- (1:pts)
meanacaddnk <- (1:pts)
meanrostock <- (1:pts)
meankiel <- (1:pts)
meanabo <- (1:pts)
meandorpat <- (1:pts)
meanamsterdam <- (1:pts)
meanAdanzig <- (1:pts)
meanGdanzig <- (1:pts)
meanedinburgh <- (1:pts)
meankonigsberg <- (1:pts)
meanstockholm <- (1:pts)
meanglasgow <- (1:pts)
meangreifswald <- (1:pts)
meanpetersburg <- (1:pts)
meanhaarlem <- (1:pts)
meanaedinburgh <- (1:pts)
meanauppsala <- (1:pts)
meanaberdeen <- (1:pts)
meanandrews <- (1:pts)
meanleiden  <- (1:pts)
meanhamburg  <- (1:pts)
meantrondheim<- (1:pts)
meanutrecht  <- (1:pts)
meansoro  <- (1:pts)





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

smpl0 <- profsgood
ns<-nrow(smpl0)

eqvec <- list()


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(titles ~ 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

meanfertility[i] <- mean(smpl$SibshipSizeMale)
sdfertility[i] <- sd(smpl$SibshipSizeMale)

nobser[i]<-eq1$nobs
  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]

  eqvec[[i]]<-eq1
  datemin[i]<- min(smpl$refdate)
datemed[i]<-median(smpl$refdate)
  datemax[i] <- max(smpl$refdate)
mediandate[i] <- round(median(smpl$refdate),0)
meanlongevity[i] <- mean(smpl$longevity)
meanageentry[i] <- mean(smpl$ageentry)
meanlaw[i] <- mean(smpl$law)
meanmed[i] <- mean(smpl$med)
meantheo[i] <- mean(smpl$theo)
meansci[i] <- mean(smpl$sci)
meanWeakLink[i] <- mean(smpl$WeakLink)
meanHalfSiblings[i] <- mean(smpl$HalfSiblings)
meanUrban[i] <- mean(smpl$Urban)
meanfatherprof[i] <- mean(smpl$fatherprof)
meangeni[i] <- mean(smpl$geni)
meansoctop[i] <- mean(smpl$soctop)
meansocmid[i] <- mean(smpl$socmid)
meancopenhagen[i] <- mean(smpl$copenhagen)
meangroningen[i] <- mean(smpl$groningen)
meanfraneker[i] <- mean(smpl$franeker)
meanalund[i] <- mean(smpl$alund)
meanulund[i] <- mean(smpl$ulund)
meanuppsala[i] <- mean(smpl$uppsala)
meanacaddnk[i] <- mean(smpl$acaddnk)
meanrostock[i] <- mean(smpl$rostock)
meankiel[i] <- mean(smpl$kiel)
meanabo[i] <- mean(smpl$abo)
meandorpat[i] <- mean(smpl$dorpat)
meanamsterdam[i] <- mean(smpl$amsterdam)
meanAdanzig[i] <- mean(smpl$Adanzig)
meanGdanzig[i] <- mean(smpl$Gdanzig)
meanedinburgh[i] <- mean(smpl$edinburgh)
meankonigsberg[i] <- mean(smpl$konigsberg)
meanstockholm[i] <- mean(smpl$stockholm)
meanglasgow[i] <- mean(smpl$glasgow)
meangreifswald[i] <- mean(smpl$greifswald)
meanpetersburg[i] <- mean(smpl$petersburg)
meanhaarlem[i] <- mean(smpl$haarlem)
meanaedinburgh[i] <- mean(smpl$aedinburgh)
meanauppsala[i] <- mean(smpl$auppsala)
meanaberdeen[i] <- mean(smpl$aberdeen)
meanandrews[i] <- mean(smpl$andrews)
meanleiden[i] <- mean(smpl$leiden)
meanhamburg[i] <- mean(smpl$hamburg)
meantrondheim[i] <- mean(smpl$trondheim)
meanutrecht[i] <- mean(smpl$utrecht)
meansoro[i] <- mean(smpl$soro)
}





################################################################################
# Figure: academic fields over time
################################################################################





# Load RColorBrewer
library(RColorBrewer)
# Build a color palette,
coul <- colorRampPalette(brewer.pal(26, "Set3"))(21)

pdf(file="fields.pdf", width = 10, height = 7)
plot(meanlaw,type="l",xlab="Year of Birth", ylab="Number of persons",lwd=2,
frame.plot=TRUE,ylim=c(0,.5),col=coul[1],
     main="Proportions of academic fields 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(mediandate[1]),mediandate[10],mediandate[20],mediandate[30],
mediandate[40],mediandate[50],mediandate[60],mediandate[pts]))
axis(2,at=seq(0,.6, by=.1),labels=c(0,.1,.2,.3,.4,.5,.6))
lines(meantheo,col=coul[6],lwd=2)
lines(meanmed,col=coul[4],lwd=2)
lines(meansci,col=coul[2],lwd=2)
legend(10, .5, legend=c("Law", "Theology","Medicine","Science"),
       col=c(coul[1],coul[6],coul[4],coul[2]),  lwd=2,box.lty=0)
dev.off()



sink("rolling.txt")
cbind(datemin,datemed,datemax,round(bet,3),round(low10,3),round(high10,3),nobser,df+1,meanfertility,sdfertility)
sink()


write_xlsx(as.data.frame(cbind(datemin,datemax,meanfertility,meanlongevity,meanageentry,meanlaw,meanmed,meantheo,meansci,
meanWeakLink,meanHalfSiblings,meanUrban,meanfatherprof,meangeni,meansoctop,meansocmid,
meancopenhagen,meangroningen,meanfraneker,meanalund,
meanulund,meanuppsala,meanacaddnk,meanrostock,meankiel,meanabo,meandorpat,meanamsterdam,meanAdanzig,
meanGdanzig,meanedinburgh,meankonigsberg,meanstockholm,meanglasgow,meangreifswald,meanpetersburg,
meanhaarlem,meanaedinburgh,meanauppsala,meanaberdeen,meanandrews,meanleiden,meanhamburg,meantrondheim,meansoro,meanutrecht)),"rolling-allvar.xlsx")

lcoef <-bet
scoef <-sbet
rsquare <- rrr

write_xlsx(as.data.frame(t(cbind(lcoef,scoef,rsquare))),"output-for-table.xlsx")

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

pdf(file="recursive-preferred-good.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Rolling Regression - Final sample" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()


#The trick with polygon is that you must provide 2 times the x coordinates in one vector, once in normal order and
# once in reverse order (with function rev) and you must provide the y coordinates as a vector of the upper bounds
# followed by the lower bounds in reverse order.
#We use the adjustcolor function to make standard colors transparent.

sink('allcoef.txt')
texreg(list(eqvec[[1]],eqvec[[10]],eqvec[[20]],eqvec[[30]],eqvec[[40]],eqvec[[50]],eqvec[[60]],eqvec[[70]]),include.ci = F,digits=3,
stars = c(0.01, 0.05, 0.1),booktabs=T,longtable=T)
sink()



##############################
#regression with quality observations - girls
##############################

smpl0 <- profsgood
ns<-nrow(smpl0)

eqvec <- list()

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(titles ~ AmongWhichGirls+ 
                                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

meanfertility[i] <- mean(smpl$SibshipSizeMale)
sdfertility[i] <- sd(smpl$SibshipSizeMale)

nobser[i]<-eq1$nobs
  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]
}

lcoef <-bet
scoef <-sbet
rsquare <- rrr


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

pdf(file="recursive-preferred-good-girls.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.13),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",nrow(smpl0)),
     ylab="Reg. coef. of Sibshipsize (female) on publications",frame.plot=TRUE,cex=1.3,
     main="(A) Rolling Regression - No. of sisters as dependent variable" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()

##############################
#Robusteness - regression with SibShipSize and quality observations
##############################

smpl0 <- profsgood
ns<-nrow(smpl0)

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(titles ~ SibshipSize + 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="robustness-sibshipsize.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.13),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",nrow(smpl0)),
     ylab="Reg. coef. of Sibshipsize on publications",frame.plot=TRUE,cex=1.3,
     main="(B) Rolling Regression - No. of sisters and brothers as dependent variable" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()


##############################
#robustness - regression with all observations
##############################


smpl0 <- profs ######################here we take all profs
ns<-nrow(smpl0)

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(titles ~ SibshipSizeMale + longevity+ageentry +law+med+theo+sci + 
                                WeakLink+HalfSiblings +Urban+fatherprof+geni+
                                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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="robustness-all-genealogies.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(A) Rolling Regression - Full sample" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()





##############################
#regression with quality observations - intensive margin of publications only
##############################


smpl0 <- profsgood[profsgood$titles>0,]
ns<-nrow(smpl0)

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(titles ~ SibshipSize + 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]]

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

cbind(datemin,datemax,round(bet,3),round(se,3),nobser,df+1)


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


pdf(file="robustness-intensive-margin.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(B) Published profs, final sample" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()



##############################
#regression with quality observations - extensive margin of publications 
##############################

smpl0 <- profsgood
ns<-nrow(smpl0)

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(extensive ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="robustness-extensive-margin.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-0.02,0.02),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",nrow(smpl0)),
     ylab="Reg. coef. of Sibshipsize (male) on probability to publish",frame.plot=TRUE,cex=1.3,
     main="(C) Final sample: probability to publish" ,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=seq(-.02,.02,by=.01),labels=c(-.02,-.01,0,.01,.02))
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)
dev.off()

##############################
#Mechanism - regression with quality observations without Weaklinks
##############################


smpl0 <- profsgood[profsgood$WeakLink==F,]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-weaklink.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.13),xlab=paste("Year of Birth. Window = 30% of sample. No. obs=",nrow(smpl0)),
     ylab="Reg. coef. of Sibshipsize (male) on publications",frame.plot=TRUE,cex=1.3,
     main="(D) Final sample. No Weaklinks." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()




##############################
#regression with quality observations without  a field
##############################


smpl0 <- profsgood[profsgood$theo==0,]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-theo.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No Theology." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()





smpl0 <- profsgood[profsgood$law==0,]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-law.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No Law." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()





smpl0 <- profsgood[profsgood$sci==0,]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-sci.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(F) Final sample. No Science." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()




smpl0 <- profsgood[profsgood$med==0,]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-med.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No Medicine." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()



##############################
#regression with quality observations - no french
##############################


smpl0 <- profsgood[profsgood$Locations_Country2016!="FRA",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-french.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(E) Final sample. No French." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()

##############################
#regression with quality observations - no swedish
##############################


smpl0 <- profsgood[profsgood$Locations_Country2016!="SWE",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-swedish.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(E) Final sample. No Swedish." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()


##############################
#regression with quality observations - no denmark
##############################

table(profsgood$Locations_Country2016)

smpl0 <- profsgood[profsgood$Locations_Country2016!="DNK",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-danish.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(E) Final sample. No Danish." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()


##############################
#regression with quality observations - no finnish
##############################

table(profsgood$Locations_Country2016)

smpl0 <- profsgood[profsgood$Locations_Country2016!="FIN",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-finnish.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(E) Final sample. No Finnish." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()


##############################
#regression with quality observations - no NA in place of birth
##############################

table(profsgood$Locations_Country2016)

smpl0 <- profsgood[profsgood$Locations_Country2016!="NA",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-NA.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(E) Final sample. No NA." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()



##############################
#regression with quality observations - no british
##############################


smpl0 <- profsgood[profsgood$Locations_Country2016!="GBR" & profsgood$Locations_Country2016!="USA",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-british.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No British." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()





##############################
#regression with quality observations - no german
##############################


smpl0 <- profsgood[profsgood$Locations_Country2016!="DEU",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-german.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No German." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()



##############################
#regression with quality observations - no dutch
##############################

smpl0 <- profsgood[profsgood$Locations_Country2016!="NLD",]
ns<-nrow(smpl0)

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(titles ~ 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]]

  df[i] <- eq1$nclusters-1  # degrees of freedom 
   high10[i]<- bet[i]+qt(0.975,df[i])*se[i]
   low10[i]<- bet[i]-qt(0.975,df[i])*se[i]
   high05[i]<- bet[i]+qt(0.95,df[i])*se[i]
   low05[i]<- bet[i]-qt(0.95,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="mechanism-no-dutch.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="Final sample. No Dutch." ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()





##############################
#regression with quality observations and lots of controls
##############################

smpl0 <- profsgood
ns<-nrow(smpl0)

eqvec <- list()


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(titles ~ SibshipSizeMale + longevity+ageentry +law+med+theo+sci + 
                                WeakLink+
                                 HalfSiblings + Urban+fatherprof+geni+soctop+socmid+
                                elder + 
                                  socbot+
                                log(1+NumberDescendants) +
                                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

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

  eqvec[[i]]<-eq1
  datemin[i]<- min(smpl$refdate)
  datemax[i] <- max(smpl$refdate)
}



lcoef <-bet
scoef <-sbet
rsquare <- rrr


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

pdf(file="recursive-preferred-good-full-control.pdf", width = 10, height = 7)
plot(t,bet,type="l",ylim=c(-.13,.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="(G) Rolling Regression - Final sample - All controls" ,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(-.1,-.05,0,.05,.1),labels=c(-.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)
dev.off()



######################################
### Correlation of fertility rates ###
######################################


correlation <- (1:pts)
smpl0 <- profsgood
ns<-nrow(smpl0)

ci_list <- vector("list", length = pts)

correlation_function <- function(data, i) {
  subset_data <- data[i, ]
  correlation <- cor(subset_data$SibshipSize, subset_data$NumberDescendants, use = "complete.obs")
  return(correlation)
}


num_bootstrap_samples <- 100

correlation <- numeric(pts)
ci_list <- vector("list", length = pts)

# Loop through subsamples
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),]


  # Perform bootstrap
  boot_results <- boot(data = smpl, statistic = correlation_function, R = num_bootstrap_samples)
  
  # Calculate confidence interval for the correlation
  ci <- quantile(boot_results$t, c(0.025, 0.975))
  low05[i] <- quantile(boot_results$t, 0.025)
  high05[i] <- quantile(boot_results$t, 0.975)
  low10[i] <- quantile(boot_results$t, 0.05)
  high10[i] <- quantile(boot_results$t, 0.95)
  
  # Store the lower and upper bounds in the ci_list
  ci_list[[i]] <- ci
  correlation[i] <- cor(smpl$SibshipSize, smpl$NumberDescendants, use = "complete.obs")
}

# Print the correlation and confidence intervals
for (i in 1:pts) {
  cat("Subsample", i, ": Correlation =", correlation[i], ", 95% CI =", ci_list[[i]], "\n")
}


corr <-correlation


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

pdf(file="recursive-correlations.pdf", width = 10, height = 7)
plot(t,corr,type="l",ylim=c(-.15,.25),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",
               nrow(smpl0)),
     ylab="Pearson correlation",frame.plot=TRUE,cex=1.3,
     main="Correlation between parental and professoral fertility" ,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(-.1,-.05,0,.05,.1,.15),labels=c(-.1,-.05,0,.05,.1,.15))
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.15,pts),lty=3,lwd=1)
lines(t,rep(0.2,pts),lty=3,lwd=1)
lines(t,rep(0.25,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()




######################################
### Correlation of fertility rates - male ###
######################################


correlation <- (1:pts)
smpl0 <- profsgood
ns<-nrow(smpl0)

ci_list <- vector("list", length = pts)

correlation_function <- function(data, i) {
  subset_data <- data[i, ]
  correlation <- cor(subset_data$SibshipSizeMale, subset_data$MaleDescendants, use = "complete.obs")
  return(correlation)
}


num_bootstrap_samples <- 100

correlation <- numeric(pts)
ci_list <- vector("list", length = pts)

# Loop through subsamples
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),]
  
  # Perform bootstrap
  boot_results <- boot(data = smpl, statistic = correlation_function, R = num_bootstrap_samples)
  
  # Calculate confidence interval for the correlation
  ci <- quantile(boot_results$t, c(0.025, 0.975))
  low05[i] <- quantile(boot_results$t, 0.025)
  high05[i] <- quantile(boot_results$t, 0.975)
  low10[i] <- quantile(boot_results$t, 0.05)
  high10[i] <- quantile(boot_results$t, 0.95)
  
  # Store the lower and upper bounds in the ci_list
  ci_list[[i]] <- ci
  correlation[i] <- cor(smpl$SibshipSizeMale, smpl$MaleDescendants, use = "complete.obs")
}

# Print the correlation and confidence intervals
for (i in 1:pts) {
  cat("Subsample", i, ": Correlation =", correlation[i], ", 95% CI =", ci_list[[i]], "\n")
}


corr <-correlation

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

pdf(file="recursive-correlations-male.pdf", width = 10, height = 7)
plot(t,corr,type="l",ylim=c(-.15,.25),xlab=paste("Year of Birth. Window = 30% of sample. Robust SE. No. obs=",
               nrow(smpl0)),
     ylab="Pearson correlation",frame.plot=TRUE,cex=1.3,
     main="Correlation between parental and professoral fertility (male)" ,axes=FALSE,lwd=3)
axis(1,at=c(1,10,20,30,40,50,60,pts),labels=as.integer(c(date[1],date[10],date[20],date[30],date[40],date[50],date[60],date[pts])))
axis(2,at=c(-.1,-.05,0,.05,.1,.15),labels=c(-.1,-.05,0,.05,.1,.15))
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.15,pts),lty=3,lwd=1)
lines(t,rep(0.2,pts),lty=3,lwd=1)
lines(t,rep(0.25,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()


