rm()

if(!require(nleqslv)){
    install.packages("nleqslv")
    library(nleqslv)}
if(!require(DEoptim)){
    install.packages("DEoptim")
    library(DEoptim)}
if(!require(latex2exp)){
    install.packages("latex2exp")
    library(latex2exp)}

if(!require(doParallel)){
    install.packages("doParallel")
    library(doParallel)}
if(!require(foreach)){
    install.packages("foreach")
    library(foreach)}
if(!require(doRNG)){
    install.packages("doRNG")
    library(doRNG)}


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


sink("report-indirect-inference.txt")

# four parameters, psi, cbar/a, gamma, eta,  are determining by perfect identification with 4 moments 
# the level of a is set to reach hbar after 4 periods, i.e. hbar = 1.3*h0

#the distribution of eta is assumed beta. Symmetry is assumed to the two parameters of the beta are equal
# eta has a minimum value which cannot be below the one generation stagant human capital for that family

# the three parameters are then estimated matching the coefficient of the regression at 7 point in time
#  restricting the simulation to match perfectly the moment at time 4

h0 <- 1   # initial condition
phi <- 1/11  # 11 sons is a maximum

gr <- (1+0.001256)^20   # from the file "growth from regression.txt"

print(c("implied growth rate (from regression)",gr))


#vector x is : psi, cbar/a, gamma, eta, a
base <- function(x) {
y <- numeric(5)
y[1] <- (x[3]-x[4])-(1+x[3])*phi   # one child when h is infinite
y[2] <- (x[3]-x[4])*x[2] - 3.316*phi*((1+x[3])*x[2]-1) # 3.316 (mean of Sibshipsize) children at h = bar h
y[3] <- phi*x[4]-(1-phi)*(x[3]-x[4])*0.073  # education spending 7.3% of gdp when h is infinite
y[4] <- x[1]*phi*x[4] - gr *(x[3]-x[4])  # growth rate from one generation to the next of 2.63%
            # = 1.0013^20
#y[5] <- 40*x[5]-15*(x[5]+h0*(1-3*phi))      # Petty's national account, share of labor in GDP is 25/40, 
#                                            #impose at h=h0
y[5] <- ((1+x[3])*x[2]-1)*x[5]-h0*1.30             # hbar is 30% higher than h0 - this is a way to peg the level of a
y}

xstart <- c(15,1.5,.3,.1,1) #initial guess
sol<- nleqslv(xstart, base)# solving the non-linear system
par<-sol$x
print("Attained value of function:")
print(round(sol$fvec,5))

print("Parameters")
cbind(c("psi","cbar/a","gamma","eta","a","cbar","hbar"),
                   c(round(par,5),round(par[2]*par[5],5),round( ((1+par[3])*par[2]-1)*par[5],5)))


psi<-par[1]
gamma<-par[3]
a<-par[5]
cbar<-par[2]*par[5]
etamean <- par[4]

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

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



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

print(paste("Bounds on eta:",round(eta1,5),round(eta2,5)))

#---------------------------------------------------------------------------------------
#this functions will have parameters as input and will deliver a vector of coefficients
#from a rolling regression on artificial data

simulation <- function(etamin, etamax, distpar1, distpar2, noise){

  nh <- nhouseholds
  T  <- horizon

  # eta: quantiles of Beta, evenly spaced across households
  eta <- etamin + (etamax - etamin) * qbeta((0:(nh-1))/(nh-1), distpar1, distpar2)

  # Human capital h: geometric growth per household across t
  # h[i,t] = h0 * ((psi*phi*eta[i]/(gamma-eta[i]))^(t-1))
  r    <- (psi * phi * eta) / (gamma - eta)            # length nh
  tpow <- 0:(T-1)
  # Use exp(outer(log(r), tpow)) for numerical stability and speed
  h <- h0 * exp(outer(log(r), tpow))                   # nh x T

  # Fertility n with bounds, fully vectorized
  # b1 = (gamma-eta)/gamma * (h+a-cbar)/(phi*h)
  # b2 = (gamma-eta)/(1+gamma) * (h+a)/(phi*h)
  term <- (gamma - eta)                                # length nh
  H    <- h
  b1   <- (term / gamma)     * (H + a - cbar) / (phi * H)
  b2   <- (term / (1 + gamma)) * (H + a) / (phi * H)
  n    <- pmax(0, pmin(b1, b2))

  # Add noise in one shot
  nint <- n + matrix(rnorm(nh * T, 0, noise), nrow = nh, ncol = T)
  nint <- pmax(nint, 0)   # forces everything to be ≥ 0

  # Regress asinh(h[, t+1]) on nint[, t] for t = 1..T-1
  # slope = cov(x,y)/var(x) computed column-wise (with intercept)
  X  <- nint[, 1:(T-1), drop = FALSE]                  # nh x (T-1)
  Y  <- asinh(h[, 2:T, drop = FALSE])                  # nh x (T-1)

  mx <- colMeans(X)
  my <- colMeans(Y)
  Xc <- sweep(X, 2, mx, FUN = "-")
  Yc <- sweep(Y, 2, my, FUN = "-")

  coef <- colSums(Xc * Yc) / colSums(Xc * Xc)

  return(coef)
}
#---------------------------------------------------------------------------------------


obj <- function(param){
dist1<-param[1]
dist2<-param[1]
noise<-param[2]
simf <- simulation(eta1,eta2,dist1,dist2,noise)[p:(p+6)]
omega <- (simf[1]-mom[1])^2+(simf[2]-mom[2])^2+(simf[3]-mom[3])^2+(simf[4]-mom[4])^2+
                   (simf[5]-mom[5])^2+(simf[6]-mom[6])^2+(simf[7]-mom[7])^2
return(omega)
}



