A Data Envelopment Analysis Toolbox for MATLAB

The Data Envelopment Analysis Toolbox is a new package for MATLAB that includes functions to calculate the main data envelopment analysis models. The package includes code for the standard radial input, output and additive measures, allowing for constant and variable returns to scale, as well as recent developments related to the directional distance function, and including both desirable and undesirable outputs when measuring eﬃciency 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 using a reliable productivity database on US agriculture to illustrate their use.


Introduction
Data envelopment analysis, DEA, has grown in importance over the past decades due to increase in the availability of data related to the performance of decision making units, DMUs, regardless their market, governmental, or non-for-profit orientation. DEA is a nonparametric method in operations research and economics that allows to practice benchmarking exercises among a set of comparable observations. This comparison is made through optimizing mathematical programming techniques that approximate the production technology and identify the most efficient facets or hyperplanes that, enveloping the data (hence the name), define the production or best-practice frontier, serving as reference for all observations. The meth-ods yield individual efficiency scores and reference benchmarks for each observation, thereby providing real-life peers that serve as example for managers and decision-makers interested in improving their performance. Because of the flexibility of DEA, researchers in a number of fields have quickly acknowledged these techniques as an excellent methodology for modeling operational processes. The empirical orientation and absence of a priori assumptions have resulted in numerous studies involving best-practice identification in a wide array of sectors -for a comprehensive introduction and a complete list of applications see Fried, Lovell, and Schmidt (2008). Indeed, the popularity of DEA has increased exponentially in these years as recent volumes devoted to specific sectors and industries bear witness; e.g., Blackburn, Brennan, and Ruggiero (2014), Ozcan (2014), Paradi, Sherman, and Tam (2018) and Khezrimotlagh and Chen (2018). The present toolbox allows interested practitioners to implement these studies in their everyday efficiency improvement strategies.
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, Inc 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; R Core Team 2020); commercial software -including trials versions, Emrouznejad and Cabanda (2014); free-ware programs - Sheel (2000); and even tutorials for spreadsheets: Sherman and Zhu (2006) and Zhu (2014). Earlier versions of these programs have been reviewed, among others, by Hollingsworth (2004) and Barr (2004). More recently, Daraio, Kerstens, Nepomuceno, and Sickles (2019) provide the most recent survey of the available options. Their eligibility criteria include only software that is diffused as self-contained programs, or packages and toolboxes for general computing environments. To complement the above references with general purpose DEA software that implements newer proposals, but excluding very specialized contributions that code specific models normally linked to particular papers, it is worth referencing Oh and Suh (2013) and Badunenko, Mozharovskyi, and Kolomiytseva (2020).
While these software packages implement the main DEA models, there is a lack of a full set of functions for MATLAB (The MathWorks, Inc. 2017), including some recent theoretical contributions that are missed in the existing software. 1 Thus, besides implementing the most popular DEA models in the MATLAB environment, the toolbox solves new models such as those corresponding to the directional and the weighted additive measures, while allowing for a greater degree of flexibility in the formulations; e.g., users can supply their preferred directional and weights vectors. Also, the models implemented, such as the directional measure with undesirable outputs, or the Malmquist-Luenberger index, correspond to the most recent theoretical proposals that solve relevant inconsistencies. As the most recent contribution to statistical software -mostly in R, the toolbox is up to date with respect to the models it includes, including bootstrapping techniques and hypotheses testing that are almost inexistent in older alternatives.
1 In principle, it would be desirable that our toolbox were compatible with the current distribution of Octave (v5.2.0; Eaton, Bateman, Hauberg, and Wehbring 2020). However, this is not the case with the out-of-the-box distribution. The reason is that the linprog function for linear programming in MATLAB is more advanced and flexible (in terms of the optimization options) than the equivalent function provided in Octave. In addition, some of the functions used are not available in Octave (such as optimoptions) and need to be replaced with other functions that do not have the same functionality as in MATLAB. Hence, it is not possible to make the toolbox compatible with Octave while maintaining some of its desirable features. However, any reader interested in solving any of our DEA models in Octave may adapt the syntax of our functions so they run in this program.
Nevertheless, it is virtually impossible to keep track and implement all current proposals in the literature. We refer the interested reader to the most recent specialized publications in DEA topics. In particular Zhu (2015) includes a series of chapters on current DEA theory dealing with most advanced topics such as stochastic nonparametric frontier approaches, efficiency measurement with fuzzy data, models with production trade-offs and weight restrictions, etc. Daraio and Simar (2007) deal with statistically robust DEA methods, including bootstrapping techniques -partially implemented in R by Simm and Besstremyannaya (2020). Finally Aparicio, Lovell, and Pastor (2016) and Hwang, Lee, and Zhu (2016) include state of the art and recent advances in the field of efficiency and productivity theory. Indeed, issues related to enhancing the discriminatory power of DEA as well as its statistical robustness are still at the forefront of current research in DEA, e.g, Angulo-Meza and Pereira Estellita Lins (2002). However while these topics merit special attention, their implementation exceeds the scope of the present software, which should be seen as the baseline tool for further specific implementations based on the proposed functions.
Arguably, a separate toolbox should be devoted to topics that are relevant by themselves and self-contained, such as weight restrictions and the inclusion of a priori information, or statistical methods such as robust, stochastic, and fuzzy DEA methods; e.g., see the recent reviews on stochastic DEA and fuzzy DEA by Olesen and Petersen (2016) and Hatami-Marbini, Emrouznejad, and Tavana (2011), respectively, and related applications by Kao and Liu (2009) and Costantino, Dotoli, Epicoco, Falagario, and Sciancalepore (2012b).
The Data Envelopment Analysis Toolbox introduces a complete set of baseline functions, covering a wide range of efficiency and productivity models, and reporting numerical results based on classical examples presented in the literature. The 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. 2 The paper is organized as follows. The following section presents the data structures characterizing the production possibility sets, the structure of the functions, results, etc., and briefly describes the data that is used to illustrate the toolbox. 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. The complementary cross efficiency model that also allows to overcome the low discriminatory power of standard DEA is presented in Section 4. Malmquist productivity indices and their decomposition into efficiency change and technical change are shown in Section 5, while Section 6 deals with the measurement of economic efficiency, and its decomposition into technical and allocative factors. Section 7 is devoted to the measurement of efficiency with undesirable outputs, most notably environmental efficiency, followed by Section 8 presenting the Malmquist-Luenberger index. Statistical analyses and hypotheses testing using bootstrapping techniques are discussed in Section 9. Advanced options, including displaying and exporting results can be found in Section 10. Section 11 concludes.

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 positive inputs x ∈ R m ++ into a vector of i = 1, 2, . . . , s positive outputs y ∈ R s ++ using the technology represented by the following constant returns to scale production possibility set: 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: 3 • 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.
• names: Names of the DMUs.

