# we will create a model for this simulation

beta0 <- 1.5
beta1 <- 2.5
beta2 <- -1.0

x <- 1:20
y <- beta0 + beta1*x + beta2*x^2 + rnorm(20,0,20)
x2 <- x^2

fit <- lm(y ~ x + x2)

# check the SE
summary(fit)

#check the variance-covariance matrix
summary(fit)$cov.unscaled*summary(fit)$s^2

# now lets do a parameteric bootstrap

n.boots <- 10000

new.betas <- matrix(0,n.boots,3)
for (i in 1:n.boots){
  y.new <- cbind(1,x,x2) %*% fit$coef + rnorm(20,0,summary(fit)$s)
  new.betas[i,] <- lm(y.new ~ x + x2)$coef
}

# what are the bootstrap-estimates of the SE

sd(new.betas[,1])
sd(new.betas[,2])
sd(new.betas[,3])


# a bootstrap estimate of the variance-covariance matrix

sum.cov <- matrix(0,dim(new.betas)[2],dim(new.betas)[2])
for (i in 1:n.boots){
  sum.cov <- sum.cov + (new.betas[i,] - fit$coef) %o%(new.betas[i,]-fit$coef)
}

sum.cov/n.boots


# the non-parameteric bootstrap

for (i in 1:n.boots){
  y.new <- cbind(1,x,x2) %*% fit$coef +sample(resid(fit),replace=T)
  new.betas[i,] <- lm(y.new ~ x + x2)$coef
}

# what are the bootstrap-estimates of the SE
sd(new.betas[,1])
sd(new.betas[,2])
sd(new.betas[,3])

# a bootstrap estimate of the variance-covariance matrix

sum.cov <- matrix(0,dim(new.betas)[2],dim(new.betas)[2])
for (i in 1:n.boots){
  sum.cov <- sum.cov + (new.betas[i,] - fit$coef) %o%(new.betas[i,]-fit$coef)
}

sum.cov/n.boots