§ INTRODUCTION

Early contributions to nonparametric regression by Hildreth, 1954, Brunk, 1955, and Grenander, 1956 built exclusively on the convexity and monotonicity constraints of the regression function. However, extending these approaches from the univariate setting to the more general multivariate regression proved a vexing challenge. Since the development of an explicit piece-wise linear characterization of the multivariate convex nonparametric least squares (CNLS) by Kuosmanen, 2008, convex regression has attracted growing interest in econometrics, statistics, operations research, machine learning and related fields (e.g., [Magnani and Boyd, 2009], [Magnani and Boyd, 2009]; [Seijo and Sen, 2011], [Seijo and Sen, 2011]; [Lim and Glynn, 2012], [Lim and Glynn, 2012]; [Hannah and Dunson, 2013], [Hannah and Dunson, 2013]; [Mazumder et al., 2019], [Mazumder et al., 2019]; [Bertsimas and Mundru, 2021], [Bertsimas and Mundru, 2021]). The recent study by Yagi et al., 2020 applies insights from convex regression to impose shape constraints on a local polynomial kernel estimator.

Convexity and monotonicity constraints are particularly relevant in the microeconomic applications where the duality theory of production and consumption directly implies certain monotonicity and convexity/concavity properties for many functions of interest (e.g., [Afriat, 1967], [Afriat, 1967], [Afriat, 1972]; [Varian, 1982], [Varian, 1982], [Varian, 1984]). For example, the cost function of the firm must be monotonic increasing and convex with respect to the input prices. Similar to the fact that a density function must be non-negative and its definite integral is equal to one, the cost function must satisfy the monotonicity and convexity properties implied by the theory, otherwise it is not really a cost function at all. The recent developments in the convex regression enable researchers to impose the concavity or convexity constraints implied by the theory to estimate the functions of interest without any parametric functional form assumptions.

In recent years, convex regression and related techniques are increasingly utilized in the estimation of frontier cost and production functions in the field of productivity and efficiency analysis, a multidisciplinary field that is widely applied in such areas as agriculture, banking, education, environment, health care, energy, manufacturing, transportation, and utilities (e.g., [Kuosmanen et al., 2015], [Kuosmanen et al., 2015]; [Johnson and Kuosmanen, 2015], [Johnson and Kuosmanen, 2015]). Traditionally, this field was divided between two competing paradigms: Data Envelopment Analysis (DEA) ([Charnes et al., 1978], [Charnes et al., 1978]) and Stochastic Frontier Analysis (SFA) ([Aigner et al., 1977], [Aigner et al., 1977]; [Meeusen et al., 1977], [Meeusen et al., 1977]). DEA is a deterministic, fully nonparametric approach whereas SFA is a probabilistic, fully parametric approach. To bridge the gap between these two paradigms, Stochastic nonparametric envelopment of data (StoNED) ([Kuosmanen, 2006], [Kuosmanen, 2006]; [Kuosmanen and Kortelainen, 2012], [Kuosmanen and Kortelainen, 2012]) was proposed as a unified framework that combines virtues of DEA and SFA, encompassing both approaches as its restricted special cases.

In practice, convex regression and StoNED are computationally demanding approaches, requiring a user to solve a mathematical programming problem subject to a large number of linear constraints. For example, the additive CNLS formulation by Kuosmanen, 2008 is a quadratic programming (QP) problem, whereas the multiplicative logarithmic formulation first considered by Kuosmanen and Kortelainen, 2012 requires solving a nonlinear programming (NLP) problem. Therefore, most empirical applications published thus far make use of commercial QP and NLP solvers, which can be coded using high-level mathematical computing languages such as GAMS ([GAMS Development Corporation, 2013], [GAMS Development Corporation, 2013]) or MATLAB ([The Mathworks, Inc., 2017], [The Mathworks, Inc., 2017]). Johnson and Kuosmanen, 2015 present detailed examples of how to implement the most basic CNLS and StoNED formulations in MATLAB and GAMS. Recently, the Benchmarking package ([Bogetoft and Otto, 2010], [Bogetoft and Otto, 2010]) implemented in R ([R Core Team, 2020], [R Core Team, 2020]) includes a new function StoNED() to estimate the CNLS/StoNED model, however, the R implementation is currently restricted to the additive CNLS QP formulation model.

The lack of a comprehensive, powerful, reliable, and fully open access computational package for the CNLS, StoNED, and related methods has slowed down the diffusion of these techniques to the empirical practice, which still heavily relies on the simple DEA and SFA techniques that either assume away noise or rely on restrictive functional form assumptions.[A plethora of computational tools, packages and resources are available for the conventional DEA and SFA. In the open access environment (R), one could use Benchmarking ([Bogetoft and Otto, 2010], [Bogetoft and Otto, 2010]), FEAR ([Wilson, 2008], [Wilson, 2008]), rDEA ([Simm and Besstremyannaya, 2020], [Simm and Besstremyannaya, 2020]), nonparaeff ([Oh, DH and Suh, D], 2013], [Oh, DH and Suh, D, 2013]), npsf ([Badunenko et al., 2020], [Badunenko et al., 2020]), and frontiles ([Daouia and Laurent, 2013], [Daouia and Laurent, 2013]), among many others, for solving different types of DEA models. For the SFA models, there exist frontier([Coelli and Henningsen, 2020], [Coelli and Henningsen, 2020]), sfa ([Straub, 2014], [Straub, 2014]), and sfaR ([Dakpo et al., 2021], [Dakpo et al., 2021]). Other free packages such as pyDEA ([Raith et al., 2021], [Raith et al., 2021]) and DataEnvelopmentAnalysis ([Barbero and Zofío, 2021], [Barbero and Zofío, 2021]) are available to implement efficiency analysis. Further, there are also a plenty of DEA and SFA codes and packages for MATLAB, GAMS, and other high-level mathematical computing languages (e.g., [Kalvelagen, 2002], [Kalvelagen, 2002]; [Ji and Lee, 2010], [Ji and Lee, 2010]; [Belotti et al., 2013], [Belotti et al., 2013]; [Álvarez et al., 2020], [Álvarez et al., 2020]).