#from the file "rolling.txt"
# assuming the following points in time:  1635, 1655, 1675, 1695, 1715, 1735, 1755
# using datemed (and cheating for the last period which has datemed=1748)

#[2,]		1514	1636	1692	0.057		0.013		0.1		694	615	3.439481	2.200275
#[7,]		1585	1657	1701	0.058		0.014		0.102		694	588	3.449568	2.223362
#[12,]	1615	1676	1708	0.056		0.01		0.102		693	554	3.460317	2.200945
#[18,]	1641	1694	1715	0.004		-0.049	0.056		694	519	3.357349	2.177339
#[33,]	1694	1715	1729	-0.061	-0.12		-0.001	694	458	3.278098	2.066474
#[55,]	1723	1735	1747	-0.083	-0.135	-0.031	694	443	3.360231	1.997814
#[70,]	1736	1748	1779	-0.056	-0.109	-0.002	694	496	3.345821	2.055859


mommean <- c(0.057,0.058,0.056,0.004,-0.061,-0.083,-0.056)
low05 <- c(0.013,0.014,0.01,-0.049,-0.12,-0.135,-0.109)
high05 <- c(0.1,0.102,0.102,0.056,-0.001,-0.031,-0.002)
momse <- (mommean-low05)/qt(.975,599)


#estimation with moments set at their mean

set.seed(10)  #set same seed to get the same result each time estimation is run 
mom<-mommean
best <- DEoptim(obj,lower=c(1,0.1),upper=c(5,2))
vecsol<-best[[1]]$bestmem


w1<-p #window
w2<-p+6
print('parameters')
round(vecsol,5)
print('fitted moments')
matrix(round(simulation(eta1,eta2,vecsol[1],vecsol[1],vecsol[2])[w1:w2],5))


###############################################################################################################
#bootstrapping (parallel computations)

cl <- makeCluster(max(1, parallel::detectCores() - 1))
on.exit(stopCluster(cl), add = TRUE)
registerDoParallel(cl)
registerDoRNG(10)  # reproducible parallel RNG

nboot <- 100

vecsolboot <- foreach(i = 1:nboot, .combine = rbind, .packages = "DEoptim") %dopar% {
  mom <- rnorm(7, mean = mommean, sd = momse)  #drawing moments

  obj_i <- function(theta) obj(theta)

  best <- DEoptim(
    fn = obj_i,
    lower = c(1, 0.1),
    upper = c(5, 2),
    control = DEoptim.control(itermax = 150, trace = FALSE)
  )
  unname(best$optim$bestmem)
}

vecsolboot <- as.data.frame(vecsolboot)
names(vecsolboot) <- c("zeta", "sigma")
paste("zeta: ",round(mean(vecsolboot$zeta),5)," (",round(sd(vecsolboot$zeta),5),")",sep="")
paste("sigma: ",round(mean(vecsolboot$sigma),5)," (",round(sd(vecsolboot$sigma),5),")",sep="")

sink()



#graph

pdf(file="simulation.pdf", width = 11, height = 7)
sims <- round(simulation(eta1,eta2,vecsol[1],vecsol[1],vecsol[2]),2)
year<-c(1475,1495,1515,1535,1555,1575,1595,1615,1635, 1655, 1675, 1695, 1715, 1735, 1755,1775,1795,1815,1835)
plot(sims~year,type = "l", lty = 1,lwd=2,ylim=c(-.5,.35),xlim=c(1500,1850),cex.axis=1.3,cex.lab=1.5,ylab=TeX("Regression coef. of $n_t$ on $h_{t+1}$"))

lines(mommean~year[p:(p+6)],lty = 3,lwd=2)
points(mommean~year[p:(p+6)],cex=1.5,lwd=2)
polygon(x=c(year[p:(p+6)],rev(year[p:(p+6)])),y=c(high05,rev(low05)), col =  adjustcolor("gray", alpha.f = 0.20), border = NA) #confidence interval

sims2 <- round(simulation(eta1,eta2,vecsol[1],vecsol[1],0),2)
lines(sims2~year,lty = 2)

legend(1500, -.3, legend=c(TeX("Regression coefficients $\\hat{\\alpha}_1$ (data)"),
                           TeX("Regression coefficients $\\beta_t(\\hat{\\zeta},\\hat{\\sigma})$ (benchmark simulation)"),
                           TeX("Regression coefficients $\\beta_t(\\hat{\\zeta},0)$ (simulation without measurement error)")),
                  lty=c(3,1,2),pch=c(1,NA,NA),lwd=c(2,2,1),bty = "n",cex=c(1.2,1.2,1.2))
dev.off()



pdf(file="simulation2.pdf", width = 11, height = 7)
sims <- round(simulation(eta1,eta2,vecsol[1],vecsol[1],vecsol[2]),2)

plot(sims~year,type = "l", lty = 1,lwd=2,ylim=c(-.35,.35),,cex.axis=1.3,cex.lab=1.5,ylab=TeX("Regression coef. of $n_t$ on $h_{t+1}$"))

for (i in (1:nboot)){
sims <- round(simulation(eta1,eta2,vecsolboot[i,1],vecsolboot[i,1],vecsolboot[i,2]),2)
lines(sims~year,type = "l",col="lightgray",add=T)
}


dev.off()



