#  written by O. Badunenko. Summer of '07

library(FEAR)
rts.test = function(x, y, alpha = 0.05, orient = 2, B = 1000){
# 
x = as.matrix(x); if(ncol(x) == 1){ x = t(x) }
y = as.matrix(y); if(ncol(y) == 1){ y = t(y) }
# 
n = ncol(x)
p = nrow(x)
q = nrow(y)

# 
# CRS Scores
# 
dhat.c = dea(XOBS = x,YOBS = y, RTS = 3, ORIENTATION = orient)
d1 = ifelse(is.na(dhat.c), 0, 1)
dhat.c = dhat.c[d1 == 1]
x = x[,d1 == 1]
y = y[,d1 == 1]
x = as.matrix(x); if(ncol(x) == 1){ x = t(x) }
y = as.matrix(y); if(ncol(y) == 1){ y = t(y) }
if (sum(d1) < n){
    cat("\n\n", "   CRS distance function estimates cannot be computed for observation(s):", "\n")
    cat("   ", c(1:length(d1))[d1 == 0], "\n\n\n")
}
n = ncol(x)
# 
# VRS Scores
# 
dhat.v = dea(XOBS = x,YOBS = y, RTS = 1, ORIENTATION = orient)
d2 = ifelse(is.na(dhat.v), 0, 1)
dhat.v = dhat.v[d2 == 1]
dhat.c = dhat.c[d2 == 1]
x = x[,d2 == 1]
y = y[,d2 == 1]
x = as.matrix(x); if(ncol(x) == 1){ x = t(x) }
y = as.matrix(y); if(ncol(y) == 1){ y = t(y) }
if (sum(d2) < n){
    cat("\n\n", "   VRS distance function estimates cannot be computed for observation(s):", "\n")
    cat("   ", c(1:length(d2))[d2 == 0], "\n\n\n")
}
n = ncol(x)
# 
# NIRS Scores
# 
dhat.n = dea(XOBS = x,YOBS = y, RTS = 2, ORIENTATION = orient)
d3 = ifelse(is.na(dhat.n), 0, 1)
dhat.n = dhat.n[d3 == 1]
dhat.v = dhat.v[d3 == 1]
dhat.c = dhat.c[d3 == 1]
x = x[,d3 == 1]
y = y[,d3 == 1]
x = as.matrix(x); if(ncol(x) == 1){ x = t(x) }
y = as.matrix(y); if(ncol(y) == 1){ y = t(y) }
if (sum(d3) < n){
    cat("\n\n", "   NIRS distance function estimates cannot be computed for observation(s):", "\n")
    cat("   ", c(1:length(d3))[d3 == 0], "\n\n\n")
}
n = ncol(x)
# 
# calculating bandwidthes
# 
h.c = eff.bw(dhat.c, method = 2)
h.n = eff.bw(dhat.n, method = 2)

se.c = dhat.c/dhat.v
se.mean.c = mean(dhat.c)/mean(dhat.v)

x.c.star = matrix(nrow = p,ncol = n)
y.c.star = matrix(nrow = q,ncol = n)

dhat.c.star = matrix(nrow = n,ncol = B)
dhat.c.v.star = matrix(nrow = n,ncol = B)

se.star.c = matrix(nrow = n,ncol = B)
se.star.mean.c = matrix(nrow = 1,ncol = B)
q2 = rep(0,B)

for(qq in 1:B){
    theta.c.star = dea.resample(dhat.c,bw = h.c)
    ratio = matrix(theta.c.star/dhat.c,nrow = 1)
    if(orient == 1){
        # input oriented
        y.c.star = y
        for(ww in 1:p){
        x.c.star[ww,]=x[ww,]*ratio
        }
    }
    else{
        # output oriented
        x.c.star = x
        for(ww in 1:q){
        y.c.star[ww,]=y[ww,]*ratio
        }       
    }
    dhat.c.star[,qq] = dea(XOBS = x,YOBS = y,
                XREF = x.c.star,YREF = y.c.star,
                RTS = 3, ORIENTATION = orient)
    dhat.c.v.star[,qq] = dea(XOBS = x,YOBS = y,
                XREF = x.c.star,YREF = y.c.star,
                RTS = 1, ORIENTATION = orient)
    se.star.c[,qq] = dhat.c.star[,qq] / dhat.c.v.star[,qq]
    se.star.mean.c[1,qq] = mean(na.omit(dhat.c.star[,qq])) / mean(na.omit(dhat.c.v.star[,qq]))
    q1 = ifelse(is.na(dhat.c.v.star[,qq]),1,0); q2[qq] = sum(q1)
}

# testing for global RTS
# real number of replications

B.c = length(na.omit(se.star.mean.c))
p.value.global.c = ifelse(orient == 1,
                sum(na.omit(se.star.mean.c) >= se.mean.c)/B.c,
                sum(na.omit(se.star.mean.c) <= se.mean.c)/B.c)

# testing for CRS for each DMU; H0: CRS == VRS

p.value.c = matrix(nrow = 1,ncol = n)
for(qqq in 1:n){
    B.c = length(na.omit(se.star.c[qqq,]))
    p.value.c[qqq] = ifelse(orient == 2,
                sum(na.omit(se.star.c[qqq,]) <= se.c[qqq])/B.c,
                sum(na.omit(se.star.c[qqq,]) >= se.c[qqq])/B.c)
}

# end CRS vs. VRS

#=====================================================================

# begin NIRS vs. VRS

se.n = dhat.n/dhat.v
se.mean.n = mean(dhat.n)/mean(dhat.v)

x.n.star = matrix(nrow = p,ncol = n)
y.n.star = matrix(nrow = q,ncol = n)

dhat.n.star = matrix(nrow = n,ncol = B)
dhat.n.v.star = matrix(nrow = n,ncol = B)

se.star.n = matrix(nrow = n,ncol = B)
se.star.mean.n = matrix(nrow = 1,ncol = B)
q4 = rep(0,B)

for(qq in 1:B){
    theta.n.star = dea.resample(dhat.n,bw = h.n)
    ratio = matrix(theta.n.star/dhat.n,nrow = 1)
    if(orient == 1){
        # input oriented
        y.n.star = y
        for(ww in 1:p){
        x.n.star[ww,]=x[ww,]*ratio
        }
    }
    else{
        # output oriented
        x.n.star = x
        for(ww in 1:q){
        y.n.star[ww,]=y[ww,]*ratio
        }       
    }
    dhat.n.star[,qq] = dea(XOBS = x,YOBS = y,
                XREF = x.n.star,YREF = y.n.star,
                RTS = 2, ORIENTATION = orient)
    dhat.n.v.star[,qq] = dea(XOBS = x,YOBS = y,
                XREF = x.n.star,YREF = y.n.star,
                RTS = 1, ORIENTATION = orient)
    se.star.n[,qq] = dhat.n.star[,qq] / dhat.n.v.star[,qq]
    se.star.mean.n[1,qq] = mean(na.omit(dhat.n.star[,qq])) / mean(na.omit(dhat.n.v.star[,qq]))
    q3 = ifelse(is.na(dhat.n.v.star[,qq]),1,0); q4[qq] = sum(q3)
}

# testing for global RTS: NIRS
# real number of replications

B.n = length(na.omit(se.star.mean.n))
p.value.global.n = ifelse(orient == 1,
                sum(na.omit(se.star.mean.n) >= se.mean.n)/B.n,
                sum(na.omit(se.star.mean.n) <= se.mean.n)/B.n)

# testing for NIRS for each DMU; H0: NIRS == VRS

p.value.n = matrix(nrow = 1,ncol = n)
for(qqq in 1:n){
    B.n = length(na.omit(se.star.n[qqq,]))
    p.value.n[qqq] = ifelse(orient == 2,
                sum(na.omit(se.star.n[qqq,]) <= se.n[qqq])/B.n,
                sum(na.omit(se.star.n[qqq,]) >= se.n[qqq])/B.n)
    
}

# end NIRS vs. VRS

#=====================================================================

rts = 0 ; #alpha = 0.05 ; # B = 10 ; orient = 1

if (p.value.global.c > alpha){
    rts = 3 # CRS
} else {
        if (p.value.global.n > alpha){
            rts = 2
        } else {
            rts = 1
             }
     }
rts

alpha.l = 1 - (1 - alpha)^(1 / n)

sc.eff = matrix(nrow = 1, ncol = n)

for(qq in 1:n){
    sc.eff[1,qq] = ifelse(p.value.c[1,qq] > alpha.l, 1, 0)
}

# number of scale efficient firms
n.sc.eff = mean(sc.eff)

# number of scale inefficient firms
n.sc.ineff = n - n * n.sc.eff

# those scale inefficient
sc.ineff = matrix(nrow = 1, ncol = n)
for(qq in 1:n){
    if (sc.eff[1,qq] == 0){
        sc.ineff[1,qq] = ifelse(p.value.n[1,qq] > alpha.l, 1, 0)
    } # 1 stands for DRS reason of scale inefficiency
}

# look at correctness
# cbind(t(p.value.c),t(sc.eff),t(p.value.n),t(sc.ineff))

# number of DMU's with DRS as a reason of scale inefficiency
if(n.sc.eff == 1) { 
    n.sc.ineff.drs = 0; n.sc.ineff.irs = 0;
} else {
    if(sum(sc.ineff,na.rm = TRUE) == 0){ 
        n.sc.ineff.drs = 0 
    } else {
        n.sc.ineff.drs = sum(sc.ineff,na.rm = TRUE)/n.sc.ineff
        }
      n.sc.ineff.irs = 1 - n.sc.ineff.drs
    }


cat("               number of observations is: ",n,"\n")
cat("         p-value that globally is CRS is: ",p.value.global.c,"\n")
if (p.value.global.c < alpha){
cat("        p-value that globally is NIRS is: ",p.value.global.n,"\n")
}
cat("           RTS is {1:VRS, 2:NIRS, 3:CRS}: ",rts,"\n")
cat("                           used alpha is: ",alpha,"\n")
cat("                  scale efficient, share: ",n.sc.eff,"\n")
cat("     scale inefficient due to DRS, share: ",n.sc.ineff.drs,"\n")
cat("     scale inefficient due to IRS, share: ",n.sc.ineff.irs,"\n")
cat("   cannot be computed in test 1, average: ",mean(q2),"\n")
cat("   cannot be computed in test 2, average: ",mean(q4),"\n")
tmp = list(sc.eff = c(sc.eff), sc.ineff = c(sc.ineff))
return(tmp)
}


#####################

# HOW TO USE

#####################

#data(ccr)
#y = t(matrix(c(ccr$y1,ccr$y2,ccr$y3),nrow = 70,ncol = 3))
#x = t(matrix(c(ccr$x1,ccr$x2,ccr$x4,ccr$x5),nrow = 70,ncol = 4))

#rts.test(x = x, y = y, orient = 1, B = 100)
#rts.test(x = x, y = y, orient = 2, B = 100)