2015-09-26 37 views
0

我正在应用CCR数据包络分析模型来基准库存数据。为此,我运行了DEA发布的here论文中的R代码。本文使用lpSolve解决线性问题。 lpSolve的文档是here如何更改lpSolve中的重量限制?线性编程,最大化CCR DEA

该文献随附步骤通过关于如何执行下面的模型在R.

步说明

背景信息:

数据包络分析(= DEA)是一个面向数据的方法来评估执行称为决策制定单元的一组实体,其将多个数据输入转换为多个数据输出.DEA用于通过考虑其与生产单位的接近度来确定生产单位的相对效率来评估和比较决策单元(DMU)的性能效率前沿。

我们假设有n个DMU需要评估。每个DMU消耗不同数量的不同输入来产生不同的输出。具体消耗输入量并产生输出量。我们进一步假设≥0,并且每个DMU至少有一个正输入和一个正输出值。在CCR模型中,目标是使用线性编程来确定权重,从而最大化每个决策单元的综合效率。这是通过比率定义为:

enter image description here

,然后将其转变成线性问题,像这样:

enter image description here

应用:

我申请这个效率模型以股票数据CAPM风险因素:

它是一个5个输入 - 输出1的情况下

enter image description here

这里的目标是找到“最好”的权重集(1,1 ... 5),用于被分析的每个资产。这里使用的术语“最佳”是指当将这些权重分配给所有DMU /股票(1 ... 10)的输入和输出时,上述比率相对于所有其他比例最大化。

这里是一个已经被标准化,使得所有的值位于0和1

> dput(testdfst) 
structure(list(Name = structure(1:10, .Label = c("Stock1", "Stock2", 
"Stock3", "Stock4", "Stock5", "Stock6", "Stock7", "Stock8", "Stock9", 
"Stock10"), class = "factor"), Date = structure(c(14917, 14917, 
14917, 14917, 14917, 14917, 14917, 14917, 14917, 14917), class = "Date"), 
    `(Intercept)` = c(0.454991569278089, 1, 0, 0.459437188169979, 
    0.520523252955415, 0.827294243132907, 0.642696631099892, 
    0.166219881886161, 0.086341470900152, 0.882092217743293), 
    rmrf = c(0.373075150411683, 0.0349067218712968, 0.895550280607866, 
    1, 0.180151549474574, 0.28669170468735, 0.0939821798173586, 
    0, 0.269645291515763, 0.0900619760898984), smb = c(0.764987877309785, 
    0.509094491489323, 0.933653313048327, 0.355340700554647, 
    0.654000372286503, 1, 0, 0.221454091364611, 0.660571586102851, 
    0.545086931342479), hml = c(0.100608151187926, 0.155064872867367, 
    1, 0.464298576152336, 0.110803875258027, 0.0720803195598597, 
    0, 0.132407005239869, 0.059742053684015, 0.0661623383303703 
    ), rmw = c(0.544512524466665, 0.0761995312858816, 1, 0, 0.507699534880555, 
    0.590607506295898, 0.460148690870041, 0.451871218073951, 
    0.801698199214685, 0.429094840372901), cma = c(0.671162426988512, 
    0.658898571758625, 0, 0.695830176886926, 0.567814542084284, 
    0.942862571603074, 1, 0.37571611336359, 0.72565234813082, 
    0.636762557753099), Returns = c(0.601347600017365, 0.806071701848376, 
    0.187500487065719, 0.602971876359073, 0.470386289298666, 
    0.655773224143057, 0.414258177255333, 0, 0.266112191477882, 
    1)), .Names = c("Name", "Date", "(Intercept)", "rmrf", "smb", 
"hml", "rmw", "cma", "Returns"), row.names = c("Stock1.2010-11-04", 
"Stock2.2010-11-04", "Stock3.2010-11-04", "Stock4.2010-11-04", 
"Stock5.2010-11-04", "Stock6.2010-11-04", "Stock7.2010-11-04", 
"Stock8.2010-11-04", "Stock9.2010-11-04", "Stock10.2010-11-04" 
), class = "data.frame") 

