R version - R 3.5.1
Windows 10


getwd()
setwd("C:/ProjectPhD")
getwd()
data <- read.csv(file="DEAdata1.csv", header=T)

names(data)
library(lpSolveAPI)
library(ucminf)
library(Benchmarking)
x = matrix(c(data$x1,data$x2,data$x3,data$x4),nrow = 2989, ncol = 4)
yd = matrix(c(data$y1,data$y),nrow = 2989, ncol = 2)
yu = matrix(c(data$Yu),nrow = 2989,ncol = 1)
te <- data$TE

#attach(data)
#xa=t(cbind(x1,x2,x3,x4))
# yb=t(cbind(y1,y2))

# is.na(te)-> detect_NA
# fix(data)


# Step 1
dhat <- dea.direct(x,yd,DIRECT=x, RTS="vrs", ORIENTATION=1)
# Step 2
n <- dim(x)[1]
h <- density(c(dhat,2-dhat))$bw # bw: bandwidth      
dhat2 <- c(dhat,2-dhat)
B <- 2000
Bm <- matrix(NA,B,n)

# Step 8
for (i in 1:B)
{
# Step 3
dstar <- sample(dhat2,n,replace=T)
eps <- rnorm(n)
dt <- dstar+h*eps
# Step 4
v <- var(dhat2)
dts <- 1+1/sqrt(1+h^2/v)*(dt-1)
# Step 5
dtss <- ifelse(dts > 1,dts,2-dts)
# Step 6
xs <- x*dhat/dtss
# Step 7
#z <- c(1,1,1,1)
Bm[i,] <- dea.direct(x,yd,DIRECT=x, XREF = xs, YREF =yd)
}


# Step 9
dhat.m <- colMeans(Bm)

dhat.bias <- dhat.m-dhat
# Step 10
dhat.bc <- dhat-dhat.bias
# Step 11
Bm.c <- t(apply(-Bm,1,function(z) z+2*dhat))
ci.low <- apply(Bm.c,2,quantile,0.05)
ci.high <- apply(Bm.c,2,quantile,0.95)
# Results
tab <- round(rbind(dhat,dhat.m,dhat.bias,
dhat.bc,ci.low,ci.high),2)



#**************************** Error**********************
Bm[i,] <- dea.direct(x,yd,DIRECT=x, XREF = xs, YREF =yd)

[1] Unit 54 is not in the technology set.  Status = 2
[1] Unit 147 is not in the technology set.  Status = 2
[1] Unit 362 is not in the technology set.  Status = 2
[1] Unit 404 is not in the technology set.  Status = 2
[1] Error in solving for unit 514 :  Status = 5
[1] Unit 736 is not in the technology set.  Status = 2
[1] Unit 737 is not in the technology set.  Status = 2
[1] Unit 738 is not in the technology set.  Status = 2
[1] Unit 739 is not in the technology set.  Status = 2
[1] Unit 740 is not in the technology set.  Status = 2
[1] Error in solving for unit 764 :  Status = 5
[1] Error in solving for unit 936 :  Status = 5
[1] Error in solving for unit 1027 :  Status = 5
[1] Unit 1627 is not in the technology set.  Status = 2
[1] Error in solving for unit 2575 :  Status = 5
Error in Bm[i, ] <- dea.direct(x, yd, DIRECT = x, XREF = xs, YREF = yd) : 
  number of items to replace is not a multiple of replacement length
                             or
 
Bm[i,] <- dea.direct(x,yd,DIRECT=x,XREF=xs,YREF=yd)$eff

Error in dea(x, yd, RTS = RTS, ORIENTATION = ORIENTATION, XREF = XREF,  : 
  Number of units must be the same in x and yd