Deamatlabtoolbox
Uploaded by
Bulent MiranDeamatlabtoolbox
Uploaded by
Bulent MiranI.S.S.
N: 1885-6888
ECONOMIC ANALYSIS
WORKING PAPER SERIES
A Data Envelopment Analysis Toolbox for MATLAB
Inmaculada C. Álvarez, Javier Barbero and José L. Zofío
Working Paper 3/2016
DEPARTAMENTO DE ANÁLISIS ECONÓMICO:
TEORÍA ECONÓMICA E HISTORIA ECONÓMICA
A Data Envelopment Analysis Toolbox for MATLAB
Inmaculada C. Álvarez Javier Barbero José L. Zofı́o
Universidad Autónoma de Madrid
Abstract
Data Envelopment Analysis Toolbox is a new package for MATLAB that includes
functions to calculate the main DEA models. The package includes code for the stan-
dard additive and radial input and output measures, allowing for constant and variable
returns to scale, as well as recent developments related to the directional distance func-
tion, and including both desirable and undesirable outputs when measuring efficiency
and productivity; i.e., Malmquist and Malmquist-Luenberger indices. Bootstrapping to
perform statistical analysis is also included. This paper describes the methodology and
implementation of the functions and reports numerical results with well-known examples
to illustrate their use.
Keywords: Data Envelopment Analysis, distance functions, technical efficiency, MATLAB
JEL codes: D20, 24, C44 .
Affiliation:
Inmaculada C.Álvarez, Javier Barbero, José L. Zofı́o
Department of Economics
Universidad Autónoma de Madrid
28049 Madrid, Spain
E-mail: [email protected], [email protected], [email protected]
URL: http://www.deatoolbox.com/
2 A Data Envelopment Analysis Toolbox for MATLAB
1. Introduction
Data Envelopment Analysis, DEA, has grown in importance over the past decades due to in-
crease in the availability of data related to the performance of Decision Making Units, DMUs,
regardless their market, governmental, or non-for-profit orientation. Basic DEA methods are
included in some standard software packages used by econometricians (e.g., Stata, StataCorp
(2015) with the user-written command by Ji and Lee (2010); LIMDEP, Econometric Software
(2009)); available in dedicated non-commercial software accompanying academic handbooks–
Cooper, Seiford, and Tone (2007), Wilson (2008), Bogetoft and Otto (2011) (these latter
two implemented in R); commercial software–including trials versions, Emrouznejad and Ca-
banda (2014); free-ware programs–Sheel (2000); and even tutorials for spreadsheets, Sherman
and Zhu (2006). Earlier versions of these programs have been reviewed, among others, by
Hollingsworth (2004) and Barr (2004), and there are continuous proposals expanding these
options. While these software packages implement the main DEA models, there is a lack
of a full set of functions for MATLAB (The MathWorks, Inc. 2015), including some recent
theoretical contributions that are missed in the existing software.
The Data Envelopment Analysis Toolbox introduces such set of functions, covering a wide
range of efficiency and productivity models, and reporting numerical results based on classical
examples presented in the literature. Data Envelopment Analysis Toolbox is available as
free software, under the GNU General Public License version 3, and can be downloaded
from http://www.deatoolbox.com, with all the supplementary material (data, examples and
source code) to replicate all the results presented in this paper. The toolbox is also hosted on
an open source repository on GitHub.1
The paper is organized as follows. The following section presents the data structures charac-
terizing the production possibility sets, the structure of the functions, results, etc... Section 3
covers the standard DEA models introduced by Charnes, Cooper, and Rhodes (1978) and
Banker, Charnes, and Cooper (1984) corresponding to the radial input and output efficiency
measures, allowing for constant and variable returns to scale, as well as newer proposals based
on the flexible directional distance function. The non-oriented additive model is also presented
as well as the super-efficiency model for all the previous efficiency measures. Malmquist pro-
ductivity indices and their decomposition into efficiency change and technical change are
shown in section 4, while section 5 deals with the measurement of economic efficiency, and
its decomposition into technical and allocative terms. Section 6 is devoted to the measure-
ment of efficiency with undesirable outputs, most notably environmental efficiency, followed
by Section 7 presenting the Malmquist-Luenberger index. Statistical analyses and hypothe-
ses testing using bootstrapping techniques are presented in Section 8. Advanced options,
including displaying and exporting results can be found in Section 9. Section 10 concludes.
2. Data structures
Data Envelopment Analysis measures productive and economic performance of a set of j =
1, 2, . . . , n observed DMUs (firms, activities, countries, individuals, etc.). These observations
transform a vector of i = 1, 2, . . . , m inputs x ∈ Rm++ into a vector of i = 1, 2, . . . , s outputs
y ∈ Rs++ using the technology represented by the following constant returns to scale produc-
tion possibility set: Pcrs = {(x, y) |x > Xλ, y 6 Y λ, λ > 0}, where X = (x)j ∈ Rs×n , Y =
1
The address of the repository is https://github.com/javierbarbero/DEAMATLAB
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 3
(y)j ∈ Rm×n and λ = (λ1 , . . . , λn )T is a semipositive vector.
Data are managed as regular MATLAB vectors and matrices, constituting the inputs of the
estimation functions. All estimation functions return a structure deaout that contains fields
with the estimation results as well as the input of the estimation function. Fields can be
accessed directly using the dot notation and the whole structure can be used as an input to
other functions that print or export results (e.g., deadisp).
Some of the fields of the deaout structure are the following:2
X, Y and Yu: contain the inputs, outputs and undesirable outputs variables, respectively.
n and neval: number of DMUs, and number of evaluated DMUs.
m, s and r: number of inputs, outputs and undesirable outputs.
model, orient, rts: strings containing the model type, the orientation, and the returns
to scale assumption.
eff: computed efficiency measure.
slackX, slackY, slackYu: computed input, output and undesirable output slacks.
names: names of the DMUs.
3. Basic DEA models
3.1. Radial input oriented model: Constant and variable returns to scale
Based on the data matrix (X, Y ), we measure the input oriented efficiency of each observation
o by solving n times the following linear programming problem–known as the ccr model:3
min θ (1)
θ,λ
subject to
θxo ≥ Xλ
Y λ ≥ yo
λ ≥ 0.
The optimal solution to this program—characterizing a technology with constant returns to
∗ . The constraints require the observation (θ
scale, is denoted by θcrs crs xo , yo ) to belong to
Pcrs , while the objective seeks the minimum θcrs that reduces the input vector xo radially
∗
to θcrs xo while remaining in Pcrs . A feasible solution signaling radial efficiency is θcrs = 1.
2
For a full list see the help of the function typing help deaout in MATLAB.
3
This program corresponds to the so-called “envelopment form” of the formulation introduced by Charnes
et al. (1978). Correspondences with the dual approaches (“multipliers form”) can be found in Cooper et al.
(2007).
4 A Data Envelopment Analysis Toolbox for MATLAB
∗
Therefore if θcrs < 1, the observation is radially inefficient and (λX, λY ) outperforms (xo , yo ).
With regard to this property, we define the additional input excesses and outputs shortfalls by
the following slack vectors: s− ∈ Rm and s+ ∈ Rs , respectively. Therefore: s− = θcrs ∗ x − Xλ,
o
and s+ = Y λ − yo with s− ≥ 0 and s+ ≥ 0 for any feasible solution (θ, λ).
To obtain the possible input excesses and output shortfalls, the following second stage program
∗
that incorporates the optimal value θcrs and corrects radial inefficiency is solved:
max ω = es− + es+ (2)
λ, s− , s+
subject to
s− = θcrs
∗
xo − Xλ
s+ = Y λ − yo
λ ≥ 0, s− ≥ 0, s+ ≥ 0,
Pm − Ps +
where e = (1, . . . , 1)T so es− = i=1 si and es+ = i=1 si .
∗ , λ∗ , s−∗ , s+∗ )
As a result, an observation is technically efficient if the optimal solution (θcrs
∗ −∗ +∗
of the two above programs satisfy θcrs = 1, s = 0, and s = 0, so no equiproportional
contraction of inputs, and individual inputs reductions and outputs increases are possible
(Pareto-Koopmans Efficiency).4
The measurement of technical efficiency assuming variables returns to scale as introduced by
Banker et al. (1984)–known as the bcc model, considers the following production possibility
set Pvrs = {(x, y) |x > Xλ, y 6 Y λ, eλ = 1, λP > 0.}. Therefore, the only difference with
the crs model is the adjunction of the condition nj=1 λj = 1. Calculating the vrs efficiency,
along with the subsequent second stage program analogous to (2), yields the corresponding
∗ , λ∗ , s−∗ , s+∗ ). As before, an observation is efficient with respect to the
optimal solution (θvrs
vrs technology if the optimal solution of the two programs satisfy θvrs ∗ = 1, s−∗ = 0, and
+∗
s = 0.
Finally, the scale efficiency of each observation is calculated as the ratio of the vrs to the
crs scores: SE = θcrs∗ /θ ∗ . As a result the radial technical efficiency of an observation can
vrs
be decomposed into its variable returns to scale efficiency (“pure” technical efficiency, P T E)
∗
and scale efficiency: T E = θcrs ∗ × SE.
= P T E × SE = θvrs
The radial input oriented model can be computed in MATLAB using the dea(X, Y, ...)
function with the orient parameter set to io (input oriented). The returns to scale assump-
tion can be specified by setting the rts parameter to crs (constant returns to scale; default)
or vrs (variable returns to scale). Results are returned in a deaout structure and can be
accessed directly (see Section 2) or displayed using the deadisp function.5
> load 'deadataFLS'
4
It is possible to solve both programs in a single stage formulation employing a “non-Archimedian” in-
finitesimal constant , e.g., (Fried, Lovell, and Schmidt 1993, 140). However, this may result in computational
inaccuracies and erroneous results.
5
See Section 9.3 for advanced uses of the deadisp function. The example for the input, output, directional
and additive models corresponds to (Fried et al. 1993, 122).
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 5
> io = dea(X, Y, 'orient', 'io');
> deadisp(io);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| Theta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.6223| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 0.8199| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 1.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.3104| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 0.5556| 4.4444| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 1.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.7577| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.8201| 1.6402| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 0.5000| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 1.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
> io_vrs = dea(X, Y, 'orient', 'io', 'rts', 'vrs');
> deadisp(io_vrs);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: io (Input oriented)
Returns to scale: vrs (Variable)
--------------------------------------------------------------------
DMU| X1| X2| Y| Theta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.8700| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 1.0000| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 1.0000| 0.0000| 0.0000| 0.0000|
6 A Data Envelopment Analysis Toolbox for MATLAB
5| 18.0000| 14.0000| 8.0000| 0.7116| 0.0000| 0.0000| 2.6984|
6| 23.0000| 6.0000| 9.0000| 1.0000| 0.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 1.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 1.0000| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 1.0000| 0.0000| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 0.5000| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 1.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
The scale efficiency can be calculated using the deascaleeff(X, Y, ...) function. The
function parameters are the same as those of the dea function, although the rts parameter
specified will be omitted since both are needed in order to compute scale efficiency.
> io_scale = deascale(X, Y, 'orient', 'io');
> deadisp(io_scale);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: io (Input oriented)
Returns to scale: scaleeff (Scale efficiency)
---------------------------------
DMU| CRS| VRS| ScaleEff|
---------------------------------
1| 1.0000| 1.0000| 1.0000|
2| 0.6223| 0.8700| 0.7153|
3| 0.8199| 1.0000| 0.8199|
4| 1.0000| 1.0000| 1.0000|
5| 0.3104| 0.7116| 0.4361|
6| 0.5556| 1.0000| 0.5556|
7| 1.0000| 1.0000| 1.0000|
8| 0.7577| 1.0000| 0.7577|
9| 0.8201| 1.0000| 0.8201|
10| 0.5000| 0.5000| 1.0000|
11| 1.0000| 1.0000| 1.0000|
---------------------------------
3.2. Radial output oriented model: Constant and variable returns to scale
It is possible to measure the output oriented technical efficiency of each observation by solving
the following linear program, counterpart to (1):
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 7
max φ (3)
φ,λ
subject to
Xλ ≤ xo
φyo ≤ Y λ
λ ≥ 0.
In this case, the optimal solution is denoted by φ∗crs with the constraints ensuring that
(xo , φ∗crs yo ) belongs to Pcrs . Now the objective seeks the maximum φcrs that increases
the output vector yo radially to φ∗crs yo while remaining in Pcrs . A feasible solution signaling
radial efficiency is φ∗crs = 1. Therefore if φ∗crs > 1, the observation is radially inefficient and
(λX, λY ) outperforms (xo , yo ). Again, there might be farther input excesses and outputs
shortfalls, with: s− = xo − Xλ, and s+ = Y λ − φ∗crs yo with s− ≥ 0 and s+ ≥ 0 for any fea-
sible solution (φ, λ). To calculate these slacks in a second stage, the corresponding program
incorporating the optimal value φ∗crs is needed:
max ω = es− + es+ (4)
λ, s− , s+
subject to
s− = xo − Xλ
s+ = Y λ − φ∗crs yo
λ ≥ 0, s− ≥ 0, s+ ≥ 0.
Finally, it is also possible to calculate the technical efficiency with respect to Pvrs by solving
the programs for the radial component andP its associated input and output slacks analogous to
(3) and (4), but adding the vrs contraint nj=1 λj = 1. If φ∗vrs = 1 the observation is radially
efficient, while it is technically efficient if s−∗ = 0 and s+∗ = 0 in the second stage. The scale
efficiency defines as SE = φ∗vrs /φ∗crs and radial efficiency can be now decomposed into pure
technical efficiency, P T E, and scale efficiency: T E = φ∗crs = P T E × SE = φ∗vrs × SE.
The radial output oriented model is computed in MATLAB using the same dea(X, Y, ...)
function with the orient parameter set to oo (output oriented). Again, the returns to scale
assumption can be specified by setting the rts parameter to crs (constant returns to scale;
default) or vrs (variable returns to scale).
> oo = dea(X, Y, 'orient', 'oo');
> deadisp(oo);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
8 A Data Envelopment Analysis Toolbox for MATLAB
Model: radial
Orientation: oo (Output oriented)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| Phi| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 1.6070| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 1.2197| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 1.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 3.2220| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 1.8000| 8.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 1.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 1.3198| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 1.2194| 2.0000| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 2.0000| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 1.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
> oo_vrs = dea(X, Y, 'orient', 'oo', 'rts', 'vrs');
> deadisp(oo_vrs);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: oo (Output oriented)
Returns to scale: vrs (Variable)
--------------------------------------------------------------------
DMU| X1| X2| Y| Phi| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 1.5075| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 1.0000| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 1.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 3.2039| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 1.0000| 0.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 1.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 1.0000| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 1.0000| 0.0000| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 1.1698| 5.0000| 11.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 1.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 9
> oo_scale = deascale(X, Y, 'orient', 'oo');
> deadisp(oo_scale);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: oo (Output oriented)
Returns to scale: scaleeff (Scale efficiency)
---------------------------------
DMU| CRS| VRS| ScaleEff|
---------------------------------
1| 1.0000| 1.0000| 1.0000|
2| 1.6070| 1.5075| 1.0660|
3| 1.2197| 1.0000| 1.2197|
4| 1.0000| 1.0000| 1.0000|
5| 3.2220| 3.2039| 1.0056|
6| 1.8000| 1.0000| 1.8000|
7| 1.0000| 1.0000| 1.0000|
8| 1.3198| 1.0000| 1.3198|
9| 1.2194| 1.0000| 1.2194|
10| 2.0000| 1.1698| 1.7097|
11| 1.0000| 1.0000| 1.0000|
---------------------------------
3.3. The directional model: Constant and variable returns to scale
Chambers, Chung, and Färe (1996) introduced a measure of efficiency that projects observa-
tion (xo , yo ) in a pre-assigned direction g = −gx− , gy+ 6= 0m+s , gx− ∈ Rm and gy+ ∈ Rs , in a
proportion β. The associated linear program is:
max β (5)
β,λ
subject to
Xλ ≤ xo − βgx−
Y λ ≥ yo + βgy+
λ ≥ 0.
In this occasion the optimal solution to this program corresponds to βcrs ∗ . Now β ∗
crs =
∗
0 signals directional inefficiency. Therefore if βcrs > 0, the observation is inefficient and
∗ g− , y + β ∗ g+ ∈ P
(λX, λY ) outperforms (xo , yo ), with xo − βcrs x o crs y crs . It is again possible
10 A Data Envelopment Analysis Toolbox for MATLAB
that further input excesses and outputs shortfalls exist. The slacks corresponding to s− =
xo − βgx− − Xλ, and s+ = Y λ − yo + βgy+ , with s− ≥ 0 and s+ ≥ 0 for any feasible solution
(β, λ). As a result, a second stage is once again needed to calculate these slacks. The next
program incorporating the optimal value βcrs∗ allows determination of these values:
max ω = es− + es+ (6)
λ, s− , s+
subject to
s− = xo − βgx− − Xλ
s+ = Y λ − yo + βgy+
λ ≥ 0, s− ≥ 0, s+ ≥ 0.
As in the input and output oriented models, one may also calculate technical efficiency with
respect toP Pvrs . This requires solving equivalent programs to (5) and (6) adding the vrs
contraint nj=1 λj = 1. If βvrs
∗ = 0 and s−∗ = 0 and s+∗ = 0, the observation is technically
∗ − β∗
efficient. Consequently, we now define scale efficiency as SE = βcrs vrs and directional
efficiency is decomposed into pure technical efficiency, P T E, and the scale efficiency term:
∗ = P T E + SE = β ∗ + SE.
T E = βcrs vrs
The directional model is oriented both in the input and output dimensions, while the choice
of directional vector corresponds to the researcher. Customarily, to keep consistency with the
radial models, the observed amounts of input and output set the direction: g = −gx− , gy+ =
(−xo , yo ), coinciding with the generalized Farrell measure, Briec (1997). In this case it can
be shown that the directional model nests the input and output oriented models. Indeed,
if −gx− , gy+ = (−xo , 0), then β ∗ = 1 − θ∗ , while if −gx− ,
y g + = (0, y ), β ∗ = φ∗ − 1.
o
However, other choices are available; particulary −gx− , gy+ = (−1, 1) or the mean of the
data: −gx− , gy+ =(−x̄o , ȳo ), which are neutral with respect to the orientation, as it does
not use the individual weights corresponding to the observed amounts of inputs and outputs.6
When deciding on the directional model, the researcher must declare whether the direction
corresponds to the observed input or outputs mixes, the unitary vectors, or her own choice of
directional vector. In this case she must introduce the directional input and output matrices.
The directional model can be computed in MATLAB using the dea(X, Y, ...) function with
the orient parameter set to ddf (directional). The input and output directions are specified
in the Gx and Gy parameters as a matrix or as a scalar (usually, 0 or 1). If omitted, X and Y
will be used for Gx and Gy respectively. The returns to scale assumption can be specified by
setting the rts parameter to crs (constant returns to scale; default) or vrs (variable returns
to scale).
> ddf = dea(X, Y, 'orient', 'ddf', 'Gx', X, 'Gy', Y);
> deadisp(ddf);
_______________________________
Data Envelopment Analysis (DEA)
6
Other directions are possible, including elaborated transformations driven by the data as proposed by
Daraio and Simar (2016), or those projecting observations to economic optima as introduced by Zofı́o, Pastor,
and Aparicio (2013), e.g. maximum profit–as shown in section 5.
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 11
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: ddf (Directional distance function)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| Beta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 0.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.2328| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 0.0990| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 0.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.5263| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 0.2857| 5.7143| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 0.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.1379| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.0988| 1.8023| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 0.3333| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 0.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
> ddf_vrs = dea(X, Y, 'orient', 'ddf', 'rts', 'vrs', 'Gx', X, 'Gy', Y);
> deadisp(ddf_vrs);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: ddf (Directional distance function)
Returns to scale: vrs (Variable)
--------------------------------------------------------------------
DMU| X1| X2| Y| Beta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 0.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.1076| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 0.0000| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 0.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.2884| 0.0000| 0.0000| 0.3915|
6| 23.0000| 6.0000| 9.0000| 0.0000| 0.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 0.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.0000| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.0000| 0.0000| 0.0000| 0.0000|
12 A Data Envelopment Analysis Toolbox for MATLAB
10| 42.0000| 25.0000| 26.5000| 0.1629| 0.0000| 5.4560| 0.0000|
11| 5.0000| 17.0000| 12.0000| 0.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
> ddf_scale = deascale(X, Y, 'orient', 'ddf', 'Gx', X, 'Gy', Y);
> deadisp(ddf_scale);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: ddf (Directional distance function)
Returns to scale: scaleeff (Scale efficiency)
---------------------------------
DMU| CRS| VRS| ScaleEff|
---------------------------------
1| 0.0000| 0.0000| 0.0000|
2| 0.2328| 0.1076| 0.1252|
3| 0.0990| 0.0000| 0.0990|
4| 0.0000| 0.0000| 0.0000|
5| 0.5263| 0.2884| 0.2379|
6| 0.2857| 0.0000| 0.2857|
7| 0.0000| 0.0000| 0.0000|
8| 0.1379| 0.0000| 0.1379|
9| 0.0988| 0.0000| 0.0988|
10| 0.3333| 0.1629| 0.1705|
11| 0.0000| 0.0000| 0.0000|
---------------------------------
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 13
3.4. The additive model
The additive model measures technical efficiency based solely on input excesses and output
shortfalls. It does not calculate efficiency scores corresponding to the radial or directional
interpretation of technical efficiency a la Farrell (1957), but characterizes efficiency in terms
of the input and output slacks: s− ∈ Rm and s+ ∈ Rs , respectively. The toolbox implements
the weighted additive program of Lovell and Pastor (1995), whose associated linear program
is:
max ω = ρ− − + +
x s + ρy s (7)
λ, s− , s+
subject to
Xλ + s− =
Y λ− s+ = yo
eλ = 1
λ ≥ 0, s− ≥ 0, s+ ≥ 0,
where (ρ− + m s
x , ρy ) ∈ R+ × R+ are the inputs and outputs weight vectors whose elements can
vary across DMUs. Therefore, assigning unitary values, program (7) corresponds to the
standard additive model, while it is worth noting that it encompasses a wide class of different
DEA models known as General Efficiency Measures (GEMs). Particularly, for the Measure
of Inefficiency Proportions (MIP): (ρ− +
x , ρy ) = (1/xo , 1/yo ); for the range adjusted measure
(RAM): (ρ− , ρ+ ) = (1/(m+s)R− , (1/(m+s)R+ ), where R− and R+ are the variables’ ranges;
while for the Bounded Adjusted Measure (BAM): (ρ− +
x , ρy ) = (1/(m + s)(xo − x), (1/(m +
s)(yo − y)), where x and y are the minimum observed values. Additionally, it is possible to
consider other fixed values with the weights representing value judgements. For observation
(xo , yo ) the objective seeks the maximum feasible reduction in its inputs and increase in
its outputs while remaining in Pvrs . An observation is technically efficient if the optimal
solution (λ∗ , s−∗ , s+∗ ) of the the program is s−∗ = 0, and s+∗ = 0. Otherwise individual input
reductions and output increases would be feasible, and the largest the sum of the slacks, ωvrs ∗ ,
the largest the inefficiency.
The function deaaddit(X, Y, ...) solves the weighted additive model in MATLAB. The
returns to scale assumption can be specified by setting the rts parameter to vrs (variable
returns to scale). Inputs and outputs weights are specified in the rhoX and rhoY parameters.
The default weights correspond to the MIP model if not included.
> add_vrs = deaaddit(X, Y, 'rts', 'vrs');
> deadisp(add_vrs);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: additive
14 A Data Envelopment Analysis Toolbox for MATLAB
Orientation: none
Returns to scale: vrs (Variable)
--------------------------------------------------------------------
DMU| X1| X2| Y| slackX1| slackX2| slackY| Eff|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 0.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.0000| 0.0000| 7.1053| 0.5075|
3| 16.0000| 26.0000| 25.0000| 0.0000| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 0.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.0000| 0.0000| 17.6316| 2.2039|
6| 23.0000| 6.0000| 9.0000| 0.0000| 0.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 0.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.0000| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.0000| 0.0000| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 17.0000| 15.0000| 0.5000| 1.0236|
11| 5.0000| 17.0000| 12.0000| 0.0000| 4.0000| 0.0000| 0.2353|
--------------------------------------------------------------------
For illustration purposes, the Bounded Adjusted Measure (BAM) model can be computed by
specifying the appropriate input and output slacks weights:
> n = size(X, 1);
> m = size(X, 2);
> s = size(Y, 2);
> rhoX = repelem(1 ./ ((m + s) * range(X, 1)), n, 1);
> rhoY = repelem(1 ./ ((m + s) * range(Y, 1)), n, 1);
> add_ram = deaaddit(X, Y, 'rts', 'vrs', 'rhoX', rhoX, 'rhoY', rhoY);
> deadisp(add_ram);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: additive
Orientation: none
Returns to scale: vrs (Variable)
--------------------------------------------------------------------
DMU| X1| X2| Y| slackX1| slackX2| slackY| Eff|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 0.0000| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.0000| 0.0000| 7.1053| 0.1030|
3| 16.0000| 26.0000| 25.0000| 0.0000| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 0.0000| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.0000| 0.0000| 17.6316| 0.2555|
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 15
6| 23.0000| 6.0000| 9.0000| 0.0000| 0.0000| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 0.0000| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.0000| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.0000| 0.0000| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 17.0000| 15.0000| 0.5000| 0.4104|
11| 5.0000| 17.0000| 12.0000| 0.0000| 4.0000| 0.0000| 0.0667|
--------------------------------------------------------------------
3.5. Super-efficiency models
One interesting model that allows to differentiate across technically efficient observations
is that proposed by Andersen and Petersen (1993). While dea assigns the same value
to all efficient units regardless of their performance, the super-efficiency scores allows dis-
criminating across them depending on their values. These scores are obtained by individ-
ually solving for each observation any of the previous models, but excluding them from
the
n reference P data set, which therefore reduces to n o− 1 observations; i.e., P̃crs (xo , yo ) =
(x, y) |x > n−1
Pn−1
j=1,6=o λj xj , y 6 j=1,6=o λj yj , λ > 0 . The magnitude of the super-efficiency
score, as well as the number of observations whose efficiency changes as a result of each in-
dividual exclusion determines the importance of each efficient observation in the complete
dataset.
Therefore, if an observation is inefficient, its super-efficiency score is the same as that previ-
ously calculated, while in the efficient case it is greater than one for the radially input oriented
model, less than one for the output counterpart, and negative for the directional model–as
inputs are to be increase and outputs reduced to reach the reference benchmarks.7
For these oriented models we show in what follows the formulation corresponding to the input
orientation under crs, while its output and directional counterparts are omitted for the sake
of brevity. For the same reason, while the MATLAB toolbox internally solves a two stage
process, the problem can be equivalently expressed according to the following single stage
non-Archimedean formulation:
min θ̃ − ε(1s− + 1s+ ) (8)
θ̃,λ,s− ,s+
subject to
n−1
X
θ̃xo = λj xj + s−
j=1,6=o
n−1
X
yo = λj yj − s+
j=1,6=o
−
λ ≥ 0, s ≥ 0, s+ ≥ 0,
where s− ∈ Rm and s+ ∈ Rs , and is an infinitesimal constant. The optimal solution–obtained
from a first stage equivalent to (1) but excluding the DMU under evaluation from the reference
7
When solving these models under the vrs assumption, it is likely that some unfeasible results are obtained.
16 A Data Envelopment Analysis Toolbox for MATLAB
∗ . In this case, if θ̃ ∗
set–is again denoted by θ̃crs crs > 1, the observation is super-efficient and
the largest the score and the increase in the values of other observation with respect to the
original calculations (1), the most relevant is the removed observation.
The radial super-efficiency model corresponds to the deasuper(X, Y, ...) function in MAT-
LAB with the orient parameter set to the desired orientation (io, oo, or ddf). Once again,
the returns to scale assumption can be specified by setting the rts parameter to crs (constant
returns to scale; default) or vrs (variable returns to scale).
> super = deasuper(X, Y, 'orient', 'io');
> deadisp(super);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial-supereff
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| Theta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.1303| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.6223| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 0.8199| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| 1.1168| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.3104| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 0.5556| 4.4444| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| 1.2449| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.7577| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.8201| 1.6402| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 0.5000| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 1.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
> superddf = deasuper(X, Y, 'orient', 'ddf', 'Gx', X, 'Gy', Y);
> deadisp(superddf);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: directional-supereff
Orientation: ddf (Directional distance function)
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 17
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| Beta| slackX1| slackX2| slackY|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| -0.0612| 0.0000| 0.0000| 0.0000|
2| 16.0000| 12.0000| 14.0000| 0.2328| 0.0000| 0.0000| 0.0000|
3| 16.0000| 26.0000| 25.0000| 0.0990| 0.0000| 0.0000| 0.0000|
4| 17.0000| 15.0000| 26.0000| -0.0552| 0.0000| 0.0000| 0.0000|
5| 18.0000| 14.0000| 8.0000| 0.5263| 0.0000| 0.0000| 0.0000|
6| 23.0000| 6.0000| 9.0000| 0.2857| 5.7143| 0.0000| 0.0000|
7| 25.0000| 10.0000| 27.0000| -0.1091| 0.0000| 0.0000| 0.0000|
8| 27.0000| 22.0000| 30.0000| 0.1379| 0.0000| 0.0000| 0.0000|
9| 37.0000| 14.0000| 31.0000| 0.0988| 1.8023| 0.0000| 0.0000|
10| 42.0000| 25.0000| 26.5000| 0.3333| 0.0000| 0.0000| 0.0000|
11| 5.0000| 17.0000| 12.0000| 0.0000| 0.0000| 4.0000| 0.0000|
--------------------------------------------------------------------
As for the super-efficiency calculations in the additive model, we follow Du, Liang, and Zhu
(2010), and solve the corresponding counterpart to (7):
min ω = es− + es+ (9)
λ, s− , s+
subject to
n−1
X
xo ≥ λj xj − s−
j=1,6=o
n−1
X
yo ≤ λj yj + s+
j=1,6=o
eλ = 1
λ ≥ 0, s− ≥ 0, s+ ≥ 0.
The constraints in the program are modified as inputs need to be increased and outputs
reduced so as to reach the production possibility set. In this occasion an observation is super-
efficient if the optimal solution (λ∗ , s−∗ , s+∗ ) yields positive slacks, and the largest their sum,
the largest the super-efficiency.
The additive super-efficiency model can be computed in MATLAB using the deaadditsuper(X,
Y, ...) function, and specifying returns to scale by setting the rts parameter to crs (con-
stant returns to scale; default) or vrs (variable returns to scale). Inputs and outputs weights
are specified in the rhoX and rhoY parameters. The default weights correspond to the MIP
model if not included.
> additsuper = deaadditsuper(X, Y);
> deadisp(additsuper);
18 A Data Envelopment Analysis Toolbox for MATLAB
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: additive-supereff
Orientation: none
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X1| X2| Y| slackX1| slackX2| slackY| Eff|
--------------------------------------------------------------------
1| 5.0000| 13.0000| 12.0000| 1.1298| 0.0000| 0.0000| 1.1298|
2| 16.0000| 12.0000| 14.0000| NaN| NaN| NaN| NaN|
3| 16.0000| 26.0000| 25.0000| NaN| NaN| NaN| NaN|
4| 17.0000| 15.0000| 26.0000| 0.0000| 0.0000| 2.7200| 2.7200|
5| 18.0000| 14.0000| 8.0000| NaN| NaN| NaN| NaN|
6| 23.0000| 6.0000| 9.0000| NaN| NaN| NaN| NaN|
7| 25.0000| 10.0000| 27.0000| 0.0000| 3.8713| 0.0000| 3.8713|
8| 27.0000| 22.0000| 30.0000| NaN| NaN| NaN| NaN|
9| 37.0000| 14.0000| 31.0000| NaN| NaN| NaN| NaN|
10| 42.0000| 25.0000| 26.5000| NaN| NaN| NaN| NaN|
11| 5.0000| 17.0000| 12.0000| NaN| NaN| NaN| NaN|
--------------------------------------------------------------------
4. Productivity change: The Malmquist index
The Malmquist index, introduced by Caves, Christensen, and Diewert (1982), measures the
change in productivity of the observation under evaluation by comparing its relative per-
formance with respect to reference technologies corresponding to two different time periods.
The constant returns to scale production possibility set corresponding to t = 1,. . . ,T peri-
ods defines as Pcrs t = (x, y) |x > X t λ, y 6 Y t λ, λ > 0 , with the variable returns to scale
counterpart denoted accordingly by Pvrs t . The standard Malmquist index relies solely on the
concept of radial efficiency and requires calculation of the input–or output–oriented scores
of (xto , yot ) observed in two periods t = 1, 2, with respect to the constant returns to scale
reference frontier of any period. Taking as the reference the first period, Pcrs1 , we denote both
1,1 2,1
scores by θcrs and θcrs , where the first superscript refers to the time period of the observation
1,1
and the second one to that of the reference technology. While θcrs is the solution to program
2,1
(1), the intertemporal score θcrs corresponds to the following program that evaluates period
2 observation (x2o , yo2 ) with respect to period 1 technology:
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 19
min θ (10)
θ,λ
subject to
θx2o ≥ X 1 λ
Y 1 λ ≥ yo2
λ ≥ 0.
1,2
Equivalently, the analogous intertemporal score θcrs using the second period technology as
reference corresponds to the same program but reversing the time superscripts of the firm
under evaluation from (x2o , yo2 ) to (x1o , yo1 ), and those of the reference technology from (X 1 , Y 1 )
to (X 2 , Y 2 ).
After the contemporary and intertemporal efficiency scores have been calculated it is pos-
2,1 1,1 2,2 1,2
sible to define the following Malmquist indices: M1 = θcrs /θcrs and M2 = θcrs /θcrs . For
both indices, if M > 1 there is productivity increase, while if M = 1 productivity remains
unchanged and M < 1 signals productivity decline. Following Färe, Grosskopf, Norris, and
Zhang (1994), productivity change can be decomposed into efficiency change and techni-
cal change.8 The former corresponds to the so-called catch-up effect; i.e., the change in
the technical efficiency of the observation under evaluation with respect to the two peri-
2,2 1,1
ods, which defines for both indices as M T EC = θcrs /θcrs . The latter corresponds to the
frontier-shift effect, i.e. the change in the reference frontier between both periods, which
2,1 2,2
defines for M1 as M T C1 = (θcrs /θcrs ) using period 2 observation as the reference benchmark
1,1 1,2
to evaluate the shift in the frontier. For M2 it defines a M T C2 = (θcrs /θcrs ). As in the
previous cases, if M T EC > 1 or M T C > 1 productivity change is driven respectively by
both technical efficiency gains and technical change improvements (technical progress); while
M T EC < 1 or M T C < 1 imply lower productivity due to greater inefficiency and technical
regress. Finally, unitary values signal that the technical efficiency and the reference frontier
remain unchanged. Therefore the decomposition of year-to-year productivity change defines
2,2 1,1 2,1 2,2
as M1 = M T EC × M T C1 = θcrs /θcrs × θcrs /θcrs –and similarly for M2 , while the geomet-
ric mean of both indices and their corresponding decompositions represents a compromise
2,2 1,1 2,1 2,2 1,1 1,2 (1/2)
between both values: M = M T EC × M T C = θcrs /θcrs × (θcrs /θcrs × θcrs /θcrs ) .
Finally, it is also possible to decompose long term productivity change from an initial period
to a final period into consecutive subperiods relying on the transitivity property of index
numbers. This allows the analysis of productivity change by subperiods. Therefore, given a
sequence of periods, i.e., t = 1, 2, 3, it is verified that the Malmquist index between the base
and final periods can be expressed in terms of its chain components: M 1,3 = M 1,2 × M 2,3 .
The toolbox calculates the sequence of year-to-year Malmquist indices and the cumulated
Malmquist index taking the based period as reference.
To compute the Malmquist indices in MATLAB one calls the deamalm(X, Y, ...) function
with the orient parameter set to any of the two orientations (io or oo). The parameters X
8
The toolbox calculates a first level decomposition that does not take into account the contribution that
scale efficiency change or returns to scale make to productivity change. Therefore, as productivity change is
measured under the crs assumptions, the input and output orientations are equal. To account for scale effects
alternative second level decompositions have been proposed in the literature, see Zofı́o (2007). These terms
can be calculated solving the necessary models under vrs.
20 A Data Envelopment Analysis Toolbox for MATLAB
and Y must be 3D arrays, with the third dimension corresponding to different time periods.
By default, Malmquist indices are computed as the geometric mean of M1 and M2 . However,
for simplicity, the indices can be computed taking the first period as the base for technical
change, M1 , by setting the parameter geomean to 0.
> X = [2; 3; 5; 4; 4];
> X(:, :, 2) = [1; 2; 4; 3; 4];
> X(:, :, 3) = [0.5; 1.5; 3; 2; 4];
> Y = [1; 4; 6; 3; 5];
> Y(:, :, 2) = [1; 4; 6; 3; 3];
> Y(:, :, 3) = [ 2; 4; 6; 3; 1];
> malmquist = deamalm(X, Y, 'orient', 'io');
> deadisp(malmquist)
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 5
Inputs: 1 Outputs: 1
Model: radial-malmquist
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
Malmquist:
Base period is previous period
-----------------------------------------------------------
DMU| M1| M2| MTEC1| MTEC2| MTC1| MTC2|
-----------------------------------------------------------
1| 2.0000| 4.0000| 1.3333| 2.0000| 1.5000| 2.0000|
2| 1.5000| 1.3333| 1.0000| 0.6667| 1.5000| 2.0000|
3| 1.2500| 1.3333| 0.8333| 0.6667| 1.5000| 2.0000|
4| 1.3333| 1.5000| 0.8889| 0.7500| 1.5000| 2.0000|
5| 0.6000| 0.3333| 0.4000| 0.1667| 1.5000| 2.0000|
-----------------------------------------------------------
M = Malmquist. MTEC = Technical Efficiency Change. MTC = Technical Change.
The aggregate Malmquist index taking the first period as the base period can be computed
by setting the fixbaset parameter to 1.
> malmquist1 = deamalm(X, Y, 'orient', 'oo', 'fixbaset', 1);
> deadisp(malmquist1)
_______________________________
Data Envelopment Analysis (DEA)
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 21
DMUs: 5
Inputs: 1 Outputs: 1
Model: radial-malmquist
Orientation: oo (Output oriented)
Returns to scale: crs (Constant)
Malmquist:
Base period is period 1
-----------------------------------------------------------
DMU| M1| M2| MTEC1| MTEC2| MTC1| MTC2|
-----------------------------------------------------------
1| 2.0000| 8.0000| 1.3333| 2.6667| 1.5000| 3.0000|
2| 1.5000| 2.0000| 1.0000| 0.6667| 1.5000| 3.0000|
3| 1.2500| 1.6667| 0.8333| 0.5556| 1.5000| 3.0000|
4| 1.3333| 2.0000| 0.8889| 0.6667| 1.5000| 3.0000|
5| 0.6000| 0.2000| 0.4000| 0.0667| 1.5000| 3.0000|
-----------------------------------------------------------
M = Malmquist. MTEC = Technical Efficiency Change. MTC = Technical Change.
5. Allocative models: Economic efficiency
Assuming an economic optimizing behaviour on the part of the observations, e.g., cost min-
imization, revenue maximization, or profit maximization, and the corresponding input and
output prices: w ∈ Rm s
++ and p ∈ R++ , it is possible to measure their economic efficiency and,
based on duality theory, decompose it into the technical efficiency and allocative efficiency
terms, Farrell (1957).
5.1. Cost efficiency
Since we are concerned with overall efficiency in the input space, we assume that obser-
vations minimize production costs. This implies that if observations succeed in using the
inputs mix (bundle) resulting in the minimum cost of producing a given output level at
the existing market prices, they are cost efficient. Let us denote by C (y, w) the mini-
mum cost of producing the output level y given the input price vector w: C (y, w) =
m
P
min wi xi |x > Xλyo 6 Y λ, λ > 0 , which considers the input possibility set capable of
i=1
producing yo .9 For the observed outputs levels we can calculate minimum cost and the asso-
ciated optimal quantities of inputs x∗ consistent with the production technology by solving
the following program:
9
For a recent discussion on the properties of the technology when decomposing economic efficiency into
technical and allocative efficiencies see Aparicio, Pastor, and Zofio (2015).
22 A Data Envelopment Analysis Toolbox for MATLAB
min C (y, w) = wx∗ (11)
x,λ
subject to
x ≥ Xλ
Y λ ≥ yo
λ ≥ 0.
Once minimum cost is calculated, cost efficiency defines as the ratio of minimum cost to
observed cost: CE = C (y, w) /wxo . Thanks to duality results—Shephard(1953), CE can
be decomposed into the technical efficiency measure associated to (1), θcrs ∗ , and the residual
difference corresponding to the allocative cost efficiency: AE = CE/θcrs ∗ . Therefore allocative
efficiency defines as the ratio between minimum cost and production cost at the technically
efficient projection of the observation: AE = C (y, w) /wx̂o –with x̂o = θcrs ∗ x and CE =
T E × AE. Consequently, if CE < 1 and the observation is technically efficient, θcrs ∗ = 1,
all cost inefficiency is allocative, while if the observation uses the right proportions of input
quantities: x̂o = x∗ , it is allocative efficient and AE = 1.
The cost efficiency model can be computed in MATLAB using the deaalloc(X, Y, ...)
function. The parameter Xprice with inputs’ prices, as a matrix or as a row vector, must be
included.10
> X = [3 2; 1 3; 4 6];
> Y = [3; 5; 6];
> W = [4 2];
> P = 6;
> cost = deaalloc(X, Y, 'Xprice', W);
> deadisp(cost);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 3
Inputs: 2 Outputs: 1
Model: allocative-cost
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
------------------------------------------------------------------------------
DMU| X1| X2| Xprice1| Xprice2| Y| TechEff| AllocEff| CostEff|
------------------------------------------------------------------------------
1| 3.0000| 2.0000| 4.0000| 2.0000| 3.0000| 0.9000| 0.4167| 0.3750|
10
If input prices are the same for all DMUs the toolbox generates a matrix with the same values. If inputs’
prices differ between DMUs a matrix with individual price information can be supplied. On this occasion we
illustrate the economic models with an example from (Cooper et al. 2007, 261).
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 23
2| 1.0000| 3.0000| 4.0000| 2.0000| 5.0000| 1.0000| 1.0000| 1.0000|
3| 4.0000| 6.0000| 4.0000| 2.0000| 6.0000| 0.6000| 0.7143| 0.4286|
------------------------------------------------------------------------------
5.2. Revenue efficiency
In an equivalent way, from an output dimension one may be interested in revenue efficiency.
Now if observations are able to produce the output mix (bundle) resulting in maximum
revenue given their inputs levels at the existing market prices, they are revenue efficient. Let
us denote by R (x, p) the maximum sfeasible revenue using inputs’ levels
x and given the
P
outputs’ prices p: R (x, p) = max pi yi |xo > Xλ, y 6 Y λ, λ > 0 ; i.e., considering the
i=1
output possibility set producible with xo .
In this case, we calculate maximum revenue along with the optimal output quantities y∗ by
solving the following program:
max R (xo , p) = py∗ (12)
y,λ
subject to
xo ≥ Xλ
Yλ ≥y
λ ≥ 0.
Revenue efficiency defines as the ratio of observed revenue to maximum revenue: RE =
pyo /R (x, p). Again, duality results from an output perspective allow us to decompose RE
into the output technical efficiency measure associated to the inverse of (3), T E = 1/φ∗crs , and
the residual difference corresponding to the allocative revenue efficiency: AE = RE × φ∗crs .
In this occasion allocative efficiency defines as the ratio between revenue at the technically
efficient projection of the observation and maximum revenue: AE = pŷo /R (x, p)–with
ŷo = φ∗crs yo , and therefore RE = T E × AE. Now, if RE < 1 and the observation is
technically efficient, φ∗crs = 1, all revenue inefficiency is allocative and the observation is
unable to produce the right output mix. Contrarily, if the observation produces the right
proportion of output quantities given their prices: ŷo = y∗ , it is allocative efficient: AE = 1.
The revenue efficiency model can be computed in MATLAB using the deaalloc(X, Y, ...)
function. The parameter Yprice with the ouputs’ prices, as a matrix or as a row vector, must
be included.
> revenue = deaalloc(X, Y, 'Yprice', P);
> deadisp(revenue);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 3
Inputs: 2 Outputs: 1
24 A Data Envelopment Analysis Toolbox for MATLAB
Model: allocative-revenue
Orientation: oo (Output oriented)
Returns to scale: crs (Constant)
---------------------------------------------------------------------
DMU| X1| X2| Y| Yprice| TechEff| AllocEff| RevEff|
---------------------------------------------------------------------
1| 3.0000| 2.0000| 3.0000| 6.0000| 0.9000| 1.0000| 0.9000|
2| 1.0000| 3.0000| 5.0000| 6.0000| 1.0000| 1.0000| 1.0000|
3| 4.0000| 6.0000| 6.0000| 6.0000| 0.6000| 1.0000| 0.6000|
---------------------------------------------------------------------
5.3. Profit efficiency: The directional approach
Profit efficiency allows studying the economic behavior of observations as profit maximizers.
In this case, if observations are able to maximize the difference between revenue and cost given
market prices by producing and using the right quantities of outputs and inputs respectively,
nP s Pm
they are profit efficient. The profit function defines as Π (w, p) = max p i yi − wi x i |
o i=1 i=1
x > Xλ, y 6 Y λ, eλ = 1, λ > 0 .
Calculating maximum revenue along with the optimal output and input quantities y∗ and x∗
requires solving:
max Π (w, p) = py∗ − wx∗ (13)
x,y,λ
subject to
x ≥ Xλ = x
y ≤Yλ=y
eλ = 1
λ ≥ 0.
Profit efficiency defines as the difference between maximum profit and observed profit. As
duality results concerning its decomposition into technical and allocative profit efficiencies rely
on the directional model (5), we can use the directional vector to define a normalized profit
efficiency measure that is homogenous of degree one in prices–invariant to proportional price
changes; i.e., P E = (Π (w, p) − (pyo − wxo ))/(pgy+ + wgx− )–see Chambers, Chung, and Färe
(1998) and Zofı́o et al. (2013) for more recent proposals. Subsequently, profit efficiency can be
decomposed into the directional technical efficiency measure associated to (5) under variable
returns to scale, T E = βvrs∗ , and the residual difference corresponding to the allocative profit
efficiency term: AE = P E − βvrs ∗ . Therefore, allocative efficiency defines as the difference
between maximum profit and profit at the technically efficient projection of the observation:
AE = (Π (w, p) − (pŷo − wx̂o ))/(pgy+ + wgx− ), with (ŷo , x̂o ) = (yo + βvrs
∗ g+ , x − β ∗ g− ),
y o vrs x
and P E = T E + AE. Now, if P E > 0 and the observation is technically efficient, βvrs ∗ = 0,
all profit efficiency is allocative, and the observation is unable to produce with the optimal
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 25
combination of outputs and inputs. Contrarily, if the observation supplies and demands the
profit maximizing quantities of outputs and inputs, (ŷo , x̂o ) = (y∗ , x∗ ), it is allocative efficient:
AE = 0.
The profit efficiency model can be computed in MATLAB using the deaalloc(X, Y, ...)
function. Both the parameters Xprice and Yprice with inputs and outputs prices must be
included. The input and output directions are specified in the Gx and Gy parameters. If
omitted, X and Y will be used for Gx and Gy respectively.
> profit = deaalloc(X, Y, 'Xprice', W, 'Yprice', P);
> deadisp(profit, 'names/X/Y/eff.T/eff.A/eff.P');
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 3
Inputs: 2 Outputs: 1
Model: allocative-profit
Orientation: ddf (Directional distance function)
Returns to scale: crs (Constant)
------------------------------------------------------------
DMU| X1| X2| Y| TechEff| AllocEff| ProfEff|
------------------------------------------------------------
1| 3.0000| 2.0000| 3.0000| 0.0000| 0.5294| 0.5294|
2| 1.0000| 3.0000| 5.0000| 0.0000| -0.0000| -0.0000|
3| 4.0000| 6.0000| 6.0000| 0.0000| 0.1875| 0.1875|
------------------------------------------------------------
5.4. Profit efficiency: The weighted additive approach
An alternative decomposition of profit efficiency based on the weighted additive measure and
duality results is possible, Cooper, Pastor, Aparicio, and Borras (2011). In contrast to the
radial and directional approaches, the technical efficiency component of the decomposition
accounts for all inefficiencies (i.e., it includes individual inputs and outputs slacks).
Profit efficiency defines in this case as the difference between maximum profit and observed
profit, normalized by the minimum among the ratios of of market prices to weights. i.e., P E =
−
(Π (w, p)−(pyo −wxo ))/min{(p1 /ρ+ + −
1 ), ..., (pm /ρs ), (w1 /ρ1 ), ..., (ρ1 /ρm )}. This normalization
ensures once again that the profit efficiency is homogenous of degree one in prices, rendering
it units invariant.
Once the weighted additive measure (7) is calculated, profit efficiency can be decomposed
into technical efficiency, T E = ωvrs∗ , and the residual difference corresponding to the alloca-
tive profit efficiency term: AE = P E − ωvrs ∗ . Consequently, allocative efficiency defines as the
difference between maximum profit and profit at the technically efficient projection of the ob-
−
servation: AE = (Π (w, p) − (pŷo − wx̂o ))/min{(p1 /ρ+ + −
1 ), ..., (pm /ρs ), (w1 /ρ1 ), ..., (ρ1 /ρm )},
+ + − −
with (ŷo , x̂o ) = (yo + ρy s , xo − ρx s ), and once again P E = T E + AE. Now, if P E > 0
∗
and the observation is technically efficient, ωvrs = 0, all profit efficiency is allocative, and the
26 A Data Envelopment Analysis Toolbox for MATLAB
observation is unable to produce with the optimal combination of outputs and inputs. Con-
trarily, if the observation supplies and demands the profit maximizing quantities of outputs
and inputs, (ŷo , x̂o ) = (y∗ , x∗ ), it is allocative efficient: AE = 0.
In this section we have focused on the decomposition of profit efficiency, while the cost and
revenue counterparts presented in the previous sections can be easily obtained by setting
outputs and inputs weights equal to zero, and using cost and revenue as support functions,
respectively.
The weighted additive profit efficiency model is computed with the deaadditprofit(X, Y,
...) function. Both the parameters Xprice and Yprice with inputs and outputs prices must
be included. The returns to scale assumption, rts parameter, must be set to to vrs (variable
returns to scale). Inputs and outputs weights are specified in the rhoX and rhoY parameters.
As before, the default weights correspond to the MIP model if not included.
> addprofit = deaadditprofit(X, Y, 'rts', 'vrs', 'Xprice', W, 'Yprice', P);
> deadisp(addprofit, 'names/X/Y/eff.T/eff.A/eff.P');
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 3
Inputs: 2 Outputs: 1
Model: additive-profit
Orientation: none
Returns to scale: vrs (Variable)
------------------------------------------------------------
DMU| X1| X2| Y| TechEff| AllocEff| ProfEff|
------------------------------------------------------------
1| 3.0000| 2.0000| 3.0000| 0.0000| 4.5000| 4.5000|
2| 1.0000| 3.0000| 5.0000| 0.0000| -0.0000| -0.0000|
3| 4.0000| 6.0000| 6.0000| 0.0000| 3.0000| 3.0000|
------------------------------------------------------------
6. Undesirable outputs
Along with desirable and market oriented products, observations may produce undesirable
or detrimental outputs as byproducts, such as pollutants or hazardous wastes from an en-
vironmental perspective. As a result, efficiency and productivity measures that do not take
into account the asymmetry between both types of production: desirable and undesirable,
will result in biased assessments of performance and erroneous calculations; e.g., when assess-
ing environmental performance and making recommendations to improve technical efficiency,
one seeks increments in desirable production but reductions in undesirable outputs. To in-
corporate undesirable outputs into the efficiency and productivity change models, we rely
on the directional measures (5) that treat both sets of outputs differently. This requires
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 27
a redefinition of the production technology where the initial vector of i = 1, 2, . . . , s out-
puts y ∈ Rs++ is partitioned into desirable and undesirable production, i.e., y = (yd , yu ),
with yd ∈ Rq++ and yu ∈ Rr++ , respectively. This translates into the corresponding reference
technology Pcrs = x, yd , yu |x > Xλ, yd 6 Y λ, yu = Y λ, λ > 0 , characterizing unde-
sirable outputs as weakly disposable, (Cooper et al. 2007, ch. 12). As the directional ef-
ficiency measure resulting in increases in desirable outputs and reductions in undesirable
outputs relative to the same amount of inputs is used in the following section to define the
Malmquist-Luenberger productivity measure, we rely on Aparicio et al. (2015) contribution
to calculate it, thereby preventing the inconsistencies of the original approach introduced by
Chung, Färe, and Grosskopf (1997). In this case, the directional efficiency measure projecting
observation xo , yod , you along the pre-assigned direction corresponding to the output vector
gy = yd , yu 6= 0m+s , corresponds to the solution of the following program:
max β (14)
β,λ
subject to
Xλ ≤ xo
Y d λ ≥ yod + βyod
Y u λ ≤ you − βyou
max {yiu } ≥ you − βyou
λ ≥ 0.
∗ , and if β ∗
Again, the optimal solution corresponds to βcrs crs = 0, with λo = 1, λj = 0 (j 6=
o), the observation is directional efficient. Otherwise, βcrs∗ > 0 signals inefficiency and
d u d u
(λX, λY , λY ) outperforms xo , yo , yo . It is possible also to calculate the non-directional
slacks, checking for further desirable output shortfalls or inputs and undesirable output ex-
cesses.
The undersirable outputs model is computed in MATLAB using the deaund(X, Y, Yu, ...)
function, where Y is now reserved for the matrix of desirable outputs–i.e., from a computational
perspective we drop the superscript d–, and Yu is the matrix of undesirable outputs.
> X0 = ones(5,1);
> Y0 = [7; 5; 1; 3; 4];
> Yu0 = [2; 5; 3; 3; 2];
> undesirable = deaund(X0, Y0, Yu0);
> deadisp(undesirable);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 5
Inputs: 1 Outputs: 1 Undesirable: 1
28 A Data Envelopment Analysis Toolbox for MATLAB
Model: directional-undesirable
Orientation: ddf (Directional distance function)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
DMU| X| Y| Yu| Beta| slackX| slackY| slackYu|
--------------------------------------------------------------------
1| 1.0000| 7.0000| 2.0000| 0.0000| 0.0000| 0.0000| 0.0000|
2| 1.0000| 5.0000| 5.0000| 0.4000| 0.0000| 0.0000| 1.0000|
3| 1.0000| 1.0000| 3.0000| 0.8261| 0.7391| 0.0000| 0.0000|
4| 1.0000| 3.0000| 3.0000| 0.5556| 0.3333| 0.0000| 0.0000|
5| 1.0000| 4.0000| 2.0000| 0.2727| 0.2727| 0.0000| 0.0000|
--------------------------------------------------------------------
7. Productivity change: The Malmquist-Luenberger index
Finally, it is possible to define the undesirable outputs productivity counterpart to the
Malmquist index. Chung et al. (1997) introduced the Malmquist-Luenberger (M L) pro-
ductivity index measuring the change in productivity of the observation under evaluation by
comparing its relative efficiency with respect to reference technologies corresponding to two
different time periods. As its radial counterpart, the standard M L uses the directional effi-
ciency score only–disregarding
second stage slacks–and requires calculation of the mix period
efficiency of xto , yot,d , yot,u observed in two periods t = 1, 2. We denote both scores by βcrs1,1
2,1
and βcrs , where once again the first superscript refers to the time period of the observation
1,1
and the second one to that of the reference technology. While βcrs is the solution to program
2,1
(14), the intertemporal
score βcrs corresponds to the following program that evaluates period
2 observation xo , yo , yo2,u with respect to period 1 technology:
2 2,d
max β (15)
β,λ
subject to
X 1 λ ≤ x2o
Y 1,d λ ≥ yo2,d + βyo2,d
Y 1,u λ ≤ yo2,u − βyo2,u
max {yit,u } ≥ yo2,u − βyo2,u
λ ≥ 0,
where max {yit,u } is the maximum observed amount of each undesirable output in all time
1,2
periods, Aparicio, Pastor, and Zofio (2013). Correspondingly, the score βcrs , using the second
period technology as reference can be calculated solving an equivalent program that evaluates
(x1o , yo1,d , yo1,u ) with respect to the reference technology (X 2 , Y 2,d , Y 2,u ).
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 29
1,1 2,1
Once these scores are calculated the M L indices define as follows: M L1 = (1+βcrs )/(1+βcrs )
1,2 2,2
and M L2 = (1 + βcrs )/(1 + βcrs ). If M L > 1 efficiency increases and the evaluated unit is
capable of producing more desirable output with less undesirable production, while if M L = 1
productivity remains unchanged, and M L < 1 signals productivity decline. Following Chung
et al. (1997) the index can be decomposed into efficiency change and technical change, which
have the same interpretation that their previous Malmquist counterparts. The change in
1,1 2,2
technical efficiency defines now as M LT EC = (1 + βcrs )/(1 + βcrs ), while the frontier-shift
2,2 2,1
effect corresponding to technical change are M LT C1 = (1 + βcrs )/(1 + βcrs ) and M LT C2 =
1,2 1,1
(1 + βcrs )/(1 + βcrs ). Again, if M LT EC > 1 or M LT C > 1 productivity change responds
to both technical efficiency gains and technical change improvements (technical progress);
while M LT EC < 1 or M LT C < 1 imply lower productivity with greater inefficiency and
technical regress. Finally, unitary values signal that the technical efficiency and the reference
frontier remain unchanged. Therefore the decomposition of productivity change defines as
1,1 2,2 2,2 2,1
M L1 = M LT EC × M LT C1 = (1 + βcrs )/(1 + βcrs ) × (1 + βcrs )/(1 + βcrs )–and similarly for
M L2 . Given a sequence a of years, the toolbox calculates the year-to-year variation of M L1
1,1 2,2
and the geometric mean of both indices: M L = M LT EC × M LT C = (1 + βcrs )/(1 + βcrs )×
2,2 2,1 1,2 1,1 (1/2)
((1 + βcrs )/(1 + βcrs ) × ((1 + βcrs )/(1 + βcrs )) .
The Malmquist-Luenberger indices are computed in MATLAB calling the deamalmluen(X,
Y, Yu,...). The parameters X, Y and Yu must be 3D arrays, with the third dimension
corresponding to different time periods. By default, Malmquist-Luenberger indices are com-
puted on a year-to-year basis taking the previous year as reference period. The aggregate
Malmquist-Luenberger index taking the first period as the base period can be computed by
setting the fixbaset parameter to 1. By default, Malmquist indices are computed as the
geometric mean of M1 and M2 . However, for simplicity, the indices can be computed taking
the first period as the base for technical change, M1 , by setting the parameter geomean to 0.
> X1 = ones(5,1);
> Y1 = [8; 5.5; 2; 2; 4];
> Yu1 = [1; 3; 2; 4; 1];
> X = X0;
> X(:, :, 2) = X1;
> Y = Y0;
> Y(:, :, 2) = Y1;
> Yu = Yu0;
> Yu(:, :, 2) = Yu1;
> ml = deamalmluen(X, Y, Yu);
> deadisp(ml);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 5
Inputs: 1 Outputs: 1 Undesirable: 1
Model: directional-malmquist-luenberger
Orientation: ddf (Directional distance function)
30 A Data Envelopment Analysis Toolbox for MATLAB
Returns to scale: crs (Constant)
Malmquist-Luenberger:
Base period is previous period
--------------------------------
DMU| ML| MLTEC| MLTC|
--------------------------------
1| 1.3702| 1.0000| 1.3702|
2| 1.1000| 0.9625| 1.1429|
3| 1.1260| 1.0272| 1.0962|
4| 0.9162| 0.8264| 1.1087|
5| 1.2792| 0.9545| 1.3401|
--------------------------------
ML: Malmquist-Luenberger. MLTEC: Technical Efficiency Change.
MLTC: Technical Change.
8. Bootstrapping DEA estimators
The efficiency scores and productivity indices calculated in the previous sections can be con-
sidered as estimates as their values are subject to uncertainty due to sampling variation.
Simar and Wilson (1998) introduce bootstrap methods that, based on resampling, provide es-
timation bias, confidence intervals and allow hypotheses testing. This toolbox implements the
algorithms presented by these authors following Bogetoft and Otto (2011), who also discuss
other non-parametric and asymptotic tests. We implement the tests that allow determining
the significance of the models covered by the toolbox: The assumption about returns to scale,
i.e., whether constant and variable returns to scale scores are significantly different–Simar
and Wilson (2002), and the existence of productivity change over time, i.e., whether the
Malmquist index and its components are different from unity, Simar and Wilson (1999).
The algorithm implemented in the toolbox is the following:
1. Selection of B independent bootstrap samples–drawn from the original dataset with
replacements;
2. Calculate an initial estimate for the efficiency score of each DMU with respect to each
bootstrapped sample and smooth their distribution by perturbating them with a random
noise generated from from a kernel density function with scale given with bandwith h;
3. Correct the original estimates for the mean and variance of the smoothed values;
4. Obtain a second set of bootstrapped samples generating inefficient DMUs, inside the
DEA attainable set and conditional on the original input–or output–mix;
5. Repeat the process, estimate the efficiency scores for each original DMU with respect
to that second set, so as to obtain a set of B bootstrap estimates; and finally,
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 31
6. Based on this distribution calculate the threshold values that truncate it according to
the pre-determined significance value α, so as to determine the confidence intervals for
the efficiency score of each DMU.
In addition, the bootstrapped scores can be used to obtain an estimate of the bias of the true
efficiency value, and thereby a bias-corrected estimator. 11
The radial bootstrap model can be computed in MATLAB using the deaboot(X, Y, ...)
function with the orient parameter set to the desired orientation (io for input oriented or
oo for output oriented). The returns to scale assumption can be specified by setting the rts
parameter to crs (constant returns to scale; default) or vrs (variable returns to scale). The
number of bootstrap replication are set in the nreps parameter (which defaults to 200), and
the significance level in the alpha parameter (0.05).
> load 'deadataFLS'
> rng(1234567); % Set seed for reproducibility
> io_b = deaboot(X, Y, 'orient', 'io', 'nreps', 200, 'alpha', 0.05);
deadisp(io_b);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial-bootstrap
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
Bootstrap replications: 200
Significance level: 0.05
-----------------------------------------
DMU| eff| effboot| effCI1| effCI2|
-----------------------------------------
1| 1.0000| 0.8245| 0.7024| 0.9853|
2| 0.6223| 0.5606| 0.4949| 0.6150|
3| 0.8199| 0.7123| 0.6289| 0.8064|
4| 1.0000| 0.8864| 0.7825| 0.9841|
5| 0.3104| 0.2793| 0.2468| 0.3063|
6| 0.5556| 0.4846| 0.4083| 0.5527|
7| 1.0000| 0.8473| 0.7345| 0.9814|
8| 0.7577| 0.6794| 0.5989| 0.7513|
9| 0.8201| 0.7008| 0.6024| 0.8159|
10| 0.5000| 0.4460| 0.3924| 0.4936|
11| 1.0000| 0.8461| 0.7024| 0.9888|
-----------------------------------------
11
A similar procedure is implemented to determine the significance of the Malmquist indices–see Simar and
Wilson (1999) for the specific algorithm.
32 A Data Envelopment Analysis Toolbox for MATLAB
The returns to scale test is compute by calling the deatestrts(X, Y, ...) function with
the orient parameter set to the desired orientation (io for input oriented or oo for output
oriented). The number of bootstrap replication are set in the nreps parameter (which defaults
to 200), and the significance level in the alpha parameter (0.05). Results of the test can be
displayed on screen by setting the parameter disp to 1.
> rng(1234567); % Set seed for reproducibility
> deatestrts(X, Y, 'orient', 'io', 'nreps', 200, 'alpha', 0.05, 'disp', 1);
_______________________________
DEA Test of RTS
H0: Globally CRS
H1: VRS
Bootstrap replications: 200
Significance level: 0.05
S statistic: 0.8318
Critical value: 0.7302
p-value: 0.2900
To compute bootstrapped Malmquist indices in MATLAB one calls the deamalmboot(X, Y,
...) function with the orient parameter set to input orientation (io). The parameters X
and Y must be 3D arrays, with the third dimension corresponding to different time periods.
The number of bootstrap replication are set in the nreps parameter (which defaults to 200),
and the significance level in the alpha parameter (0.05). Again, by default, Malmquist indices
are computed as the geometric mean of M1 and M2 . However, for simplicity, the indices can
be computed taking the first period as the base for technical change, M1 , by setting the
parameter geomean to 0.
> X = [2; 3; 5; 4; 4];
> X(:, :, 2) = [1; 2; 4; 3; 4];
> Y = [1; 4; 6; 3; 5];
> Y(:, :, 2) = [1; 4; 6; 3; 3];
> rng(1234567); % Set seed for reproducibility
> malmquist = deamalmboot(X, Y, 'orient', 'io', 'nreps', 200, 'alpha', 0.05);
deadisp(malmquist)
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 5
Inputs: 1 Outputs: 1
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 33
Model: radial-malmquist-bootstrap
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
Bootstrap replications: 200
Significance level: 0.05
Base period is previous period
-----------------------------------------
DMU| M| Mboot| McLow| McUpp|
-----------------------------------------
1| 2.0000| 1.7355| 1.3193| 2.1003|
2| 1.5000| 1.5191| 1.2958| 1.7594|
3| 1.2500| 1.3943| 1.2383| 1.5831|
4| 1.3333| 1.4445| 1.2798| 1.6690|
5| 0.6000| 0.8266| 0.7780| 0.8941|
-----------------------------------------
9. Advanced options, displaying and exporting results
9.1. Specifying DMU names
When using the functions to compute DEA models, by default observations are numbered
from 1 to n. However, custom names can be assigned to DMUs by adding the parameter
names and a cell string with the desired names when calling the functions. Using data from
the example in section 6:
> X0 = ones(5,1);
> Y0 = [7; 5; 1; 3; 4];
> Yu0 = [2; 5; 3; 3; 2];
> names = {'A','B','C','D','E'}';
> undesirable = deaund(X0, Y0, Yu0, 'names', names);
> deadisp(undesirable);
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 5
Inputs: 1 Outputs: 1 Undesirable: 1
Model: directional-undesirable
Orientation: ddf (Directional distance function)
Returns to scale: crs (Constant)
--------------------------------------------------------------------
34 A Data Envelopment Analysis Toolbox for MATLAB
DMU| X| Y| Yu| Beta| slackX| slackY| slackYu|
--------------------------------------------------------------------
A| 1.0000| 7.0000| 2.0000| 0.0000| 0.0000| 0.0000| 0.0000|
B| 1.0000| 5.0000| 5.0000| 0.4000| 0.0000| 0.0000| 1.0000|
C| 1.0000| 1.0000| 3.0000| 0.8261| 0.7391| 0.0000| 0.0000|
D| 1.0000| 3.0000| 3.0000| 0.5556| 0.3333| 0.0000| 0.0000|
E| 1.0000| 4.0000| 2.0000| 0.2727| 0.2727| 0.0000| 0.0000|
--------------------------------------------------------------------
Another approach is to change the field names of the returned deaout structure before dis-
playing or exporting the results.
9.2. Changing the reference set
Sometimes we need to evaluate some DMUs with respect to a reference set that differs from
the observed input and output data. By default, all functions solve all models using the data
corresponding to X, Y and Yu –the observed DMUs, which also represents the reference set.
However, we can change the data corresponding to the DMU to be evaluated with the Xeval,
Yeval and Yueval optional parameters.12
As an example, using data from section 3, if we want to evaluate the first DMU with respect
to a reference set including all the other DMUs but not itself:
> load 'deadataFLS'
> Xref = X(2:end, :);
> Yref = Y(2:end, :);
> Xeval = X(1, :);
> Yeval = Y(1, :);
> io1 = dea(Xref, Yref, 'orient', 'io', 'Xeval', Xeval, 'Yeval', Yeval);
> disp(io1.eff);
1.1303
9.3. Custom display
When calling the deadisp(out, dispstr) function after computing a DEA model, appro-
priate information depending on the estimated model will be displayed on the screen. This
setting can be changed to display the desired information by specifying in the deadisp func-
tion the string dispstr (display string) as a second parameter.
The default dispstr after using the dea function is names/X/Y/eff/slack.X/slack.Y. Each
of the fields to be displayed in the output table must be separated with a / including the
names corresponding to the field names of the deaout structure. The available fields are
presented in Table 1.
As an example, we can display DMU names, the efficiency measure and the efficient X’s and
Y’s of the input-oriented model of Section 3 using names/eff/Xeff/Yeff in dispstr:
12
This is the procedure used internally by the toolbox to compute super-efficiency models and the produc-
tivity indices.
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 35
Common fields
Field Data
names DMU names
X Inputs
Y Outputs
eff Efficiency measure
slack.X Inputs slacks
slack.Y Outputs slacks
lambda Computed λ’s
Xeff Efficient X’s
Yeff Efficient Y’s
exitflag Exit flag of the optimization
Scale efficiency models
eff.crs CRS efficiency
eff.vrs VRS efficiency
eff.scale Scale efficiency
Malmquist index
eff.M Malmquist index
eff.MTEC Technical efficiency change
eff.MTC Technical change
Allocative efficiency model
Xprice Inputs’ prices
Yprice Outputs’ prices
eff.C Cost efficiency
eff.R Revenue efficiency
eff.P Profit efficiency
eff.A Allocative efficiency
eff.T Technical efficiency
Undesirable outputs model
Yu Undesirable outputs
slack.Yu Undesirable outputs’ slacks
Yueff Efficient Yu’s
Malmquist-Luenberger index
eff.ML Malmquist-Luenberger index
eff.MLTEC Technical efficiency change
eff.MLTC Technical change
Table 1: Fields of the deaout structure available for the dispstr string
36 A Data Envelopment Analysis Toolbox for MATLAB
> load 'deadataFLS'
> io = dea(X, Y, 'orient', 'io');
> deadisp(io, 'names/eff/Xeff/Yeff');
_______________________________
Data Envelopment Analysis (DEA)
DMUs: 11
Inputs: 2 Outputs: 1
Model: radial
Orientation: io (Input oriented)
Returns to scale: crs (Constant)
-----------------------------------------
DMU| Theta| Xeff1| Xeff2| Yeff|
-----------------------------------------
1| 1.0000| 5.0000| 13.0000| 12.0000|
2| 0.6223| 9.9566| 7.4675| 14.0000|
3| 0.8199| 13.1177| 21.3163| 25.0000|
4| 1.0000| 17.0000| 15.0000| 26.0000|
5| 0.3104| 5.5867| 4.3452| 8.0000|
6| 0.5556| 8.3333| 3.3333| 9.0000|
7| 1.0000| 25.0000| 10.0000| 27.0000|
8| 0.7577| 20.4571| 16.6687| 30.0000|
9| 0.8201| 28.7037| 11.4815| 31.0000|
10| 0.5000| 21.0000| 12.5000| 26.5000|
11| 1.0000| 5.0000| 13.0000| 12.0000|
-----------------------------------------
9.4. Exporting results
Results of the computed DEA models can be easily exported to diverse file formats for pos-
terior analysis and sharing. First the deaout structure should be converted to a MATLAB
table data type using the dea2table(out, dispstr) function using the desired dispstr.13 .
> load 'deadataFLS'
> io = dea(X, Y, 'orient', 'io');
> T = dea2table(io);
Then, the table can be exported using the MATLAB function writetable.14
> writetable(T, 'ioresults.csv');
13
If the dispstr parameter is omitted, the default one is used.
14
See the official MATLAB documentation for this function: http://www.mathworks.com/help/matlab/ref/
writetable.html
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 37
10. Conclusions
The new Data Envelopment Analysis Toolbox covers a wide variety of models calculating
efficiency and productivity measures in an organized environment for MATLAB. The models
implemented correspond to the classic radially oriented, the directional model and the additive
formulation. Both constant and variable returns to scale technical efficiency measures are
calculated, which allows the calculation of scale efficiency. The economic performance of
firms in terms of technical and allocative criteria is also presented, along with efficiency
models including undesirable outputs. Productivity indices are also implemented, both the
standard Malmquist index based on radial efficiency, and the Malmquist-Luenberger defined
in terms of the directional distance function. Finally, statistical analyses and hypotheses
testing using bootstrapping techniques are also available
We show how to organize the data, use the available functions and interpret results. To
illustrate the toolbox we solve several examples that are presented in reference dea texts
and handbooks. This positions the new toolbox as a valid self-contained package for these
evaluating techniques in MATLAB.
Since the code is freely available in an open source repository on GitHub, under the GNU
General Public License version 3, users will benefit from the collaboration and review of the
community, and can check the code to learn how dea optimizing programs are translated into
suitable code.15
Acknowledgments
We are grateful to Juan Aparicio for useful comments and suggestions. This research was
supported by the Spanish Ministry for Economy and Competitiveness under grant MTM2013-
43903-P.
References
Andersen P, Petersen NC (1993). “A Procedure for Ranking Efficient Units in Data Envelop-
ment Analysis.” Management Science, 39(10), 1261–1264.
Aparicio J, Pastor JT, Zofio JL (2013). “On the inconsistency of the Malmquist-Luenberger
index.” European Journal of Operational Research, 229(3), 738 – 742.
Aparicio J, Pastor JT, Zofio JL (2015). “How to properly decompose economic efficiency
using technical and allocative criteria with non-homothetic DEA technologies.” European
Journal of Operational Research, 240(3), 882 – 891.
Banker R, Charnes A, Cooper W (1984). “Some Models for Estimating Technical and Scale
Inefficiencies in Data Envelopment Analysis.” Management Science, 30(9), 1078–1092.
Barr RS (2004). Handbook on Data Envelopment Analysis, chapter Dea Software Tools and
Technology, pp. 539–566. Springer US, Boston, MA.
15
The address of the repository is https://github.com/javierbarbero/DEAMATLAB
38 A Data Envelopment Analysis Toolbox for MATLAB
Bogetoft P, Otto L (2011). Benchmarking with DEA, SFA, and R. Springer-Verlag New York,
New York: Springer.
Briec W (1997). “A Graph-Type Extension of Farrell Technical Efficiency Measure.” Journal
of Productivity Analysis, 8(1), 95–110. ISSN 1573-0441.
Caves DW, Christensen LR, Diewert WE (1982). “The Economic Theory of Index Numbers
and the Measurement of Input, Output, and Productivity.” Econometrica, 50(6), 1393–
1414.
Chambers RG, Chung Y, Färe R (1996). “Benefit and Distance Functions.” Journal of
Economic Theory, 70(2), 407 – 419.
Chambers RG, Chung Y, Färe R (1998). “Profit, Directional Distance Functions, and Nerlo-
vian Efficiency.” Journal of Optimization Theory and Applications, 98(2), 351–364.
Charnes A, Cooper W, Rhodes E (1978). “Measuring the efficiency of decision making units.”
European Journal of Operational Research, 2(6), 429 – 444.
Chung Y, Färe R, Grosskopf S (1997). “Productivity and Undesirable Outputs: A Directional
Distance Function Approach.” Journal of Environmental Management, 51(3), 229 – 240.
Cooper W, Pastor JT, Aparicio J, Borras F (2011). “Decomposing profit inefficiency in DEA
through the weighted additive model.” European Journal of Operational Research, 212(2),
411–416.
Cooper WW, Seiford LM, Tone K (2007). Data Envelopment Analysis: A Comprehensive Text
with Models, Applications, References and DEA-Solver Software. NJ: Princeton University
Press., Princeton.
Daraio C, Simar L (2016). “Efficiency and benchmarking with directional distances: a data-
driven approach.” Journal of the Operations Research Society. Forthcoming.
Du J, Liang L, Zhu J (2010). “A slacks-based measure of super-efficiency in data envelopment
analysis: A comment.” European Journal of Operational Research, 204(3), 694 – 697.
Econometric Software I (2009). LIMDEP. Plainview, NY. URL http://www.limdep.com/.
Emrouznejad A, Cabanda E (eds.) (2014). Managing Service Productivity: Using Frontier
Efficiency Methodologies and Multicriteria Decision Making for Improving Service Perfor-
mance. Springer-Verlag, Berlin.
Färe R, Grosskopf S, Norris M, Zhang Z (1994). “Productivity Growth, Technical Progress,
and Efficiency Change in Industrialized Countries.” American Economic Review, 84(1),
66–83.
Farrell M (1957). “The Measurement of Productive Efficiency.” Journal of the Royal Statistical
Society, Series A, General, 120(3), 253–81.
Fried HO, Lovell CAK, Schmidt SD (1993). The Measurement of Productivity Efficiency.
Techniques and Applications. Oxford University Press, New York.
Inmaculada C. Álvarez, Javier Barbero, José L. Zofı́o 39
Hollingsworth B (2004). “Non Parametric Efficiency Measurement.” The Economic Journal,
114(496), F307–F311.
Ji Y, Lee C (2010). “Data envelopment analysis.” Stata Journal, 10(2), 267–280(14).
Lovell CK, Pastor JT (1995). “Units invariant and translation invariant {DEA} models.”
Operations Research Letters, 18(3), 147 – 151.
Sheel H (2000). EMS: Efficiency Measurement System User’s Manual. University of Dort-
mund, Dortmund.
Sherman H, Zhu J (2006). Service Productivity Management: Improving Service Productivity
using Data Envelopment Analysis (DEA). Springer, New York.
Simar L, Wilson PW (1998). “Sensitivity Analysis of Efficiency Scores: How to Bootstrap in
Nonparametric Frontier Models.” Management Science, 44(1), 49–61.
Simar L, Wilson PW (1999). “Estimating and bootstrapping Malmquist indices.” European
Journal of Operational Research, 115(3), 459–471.
Simar L, Wilson PW (2002). “Non-parametric tests of returns to scale.” European Journal of
Operational Research, 139(1), 115 – 132.
StataCorp (2015). Stata Statistical Software: Release 14. StataCorp LP, College Station, TX.
URL http://www.stata.com.
The MathWorks, Inc (2015). MATLAB — The Language of Technical Computing, Ver-
sion R2015a (8.5). Natick, Massachusetts. URL http://www.mathworks.com/products/
matlab/.
Wilson PW (2008). “FEAR: A software package for frontier efficiency analysis with R.” Socio-
Economic Planning Sciences, 42(4), 247–254.
Zofı́o JL (2007). “Malmquist productivity index decompositions: a unifying framework.”
Applied Economics, 39(18), 2371–2387.
Zofı́o JL, Pastor JT, Aparicio J (2013). “The directional profit efficiency measure: on why
profit inefficiency is either technical or allocative.” Journal of Productivity Analysis, 40(3),
257–266.
You might also like
- Implementing DEA Models in The R Program: January 2013No ratings yetImplementing DEA Models in The R Program: January 20137 pages
- Stata DEA Command for Efficiency AnalysisNo ratings yetStata DEA Command for Efficiency Analysis14 pages
- A Review On The 40 Years of Existence of Data Envelopment Analysis Models: Historic Development and Current TrendsNo ratings yetA Review On The 40 Years of Existence of Data Envelopment Analysis Models: Historic Development and Current Trends30 pages
- A Review On The 40 Years of Existence of Data Envelopment Analysis Models - Historic Development and Current TrendsNo ratings yetA Review On The 40 Years of Existence of Data Envelopment Analysis Models - Historic Development and Current Trends30 pages
- Understanding Data Envelopment AnalysisNo ratings yetUnderstanding Data Envelopment Analysis14 pages
- Data Envelopment Analysis For Performance MeasurementNo ratings yetData Envelopment Analysis For Performance Measurement44 pages
- (International Series in Operations Research &Amp_ Management Science 239) Shiuh-Nan Hwang, Hsuan-Shih Lee, Joe Zhu (Eds.)-Handbook of Operations Analytics Using Data Envelopment Analysis-Springer US100% (1)(International Series in Operations Research &Amp_ Management Science 239) Shiuh-Nan Hwang, Hsuan-Shih Lee, Joe Zhu (Eds.)-Handbook of Operations Analytics Using Data Envelopment Analysis-Springer US511 pages
- An Introduction To Data Envelopment Analysis in Technology ManageNo ratings yetAn Introduction To Data Envelopment Analysis in Technology Manage7 pages
- Annals of Operations Research (In Press) : AcknowledgementsNo ratings yetAnnals of Operations Research (In Press) : Acknowledgements34 pages
- Emmanuel Thanassoulis Auth. Introduction To The Theory and Application of Data Envelopment Analysis A Foundation Text With Integrated Software33% (3)Emmanuel Thanassoulis Auth. Introduction To The Theory and Application of Data Envelopment Analysis A Foundation Text With Integrated Software296 pages
- The Final Frontier A SAS Approach To Data Envelopment AnalysisNo ratings yetThe Final Frontier A SAS Approach To Data Envelopment Analysis7 pages
- Model Cluster Limitations Using Standard Excel Solver: Two-Stage Network DEANo ratings yetModel Cluster Limitations Using Standard Excel Solver: Two-Stage Network DEA6 pages
- Data Envelopment Analysis Approach and Its Application in Information and Communication TechnologiesNo ratings yetData Envelopment Analysis Approach and Its Application in Information and Communication Technologies10 pages
- Data Envelopment Analysis (Dea) : A Tool For Measuring Efficiency and PerformanceNo ratings yetData Envelopment Analysis (Dea) : A Tool For Measuring Efficiency and Performance6 pages
- Data Envelo P Ment Analy Si S - Basi C Mo Dels and Thei R Uti Li Zati o NNo ratings yetData Envelo P Ment Analy Si S - Basi C Mo Dels and Thei R Uti Li Zati o N7 pages
- Journal of Interdisciplinary MathematicsNo ratings yetJournal of Interdisciplinary Mathematics12 pages
- A Guide To DEAP Version 2.1: A Data Envelopment Analysis (Computer) ProgramNo ratings yetA Guide To DEAP Version 2.1: A Data Envelopment Analysis (Computer) Program49 pages
- Jablonski (2009) - Software Support For Multiple Criteria Decision Making ProblemsNo ratings yetJablonski (2009) - Software Support For Multiple Criteria Decision Making Problems6 pages
- Data Envelopment Analysis DEA Thirty YeaNo ratings yetData Envelopment Analysis DEA Thirty Yea17 pages
- Data Envelopment Analysis: December 2017No ratings yetData Envelopment Analysis: December 20176 pages
- Inverse Data Envelopment Analysis For Parallel Production Systems: An Analysis in Four Types of Chinese BanksNo ratings yetInverse Data Envelopment Analysis For Parallel Production Systems: An Analysis in Four Types of Chinese Banks9 pages
- Banker Charnes Et Al 1989 CCR BCC DEA Model An Introduction To Data Envelopment Analysis With Some of Its Model and Their UseNo ratings yetBanker Charnes Et Al 1989 CCR BCC DEA Model An Introduction To Data Envelopment Analysis With Some of Its Model and Their Use20 pages
- Mercedes-Maybach GLS 600 4MATIC First Class Night Series MUXH22F6No ratings yetMercedes-Maybach GLS 600 4MATIC First Class Night Series MUXH22F65 pages
- Cambridge First Certificate in English3 For Updated Exam Upper Intermediate Students Book With Answers Frontmatter PDF100% (1)Cambridge First Certificate in English3 For Updated Exam Upper Intermediate Students Book With Answers Frontmatter PDF4 pages
- Test Your Understanding - Constructor (Copy) - Attempt Review100% (1)Test Your Understanding - Constructor (Copy) - Attempt Review7 pages
- Is - 4503 Specification For Shell and Tube Type Heat Exchangers100% (3)Is - 4503 Specification For Shell and Tube Type Heat Exchangers69 pages
- Chapter 6 - Big Data Architecture Part 1No ratings yetChapter 6 - Big Data Architecture Part 141 pages
- Safe Chemical Development Practices: Risks From Rising TemperatureNo ratings yetSafe Chemical Development Practices: Risks From Rising Temperature14 pages
- HBL632RT2: Construction Electrical Optics Specification FeaturesNo ratings yetHBL632RT2: Construction Electrical Optics Specification Features2 pages
- Mirpuri Vs Court of Appeals, 318 SCRA 516, G.R. No. 114508, November 19, 1999No ratings yetMirpuri Vs Court of Appeals, 318 SCRA 516, G.R. No. 114508, November 19, 199951 pages