To lower the barrier for applied researchers and practitioners to apply more advanced techniques that help to relax unnecessarily restrictive assumption, the Python package pyStoNED was first introduced in April 2020 to prove a freely available and user-friendly tool for the multivariate CNLS and StoNED methods. Its latest edition 0.5.2 also includes modules for the convex quantile regression, convex expectile regression, isotonic regression, and graphical illustration. It also facilitates efficiency measurement using the conventional DEA and free disposable hull (FDH) approaches. The pyStoNED package allows practitioners to estimate these models in an open access environment under a GPL-3.0 License. The project, including source code, internal data, notebook tutorials, and web documentation, is publicly available at GitHub.[
The pyStoNED package GitHub repository: <https://github.com/ds2010/pyStoNED>.

The purpose of this paper is to present a tutorial of the pyStoNED package, briefly review the alternative models supported, and illustrate its application. We focus on the estimation of frontier cost and production functions, which currently forms the main application area of these techniques, emphasizing that the various modules of the pyStoNED package are directly applicable for (semi-)nonparametric regression analysis in any other contexts as well. We emphasize that the pyStoNED package is an ongoing development by the users for the users: further model specifications and methodological advances will be implemented and added to the pyStoNED package on a continuous basis.

The paper is organized as follows. Section <ref> describes the basic setups of the pyStoNED package, and Section <ref> introduces the structures of example data and the attributes of different models. Section <ref> describes the first step of the StoNED model (e.g., the CNLS estimator) to estimate the conditional mean and some other commonly used extensions. The Python code for the CNLS estimator and these extensions are included. Section <ref> demonstrates the rest of the steps of the StoNED model and related codes. Section <ref> illustrates how to implement the CNLS-G algorithm to calculate the CNLS estimator using the developed pyStoNED package. The plot of the estimated function can be found in Section <ref>. Section <ref> concludes this paper. The list of acronyms is presented in Appendix <ref>. Appendices <ref> and <ref> present the estimated residuals for a same CNLS model solved by GAMS and pyStoNED, respectively.

§ SETUP

§.§ Installation

The pyStoNED package supports Python 3.8 or later versions on Linux, macOS, and Windows, and is freely available on the Python Package Index (PyPI) at <https://pypi.org/project/pystoned>. The package can be installed by any of the following three approaches:

* Install from the PyPI

pip install pystoned

* Install from the GitHub repository

pip install -U git+https://github.com/ds2010/pyStoNED

* Clone developing branches and install with Python using the setup.py script (for developers)

git clone https://github.com/ds2010/pyStoNED.git

python setup.py install

The pyStoNED package is built based on a few existing dependencies.[
List of dependencies used in pyStoNED: NumPy ([Harris et al., 2020], [Harris et al., 2020]), pandas ([The pandas development team], 2020], [The pandas development team, 2020]), SciPy ([Virtanen et al., 2020], [Virtanen et al., 2020]), Matplotlib ([Hunter, 2007], [Hunter, 2007]), and Pyomo ([Bynum et al., 2021], [Bynum et al., 2021]). All the dependencies will be automatically added when installing the pyStoNED package.

It is worth highlighting that Pyomo is a full-featured high-level programming language that provides a rich set of supporting libraries to program  CNLS,  StoNED, and various extensions. The other dependencies are also essential for pyStoNED. Specifically, NumPy and pandas are used to import input-output data, manipulate data, and export the estimation results. SciPy provides the algorithm to optimize the quasi-likelihood function. Matplotlib is then associated to plot the estimated functions.

§.§ Solvers

In the CNLS/StoNED framework, the existing models are either additive or multiplicative models, depending on the specification of the error term's structure. In the context of the optimization problems, the additive models are usually the QP problem with an exception of the CQR model, which is a linear programming (LP) problem, and all multiplicative models are the NLP problem. The attribute of models determines which type of solver will meet our needs, i.e., QP/LP-solver or NLP-solver. In the pyStoNED package, we can import these two kinds of solvers locally or remotely.

§.§.§ Remote Optimization

With the help of Network-Enabled Optimization System (NEOS) Server that freely provides a larger number of academic solvers such as CPLEX ([Cplex, IBM ILOG, 2009], [Cplex, IBM ILOG, 2009]), MOSEK ([MOSEK ApS, 2021], [MOSEK ApS, 2021]), MINOS ([Murtagh and Saunders, 2003], [Murtagh and Saunders, 2003]), and KNITRO ([Byrd et al., 2006], [Byrd et al., 2006]),[
NEOS Server: State-of-the-Art Solvers for Numerical Optimization, <https://neos-server.org/neos>
pyStoNED has the ability to calculate the StoNED related models without installing the local solvers.  Here we have a model optimizing remotely:

>>> model.optimize(email="email@address", solver='mosek')

The use of the remote optimization requiring email and solver as parameters.[
Replace with your own email address required by NEOS server, see <https://neos-guide.org/content/FAQ#email>
If the user leaves the parameter solver out (i.e., model.optimize(email="email@address")), the optimization job will be calculated by the default solver. That is, the additive and multiplicative models will be solved by MOSEK and KNITRO, respectively. The remote optimization is highly recommended for light computing jobs.

§.§.§ Local Optimization

We can resort to the local solver pre-installed by ourselves to estimate the StoNED related models (e.g., MOSEK and CPLEX). The parameter OPT_LOCAL is added in the function .optimize(...) to indicate that the model optimizes locally. Here is an example of optimizing the additive model with local solver MOSEK.[
The tutorial of the MOSEK installation is available at <https://pystoned.readthedocs.io/en/latest/install/index.html>.

>>> model.optimize(OPT_LOCAL, solver='mosek')

Since the free and stable NLP solver remains scant, the pyStoNED package does not support the multiplicative model estimation locally at present. Hence, we recommend the user solve their multiplicative models remotely. We will update the package immediately when the free and stable NLP solver is available.

§ DATA STRUCTURES AND DATASET

§.§ Data structures

Data pre-processing is the first step to use the developed package, and thus the user must prepare the dataset based on the data structures of pyStoNED. In all CNLS/StoNED models, two common vectors or matrix are required: input variables $\bm{x}_i$ and output variables $\bm{y}_i$ for observed Decision Making Units (DMU) $i = 1, 2, \cdots, n$.  For a model considering operational conditions and practices, contextual variables $\bm{z}_i$ are required. The directional distance function (DDF) based models handle the multi-dimensional inputs and outputs with given directional vectors $\bm{g}^x$ and $\bm{g}^y$, respectively. The DDF based models also take undesired outputs $\bm{b}_i$ and its directional vector $\bm{g}^b$ into consideration. Table <ref> summarizes the data structures used in all the CNLS/StoNED models.

Symbol         Model                       Description
$\bm{x}$       All models                  Input variables

$\bm{y}$       All models                  Output variables

$\bm{z}$       Contextual based models     Contextual variables

$\bm{b}$       DDF based models            Undesirable outputs

$\bm{g}^x$     DDF based models            The direction of inputs

$\bm{g}^y$     DDF based models            The direction of outputs

$\bm{g}^b$     DDF based models            The direction of undesirable outputs

Data structures.

The inputs of the estimation functions support two forms: matrix from NumPy and list from Python. The outputs (i.e., final estimates) are retrieved by using the corresponding functions .get_xxx() provided by pyStoNED, all of which are in the form of ndarray from NumPy.

§.§ Internal data

To illustrate the application of the pyStoNED package, four commonly used datasets are attached:

* Regulation of Finnish electricity distribution firms (load_Finnish_electricity_firm)

The data of the regulation of Finnish electricity distribution firms are collected from Kuosmanen, 2012 and Kuosmanen et al., 2013. The data consist of seven variables: three different expenditures are used as inputs (i.e., OPEX, CAPEX, and TOTEX);[
Note that TOTEX = OPEX + CAPEX. It is possible to use TOTEX as an aggregate input, or model OPEX and CAPEX as two separate input variables.
Energy, Length, and Customers are considered as outputs; Further, PerUndGr is denoted as the contextual variable. Table <ref> presents the description of the dataset.

Variable             Unit                  Description
OPEX          Thousand Euro         Controllable operational expenditure

CAPEX         Thousand Euro         Total capital expenditure

TOTEX         Thousand Euro         Total expenditure

Energy        Gigawatt Hours        Weighted amount of energy transmitted

Length        Kilometer             Length of the network

Customers     Person                Customers connected to the network

PerUndGr      Percentage            Proportion of underground cabling

The regulation of Finnish electricity distribution firms.

* GHG abatement cost of OECD countries (load_GHG_abatement_cost)

The data of the Greenhouse gas (GHG) abatement cost of OECD countries are provided by Kuosmanen et al., 2020. The data contain two input variables (i.e., CPNK and HRSN), one good output variable (i.e., VALK), and one undesirable output variable (i.e., GHG). Table <ref> describes the dataset in detail.

Variable        Unit                                   Description
CPNK     Billion Euro$^{2010}$                  Net capital stock

HRSN     Billion hours                          Hours worked by total engaged

VALK     Billion Euro$^{2010}$                  Value added

GHG      Million tons of CO$_2$ equivalents     Total GHG emissions

The GHG abatement cost of OECD countries.

* Data provided with Tim Coelli's Frontier 4.1 (load_Tim_Coelli_frontier)

The classic 60 firms dataset attached in Frontier 4.1 ([Coelli, 1996], [Coelli, 1996]) includes two input variables (i.e., capital and labour) and one output variable (i.e., output) (see Table <ref>).

Variable           Unit               Description
firm        Quantity           Firm ID

output      Quantity Index     Output quantity

capital     Quantity Index     Capital input

labour      Quantity Index     Labour input

The data provided with Tim Coelli's Frontier 4.1.

* Rice Production in the Philippines (load_Philipines_rice_production)

The Rice Production in the Philippines dataset collected from Coelli et al., 2005 consists of 17 different variables. The different variables can be organized into diversified combinations for target models. Table <ref> summarizes the variables of the data.

Variable        	    Unit              Description
YEARDUM          Year              Time period

FMERCODE         Quantity          Farmer code

PROD     	    Tonnes            Toones of freshly threshed rice

AREA    		    Hactares          Area planted

LABOR    	    Mandays           Labour used

NPK      	    Kilogram          Fertiliser used

OTHER    	    Laspeyres index     Other inputs used

PRICE   		    Pesos/kilogram       Output price

AREAP   		    Pesos/hectare        Rental price of land

LABORP   	    Pesos/day            Labour price

NPKP             Pesos/kilogram       Fertiliser price

OTHERP           Implicit price index     Price of other inputs

AGE              Years                      Age of the household head

EDYRS            Years                      Education of the household head

HHSIZE           Quantity                   Household size

NADULT           Quantity                   Number of adults in the household

BANRAT           Percentage                 Percentage of area classified as upland fields

The rice production in the Philippines.

These datasets can be imported through the module dataset of pyStoNED. The variables in Tables <ref>–<ref> are used as parameters of dataset to load the input data. The following example demonstrates how to load the input-output data from the regulation of Finnish electricity distribution firms dataset. We first import the dataset module using the function load_Finnish_electricity_firm (i.e., Line 1),[
In addition to the presented example, we can load other datasets using the corresponding dataset module, e.g., loading the GHG data: from pystoned.dataset import load_GHG_abatement_cost.
then define the $\bm{x}$ and $y$ according to the imported dataset (Line 2). We can check the input and output data using the function print(...) (Lines 3 and 4). The parameters x_select, y_select, and z_select in load_Finnish_electricity_firm(x_select, y_select, z_select) are used to select the inputs, outputs, and  contextual variables, respectively.

>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'], z_select=['PerUndGr'])
>>> print(data.x)
>>> print(data.y)
>>> print(data.z)

[[    75    878   4933]
[   105    575   9084]]

[[  1612]
[  1776]]

Note that the parameters in the module dataset can be defined according to the user's purpose. For example, if the target model only consists of two inputs (e.g., Energy and Customers) and one output (e.g., TOTEX), then the data selection in Line 2 should be

>>> data=load_Finnish_electricity_firm(x_select=['Energy','Customers'],
...     y_select=['TOTEX'])

§.§ External data

In practice, the user's own dataset is the main input of pyStoNED. We present an example to show how to import the user's own data. Assume that Table <ref> is the input-output data stored in Excel file table1.xlsx, the following code utilizes the pandas to read the Excel file and organize the data with NumPy. The input variable x is a matrix and the output variable y is then an array.

id      output     input1     input2
1       120        10         55

2        80        30         49

.         .         .          .

.         .         .          .

100      90        25         72

An example of user's own dataset.

In the following case, we first import the packages NumPy and pandas in Lines 1 and 2. Line 3 is used to read the Excel data file and the rest of the lines define the input and output variables. See more similar examples in Section <ref>.

>>> import numpy as np
>>> import pandas as pd
>>> df=pd.read_excel("table1.xlsx")
>>> y=df['output']
>>> x1=df['input1']
>>> x1=np.asmatrix(x1).T
>>> x2= df['input2']
>>> x2=np.asmatrix(x2).T
>>> x=np.concatenate((x1, x2), axis=1)

§ SHAPE-CONSTRAINED NONPARAMETRIC REGRESSION

Consider a standard multivariate, cross-sectional model in production economics:

y_i     = f(x_i) + ε_i

= f(x_i) + v_i - u_i           ∀i
where $y_i$ is the output of the DMU $i$, $f: R_+^m \rightarrow R_+$ is the production (cost) function that characterizes the production (cost) technology, and $\bm{x}_i = (x_{i1}, x_{i2}, \cdots, x_{im})^{'}$ denotes the input vector of unit $i$. Similar to the literature in SFA, for production function, the presented composite error term $\varepsilon_i = v_i - u_i$ consists of the inefficiency term $u_i>0$ and stochastic noise term $v_i$. Note that by setting $u_i = 0$ we have the standard nonparametric regression model as a special case. To estimate the function $f$, one could resort to the parametric and nonparametric methods or neoclassical and frontier models, which are classified based on the specifications of $f$ and error term $\varepsilon$ (see [Kuosmanen and Johnson, 2010], [Kuosmanen and Johnson, 2010]). In this paper, we assume certain axiomatic properties (e.g., monotonicity, concavity) instead of a priori functional form for the function $f$ and apply the following nonparametric methods to estimate the function $f$.[
In this paper, we focus on introducing the CNLS/StoNED and other related models and illustrating how to apply the pyStoNED to solve these models. Of course, one can use this package to solve other nonparametric models such as DEA and FDH. See more tutorials about DEA and FDH estimations at <https://pystoned.readthedocs.io>.

§.§ Convex Nonparametric Least Square

§.§.§ Additive CNLS model

Hildreth, 1954 is the first to consider the nonparametric regression subject to monotonicity and concavity constraints in the case of a single input variable $x$. Afriat, 1972 also proposes methods to impose convexity on the estimation of a production function. Kuosmanen, 2008 extends Hildreth's approach to the multivariate setting with the multidimensional input $\bm{x}$, and refers it to as the CNLS. CNLS builds upon the assumption that the true but unknown production function $f$ belongs to the set of continuous, monotonic increasing and globally concave (convex) functions, imposing the same production axioms as standard DEA (see further discussion in [Kuosmanen and Johnson, 2010], [Kuosmanen and Johnson, 2010]). The additive multivariate CNLS formulations are defined as

* Estimating production function (i.e., regression function $f$ is concave and increasing)

α, β, ε min   ∑_i=1^nε_i^2

y_i = α_i + β_i^'x_i + ε_i          ∀i

α_i + β_i^'x_i ≤α_j + β_j^'x_i            ∀i, j ,  and   i ≠j

β_i ≥0          ∀i

* Estimating cost function (i.e., regression function $f$ is convex and increasing)

α, β, ε min    ∑_i=1^nε_i^2

y_i = α_i + β_i^'x_i + ε_i          ∀i

α_i + β_i^'x_i ≥α_j + β_j^'x_i            ∀i, j ,  and   i ≠j

β_i ≥0          ∀i

where $\alpha_i$ and $\bm{\beta}_i$ define the intercept and slope parameters of tangent hyperplanes that characterize the estimated piecewise linear frontier, respectively. $\varepsilon_i$ denotes the regression residual. The first constraint can be interpreted as a multivariate regression equation, the second constraint imposes convexity (concavity), and the third constraint imposes monotonicity. Similar to the DEA specification, other standard specifications of returns to scale can be imposed by an additional constraint on the intercept term $\alpha_i$. If $\alpha_i=0$, then Problems (<ref>) and (<ref>) are the constant returns to scale (CRS) model, otherwise they are the variable returns to scale (VRS) model. Note that both Problems (<ref>) and (<ref>) are the QP problem and hence can be solved by MOSEK and CPLEX.

The basic additive CNLS model can be estimated in pyStoNED using the module CNLS(y, x, ...) with the contextual variable z parameter set to None (default) and the type of model cet parameter set to CET_ADDI (additive model; default). The type of estimated function can be classified by setting the fun parameter to FUN_PROD (production function; default) or FUN_COST (cost function). The returns to scale assumption can be specified by setting the rts parameter to RTS_VRS (VRS model; default) or RTS_CRS (CRS model). The results can be displayed in the screen directly using the .display_alpha() (i.e., display the coefficients $\hat{\alpha}_i$) or stored in the memory using the .get_alpha().

We first present an example to solve the VRS production model and store the estimates as the following Lines 12-14.

>>> from pystoned import CNLS
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data = load_Finnish_electricity_firm(x_select=['OPEX','CAPEX'],
...     y_select=['Energy'])
>>> model=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> model.display_alpha()
>>> model.display_beta()
>>> model.display_residual()
>>> alpha = model.get_alpha()
>>> beta = model.get_beta()
>>> residuals = model.get_residual()

alpha : alpha
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -22.935939927355427 :  None : False : False :  Reals
1 :  None :  -22.98683384275405 :  None : False : False :  Reals
88 :  None : -22.860373724764905 :  None : False : False :  Reals

beta : beta
Size=178, Index=beta_index
Key     : Lower : Value                  : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 :    0.13585845731547977 :  None : False : False :  Reals
(0, 1) :   0.0 :   0.011273984231181894 :  None : False : False :  Reals
(1, 0) :   0.0 :    0.13637662344803758 :  None : False : False :  Reals
(1, 1) :   0.0 :   0.010903786808952633 :  None : False : False :  Reals
(88, 0) :   0.0 :    0.13561549262636433 :  None : False : False :  Reals
(88, 1) :   0.0 :   0.011417220278201383 :  None : False : False :  Reals

epsilon : residual
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -2.8024040090178914 :  None : False : False :  Reals
1 :  None :  1.4140528128759229 :  None : False : False :  Reals
88 :  None :  -0.885163858119256 :  None : False : False :  Reals

In this example, Lines 1-3 import the CNLS module, parameter setting modules, and dataset module. Line 4 defines the input and output variables using the Finnish electricity distribution firm data. Lines 5-6 define the CNLS production model with a user-defined name (e.g., model) and solve the production model using the local off-the-shelf solver (i.e., MOSEK).[
The name of estimated model is free to be assigned, e.g., CNLS_model=CNLS.CNLS(y=data.y, x=data.x, z=None, cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS). Correspondingly, the Line 6 in this example should be replaced by CNLS_model.optimize(OPT_LOCAL), and Line 8 now is CNLS_model.display_alpha().
Lines 8-10 directly display the estimated coefficients (i.e., $\hat{\alpha}_i$, $\hat{\beta}_{ij}$, and $\hat{\varepsilon}_i$) on screen, and Lines 12-14 store the estimates in memory with a special variable name (e.g, alpha). Appendix <ref> presents the full estimated CNLS residuals that are equivalent to those in GAMS (cf. Appendix <ref>). Note that the estimated alpha and beta in pyStoNED may be slightly different from those in GAMS due to non-uniqueness in CNLS estimation (see [Kuosmanen, 2008], [Kuosmanen, 2008]; [Dai, 2021], [Dai, 2021]).

§.§.§ Multiplicative CNLS model

Similar to the existing SFA literature, the Cobb-Douglas and translog functions are commonly assumed to be the functional form for the function $f$, where inefficiency $u$ and noise $v$ affect production in a multiplicative fashion. We, thus, further consider the multiplicative specification in the nonparametric models. Note that the assumption of CRS would also require a multiplicative error structure. Under the multiplicative composite error structure, the function (<ref>) is rephrased as

\begin{align}
y_i = f(\bm{x}_i) \cdot \exp{(\varepsilon_i)} = f(\bm{x}_i) \cdot \exp{( v_i - u_i)}  \label{eq:eq4}
\end{align}

Applying the log-transformation to Eq.(<ref>), we obtain

\begin{align}
\ln y_i = \ln f(\bm{x}_i) + v_i - u_i \label{eq:eq5}
\end{align}

To estimate Eq.(<ref>), we reformulate the additive production model (<ref>) and obtain the following log-transformed CNLS formulation:

α, β, ϕ, ε min    ∑_i=1^nε_i^2

lny_i = ln(ϕ_i+1) + ε_i           ∀i

ϕ_i  = α_i+β_i^'x_i -1          ∀i

α_i + β_i^'x_i ≤α_j + β_j^'x_i            ∀i, j ,  and   i ≠j

β_i ≥0           ∀i
where $\phi_i+1$ is the CNLS estimator of $\E[y_i \, | \, \bm{x}_i]$. The value of one is added here to make sure that the computational algorithms do not take the logarithm of zero. The first equality can be interpreted as the log-transformed regression equation (using the natural logarithm function $\ln(\cdot)$). The rest of the constraints are the same as those of the additive models. The use of $\phi_i$ allows the estimation of a multiplicative relationship between output and input while assuring convexity of the production possibility set in the original input-output space. Note that one could not apply the log transformation directly to the input data $\bm{x}$ due to the fact that the piece-wise log-linear frontier does not satisfy the axiomatic property (i.e., concavity or convexity) of the function $f$. Since the multiplicative model (<ref>) is a nonlinear programming (NLP) problem, we need to use the nonlinear solvers, e.g., MINOS, KNITRO.

We next demonstrate the estimation of both VRS and CRS multiplicative cost function models. Let the type of model cet parameter be CET_MULT (multiplicative model). Note that the following NLP models are remotely solved by KNITRO solver via the NEOS Server.

>>> from pystoned import CNLS
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RTS_CRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data = load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model1=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model1.optimize('email@address')
>>> model2=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_CRS)
>>> model2.optimize('email@address')
>>> model1.display_residual()
>>> model2.display_residual()

epsilon : residual
Size=89, Index=I
Key : Lower : Value                 : Upper : Fixed : Stale : Domain
0 :  None :   0.03795367428994332 :  None : False : False :  Reals
1 :  None :  0.030099796252526297 :  None : False : False :  Reals
88 :  None :   0.04860152062781719 :  None : False : False :  Reals

epsilon : residual
Size=89, Index=I
Key : Lower : Value                  : Upper : Fixed : Stale : Domain
0 :  None :    0.03591704717267953 :  None : False : False :  Reals
1 :  None :   0.026035484542682643 :  None : False : False :  Reals
88 :  None :    0.07738478945809676 :  None : False : False :  Reals

§.§.§ Corrected CNLS model

Corrected convex nonparametric least squares (C$^2$NLS) is a variant of the corrected ordinary least squares (COLS) model, in which nonparametric least squares subject to monotonicity and concavity constraints replace the first-stage parametric ordinary least squares (OLS) regression. To estimate the production function, the C$^2$NLS model assumes that the regression $f$ is monotonic increasing and globally concave production function, the inefficiencies $\varepsilon$ are identically and independently distributed with mean $\mu$ and a finite variance $\sigma^2$, and that the inefficiencies $\varepsilon$ are uncorrelated with inputs $\bm{x}$.

Similar to COLS, the C$^2$NLS method includes two stages, which can be stated as follows:

* Estimate $\E[y_i \, | \, x_i] $ by solving the CNLS model, e.g., Problem (<ref>). Denote the CNLS residuals by $\varepsilon^{CNLS}_i$.

* Shift the residuals analogous to the COLS procedure; the C$^2$NLS efficiency estimator is

$$\hat{\varepsilon}_i^{C2NLS}= \varepsilon_i^{CNLS} - \max_j \varepsilon_j^{CNLS} $$

where values of $\hat{\varepsilon_i}^{C2NLS}$ range from $[0, +\infty]$ with 0 indicating efficient performance. Similarly, we adjust the CNLS intercepts $\alpha_i$ as

$$ \hat{\alpha}_i^{C2NLS}= \alpha_i^{CNLS} + \max_j \varepsilon_j^{CNLS} $$

where $\alpha^{CNLS}_i$ is the optimal intercept for firm $i$ in above CNLS problem and $\hat{\alpha}_i^{C2NLS}$ is the C$^2$NLS estimator. Slope coefficients $\bm{\beta}_i$ for C$^2$NLS are obtained directly as the optimal solution to the CNLS problem.

To solve the C$^2$NLS model (e.g., additive production model), we only need to introduce two new functions: .get_adjusted_residual() and .get_adjusted_alpha(). After completing the first-stage estimation, we use these two functions to obtain the adjusted residuals and intercept terms.

>>> from pystoned import CNLS
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data = load_Finnish_electricity_firm(x_select=['OPEX', 'CAPEX'],
...     y_select=['Energy'])
>>> model=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> print(model.get_adjusted_residual())
>>> print(model.get_adjusted_alpha())

[ -682.18627782  -677.969821    -701.60761871 -1030.2948489
-693.08377375  ...            -685.44491325 -680.26903767]

[ 656.44793389  656.39703997  656.51998647  712.844865
656.44540447  ...           656.38925598  656.52350009]

§.§ Convex Quantile and Expectile Regression

While CNLS estimates the conditional mean $\E(y_i \, | \, x_i)$, quantile regression aims at estimating the conditional median or other quantiles of the response variable ([Koenker and Bassett, 1978], [Koenker and Bassett, 1978]; [Koenker, 2005], [Koenker, 2005]) and provides an overall picture of the conditional distributions. The quantile $\tau$ splits the observations $\tau$% above and (1-$\tau$)% below. In this subsection, we extend CNLS problem to estimate CQR ([Wang et al., 2014], [Wang et al., 2014]; [Kuosmanen et al., 2015], [Kuosmanen et al., 2015]) and convex expectile regression (CER) ([Kuosmanen et al., 2020], [Kuosmanen et al., 2020]; [Dai et al., 2020], [Dai et al., 2020]; [Kuosmanen and Zhou, 2021], [Kuosmanen and Zhou, 2021]). Note that both quantile and expectile estimations are more robust to outliers and heteroscedasticity than the CNLS estimation.

§.§.§ Convex quantile regression

Given a pre-specified quantile $\tau \in (0, 1)$, the CQR model is formulated as

y_i=α_i+ β_i^'x_i+ε_i^+ - ε_i^-        ∀i

α_i+β_i^'x_i ≤α_j+β_j^'x_i         ∀i, j ,  and   i ≠j

β_i≥0        ∀i

ε_i^+≥0, ε_i^- ≥0        ∀i
where $\varepsilon^{+}_i$ and $\varepsilon^{-}_i$ denote the two non-negative components. The objective function in (<ref>) minimizes the asymmetric absolute deviations from the function rather than the symmetric quadratic deviations. The last set of constraints is the sign constraint of the error terms. The other constraints are the same as those of the CNLS problem (<ref>).

§.§.§ Convex expectile regression

Convex quantile regression (<ref>) may suffer from non-uniqueness due to that Problem (<ref>) is a linear programming (LP) problem ([Kuosmanen et al., 2015], [Kuosmanen et al., 2015]). To address this problem, Kuosmanen et al., 2015 purpose a CER approach, where a quadratic objective function is used to ensure unique estimates of the quantile functions. Consider the following QP problem

τ̃ ∑_i=1^n(ε_i^+)^2+(1-τ̃ )∑_i=1^n(ε_i^-)^2

y_i=α_i+ β_i^'x_i+ε_i^+-ε_i^-          ∀i

α_i+β_i^'x_i ≤α_j+β_j^'x_i         ∀i, j ,  and   i ≠j

β_i≥0          ∀i

ε_i^+≥0, ε_i^- ≥0         ∀i
where the expectile $\tilde{\tau} \in (0, 1)$ can be converted from/to the quantile $\tau$ ([Kuosmanen and Zhou, 2021], [Kuosmanen and Zhou, 2021]).

The alternative module CQER includes the function CQR(y, x, ...) that is designed for solving CQR problem and the function CER(y, x, ...) that is for CER problem. Therefore, we use the CQER.CQR() and CQER.CER() to define the CQR problem (see Line 5 in following example) and the CER problem, respectively. The other parameters settings are similar to those in module CNLS(y, x, ...). To display the estimated $\varepsilon_i^{+}$ and $\varepsilon_i^{-}$, the functions .display_positive_residual() and .display_negative_residual() are designed in the new module CQER.[Note that as shown in Problem <ref>, both residuals $\varepsilon _i^{+}$ and $\varepsilon_i^{-}$ are larger than 0. The function names, positive and negative residuals, are simply used to keep the consistency with the math notations in Problem <ref>.
The following additive CQR model is presented to estimate a quantile production function.

>>> from pystoned import CQER
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned import dataset as dataset
>>> data=dataset.load_GHG_abatement_cost(x_select=['HRSN','CPNK','GHG'],
...     y_select=['VALK'])
>>> model=CQER.CQR(y=data.y, x=data.x, tau=0.5, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> model.display_alpha()
>>> model.display_beta()
>>> model.display_positive_residual()
>>> model.display_negative_residual()

alpha : alpha
Size=168, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -1306.0485023191532 :  None : False : False :  Reals
1 :  None : -1306.0485023191268 :  None : False : False :  Reals
167 :  None :   170680.2112574062 :  None : False : False :  Reals

beta : beta
Size=504, Index=beta_index
Key      : Lower : Value                   : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 :      21.043360329998244 :  None : False : False :  Reals
(0, 1) :   0.0 :      0.1758315887019594 :  None : False : False :  Reals
(0, 2) :   0.0 :                     0.0 :  None : False : False :  Reals
(167, 0) :   0.0 :       13.22794884204182 :  None : False : False :  Reals
(167, 1) :   0.0 :     0.11023334987523721 :  None : False : False :  Reals
(167, 2) :   0.0 :       645.5693304144419 :  None : False : False :  Reals

epsilon_plus : positive error term
Size=168, Index=I
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :   0.0 :                0.0 :  None : False : False :  Reals
1 :   0.0 :                0.0 :  None : False : False :  Reals
167 :   0.0 :  98359.68518583104 :  None : False : False :  Reals

epsilon_minus : negative error term
Size=168, Index=I
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :   0.0 :  16157.82681402145 :  None : False : False :  Reals
1 :   0.0 :  7158.090538541088 :  None : False : False :  Reals
167 :   0.0 :                0.0 :  None : False : False :  Reals

§.§ Contextual Variables

A firm's ability to operate efficiently often depends on operational conditions and practices, such as the production environment and the firm-specific characteristics (i.e., technology selection and managerial practices). Johnson and Kuosmanen, 2011, Johnson and Kuosmanen, 2012 refer to both variables that characterize operational conditions and practices as contextual variables.

* Contextual variables are often (but not always) external factors that are beyond the control of firms

* Examples: competition, regulation, weather, location

* Need to adjust efficiency estimates for the operating environment

* Policymakers may influence the operating environment

* Contextual variables can also be internal factors

* Examples: management practices, ownership

* Better understanding of the impacts of internal factors can help the firm to improve performance

By introducing the contextual variables $\bm{z}_i$, the multiplicative model (<ref>) is reformulated as a partial log-linear model to take the operational conditions and practices into account.

\begin{align}
\ln y_i = \ln f(\bm{x_i}) + \bm{\lambda}^{'}\bm{z}_i + v_i - u_i
\end{align}

where parameter vector $\bm{\lambda}=(\lambda_1, \cdots, \lambda_r)$ represents the marginal effects of contextual variables on output. All other variables maintain their previous definitions. Similarly, we can also introduce the contextual variables to the additive model. In this subsection, we consider the multiplicative production model as our starting point.

§.§.§ CNLS with z variables

Following Johnson and Kuosmanen, 2011, we incorporate the contextual variables in the multiplicative CNLS model and redefine it as follows:

α, β, λ, ε min ∑_i=1^nε_i^2

lny_i = ln(ϕ_i+1) + λ^'z_i + ε_i          ∀i

ϕ_i  = α_i + β_i^'x_i -1         ∀i

α_i + β_i^'x_i ≤α_j + β_j^'x_i           ∀i, j ,  and   i ≠j

β_i ≥0          ∀i
Denote by $\hat{\bm{\lambda}}$ the coefficients of the contextual variables obtained as the optimal solution to the above nonlinear problem.  Johnson and Kuosmanen, 2011 examine the statistical properties of this estimator in detail, showing its unbiasedness, consistency, and asymptotic efficiency.

The contextual variables $\bm{z}$ have been integrated into the modules CNLS() and CQER(). Further, the function .display_lamda() is used to display the marginal effect of contextual variable. In the following example, we estimate a log-transformed cost function model with z-variable. Note that we assume that the firms are constant returns to scale in this example.

>>> from pystoned import CNLS
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_CRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'], z_select=['PerUndGr'])
>>> model=CNLS.CNLS(y=data.y, x=data.x, z=data.z,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_CRS)
>>> model.optimize('email@address')
>>> model.display_lamda()

lamda : z coefficient
Size=1, Index=K
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :  None : 0.3609227173656651 :  None : False : False :  Reals

§.§.§ CER with z variables

Following Kuosmanen et al., 2021, we also incorporate the contextual variable in the multiplicative CER estimation. The reformulation of CER model is

α,β,λ, ε^+,ε^-min
τ̃ ∑_i=1^n(ε_i^+)^2+(1-τ̃ )∑_i=1^n(ε_i^-)^2

lny_i = ln(ϕ_i+1) + λ^'z_i + ε_i^+-ε_i^-          ∀i

ϕ_i  = α_i + β_i^'x_i -1         ∀i

α_i + β_i^'x_i ≤α_j + β_j^'x_i           ∀i, j ,  and   i ≠j

β_i ≥0          ∀i

ε_i^+≥0, ε_i^- ≥0         ∀i
The following code is prepared to solve the CER model with z-variable. We now use the function CQER.CER() to model the CER with z-variable problem and assume that the expectile $\tilde{\tau}$ is equal to 0.5. In this example, we also estimate a CRS cost function.

>>> from pystoned import CQER
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_CRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'], z_select=['PerUndGr'])
>>> model=CQER.CER(y=data.y, x=data.x, z=data.z, tau=0.5,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_CRS)
>>> model.optimize('email@address')
>>> model.display_lamda()

lamda : z coefficient
Size=1, Index=K
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :  None : 0.3609226375022837 :  None : False : False :  Reals

§.§ Multiple Outputs (DDF Formulation)

§.§.§ CNLS with multiple outputs

Until now, the convex regression approaches have been presented within the single output, multiple input framework. In this subsection, we describe the CNLS/CQR/CER approaches with the DDF to handle multiple-input multiple-output data ([Chambers et al., 1996], [Chambers et al., 1996], [Chambers et al., 1998]).

Consider the following QP problem ([Kuosmanen and Johnson, 2017], [Kuosmanen and Johnson, 2017])

α, β, γ, εmin   ∑_i=1^nε_i^2

γ_i^'y_i = α_i + β_i^'x_i - ε_i          ∀i

α_i + β_i^'x_i -γ_i^'y_i ≤α_j + β_j^'x_i -γ_j^'y_i            ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x  = 1           ∀i

β_i ≥0, γ_i ≥0          ∀i
where the residual $\varepsilon_i$ represents the estimated value of $d$ ($\vec{D}(x_i,y_i,g^x,g^y)+u_i$). In addition to the same notations as the CNLS estimator, we also introduce firm-specific coefficients $\bm{\gamma_i}$ that represent marginal effects of outputs to the DDF.

The first constraint defines the distance to the frontier as a linear function of inputs and outputs. The linear approximation of the frontier is based on the tangent hyperplanes, analogous to the original CNLS formulation. The second set of constraints is the system of Afriat inequalities that impose global concavity. The third constraint is a normalization constraint that ensures the translation property. The last two constraints impose monotonicity in all inputs and outputs. It is straightforward to show that the CNLS estimator of function $d$ satisfies the axioms of free disposability, convexity, and translation property.

To perform the DDF related models, the pyStoNED includes the modules CNLSDDF(y, x, ...) and CQERDDF(y, x, ...), and takes the undesirable outputs into account. We first present the CNLS-DDF models in the following two examples and then describe the CQR/CER-DDF models in next subsection. To apply the module CNLSDDF(y, x, ...), we have to pre-define the directional vector for the parameters gx, gb (None; default)  and gy (i.e., Line 5: gx=[1.0, 0.0], gb=None, gy=[0.0, 0.0, 0.0]). The module also reports the estimates using .display_alpha(), .display_beta(), .display_gamma(), and .display_residual().

>>> from pystoned import CNLSDDF
>>> from pystoned.constant import FUN_PROD, OPT_LOCAL
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['OPEX','CAPEX'],
...     y_select=['Energy','Length','Customers'])
>>> model=CNLSDDF.CNLSDDF(y=data.y, x=data.x, b=None, fun=FUN_PROD,
...     gx=[1.0, 0.0], gb=None, gy=[0.0, 0.0, 0.0])
>>> model.optimize(OPT_LOCAL)
>>> model.display_alpha()
>>> model.display_beta()
>>> model.display_gamma()
>>> model.display_residual()

alpha : alpha
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -231.06893282325376 :  None : False : False :  Reals
1 :  None :  -212.5240841137753 :  None : False : False :  Reals
88 :  None :  -330.8605696176572 :  None : False : False :  Reals

beta : beta
Size=178, Index=beta_index
Key     : Lower : Value                  : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 :                    1.0 :  None : False : False :  Reals
(0, 1) :   0.0 : 0.00027176268397506657 :  None : False : False :  Reals
(88, 0) :   0.0 :                    1.0 :  None : False : False :  Reals
(88, 1) :   0.0 :     0.0526234540690525 :  None : False : False :  Reals

gamma : gamma
Size=267, Index=gamma_index
Key     : Lower : Value                  : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 :      3.806267544117393 :  None : False : False :  Reals
(0, 1) :   0.0 :     0.2311214613059167 :  None : False : False :  Reals
(0, 2) :   0.0 :   0.004895753020650473 :  None : False : False :  Reals
(88, 0) :   0.0 :   0.015392372124732883 :  None : False : False :  Reals
(88, 1) :   0.0 :  0.0029598467042669355 :  None : False : False :  Reals
(88, 2) :   0.0 :    0.06644931421296872 :  None : False : False :  Reals

epsilon : residuals
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None :  -62.41627631290396 :  None : False : False :  Reals
1 :  None : -176.54976937092056 :  None : False : False :  Reals
88 :  None : -10.336197186462073 :  None : False : False :  Reals

When considering undesirable outputs, the CNLS-DDF problem (<ref>) can be reformulated as

α, β, γ, δ, εmin   ∑_i=1^nε_i^2

γ_i^'y_i = α_i + β_i^'x_i + δ_i^'b_i - ε_i          ∀i

α_i + β_i^'x_i + δ_i^'b_i -γ_i^'y_i ≤α_j + β_j^'x_i + δ_j^'b_i -γ_j^'y_i           ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x + δ_i^'g^b = 1           ∀i

β_i ≥0, δ_i ≥0, γ_i ≥0           ∀i
where the coefficients $\bm{\delta_i}$ denote marginal effects of undesirable outputs to the DDF.

We can also model the undesirable outputs ($\bm{b}$) in the framework of DDF via the module CNLSDDF (Line 5; cf., the  example above). The estimated coefficients can be displayed by using the function .display_delta().

>>> from pystoned import CNLSDDF
>>> from pystoned.constant import FUN_PROD, OPT_LOCAL
>>> from pystoned import dataset as dataset
>>> data = dataset.load_GHG_abatement_cost()
>>> model = CNLSDDF.CNLSDDF(y=data.y, x=data.x, b=data.b,
...     fun=FUN_PROD, gx=[0.0, 0.0], gb=-1.0, gy=1.0)
>>> model.optimize(OPT_LOCAL)
>>> model.display_delta()

delta : delta
Size=168, Index=delta_index
Key      : Lower : Value               : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 :  27.122300643881694 :  None : False : False :  Reals
(1, 0) :   0.0 :  31.329802702797316 :  None : False : False :  Reals
(167, 0) :   0.0 :    6.10497985046203 :  None : False : False :  Reals

§.§.§ CQR/CER with multiple outputs

Similar to CNLS with DDF, we present another two approaches integrating DDF to convex quantile/expectile regression, and we also consider modeling the undesirable outputs.

* without undesirable outputs

* CQR-DDF model

γ_i^'y_i = α_i + β_i^'x_i + ε^+_i - ε^-_i          ∀i

α_i + β_i^'x_i -γ_i^'y_i ≤α_j + β_j^'x_i -γ_j^'y_i           ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x  = 1           ∀i

β_i ≥0 , γ_i ≥0           ∀i

ε_i^+≥0, ε_i^- ≥0          ∀i

* CER-DDF model

τ̃ ∑_i=1^n(ε_i^+)^2+(1-τ̃ )∑_i=1^n(ε_i^-)^2

γ_i^'y_i = α_i + β_i^'x_i + ε^+_i - ε^-_i          ∀i

α_i + β_i^'x_i -γ_i^'y_i ≤α_j + β_j^'x_i -γ_j^'y_i           ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x  = 1           ∀i

β_i ≥0 , γ_i ≥0          ∀i

ε_i^+≥0, ε_i^- ≥0          ∀i

Similar to the module CNLSDDF, the module CQERDDF uses the same parameters and is applied to estimate the CQR/CER-DDF model. The results can be reported using the functions that have been introduced in module CQER.CQR(y, x, ...)/.CER(y,x,...). For instance, the functions .display_positive_residual() and .display_negative_residual() are used to display the estimated residuals, respectively (Lines 8 and 9).

>>> from pystoned import CQERDDF
>>> from pystoned.constant import FUN_PROD, OPT_LOCAL
>>> from pystoned import dataset as dataset
>>> data = dataset.load_Finnish_electricity_firm(x_select=['OPEX', 'CAPEX'],
...     y_select=['Energy', 'Length', 'Customers'])
>>> model = CQERDDF.CQRDDF(y=data.y, x=data.x, b=None, tau=0.9,
...     fun = FUN_PROD, gx= [1.0, 0.0], gb=None, gy= [0.0, 0.0, 0.0])
>>> model.optimize(OPT_LOCAL)
>>> model.display_positive_residual()
>>> model.display_negative_residual()

epsilon_plus : positive error term
Size=89, Index=I
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :   0.0 :                0.0 :  None : False : False :  Reals
1 :   0.0 :                0.0 :  None : False : False :  Reals
88 :   0.0 :                0.0 :  None : False : False :  Reals

epsilon_minus : negative error term
Size=89, Index=I
Key : Lower : Value              : Upper : Fixed : Stale : Domain
0 :   0.0 : 119.91091883695128 :  None : False : False :  Reals
1 :   0.0 :               -0.0 :  None : False : False :  Reals
88 :   0.0 : 224.05081364484738 :  None : False : False :  Reals

* with undesirable outputs

* CQR-DDF-b model

α,β,γ,δ, ε^+,ε^-min

γ_i^'y_i = α_i + β_i^'x_i + δ_i^'b_i + ε^+_i - ε^-_i          ∀i

α_i + β_i^'x_i + δ_i^'b_i -γ_i^'y_i ≤α_j + β_j^'x_i + δ_j^'b_i -γ_j^'y_i          ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x + δ_i^'g^b = 1          ∀i

β_i ≥0, δ_i ≥0, γ_i ≥0           ∀i

ε_i^+≥0, ε_i^- ≥0          ∀i

* CER-DDF-b model

α,β, γ,δ, ε^+,ε^-min
τ̃ ∑_i=1^n(ε_i^+)^2+(1-τ̃ )∑_i=1^n(ε_i^-)^2

γ_i^'y_i = α_i + β_i^'x_i + δ_i^'b_i + ε^+_i - ε^-_i           ∀i

α_i + β_i^'x_i + δ_i^'b_i -γ_i^'y_i ≤α_j + β_j^'x_i + δ_j^'b_i -γ_j^'y_i            ∀i, j ,  and   i ≠j

γ_i^' g^y  + β_i^' g^x + δ_i^'g^b = 1            ∀i

β_i ≥0, δ_i ≥0, γ_i ≥0            ∀i

ε_i^+≥0, ε_i^- ≥0           ∀i

The following code shows how to solve the CQR-DDF model when the data include the undesirable output (i.e., GHG emissions).

>>> from pystoned import CQERDDF
>>> from pystoned.constant import FUN_PROD, OPT_LOCAL
>>> from pystoned import dataset as dataset
>>> data = dataset.load_GHG_abatement_cost()
>>> model = CQERDDF.CQRDDF(y=data.y, x=data.x, b=data.b, tau=0.9,
...     fun = FUN_PROD, gx= [0.0, 0.0], gb=[-1], gy=[1])
>>> model.optimize(OPT_LOCAL)
>>> model.display_delta()

delta : delta
Size=168, Index=delta_index
Key      : Lower : Value                   : Upper : Fixed : Stale : Domain
(0, 0) :   0.0 : -3.1705034323517883e-15 :  None : False : False :  Reals
(1, 0) :   0.0 :  3.2421934988654777e-15 :  None : False : False :  Reals
(167, 0) :   0.0 :      0.1270833999407421 :  None : False : False :  Reals

§.§ Relaxing convexity

§.§.§ Isotonic CNLS

This section introduces a variant of the CNLS estimator, Isotonic CNLS that relies only on the monotonicity assumption. To relax the concavity assumption in CNLS estimation (i.e., estimating a production function), we have to rephrase the Afriat inequality constraint in Problem (<ref>).

Define a binary matrix $\bm{P}=\big[p_{ij} \big]_{n \times n}$ to represent isotonicity ([Keshvari and Kuosmanen, 2013], [Keshvari and Kuosmanen, 2013]) as follows.

p_ij =

1       if  x_i ≼x_j

0       otherwise

Applying the enumeration method to define the elements of matrix $\bm{P}$, we then solve the following QP problem

α, β, ε min   ∑_i=1^nε_i^2

y_i = α_i + β_i^'x_i + ε_i          ∀i

p_ij(α_i + β_i^'x_i) ≤p_ij(α_j + β_j^'x_i)            ∀i, j ,  and   i ≠j

β_i ≥0            ∀i
Note that the concavity constraints between units $i$ and $j$ are relaxed by the matrix $P_{ij}=0$. If the $P_{ij}=1$ for all $i$ and $j$, then the above QP problem (i.e., ICNLS problem) reduces to the CNLS problem.

To calculate the monotonic models, the pyStoNED package provides the modules ICNLS(y, x, ...) and ICQER(y, x, ...), of which inherit the parameter settings from the modules CNLS(y, x, ...) and CQER(y, x, ...). Therefore, the usages of CNLS() and CQER() are similar to those of CNLS() and CQER(). Note that the matrix $\bm{P}$ is added in the module as the internal function. The readers could check the source code from the GitHub repository.

>>> from pystoned import ICNLS
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['OPEX','CAPEX'],
...     y_select=['Energy'])
>>> model=ICNLS.ICNLS(y=data.y, x=data.x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> model.display_residual()

epsilon : residual
Size=89, Index=I
Key : Lower : Value                   : Upper : Fixed : Stale : Domain
0 :  None :   0.0006849098000571985 :  None : False : False :  Reals
1 :  None :      -3.833461767868812 :  None : False : False :  Reals
88 :  None :      7.8890389440795445 :  None : False : False :  Reals

§.§.§ Isotonic CQR/CER

Similar to ICNLS, the Isotonic CQR and CER approaches are defined as follows:

* ICQR estimator

α, β, ε^+,ε^-min

y_i=α_i+ β_i^'x_i+ε_i^+-ε_i^-           ∀i

p_ij(α_i+β_i^'x_i) ≤p_ij(α_j+β_j^'x_i)            ∀i, j ,  and   i ≠j

β_i≥0           ∀i

ε_i^+≥0, ε_i^- ≥0           ∀i

>>> from pystoned import ICQER
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['OPEX','CAPEX'],
...     y_select=['Energy'])
>>> model = ICQER.ICQR(y=data.y, x=data.x, tau = 0.9, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> model.display_residual()

epsilon : error term
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -2.0000000585763473 :  None : False : False :  Reals
1 :  None : -15.000000004675286 :  None : False : False :  Reals
88 :  None :  -25.00000000488575 :  None : False : False :  Reals

* ICER estimator

α, β, ε^+,ε^-min
τ̃ ∑_i=1^n(ε_i^+)^2+(1-τ̃ )∑_i=1^n(ε_i^-)^2

y_i=α_i+ β_i^'x_i+ε_i^+-ε_i^-          ∀i

p_ij(α_i+β_i^'x_i) ≤p_ij(α_j+β_j^'x_i)            ∀i, j ,  and   i ≠j

β_i≥0            ∀i

ε_i^+≥0, ε_i^- ≥0            ∀i

>>> from pystoned import ICQER
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> data=load_Finnish_electricity_firm(x_select=['OPEX','CAPEX'],
...     y_select=['Energy'])
>>> model=ICQER.ICER(y=data.y, x=data.x, tau=0.9, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model.optimize(OPT_LOCAL)
>>> model.display_residual()

epsilon : error term
Size=89, Index=I
Key : Lower : Value                : Upper : Fixed : Stale : Domain
0 :  None :  -17.735777393077573 :  None : False : False :  Reals
1 :  None :  -30.733601094693068 :  None : False : False :  Reals
88 :  None :   -9.466741932963199 :  None : False : False :  Reals

§ STOCHASTIC NONPARAMETRIC ENVELOPMENT OF DATA

Combining virtues of SFA and DEA in a unified framework, StoNED ([Kuosmanen, 2006], [Kuosmanen, 2006]; [Kuosmanen and Kortelainen, 2012], [Kuosmanen and Kortelainen, 2012]) uses a composed error term to model both inefficiency $u$ and noise $v$ without assuming a functional form of $f$. Analogous to the COLS/C$^2$NLS estimators, the StoNED estimator consists of the following four steps:

* Step 1: Estimating the conditional mean $\E[y_i \,| \, \bm{x}_i]$ using CNLS estimator

* Step 2: Estimating the expected inefficiency $\mu$ based on the residual $\varepsilon_i^{CNLS}$

* Step 3: Estimating the StoNED frontier $\hat{f}^{StoNED}$ based on the $\hat{\mu}$

* Step 4: Estimating firm-specific inefficiencies $\E[u_i \mid \varepsilon_i^{CNLS}]$

Besides the CNLS estimator, we can apply other convex regression approaches such as ICNLS and CNLS-DDF to estimate the conditional mean in the first step (see [Keshvari and Kuosmanen, 2013], [Keshvari and Kuosmanen, 2013]; [Kuosmanen and Johnson, 2017], [Kuosmanen and Johnson, 2017]). However, the quantile and expectile related estimators introduced in Section <ref> can not be integrated into the StoNED framework at present.

§.§ Estimating the expected inefficiency

After obtaining the residuals (e.g., $\hat{\varepsilon}_i^{CNLS}$) from the convex regression approaches, one can estimate the expected value of the inefficiency term $\mu = E(u_i)$. In practice, three commonly used methods are available to estimate the expected inefficiency $\mu$: method of moments ([Aigner et al., 1977], [Aigner et al., 1977]), quasi-likelihood estimation ([Fan et al., 1996], [Fan et al., 1996]), and kernel deconvolution estimation ([Hall and Simar, 2002], [Hall and Simar, 2002]). We next briefly review these three approaches and demonstrate the application of pyStoNED; see more detailed theoretical introduction in [Kuosmanen et al., 2015].

§.§.§ Method of moments

The method of moments requires some additional parametric distributional assumptions. Following Kuosmanen et al., 2015, under the maintained assumptions of half-normal inefficiency (i.e., $u_i \sim N^+(0, \sigma_u^2)$) and normal noise (i.e., $v_i \sim N(0, \sigma_v^2)$), the second and third central moments of the composite error (i.e., $\varepsilon_i$) distribution are given by

\begin{align}
M_2 &= \bigg[\frac{\pi-2}{\pi}\bigg] \sigma_u^2 + \sigma_v^2  \notag \\
M_3 &= \bigg(\sqrt{\frac{2}{\pi}}\bigg)\bigg[1-\frac{4}{\pi}\bigg]\sigma_u^2 \notag
\end{align}

The second and third central moments can be estimated by using the CNLS residuals, i.e.,  $\hat{\varepsilon}_i^{CNLS}$:

\begin{align}
\hat{M_2} &= \sum_{i=1}^{n}(\hat{\varepsilon}_i-\bar{\varepsilon})^{2}/n  \notag \\
\hat{M_3} &= \sum_{i=1}^{n}(\hat{\varepsilon}_i-\bar{\varepsilon})^{3}/n  \notag
\end{align}

Note that the third moment $M_3$ (which measures the skewness of the distribution) only depends on the standard deviation parameter $\sigma_u$ of the inefficiency distribution. Thus, given the estimated $\hat{M}_3$ (which should be positive in the case of a cost frontier), we can estimate the parameters $\sigma_u$ and $\sigma_v$ by

\begin{align}
\hat{\sigma}_u &= \sqrt[3]{\frac{\hat{M_3}}{\bigg(\sqrt{\frac{2}{\pi}}\bigg)\bigg[1-\frac{4}{\pi}\bigg]}} \notag \\
\hat{\sigma}_v &= \sqrt[2]{\hat{M_2}-\bigg[\frac{\pi-2}{\pi}\bigg] \hat{\sigma}_u^2 } \notag
\end{align}

To estimate the unconditional expected inefficiency $\mu$ and firm-specific inefficiencies, we provide a module StoNED() that inherits the parameter setting from the module CNLS(). We estimate the conditional expected inefficiency through the presented example below. Following the calculation procedure of StoNED estimator, we first utilize the module CNLS() to estimate the conditional mean $\E[y_i \,| \, \bm{x}_i]$ (i.e., Lines 5-6), then apply the module StoNED(model), where the parameter model is defined as the name of CNLS model (cf. footnote <ref>), to decompose the estimated residuals (i.e., Line 8). We develop and import the parameter RED_MOM (i.e., Line 3) to decompose the residuals by using the method of moment. We finally resort to the function .get_unconditional_expected_inefficiency(RED_MOM)) included in the module StoNED() to retrieve the expected inefficiency $\hat{\mu}$.

>>> from pystoned import CNLS, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_MOM
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model=CNLS.CNLS(data.y, data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_unconditional_expected_inefficiency(RED_MOM))

§.§.§ Quassi-likelihood estimation

Quasi-likelihood approach is an alternative approach to decomposing $\sigma_u$ and $\sigma_v$ suggested by Fan et al., 1996. Given the shape of the CNLS curve, we apply the standard maximum likelihood method to estimate the parameters $\sigma_u$ and $\sigma_v$. The quasi-likelihood function is formulated as

\begin{align}
\ln L(\lambda) & = -n\ln(\hat{\sigma}) + \sum \ln\Phi\bigg[\frac{-\hat{\varepsilon}_i \lambda}{\hat{\sigma}}\bigg] - \frac{1}{2\hat{\sigma}^2}\sum\hat{\varepsilon}_i^2  \notag
\end{align}

\begin{align}
\hat{\varepsilon}_i &= \hat{\varepsilon}_i^{CNLS}-(\sqrt{2}\lambda\hat{\sigma})/[\pi(1+\lambda^2)]^{1/2}   \notag  \\
\hat{\sigma} &= \Bigg\{\frac{1}{n}\sum(\hat{\varepsilon}_i^{CNLS})^2 / \bigg[1 - \frac{2\lambda^2}{\pi(1+\lambda^2)}\bigg]  \Bigg\}.   \notag
\end{align}

Note that the quasi-likelihood function only consists of a single parameter $\lambda$ (i.e., the signal-to-noise ratio $\lambda = \sigma_u/\sigma_v$).[Note the notation difference between signal-to-noise ratio $\lambda$ and the marginal effect of contextual variable $\lambda$ in Problems (<ref>) and (<ref>).
The symbol $\Phi$ represents the cumulative distribution function of the standard normal distribution. In the pyStoNED package, we use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm provided by SciPy to estimate the maximum likelihood function.

Since we apply the quasi-likelihood estimation to decompose the residuals, we import the parameter RED_QLE in the module StoNED() (see Line 3).

>>> from pystoned import CNLS, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_QLE
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model = CNLS.CNLS(data.y, data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_unconditional_expected_inefficiency(RED_QLE))

§.§.§ Kernel deconvolution estimation

While the method of moments and quasi-likelihood approaches require additional distributional assumptions for the inefficiency and noise terms, an alternative nonparametric estimation of the expected inefficiency $\mu$ is available by applying nonparametric kernel deconvolution proposed by Hall and Simar, 2002. Note that the residual $\hat{\varepsilon}_i^{CNLS}$ is a consistent estimator of $e^o = \varepsilon_i + \mu$ (for production model). Following [Kuosmanen and Johnson, 2017], the density function of ${e^o}$ is

\begin{align}
\hat{f}_{e^o}(z) = (nh)^{-1} \sum_{i=1}^{n}K\bigg(\frac{z-e_i^o}{h} \bigg), \notag
\end{align}

where $K(\cdot)$ is a compactly supported kernel, and $h$ is a bandwidth. Hall and Simar, 2002 show that the first derivative of the density function of the composite error term ($f_\varepsilon^{'}$) is proportional to that of the inefficiency term ($f_u^{'}$) in the neighborhood of $\mu$. Therefore, a nonparametric estimator of expected inefficiency $\mu$ is obtained as

\begin{align}
\hat{\mu} = \arg \max_{z \in C}(\hat{f^{'}}_{e^o}(z)), \notag
\end{align}

where $C$ is a closed interval in the right tail of $f_{e^o}$.

We then import the parameter RED_KDE to decompose the estimated residuals, i.e., $\hat{\varepsilon}_i^{CNLS}$ (i.e., Line 3).

>>> from pystoned import CNLS, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_KDE
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model = CNLS.CNLS(data.y, data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_unconditional_expected_inefficiency(RED_KDE))

§.§ Estimating firm-specific inefficiencies

After estimating the expected inefficiency $\mu$ using the methods of moment (MOM) or quasi-likelihood estimation (QLE),[
For the expected inefficiency $\mu$ estimated by kernel deconvolution, Dai, 2016 proposes a non-parametric strategy where the Richardson–Lucy blind deconvolution algorithm is used to identify firm-specific inefficiencies. However, the pyStoNED package only supports the parametric estimation of firm-specific inefficiencies due to the fact that the parametric method is more widely used in efficiency analysis literature.
we then employ JLMS estimator proposed by Jondrow et al., 1982 to estimate the firm-specific inefficiencies ([Johnson and Kuosmanen, 2015], [Johnson and Kuosmanen, 2015]). Under the assumption of a normally distributed error term and a half-normally distributed inefficiency term, JLMS formulates the conditional distribution of inefficiency $u_i$, given $\hat{\varepsilon}_i$, and proposes the inefficiency estimator as the conditional mean $\E[u_i \, | \, \hat{\varepsilon}_i]$.

Following Kumbhakar and Lovell, 2000, the conditional expected value of inefficiency $\E[u_i  \, | \, \hat{\varepsilon}_i]$ for production function and cost function are shown as follows, respectively:

* Production function

\begin{align}
\E[u_i  \, | \, \hat{\varepsilon}_i] = \mu_{*i} + \sigma_* \Bigg[ \frac{\phi(-\mu_{*i}/\sigma_*)}{1-\Phi(-\mu_{*i}/\sigma_*)} \Bigg] = \sigma_* \Bigg[ \frac{\phi(\varepsilon_i \lambda/\sigma)}{1-\Phi(\varepsilon_i \lambda/\sigma)} - \frac{\varepsilon_i \lambda}{\sigma} \Bigg]  \notag
\end{align}

where $\mu_*= -\varepsilon \sigma_u^2/\sigma^2$, $\sigma_*^2 = \sigma_u^2\sigma_v^2/\sigma^2$, $\lambda = \sigma_u/\sigma_v$, and $\sigma^2 = \sigma_u^2 +\sigma_v^2$. The symbol $\phi$ is the standard normal density function, and the symbol $\Phi$ denotes the cumulative distribution function of the standard normal distribution.

* Cost function

\begin{align}
\E[u_i \mid \hat{\varepsilon}_i] = \mu_{*i} + \sigma_* \Bigg[ \frac{\phi(-\mu_{*i}/\sigma_*)}{1-\Phi(-\mu_{*i}/\sigma_*)} \Bigg] = \sigma_* \Bigg[ \frac{\phi(\varepsilon_i \lambda/\sigma)}{1-\Phi(-\varepsilon_i \lambda/\sigma)} + \frac{\varepsilon_i \lambda}{\sigma} \Bigg]  \notag
\end{align}

where $\mu_*= \varepsilon \sigma_u^2/\sigma^2$, $\sigma_*^2 = \sigma_u^2\sigma_v^2/\sigma^2$, $\lambda = \sigma_u/\sigma_v$, and $\sigma^2 = \sigma_u^2 +\sigma_v^2$.

The firm-level technical efficiency (TE) is then measured based on the estimated conditional mean. For other models, the technical efficiency can be estimated as follows.

* Production function

* Multiplicative model: $\text{TE} = \exp(-\E[u_i \mid  \varepsilon_i])$

* Additive model: $\text{TE} = \frac{y - \E[u_i \mid  \varepsilon_i]}{y}$

* Cost function

* Multiplicative model: $\text{TE} = \exp(\E[u_i \mid  \varepsilon_i])$

* Additive model: $\text{TE} = \frac{y+ \E[u_i \mid  \varepsilon_i]}{y}$

To calculate the firm-level technical efficiency, we resort to the function .get_technical_
inefficiency(). We either set the parameter in function .get_technical_inefficiency() to RED_MOM (i.e., using the MOM approach to calculate the efficiency; Line 9) or to RED_QLE (i.e., using the quassi-likelihood estimation approach).

>>> from pystoned import CNLS, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_MOM
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model = CNLS.CNLS(data.y, data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_technical_inefficiency(RED_MOM))

[1.02974326 1.029505   1.03407182 1.03011172 1.03285089 1.02624282
1.02854464 1.03305965 ...        1.02883011 1.03059918 1.03003701]

Further, we can consider the effect of contextual variables ($\bm{z}$) in the first-step estimation of the StoNED estimator (i.e., Line 5) when calculating the technical efficiency.

>>> from pystoned import CNLS, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_MOM
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'], z_select=['PerUndGr'])
>>> model = CNLS.CNLS(y=data.y, x=data.x, z=data.z,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_technical_inefficiency(RED_MOM))

[1.03110069 1.02940641 1.02886856 1.02977033 1.03447988 1.02600838
1.02809928 1.03574493 ...        1.02743418 1.03236693 1.02948783]

§ CNLS-G ALGORITHM

Since convex regression approaches shape the convexity (concavity) of function using the Afrait inequality, the estimation becomes excessively expensive due to the $O(n^2)$ linear constraints. For example, if the data samples have 500 observations, the total number of linear constraints is equal to 250,000. To speed up the computational time, Lee et al., 2013 propose a more efficient generic algorithm, CNLS-G, which uses the relaxed Afriat constraint set and iteratively adds violated constraints to the relaxed model as necessary. See further discussions in Lee et al., 2013.

To illustrate the CNLS-G algorithm, we follow Lee et al., 2013 to generate the input and output variables. In this section, we assume an additive production function with two-input and one-output, $y=x_1^{0.4} \times x_2^{0.4}+u$. We randomly draw the inputs $x_1$ and $x_2$ from a uniform distribution, $x \sim U[1, 10]$, and the error term $u$ from a normal distribution, $u \sim N(0, 0.7^2)$. Based on these settings, we generate 500 artificial observations, estimate the CNLS problem (<ref>) and the CER problem (<ref>), and calculate the firm-level technical efficiency using the CNLS-G algorithm.

§.§ Solving CNLS model

We first compare the running time of original CNLS and CNLS-G algorithm in the same computation environment.  Line 1 imports the modules CNLSG() that is designed to perform the CNLS-G algorithm and CNLS() from the package pyStoNED. Note that the module CNLSG() has the same parameters as the module CNLS(). Line 3 imports the NumPy to provide multidimensional array and functions for linear algebra. We use the time to count the running time for CNLS estimation (i.e., Lines 4, 9,  12) and employ the function .get_runningtime() to directly obtain the running time for CNLS-G algorithm (i.e., Line 16). To replicate the experiment, we set a random seed using the np.random.seed(). Lines 6-8 generate the variables $\bm{x}$, $u$, and $y$. model1 and model2 are the normal CNLS model and CNLS-G model, respectively. To count the number of constraints included in the CNLS-G algorithm, the module CNLSG() provides an internal function .get_totalconstr() (Line 19).

>>> from pystoned import CNLSG, CNLS
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> import numpy as np
>>> import time
>>> np.random.seed(0)
>>> x=np.random.uniform(low=1, high=10, size=(500, 2))
>>> u=np.random.normal(loc=0, scale=0.7, size=500)
>>> y=x[:, 0]**0.4*x[:, 1]**0.4+u
>>> t1=time.time()
>>> model1=CNLS.CNLS(y, x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model1.optimize(OPT_LOCAL)
>>> CNLS_time = time.time()-t1
>>> model2=CNLSG.CNLSG(y, x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model2.optimize(OPT_LOCAL)
>>> print("The running time with algorithm is ", model2.get_runningtime())
>>> print("The running time without algorithm is ", CNLS_time)
>>> print("The total number of constraints is ", model2.get_totalconstr())

The running time with algorithm is  20.35739278793335
The running time without algorithm is  30.252474308013916
The total number of constraints is  9700.0

§.§ Solving CER model

We next demonstrate a CER model solved by the CNLS-G algorithm prepared in module CQERG(). The other experimental settings are similar to those in Section <ref>. Note that the CQR model can also be solved by the CNLS-G algorithm via the function CQERG.CQRG(y, x, ...) (cf. Line 14).

>>> from pystoned import CQERG, CQER
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> import numpy as np
>>> import time
>>> np.random.seed(0)
>>> x=np.random.uniform(low=1, high=10, size=(500, 2))
>>> u=np.random.normal(loc=0, scale=0.7, size=500)
>>> y=x[:, 0]**0.4*x[:, 1]**0.4+u
>>> tau=0.5
>>> t1=time.time()
>>> model1=CQER.CER(y, x, tau, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model1.optimize(OPT_LOCAL)
>>> CER_time=time.time()-t1
>>> model2=CQERG.CERG(y, x, tau, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> model2.optimize(OPT_LOCAL)
>>> print("The running time with algorithm is ", model2.get_runningtime())
>>> print("The running time without algorithm is ", CER_time)
>>> print("The total number of constraints in CER model is ",
...     model2.get_totalconstr())

The running time with algorithm is  22.845768451690674
The running time without algorithm is  29.961899280548096
The total number of constraints in CER model is  9701.0

§.§ Calculating firm-level efficiency

We can apply the CNLS-G algorithm to calculate the firm-level efficiency more efficiently. In the first step of StoNED() estimator, we use the module CNLSG(y, x, ...) as a substitute for the module CNLS(y, x, ...) (Lines 1 and 5). The other settings are similar to those in Section <ref>.

>>> from pystoned import CNLSG, StoNED
>>> from pystoned.dataset import load_Finnish_electricity_firm
>>> from pystoned.constant import CET_MULT, FUN_COST, RTS_VRS, RED_MOM
>>> data=load_Finnish_electricity_firm(x_select=['Energy','Length','Customers'],
...     y_select=['TOTEX'])
>>> model=CNLSG.CNLSG(data.y, data.x, z=None,
...     cet=CET_MULT, fun=FUN_COST, rts=RTS_VRS)
>>> model.optimize('email@address')
>>> rd = StoNED.StoNED(model)
>>> print(rd.get_technical_inefficiency(RED_MOM))

[1.09156319 1.08650469 1.14247297 1.09270023 1.12171169 1.06035128
1.07631848 1.13245426 ...        1.08719327 1.09467879 1.09408226]

§ GRAPHICAL ILLUSTRATION OF ESTIMATED FUNCTION

To illustrate how the estimated function looks like, the pyStoNED package provides the functions plot.plot2d() and plot.plot3d() to plot two- and three-dimensional (2D and 3D) estimated functions. Similar to the usage of the module StoNED(y, x, ...), we first estimate the nonparametric regression such as CNLS, CQR, and ICNLS  and then apply the plot function to draw the figures. In this section, we use the internal data provided with Tim Coelli's Frontier 4.1 to demonstrate the plotting process.

§.§ One-input and One-output

In the one-input and one-output case, we present two different estimated functions: the CNLS function and the CQR function. Therefore, we import the modules CNLS and CQER in Line 1 as the estimators, the plotting module plot2d (Line 2), and the example dataset (Line 4). Lines 6 and 7 and Lines 10 and 11 define and solve the CNLS and CQR models. Lines 8 and 12 are used to plot the estimated function. There are four parameters in the module plot2d(...): the first is the model's name (see footnote <ref>), the second parameter x_select defines which is the selected input $x$, the third and last are the given names of the legend and generated picture, respectively.

>>> from pystoned import CNLS, CQER
>>> from pystoned.plot import plot2d
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Tim_Coelli_frontier
>>> data=load_Tim_Coelli_frontier(x_select=['labour'],y_select=['output'])
>>> CNLS_model=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> CNLS_model.optimize(OPT_LOCAL)
>>> plot2d(CNLS_model, x_select=0, label_name="CNLS estimated function",
...     fig_name="CNLS_2d")
>>> CQR_model=CQER.CQR(y=data.y, x=data.x, tau=0.5, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> CQR_model.optimize(OPT_LOCAL)
>>> plot2d(CQR_model, x_select=0, label_name="CQR estimated function",
...     fig_name="CQR_2d")

Figure <ref> depicts the functions estimated by the CNLS and CQR model (with $\tau=0.5$).

[CNLS model]

[CQR model ($\tau=0.5$)]
Estimated function by CNLS/CQR model.

§.§ Two-input and One-output

In the two-input and one-output case, we use linear interpolation to obtain the 3D surface due to the fact that CNLS estimates hyperplanes at the observation points.[
Note that the CNLS estimator includes linear programming as the second stage estimation to find the minimum extrapolated production function ([Kuosmanen and Kortelainen, 2012], [Kuosmanen and Kortelainen, 2012]).
The function plot3d(model, x_select_1, x_select_2, fig_name=None, line_transparent=False,
pane_transparent=False) includes six parameters: the first is the name of the estimated model, the second and third define the selected input, the fourth is the name of generated figure, the last two are the basic settings for the figure (False; default). We import function plot3d(...) in Line 2 and plot the figure in Line 8. Figure <ref> presents the estimated 3D estimated function.

>>> from pystoned import CNLS
>>> from pystoned.plot import plot3d
>>> from pystoned.constant import CET_ADDI, FUN_PROD, OPT_LOCAL, RTS_VRS
>>> from pystoned.dataset import load_Tim_Coelli_frontier
>>> data=load_Tim_Coelli_frontier(x_select=['capital', 'labour'],
...     y_select=['output'])
>>> CNLS_model=CNLS.CNLS(y=data.y, x=data.x, z=None,
...     cet=CET_ADDI, fun=FUN_PROD, rts=RTS_VRS)
>>> CNLS_model.optimize(OPT_LOCAL)
>>> plot3d(CNLS_model,x_select_1=0, x_select_2=1, fig_name="CNLS_3d")

Estimated function by CNLS model.

§ CONCLUSIONS

Convex regression and related methods provide an appealing way to impose shape constraints implied by the theory without making restrictive functional form assumptions. As these techniques are becoming increasingly popular, there is a great demand for a powerful, reliable, and fully open access computational package. The Python package pyStoNED aims to address this need, providing a comprehensive set of functions for estimating CNLS, StoNED, and their numerous variants. For more information, see its developed repository on GitHub at <https://github.com/ds2010/pyStoNED> and documentation at <https://pystoned.readthedocs.io>.

This paper has reviewed the models and specifications currently supported by pyStoNED, and demonstrated its use by empirical examples taken from the literature of productivity and efficiency analysis. All modules are implemented in a fully open access environment. We encourage the users to utilize pyStoNED to further develop their own packages for specific estimation purposes.

In the future, our plan is to include further extensions to pyStoNED on a continuous basis. One interesting avenue of ongoing research is to include additional penalty terms to the objective function of the CNLS/CQR/CER problems to alleviate overfitting and the curse of dimensionality (e.g.,  [Lee and Cai, 2020], [Lee and Cai, 2020]). More efficient computational algorithms such as the CNLS-A proposed by Dai, 2021 are under active development, and will be included in the future editions of pyStoNED.

We hope that the pyStoNED package could contribute to raising the standards of empirical applications of productivity and efficiency analysis, and facilitate more meaningful and relevant applications that influence the managerial and policy decisions. The basic idea of the StoNED approach is to enable applied researchers and practitioners of efficiency analysis to incorporate existing tools and techniques from different domains such as econometrics, statistics, operational research, and machine learning to a logically consistent unified framework, but also to facilitate further methodological development in a multi-disciplinary environment. The purpose of the pyStoNED package is to support this development.

While we have phrased this review and the pyStoNED modules in terms of the cost and production function, most of the modules are readily applicable to nonparametric regression analysis in any other context as well. We hope that the convex regression and related techniques could also prove useful in other application areas where shape constraints play an essential role. For example, optimization behavior also implies specific convexity constraints in the context of consumer demand analysis (see, e.g., [Afriat, 1967], [Afriat, 1967]; [Varian, 1982], [Varian, 1982]).

§ ACKNOWLEDGMENTS

Sheng Dai gratefully acknowledges financial support from the Foundation for Economic Education (Liikesivistysrahasto) [No. 180019, 190073] and the HSE Support Foundation [No. 11–2290].

[Afriat, 1972]
Afriat SN (1972).
Efficiency Estimation of Production Functions.
International Economic Review, 13, 568–598.

[Afriat, 1967]
Afriat SN (1967).
The Construction of Utility Functions from Expenditure
International Economic Review, 8, 67.

[Aigner et al., 1977]
Aigner D, Lovell CK, Schmidt P (1977).
Formulation and Estimation of Stochastic Frontier
Production Function Models.
Journal of Econometrics, 6, 21–37.

[Álvarez et al., 2020]
Álvarez IC, Barbero J, Zofío JL (2020).
A Data Envelopment Analysis Toolbox for
Journal of Statistical Software, 95, 1–49.

[Badunenko et al., 2020]
Badunenko O, Mozharovskyi P, Kolomiytseva Y (2020).
npsf: Nonparametric and Stochastic Efficiency and
Productivity Analysis.
R package version 0.8-0,

[Barbero and Zofío, 2021]
Barbero J, Zofío J (2021).
DataEnvelopmentAnalysis: A Julia Package
for Efficiency and Productivity Measurement using Data Envelopment Analysis

[Belotti et al., 2013]
Belotti F, Daidone S, Ilardi G, Atella V (2013).
Stochastic Frontier Analysis using Stata.
The Stata Journal, 13, 719–758.

[Bertsimas and Mundru, 2021]
Bertsimas D, Mundru N (2021).
Sparse convex regression.
INFORMS Journal on Computing, 33, 262–279.

[Bogetoft and Otto, 2010]
Bogetoft P, Otto L (2010).
Benchmarking With DEA, SFA, and R.
Springer-Verlag, New York.

[Brunk, 1955]
Brunk HD (1955).
Maximum Likelihood Estimates of Monotone Parameters.
The Annals of Mathematical Statistics, 26, 607–616.

[Bynum et al., 2021]
Bynum ML, Hackebeil GA, Hart WE, Laird CD, Nicholson BL, Siirola JD, Watson JP,
Woodruff DL (2021).
Pyomo–Optimization Modeling in Python,
volume 67.
Third edition. Springer, Berlin.

[Byrd et al., 2006]
Byrd RH, Nocedal J, Waltz RA (2006).
KNITRO: An Integrated Package for Nonlinear
In Large-scale nonlinear optimization, pp. 35–59. Springer.

[Chambers et al., 1996]
Chambers RG, Chung Y, Färe R (1996).
Benefit and Distance Functions.
Journal of Economic Theory, 70, 407–419.

[Chambers et al., 1998]
Chambers RG, Chung Y, Färe R (1998).
Profit, Directional Distance Functions, and Nerlovian
Journal of Optimization Theory and Applications, 98,

[Charnes et al., 1978]
Charnes A, Cooper W, Rhodes E (1978).
Measuring the Efficiency of Decision Making Units.
European Journal of Operational Research, 2,

[Coelli, 1996]
Coelli T (1996).
A Guide to FRONTIER Version 4.1: A Computer Program for
Stochastic Frontier Production and Cost Function Estimation.
CEPA working paper 96/08, University of New England.

[Coelli and Henningsen, 2020]
Coelli T, Henningsen A (2020).
frontier: Stochastic Frontier Analysis.
R package version 1.1-8,

[Coelli et al., 2005]
Coelli T, Rao DSP, O'Donnell CJ, Battese GE (2005).
An Introduction to Efficiency and Productivity Analysis.
Springer, New York.

[Cplex, IBM ILOG, 2009]
Cplex, IBM ILOG (2009).
V12. 1: User’s Manual for CPLEX.
International Business Machines Corporation, 46(53),

[Dai, 2021]
Dai S (2021).
Variable Selection in Convex Quantile Regression: L1-norm
or L0-norm Regularization?

[Dai et al., 2020]
Dai S, Zhou X, Kuosmanen T (2020).
Forward-looking Assessment of the GHG Abatement Cost:
Application to China.
Energy Economics, 88, 104758.

[Dai, 2016]
Dai X (2016).
Non-parametric Efficiency Estimation Using Richardson-Lucy
Blind Deconvolution.
European Journal of Operational Research, 248,

[Dakpo et al., 2021]
Dakpo KH, Desjeux D, Latruffe L (2021).
sfaR: Stochastic Frontier Analysis Using
R package version 0.1-0,

[Daouia and Laurent, 2013]
Daouia A, Laurent T (2013).
frontiles: Partial Frontier Efficiency Analysis.
R package version 1.2,

[Fan et al., 1996]
Fan Y, Li Q, Weersink A (1996).
Semiparametric Estimation of Stochastic Production Frontier
Journal of Business & Economic Statistics, 14,

[GAMS Development Corporation, 2013]
GAMS Development Corporation (2013).
General Algebraic Modeling System (GAMS) Release
Washington, DC, USA.

[Grenander, 1956]
Grenander U (1956).
On the theory of mortality measurement.
Scandinavian Actuarial Journal, 1956(2), 125–153.

[Hall and Simar, 2002]
Hall P, Simar L (2002).
Estimating a Change Point, Boundary, or Frontier in the
Presence of Observation Error.
Journal of the American Statistical Association, 97,

[Hannah and Dunson, 2013]
Hannah LA, Dunson DB (2013).
Multivariate convex regression with adaptive
Journal of Machine Learning Research, 14, 3153–3188.

[Harris et al., 2020]
Harris CR, Millman KJ, van der Walt SJ, Gommers R, Virtanen P, Cournapeau D,
Wieser E, Taylor J, Berg S, Smith NJ, Kern R, Picus M, Hoyer S, van Kerkwijk
MH, Brett M, Haldane A, del Río JF, Wiebe M, Peterson P,
Gérard-Marchant P, Sheppard K, Reddy T, Weckesser W, Abbasi H, Gohlke
C, Oliphant TE (2020).
Array programming with NumPy.
Nature, 585, 357–362.

[Hildreth, 1954]
Hildreth C (1954).
Point Estimates of Ordinates of Concave Functions.
Journal of the American Statistical Association, 49,

[Hunter, 2007]
Hunter JD (2007).
Matplotlib: A 2D Graphics Environment.
Computing in Science & Engineering, 9, 90–95.

[Ji and Lee, 2010]
Ji Yb, Lee C (2010).
Data Envelopment Analysis.
The Stata Journal, 10, 267–280.

[Johnson and Kuosmanen, 2011]
Johnson AL, Kuosmanen T (2011).
One-stage Estimation of the Effects of Operational
Conditions and Practices on Productive Performance: Asymptotically Normal and
Efficient, Root-n Consistent StoNEZD Method.
Journal of Productivity Analysis, 36, 219–230.

[Johnson and Kuosmanen, 2012]
Johnson AL, Kuosmanen T (2012).
One-stage and two-stage DEA estimation of the effects of
contextual variables.
European Journal of Operational Research, 220,

[Johnson and Kuosmanen, 2015]
Johnson AL, Kuosmanen T (2015).
An Introduction to CNLS and StoNED Methods for Efficiency
Analysis: Economic Insights and Computational Aspects.
In SC Ray, SC Kumbhakar, P Dua (eds.), Benchmarking for
Performance Evaluation: A Production Frontier Approach, chapter 3, pp.
117–186. Springer.

[Jondrow et al., 1982]
Jondrow J, Lovell CK, Materov IS, Schmidt P (1982).
On the Estimation of Technical Inefficiency in the
Stochastic Frontier Production Function Model.
Journal of Econometrics, 19, 233–238.

[Kalvelagen, 2002]
Kalvelagen E (2002).
Efficiently Solving DEA Models with GAMS.
GAMS Development Corporation, Washington, DC.

[Keshvari and Kuosmanen, 2013]
Keshvari A, Kuosmanen T (2013).
Stochastic Non-convex Envelopment of Data: Applying
Isotonic Regression to Frontier Estimation.
European Journal of Operational Research, 231,

[Koenker, 2005]
Koenker R (2005).
Quantile Regression.
Cambridge University Press, Cambridge.

[Koenker and Bassett, 1978]
Koenker R, Bassett G (1978).
Regression Quantiles.
Econometrica, 46, 33–50.

[Kumbhakar and Lovell, 2000]
Kumbhakar SC, Lovell CK (2000).
Stochastic Frontier Analysis.
Cambridge University Press, Cambridge.

[Kuosmanen, 2006]
Kuosmanen T (2006).
Stochastic Nonparametric Envelopment of Data: Combining
Virtues of SFA and DEA in a Unified Framework.
MTT Discussion Paper No.3/2006.

[Kuosmanen, 2008]
Kuosmanen T (2008).
Representation Theorem for Convex Nonparametric Least
Econometrics Journal, 11, 308–325.

[Kuosmanen, 2012]
Kuosmanen T (2012).
Stochastic Semi-Nonparametric Frontier Estimation of
Electricity Distribution Networks: Application of the StoNED method in the
Finnish regulatory model.
Energy Economics, 34, 2189–2199.

[Kuosmanen et al., 2015]
Kuosmanen T, Johnson A, Saastamoinen A (2015).
Stochastic Nonparametric Approach to Efficiency Analysis: A
Unified Framework.
In J Zhu (ed.), Data Envelopment Analysis, chapter 7, pp.
191–244. Springer.

[Kuosmanen and Johnson, 2010]
Kuosmanen T, Johnson AL (2010).
Data Envelopment Analysis as Nonparametric Least-Squares
Operations Research, 58, 149–160.

[Kuosmanen and Johnson, 2017]
Kuosmanen T, Johnson AL (2017).
Modeling Joint Production of Multiple Outputs in StoNED:
Directional Distance Function Approach.
European Journal of Operational Research, 262,

[Kuosmanen and Kortelainen, 2012]
Kuosmanen T, Kortelainen M (2012).
Stochastic Non-smooth Envelopment of Data: Semi-parametric
Frontier Estimation Subject to Shape Constraints.
Journal of Productivity Analysis, 38, 11–28.

[Kuosmanen et al., 2013]
Kuosmanen T, Saastamoinen A, Sipiläinen T (2013).
What is the Best Practice for Benchmark Regulation of
Electricity Distribution? Comparison of DEA, SFA and StoNED Methods.
Energy Policy, 61, 740–750.

[Kuosmanen et al., 2021]
Kuosmanen T, Tan A, Dai S (2021).
Performance of English NHS Hospitals During the First and
Second Waves of the COVID-19 Pandemic.
working paper.
Available online at

[Kuosmanen and Zhou, 2021]
Kuosmanen T, Zhou X (2021).
Shadow Prices and Marginal Abatement Costs: Convex Quantile
Regression Approach.
European Journal of Operational Research, 289(2),

[Kuosmanen et al., 2020]
Kuosmanen T, Zhou X, Dai S (2020).
How Much Climate Policy Has Cost for OECD Countries?
World Development, 125, 104681.

[Lee and Cai, 2020]
Lee CY, Cai JY (2020).
LASSO Variable Selection in Data Envelopment Analysis with
Small Datasets.
Omega, 91, 102019.
ISSN 03050483.

[Lee et al., 2013]
Lee CY, Johnson AL, Moreno-Centeno E, Kuosmanen T (2013).
A More Efficient Algorithm for Convex Nonparametric Least
European Journal of Operational Research, 227,

[Lim and Glynn, 2012]
Lim E, Glynn PW (2012).
Consistency of multidimensional convex regression.
Operations Research, 60, 196–208.

[Magnani and Boyd, 2009]
Magnani A, Boyd SP (2009).
Convex piecewise-linear fitting.
Optimization and Engineering, 10, 1–17.

[Mazumder et al., 2019]
Mazumder R, Choudhury A, Iyengar G, Sen B (2019).
A computational framework for multivariate convex
regression and its variants.
Journal of the American Statistical Association, 114,

[Meeusen et al., 1977]
Meeusen W, Van J, Broeck D (1977).
Efficiency Estimation from Cobb-Douglas Production
Functions with Composed Error.
International Economic Review, 18, 435–444.

[MOSEK ApS, 2021]
MOSEK ApS (2021).
The MOSEK Optimization Toolbox for Python
Manual. Version 9.2.47.

[Murtagh and Saunders, 2003]
Murtagh B, Saunders M (2003).
MINOS 5.51 User’s Guide.

[Oh, DH and Suh, D, 2013]
Oh, DH and Suh, D (2013).
nonparaeff: Nonparametric Methods for Measuring
Efficiency and Productivity.
R package version 0.5-8,

[R Core Team, 2020]
R Core Team (2020).
R: A Language and Environment for Statistical
R Foundation for Statistical Computing, Vienna, Austria.

[Raith et al., 2021]
Raith A, Perederieieva O, Fauzi F, Harton K, Lee A, Lin K, Priddey H, Rouse M
pyDEA: Package for Conducting Data Envelopment
Python package version 1.6,

[Seijo and Sen, 2011]
Seijo E, Sen B (2011).
Nonparametric least squares estimation of a multivariate
convex regression function.
Annals of Statistics, 39, 1633–1657.

[Simm and Besstremyannaya, 2020]
Simm J, Besstremyannaya G (2020).
rDEA: Robust Data Envelopment Analysis (DEA) for
R package version 1.2-6,

[Straub, 2014]
Straub A (2014).
sfa: Stochastic Frontier Analysis.
R package version 1.0-1,

[The Mathworks, Inc., 2017]
The Mathworks, Inc (2017).
MATLAB — The Language of Technical Computing,
Version R2015b.
Natick, Massachusetts.

[The pandas development team, 2020]
The pandas development team (2020).
Pandas-dev/pandas: Pandas.

[Varian, 1982]
Varian HR (1982).
The Nonparametric Approach to Demand Analysis.
Econometrica, 50, 945.

[Varian, 1984]
Varian HR (1984).
The Nonparametric Approach to Production Analysis.
Econometrica, 52, 579.

[Virtanen et al., 2020]
Virtanen P, Gommers R, Oliphant TE, Haberland M, Reddy T, Cournapeau D,
Burovski E, Peterson P, Weckesser W, Bright J, van der Walt SJ, Brett M,
Wilson J, Millman KJ, Mayorov N, Nelson ARJ, Jones E, Kern R, Larson E, Carey
CJ, Polat İ, Feng Y, Moore EW, VanderPlas J, Laxalde D, Perktold J,
Cimrman R, Henriksen I, Quintero EA, Harris CR, Archibald AM, Ribeiro AH,
Pedregosa F, van Mulbregt P, SciPy 10 Contributors (2020).
SciPy 1.0: Fundamental Algorithms for Scientific
Computing in Python.
Nature Methods, 17, 261–272.

[Wang et al., 2014]
Wang Y, Wang S, Dang C, Ge W (2014).
Nonparametric Quantile Frontier Estimation Under Shape
European Journal of Operational Research, 232,

[Wilson, 2008]
Wilson PW (2008).
FEAR: A Software Package for Frontier Efficiency
Analysis With R.
Socio-economic Planning Sciences, 42, 247–254.

[Yagi et al., 2020]
Yagi D, Chen Y, Johnson AL, Kuosmanen T (2020).
Shape-Constrained Kernel-Weighted Least Squares: Estimating
Production Functions for Chilean Manufacturing Industries.
Journal of Business and Economic Statistics, 38,

§ LIST OF ACRONYMS

COLS: Corrected Ordinary Least Squares
Corrected Convex Nonparametric Least Squares
CER: Convex Expectile Regression
CNLS: Convex Nonparametric Least Squares
CNLS-G: Convex Nonparametric Least Squares Generic Algorithm
CQR: Convex Quantile Regression
CRS: Constant Returns to Scale
DDF: Directional Distance Function
DEA: Data Envelopment Analysis
DMU: Decision-Making Unit
FDH: Free Disposal Hull

MOM: Method of Moments

QLE: Quasi-likelihood Estimation

ICNLS: Isotonic Convex Nonparametric Least Square
ICER: Isotonic Convex Quantile Regression
ICQR: Isotonic Convex Expectile Regression
SFA: Stochastic Frontier Analysis
StoNED: Stochastic Nonparametric Envelopment of Data
StoNEZD: Stochastic Semi-nonparametric Envelopment of Z Variables Data
KDE: Kernel Density Estimation
VRS: Variable Returns to Scale

§ ESTIMATED CNLS RESIDUALS: AN ADDITIVE VRS MODEL (GAMS)

—-    160 VARIABLE e.L  error terms
1    -2.802,    2     1.414,    3   -22.224,    4  -350.911,    5   -13.700
6   101.015,    7   -28.872,    8   -14.040,    9    -0.847,    10   56.894
11  285.507,    12  679.384,    13  -20.230,    14  -70.008,    15   10.507
16   74.597,    17   -6.539,    18  -30.076,    19  -40.135,    20  -27.778
21   48.760,    22   87.008,    23   22.481,    24   23.951,    25   -2.142
26  -22.837,    27  -37.861,    28 -351.071,    29   66.858,    30  -21.703
31   -8.951,    32  216.006,    33   37.597,    34  -31.811,    35  -23.785
36 -201.341,    37  -69.290,    38    6.269,    39    3.756,    40  -74.607
41   -1.975,    42  -11.971,    43 -236.746,    44    6.575,    45   11.657
46  -20.004,    47  -67.215,    48   -5.239,    49  -55.947,    50  265.514
51    1.924,    52  -17.651,    53    4.656,    54   38.854,    55   -2.903
56  349.528,    57  163.991,    58   35.001,    59   28.393,    60  -99.132
61   32.789,    62 -604.040,    63  197.101,    64   21.411,    65  -26.755
66  -15.130,    67   29.384,    68 -128.271,    69  -16.445,    70  -13.386
71 -293.862,    72   -2.551,    73  476.620,    74   -9.804,    75   32.090
76  -17.241,    77  114.768,    78   -3.852,    79   35.396,    80  -62.937
81    8.919,    82  349.397,    83  -80.700,    84 -604.758,    85  -59.840
86    6.718,    87    6.174,    88   -6.061,    89   -0.885

§ ESTIMATED CNLS RESIDUALS: AN ADDITIVE VRS MODEL (PYSTONED)

epsilon : residual
Size=89, Index=I
Key : Lower : Value               : Upper : Fixed : Stale : Domain
0 :  None : -2.8024040090178914 :  None : False : False :  Reals
1 :  None :  1.4140528128759229 :  None : False : False :  Reals
2 :  None : -22.223744892648355 :  None : False : False :  Reals
3 :  None :    -350.91097508901 :  None : False : False :  Reals
4 :  None : -13.699899937048826 :  None : False : False :  Reals
5 :  None :  101.01484316190096 :  None : False : False :  Reals
6 :  None : -28.872348336757398 :  None : False : False :  Reals
7 :  None : -14.039749134653078 :  None : False : False :  Reals
8 :  None : -0.8474521090027167 :  None : False : False :  Reals
9 :  None :   56.89430545783776 :  None : False : False :  Reals
10 :  None :   285.5068982735597 :  None : False : False :  Reals
11 :  None :   679.3838738157001 :  None : False : False :  Reals
12 :  None : -20.230026149875016 :  None : False : False :  Reals
13 :  None :  -70.00793999792549 :  None : False : False :  Reals
14 :  None :  10.506530291158015 :  None : False : False :  Reals
15 :  None :   74.59733583012911 :  None : False : False :  Reals
16 :  None :  -6.538916384385139 :  None : False : False :  Reals
17 :  None : -30.076412673918952 :  None : False : False :  Reals
18 :  None :  -40.13473761433872 :  None : False : False :  Reals
19 :  None : -27.777867439090812 :  None : False : False :  Reals
20 :  None :   48.75945800794108 :  None : False : False :  Reals
21 :  None :   87.00800578768963 :  None : False : False :  Reals
22 :  None :  22.480756177157602 :  None : False : False :  Reals
23 :  None :    23.9507283965699 :  None : False : False :  Reals
24 :  None : -2.1416654007527427 :  None : False : False :  Reals
25 :  None : -22.836604777669237 :  None : False : False :  Reals
26 :  None : -37.861048214358334 :  None : False : False :  Reals
27 :  None :  -351.0709400239564 :  None : False : False :  Reals
28 :  None :    66.8575905653208 :  None : False : False :  Reals
29 :  None : -21.702814864843447 :  None : False : False :  Reals
30 :  None :  -8.951285786354433 :  None : False : False :  Reals
31 :  None :   216.0059312936645 :  None : False : False :  Reals
32 :  None :  37.597340051237836 :  None : False : False :  Reals
33 :  None : -31.811496011595466 :  None : False : False :  Reals
34 :  None : -23.784708760965316 :  None : False : False :  Reals
35 :  None : -201.34125864236086 :  None : False : False :  Reals
36 :  None :  -69.29033394425426 :  None : False : False :  Reals
37 :  None :  6.2687614525916615 :  None : False : False :  Reals
38 :  None :   3.755841883998798 :  None : False : False :  Reals
39 :  None :  -74.60739790684079 :  None : False : False :  Reals
40 :  None :  -1.974792422660812 :  None : False : False :  Reals
41 :  None : -11.970877239816076 :  None : False : False :  Reals
42 :  None :  -236.7462774012531 :  None : False : False :  Reals
43 :  None :   6.574582643949043 :  None : False : False :  Reals
44 :  None :  11.656881908592197 :  None : False : False :  Reals
45 :  None : -20.003997143376267 :  None : False : False :  Reals
46 :  None :  -67.21455740523754 :  None : False : False :  Reals
47 :  None :    -5.2391132374253 :  None : False : False :  Reals
48 :  None :  -55.94694781186843 :  None : False : False :  Reals
49 :  None :  265.51412810307147 :  None : False : False :  Reals
50 :  None :  1.9241988681286273 :  None : False : False :  Reals
51 :  None : -17.650818142520706 :  None : False : False :  Reals
52 :  None :   4.656303515995944 :  None : False : False :  Reals
53 :  None :   38.85434607338044 :  None : False : False :  Reals
54 :  None : -2.9033063476584573 :  None : False : False :  Reals
55 :  None :   349.5282971070195 :  None : False : False :  Reals
56 :  None :  163.99050232155884 :  None : False : False :  Reals
57 :  None :   35.00136220662897 :  None : False : False :  Reals
58 :  None :  28.392936012634948 :  None : False : False :  Reals
59 :  None :  -99.13217033642508 :  None : False : False :  Reals
60 :  None :   32.78887118253597 :  None : False : False :  Reals
61 :  None :   -604.039828622092 :  None : False : False :  Reals
62 :  None :  197.10103085237904 :  None : False : False :  Reals
63 :  None :   21.41105729756117 :  None : False : False :  Reals
64 :  None : -26.755385499560376 :  None : False : False :  Reals
65 :  None : -15.129727585118061 :  None : False : False :  Reals
66 :  None :  29.384158404973107 :  None : False : False :  Reals
67 :  None : -128.27098971868338 :  None : False : False :  Reals
68 :  None :  -16.44516975139834 :  None : False : False :  Reals
69 :  None : -13.386282784154389 :  None : False : False :  Reals
70 :  None :  -293.8615590996965 :  None : False : False :  Reals
71 :  None : -2.5506262260430503 :  None : False : False :  Reals
72 :  None :  476.61961342030463 :  None : False : False :  Reals
73 :  None :  -9.803554843610236 :  None : False : False :  Reals
74 :  None :  32.089472403185525 :  None : False : False :  Reals
75 :  None :  -17.24135562459952 :  None : False : False :  Reals
76 :  None :  114.76795801536423 :  None : False : False :  Reals
77 :  None :  -3.851436745667538 :  None : False : False :  Reals
78 :  None :   35.39622164193156 :  None : False : False :  Reals
79 :  None :  -62.93737365034673 :  None : False : False :  Reals
80 :  None :    8.91934740692409 :  None : False : False :  Reals
81 :  None :   349.3972830690898 :  None : False : False :  Reals
82 :  None :  -80.69964425322985 :  None : False : False :  Reals
83 :  None :   -604.758144376533 :  None : False : False :  Reals
84 :  None :  -59.83972941268223 :  None : False : False :  Reals
85 :  None :  6.7184004638272405 :  None : False : False :  Reals
86 :  None :   6.173694361311188 :  None : False : False :  Reals
87 :  None :  -6.061039431942607 :  None : False : False :  Reals
88 :  None :  -0.885163858119256 :  None : False : False :  Reals