Dataset and statistical sources
We illustrate all the models presented in this toolbox resorting to a single dataset. The data has been collected and tabulated by the Economic Research Section (ERS) of the United States Department of Agriculture (USDA) in an effort to study long term productivity trends using prices-based indices. It corresponds to a subset of the state-level tables including price indices and implicit quantities of farm outputs and inputs (Table  23), which can be downloaded in full from https://www.ers.usda.gov/data-products/ agricultural-productivity-in-the-us/. More information on the data, including a discussion on the agricultural productivity growth in the US can be found in the above link. Due to its comprehensiveness and reliability, this dataset has been used over the years by many authors, whose results are complemented with those obtained by solving the models included in this toolbox -see, e.g., Färe, Grosskopf, and Margaritis (2008), Ball, Färe, Grosskopf, and Nehring (2001) and Zofío and Lovell (2001).

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 Charnes,Cooper,and Rhodes,ccr,model: 4 The optimal solution to this program -characterizing a technology with constant returns to scale, crs, is denoted by θ * crs . The constraints require the observation (θ crs x o , y o ) to belong to P crs , while the objective seeks the minimum θ crs that reduces the input vector x o radially to θ crs x o while remaining in P crs ; i.e., that projects the observation to the production frontier. The frontier corresponds to the supporting hyperplane defined by the linear combination of the observations that serve as reference benchmark for the evaluated unit. Those observations whose associated λ multipliers are greater than zero define the enveloping hyperplane. A feasible solution signaling radial efficiency is θ * crs = 1. Therefore if θ * crs < 1, the observation is radially inefficient and (λX, λY ) outperforms (x o , y o ). With regard to this property, we define the additional input excesses and output shortfalls by the following slack vectors: s − ∈ R m and s + ∈ R s , respectively. Therefore: s − = θ * crs x o − Xλ, and s + = Y λ − y o 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 subject to As a result, an observation is technically efficient if the optimal solution (θ * crs , λ * , s − * , s + * ) of the two above programs satisfies θ * 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). 5 The measurement of technical efficiency assuming variable returns to scale, vrs, as introduced by Banker et al. (1984) -known as the Banker, Charnes and Cooper, bcc, model, considers the following production possibility set P vrs = {(x, y) |x Xλ, y Y λ, eλ = 1, λ 0.}. Therefore, the only difference with the crs model is the adjunction of the condition n j=1 λ j = 1. Calculating the vrs efficiency, along with the subsequent second stage program analogous to (2), yields the corresponding optimal solution (θ * vrs , λ * , s − * , s + * ). As before, an observation is efficient with respect to the vrs technology if the optimal solution of the two programs satisfies θ * vrs = 1, s − * = 0, and s + * = 0. Finally, the scale efficiency, SE, of each observation is calculated as the ratio of the vrs to the crs scores: SE = θ * crs /θ * vrs . As a result the radial technical efficiency of an observation can be decomposed into its variable returns to scale efficiency ("pure" technical efficiency, PTE) and scale efficiency: 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 assumption can be specified by setting the rts parameter to crs (constant returns to scale; default) or vrs (variable returns to scale). With the optional parameter names we can specify a cell string with the name of the decision making units, which in this case are the names of the American states. 6 Results are returned in a 'deaout' structure and can be accessed directly (see Section 2) or displayed using the deadisp function. The second parameter of the deadisp function allows to specify the information that is displayed on screen. If not specified, the information displayed will depend on the calculated model. 7 load DataAgriculture X = [CAPITAL04, LAND04, LABOR04, INTINP04]; Y = [LS04, CROPS04, FARMOUT04]; X = X ./ 1000; Y = Y ./ 1000; statenames = STATE_NAME04; io = dea(X, Y, orient , io , names , statenames); deadisp(io, names/eff/slack.X/slack.Y ));

_______________________________ Data Envelopment Analysis (DEA)
DMUs: 48 5 It is possible to solve both programs in a single stage formulation employing a "non-Archimedian" infinitesimal constant , e.g., Ali and Seiford (1993, p. 140). However, this may result in computational inaccuracies and erroneous results.
6 If the optional parameter names is not used, observations (DMUs) are numbered from 1 to n. 7 See Section 10.3 for advanced uses of the deadisp function. Besides the US agricultural data, a second set of data used by Ali and Seiford (1993) is available at the toolbox website. A previous version of this toolbox illustrated the input, output, directional and additive models with these data.

