2015-09-26 42 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輸入包含在計算中。所以我開始相信這個問題可能不是代碼,但它對於計算來說根本不包括這個變量是最合適的。你認爲這可能是真的嗎?我再次更新了問題以包含此信息。謝謝! –