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


# to import xl files
library(readxl)
#to write excel files
require(writexl)
# to output towards latex
#require(stargazer)
#to do clustered SE
library(estimatr)
library(texreg) #alternative to stargazer


# 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("sample-for-balance.xlsx"))



sink("balance.txt", split = TRUE)


# 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 <- 0.30 # range of quantiles at each point

profs$longevity <- profs$YearDied - profs$YearBorn


##############################
# regression with quality observations
##############################

smpl0 <- profs[!is.na(profs$longevity), ]
ns <- nrow(smpl0)

models_to_print <- c(1, 10, 20, 30, 40, 50, 60, 70)


# ============================================================
# Table 1: Outcome = titles
# ============================================================

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 ~ exist +
      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$PKey
  )
  
  eqvec[[i]] <- eq1
}

cat("\n\n")
cat("============================================================\n")
cat("Outcome: Titles\n")
cat("============================================================\n\n")

texreg(
  eqvec[models_to_print],
  include.ci = FALSE,
  stars = c(0.01, 0.05, 0.1),
  booktabs = TRUE,
  longtable = TRUE
)


# ============================================================
# Table 2: Outcome = longevity
# ============================================================

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(
    longevity ~ exist +
      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$PKey
  )
  
  eqvec[[i]] <- eq1
}

cat("\n\n")
cat("============================================================\n")
cat("Outcome: Longevity\n")
cat("============================================================\n\n")

texreg(
  eqvec[models_to_print],
  include.ci = FALSE,
  stars = c(0.01, 0.05, 0.1),
  booktabs = TRUE,
  longtable = TRUE
)


# ============================================================
# Table 3: Outcome = sci
# ============================================================

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(
    sci ~ exist +
      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$PKey
  )
  
  eqvec[[i]] <- eq1
}

cat("\n\n")
cat("============================================================\n")
cat("Outcome: Scientific output\n")
cat("============================================================\n\n")

texreg(
  eqvec[models_to_print],
  include.ci = FALSE,
  stars = c(0.01, 0.05, 0.1),
  booktabs = TRUE,
  longtable = TRUE
)


# ============================================================
# Stop saving output to the text file
# ============================================================

sink()

cat("\nAll regression tables have been saved\n")




