I was researching methods to rank forecasts and found this paper, A novel ranking procedure for forecasting approaches using Data Envelopment Analysis, I've been working through the article and setting up my data, but I can't seem to replicate their LP formulation in R?
This is the formulation I'm referring to: LP Formulation
Here is an example I'm working through to try and replicate their results with the above formulation as a reference. The data is based off of 'Table 2. Log values of data of illustrative example' from the above mentioned article.
library(lpSolve)
library(nonparaeff)
DMU = c("FOR01", "FOR02", "FOR03", "FOR04", "FOR05")
log.data = matrix(data = as.numeric(c("1.794","1.575","3.576"
,"2.228","2.106","6.628"
,"2.399","1.871","6.354"
,"2.619","1.983","5.849"
,"2.559","1.541","5.676")), ncol = 3, byrow = TRUE)
colnames(log.data) = c("M1", "M2", "M3")
rownames(log.data) = DMU
THETA = c(-1,-1,-1,0)
add.to.one = c(1,1,1,1,1) # Constraint so each lambda adds up to one.
f.obj = c(1)
f.con = cbind(THETA, rbind(t(log.data), add.to.one))
f.dir = c("<=","<=","<=","=")
f.rhs = c(1.794,1.575,3.576,1)
lp2(direction = "min", f.obj, f.con, f.dir, f.rhs, free.var = c(1))
I'm using the packages lpsolve and nonparaeff; nonparaeff extends the lp() function so it can handle free variables.
Using this code I end up getting the error:
"Error: no feasible solution found".
However, in the article they end up with theta as 0 or 1 and the first lambda returning as 1. So I must be doing something wrong.
Did I apply the fourth constraint correctly ('add.to.one')? Also, lp() already assumes each variable to be >= to zero, but does it assume anything else?
Did I correctly translate the formulation to R? Am I using the lp2 function correctly?
I've looked through other similar lp questions on here, but I don't see very many on free variables. But please link me to other questions if you think otherwise.
Thanks in advance.
lpSolve is somewhat of a pain to use. However, it should not be too difficult to implement a free variable using a technique called variable splitting. I.e. replace the free variable x by xplus-xmin with xplus,xmin>=0. Not both can be non-zero (then the basis matrix would be singular).
I think this is the correct updated code:
library(lpSolve)
DMU = c("FOR01", "FOR02", "FOR03", "FOR04", "FOR05")
log.data = matrix(data = as.numeric(c("1.794","1.575","3.576"
,"2.228","2.106","6.628"
,"2.399","1.871","6.354"
,"2.619","1.983","5.849"
,"2.559","1.541","5.676")), ncol = 3, byrow = TRUE)
colnames(log.data) = c("M1", "M2", "M3")
rownames(log.data) = DMU
theta = c(-1,-1,-1,0)
add.to.one = c(1,1,1,1,1) # Constraint so each lambda adds up to one.
f.obj = c(0,0,0,0,0,1,-1)
f.con = cbind(rbind(t(log.data), add.to.one),theta,-theta)
f.dir = c("<=","<=","<=","=")
f.rhs = c(1.794,1.575,3.576,1)
r <- lp(direction = "min", f.obj, f.con, f.dir, f.rhs)
Better would be to use a more capable tool such as CVXR.