"bootstrap" <-
function(effa,effb,effcp,n,nrep,nboot,alpha)
{

 #emp_power.r by Matt Fritz, 2005, Arizona State University
 #This program calculates power for the percentile and bias-
 #corrected bootstrap for a specified sample size.
 #To use, import bootstrap_power into R using the command
 #		source("C:/bootstrap_power.r")
 #where C:/ should be replaced with the address of the 
 #bootstrap_power.r file which will vary by computer.  Then, to 
 #run the program use the command 
 #		bootstrap(effa, effb, effcp, n, nrep, nboot,alpha)
 #inputting the values of interest, where effa=alpha path, 
 #effb=beta path, effcp=direct effect(c prime path),
 #sampsize=sample size, nrep=number of replications,
 #nboot=number of bootstrap samples, and alpha=the Type I
 #error rate.

percent <- rep(0,nrep)
bc.old <- rep(0,nrep)
numcases <- seq(1:n)

for (i in 1:nrep) {
  x <- rnorm(n)
  x2 <- rnorm(n)
  x3 <- rnorm(n)
  m <- effa * x + x2
  y <- effb * m + effcp * x + x3

  fit.mx <- lm(m~x)
  sum.mx <- summary(fit.mx)
  a <- sum.mx$coefficients[2,1]
  sea <- (sum.mx$coefficients[2,2])
	
  fit.yxm <- lm(y~x+m)
  sum.yxm <- summary(fit.yxm)
  b <- sum.yxm$coefficients[3,1]
  seb <- (sum.yxm$coefficients[3,2])

  ab <- a*b

  #Start Bootstrap Loop j= bootstrap number
  boot.percent<-rep(0,nboot)

   for (k in 1:nboot) {
    case <- sample(numcases, n, replace = TRUE)
    xnew <- rep(0,n)
    mnew <- rep(0,n)
    ynew <- rep(0,n)	
	   
    for (l in 1:n) {
     xnew[l] <- x[case[l]]
     mnew[l] <- m[case[l]]
     ynew[l] <- y[case[l]] }	

    bfit.mx<-lm(mnew~xnew)
    bsum.mx<-summary(bfit.mx)
    boot.a<-bsum.mx$coefficients[2,1]
	   
    bfit.yxm<-lm(ynew~xnew+mnew)
    bsum.yxm<-summary(bfit.yxm)
    boot.b<-bsum.yxm$coefficients[3,1]
		
    boot.percent[k] <- boot.a*boot.b
   }
  #End Bootstrap Loop	

  ##Percentile Bootstrap
  sort.boot.percent<-sort(boot.percent)
  upper<-round((1-(alpha/2))*nboot)
  lower<-round((alpha/2)*nboot)
  ifelse(lower==0, lower<-1, lower<-lower)
  uperboot<-sort.boot.percent[upper]
  lperboot<-sort.boot.percent[lower]
  ifelse(lperboot<=0 && uperboot>=0, percent[i]<-0, percent[i]<-1)

  #Bias Corrected
  bias<-rep(0,n)
  for (q in 1:n) {
   ifelse(boot.percent[q]<=ab, bias[q]<-1, bias[q]<-0) }
  z0<-mean(bias)
  zu.bias.old<-round((pnorm(2*qnorm(z0)+qnorm(1-(alpha/2))))*nboot)
  zl.bias.old<-round((pnorm(2*qnorm(z0)-qnorm(1-(alpha/2))))*nboot)
  ifelse(zl.bias.old==0, zl.bias.old<-1, zl.bias.old<-zl.bias.old)
  ubias.old<-sort.boot.percent[zu.bias.old]
  lbias.old<-sort.boot.percent[zl.bias.old]
  ifelse(lbias.old<=0 && ubias.old>=0, bc.old[i]<-0, bc.old[i]<-1)	
}

mpercent <- mean(percent)
mbc.old <- mean(bc.old)

out<-list(	'Percentile Bootstrap'=mpercent,
		'Old Bias Corrected'=mbc.old)
out

}