数据包络分析DEA

测量一些决策部门的生产效率的方法。

简单说，现在有三个人，你可以叫他们DMU，a用一个甲生产了一个乙，b用两个甲生产了一个乙，c用三个甲生产了一个乙，显而易见a的效率最高，我们就把他定义为技术前沿面，他的效率为1。那么b就是1/2,c就是1/3。这样我们就可以计算出每个人的效率水平了。

1. 基本概念

数据包络分析（Data envelopment analysis，DEA）是运筹学和研究经济生产边界的一种方法。该方法一般被用来测量一些决策部门的生产效率。

关于DEA的概念和使用入门，可以参见

资料包络分析法 - 百度文库​wenku.baidu.com/link?url=5QTVRpdlhN82xkuKnmBox0csPCggStzVo7SEqakP-ipLwu5757VouBDyQGUl0skclxgGctpirFukXA4pi_bWbwh5KTVmMwFbKw1GSk2SspK​编辑

，其中基本概念和应用举例讲解的比较通俗。

为了更详细的了解DEA，以下几个基本概念有必要了解一下。

1. 决策单元

一个经济系统通常可以看成是一个“公司”,通过投入一定数量的生产要素并产出一定数量的“产品”，尽管这种活动的具体内容各不相同,但其目的都是尽可能地使这一活动取得最大的“效益”。由于从“投入”到“产出”需要经过一系列决策才能实现,或者说,由于“产出”是决策的结果,所以这样的“公司”被称为决策单元（decision making unit，DMU）。所以,可以认为,每个决策单元都代表一定的经济意义,它的基本特点是具有一定的输入和输出,并且将输入转化成输出的过程中,努力实现自身的决策目标。

3. 有效生产（前沿）

对于生产可能集， 如果不存在(�,�)∈�,如果不存在�′≥�,(�,�′)∈� ,则称（X,Y）为有效生产活动，此投入产出对应一个前沿，由众多“有效生产”构成的凸包即为前沿。

2. DEA方法的演化。

Farrell （1957）[1]首先提出了效率概念，将生产效率分为技术效率和配置效率来衡量生产前沿。此后，Charnes, Cooper和Rhodes（1978）[2]开发了CCR模型。Banker, Charnes和Cooper（1984）[3]进一步将CCR模型修正为BCC模型。CCR和BCC模型都测量径向效率，允许输入或输出项按相等比例增加或减少。然而，这一假设并非适用于所有情况。Tone[4]在2001年提出了SBM模型，一种以差分变量为测量基础，考虑输入和输出之间的差异(slack)，并利用非径向估计方法和标量来表示0到1之间的SBM效率值。Malmquist(1953)[5]首先提出了Malmquist指数来解释动态效率。Fare(1994)[6] 将该概念扩展到测量跨时间效率变化。然而，这些模型没有考虑跨期持续活动的影响，因此不太适合测量长期效率。Fare和Grosskopf(1996)[7]首先建立了动态DEA概念，设计了一种动态分析形式，然后提出了动态模型的延滞（carryover）变量。Tone and Tsutsui(2010)[8] 随后扩展到基于加权松弛的动态DEA方法，包括四种类型的连接活动：（1）期望（好的）；（2）非期望（坏的）；（3）可任意支配的（自由的）；（4）不可任意支配的（固定的）。Ruttan, Binswanger,Hayami, Wade, and Weber (1978) [9]将metafrontier (MF)定义为一个包络线，其中包含所有组的生产前沿，可以在一个公共基准下对不同组的效率进行测量。Battese and Rao (2002) [10]and Battese, Rao, and O’donnell (2004) [11]然后证明了利用MF模型可以比较不同组的技术效率。Portela和Thanassoulis(2008)[12]提出了一个凸MF概念，可以考虑所有组的技术，最先进的技术生产水平，以及组之间的沟通，并可以进一步扩展，以提高业务绩效。O’Donnell, Rao, and Battese (2008) [13]提出了一种利用输出距离函数定义技术效率的MF模型，该模型能够准确地计算组和MF技术效率，发现所有组的技术水平都优于任何一个组的技术水平。

3. CCR模型

3.1 模型的基本推导

CCR模型是最早也最经典的DEA模型，模型如下：

经过Charnes-Cooper变换：

得到以下模型：

对于等式的处理，通常为：

则，可以进一步将上述模型写为矩阵形式：

则其对偶形式可以表示为如下形式：

简写为如下形式：

各类论文中最常见的CCR模型形式

3.2 代码

我们使用较成熟的MATLAB中linprog()函数作为线性规划的求解器开展代码编写.

为了更好的理解代码，我把linprog()函数官方说明中比较重要的几个基本点先列出来。

linprog Linear programming.

for input:

X = linprog(f,A,b) attempts to solve the linear programming problem:

min f'*x subject to: A*x <= b

X = linprog(f,A,b,Aeq,beq,LB,UB) defines a set of lower and upper bounds on the design variables, X, so that the solution is in the range LB <= X <= UB. Use empty matrices for LB and UB. if no bounds exist. Set LB(i) = -Inf if X(i) is unbounded below; set UB(i) = Inf if X(i) is unbounded above.

for output:

X = linprog(PROBLEM) finds the minimum for PROBLEM. PROBLEM is a structure with the vector 'f' in PROBLEM.f, the linear inequality constraints in PROBLEM.Aineq and PROBLEM.bineq, the linear equality constraints in PROBLEM.Aeq and PROBLEM.beq, the lower bounds in PROBLEM.lb, the upper bounds in PROBLEM.ub, the start point in PROBLEM.x0, the options structure in PROBLEM.options, and solver name 'linprog' in PROBLEM.solver. Use this syntax to solve at the command line a problem exported from OPTIMTOOL.

[X,FVAL] = linprog(f,A,b) returns the value of the objective function at X:

FVAL = f'*X.

[X,FVAL,EXITFLAG] = linprog(f,A,b) returns an EXITFLAG that describes the exit condition. Possible values of EXITFLAG and the corresponding exit conditions are

1 linprog converged to a solution X.

0 Maximum number of iterations reached.

-2 No feasible point found.

-3 Problem is unbounded.

-4 NaN value encountered during execution of algorithm.

-5 Both primal and dual problems are infeasible.

-7 Magnitude of search direction became too small; no further

progress can be made. The problem is ill-posed or badly

conditioned.

[X,FVAL,EXITFLAG,OUTPUT] = linprog(f,A,b) returns a structure OUTPUT with the number of iterations taken in OUTPUT.iterations, maximum of constraint violations in OUTPUT.constrviolation, the type of algorithm used in OUTPUT.algorithm, the number of conjugate gradient iterations in OUTPUT.cgiterations (= 0, included for backward compatibility), and the exit message in OUTPUT.message.

[X,FVAL,EXITFLAG,OUTPUT,LAMBDA] = linprog(f,A,b) returns the set of Lagrangian multipliers LAMBDA, at the solution: LAMBDA.ineqlin for the linear inequalities A, LAMBDA.eqlin for the linear equalities Aeq, LAMBDA.lower for LB, and LAMBDA.upper for UB.

function [eff,slackX,slackY,Xeff,Yeff ] = dea_ccr( X, Y)
%DEA Data envelopment analysis radial model
%   X表示投入变量矩阵
%   Y表示产出变量矩阵
    
    %1. 提取变量的维度
    [n, m] = size(X);
    s = size(Y,2);
    
    %2. 依据提取的维度信息，创建求解变量的存储矩阵
   
    slackX = nan(n, m);%投入的松弛变量
    slackY = nan(n, s);%产出的松弛变量
    eff = nan(n, 1);%效率值，第一阶段可计算得的
    Eflag = nan(n, 2);%求解的过程信息
    lambda = nan(n, n);% 每个单元的intensity variable
    Xeff = nan(n, m);%投入对应的前沿值
    Yeff = nan(n, s);%产出对应的前沿值
    
    %3. 参数准备及求解线性规划过程
   
        % 对每个 DMU
        for j=1:n
            
            %第一阶段，求解效率θ
                % 目标函数
                f = [zeros(1,n), 1];
                
                % 约束条件
                A = [ X', -X(j,:)';
                     -Y', zeros(s,1)];
                b = [zeros(m,1);
                    -Y(j,:)'];
                Aeq = [];
                beq = [];
                lb = zeros(1, n + 1);
                ub = [];                
                % linprog（）求解
                [z, ~, exitflag] = linprog(f, A, b, Aeq, beq, lb, ub);
                if exitflag ~= 1
                    fprintf(1,'DMU %6d 第一阶段的exitflag: %6d', j, exitflag)
                end
                if isempty(z)
                    fprintf(1,'DMU %6d 第一阶段的解不存在.', j)
                    z = nan(n + 1, 1);
                end
                
                % 提取计算结果
                theta = z(end);
                Eflag(j, 1) = exitflag;
                eff(j) = theta;
    
                
           %第二阶段，计算松弛变量及每个单元各变量对应的前沿值
                
                if ~isnan(theta)
                
                    % 目标函数
                    f = [zeros(1, n), -ones(1, m + s)];

                    % 约束条件
                    Aeq = [ X', eye(m,m)  , zeros(m,s);
                            Y', zeros(s,m), -eye(s,s)];
                    beq = [theta .* Xeval(j,:)';
                            Yeval(j,:)'];
                    lb = zeros(n + s + m, 1);

                    % 求解
                    [z, ~, exitflag] = linprog(f, [], [], Aeq, beq, lb, []);
                    if exitflag ~= 1
                        fprintf(1, 'DMU %6d 第二阶段的exitflag: %6d: %i', j, exitflag)
                    end
                    if isempty(z)
                        fprintf(1,'DMU %6d 第二阶段的解不存在.', j)
                        z = nan(n + m + s, 1);                   
                    end

                    % 提取结果
                    lambda(j,:) = z(1:n);
                    slackX(j,:) = z(n + 1 : n + m);
                    slackY(j,:) = z(n + m + 1 : n + m + s);                
                    Eflag(j, 2) = exitflag;

                    % 计算每个变量的前沿值
                    Xeff(j,:) = repmat(eff(j), 1, m) .* X(j,:) - slackX(j,:);
                    Yeff(j,:) = Y(j,:) + slackY(j,:);

                end
                
        end
end