y=matrix(c(3,12,5,4,17),byrow=T,5,1)
y
x3=matrix(c(3,2,5,1,6),byrow=T,5,1)
x2=matrix(c(4,2,6,2,7),byrow=T,5,1)
model=lm(y~x2+x3)
summary(model)
u=resid(model)
s2=(t(u)%*%u)/(5-3)
s2
bigb=999
numb=2345
nobs=5
b=model$coef
bb=b #first row of bb has regression coefficients from original data
i=1
jj=c(0,0,0)
while (i <=bigb) {
set.seed(numb+i)
#now create y* data for parametric boot using s2= (residSS/df) and normal errors
ystar=b[1]+ b[2]*x2 +b[3]*x3 + rnorm(nobs,0,s2)
modelstar=lm(ystar~x2+x3)
bbs=modelstar$coeff
for (j in 1:3) if (bbs[j]>b[j])  jj[j]=jj[j]+1
bb=rbind(bb, modelstar$coeff) #insert rows to bb
# number of times simulated coefficient exceeds the original
i=i+1}
print(jj)
#now P values for 3 coefficients from parametric bootstrap
print(jj/bigb)


#Now use the procedure in R
library(simpleboot)
set.seed(34)
yx=matrix(rnorm(120,5,2),40,3)
y=yx[,1]
x2=yx[,2]
x3=yx[,3]
model=lm(y~x2+x3)
summary(model)
bigr=999
boot.model = lm.boot(model, R=999)
summary(boot.model)
b1star =10  #dummy starting value for the object
for (i in 1:bigr) b1star[i] = boot.model$boot.list[[i]]$coef[1]  #correctly extracts first regr coeff.
b1s=sort(b1star)

b2star =10  #dummy starting value for the object
for (i in 1:bigr) b2star[i] = boot.model$boot.list[[i]]$coef[2]  #correctly extracts second regr coeff.
b2s=sort(b2star)

b3star =10  #dummy starting value for the object
for (i in 1:bigr) b3star[i] = boot.model$boot.list[[i]]$coef[3]  #correctly extracts third regr coeff.
b3s=sort(b3star)

bootans=cbind(b1s,b2s,b3s)
bootans[25,]#lower limit of 95 conf interval for regr coefficients
bootans[975,]#upper limit