#-------------------------------
#     SINGLE VERSION
#-------------------------------
N=dim(labor1)[1]
set.seed(1)

## reduced.prop is a matrix that stores the proportion of targeting the original model when we reduce the sample size. The first row stores the proportions of BIC and the second row stores the proportions of AIC
reduced.prop=matrix(NA, 2, 28)
judge.BIC=rep(FALSE,100); judge.AIC=rep(FALSE,100)
for (T in seq(285, 420, 5))#
{
  for (i in 1:100)
  {
    ind=sample(N, T)
    train=labor1[ind,]
    reg.all=lm(lnearnings ~ WA+I(WA^2)+WE+kids+HA+HE+WMED+WFED+UN+CIT+AX+MTR+lnearning.H, data=train)

    train.BIC=step(reg.all,scope=list(upper=~.,lower=~1),k=log(T),trace=0)
    judge.BIC[i]=(train.BIC$terms[[3]]==m.BIC$terms[[3]])

    train.AIC=step(reg.all,scope=list(upper=~.,lower=~1),trace=0)
    judge.AIC[i]=(train.AIC$terms[[3]]==m.AIC$terms[[3]])
  }
  reduced.prop[1,T/5 - 56]=sum(judge.BIC)/100
  reduced.prop[2,T/5 - 56]=sum(judge.AIC)/100
}
reduced.prop

#---------------------------------------------
#     PARALLEL COMPUTING
#---------------------------------------------
library(snowfall)
sfInit(parallel = TRUE, cpus = 7, type="SOCK")
sfExport('reg.all','N','labor1','m.BIC','m.AIC')

inner.loop = function(T){
  reg.all = apply(matrix(1:100,100,1),1,function(i){
    ind = sample(N,T)
    train=labor1[ind,]
    reg.all=lm(lnearnings ~ WA+I(WA^2)+WE+kids+HA+HE+WMED+WFED+UN+CIT+AX+MTR+lnearning.H, data=train)

    train.BIC=step(reg.all,scope=list(upper=~.,lower=~1),k=log(T),trace=0)
    judge.BIC=(train.BIC$terms[[3]]==m.BIC$terms[[3]])

    train.AIC=step(reg.all,scope=list(upper=~.,lower=~1),trace=0)
    judge.AIC=(train.AIC$terms[[3]]==m.AIC$terms[[3]])
    return(list(judge.BIC,judge.AIC))

  })
  reg.all = matrix(unlist(reg.all),2,100)
  return(list(mean(reg.all[1,]),mean(reg.all[2,])))
}

reduced.prop <- sfSapply(seq(285,420,5), inner.loop)