我现在通过从在所提到的文件采取将R代码运行上述数据之间的原始样值Z得分数据开始。请注意,从计算中排除testdfst数据帧的截取变量。代码

require(lpSolve) 

dea_results <- list() 

namesDMU <- testdfst[1] 
inputs <- testdfst[c(4,5,6,7,8)] 
outputs <- testdfst[9] 

N <- dim(testdfst)[1] # number of DMU 
s <- dim(inputs)[2] # number of inputs 
m <- dim(outputs)[2] # number of outputs 

f.rhs <- c(rep(0,N),1) # RHS constraints 
f.dir <- c(rep("<=",N),"=") # directions of the constraints 
aux <- cbind(-1*inputs,outputs) # matrix of constraint coefficients in (6) 
for (i in 1:N) { 
    f.obj <- c(0*rep(1,s),outputs[i,]) # objective function coefficients 
    f.con <- rbind(aux, c(as.matrix(inputs[i,]), rep(0, m))) # add LHS of bTz=1 
    results <-lp("max",f.obj,f.con,f.dir,f.rhs,scale=1,compute.sens=TRUE) # solve LPP 
    multipliers <- results$solution # input and output weights 
    efficiency <- results$objval # efficiency score 
    duals <- results$duals # shadow prices 
    if (i==1) { 
    weights <- multipliers 
    effcrs <- efficiency 
    lambdas <- duals [seq(1,N)] 
    } else { 
    weights <- rbind(weights,multipliers) 
    effcrs <- rbind(effcrs , efficiency) 
    lambdas <- rbind(lambdas,duals[seq(1,N)]) 
    } 
} 

report <- as.data.frame(cbind(effcrs,weights)) 
colnames(report) <- c('efficiency',names(inputs), 
         names(outputs)) # header 
rownames(report) <- namesDMU[,1] 

说明:

enter image description here

enter image description here

可能与我的问题: 音符的Z> = 0的一部分?纠正我,如果我错了,但lpSolve必须采取这种默认,因为这个限制不是直接包含在代码中?如果我能正确识别这部分我可以调整它,以便它解决我的问题。

该问题:

这是最终的结果:

> report 

     efficiency rmrf  smb  hml  rmw  cma Returns 
Stock1 0.5674100 0 0.000000 0.1769187 0.0000000 1.4634319 0.9435640 
Stock2 1.0000000 0 0.000000 0.0000000 0.7713486 1.4284803 1.2405844 
Stock3 1.0000000 0 1.071061 0.0000000 0.0000000 7.4588210 5.3333195 
Stock4 1.0000000 0 1.930968 0.0000000 0.7427269 0.4510419 1.6584521 
Stock5 0.5218197 0 0.000000 0.2080023 0.0000000 1.7205486 1.1093429 
Stock6 0.5498426 0 0.000000 9.3299443 0.0000000 0.3473408 0.8384645 
Stock7 1.0000000 0 3.260381 0.0000000 0.0000000 1.0000000 2.4139536 
Stock8 0.0000000 0 0.000000 0.0000000 2.2130199 0.0000000 0.0000000 
Stock9 0.2756548 0 0.000000 11.5264372 0.0000000 0.4291132 1.0358592 
Stock10 1.0000000 0 0.000000 0.1875005 0.0000000 1.5509620 1.0000000 

的RMRF输入已被分配在所有情况下为0的权重。为了使这种效率度量在我应用它的情况下是相关的,有必要在每次计算中考虑所有5个输入。

我想解决的问题是进一步限制权重。我要让所有5个输入被纳入这个问题,所以我想限制“V”重像这样:

0.2 <= V(n)/V(n+1) <=5 

这是一样的

1/5 <= V(n)/V(n+1) <= 5/1 

