##########  PS-239 - FALL 2006 #############################################
##########  SECTION 10      #############################################
# Generate data
# Draw covariates
set.seed(101)
n<-1000
X1<-rnorm(n,mean=10,sd=3)
# Draw disturbance
u<-rnorm(n,mean=0,sd=1)
# Generate Y
# Y<-10+0.005*X1++u
Y<-10+0.2*X1++u
# Null hypothesis: beta1=0
# Alternative: beta1>0

# Full sample and parameters
fsbeta0<-lm(Y~X1)$coefficients[1]
fsbeta1<-lm(Y~X1)$coefficients[2]
fse<-lm(Y~X1)$residuals
fsbeta1.sd<-summary(lm(Y~X1))$coefficients[2,1]/summary(lm(Y~X1))$coefficients[2,2]

nboot<-1000
dta<-data.frame(Y,X1)
obs<-length(Y)
set.seed(1254)
count1<-0
count2<-0
count3<-0
count4<-0
count5<-0


for (i in 1:nboot) {

  # Bootstrap method 1 ==> permutation test
  bx<-X1
  by<-Y[sample(1:obs,replace=FALSE)]
  bb<-lm(by~bx)$coefficients[2]

  if(bb>=fsbeta1) count1=count1+1

  # Bootstrap method 2 ==> similar to permutation test
  bx<-X1[sample(1:obs,replace=TRUE)]
  by<-Y[sample(1:obs,replace=TRUE)]
  bb<-lm(by~bx)$coefficients[2]

  if(bb>=fsbeta1) count2=count2+1
  
  # Bootstrap method 3 ==> use the model structure
  bx<-X1
  be<-sample(fse,replace=TRUE)
  by<-fsbeta0+be # plug in the null beta1=0
  bb<-lm(by~bx)$coefficients[2]
  
  if(bb>=fsbeta1) count3=count3+1
  
  # Bootstrap method 4 ==> studentized
  be<-sample(fse,replace=TRUE)
  by<-fsbeta0+fsbeta1*X1+be # plug in the null beta1=0
  bls<-summary(lm(by~X1))
  bb<-bls$coefficients[2,1]
  bse<-bls$coefficients[2,2]
  bbeta1.sd<-(bb-fsbeta1)/bse
  
  if(bbeta1.sd>=fsbeta1.sd) count4=count4+1
  
  # Bootstrap method 5 ==> count #(beta.boot-beta.fs>beta.fs)
  if(bb-fsbeta1>=fsbeta1) count5=count5+1

}
  
count1/nboot
count2/nboot
count3/nboot
count4/nboot
count5/nboot
  
  
# Monte Carlo
mc<-200
nboot<-1000
obs<-1000
mc.results<-matrix(nrow=mc,ncol=5)
set.seed(101)

for (m in 1:mc) {
  # Generate data

  X1<-rnorm(obs,mean=10,sd=3)
  # Draw disturbance
  u<-rnorm(obs,mean=0,sd=1)
  # Generate Y
  Y<-10+0*X1++u
  #Y<-10+0.003*X1++u

  #Full sample and parameters
  fsbeta0<-lm(Y~X1)$coefficients[1]
  fsbeta1<-lm(Y~X1)$coefficients[2]
  fse<-lm(Y~X1)$residuals
  fsbeta1.sd<-summary(lm(Y~X1))$coefficients[2,1]/summary(lm(Y~X1))$coefficients[2,2]

  dta<-data.frame(Y,X1)
  count1<-0
  count2<-0
  count3<-0
  count4<-0


  for (i in 1:nboot) {

    # Bootstrap method 1 ==> permutation test
    bx<-X1
    by<-Y[sample(1:obs,replace=FALSE)]
    bb<-lm(by~bx)$coefficients[2]

    if(bb>=fsbeta1) count1=count1+1

    # Bootstrap method 2 ==> similar to permutation test
    bx<-X1[sample(1:obs,replace=TRUE)]
    by<-Y[sample(1:obs,replace=TRUE)]
    bb<-lm(by~bx)$coefficients[2]

    if(bb>=fsbeta1) count2=count2+1

    # Bootstrap method 3 ==> use the model structure
    bx<-X1
    be<-sample(fse,replace=TRUE)
    by<-fsbeta0+be # plug in the null beta1=0
    bb<-lm(by~bx)$coefficients[2]

    if(bb>=fsbeta1) count3=count3+1

    # Bootstrap method 4 ==> studentized
    be<-sample(fse,replace=TRUE)
    by<-fsbeta0+fsbeta1*X1+be
    bls<-summary(lm(by~X1))
    bb<-bls$coefficients[2,1]
    bse<-bls$coefficients[2,2]
    bbeta1.sd<-(bb-fsbeta1)/bse

    if(bbeta1.sd>=fsbeta1.sd) count4=count4+1
    
    # Bootstrap method 5 ==> count #(beta.boot-beta.fs>beta.fs)
    if(bb-fsbeta1>=fsbeta1) count5=count5+1
    


  } # close bootstrap
  cat("Round: ",m,"\n")
  mc.results[m,1] = (count1/nboot)
  mc.results[m,2] = (count2/nboot)
  mc.results[m,3] = (count3/nboot)
  mc.results[m,4] = (count4/nboot)
  mc.results[m,5] = (count5/nboot)
} # close mc

cat("\n")
alpha = .1
cat(alpha,"NOMINAL COVERAGE RESULTS:\n")
cat("Method 1:", mean(mc.results[,1] < .01), "\n")
cat("Method 2:", mean(mc.results[,2] < .01), "\n")
cat("Method 3:", mean(mc.results[,3] < .01), "\n")
cat("Method 4:", mean(mc.results[,4] < .01), "\n")
cat("Method 5:", mean(mc.results[,5] < .01), "\n")
mean(mc.results[,1]<0.01)
mean(mc.results[,2]<0.01)
save(list = ls(all=TRUE),file="W:/Teaching/Fall2006_PolSci239/Section_Notes/Section10/BootstrapOLS.RData")
load(file="W:/Teaching/Fall2006_PolSci239/Section_Notes/Section10/BootstrapOLS")