rm(list = ls())
# simple data
dat = data.frame(
y = c(1,4,3,3,5,6,9,8,11,12,10),
x = c(1,2,3,5,6,4,5,7,7,8,9)
)
dat = dat[order( dat$x ), ]

# plot
plot(dat$x,dat$y, xlim = c(0,10), ylim = c(0, 18), pch = 19, lwd = 8)

x = matrix(dat$x, nrow = 1)
y = matrix(dat$y, nrow = 1)

# FDH
library(FEAR)
eff_fdh = fdh(XOBS = x, YOBS = y, ORIENTATION = 2)
dat$eff_fdh = c( eff_fdh[1,] )
dat$y_fdh = dat$y / dat$eff_fdh
# FHD frotnier
plot(dat$x,dat$y, xlim = c(0,10), ylim = c(0, 18), pch = 1, lwd = 8)
points(dat$x,dat$y_fdh, pch = 19, lwd = 2, col = "red")

# DEA
dat$eff_dea_vrs = dea(XOBS = x, YOBS = y, RTS = 1, ORIENTATION = 2)
dat$eff_dea_nrs = dea(XOBS = x, YOBS = y, RTS = 2, ORIENTATION = 2)
dat$eff_dea_crs = dea(XOBS = x, YOBS = y, RTS = 3, ORIENTATION = 2)

# DEA frotnier(s)
dat$y_dea_vrs = dat$y / dat$eff_dea_vrs
dat$y_dea_nrs = dat$y / dat$eff_dea_nrs
dat$y_dea_crs = dat$y / dat$eff_dea_crs

plot(dat$x,dat$y, xlim = c(0,10), ylim = c(0, 18), pch = 19, lwd = 8)
lines(dat$x,dat$y_dea_vrs, col = "blue", lty = 1, lwd = 1)
lines(dat$x,dat$y_dea_nrs, col = "red", lty = 2, lwd = 6)
lines(dat$x,dat$y_dea_crs, col = "green", lty = 1, lwd = 2)

# aggregate
colMeans(dat)

# Weights
weights.price.independent=function(y){
n=ncol(y)
M=nrow(y)
ratio=matrix(,M,n)
for (i in 1:n){
 for (j in 1:M){
  ratio[j,i]=y[j,i]/(sum(y[j,]))
}
}
w=c(1:n)
for (i in 1:n){
 w[i]=sum(ratio[,i])/M
}
w
}
dat$w = weights.price.independent(y)
eff_dea_crs_ave = weighted.mean(dat$eff_dea_crs,dat$w)
eff_dea_nrs_ave = weighted.mean(dat$eff_dea_nrs,dat$w)
eff_dea_vrs_ave = weighted.mean(dat$eff_dea_vrs,dat$w)

aggr_dea = cbind(
unweighted = c(mean(dat$eff_dea_crs), mean(dat$eff_dea_nrs), mean(dat$eff_dea_vrs)),
weighted = c(eff_dea_crs_ave, eff_dea_nrs_ave, eff_dea_vrs_ave)
)
aggr_dea

# bootstrap VRS

t1 = boot.sw98(XOBS = x, YOBS = y, RTS = 1, ORIENTATION = 2,
	NREP = 2000, alpha = 0.05)

t1$bias
t1$var
t1$conf.int
t1$dhat
t1$dhat.bc
t1$bias.flag
t1$bias / sqrt(t1$var)

# plot VRS case

# DEA frotnier(s)
dat$y_dea_vrs = dat$y / dat$eff_dea_vrs
dat$y_dea_vrs1 = dat$y / t1$conf.int[,1]
dat$y_dea_vrs2 = dat$y / t1$conf.int[,2]

plot(dat$x,dat$y, xlim = c(0,10), ylim = c(0, 18), pch = 19, lwd = 8)
lines(dat$x,dat$y_dea_vrs, col = "blue", lty = 1, lwd = 3)
lines(dat$x,dat$y_dea_vrs1, col = "red", lty = 2, lwd = 4)
lines(dat$x,dat$y_dea_vrs2, col = "red", lty = 2, lwd = 4)

# table with results
cbind(
x = dat$x,
bias = t1$bias,
CI.lower = t1$conf.int[,1],
Bias.corrected = t1$dhat.bc,
CI.upper = t1$conf.int[,2],
if.to.do = t1$bias / sqrt(t1$var)
)

# scale efficiency

plot(dat$x,dat$y, xlim = c(0,10), ylim = c(0, 18), pch = 19, lwd = 8)
lines(dat$x,dat$y_dea_vrs, col = "blue", lty = 1, lwd = 1)
lines(dat$x,dat$y_dea_nrs, col = "red", lty = 2, lwd = 6)
lines(dat$x,dat$y_dea_crs, col = "green", lty = 1, lwd = 2)

cbind(
x = dat$x,
se = dat$y_dea_vrs / dat$y_dea_crs
)

# Sourse of SE
cbind(
x = dat$x,
se = dat$y_dea_vrs / dat$y_dea_crs,
se.ineff = dat$y_dea_vrs / dat$y_dea_nrs
)

# load "RTS_test_Simar_Wilson_02.R" file
rts.test(x = x, y = y, orient = 2, B = 2000, alpha = 0.05)