Inputs: 4
Outputs: 3 Model: radial Orientation: io (Input oriented) Returns to scale: crs (Constant)  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.
Results show the relative difference between the constant returns to scale (crs) and variable returns to scale (vrs) efficiency scores. Since the vrs model envelops the data more tightly as a result of the additional constraint eλ = 1, these scores are greater than their crs counterparts, and therefore scale efficiency is also equal or less than one. For the selected states we observe that while none of them is technically efficient in the constant returns to scale model, with scores less than one, Wisconsin (WI) is efficient under the vrs assumption with a unitary score, and a relative scale efficiency of SE = 0.8452. Therefore, this state presents a suboptimal scale in terms of relative input and output quantities when compared to other counterparts that are almost efficient under crs; e.g., Arkansas (AR). As for the input and output slacks, corresponding to no radial reductions and increased, we observe that there is an excessive usage of the second input (land), whose values are greater than zero across the majority of states. Other intermediate inputs and livestock output do not exhibit slacks at all.

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): In this case, the optimal solution is denoted by φ * crs with the constraints ensuring that (x o , φ * crs y o ) belongs to P crs . Now the objective seeks the maximum φ crs that increases the output vector y o radially to φ * crs y o while remaining in P crs . A feasible solution signaling radial efficiency is φ * crs = 1. Therefore if φ * crs > 1, the observation is radially inefficient and (λX, λY ) outperforms (x o , y o ). Again, there might be further input excesses and output shortfalls: s − = x o − Xλ, and s + = Y λ − φ * crs y o with s − ≥ 0 and s + ≥ 0 for any feasible solution (φ, λ). To calculate these slacks in a second stage, the corresponding program incorporating the optimal value φ * crs is needed: subject to Finally, it is also possible to calculate the technical efficiency with respect to P vrs by solving the programs for the radial component and its associated input and output slacks analogous to (3) and (4), but adding the vrs constraint n j=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 in this case as SE = φ * crs /φ * vrs and radial efficiency can be now decomposed into pure technical efficiency, PTE, and scale efficiency: 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 , names , statenames); deadisp(oo, names/eff/slack.X/slack.Y );

_______________________________
Regarding the output oriented results, we first highlight that from DEA theory the output constant returns to scale scores are the inverse of their input counterparts; i.e., φ * crs = 1/θ * crs , and therefore none of our selected states are efficient. Additionally, it is also the case that if an observation is efficient under variable returns to scale from an input orientation, it is efficient from the output perspective (and for the forthcoming directional and additive models). Consequently Wisconsin (WI) is also vrs efficient from an output perspective. However scale efficiency can be substantially different depending on the orientation. The consistency of the results between the input and output models is also inferred from the positive values of the slacks corresponding to the second input (land). Chambers, Chung, and Färe (1996) introduced a measure of efficiency that projects observa-

The directional model: Constant and variable returns to scale
In this occasion the optimal solution to this program corresponds to β * crs . Now β * crs = 0 signals directional efficiency. Therefore if β * crs > 0, the observation is inefficient and (λX, λY ) It is again possible that further input excesses and output shortfalls exist. The slacks correspond to y , respectively, 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 subject to As in the inputs and outputs oriented models, one may also calculate technical efficiency with respect to P vrs . This requires solving equivalent programs to (5) and (6) adding the vrs constraint n j=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, PTE, and the scale efficiency term: 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 inputs and outputs set the direction: , 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, However, other choices are available; particularly −g − x , g + y = (−1, 1) or the mean of the , 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. 8 When deciding on the directional model, the researcher must declare whether the direction corresponds to the observed input or output 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, names , statenames); deadisp(ddf, names/eff/slack.X/slack.Y );

_______________________________
The results corresponding to the directional distance function are not bounded by one as in the previous input and output models. Nevertheless, if any observation is efficient in the former case with a score equal to one, either under crs or vrs, it is also efficient under the directional approach because it defines the production frontier; i.e., inputs and outputs cannot be simultaneously reduced or increased, respectively. That is why Wisconsin (WI) presents a zero valued score under vrs. These scores can be compared among themselves, and in this case the greater the value, the more inefficient is the state. Correspondingly, as in the previous models, the most inefficient state among the selected ones is Wyoming (WY). The values of the input and output slacks also present a similar pattern to that previously commented.

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 − ∈ R m and s + ∈ R s , respectively. The toolbox implements the weighted additive formulation of Lovell and Pastor (1995) and Pastor, Lovell, and Aparicio (2011), whose associated linear program is: subject to 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): other fixed values with the weights representing value judgments. For observation (x o , y o ) the objective seeks the maximum feasible reduction in its inputs and increase in its outputs while remaining in P vrs . 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 larger the sum of the slacks, ω * vrs , the larger the inefficiency. The relevance of these transformations is that they make the additive measures independent of the units of measurement, which is a desirable property.
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.
For illustration purposes, the range adjusted measure (RAM) 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, ... names , statenames); deadisp(add_ram, names/eff/slack.X/slack.Y ); Again, if an observation is identified as efficient by belonging to the the production frontier under any of the previous models, including zero valued input or output slacks, it shows an inefficiency score equal to zero in the additive models. This is the case of Wisconsin (WI). Note also that since the slacks weighting scheme is introduced in the objective function, the value of the optimal input and output slacks may remain the same. This is the case if the evaluated state is projected to the strongly efficient frontier, i.e., once slacks are accounted for, as it is the case of Arkansas (AR) in the above results. This not the case however for the scores of inefficient states. Taking Arizona (AZ) as the exemplifying case, its score is ω = 0.1851 if the measure of inefficiency proportions (MIP) is calculated, while if the ranges of the variables are used to weight the slacks (RAM), the score reduces to ω = 0.0149. It is relevant to remark that the use of alternative weights does not only change the value of the inefficiency scores, but may also change the relative rankings, even if the weights are the same for all observations, i.e., as in the RAM case. Hence, researchers should bear in mind that the choice of weights is not neutral, because it may result in alternative efficiency rankings.

