################
### Preparation
################

rm(list = ls())

require(foreign)
require(Hmisc)
require(combinat)
require(simex)

EULIS_EUP <- read.dta("EULIS_EUP.dta", convert.factors=FALSE)

nman <- nrow(EULIS_EUP)
nrepl <- 2000

actorBSn <- actorBSnRand <- array(as.matrix(EULIS_EUP[,17:56]),
                 c(nman, 40, nrepl+1),
                 dimnames=list(1:nman, names(EULIS_EUP[,17:56]), 0:nrepl))
				 
###############################
## Bootstrap (multinomial draw)
###############################
require(combinat)
n <- apply(actorBSn[1:nrow(actorBSn),,1],1,sum)
p <- actorBSn[,,1]/n
cat("Resample number: ")
 for (i in 1:nrepl) {
  cat(i," ")
  actorBSn[,,i] <- rmultinomial(n, p)
}
cat("\n")

############
## Positions
############

pnsBS <-
  ((log(actorBSn[,"p101s_107",]+0.5)-log(actorBSn[,"p101n_109",]+0.5)) +
   (log(actorBSn[,"p102s_108",]+0.5)-log(actorBSn[,"p102n_110",]+0.5)) +
   (log(actorBSn[,"p103s_602",]+0.5)-log(actorBSn[,"p103n_601",]+0.5))) /3
pnsMean_rep <- apply(pnsBS,1,mean)
pnsSE_rep <- apply(pnsBS,1,sd)


plrBS <-
 ((log(actorBSn[,"p201r_104",]+0.5)-log(actorBSn[,"p201l_105",]+actorBSn[,"p201l_106",]+0.5)) +
  (log(actorBSn[,"p202r_605",]+0.5)-log(actorBSn[,"p202l_201",]+actorBSn[,"p202l_202",]+0.5)) +
  (log(actorBSn[,"p203r_303",]+0.5)-log(actorBSn[,"p203l_404",]+actorBSn[,"p203l_405",]+0.5)) +
  (log(actorBSn[,"p204r_401",]+0.5)-log(actorBSn[,"p204l_412",]+actorBSn[,"p204l_413",]+0.5)) +
  (log(actorBSn[,"p205r_402",]+0.5)-log(actorBSn[,"p205l_403",]+0.5)) +
  (log(actorBSn[,"p206r_407",]+0.5)-log(actorBSn[,"p206l_406",]+0.5)) +
  (log(actorBSn[,"p207r_414",]+0.5)-log(actorBSn[,"p207l_409",]+0.5)) +
  (log(actorBSn[,"p208r_410",]+0.5)-log(actorBSn[,"p208l_416",]+actorBSn[,"p208l_501",]+0.5)) +
  (log(actorBSn[,"p209r_505",]+0.5)-log(actorBSn[,"p209l_503",]+actorBSn[,"p209l_504",]+0.5)) +
  (log(actorBSn[,"p210r_603",]+0.5)-log(actorBSn[,"p210l_604",]+0.5)) +
  (log(actorBSn[,"p211r_608",]+0.5)-log(actorBSn[,"p211l_607",]+0.5)) +
  (log(actorBSn[,"p212r_702",]+0.5)-log(actorBSn[,"p212l_701",]+0.5)) +
  (log(actorBSn[,"p213r_704",]+0.5)-log(actorBSn[,"p213l_705",]+0.5))) /13
plrMean_rep <- apply(plrBS,1,mean)
plrSE_rep <- apply(plrBS,1,sd)


eulisnsBS <-
  (((log(actorBSn[,"p101s_107",]+0.5)-log(actorBSn[,"p101n_109",]+0.5))*EULIS_EUP$s101) +
   ((log(actorBSn[,"p102s_108",]+0.5)-log(actorBSn[,"p102n_110",]+0.5))*EULIS_EUP$s102) +
   ((log(actorBSn[,"p103s_602",]+0.5)-log(actorBSn[,"p103n_601",]+0.5))*EULIS_EUP$s103)) /100
eulisnsMean_rep <- apply(eulisnsBS,1,mean)
eulisnsSE_rep <- apply(eulisnsBS,1,sd)


eulislrBS <-
 (((log(actorBSn[,"p201r_104",]+0.5)-log(actorBSn[,"p201l_105",]+actorBSn[,"p201l_106",]+0.5))*EULIS_EUP$s201) +
  ((log(actorBSn[,"p202r_605",]+0.5)-log(actorBSn[,"p202l_201",]+actorBSn[,"p202l_202",]+0.5))*EULIS_EUP$s202) +
  ((log(actorBSn[,"p203r_303",]+0.5)-log(actorBSn[,"p203l_404",]+actorBSn[,"p203l_405",]+0.5))*EULIS_EUP$s203) +
  ((log(actorBSn[,"p204r_401",]+0.5)-log(actorBSn[,"p204l_412",]+actorBSn[,"p204l_413",]+0.5))*EULIS_EUP$s204) +
  ((log(actorBSn[,"p205r_402",]+0.5)-log(actorBSn[,"p205l_403",]+0.5))*EULIS_EUP$s205) +
  ((log(actorBSn[,"p206r_407",]+0.5)-log(actorBSn[,"p206l_406",]+0.5))*EULIS_EUP$s206) +
  ((log(actorBSn[,"p207r_414",]+0.5)-log(actorBSn[,"p207l_409",]+0.5))*EULIS_EUP$s207) +
  ((log(actorBSn[,"p208r_410",]+0.5)-log(actorBSn[,"p208l_416",]+actorBSn[,"p208l_501",]+0.5))*EULIS_EUP$s208) +
  ((log(actorBSn[,"p209r_505",]+0.5)-log(actorBSn[,"p209l_503",]+actorBSn[,"p209l_504",]+0.5))*EULIS_EUP$s209) +
  ((log(actorBSn[,"p210r_603",]+0.5)-log(actorBSn[,"p210l_604",]+0.5))*EULIS_EUP$s210) +
  ((log(actorBSn[,"p211r_608",]+0.5)-log(actorBSn[,"p211l_607",]+0.5))*EULIS_EUP$s211) +
  ((log(actorBSn[,"p212r_702",]+0.5)-log(actorBSn[,"p212l_701",]+0.5))*EULIS_EUP$s212) +
  ((log(actorBSn[,"p213r_704",]+0.5)-log(actorBSn[,"p213l_705",]+0.5))*EULIS_EUP$s213)) /100
eulislrMean_rep <- apply(eulislrBS,1,mean)
eulislrSE_rep <- apply(eulislrBS,1,sd)


##########################
## Write out Stata dataset
##########################
save.image(file="EULIS_EUP_rep.Rdata", safe=FALSE)
EULIS_EUPbs <- data.frame(cbind(EULIS_EUP[,1:134],
							   pnsMean_rep,
							   pnsSE_rep,
							   plrMean_rep,
                               plrSE_rep,
							   eulisnsMean_rep,
							   eulisnsSE_rep,
							   eulislrMean_rep,
							   eulislrSE_rep))
write.dta(EULIS_EUPbs, "EULIS_EUP_rep.dta", version=7)