rm()



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


#to be copied from the output of simulation-new: "report-indirect-inference.txt"


h0 <- 1   # initial condition
pop0 <-1
phi <- 1/11  # 11 sons is a maximum
psi<-15.45156
gamma<-0.18662
a<-3.0108
cbar<-1.2066*a
etamean <- 0.07875

zeta<-3.078
sigma<-0.34623

etalow <- gamma/(1+phi*psi)

nhouseholds<-600
horizon<-20
time <- 1:horizon



p<-9  #we set model date p = 1635
etamin<-etalow
etamax<-etamean+(etamean-etalow)

n <- matrix(0, nrow = nhouseholds, ncol = horizon)
h <- matrix(0, nrow = nhouseholds, ncol = horizon+1)
p <- matrix(0, nrow = nhouseholds, ncol = horizon+1)
y <- matrix(0, nrow = nhouseholds, ncol = horizon+1)
eta <- rep(0,nhouseholds)

#we assume a beta distribution for eta, and the i families receive
#                  the (i-1)/(nhouseholds-1) quantile of the distribution
for (i in (1:nhouseholds)){
eta[i] <- etamin + (etamax-etamin) * qbeta((i-1)/(nhouseholds-1), zeta, zeta)
}

ilow <- which.min(abs(eta - quantile(eta,.25)))   # index of eta closest to y
imed <- which.min(abs(eta - quantile(eta,.5)))   # index of eta closest to y
ihigh <- which.min(abs(eta - quantile(eta,.75)))   # index of eta closest to y



for (t in (1:horizon)){
for (i in (1:nhouseholds)){

h[i,t] <- (psi*phi*eta[i]/(gamma-eta[i]))^(t-1) *h0

n[i,t] <- max(0,min( (gamma-eta[i])/gamma * (h[i,t]+a-cbar)/phi/h[i,t] ,
          (gamma-eta[i])/(1+gamma) * (h[i,t]+a)/phi/h[i,t]      ))

p[i,t] <- ifelse(t==1, pop0, p[i,t-1]*n[i,t])

y[i,t] <- (1-phi*n[i,t])*h[i,t]+a

}}





nlow <- rep(1,horizon-1)
nhigh <- rep(1,horizon-1)
nmed <- rep(1,horizon-1)
incmed <-  rep(1,horizon-1)
inchigh <-  rep(1,horizon-1)
inclow <-  rep(1,horizon-1)
grmed <-  rep(1,horizon-1)
grhigh <-  rep(1,horizon-1)
grlow <-  rep(1,horizon-1)

popmed <-  rep(1,horizon-1)
pophigh <-  rep(1,horizon-1)
poplow <-  rep(1,horizon-1)

educmed <-  rep(1,horizon-1)
educhigh <-  rep(1,horizon-1)
educlow <-  rep(1,horizon-1)

pop<-  rep(1,horizon-1)
popmed[1]<-1
pophigh[1]<-1
poplow[1]<-1

for (t in (1:(horizon-1))){
nlow[t] <- n[ilow,t]
nhigh[t] <- n[ihigh,t]
nmed[t] <- n[imed,t]

incmed[t] <- y[imed,t]
inchigh[t] <- y[ihigh,t]
inclow[t] <- y[ilow,t]

educlow[t] <- phi*eta[ilow]/(gamma-eta[ilow])*h[ilow,t]*n[ilow,t]/inclow[t]
educhigh[t] <- phi*eta[ihigh]/(gamma-eta[ihigh])*h[ihigh,t]*n[ihigh,t]/inchigh[t]
educmed[t] <- phi*eta[imed]/(gamma-eta[imed])*h[imed,t]*n[imed,t]/incmed[t]

grmed[t] <- (((1-phi*nmed[t+1])*h[imed,t+1]+a)/((1-phi*nmed[t])*h[imed,t]+a))^(1/20)
grhigh[t] <- (((1-phi*nhigh[t+1])*h[ihigh,t+1]+a)/((1-phi*nhigh[t])*h[ihigh,t]+a))^(1/20)
grlow[t] <- (((1-phi*nlow[t+1])*h[ilow,t+1]+a)/((1-phi*nlow[t])*h[ilow,t]+a))^(1/20)

pop[t]<- sum(p[1:nhouseholds,t]) 
popmed[t] <- p[imed,t]/pop[t]
pophigh[t] <- p[ihigh,t]/pop[t]
poplow[t] <- p[ilow,t]/pop[t]
}

#aggregates
navg <- rep(1,horizon-1)
incavg <- rep(1,horizon-1)

for (t in (1:(horizon-1))){
navg[t] <- sum(n[1:nhouseholds,t]*p[1:nhouseholds,t])/sum(p[1:nhouseholds,t])
incavg[t] <- sum(y[1:nhouseholds,t]*p[1:nhouseholds,t])/sum(p[1:nhouseholds,t])
}

year <- rep(1,horizon-1)
for (t in (1:(horizon-1))){
year[t] <- 1475+20*(t-1)
}

pdf("macro-sim.pdf", width = 7, height = 6)  # save to PDF

par(mfrow = c(3, 2), mar = c(4, 4, 1.5, 1), mgp = c(2, 0.7, 0))

plot(nhigh~year, type = "l", xlab = "time", ylab = "no. of children per adult",ylim=c(2,4),main="Fertility")
lines( nmed~year, lty = 2)
lines( nlow~year, lty = 4)
legend("topleft", legend = c("Q75","Q50","Q25"), lty = c(1,2,4), bty = "n")

plot((pophigh)~year, type = "l", xlab = "time", ylab = "share",ylim=c(0,.004),main="Population  in % total (6oo households)")
lines( (popmed)~year, lty = 2)
lines( (poplow)~year, lty = 4)
legend("topleft", legend = c("Q75","Q50","Q25"), lty = c(1,2,4), bty = "n")


plot((educhigh)~year, type = "l", xlab = "time", ylab = "% of income",ylim=c(.04,.085),main="Education spending in % of income")
lines( (educmed)~year, lty = 2)
lines( (educlow)~year, lty = 4)
legend("topleft", legend = c("Q75","Q50","Q25"), lty = c(1,2,4), bty = "n")



plot((grhigh)~year, type = "l", xlab = "time", ylab = "growth factor",ylim=c(1,1.006),main="Annualized growth factor of income")
lines( (grmed)~year, lty = 2)
lines( (grlow)~year, lty = 4)
legend("topleft", legend = c("Q75","Q50","Q25"), lty = c(1,2,4), bty = "n")

plot((navg)~year, type = "l", xlab = "time", ylab = "no. of children per adult",main="Aggregate fertility of academics")

plot((incavg)~year, type = "l", xlab = "time", ylab = "",main="Aggregate income of academics")


dev.off()