Super-efficiency models
One interesting model that allows to differentiate across technically efficient observations is that proposed by Andersen and Petersen (1993). As Angulo-Meza and Pereira Estellita Lins (2002) highlight, one of the relevant drawbacks of standard DEA models is their inability to rank weakly efficient units differently as they all are assigned the same unitary efficiency score. The super-efficiency scores allow discriminating across them by calculating an individual score that is different across observations. These scores are obtained by individually solving for each observation any of the previous models, but excluding them from the reference dataset, which therefore reduces to n − 1 observations; i.e.,P crs ( The magnitude of the super-efficiency score, as well as the number of observations whose efficiency changes as a result of each individual exclusion determine 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 previously 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 increased and outputs reduced to reach the reference benchmarks. 9 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-Archimedian formulation: subject toθ where s − ∈ R m and s + ∈ R s , 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 set -is again denoted byθ * crs . In this case, ifθ * crs > 1, the observation is super-efficient and the larger the score and the increase in the values of the remaining observations with respect to the original calculations in (1), the more relevant is the observation that has been removed when defining the production frontier. Also, if an observation is inefficient when solving (1). its efficiency score does not change when removing it from the production possibility set, simply because it never defined the production frontier in the first place, and its original benchmarks remain the same. Implementing the super-efficiency method in the directional distance function model follows similar steps, while variable returns to scale can be computed by adding the corresponding constraint: n−1 j=1, =o λ j = 1. 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).
The results obtained for the super-efficient radial input oriented model under vrs show that once the observations are removed from the production possibility set, their efficiency scores may exceed the unitary value. This is the case of Wisconsin (WI), withθ * vrs = 1.0273 > 1, showing that input quantities must be actually increased to reach the production frontier. On the contrary, for the rest of the selected states, their efficiency scores remain unchanged from those reported for (1) since they are inefficient. Similarly, for the super-efficient directional distance function model. The optimal inefficiency score for Wisconsin (WI) is now negative: β * vrs = −0.0104 < 0, and therefore inputs must be increased and outputs reduced to reach the production frontier. Again, the values for the remaining states do not change from those obtained when solving (5). Note also that the super-efficiency model yields its own set of optimal input and output slacks.
As for the super-efficiency calculations in the additive model, we follow Du, Liang, and Zhu (2010), and solve the corresponding counterpart to (7) -for the case of the standard additive model: subject to The constraints in the program are modified as inputs need to be increased and outputs reduced so as to reach the production possibility set. Again, (ρ − x , ρ + y ) ∈ R m + × R s + are the corresponding weight vectors, and changing their value allows to calculate a wide range of efficiency scores, including the MIP, RAM and BAM models. By default, as in (7), program (9) calculates the measure of inefficiency proportions (MIP) with (ρ − x , ρ + y ) = (1/x o , 1/y o ). However one may set the individual weights to any value, including the unitary values corresponding to the standard additive model. On 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 using the function deaadditsuper(X, Y, ...), and specifying returns to scale by setting the rts parameter to crs (constant returns to scale; default) or vrs (variable returns to scale). Inputs and outputs weights are specified in the rhoX and rhoY parameters. If omitted, the weights corresponding to the MIP model are used.  (9) solves the super-efficiency model only for the states that are efficient in its standard additive counterpart (7), as inefficient observations have no interpretable solution. In this case, the MATLAB output returns NaN for these observations. That is why previous inefficient states are not reported in the above results (e.g., Alabama (AL), Arkansas (AR), etc.). Only Wisconsin (WI) is shown with a super-efficiency score ofω * vrs = 0.0168 > 0. In this case, the super-efficiency slack can be found on the output side, implying that to reach the production frontier from which this state is excluded, only the first output (livestock) must be reduced. This is the reason why its super-efficiency score is rather low when compared to other efficient states such as California (CA), withω * vrs = 1.5210. California presents both a greater slack value in the first output, as well as positive values in the remaining two outputs (crops and other farm-related output). Other states, e.g., Delaware (DE) withω * vrs = 1.0445, present super-efficiency slacks both in the input and output dimensions, but their absolute values are less than those of California, and therefore this state ranks in a lower position. From these results it is worth remarking that a livestock slack is consistently found for the majority of states, showing its relevance when defining the production frontier in the standard additive model.