这样,所有的权重应该是限制为最大5倍,比其他任何重量还要小5倍。

谢谢任何​​能提供一些意见的人。随着我取得更多进展,我将继续更新这篇文章。

UPDATE

我一直在检查从模型中每个变量,也尝试了一些不同的数据集。我还在数据框中交换了rmrf和smb的位置,试图查看它是错误的代码还是数据。

 efficiency  smb  rmrf  hml  rmw  cma Returns 
Stock1 1.0000000 0.2540156 0.0000000 1.1812905 1.043781 0.03466956 1.000000 
Stock2 1.0000000 1.2150960 0.0000000 0.5443910 2.854127 0.00000000 2.422061 
Stock3 1.0000000 0.0000000 0.0000000 1.0000000 0.000000 1.48815164 1.104670 
Stock4 1.0000000 1.2348830 0.0000000 0.0000000 4.088438 0.89216307 3.674087 
Stock5 0.8151261 0.0000000 0.0000000 0.5921236 1.173539 0.61261329 1.231220 
Stock6 0.1491688 0.0000000 0.0000000 2.8715984 1.262761 0.00000000 1.148347 
Stock7 1.0000000 2.0599619 0.0000000 0.0000000 0.000000 1.03785902 1.520484 
Stock8 1.0000000 0.8346388 0.1206426 0.0000000 1.862840 0.00000000 1.535768 
Stock9 0.0000000 0.0000000 0.0000000 0.0000000 1.167498 0.00000000 0.000000 
Stock10 0.8065724 0.0000000 0.6973118 2.9529776 1.318778 0.00000000 1.357774 

在尝试了不同的数据集之后,似乎有包括RMRF在内的例外(如上所述)。我开始相信我的问题的答案可能非常简单。程序发现它最适合简单地不包括我的RMRF输入到计算中。解决这个问题的办法就是对权重进行额外限制。

这里是具有非标准化的输入中的溶液(在数据完全不同的分布号码)

  efficiency smb  rmrf hml rmw cma Returns 
Stock1   1 0 1.0775390 0 0 0 346.1943 
Stock2   0 0 1.3576882 0 0 0 0.0000 
Stock3   0 0 0.7443477 0 0 0 0.0000 
Stock4   0 0 0.6879011 0 0 0 0.0000 
Stock5   0 0 1.2115507 0 0 0 0.0000 
Stock6   0 0 1.1305046 0 0 0 0.0000 
Stock7   0 0 1.2472639 0 0 0 0.0000 
Stock8   0 0 1.4552312 0 0 0 0.0000 
Stock9   0 0 1.1937843 0 0 0 0.0000 
Stock10   0 0 1.3569824 0 0 0 0.0000 

将继续更新的实例。

回答

1

寻址第二个问题的第一,通过V(n+1)约束0.2 <= V(n)/V(n+1)相乘可以转化为:

0.2 * V(n+1) <= V(n) 

由于这是一个线性约束,它可直接加入到模型中。同样,不平等的另一侧可以模拟成线性约束

V(n) <= 5 * V(n+1) 

关于你的第一个问题,因为你还没有真正描述的其他算法不是链接到文件链接到另一个文件描述变量,你至少没有给我足够的上下文来回答。如果您期望rmrf系数取正值,则可以添加一个约束,要求rmrf系数的总和超过您选择的某个正值;这将迫使系数中的至少一个是正数。

+0

嗨josilber,我编辑了这个问题,并添加了更多关于问题实际上试图解决的信息。当我第一次写这个问题时,我没有意识到我错过了主要想法。感谢您提出。至于额外的限制,我想补充说,将它分成两个绝对有帮助。 –

+0

嗨josil,通过试验和错误,我发现了一些情况下,我的RMRF输入包含在计算中。所以我开始相信这个问题可能不是代码,但它对于计算来说根本不包括这个变量是最合适的。你认为这可能是真的吗?我再次更新了问题以包含此信息。谢谢! –