set.seed(1234)
library ('arm')
lalondedat <- data('lalonde')
datacont <- subset(lalonde, treat == '0')

contmodel <- lm(datacont$re78 ~ datacont$age+datacont$educ+datacont$re74+datacont$re75+(datacont$educ*datacont$re74)+(datacont$educ*datacont$re75)+(datacont$age*datacont$re74)+(datacont$age*datacont$re75)+(datacont$re74*datacont$re75))
summary(contmodel)
simulations <- sim(contmodel, n.sims = 10000)

storage_array <- rep(0:10000)
result_matrix <- matrix(rep(1:10000*39), nrow = 39, ncol = 10000)

for (i in 17:55){
  matrix_1 <- matrix(c(1,i,quantile(datacont$educ, 0.9),quantile(datacont$re74, 0.9),quantile(datacont$re75, 0.9),quantile(datacont$educ, 0.9)*quantile(datacont$re74, 0.9),quantile(datacont$educ, 0.9)*quantile(datacont$re75, 0.9),i*quantile(datacont$re74, 0.9),i*quantile(datacont$re75, 0.9),quantile(datacont$re74, 0.9)*quantile(datacont$re75, 0.9)))
  storage_array <- simulations@coef %*% matrix_1 + rnorm(10000,0,simulations@sigma)
  result_matrix[i -16,] <- storage_array
}

q_lower_quantile <- rep(0:39)
q_upper_quantile <- rep(0:39)
age <- rep(0:39)

for (i in 1:39){
  q_lower_quantile[i] <- quantile(result_matrix[i,], 0.025)
  q_upper_quantile[i] <- quantile(result_matrix[i,], 0.975)
  age[i] <- i+16
}

pointestimates <- data.frame(age,q_lower_quantile,q_upper_quantile)
write.csv(pointestimates, file = "Quantile Point Estimates.csv")

plot(c(1:100), c(1:100), type = 'n', ylim = c(-7000, 20000), xlim = c(17,55),
     main = "Scatter Plot While Holding Predictors at 90% Quantile", xlab = "Age", ylab = "Real Earning Income 1978")
for (i in 17:55) {
  segments(
    x0 = i,
    y0= q_lower_quantile[i - 16],
    x1 = i,
    y1= q_upper_quantile[i - 16])
}