Cross efficiency
As the previous approach, cross efficiency models allow to overcome the low discriminatory power of standard DEA, incapable of ranking the subset of efficient firms. Ultimately, the flexibility of DEA when searching for the most favorable weights may turn against the method itself by hampering its discriminatory power. A problem that is aggravated when the degrees of freedom are limited, with a small number of observations relative to the number of inputs and outputs. Initially developed by Sexton, Silkman, and Hogan (1986), cross efficiency methods rely on peer appraisals instead of self-evaluations. Here we follow Doyle and Green (1994) and implement the cross efficiency model based on the multiplier formulations of the radial input measures under constant returns to scale, (1). The basic cross efficiency concept uses the obtained weights for observation j, (u * j i , v * j i ), and multiplies them by the input and output quantities of the unit under evaluation (x o , y o ), thereby generating a peer-appraisal (cross efficiency) score for every pair (j, o). This peer-appraisal efficiency score corresponds to: After performing these calculations, one is left with a n × n matrix of peer-appraisal scores, from which the cross efficiency scores for every observation can be derived. θ j,o if they are excluded. However, because the optimal weights (u * j i , v * j i ) obtained from the first stage DEA run are not unique, one may obtain different results for e o in successive implementations. To remedy this situation a secondary goal is defined. The general idea of the benevolent and aggressive approaches is to solve the cross efficiency model with an objective to maximize or minimize the sum of all peer-appraisal scores, subject to the restrictions that: a) the self-appraisal scores remain equal to the results of the standard first stage DEA run, and b) no peer-appraisal score is greater than one. Because optimizing over the sum of ratios results in a non-linear problem, Sexton et al. (1986) propose a "linear surrogate", which uses the sum of all observed j inputs and outputs, j Program (10) corresponds to the benevolent approach. Minimizing of the objective function results in the aggressive formulation. One may average the results obtained from both approaches, for a less extreme set of cross efficiencies. 10 The cross efficiency model can be computed in MATLAB using the deacross(X, Y, ...) function with the objective parameter set to benevolent for the benevolent formulation (maximization) or to aggressive for the aggressive formulation (minimization). If the parameter mean is set to inclusive the cross efficiency score is calculated including the evaluated DMU, whereas if it is set to exclusive the evaluated DMU is excluded. Results for the benevolent and aggressive models show that, when compared to their standard radial input measure (1) counterparts, cross efficiency scores have smaller values. This is consistent with the fact that the optimal solution to the standard DEA corresponds to the most favorable weights (u * o i , v * o i ) that the program can find, and therefore assessing efficiency based on any other pair of weights (u * j i , v * j i ) results in smaller individual cross efficiencies; i.e., θ o,o > θ j,o , and therefore their inclusive or exclusive averages. For illustrative purposes, and taking as example West Virginia (WV), its standard efficiency score (1) is θ * crs = 0.6612, greater than its benevolent cross efficient counterpart, CrossEff = 0.4597. In turn, this latter score is larger than that corresponding to the aggressive approach 0.4005. Similar remarks can be made for the other selected states. Note also that for efficient states, their cross efficiency score would be less than one, and a comprehensive ranking could be obtained using these approaches.

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 performance with respect to reference technologies corresponding to two different time periods. The constant returns to scale production possibility set corresponding to t = 1,. . . ,T periods defines as P t crs = (x, y) |x X t λ, y Y t λ, λ 0 , with the variable returns to scale counterpart denoted accordingly by P t vrs . The standard Malmquist index relies solely on the concept of radial efficiency and requires calculation of the input -or output -oriented scores of (x t o , y t o ) 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, P 1 crs , we denote both scores by θ 1,1 crs and θ 2,1 crs , where the first superscript refers to the time period of the observation and the second one to that of the reference technology. While θ 1,1 crs is the solution to program (1), the intertemporal score θ 2,1 crs corresponds to the following program that evaluates period 2 observation (x 2 o , y 2 o ) with respect to period 1 technology: min Equivalently, the analogous intertemporal score θ 1,2 crs using the second period technology as reference corresponds to the same program but reverses the time superscripts of the firm under evaluation from (x 2 o , y 2 o ) to (x 1 o , y 1 o ), 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 possible to define the following Malmquist indices: M 1 = θ 2,1 crs /θ 1,1 crs and M 2 = θ 2,2 crs /θ 1,2 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 technical change. 11 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 periods, which defines for both indices as MTEC = θ 2,2 crs /θ 1,1 crs . The latter corresponds to the frontier-shift effect, i.e., the change in the reference frontier between both periods, which defines for M 1 as MTC 1 = (θ 2,1 crs /θ 2,2 crs ) using period 2 observation as the reference benchmark to evaluate the shift in the frontier. For M 2 it defines a MTC 2 = (θ 1,1 crs /θ 1,2 crs ). As in the previous cases, if MTEC > 1 or MTC > 1 productivity change is driven respectively by both technical efficiency gains and technical change improvements (technical progress); while MTEC < 1 or MTC < 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 as M 1 = MTEC × MTC 1 = θ 2,2 crs /θ 1,1 crs × θ 2,1 crs /θ 2,2 crs -and similarly for M 2 , while the geometric mean of both indices and their corresponding decompositions represents a compromise between both values: M = MTEC × MTC = θ 2,2 crs /θ 1,1 crs × (θ 2,1 crs /θ 2,2 crs × θ 1,1 crs /θ 1,2 crs ) (1/2) . 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 cumulative 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 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 M 1 and M 2 . However, for simplicity, the indices can be computed taking the first period as the base for technical change, M 1 , by setting the parameter geomean to 0. Also, in case of having more than two time periods, the transitive Malmquist index taking always the first period as the reference technology can be computed by setting the fixbaset parameter to 1. Calculating the Malmquist productivity index for the selected states results in quite different trends and sources of productivity change in the five year period between 2000-2004. While Alabama (AL), Arkansas (AR), Wisconsin (WI) and West Virginia (VW) experience productivity growth, M > 1, Arizona (AR) presents productivity stagnation, M ≈ 1, and Wyoming (WY) exhibits productivity decline, M < 1. Moreover, looking at the sources of productivity change, technological progress drives productivity growth, MTC > 1, for all states but West Virginia (VW), MTC < 1. Efficiency change is in line with productivity change, contributing to its growth in the four states where progress is observed, MTEC > 1, and causing productivity decline in the remaining two states, MTEC < 1. Indeed, in the case of Wyoming (WY), falling behind the production frontier is the main cause of productivity decline: MTEC = 0.9115, while observed technological progress, MTC > 1.0061, cannot counterbalance this declining trend.

Allocative models: Economic efficiency
Assuming an economic optimizing behavior on the part of the observations, e.g., cost minimization, revenue maximization, or profit maximization, and the corresponding input and output positive prices: w ∈ R m ++ and p ∈ R s ++ , 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).

