##delta method and bootstrap for minimum of quadratic function

rm(list= ls())


pdf("delta.pdf")
set.seed(89777)
n = 100
X = rnorm(n)
Y = 2*X^2 + rnorm(n)

X2 = X^2
out = lm(Y ~ X + X2)
## print(summary(out))

V = vcov(out)
beta = out$coef
psi = -.5*beta[2]/beta[3]
g = c(0,-.5/beta[3],.5*beta[2]/beta[3]^2)
g = as.matrix(g)
b = sqrt(t(g) %*% V %*% g)
## print(psi)

grid = seq(-3,3,length=100)
mu = beta[1] + beta[2]*grid + beta[3]*grid^2
plot(X,Y,xlab="",ylab="", pch=21,cex=2.,bg="pink",cex.axis=1.4)
lines(grid,mu,col="purple",lwd=5)
abline(v=psi,col="red",lwd=3)
title(xlab="X", line=2.5, cex.lab=1.8)
title(ylab="Y", line=2.3, cex.lab=1.8)

## print(b)
## [1,] 0.0283808  
## print(c(psi - 1.96*b,psi + 1.96*b))

dev.off()

pdf("boot.pdf")

hat = psi
B = 10000
psi = rep(0,B)
for(i in 1:B){
     I = sample(1:n,size=n,replace=TRUE)
     XX = X[I]
     YY = Y[I]
     XX2 = XX^2
     out = lm(YY ~ XX + XX2)
     beta = out$coef
     psi[i] = -.5*beta[2]/beta[3]
     }

hist(psi,main=NULL,col="thistle",border="purple",,xlab="",
         ylab="",cex.axis=1.4)
title(xlab=expression(paste("Histogram of ", psi)), line=2.7, cex.lab=1.8)
title(ylab="Frequency", line=2.3, cex.lab=1.8)

b = sqrt(var(psi))
print(b)
## [1] 0.02223679
## print(c(hat - 1.96*b,hat + 1.96*b))
dev.off()