Cost efficiency
Since we are concerned with overall efficiency in the input space, we assume that observations 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 minimum cost of producing the output level y given the input price vector w: C (y, w) = min m i=1 w i x i |x X, λy o Y λ, λ 0 , which considers the input possibility set capable of producing y o . 12 For the observed output levels we can calculate minimum cost and the associated optimal quantities of inputs consistent with the production technology by solving the following program: subject to Once minimum cost is calculated, cost efficiency defines as the ratio of minimum cost to observed cost: CE = C (y, w) /wx o . 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 -withx o = θ * crs x and CE = TE × AE. Consequently, if CE < 1 and the observation is technically efficient, θ * crs = 1, all cost inefficiency is allocative, while if the projected benchmark uses the right proportions of the optimal input quantities, which we denote by x * , thenx o = x * , and the observation 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 input prices, as a matrix or as a row vector, must be included. 13 If input prices are the same for all DMUs the toolbox generates a matrix with the same values. If input prices differ between DMUs a matrix with individual price information can be supplied. In a previous version of the toolbox the economic models in this section were illustrated using an example from Cooper et al. (2007, p. 261).

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 input levels at the existing market prices, they are revenue efficient. Let us denote by R (x, p) the maximum feasible revenue using input levels x and given the output prices p: Xλ, y Y λ, λ 0 ; i.e., considering the output possibility set producible with x o .
In this case, we calculate maximum revenue along with the optimal output quantities by solving the following program: subject to Revenue efficiency defines as the ratio of observed revenue to maximum revenue: RE = py o /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), TE = 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) -witĥ y o = φ * crs y o , and therefore RE = TE ×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 projected benchmark produces the right mix of optimal output quantities given their prices, which we denote by y * , thenŷ o = y * , and the observation 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 output prices, as a matrix or as a row vector, must be included.
The cost and revenue efficiency models show that none of the selected states minimize the cost of production or maximize the revenue associated to agricultural activities. None of them present a score equal to one, and indeed they are far away from it except Arkansas (AR), with CE = 0.9187 and RE = 0.9504, and Arizona (AR), with CE = 0.8963 and RE = 0.9425. These two states present relatively high values of technical efficiency and allocative efficiency, but are inefficient in both dimensions, either from a cost or revenue perspective. The states that are cost efficient (not reported here) are Delaware (DE), Florida (FL), Georgia (GA), Iowa (IA), Idaho (ID) and Illinois (IL) -and similarly from a revenue efficiency perspective. Therefore all of them are both technical and allocative efficient. On the contrary West Virginia (VW) and Wyoming (WY) have rather low cost and revenue efficiency scores. In general, both sources of inefficiency are balanced, except in the case of Wyoming (WY), whose technical inefficiency is mainly responsible for all revenue inefficiency: TE = 0.5911 vs. AE = 0.9536.

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, they are profit efficient. The profit function defines as Π (w, p) = max Xλ, y Y λ, eλ = 1, λ 0 . Calculating maximum profit along with the optimal output and input quantities and requires solving: subject to 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 homogeneous of degree one in prices -invariant to proportional price changes; i.e., PE = (Π (w, p)−(py o −wx o ))/(pg + y + wg − x ) -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, TE = β * vrs , and the residual difference corresponding to the allocative profit efficiency term: AE = PE − β * vrs . Therefore, allocative efficiency defines as the difference between maximum profit and profit at the technically efficient projection of the observation: , and PE = TE + AE. Now, if PE > 0 and the observation is technically efficient, β * vrs = 0, all profit inefficiency is allocative, and the observation is unable to produce with the optimal combination of outputs and inputs. Contrarily, if the benchmark observation supplies and demands the profit maximizing quantities of optimal outputs and inputs, which we denote by y * and x * , then (ŷ o ,x o ) = (y * , x * ), and the observation is allocative efficient: AE = 0. The profit inefficiency 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, names , statenames); deadisp(profit, names/eff.T/eff.A/eff.P );

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 input and output 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., This normalization ensures once again that the profit efficiency is homogeneous 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, TE = ω * vrs , and the residual difference corresponding to the allocative profit efficiency term: AE = PE − ω * vrs . Consequently, allocative efficiency defines as the difference between maximum profit and profit at the technically efficient projection of the observation: , and once again PE = TE +AE. Now, if PE > 0 and the observation is technically efficient, ω * vrs = 0, all profit inefficiency is allocative, and the observation is unable to produce with the optimal 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.
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, is set 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 measure of inefficiency proportions (MIP) model if not included. addprofit = deaadditprofit(X, Y, rts , vrs , Xprice , W, Yprice , P, ... names , statenames); deadisp(addprofit, names/eff.T/eff.A/eff.P ); Looking at the profit efficiency results, either directional or additive, the state that maximizes profit by not incurring in technical and allocative inefficiencies is California (CA). Consequently this state defines the benchmark against which the remaining states compare themselves. In the profit efficiency models the larger the value of the score, the greater the inefficiency. In the results for the directional distance function model West Virginia (VW) and Wyoming (WY) rank again amongst the worse states, with efficiency scores equal to PE = 6.4673 and PE = 6.3739, respectively. Arkansas (AR) and Wisconsin (WI) present values that are relatively low, PE = 0.5104 and PE = 0.3881. As this latter state is efficient technically, TE = 0, all inefficiency is allocative, PE = AE = 0.3881, which is due to a wrong choice of input and output mixes; i.e., farms are not demanding and supplying the optimal amounts of inputs and outputs given their prices. The results are similar in qualitative terms if the weighted additive approach is used to measure profit efficiency, although the values are substantially different as a result of the change in the objective function. Nevertheless, we can observe that, once again, the change of model does not have a neutral effect on the efficiency rankings. For example, based on the results for the directional distance function, Alabama (AL) is ranked ahead of West Virginia (VW): PE = 1.3481 vs. PE = 6.4673, but this is reversed in the weighted additive model: PE = 1312.8669 vs. PE = 1222.2983. Therefore, considering an efficiency measure that takes slacks into account by projecting observations to the strongly efficient frontier, rather than to the weakly efficient frontier, where slacks may still exist, may have relevant implications.

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 environmental 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 assessing environmental performance and making recommendations to improve technical efficiency, one seeks increments in desirable production but reductions in undesirable outputs. To incorporate 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 a redefinition of the production technology where the initial vector of i = 1, 2, . . . , s outputs y ∈ R s ++ is partitioned into desirable and undesirable production, i.e., y = (y d , y u ), with y d ∈ R q ++ and y u ∈ R r ++ , respectively. This translates into the corresponding reference technology P crs = x, y d , y u |x Xλ, y d Y λ, y u = Y λ, λ 0 , characterizing undesirable outputs as weakly disposable (Cooper et al. 2007, Chapter 12). As the directional efficiency 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  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 x o , y d o , y u o along the pre-assigned direction corresponding to the output vector g y = y d , y u = 0 m+s , corresponds to the solution of the following program: Again, the optimal solution corresponds to β * crs , and if β * crs = 0, with λ o = 1, λ j = 0 (j = o), the observation is directional efficient. Otherwise, β * crs > 0 signals inefficiency and (λX, λY d , λY u ) outperforms x o , y d o , y u o . It is possible also to calculate any additional slacks associated to further increases of individual outputs as well as reductions in individual inputs and undesirable outputs. The undesirable 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. The environmental or eco-efficiency model based on the directional distance function takes into consideration the possibility of reducing undesirable outputs when assessing productive performance; in this case the possibility of reducing the number of pesticide exposures among the population as a results of the agricultural activity. It is therefore relevant to compare the results of this model to those obtained for the standard directional model (5). For example, focusing on Arkansas (AR), which is almost efficient in the variable returns to scale model ignoring undesirable outputs β * vrs = 0.0055, it becomes fully efficient when they are taken into consideration (even under the chosen constant returns to scale model, which implies that it is also efficient if variable returns to scale were considered). Moreover, this state defines the weakly efficient production frontier, since no slacks are observed in that dimension. This is not the case however for the rest of the selected states, which are environmentally inefficient. Also, it is worth remarking that Arizona (AZ), besides reducing inputs and increasing outputs in the magnitude represented by the efficiency score times its observed quantities, could further reduce the number of pesticide exposures in 1, 042 cases, as measured by its associated slack, showing that additional gains in efficiency could obtained in that individual dimension.

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 (ML) productivity 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 ML uses the directional efficiency score only -disregarding second stage slacks -and requires calculation of the mix period efficiency of x t o , y t,d o , y t,u o observed in two periods t = 1, 2. We denote both scores by β 1,1 crs and β 2,1 crs , where once again the first superscript refers to the time period of the observation and the second one to that of the reference technology. While β 1,1 crs is the solution to program (15), the intertemporal score β 2,1 crs corresponds to the following program that evaluates period 2 observation x 2 o , y 2,d o , y 2,u o with respect to period 1 technology: subject to where max {y t,u i } is the maximum observed amount of each undesirable output in all time periods . Correspondingly, the score β 1,2 crs , using the second period technology as reference can be calculated solving an equivalent program that evaluates (x 1 o , y 1,d o , y 1,u o ) with respect to the reference technology (X 2 , Y 2,d , Y 2,u ).
Once these scores are calculated the ML indices define as follows: ML 1 = (1 + β 1,1 crs )/(1 + β 2,1 crs ) and ML 2 = (1 + β 1,2 crs )/(1 + β 2,2 crs ). If ML > 1 efficiency increases and the evaluated unit is capable of producing more desirable output with less undesirable production, while if ML = 1 productivity remains unchanged, and ML < 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 than their previous Malmquist counterparts. The change in technical efficiency defines now as MLTEC = (1 + β 1,1 crs )/(1 + β 2,2 crs ), while the frontier-shift effects corresponding to technical change are MLTC 1 = (1 + β 2,2 crs )/(1 + β 2,1 crs ) and MLTC 2 = (1 + β 1,2 crs )/(1 + β 1,1 crs ). Again, if MLTEC > 1 or MLTC > 1 productivity change responds to both technical efficiency gains and technical change improvements (technical progress); while MLTEC < 1 or MLTC < 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 ML 1 = MLTEC × MLTC 1 = (1 + β 1,1 crs )/(1 + β 2,2 crs ) × (1 + β 2,2 crs )/(1 + β 2,1 crs ) -and similarly for ML 2 . Given a sequence of years, the toolbox calculates the year-to-year variation of ML 1 and the geometric mean of both indices: ML = MLTEC × MLTC = (1 + β 1,1 crs )/(1 + β 2,2 crs ) × ((1 + β 2,2 crs )/(1 + β 2,1 crs ) × ((1 + β 1,2 crs )/(1 + β 1,1 crs )) (1/2) . The Malmquist-Luenberger indices are computed in MATLAB calling the deamalmluen(X, Y, Yu, ...) function. 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 computed 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-Luenberger indices are computed as the geometric mean of M L 1 and M L 2 . However, for simplicity, the indices can be computed taking the first period as the base for technical change, M L 1 , by setting the parameter geomean to 0. Just as the undesirable output dimension can be incorporated to the directional efficiency cross-sectional model, the same can be done from a panel data perspective, redefining the notion of productivity to account for undesirable outputs. For this reason it is relevant to compare the productivity trends previously discussed in Section 5 for the Malmquist index (M ) with those reported here for the Malmquist-Luenberger indices (ML). In the latter case, the inclusion of the number of pesticide exposures as undesirable output may change the nature of productivity change. A good example is Wisconsin (WI), which shows productivity growth according to the Malmquist index, M = 1.0469, but exhibits productivity decline in the case of the Malmquist-Luenberger index ML = 0.9172. Moreover, compared with the differences observed for other states, this gap between both indices is remarkable, showing that if one is willing to adopt a definition that accounts for the externalities of production processes, the results can vary dramatically. Also, in the Malmquist case, while Arizona (AR) leads the ranking of productivity growth by far, with Wisconsin (WI) in second place, this is not the case for the ML index, whose value ML = 1.0389 is closely followed by that of Alabama (AL), ML = 1.0384, with Wisconsin (WI) falling from the second to the fifth place. As for the sources of productivity change, while technological progress, MTC > 1, was the main driver of productivity change in the Malmquist case for all but one of the selected states (West Virginia, WV), in the Malmquist-Luenberger results this is only observed for the first three states, MLTC > 1. Similarly, the values for the technical efficiency change component are different depending on the definition of productivity. These results simply show that productivity change is much lower when undesirable outputs are taken into account, and therefore its components also show smaller values. These trends suggest that the effort exerted in reducing pesticide run-offs, so as to internalize the damage that they produce on humans, is lower than that aimed at increasing marketable outputs and reducing costly inputs. Hence the difference in productivity growth in favor of the Malmquist values.

Bootstrapping DEA estimators
The efficiency scores and productivity indices calculated in the previous sections can be considered as estimates as their values are subject to uncertainty due to sampling variation. Indeed, because of the uncertainty in many real situations, input and/or output data could be in stochastic form (e.g., historical data) or arising from qualitative judgments (e.g., decision maker's preferences). 14 Simar and Wilson (1998) introduce bootstrap methods that, based on resampling, provide estimation 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 perturbing them with a random noise generated from a kernel density function with scale given with bandwidth 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; 14 For recent reviews on alternative methods to overcome the deterministic nature of DEA we refer the reader to Olesen and Petersen (2016) and Hatami-Marbini et al. (2011), who review the stochastic and fuzzy approaches to DEA, respectively. Stochastic DEA complements the usual technological axioms characterizing the deterministic production technology with a set of distributional assumptions that provide a statistical foundation for the model. The former authors discuss three recent extensions of deterministic DEA: (i) deviations from the standard frontier are modeled as stochastic variables, (ii) random noise in terms of measurement errors, sample noise, and specification errors is made an integral part of the model, and (iii) the frontier is stochastic as is the underlying production possibility set (PPS). These approaches allow for an estimation of stochastic inefficiency and for statistical inference, while maintaining an axiomatic foundation. The latter authors also review recent advances in fuzzy DEA by presenting a taxonomy of methods. In particular they propose a classification scheme based on the following typology: the tolerance approach, the α-level based approach, the fuzzy ranking approach, and the possibility approach. The interested reader is also referred to empirical implementations of these methods by Kao and Liu (2009) and Costantino et al. (2012b). 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, 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. 15 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 A similar procedure is implemented to determine the significance of the Malmquist indices -see Simar and Wilson (1999) for the specific algorithm.
16 Calculation of the bootstrapped radially input model with the default options takes around 26 seconds on a PC with an i7-8700K CPU (3.70GHz, 3696 Mhz), 32 GB DDR4 (2666 MH) of RAM, and making use of the parallel computing toolbox connected to 6 workers. Computing time rises to 37 seconds when running the returns to scale test -since both the constant and variables returns to scale need to be calculated, and up to 70 seconds for the bootstrapped Malmquist (including the calculation of several cross-period efficiency scores). The rest of the models in the toolbox are solved within one or two seconds.
The returns to scale test is computed 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. 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 M 1 and M 2 . However, for simplicity, the indices can be computed taking the first period as the base for technical change, M 1 , by setting the parameter geomean to 0. Bootstrapping allows to test whether the standard efficiency scores, the existence of constant returns to scale, and productivity change are statistically significant or not. Focusing first on the bootstrapped scores calculated for the radially oriented input measure in Section 3.1, one sees that at the 5% default significance level, all original scores exceed the upper confidence interval threshold, and therefore none of them are reliable estimates of efficiency change. Taking Alabama (AL) as example, its standard score from program (1), θ * crs = 0.9077, is above the confidence interval (effCI2 = 0.9042), while the bootstrapped index, effboot = 0.8471, falls within the confidence interval defined by the lower and upper bounds. Consequently, bias corrected scores provide a robust and statistically sound alternative. Note however that the ranking among the selected states does not change for this particular database and default settings, with Arkansas (AR) and Wyoming (WY) ranking in the first and sixth place, respectively. Regarding returns to scale, for the default settings the null hypotheses of constant returns to scale cannot be rejected, suggesting that these property should be considered when solving the different models. Turning now to the standard and bootstrapped Malmquist indices, all of them fall within the confidence interval defined by the lower and upper levels, which also comprise the unitary value, except in the case of Wyoming (WY), whose bounds are McLow = 0.8866 and McUpp = 0.9903. Therefore, in this particular case, it can be concluded that the bootstrapped value equal to Mboot = 0.9499, signaling productivity decline, is indeed different from one, which allows rejecting the hypothesis that productivity remains unchanged.

Advanced optimization options
As recommended by The MathWorks, Inc. (2017), by default all DEA programs are solved using the dual-simplex algorithm, with a constraint tolerance of 1e-7 and an optimality tolerance of 1e-10. However, if needed, optimization options can be changed by passing an 'optimoptions' structure to the optimopts optional parameter. 17 As an example, if we want to solve the input oriented DEA program using the interior-point algorithm: opts = optimoptions( linprog , display , off , Algorithm , ... interior-point ); io_ip = dea(X, Y, orient , io , names , statenames, secondstep , 0, ... optimopts , opts);

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 represent the reference set. However, we can change the data corresponding to the DMU to be evaluated with the Xeval, Yeval and Yueval optional parameters. 18 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: Xref = X(2:end, :); Yref = Y(2:end, :); Xeval = X(1, :); Yeval = Y(1, :); 17 See the official MATLAB documentation to create an optimization options structure at https://www. mathworks.com/help/optim/ug/optimoptions.html. 18 This is the procedure used internally by the toolbox to compute super-efficiency models and the productivity indices. io1 = dea (Xref, Yref, orient , io , Xeval , Xeval, Yeval , Yeval); disp(io1.eff); 0.9077

Custom display
When calling the deadisp(out, dispstr) function after computing a DEA model, appropriate 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 function 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.

Exporting results
Results of the computed DEA models can be easily exported to diverse file formats for posterior 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. 19 io = dea(X, Y, orient , io , names , statenames); T = dea2table(io); Then, the table can be exported using the MATLAB function writetable. 20 writetable(T, ioresults.csv );

Conclusions
The 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. Models that overcome the low discriminatory power of the standard DEA models such as the super-efficiency or cross efficiency model are also included. 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 also show how to organize the data, use the available functions and interpret results. To illustrate all the models in the toolbox we use a reliable database on U.S. agriculture. 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. 21