abbr

# Variable selection in convex quantile regression: $\mathcal{L}_{1}$-norm or
$\mathcal{L}_{0}$-norm regularization?

Sheng Dai 111 Department of Information and Service Management, Aalto
University School of Business, 02150 Espoo, Finland. Email:
<EMAIL_ADDRESS>
(July 2021)

###### Abstract

The curse of dimensionality is a recognized challenge in nonparametric
estimation. This paper develops a new $\mathcal{L}_{0}$-norm regularization
approach to the convex quantile and expectile regressions for subset variable
selection. We show how to use mixed integer programming to solve the proposed
$\mathcal{L}_{0}$-norm regularization approach in practice and build a link to
the commonly used $\mathcal{L}_{1}$-norm regularization approach. A Monte
Carlo study is performed to compare the finite sample performances of the
proposed $\mathcal{L}_{0}$-penalized convex quantile and expectile regression
approaches with the $\mathcal{L}_{1}$-norm regularization approaches. The
proposed approach is further applied to benchmark the sustainable development
performance of the OECD countries and empirically analyze the accuracy in the
dimensionality reduction of variables. The results from the simulation and
application illustrate that the proposed $\mathcal{L}_{0}$-norm regularization
approach can more effectively address the curse of dimensionality than the
$\mathcal{L}_{1}$-norm regularization approach in multidimensional spaces.
Keywords: Variable selection, Convex quantile regression, Regularization, SDG
evaluation

##  1 Introduction

The curse of dimensionality, a well-known problem in statistics and
econometrics, refers to the poor performance of nonparametric methods in a
high-dimensional data space (Lavergne2008, *Lavergne2008). The nonparametric
estimators are weakened in prediction accuracy and exploratory power as the
dimension increases due to the sparsity of data (Stone2007, *Stone2007). This
problem has motivated a wide range of literature over the past decades
addressing different dimensionality-reduction methods, which can be utilized
to reduce the effects of dimensionality (e.g., Sinha2015, *Sinha2015;
Wilson2018, *Wilson2018; Lee2020, *Lee2020; and references therein). However,
in the production efficiency analysis area, the curse of dimensionality
remains an unresolved issue for nonparametric estimators such as data
envelopment analysis (DEA), free disposal hull (FDH), and sign-constrained
convex nonparametric least squares (SCNLS) (Lee2020, *Lee2020).

These nonparametric estimators would lead to an overfitting problem when faced
with a limited data sample size in multidimensional spaces, especially when
the estimated production function approaches the boundary of the convex hull
of the production set (Mazumder2019, *Mazumder2019). Accordingly, observations
are likely to lie close to or on the production frontier and become more
efficient, presumably jeopardizing the accuracy of production function
estimation. To date, no prominent method has been made available to directly
circumvent these problems caused by the curse of dimensionality.

To address the curse of dimensionality in production models, we could consider
either increasing sample sizes or reducing dimensions. Increasing sample sizes
aims to resatisfy a “rules of thumb” relation among the number of observations
($n$) and the number of inputs ($d$) and outputs ($q$), e.g., $n\geq 2(d+q)$
proposed in Homburg2001 and $n\geq 2d*q$ proposed in Dyson2001. However, we
cannot always enlarge sample sizes in actual applications. Moreover,
increasing sample sizes could cause a large computational burden (see, e.g.,
Dula2011, *Dula2011; Lee2013, *Lee2013). For instance, for the evaluation of
the Sustainable Development Goals (SDGs), even when using classical production
models, the input-output variables have substantial multicollinearity (e.g.,
Lamichhane2020, *Lamichhane2020; Dai2020, *Dai2020). It is not applicable to
reduce multicollinearity by increasing the number of observations. Thus,
increasing sample sizes is not an optimal choice.

By contrast, dimensionality reduction, which includes feature extraction
(i.e., variable extraction) and feature selection (i.e., variable selection)
techniques, is a more appealing approach. A vast body of literature on feature
extraction for DEA has emerged. One typical choice is to use principal
component analysis (PCA) to reduce the dimensionality of DEA production models
(e.g., Adler2010a, *Adler2010a; Nataraja2011, *Nataraja2011; Benitez-Pena2019,
*Benitez-Pena2019). Another commonly used method is to utilize the
multicollinearity among inputs or outputs (Wilson2018, *Wilson2018). However,
using the new artificial variables created by the variable extraction
technique instead of the original variables cannot lead to the formation of a
meaningful DEA-estimated frontier (Lee2020, *Lee2020). Furthermore, it is
difficult to identify the production transformation process from the new
artificial inputs to the output. To better understand the production
transformation process, variable extraction is obviously not a main technique
to reduce the dimensionality in production models.

In addition to variable extraction, variable selection is commonly seen in the
context of regression. According to a review on variable selection methods by
Sinha2015, the mainstream approaches in the statistical modeling and machine
learning fields include regularization, stepwise selection, best subset,
heuristic algorithms, and Bayesian model averaging. The pros and cons of each
approach are well addressed in Sinha2015 in detail. More importantly, most of
these approaches have been widely applied in production efficiency analysis to
handle dimensionality reduction and variable selection.

Stepwise selection (i.e., forward selection and backward elimination) is most
frequently used as a variable selection technique in production models,
implementing an ex post analysis of the sensitivity of the production frontier
to additional variables and detecting whether the relevant variables should be
included or removed (e.g., Pastor2002, *Pastor2002; Wagner2007, *Wagner2007).
However, the relevance among variables in stepwise selection will affect the
final efficiency of the frontier estimation (Li2017, *Li2017). The least
absolute shrinkage and selection operator (Lasso)-type regularization (i.e.,
$\mathcal{L}_{1}$-norm regularization) is another common technique for
choosing a subset of variables (e.g., Qin2014, *Qin2014; Lee2020, *Lee2020;
Chen2021, *Chen2021). A few recent attempts have been made to apply integer
programming (i.e., a variant of $\mathcal{L}_{0}$-norm regularization) to
achieve direct control over the number of subset variables instead of seeking
sparsity in $\mathcal{L}_{1}$-norm regularization (Keshvari2018,
*Keshvari2018; Benitez-Pena2019, *Benitez-Pena2019). However, the existing
literature on production models does not address the question of which
regularization would be a better subset selection approach via, e.g., an
effectiveness comparison. Therefore, the effectiveness comparison between the
$\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm regularization approaches in
the dimensionality reduction of variables would warrant further research.

This study is also motivated by the application to SDG evaluation. The 2030
Agenda for Sustainable Development, adopted by all 193 United Nations member
countries in 2015, is a widely acknowledged urgent call for actions with the
sustainable aims to end poverty, reduce inequality, boost economic growth and
tackle climate change.111United Nations, The 17 Goals. See more details
regarding the SDGs at https://sdgs.un.org/goals.  The core of Agenda 2030
comprises the 17 SDGs and 169 targets, representing a global consensus. In
addition, under the global indicator framework,222 United Nations, Global
indicator framework for the Sustainable Development Goals and targets of the
2030 Agenda for Sustainable Development. Available at
https://unstats.un.org/sdgs/indicators/indicators-list/.  231 unique
indicators have been designed to measure these targets and monitor progress
toward the achievement of the SDGs. Consequently, the integrated assessment
approaches based on this high number of indicators are immensely challenging
to implement due to the correlation among the indicators and the unit of
measurement differences (Lamichhane2020, *Lamichhane2020). The need and
challenge to better benchmark the degree of sustainable development of
countries triggered the emerging integration of SDG assessment and production
economics (e.g., Singpai, *Singpai).

Furthermore, the heterogeneity regarding the SDG assessment has become another
pervasive concern throughout the production function estimation process. Very
large differences can be noted in the level, the stage, and the goal of
development worldwide. It is thus improper to evaluate the degree of
sustainable development of countries following a uniform standard. However,
convex quantile regression (CQR) (Wang2014c, *Wang2014c; Kuosmanen2015d,
*Kuosmanen2015d) is an alternative that can be used to overcome heterogeneity
and that is more robust to outliers, providing an overall picture of the
conditional distributions at any given quantile. In contrast to the
conventional full frontier estimation approach, CQR can evaluate the actual
level of sustainable development locally in the interior of the production
possibility set.

In this paper, we have two main purposes: 1) We aim to compare the
$\mathcal{L}_{1}$-norm regularization with $\mathcal{L}_{0}$-norm
regularization and determine which would outperform in subset selection in the
context of the production model. In doing so, a Monte Carlo (MC) study and an
empirical application are performed to compare the effectiveness of methods in
reducing dimensionality. 2) We aim to evaluate the progress of SDGs for OECD
countries and to determine what causes the inequality in sustainable
development. Hence, extensions to convex quantile and expectile regression
using the $\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm regularization have
been developed and applied. The main contributions are thus summarized below.

  1. 1)

$\mathcal{L}_{0}$-penalized CQR. To reduce the dimensionality of the
production model and increase the robustness to heterogeneity in SDG
assessment, we develop a new $\mathcal{L}_{0}$-norm regularization approach to
convex quantile and expectile regression (Wang2014c, *Wang2014c;
Kuosmanen2015d, *Kuosmanen2015d), show how to solve it in practice using the
mixed integer programming (MIP) method, and build a link with the commonly
used $\mathcal{L}_{1}$-norm regularization approach. By taking advantage of
$\mathcal{L}_{0}$-norm regularization and quantile estimation, we can make
full use of each observation and fairly evaluate the degree of SDGs for OECD
countries even in a small sample size.

  2. 2)

Testing the effectiveness of $\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm
regularization in dimensionality reduction within the context of the
production model. $\mathcal{L}_{0}$-norm regularization directly controls the
number of subset variables rather than the level of sparsity uniformly. Does
this mean that $\mathcal{L}_{0}$-norm regularization is a better choice for
subset selection? It is thus natural to compare the performance of the
proposed $\mathcal{L}_{0}$-norm regularization with that of the competing
$\mathcal{L}_{1}$-norm regularization. In the MC and application sections, we
show that the proposed $\mathcal{L}_{0}$-norm regularization outperforms in
terms of dimensionality reduction and that $\mathcal{L}_{1}$-norm
regularization cannot implicitly eliminate the input variables in the
regression-based production function estimation.

  3. 3)

Application to SDG assessment in OECD countries. There are few SDG evaluations
in the context of production economics. We present an empirical assessment in
Section 4 from the perspective of production function estimation using the
OECD country data from the SDG indicators database.

Furthermore, the penalized convex quantile and expectile regression approaches
would be greatly subject to the computational burden from the O($n^{2}$)
linear constraints as the number of observations increases. To mitigate the
computational burden, we improve upon the base of the cutting-plane algorithm
(Bertsimas2020, *Bertsimas2020) to solve the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-penalized convex quantile and expectile regression
approaches. An additional experiment is performed to justify the superiority
in performance of the new algorithm over the CNLS-G algorithm proposed by
Lee2013.

The rest of the paper is organized as follows. Section 2 introduces the
$\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-penalized convex quantile and
expectile regression approaches. Section 3 performs an MC study to compare the
finite sample performances of the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-norm regularization approaches from two metrics. An
empirical application to SGD assessment is prepared in Section 4. Section 5
concludes this paper. An introduction to the CNLS-A algorithm and additional
tables and figures are provided in the Appendix A.

##  2 Methodology

###  2.1 Penalized quantile production function

Consider the nonparametric conditional quantile production functions (e.g.,
Dai2020, *Dai2020),

$Q_{y_{i}}(\tau\,|\,{\boldsymbol{x}}_{i})=f_{\tau}({\boldsymbol{x}}_{i})+F^{-1}_{\varepsilon_{i}}(\tau),\quad\text{for}\quad
i=1,\cdots,n$ (1)

where $\tau\in(0,1)$ denotes the quantile order, $y_{i}\in{\mathbb{R}}$ and
${\boldsymbol{x}}_{i}\in{\mathbb{R}}^{d}$ are the output and input variables,
respectively, $\varepsilon_{i}$ denotes the independent and identically
distributed errors, and $F_{\varepsilon_{i}}$ represents the cumulative
distribution function of the errors.

The nonparametric quantile production model (1) does not require the a priori
functional form for
$f_{\tau}(\cdot):{\mathbb{R}}^{d}\xrightarrow{}{\mathbb{R}}$ but instead
assumes that $f_{\tau}(\cdot)$ satisfies certain axiomatic properties (e.g.,
monotonicity, concavity). As such, one can readily resort to CQR (Wang2014c,
*Wang2014c; Kuosmanen2015d, *Kuosmanen2015d) or the more appealing convex
expectile regression (CER) (Kuosmanen2020, *Kuosmanen2020; Kuosmanen2020b,
*Kuosmanen2020b) to estimate the quantile production functions.

To reduce the dimensions of input variables in the production model, as in
Lee2020, we introduced the $\mathcal{L}_{1}$-norm regularization, known as
Lasso (Tibshirani1996, *Tibshirani1996), for the nonparametric quantile
production function.

$\hat{Q}(\tau\,|\,{\boldsymbol{x}}_{i})=\operatorname*{arg\,min}_{f_{\tau}\in{\mathcal{F}}}\sum^{n}_{i=1}\rho_{\tau}(y_{i}-f_{\tau}({\boldsymbol{x}}_{i}))+\lambda\left\|P(f_{\tau})\right\|_{1}$
(2)

where $\rho_{\tau}(t)=(\tau-1\\{t\leq 0\\})t$ is the check function
(Koenker1978, *Koenker1978), $\lambda\geq 0$ is the tuning parameter, and
$\left\|P(f_{\tau})\right\|_{1}$ denotes the Lasso regularization term.

To date, although the Lasso regularization has been predominant in the
literature because of its computational advantage, several recent studies have
pointed out that Lasso regularization is more of a variable screening than a
variable selection procedure (e.g., Fan2010, *Fan2010; Su2017, *Su2017;
Bertsimas2017, *Bertsimas2017). To effectively select the subset, another
regularization $\mathcal{L}_{0}$ has been applied and proven to outperform the
Lasso regularization, especially in a high signal-to-noise (SNR) regime (see,
e.g., Bertsimas2016, *Bertsimas2016; Hastie2020, *Hastie2020). Compared with
Lasso regularization, which uniformly shrinks the variables,
$\mathcal{L}_{0}$-norm regularization sparsifies the variables only with the
desired shrinking (Bertsimas2020b, *Bertsimas2020b).

We therefore further considered the $\mathcal{L}_{0}$-norm regularization for
the nonparametric quantile production function. The $\mathcal{L}_{0}$-norm
regularization formulation

$\displaystyle\hat{Q}(\tau\,|\,{\boldsymbol{x}}_{i})$
$\displaystyle=\operatorname*{arg\,min}_{f_{\tau}\in{\mathcal{F}}}\sum^{n}_{i=1}\rho_{\tau}(y_{i}-f_{\tau}({\boldsymbol{x}}_{i}))$
(3) s.t. $\displaystyle\left\|P(f_{\tau})\right\|_{0}\leq k$

where $k$ represents the subset size and $\left\|P(f_{\tau})\right\|_{0}$
denotes the $\mathcal{L}_{0}$-pseudonorm regularization term. Alternatively,
one could apply a Lagrangian form similar to that of an Lasso regularization,
that is,
$\hat{Q}(\tau\,|\,{\boldsymbol{x}}_{i})=\operatorname*{arg\,min}_{f_{\tau}\in{\mathcal{F}}}\sum^{n}_{i=1}\rho_{\tau}(y_{i}-f_{\tau}({\boldsymbol{x}}_{i}))+\lambda\left\|P(f_{\tau})\right\|_{0}$.
In this paper, we used an explicit formulation to ensure that the size of the
support remains bounded to decrease the computational cost. To estimate the
penalized quantile production models (2) and (3), we propose two new penalized
CQR approaches by imposing the $\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm
regularization procedures and elaborate on both methods in the following
sections.

###  2.2 $\mathcal{L}_{1}$-penalized convex quantile regression

While the CQR estimates the full set of variables, the
$\mathcal{L}_{1}$-penalized CQR ($\mathcal{L}_{1}$-CQR) aims to retain the
true variables and estimate the subset selection by shrinking some
coefficients and setting the others to 0. Given the prespecified tuning
parameter $\lambda$, we extended the CQR problem to estimate the
$\mathcal{L}_{1}$-CQR as follows:

$\displaystyle\underset{\mathbf{\alpha},\mathbf{{\boldsymbol{\beta}}},{{\mathbf{\varepsilon}}^{\text{+}}},{{\mathbf{\varepsilon}}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}+\lambda\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}$
(4) s.t. $\displaystyle
y_{i}=\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon^{+}_{i}-\varepsilon^{-}_{i}$
$\displaystyle\forall i$
$\displaystyle\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{{{\boldsymbol{x}}}_{i}}\leq\mathbf{\alpha}_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

The $\mathcal{L}_{1}$-CQR problem differs from the CQR problem in that the
objective function adds an extra $\mathcal{L}_{1}$-norm regularization to
$\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}$, where
$\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}=\sum^{d}_{j}\sum^{n}_{i}\left|\beta_{j,i}\right|$.
Similar to the Lasso regression, we can recover the original CQR problem as
$\lambda\xrightarrow{}0$ and remove the corresponding input variables with
$\beta_{j,i}=0$ for all observations as $\lambda\xrightarrow{}\infty$. The
first constraint in Problem (4) can be interpreted as a multivariate
regression equation. The second constraint, i.e., the system of Afriat
inequality, imposes the convexity. The third constraint imposes the
monotonicity, and the last is the sign constraint of the error terms.

Similar to the connection between CQR and CER problems, the
$\mathcal{L}_{1}$-penalized CER ($\mathcal{L}_{1}$-CER) approach is a more
compelling alternative to the $\mathcal{L}_{1}$-CQR problem due to the
nonuniqueness of linear programming problem (4). Compared with the
$\mathcal{L}_{1}$-CQR problem, the $\mathcal{L}_{1}$-CER problem uses an
alternative quadratic objection function to ensure a unique solution to the
following problem.

$\displaystyle\underset{\mathbf{\alpha},\mathbf{{\boldsymbol{\beta}}},{{\mathbf{\varepsilon}}^{\text{+}}},{\mathbf{\varepsilon}}^{-}}{\mathop{\min}}$
$\displaystyle\,\tilde{\tau}\sum\limits_{i=1}^{n}(\varepsilon_{i}^{+})^{2}+(1-\tilde{\tau})\sum\limits_{i=1}^{n}(\varepsilon_{i}^{-})^{2}+\lambda\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}$
(5)

where the $\mathcal{L}_{1}$-CER problem is subject to the same constraints as
those of the $\mathcal{L}_{1}$-CQR problem. The expectile $\tilde{\tau}$ in
(5) and quantile $\tau$ in (4) can be converted to each other; see more
discussions in Kuosmanen2020. In practice, the $\mathcal{L}_{1}$-CQR problem
(4) and $\mathcal{L}_{1}$-CER (5) problem can be solved directly by commercial
off-the-shelf solvers (e.g., Cplex, Mosek, and Gurobi) or standalone
algorithms (e.g., Mazumder2019, *Mazumder2019; Lin2020 *Lin2020).

We notice that Lee2020 propose a similar Lasso-penalized DEA model (i.e.,
Lasso-SCNLS), which is, in essence, a special variant of the
$\mathcal{L}_{1}$-CER model, i.e., expectile $\tilde{\tau}=0.5$. However, in
their study, it is unclear whether this Lasso-style penalty would necessarily
reduce the dimensionality (set all $\beta_{k,i}=0$ for some $k$) (Chen2021,
*Chen2021). Although Lee2020 present some simulations that suggest improvement
in the mean squared error (MSE) compared with that obtained from other DEA-
related approaches, they fail to justify with convincing evidence whether this
Lasso-style penalty effectively helps to reduce dimensionality. However,
according to this study, the Lasso-SCNLS approach could eliminate the true
irrelevant variables to some extent. Yet, the complete evidences (e.g.,
prediction error metric or accuracy metric) to support this argument would
warrant further research.

Note further that the quadratic objective function in the
$\mathcal{L}_{1}$-CER problem can ensure a unique estimated quantile function
($\hat{Q}$) but cannot guarantee the uniqueness of the subgradients (i.e., the
estimated coefficients $\hat{\beta}_{i,j}$). To obtain a unique minimizer,
denoted by ($\hat{Q}$, $\hat{\beta}_{i,j}$) in Problem (4) or (5), one could
apply an alternative $\mathcal{L}_{1}$-norm squared regularization to the
subgradients for the original CQR problem. Note that in this case, the
convexity of the objective function becomes stronger, and both estimated
quantile functions and estimated subgradients are then unique, a detailed
proof of which would be an interesting topic for future research.

###  2.3 $\mathcal{L}_{0}$-penalized convex quantile regression

Considering the potential advantages of $\mathcal{L}_{0}$-norm regularization
in subset selection, we extended the $\mathcal{L}_{1}$-CQR problem to the
cardinality penalized CQR problem. The proposed $\mathcal{L}_{0}$-penalized
CQR approach ($\mathcal{L}_{0}$-CQR) is defined as

$\displaystyle\underset{\mathbf{\alpha},\mathbf{{\boldsymbol{\beta}}},{\mathbf{\varepsilon}^{\text{+}}},{\mathbf{\varepsilon}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(6) s.t. $\displaystyle
y_{i}=\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon^{+}_{i}-\varepsilon^{-}_{i}$
$\displaystyle\forall i$
$\displaystyle\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{{{\boldsymbol{x}}}_{i}}\leq\mathbf{\alpha}_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$
$\displaystyle\left\|{\boldsymbol{\beta}}_{i}\right\|_{0}\leq k$
$\displaystyle\forall i$

where $\left\|{\boldsymbol{\beta}}_{i}\right\|_{0}:=\left|j:\beta_{j,i}\neq
0\right|,\forall j\in[d]$, denotes the $\mathcal{L}_{0}$ pseudonorm,
indicating the number of nonzero coefficients of ${\boldsymbol{\beta}}_{i}$.
Accordingly, the last constraint in Problem (6) limits the number of input
variables to at most $k$. If $k\geq d$, the optimal solution to the
$\mathcal{L}_{0}$-CQR problem (6) is equivalent to the optimal solution to the
original CQR problem because the cardinality constraint is not binding; if
$k<d$, there are exactly $k$ input variables that contribute to the estimated
quantile function.

Although Problem (6) is theoretically appealing, it is an NP-hard optimization
problem (Natarajan1995, *Natarajan1995; Bertsimas2016, *Bertsimas2016). To
solve Problem (6) in practice, we reformulated it into a mixed integer linear
programming (MILP) problem. We defined the binary variable ${\boldsymbol{z}}$
to determine what the hyperplanes denoted by the pair of coefficients
($\alpha_{i}$, ${\boldsymbol{\beta}}_{i}$) look like. Given an a priori-
specified $k$, the MILP formulation for the $\mathcal{L}_{0}$-CQR problem is

$\displaystyle\underset{\mathbf{\alpha},\mathbf{{\boldsymbol{\beta}}},{{\mathbf{\varepsilon}}^{\text{+}}},{\mathbf{\varepsilon}}^{-},{\boldsymbol{z}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(7) s.t. $\displaystyle
y_{i}=\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon^{+}_{i}-\varepsilon^{-}_{i}$
$\displaystyle\forall i$
$\displaystyle\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{{{\boldsymbol{x}}}_{i}}\leq\mathbf{\alpha}_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$
$\displaystyle\left|{\boldsymbol{\beta}}_{i}\right|\leq M{\boldsymbol{z}}$
$\displaystyle\forall i$ $\displaystyle{\boldsymbol{z}}\in\\{0,1\\}^{d}$
$\displaystyle\sum^{d}_{j=1}z_{j}\leq k$

where $M$ is a predetermined positive number that can be either estimated from
the variants of the CQR problem (c.f., Bertsimas2020, *Bertsimas2020) or
enumerated from a range of candidates. Note that a smaller $M$ could make the
solution too constrained and that a larger $M$ may lead to overfitting. For
the $\mathcal{L}_{0}$-CER problem, we can replace the objective function of
(7) with the quadratic objective function
$\tilde{\tau}\sum\limits_{i=1}^{n}(\varepsilon_{i}^{+})^{2}+(1-\tilde{\tau})\sum\limits_{i=1}^{n}(\varepsilon_{i}^{-})^{2}$
to ensure the uniqueness, which is now a mixed integer quadratic programming
(MIQP) problem. We then employed a mixed integer optimization solver (e.g.,
Cplex and Gurobi) to solve the $\mathcal{L}_{0}$-CQR/CER problem.

######  Theorem 1.

Problems (6) and (7) are equivalent.

###### Proof.

Let ${\boldsymbol{\beta}}_{i}^{*}$ be the optimal solution to Problem (6) that
satisfies $\left\|{\boldsymbol{\beta}}_{i}^{*}\right\|_{\infty}\leq M$. After
introducing the binary variable ${\boldsymbol{z}}$, we can obtain a linear
constraint on ${\boldsymbol{z}}$ by imposing the cardinality constraint on
${\boldsymbol{\beta}}_{i}$,
${\boldsymbol{z}}:=\\{z\in\\{0,1\\}^{d}:{\boldsymbol{e}}^{{}^{\prime}}z\leq
k\\}$, where all components in vector ${\boldsymbol{e}}$ are equal to one.
Furthermore,
$M{\boldsymbol{z}}=M\sum^{d}_{j}z_{j}=M\sum^{d}_{j:z_{j}=1}z_{j}=M{\boldsymbol{z}}_{z_{j}}$.
Thus, $\left|{\boldsymbol{\beta}}_{i}\right|\leq M$ if $z_{j}=1$; otherwise,
${\boldsymbol{\beta}}_{i}=0$. The last constraint of Problem (6) can be
equivalently expressed as the last three constraints of Problem (7). ∎

Apart from the computational convenience, the explicit constraint formulation
(7) provides a link to the $\mathcal{L}_{1}$-CQR problem (4). We described
here the relation between the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
problems. Let Conv($A$) denote the convex hull of the set of $A$. Similarly to
Kuosmanen2008, we set the $\mathcal{L}_{0}$-CQR problem to the function
$\hat{Q}:{\mathbb{R}}^{d}\xrightarrow{}{\mathbb{R}}$:

$\displaystyle\hat{Q}(\tau\,|\,{\boldsymbol{x}}):=\sup\Big{\\{}\sum^{n}_{i=1}a_{i}\hat{y}_{\tau,i}:\sum^{n}_{i=1}a_{i}=1,\sum^{n}_{i=1}a_{i}{\boldsymbol{x}}_{i}={\boldsymbol{x}},a_{i}\geq
0,$ $\displaystyle\left|{\boldsymbol{\beta}}_{i}\right|\leq
M{\boldsymbol{z}},{\boldsymbol{z}}\in\\{0,1\\},\sum^{d}_{j=1}z_{j}\leq
k\Big{\\}}$ (8)

where $a_{i}$ is the weight assigned to observation $i$ and
$\sup(\emptyset)=-\infty$. $\hat{Q}$ is thus well defined and finite on
Conv(${\boldsymbol{x}}$). Using these settings, Theorem 2 shows that the
optimum objective value of the $\mathcal{L}_{0}$-CQR problem is low-bounded by
that of the $\mathcal{L}_{1}$-CQR problem.

######  Theorem 2.

The optimal solution to Problem (4) is greater than or equal to that of
Problem (7).

###### Proof.

Similar to Bertsimas2016, we first reformulated Conv(${\boldsymbol{x}}$) (8)
into a relaxation form:

$\displaystyle\hat{Q}_{1}(\tau\,|\,{\boldsymbol{x}}):$
$\displaystyle=\sup\Big{\\{}\sum^{n}_{i=1}a_{i}\hat{y}_{\tau,i}:\sum^{n}_{i=1}a_{i}=1,\sum^{n}_{i=1}a_{i}{\boldsymbol{x}}_{i}={\boldsymbol{x}},a_{i}\geq
0,\left|{\boldsymbol{\beta}}_{i}\right|\leq
M{\boldsymbol{z}},{\boldsymbol{z}}\in\\{0,1\\},\sum^{d}_{j=1}z_{j}\leq
k\Big{\\}}$
$\displaystyle=\sup\Big{\\{}\sum^{n}_{i=1}a_{i}\hat{y}_{\tau,i}:\sum^{n}_{i=1}a_{i}=1,\sum^{n}_{i=1}a_{i}{\boldsymbol{x}}_{i}={\boldsymbol{x}},a_{i}\geq
0,\left\|{\boldsymbol{\beta}}_{i}\right\|_{\infty}\leq
M,\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}\leq Mk\Big{\\}}$ (9)

We then rewrote Problem (4) in constrained form as the convex hull:

$\displaystyle\hat{Q}_{2}(\tau\,|\,{\boldsymbol{x}}):=\sup\Big{\\{}\sum^{n}_{i=1}a_{i}\hat{y}_{\tau,i}:\sum^{n}_{i=1}a_{i}=1,\sum^{n}_{i=1}a_{i}{\boldsymbol{x}}_{i}={\boldsymbol{x}},a_{i}\geq
0,\left\|{\boldsymbol{\beta}}_{i}\right\|_{1}\leq Mk\Big{\\}}$ (10)

It is straightforward to show the link between $\mathcal{L}_{1}$-CQR and
$\mathcal{L}_{0}$-CQR: $\hat{Q}_{1}\subseteq\hat{Q}_{2}$. We thus find the
inequality below: $\hat{Q}_{1}\geq\hat{Q}_{2}$. That is, Problem (7) provides
lower bounds for Problem (4). Note that we aim to find the minimum objective
function values of the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
problems. ∎

In practice, there is one concern rooted in the original CQR problem when
using the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR approaches.
Estimating a larger sample becomes excessively expensive in both the CQR and
CER approaches due to the O($n^{2}$) linear constraints (Lee2013, *Lee2013;
Mazumder2019, *Mazumder2019). For example, if the data samples contain 300
observations, the total number of linear constraints amounts to 90000.
Undoubtedly, the computational burden becomes a barrier to the application of
the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR approaches under larger
sample sizes and calls for a more efficient algorithm.

Recent works are showing several promising algorithms from the perspective of
computational burden (see, e.g., Lee2013, *Lee2013; Balazs2015,*Balazs2015;
Mazumder2019, *Mazumder2019; Bertsimas2020, *Bertsimas2020; Lin2020 *Lin2020).
It is worth highlighting two of them here, namely, the CNLS-G algorithm
proposed by Lee2013 and the cutting-plane algorithm proposed by Balazs2015 and
extended by Bertsimas2020. Although both algorithms use the relaxed Afriat
inequality constraint set and iteratively introduce new ones as necessary, the
cutting-plane algorithm is more efficient to solve the CQR and CER problems
due to its fast identification of violating constraints. We thus improved on
the base of the cutting-plane algorithm in Bertsimas2020 to solve Problems (4)
and (7) and referred to it as the CNLS-adapted cutting-plane (CNLS-A)
algorithm. In the Appendix A, we introduce the CNLS-A algorithm in detail and
perform a computational test to compare the performances of the CNLS-G and
CNLS-A algorithms in solving the CQR and CER problems.

##  3 Monte Carlo study

###  3.1 Setup

We present an MC study to compare the finite sample performances in reducing
the dimensionality between the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-penalized CQR/CER models. Following Lee2013, we considered
an additive Cobb–Douglas production function with $d$ input variables and one
output variable. Given $n$ (number of observations), $d$ (number of input
variables), $k$ (number of true input variables), $\rho$ (the SNR level), and
$\tau$ (the quantile level), our experimental process was as illustrated
below:

  1. 1)

we defined the following parameters to generate the artificial data:
$n\in\\{100,500\\}$, $d\in\\{6,8,10,12\\}$, $\rho\in\\{0.5,2,10\\}$, and
$\tau\in\\{0.1,0.3,0.5,0.7,0.9\\}$;

  2. 2)

we randomly sampled a support set ($\omega$) of size $k$ from
$\\{1,2,\cdots,d\\}$ and considered that the sizes of the true support set
($\omega^{*}$) are 2 and 4, respectively;

  3. 3)

we drew the input variables ${\boldsymbol{x}}_{i}\in{\mathbb{R}}^{n\times d}$
i.i.d. from a uniform distribution, $U[1,10]$;

  4. 4)

the observed output
$y_{i}=\prod_{d=1}^{D\in\omega^{*}}{\boldsymbol{x}}^{\frac{0.8}{d}}_{d,i}+v_{i}$.
For each noise, $v_{i}$ was randomly sampled from the normal distribution,
$N(0,\sigma^{2})$. $\sigma^{2}$ was calculated based on the SNR level,
$\sigma^{2}=var(\prod_{d=1}^{D\in\omega^{*}}{\boldsymbol{x}}^{\frac{0.8}{d}}_{d,i})/\rho$
(Bertsimas2016, *Bertsimas2016);

  5. 5)

we conducted the $\mathcal{L}_{1}$-CQR, $\mathcal{L}_{0}$-CQR,
$\mathcal{L}_{1}$-CER, and $\mathcal{L}_{0}$-CER on the data
${\boldsymbol{x}}_{i}$, $y$, with the optimal tuning parameters; for each
approach, the tuning parameter was selected by a 5-fold cross-validation
procedure from a wide range of values (Bertsimas2020, *Bertsimas2020;
chen2020c, *chen2020c; Chen2021a, *Chen2021a);

  6. 6)

we measured the finite sample performance of each approach using the following
two statistics: prediction error and accuracy;

  7. 7)

we replicated steps 1)-6) 10 times and averaged the statistics.

To assess the finite sample performance of each method, we utilized the
prediction error and accuracy statistics, which are, respectively, the level
of good data fidelity and the proportion of true variables that were selected.
For a penalized quantile function estimate, the prediction error statistic is
defined as (Bertsimas2016, *Bertsimas2016)

$\text{Prediction
error}=\left\|\hat{Q}_{i}-Q^{\omega^{*}}_{i}\right\|^{2}_{2}/\left\|Q^{\omega^{*}}_{i}\right\|^{2}_{2}$

where $\hat{Q}_{i}$ denotes the estimated conditional penalized quantile
function, identified by a particular approach, and $Q^{\omega^{*}}_{i}$
indicates the true conditional penalized quantile function. Note that for the
prediction error statistic, an “in-sample” version is used in this paper,
where the smaller the prediction error is, the better the performance will be.
Furthermore, we reported the proportion of true variables that were selected,
that is, the accuracy (%) (Bertsimas2020, *Bertsimas2020)

$\text{Accuracy}=\frac{\left|\hat{\omega}\cap\omega^{*}\right|}{k}\times
100\%$

Regarding the tuning parameter selection, by using a 5-fold cross-validation
procedure (see, e.g., Bertsimas2020, *Bertsimas2020; chen2020c, *chen2020c;
Chen2021a, *Chen2021a), we calibrated the parameter $\lambda$ over a sequence
of 100 values for the $\mathcal{L}_{1}$-CQR/CER approach. For the
$\mathcal{L}_{0}$-norm regularization, we varied the tuning value $M$ from the
multiplier set $\\{0.1,0.5,1,1.5,2,3,4,5\\}$ for the $\mathcal{L}_{0}$-CQR/CER
approach, and we chose the subset size $k$ from the set $\\{1,\cdots,d-1\\}$
(i.e., at least one input). Note that the 5-fold cross-validation technique
needs to solve each $\mathcal{L}_{1}$/$\mathcal{L}_{0}$-CQR/CER problem five
times for each candidate value; a faster procedure for finding the tuning
parameters for a similar optimization problem can be seen in chen2020c.

In all the experiments that follow, we resorted to Julia/JuMP (Dunning2017,
*Dunning2017) to solve the $\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-CQR/CER
models based on the CNLS-A algorithm. We used the commercial off-the-shelf
solver Gurobi (9.1.2) to solve the CNLS-A algorithm. All experiments were
performed on Aalto University’s high-performance computing cluster Triton with
Xeon @2.8 GHz processors, 5 CPUs, and 5 GB RAM per CPU. Since the MC study is
naturally parallel, we resorted to the multiple threads in Julia to speed up
the computation. The skeleton codes are publicly available at
https://github.com/ds2010/penalizedCQR.

###  3.2 Prediction error

We first explored the prediction error of the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-norm regularization approaches solved by the CNLS-A
algorithm. For the comparison of the quantile and expectile approaches, the
expectile $\tilde{\tau}$ was transformed into the corresponding quantile
$\tau$ based on the empirical inverse quantile functions of the noise
($\varepsilon_{i}$).

Figs. 1 and 2 depict how the prediction error scales as the sample sizes,
input dimensions, and quantiles increase when estimating the
$\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CER problems with the observation
$n=100$ but different sizes of the true support set (i.e., $k=2$ and $k=4$).
Recall that the smaller the prediction error is, the better the performance.
Combined with Figs. B1 and B2, $\mathcal{L}_{0}$-CQR shows a major advantage
over $\mathcal{L}_{1}$-CQR in quantile production function estimation when
applying the subset selection procedure. The main findings are summarized as
follows.

Figure 1: Prediction errors of the $\mathcal{L}_{1}$-CQR and
$\mathcal{L}_{0}$-CQR approaches with $n=100$ and $k_{\text{true}}=2$. Figure
2: Prediction errors of the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
approaches with $n=100$ and $k_{\text{true}}=4$.

  * •

The $\mathcal{L}_{0}$-CQR approach generally outperformed the
$\mathcal{L}_{1}$-CQR approach based on the comparison of the prediction
errors, especially in a high-SNR regime. However, the $\mathcal{L}_{1}$-CQR
performed better in certain cases in a low-SNR regime, which confirms the
conclusion by Hastie2020.

  * •

The quantile did not influence the prediction error systematically. As shown
in Figs. 1, 2, B1, and B2, there was a “U-shaped” relationship between the
prediction error and the quantile for both the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-CQR approaches. The lowest bias was observed in the median
estimation; the higher bias was observed in the lower-quantile or the higher-
quantile estimation.

  * •

The $\mathcal{L}_{0}$-CQR and $\mathcal{L}_{1}$-CQR approaches tended to yield
larger prediction errors as the number of dimensions increased, especially for
a large true support set. Lee2020 also report similar conclusions for the
Lasso-SCNLS in terms of the MSE change; that is, the MSE increases as more
inputs are included. This is because the high-dimensional data space usually
has a large sparsity, which undermines the prediction accuracy of both
approaches.

  * •

The noise term obviously affected the precision of the estimated quantile
function. As the SNR increased, the prediction error became smaller for both
the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR approaches. Furthermore,
when the size of the true support set increased, the prediction error also
decreased.

  * •

The sample size effect was also observed in the experiments. Compared to
estimating a larger sample size (Figs. B1 and B2), the results of Figs. 1 and
2 suggest that a large sample would decrease the prediction error in both
approaches. A similar conclusion is commonly seen in most of the existing
literature (see, e.g., Lee2020, *Lee2020; Chen2021, *Chen2021).

Because of the nonuniqueness of estimated objective function $\hat{Q}$ in the
penalized quantile approaches, it is unfair to compare the performances of the
penalized quantile and expectile approaches in estimating the quantile
production function. In fact, the prediction errors (which are available upon
request) show that neither the quantile approaches nor the expectile
approaches clearly outperform the others. Within this context of subset
selection, the core competitiveness of the methods relies on whether they can
accurately eliminate the irrelevant variables. The accuracy metric would thus
be a better measurement approach.

###  3.3 Accuracy

To provide direct evidence of the performance of the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-norm regularization approaches in subset variable selection,
we next evaluated the accuracy statistic by counting the number of nonzero
estimated ${\boldsymbol{\beta}}_{i}$, which is compared with the true support
set $\omega^{*}$. Recall that the higher the accuracy is, the better the
performance.

Figs. 3 and B3 plot the accuracy results with $n=500$ and $n=100$,
respectively, but yield somewhat similar results. In most cases, the
$\mathcal{L}_{0}$-CQR approach outperformed the $\mathcal{L}_{1}$-CQR
approach. As expected, it was observed that the accuracy decreased in all
scenarios with increasing sample size and that $\mathcal{L}_{0}$-CQR achieved
the highest accuracy. In certain cases (e.g., low quantile with $\tau=0.1$),
the accuracy of the $\mathcal{L}_{1}$-CQR approach was zero, indicating that
the $\mathcal{L}_{1}$-CQR approach was unable to completely eliminate the
irrelevant variables compared with the true support set. However, the
$\mathcal{L}_{1}$-CQR approach can successfully select the true variables, at
least partially. Therefore, for accuracy comparison, the $L_{0}$-CQR approach
is a better choice for subset selection.

Figure 3: Accuracy of the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
approaches with $n=500$ and $k=4$.

Furthermore, similar to the prediction error comparison, several other main
findings based on the accuracy benchmarking are summarized as follows:

  * •

Fig. 3 also shows no clear deterministic trend with decreasing or increasing
quantile. For the $\mathcal{L}_{1}$-CQR approach, the lower-quantile
estimation yielded poor accuracy and even could be zero. By contrast, the
$\mathcal{L}_{0}$-CQR approach achieved a higher accuracy at the median
quantile. Furthermore, the accuracy of $\mathcal{L}_{0}$-CQR was more than 50%
in almost all scenarios and close to 100% in large parts. Note that the
accuracy difference in each quantile estimation provides direct evidence of
the superiority of the $\mathcal{L}_{0}$-CQR approach, especially in a high-
SNR regime.

  * •

As more input variables were included, using both methods decreased the
accuracy in terms of subset variable selection. This is due to the fact that
the higher dimension makes the problem more sparse, inducing difficulty in
true variable identification and selection. Although the $\mathcal{L}_{1}$-CQR
approach achieved lower accuracy when $d=12$ with the lower SNR,
$\mathcal{L}_{1}$-CQR generally performed better in most scenarios.

  * •

The noise term also affected the accuracy of subset variable selection. The
larger noise amplitude in the data resulted in unstable accuracy estimation
and normally yielded lower accuracy results. Furthermore, considering the
results demonstrated in Figs. B3-B5, the smaller size of the true support set
made the subset selection easier for both approaches, as reflected by the low
accuracy.

We further investigated the accuracy of the penalized quantile and expectile
approaches for different input dimensions, true support set sizes, and SNR
levels. Table 1 demonstrates the simulation results with sample size $n=100$
and quantile $\tau=0.9$. As shown in Table 1, the expectile estimation
approach generally outperformed the quantile approach. That is, the expectile
approach can select the subset variable more effectively. This is because the
expectile approach has a stronger convex objective function. However, the
expectile approach performed worse compared to the quantile approach in
certain cases, possibly due to the well-known small sample bias. Since the
performances of the $\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm
regularization approaches show a nonsystematic trend as the quantile increases
or decreases, the comparison for other quantiles also presents similar results
to those in Table 1.

Table 1: Comparison of the accuracy of quantile and expectile estimation for $\mathcal{L}_{0}$-norm regularization. $d$ | $k$ | $\text{SNR}=0.5$ | $\text{SNR}=2$ | $\text{SNR}=10$
---|---|---|---|---
quantile | expectile | quantile | expectile | quantile | expectile
6 | 2 | 42 | 42 | 43 | 52 | 56 | 63
| 4 | 76 | 78 | 74 | 73 | 69 | 79
8 | 2 | 39 | 34 | 39 | 39 | 46 | 39
| 4 | 54 | 51 | 53 | 61 | 47 | 63
10 | 2 | 22 | 24 | 45 | 34 | 47 | 48
| 4 | 44 | 47 | 38 | 49 | 65 | 56
12 | 2 | 24 | 21 | 39 | 25 | 37 | 40
| 4 | 29 | 47 | 31 | 46 | 44 | 46

In conclusion, the experiments seem informative enough to indicate that
$\mathcal{L}_{0}$-CQR can outperform $\mathcal{L}_{1}$-CQR in terms of
prediction error and accuracy. This is particularly true when the
$\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-CQR associated with a large sample
size or higher-noise data are estimated. The experiments also reveal that the
expectile approaches generally perform better than the quantile approaches.
Thus, we apply the expectile-based regularization approaches
($\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-CER) to select the subset variables
and analyze the inequality in sustainable development in OECD countries.

##  4 Application to SDG evaluation

In this section, we explore the effectiveness of the proposed
$\mathcal{L}_{1}$\- and $\mathcal{L}_{0}$-norm regularization approaches in
reducing dimensionality empirically and evaluated the performance of
sustainability in 35 OECD countries using the SDG indicators.

###  4.1 Data and variables

To benchmark the degree of sustainable development for OECD countries, we
estimated the quantile production function based on the panel data of 35 OECD
countries during 2017, 2019, and 2020.333 Colombia and Lithuania are excluded
from the list of OECD countries due to insufficient data. The data for 2018
are also excluded due to the inconsistent indicators from the same database.
We collected the raw data from the UN SDG Index and Dashboard Reports
(Sachs2017, *Sachs2017, *Sachs2019, *Sachs2020), which provide a substantial
range of indicators and detailed descriptions of the indicator definition,
indicator measurement, and indicator source. The same SDG data have already
been applied in other empirical applications of sustainability evaluation
(see, e.g., Lamichhane2020, *Lamichhane2020; Huan2021, *Huan2021).

The input-output variables in this paper were selected based on three pillars
of sustainability: social, economic, and environmental. In the framework of
sustainable development, the 17 SDGs can be divided into these three aspects
but with different classifications (Dalampira2020, *Dalampira2020). For
example, environmental sustainability can be achieved through Goals 6, 13, 14,
and 15; economic sustainability through Goals 7, 8, 9, 11, and 12; and social
sustainability through Goals 1, 2, 3, 4, 5, 10, 16, and 17 (Costanza2016,
*Costanza2016). Using this criterion combined with the UN SDG Index and
Dashboard Reports, we retrieved the following variables to characterize the
production of sustainability (Table 2).

Table 2: Description of the selected indicators for the sustainability evaluation. Pillar | Goal | Indicator
---|---|---
Economic | SDG8 | ($y$) Adjusted GDP growth (%)
Social | SDG1 | ($x_{1}$) Poverty headcount ratio at US$ 1.90/day (%)
| SDG3 | ($x_{2}$) Life expectancy at birth (years)
|  | ($x_{3}$) Subjective well-being (average ladder score)
| SDG8 | ($x_{4}$) Employment-to-population ratio (%)
| SDG10 | ($x_{5}$) Gini coefficient adjusted for top income (1-100)
Environmental | SDG11 | ($x_{6}$) Annual mean concentration of PM2.5 ($\mu$g/m3)
| SDG12 | ($x_{7}$) Municipal solid waste (kg/day/capita)
|  | ($x_{8}$) E-waste generated (kg/capita)
|  | ($x_{9}$) Production-based SO2 emissions (kg/capita)
|  | ($x_{10}$) Nitrogen production footprint (kg/capita)
| SDG13 | ($x_{11}$) Energy-related CO2 emissions per capita (tCO2/capita)
|  | ($x_{12}$) Effective carbon rate (EUR/tCO2)

To integrate the sustainability evaluation (or SDG assessment) into the
production function, the aggregated variable, GDP, from the economic dimension
is treated as the output, and the other variables from the social and
environmental dimensions are considered as the inputs. In the estimation, we
used the real GDP collected from the World Bank (the same data source) instead
of the adjusted GDP growth. Regarding the indicator employment-to-population
ratio, we transformed it into the actual population according to its
definition in the dashboard reports. Table B1 summarizes the descriptive
statistics for the input and output variables. Note that one can, of course,
introduce such sustainable indicators as the elderly poverty rate, inequality
of income, and freshwater withdrawal described in Lamichhane2020, Singpai, and
the references therein, but should mind the computational burden in the
nonparametric model.

###  4.2 Subset selection and estimates

Since the expectile approaches generally outperform the quantile approaches in
terms of the accuracy comparison and can ensure uniqueness of the estimated
quantile functions, we applied the $\mathcal{L}_{1}$\- and
$\mathcal{L}_{0}$-CER approaches elaborated in Section 2 to calculate the
averaged estimates and examine the effectiveness in subset variable selection.
To obtain a more complete picture of the distribution characteristics of
sustainability, we considered 10 expectile functions with
$\tilde{\tau}\in\\{0.05,0.15,\cdots,0.85,0.95\\}$ by solving the
$\mathcal{L}_{1}$-CER and $\mathcal{L}_{0}$-CER problems. The corresponding
quantile $\tau$ was reported to connect with the usual quantile estimation by
counting the number of negative residuals $\hat{\varepsilon}_{i}^{-}$ that are
greater than zero (Efron1991, *Efron1991).

For the SDG application, we also used the 5-fold cross-validation procedure to
determine the optimal tuning parameters $\lambda$, $k$, and $M$. Specifically,
we calibrated $\lambda$ over 100 values ranging from 0.1 to 3 for the
$\mathcal{L}_{1}$-CER approach. For the $\mathcal{L}_{0}$-CER approach, the
subset size $k$ and constant $M$ were chosen from a set of integers ranging
from 1 to 11 and the set $\\{0.1,0.5,0.8,1,1.5,1.8,2,2.5,3,5\\}$,
respectively. Table 3 demonstrates the averaged estimates of each input
variable on sustainability for the selected expectiles (i.e., 5th, 35th, 65th,
and 95th) in $\mathcal{L}_{1}$-CER and $\mathcal{L}_{0}$-CER estimation.

As shown in Table 3, the $\mathcal{L}_{0}$-norm regularization is a better
choice to select the subset variables, which further confirms the conclusion
of the MC study. The blanks in Table 3 denote that the corresponding variable
is eliminated by the $\mathcal{L}_{1}$-norm or $\mathcal{L}_{0}$-norm
regularization. The number of eliminated irrelevant variables in the
$\mathcal{L}_{0}$-CER estimation is larger than that in the
$\mathcal{L}_{1}$-CER estimation for all expectiles. For example, in the 95th
expectile (i.e., 63rd and 71st quantiles for the $\mathcal{L}_{0}$ penalty and
$\mathcal{L}_{1}$ penalty, respectively), $\mathcal{L}_{0}$-CER selects four
principal input variables, whereas $\mathcal{L}_{1}$-CER estimates eight input
variables. Note that the subset size of actual variables used in the
$\mathcal{L}_{0}$-CER estimation is less than or equal to the optimal
$\hat{k}$.

Table 3: Average estimates of the subset variables for the selected expectiles. |  $\mathcal{L}_{0}$-CER |  $\mathcal{L}_{1}$-CER |  |  $\mathcal{L}_{0}$-CER |  $\mathcal{L}_{1}$-CER |  |  $\mathcal{L}_{0}$-CER |  $\mathcal{L}_{1}$-CER |  |  $\mathcal{L}_{0}$-CER |  $\mathcal{L}_{1}$-CER
---|---|---|---|---|---|---|---|---|---|---|---
E ($\tilde{\tau}$) | 0.05 |  | 0.35 |  | 0.65 |  | 0.95
Q ($\tau$) | 0.14 | 0.16 |  | 0.29 | 0.36 |  | 0.46 | 0.42 |  | 0.63 | 0.71
$\hat{\beta}_{1}$ |  |  |  |  |  |  |  |  |  |  |
$\hat{\beta}_{2}$ | 0.084 | 0.0006 |  | 0.104 | 0.0011 |  |  | 0.0007 |  | 0.046 | 0.0008
$\hat{\beta}_{3}$ | 0.519 |  |  | 0.370 |  |  |  |  |  |  |
$\hat{\beta}_{4}$ | 0.213 | 0.1039 |  | 0.251 | 0.1197 |  | 0.168 | 0.1203 |  | 0.181 | 0.1241
$\hat{\beta}_{5}$ |  | 0.0001 |  |  | 0.0010 |  |  | 0.0001 |  |  | 0.0002
$\hat{\beta}_{6}$ |  | 0.0003 |  |  | 0.0014 |  |  | 0.0007 |  |  | 0.0005
$\hat{\beta}_{7}$ |  |  |  | 0.766 |  |  | 0.417 |  |  |  |
$\hat{\beta}_{8}$ | 0.145 | 0.0008 |  | 0.000 | 0.0031 |  | 0.000 | 0.0009 |  | 0.079 | 0.0011
$\hat{\beta}_{9}$ |  | 0.0033 |  |  | 0.0136 |  |  | 0.0061 |  |  | 0.0086
$\hat{\beta}_{10}$ |  | 0.0122 |  |  | 0.0266 |  |  | 0.0184 |  |  | 0.0246
$\hat{\beta}_{11}$ |  |  |  |  | 0.0052 |  | 0.172 |  |  | 0.161 |
$\hat{\beta}_{12}$ |  | 0.0036 |  |  | 0.0054 |  |  | 0.0012 |  |  | 0.0010
$\hat{\alpha}$ | -10.249 | -0.2700 |  | -10.300 | -0.8202 |  | -1.242 | -0.2407 |  | -4.356 | -0.3166
$\hat{k}$ | 7 | – |  | 7 | – |  | 3 | – |  | 6 | –
$\hat{M}$ | 2.5 | – |  | 2.5 | – |  | 1.5 | – |  | 1.5 | –
$\hat{\lambda}$ | – | 0.95 |  | – | 0.33 |  | – | 2.41 |  | – | 0.31

The best subset variables differ across different expectiles (quantiles). It
can be seen that, for example, variable $x_{2}$, life expectancy, is
eliminated at the 65th expectile in $\mathcal{L}_{0}$-CER but not at the 5th,
35th, or 95th expectiles. Note that the CER approach, a novel data-driven
approach, estimates the production function locally and is robust to
heterogeneity (Kuosmanen2020b, *Kuosmanen2020b; Dai2020, *Dai2020). The
approaches proposed in this paper inherit the good properties from the CER.
Utilizing the full information of each observation, the
$\mathcal{L}_{1}$-/$\mathcal{L}_{0}$-CER approach estimates the expectile
production function based on the actual level of sustainability, taking both
inefficiency and noise into account, but not the full production frontier. The
subset selection is thus not exactly the same in each expectile estimation.

The average estimates of subset variables reveal the heterogeneous effects on
the expected sustainability from the expectile estimation. Note that the
estimated coefficients $\hat{\beta}$ differ across observations. In the
$\mathcal{L}_{0}$-CER estimation, we can see that the impact of the employment
population in the lower 5th and 35th expectile countries are higher than those
in other expectile countries. From the policy-recommendation perspective,
those countries located at the lower expectile can take further actions to
promote employment and even sustainable development. The heterogeneous
policies could be proposed according to the average estimates.

Table 3 shows the strong evidence of economies of scale (i.e., increasing
returns to scale) for OECD countries. We find that the estimated intercepts
$\hat{\alpha}$ are systematically less than zero for all expectiles,
exhibiting increasing returns to scale. In terms of each expectile, all
countries are operating at increasing returns to scale in the lower expectiles
(e.g., 5th and 15th expectiles), while few countries are operating at
decreasing returns to scale in the higher expectiles (e.g., only three units
with positive $\hat{\alpha}$ in the 95th expectile). In the context of
sustainability production, increasing returns suggest that the OECD countries
could improve sustainability through further expansion. The OECD countries,
especially for the smaller counties, for example, could increase inputs in the
social aspect (e.g., increasing the employment population) and decrease
pollutant emissions in the environmental dimension to accelerate sustainable
development.

###  4.3 Performance benchmarking

While the average estimates of subset variables are very interesting and
relevant as such, another main purpose of this paper is to benchmark the
performance of sustainable development for OECD countries. To this end, Table
4 reports the productive performance of OECD countries relative to the
quantile frontier based on the $\mathcal{L}_{0}$-CER estimation. For the sake
of illustration, following Kuosmanen2020, Table 4 presents the full list of
OECD countries under four groups: EU-15, EU transition economies, European
Free Trade Association (EFTA), and Non-European OECD. Note that the higher the
expectile is, the better the relative performance.

Recall that the $\mathcal{L}_{0}$-CER approach estimates the quantile
(expectile) production function instead of the full frontier. This would
suggest that the quantile production function estimation can serve as the
benchmarking for structure analysis. There is thus a good connection between
the quantile benchmarking and metafrontier benchmarking. One could build a
quantile metafrontier to analyze a group of units at the same quantile level
(see Lai2018, *Lai2018). In the present context, the quantile output technical
efficiency can be measured by
$\text{TE}(\tau)=y/\hat{Q}(\tau\,|\,{\boldsymbol{x}})$. However, for the
purpose of performance benchmarking, one can directly compare the relative
locations of the quantiles (e.g., Kuosmanen2020, *Kuosmanen2020; this paper)
instead of the quantile-based technical efficiency (e.g., Behr2010, *Behr2010;
Lai2018, *Lai2018).

Table 4: Productive performances of countries relative to the quantile
frontiers.

| 2017 | 2019 | 2020 |  | 2017 | 2019 | 2020
---|---|---|---|---|---|---|---
EU-15 | EU transition economies
Luxembourg | 10 | 10 | 10 | Estonia | 10 | (6, 7) | 10
United Kingdom | 1 | 1 | (1, 2) | Slovenia | 10 | 10 | 10
Denmark | 10 | (9, 10) | 10 | Poland | (5, 6) | (6, 7) | (2, 3)
Sweden | (9, 10) | (9, 10) | (9, 10) | Hungary | (6, 7) | (1, 2) | (5, 6)
Germany | (2, 3) | (1, 2) | 1 | Latvia | 10 | 10 | 10
Finland | (8, 9 | (8, 9) | (9, 10) | Czech Republic | (9, 10) | (9, 10) | (8, 9)
Belgium | 10 | (9, 10) | (9, 10) | Slovakia | (9, 10) | (9, 10) | (7, 8)
Netherlands | (6, 7) | (6, 7) | (6, 7) |  |  |  |
France | (3, 4) | (4, 5) | (7, 8) | Non-European OECD
Italy | (4, 5) | (4, 5) | (3, 4) | Canada | 1 | 1 | (1, 2)
Austria | (7, 8) | (8, 9) | (7, 8) | Israel | (6, 7) | (6, 7) | (6, 7)
Greece | (6, 7) | (7, 8) | (7, 8) | United States | 10 | 10 | 10
Ireland | 10 | 10 | (5, 6) | Australia | (1, 2) | (1, 2) | (2, 3)
Spain | (1, 2) | (1, 2) | (2, 3) | New Zealand | (1, 2) | (1, 2) | (6, 7)
Portugal | (5, 6) | (4, 5) | (2, 3) | South Korea | (4, 5) | (5, 6) | (6, 7)
|  |  |  | Chile | (6, 7) | (6, 7) | (6, 7)
EFTA | Japan | 1 | 1 | (5, 6)
Norway | (9, 10) | (9, 10) | (9, 10) | Mexico | (1, 2) | (4, 5) | 1
Iceland | 10 | 10 | 10 | Turkey | (7, 8) | (8, 9) | (6, 7)
Switzerland | 10 | 10 | 10 |  |  |  |

  * •

Note: 1) Legend: 10 = above 95th quantile; (9,10) = between 85th and 95th
quantile; (8,9) = between 75th and 85th quantile; (7,8) = between 65th and
75th quantile; (6,7) = between 55th and 65th quantile; (5,6) = between 45th
and 55th quantile; (4,5) = between 35th and 45th quantile; (3,4) = between
25th and 35th quantile; (2,3) = between 15th and 25th quantile; (1,2) =
between 5th and 15th quantile; 1 = below 5th quantile. 2) Although the United
Kingdom left the EU in January 2020, we still treat it as an EU country owing
to the sample period.

Several interesting findings can be drawn from Table 4. First, in general,
EFTA countries performed best in sustainable development, while the ranking of
Non-European countries was relatively lower than that of the other groups.
With respect to time, there was an increase in sustainability (i.e., the level
of performance) for most countries and a decrease for a few other countries.
Although Table 4 provides a similar ranking to the overall performance from
the SDG Dashboard Reports (Sachs2017, *Sachs2017, *Sachs2019, *Sachs2020) and
Lamichhane2020, there was a minor difference in ranking for certain countries
(c.f., the ranking of Germany) due to the fact that all 17 SDGs with 231
unique indicators were evaluated in their studies. Instead, we consider a
large portion of environmental indicators as input variables (see Table 2).

Second, countries from the EU transition economies can serve as the benchmark.
As shown in Table 4, the developed countries did not necessarily perform best
in sustainable development, whereas the less-developed countries had a chance
to become better (e.g., Estonia and Slovenia). This partly demonstrates the
advantages of the CQR approach over the full frontier estimation in structure
analysis. If a unit is located in the interior of a production probability
set, it will be treated as inefficient in the traditional frontier, in
contrast to quantile-based efficiency analysis, where we treat it as an
efficient unit. Therefore, less-developed countries, such as the Czech
Republic and Slovenia, could perform best in sustainable development. As
stated earlier, this might be related to the selection of input-output
variables.

Third, Table 4 reveals the inequality of the relative performances of
sustainability among the OECD countries. The large variation in expectiles
would suggest that the plan, action, and implementation of facilitating
sustainable development in OECD counties are inefficient. The inequality of
sustainability has become a barrier to the 2030 Agenda for Sustainable
Development for OECD countries, even to all countries worldwide. Consequently,
for the upper-expectile countries, they can achieve the 17 SDGs easily, but it
is difficult for lower-expectile countries to achieve those goals. Table 4
also suggests that there is large room for OECD countries to improve in terms
of sustainability. This considerable room for improvement calls for ambitious
plans and actions to achieve the 17 SDGs. A whole plan, for example, in terms
of SO2 emissions, similar to the Paris Agreement, could be considered and
implemented in OECD countries.

##  5 Conclusions

Addressing the curse of dimensionality under a limited sample size and sparse
data space remains a challenge in production efficiency analysis. In this
paper, we developed a new $\mathcal{L}_{0}$-norm regularization approach to
CQR and CER for subset variable selection. The paper investigated the finite
sample performance of the proposed $\mathcal{L}_{0}$-norm regularization in
contrast to the commonly used $\mathcal{L}_{1}$-norm regularization via an MC
study (Section 3). The proposed $\mathcal{L}_{0}$-norm regularization approach
was further applied to the SDG data of OECD countries for the years 2017, 2019
and 2020 to evaluate the performance of sustainability and compare the
methods’ performances empirically (Section 4). The proposed
$\mathcal{L}_{0}$-norm regularization approach can more effectively address
the curse of dimensionality via subset variable selection in multidimensional
spaces.

The evidence from the MC study suggests that $\mathcal{L}_{0}$-CQR has a major
advantage over $\mathcal{L}_{1}$-CQR in subset variable selection, especially
in the high-SNR regime. $\mathcal{L}_{1}$-CQR is unable to completely
eliminate the irrelevant variables compared with the true support set in
certain cases, whereas $\mathcal{L}_{0}$-CQR can successfully select the true
relevant variables at least partially. There is no systematic increasing or
decreasing trend as the quantile varies. The effects of the SNR and sample
size are also observed: the higher the SNR level is, the lower the prediction
error and accuracy; the larger the sample size is, the lower the prediction
error and accuracy. Furthermore, we find that the expectile approach generally
performs better than the quantile approach.

The SDG application results demonstrate that the $\mathcal{L}_{0}$-CER
approach can eliminate more irrelevant variables than the
$\mathcal{L}_{1}$-CER approach, indicating a better performance in subset
selection. The relative expectile ranking also helps reveal the heterogeneity
of performance in sustainable development. The results also suggest that the
countries in the EFTA group outperform those in other groups, that the
countries in the EU transition economies can serve as the benchmark, and that
there is large inequality of the relative performance of sustainability among
the OECD countries.

The findings drawn from this study can provide insight into the relationship
between $\mathcal{L}_{1}$-norm regularization and $\mathcal{L}_{0}$-norm
regularization in the context of production economics, especially in the era
of big data. Furthermore, the proposed approach is applicable to other
relative performance evaluations (e.g., schools, hospitals, and electricity
distribution system operators). Finding a more efficient procedure to
determine the tuning parameter is one of the fascinating avenues for future
research. The other, as noted in Section 2.2, is to provide a formal proof
that application of the $\mathcal{L}_{1}$-norm squared regularization to the
CQR problem can ensure the uniqueness of the subgradients.

## Acknowledgments

The author is indebted to Timo Kuosmanen for his valuable guidance, advice,
and comments. I also gratefully acknowledge the computational resources
provided by the Aalto Science-IT project and financial support from the
Foundation for Economic Education (Liikesivistysrahasto) (nos. 180019, 190073)
and the HSE Support Foundation (no. 11–2290).

## References

  * [1] Adler Yazhemsky2010Adler2010a Adler, N. Yazhemsky, E. 2010, Improving discrimination in data envelopment analysis: PCA-DEA or variable reduction, European Journal of Operational Research 202, 273–284.
  * [2] [Balázs et al.]Balázs, György Szepesvári2015Balazs2015 Balázs, G., György, A. Szepesvári, C. 2015, Near-optimal max-affine estimators for convex regression, 18th Artificial Intelligence and Statistics, PMLR, p. 38:56–64.
  * [3] Behr2010Behr2010 Behr, A. 2010, Quantile regression for robust bank efficiency score estimation, European Journal of Operational Research 200, 568–581.
  * [4] [Benítez-Peña et al.]Benítez-Peña, Bogetoft Romero Morales2020Benitez-Pena2019 Benítez-Peña, S., Bogetoft, P. Romero Morales, D. 2020, Feature selection in data envelopment analysis: A mathematical optimization approach, Omega 96, 102068.
  * [5] [Bertsimas et al.]Bertsimas, King Mazumder2016Bertsimas2016 Bertsimas, D., King, A. Mazumder, R. 2016, Best subset selection via a modern optimization lens, Annals of Statistics 44, 813–852.
  * [6] Bertsimas Mundru2021Bertsimas2020 Bertsimas, D. Mundru, N. 2021, Sparse convex regression, INFORMS Journal on Computing 33, 262–279.
  * [7] [Bertsimas et al.]Bertsimas, Pauphilet Van Parys2017Bertsimas2017 Bertsimas, D., Pauphilet, J. Van Parys, B. 2017, Sparse Classification: A scalable discrete optimization perspective. https://arxiv.org/abs/1710.01352
  * [8] Bertsimas van Parys2020Bertsimas2020b Bertsimas, D. van Parys, B. 2020, Sparse high-dimensional regression: Exact scalable algorithms and phase transitions, Annals of Statistics 48, 300–323.
  * [9] Chen Lee2021Chen2021a Chen, L.-Y. Lee, S. 2021, Binary classification with covariate selection through L0-penalised empirical risk minimisation, Econometrics Journal 24, 103–120.
  * [10] [Chen et al.]Chen, Lin Sen2020chen2020c Chen, X., Lin, Q. Sen, B. 2020, On degrees of freedom of projection estimators with applications to multivariate nonparametric regression, Journal of the American Statistical Association 115(529), 173–186.
  * [11] [Chen et al.]Chen, Tsionas Zelenyuk2021Chen2021 Chen, Y., Tsionas, M. G. Zelenyuk, V. 2021, LASSO+DEA for small and big wide data, Omega 102, 102419.
  * [12] [Costanza et al.]Costanza, Daly, Fioramonti, Giovannini, Kubiszewski, Mortensen, Pickett, Ragnarsdottir, De Vogli Wilkinson2016Costanza2016 Costanza, R., Daly, L., Fioramonti, L., Giovannini, E., Kubiszewski, I., Mortensen, L. F., Pickett, K. E., Ragnarsdottir, K. V., De Vogli, R. Wilkinson, R. 2016, Modelling and measuring sustainable wellbeing in connection with the UN Sustainable Development Goals, Ecological Economics 130, 350–355.
  * [13] [Dai et al.]Dai, Zhou Kuosmanen2020Dai2020 Dai, S., Zhou, X. Kuosmanen, T. 2020, Forward-looking assessment of the GHG abatement cost: Application to China, Energy Economics 88, 104758.
  * [14] Dalampira Nastis2020Dalampira2020 Dalampira, E. S. Nastis, S. A. 2020, Mapping sustainable development goals: A network analysis framework, Sustainable Development 28, 46–55.
  * [15] Dulá2011Dula2011 Dulá, J. H. 2011, An algorithm for data envelopment analysis, INFORMS Journal on Computing 23, 284–296.
  * [16] [Dunning et al.]Dunning, Huchette Lubin2017Dunning2017 Dunning, I., Huchette, J. Lubin, M. 2017, JuMP: A modeling language for mathematical optimization, SIAM review 59, 295–320.
  * [17] [Dyson et al.]Dyson, Allen, Camanho, Podinovski, Sarrico Shale2001Dyson2001 Dyson, R. G., Allen, R., Camanho, A. S., Podinovski, V. V., Sarrico, C. S. Shale, E. A. 2001, Pitfalls and protocols in DEA, European Journal of Operational Research 132, 245–259.
  * [18] Efron1991Efron1991 Efron, B. 1991, Regression percentiles using asymmetric squared error loss, Statistica Sinica 1, 93–125.
  * [19] Fan Song2010Fan2010 Fan, J. Song, R. 2010, Sure independence screening in generalized linear models with NP-dimensionality, Annals of Statistics 38, 3567–3604.
  * [20] [Hastie et al.]Hastie, Tibshirani Tibshirani2020Hastie2020 Hastie, T., Tibshirani, R. Tibshirani, R. 2020, Best subset, forward stepwise or Lasso? Analysis and recommendations based on extensive comparisons, Statistical Science 35, 579–592.
  * [21] Homburg2001Homburg2001 Homburg, C. 2001, Using data envelopment analysis to benchmark activities, International Journal of Production Economics 73, 51–58.
  * [22] [Huan et al.]Huan, Yu, Liang Burgman2021Huan2021 Huan, Y., Yu, Y., Liang, T. Burgman, M. 2021, A method for assessing the impacts of an international agreement on regional progress towards Sustainable Development Goals, Science of the Total Environment 785, 147336.
  * [23] Keshvari2018Keshvari2018 Keshvari, A. 2018, Segmented concave least squares: A nonparametric piecewise linear regression, European Journal of Operational Research 266, 585–594.
  * [24] Koenker Bassett1978Koenker1978 Koenker, R. Bassett, G. 1978, Regression quantiles, Econometrica 46, 33–50.
  * [25] Kuosmanen2008Kuosmanen2008 Kuosmanen, T. 2008, Representation theorem for convex nonparametric least squares, Econometrics Journal 11, 308–325.
  * [26] [Kuosmanen et al.]Kuosmanen, Johnson Saastamoinen2015Kuosmanen2015d Kuosmanen, T., Johnson, A. Saastamoinen, A. 2015, Stochastic nonparametric approach to efficiency analysis: A unified framework, in J. Zhu (ed.), Data Envelopment Analysis, Springer, Boston, MA, chapter 7, pp. 191–244.
  * [27] Kuosmanen Zhou2021Kuosmanen2020b Kuosmanen, T. Zhou, X. 2021, Shadow prices and marginal abatement costs: Convex quantile regression approach, European Journal of Operational Research 289, 666–675.
  * [28] [Kuosmanen et al.]Kuosmanen, Zhou Dai2020Kuosmanen2020 Kuosmanen, T., Zhou, X. Dai, S. 2020, How much climate policy has cost for OECD countries?, World Development 125, 104681.
  * [29] [Lai et al.]Lai, Huang Fu2018Lai2018 Lai, H. p., Huang, C. J. Fu, T. T. 2018, Estimation of the production profile and metafrontier technology gap: A quantile approach, Empirical Economics 58, 2709–2731.
  * [30] [Lamichhane et al.]Lamichhane, Eğilmez, Gedik, Bhutta Erenay2020Lamichhane2020 Lamichhane, S., Eğilmez, G., Gedik, R., Bhutta, M. K. S. Erenay, B. 2020, Benchmarking OECD countries’ sustainable development performance: A goal-specific principal component analysis approach, Journal of Cleaner Production 287, 125040.
  * [31] Lavergne Patilea2008Lavergne2008 Lavergne, P. Patilea, V. 2008, Breaking the curse of dimensionality in nonparametric testing, Journal of Econometrics 143, 103–122.
  * [32] Lee Cai2020Lee2020 Lee, C. Y. Cai, J. Y. 2020, LASSO variable selection in data envelopment analysis with small datasets, Omega 91, 102019.
  * [33] [Lee et al.]Lee, Johnson, Moreno-Centeno Kuosmanen2013Lee2013 Lee, C. Y., Johnson, A. L., Moreno-Centeno, E. Kuosmanen, T. 2013, A more efficient algorithm for Convex Nonparametric Least Squares, European Journal of Operational Research 227, 391–400.
  * [34] [Li et al.]Li, Shi, Yang Liang2017Li2017 Li, Y., Shi, X., Yang, M. Liang, L. 2017, Variable selection in data envelopment analysis via Akaike’s information criteria, Ann Oper Res 253, 453–476.
  * [35] [Lin et al.]Lin, Sun Toh2020Lin2020 Lin, M., Sun, D. Toh, K.-C. 2020, An augmented Lagrangian method with constraint generations for shape-constrained convex regression problems. http://arxiv.org/abs/2012.04862
  * [36] [Mazumder et al.]Mazumder, Choudhury, Iyengar Sen2019Mazumder2019 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, 318–331.
  * [37] Nataraja Johnson2011Nataraja2011 Nataraja, N. R. Johnson, A. L. 2011, Guidelines for using variable selection techniques in data envelopment analysis, European Journal of Operational Research 215, 662–669.
  * [38] Natarajan1995Natarajan1995 Natarajan, B. K. 1995, Sparse approximate solutions to linear systems, SIAM Journal on Computing 24, 227–234.
  * [39] [Pastor et al.]Pastor, Ruiz Sirvent2002Pastor2002 Pastor, J. T., Ruiz, J. L. Sirvent, I. 2002, A statistical test for nested radial DEA models, Operations Research 50, 728–735.
  * [40] Qin Song2014Qin2014 Qin, Z. Song, I. 2014, Joint variable selection for data envelopement analysis via group sparsity. http://arxiv.org/abs/1402.3740
  * [41] [Sachs et al.]Sachs, Schmidt-Traub, Kroll, Durand-Delacre Teksoz2017Sachs2017 Sachs, J., Schmidt-Traub, G., Kroll, C., Durand-Delacre, D. Teksoz, K. 2017, SDG index and dashboards report 2017, Bertelsmann Stiftung and Sustainable Development Solutions Network (SDSN), New York.
  * [42] [Sachs et al.]Sachs, Schmidt-Traub, Kroll, Lafortune Fuller2019Sachs2019 Sachs, J., Schmidt-Traub, G., Kroll, C., Lafortune, G. Fuller, G. 2019, Sustainable development report 2019, Bertelsmann Stiftung and Sustainable Development Solutions Network (SDSN), New York.
  * [43] [Sachs et al.]Sachs, Schmidt-Traub, Kroll, Lafortune, Fuller Woelm2020Sachs2020 Sachs, J., Schmidt-Traub, G., Kroll, C., Lafortune, G., Fuller, G. Woelm, F. 2020, The Sustainable Development Goals and COVID-19. Sustainable Development Report 2020, Cambridge University Press, Cambridge.
  * [44] Singpai Wu2020Singpai Singpai, B. Wu, D. 2020, Using a DEA-AutoML approach to track SDG achievements, Sustainability 12, 1–26.
  * [45] [Sinha et al.]Sinha, Malo Kuosmanen2015Sinha2015 Sinha, A., Malo, P. Kuosmanen, T. 2015, A multiobjective exploratory procedure for regression model selection, Journal of Computational and Graphical Statistics 24(1), 154–182.
  * [46] Stone1980Stone2007 Stone, C. J. 1980, Optimal rates of convergence for nonparametric estimators, The Annals of Statistics 8, 1348–1360.
  * [47] [Su et al.]Su, Bogdan Candès2017Su2017 Su, W., Bogdan, M. Candès, E. 2017, False discoveries occur early on the lasso path, Annals of Statistics 45, 2133–2150.
  * [48] Tibshirani1996Tibshirani1996 Tibshirani, R. 1996, Regression shrinkage and selection via the Lasso, Journal of the Royal Statistical Society. Series B (Methodological) 58, 267–288.
  * [49] Wagner Shimshak2007Wagner2007 Wagner, J. M. Shimshak, D. G. 2007, Stepwise selection of variables in data envelopment analysis: Procedures and managerial perspectives, European Journal of Operational Research 180, 57–67.
  * [50] [Wang et al.]Wang, Wang, Dang Ge2014Wang2014c Wang, Y., Wang, S., Dang, C. Ge, W. 2014, Nonparametric quantile frontier estimation under shape restriction, European Journal of Operational Research 232, 671–678.
  * [51] Wilson2018Wilson2018 Wilson, P. W. 2018, Dimension reduction in nonparametric models of production, European Journal of Operational Research 267, 349–367.
  * [52]

## Appendix

###  A The CNLS-A algorithm

The CNLS-A algorithm, an extension of Bertsimas2020, is designed for improving
the computational efficiency in solving the penalized CQR and CER problems.
Specifically, to speed up the calculation, the following strategy is applied:
solving the reduced master problem that contains only a few constraints and
then iteratively adding the violated constraints in a delayed manner.111As in
Bertsimas2020, we use the same notations: the master problem is the original
problem with $n^{2}$ linear constraints (i.e., the formulation (A.1)); the
reduced master problem is a problem with the same objective function and
decision variables but with a subset of constraints.

For the sake of elaboration, the CQR problem is taken as an example to compare
the computational efficiency between the CNLS-A and CNLS-G algorithms.
Therefore, the CQR problem is defined as (Kuosmanen2015d, *Kuosmanen2015d)

$\displaystyle\underset{\mathbf{\alpha},\mathbf{{\boldsymbol{\beta}}},{{\mathbf{\varepsilon}}^{\text{+}}},{{\mathbf{\varepsilon}}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(A.1) s.t. $\displaystyle
y_{i}=\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}+\varepsilon^{+}_{i}-\varepsilon^{-}_{i}$
$\displaystyle\forall i$
$\displaystyle\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{{{\boldsymbol{x}}}_{i}}\leq\mathbf{\alpha}_{h}+{\boldsymbol{\beta}}_{h}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
$\displaystyle\forall i,h$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

To adapt the existing cutting-plan algorithm (Bertsimas2020, *Bertsimas2020),
the CQR problem first needs to be reformulated. The alternative CQR problem is

$\displaystyle\underset{\hat{y},\mathbf{\beta},{{\mathbf{\varepsilon}}^{\text{+}}},{{\mathbf{\varepsilon}}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(A.2) s.t. $\displaystyle
y_{i}-\hat{y}_{i}=\varepsilon_{i}^{+}-\varepsilon_{i}^{-}$
$\displaystyle\forall i$
$\displaystyle\hat{y}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{j}-{\boldsymbol{x}}_{i})\geq\hat{y}_{j}$
$\displaystyle\forall i,j$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

where $\hat{y}_{i}$, a new decision variable, is the estimated quantile
function. Note that the alternative formulation (A.2) is derived from Problem
(A.1) by applying a mathematical transformation. For the formulation (A.1),
$\hat{y}_{i}$ is used to replace
$\mathbf{\alpha}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}{\boldsymbol{x}}_{i}$
in the first constraint, and the Afriat equality constraint is reconstructed
by substituting $\hat{y}_{i}$ for the left-hand side and substituting out
$\mathbf{\alpha}_{j}$ using
$\mathbf{\alpha}_{j}=\hat{y}_{j}-{\boldsymbol{\beta}}_{j}^{{}^{\prime}}{\boldsymbol{x}}_{j}$
for the right-hand side.

There are two commonly used approaches to form the initial $n-1$ constraints:
the spanning path (SP) based on the Euclidean distances between observations
and the minimum spanning tree (MST) among observations. The initial reduced
master problem of CQR is:

$\displaystyle\underset{\hat{y},\mathbf{\beta},{{\mathbf{\varepsilon}}^{\text{+}}},{{\mathbf{\varepsilon}}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(A.3) s.t.
$\displaystyle\hat{y}_{i_{1}}+{\boldsymbol{\beta}}_{i_{1}}^{{}^{\prime}}({\boldsymbol{x}}_{i_{2}}-{\boldsymbol{x}}_{i_{1}})\geq\hat{y}_{i_{2}}$
$\displaystyle\hat{y}_{i_{2}}+{\boldsymbol{\beta}}_{i_{2}}^{{}^{\prime}}({\boldsymbol{x}}_{i_{3}}-{\boldsymbol{x}}_{i_{2}})\geq\hat{y}_{i_{3}}$
$\displaystyle\vdots$
$\displaystyle\hat{y}_{i_{n-1}}+{\boldsymbol{\beta}}_{i_{n-1}}^{{}^{\prime}}({\boldsymbol{x}}_{i_{n}}-{\boldsymbol{x}}_{i_{n-1}})\geq\hat{y}_{i_{n}}$
$\displaystyle y_{i}-\hat{y}_{i}=\varepsilon_{i}^{+}-\varepsilon_{i}^{-}$
$\displaystyle\forall i$
$\displaystyle{\boldsymbol{\beta}}_{i}\geq{\boldsymbol{0}}$
$\displaystyle\forall i$ $\displaystyle\varepsilon_{i}^{+}\geq 0,\
\varepsilon_{i}^{-}\geq 0$ $\displaystyle\forall i$

For a given solution (i.e.,
$\hat{y},{\boldsymbol{\beta}}^{{}^{\prime}}_{i},\varepsilon_{i}^{-},\varepsilon_{i}^{+}$)
to the reduced master problem (A.3), we need to check whether it is a feasible
solution to the master problem (A.2). If yes, it is also the optimal solution
for the master problem (A.2); otherwise, we have to efficiently find a
violated Afriat constraint. Finding a violated constraint is regarded as a
separation problem (Bertsimas2020, *Bertsimas2020). Therefore, the separation
problem finds the minimal index $m(i)$ for each observation $i$ and checks if
the corresponding smallest value is less than zero.

$\displaystyle m(i)=\operatorname*{arg\,min}_{1\leq m\leq
n}\\{\hat{y}_{i}-\hat{y}_{m}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{m}-{\boldsymbol{x}}_{i})\\}$
(A.4)

Furthermore, the following violated constraint is added to the reduced master
problem (A.3) only if it is smaller than a given tolerance.

$\displaystyle\hat{y}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{m(i)}-{\boldsymbol{x}}_{i})\geq\hat{y}_{m(i)}$
(A.5)

The $m$th iteration by $T_{m}$ is then solved, and if $\min_{1\leq m\leq
n}\\{\hat{y}_{i}-\hat{y}_{m}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{m}-{\boldsymbol{x}}_{i})\\}\geq-
tol$, a final optimal solution to the master problem (A.2) is found. The $m$th
iteration problem is as follows:

$\displaystyle\underset{\hat{y},\mathbf{\beta},{{\mathbf{\varepsilon}}^{\text{+}}},{{\mathbf{\varepsilon}}^{-}}}{\mathop{\min}}$
$\displaystyle\,\tau\sum\limits_{i=1}^{n}{\varepsilon_{i}^{+}}+(1-\tau)\sum\limits_{i=1}^{n}{\varepsilon_{i}^{-}}$
(A.6) s.t.
$\displaystyle\hat{y}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{j}-{\boldsymbol{x}}_{i})\geq\hat{y}_{j}$
$\displaystyle\forall(i,j)\in T_{0}$
$\displaystyle\hat{y}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{j}-{\boldsymbol{x}}_{i})\geq\hat{y}_{j}$
$\displaystyle\forall(i,j)\in T_{1}$ $\displaystyle\vdots$
$\displaystyle\hat{y}_{i}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{j}-{\boldsymbol{x}}_{i})\geq\hat{y}_{j}$
$\displaystyle\forall(i,j)\in T_{m}$

The CNLS-A algorithm procedure is summarized as:

Data:
$\\{({\boldsymbol{x}}_{i},y_{i})\in{\mathbb{R}}^{d}\times{\mathbb{R}}:i=1,\cdots,n\\}$
and $tol=0.01$

1 $\mathrm{out}=0$ and $m=0$;

2 Solve the reduced master problem (A.6);

3 while  _$\mathrm{out}=0$_ do

4 for  _$1\leq i\leq n$_ do

5 Solve the separation problem (A.4) to find a minimal index $m(i)$;

6 Add the corresponding violated constraint (A.5) to Problem (A.6);

7

8 if  _$\min_{1\leq m\leq
n}\\{\hat{y}_{i}-\hat{y}_{m}+{\boldsymbol{\beta}}_{i}^{{}^{\prime}}({\boldsymbol{x}}_{m}-{\boldsymbol{x}}_{i})\\}
<-tol$_ then

9 Resolve the reduced master problem (A.6) with the newly added constraint set
$T_{m+1}$;

10

11 else

12 $\mathrm{out}=1$;

13

14 $m=m+1$;

15

Result:
$\hat{y},{\boldsymbol{\beta}}^{{}^{\prime}}_{i},\varepsilon_{i}^{-},\varepsilon_{i}^{+}$

Algorithm 1 CNLS-A algorithm for solving the CQR problem (A.2).

Fig. A1 depicts a performance comparison between the CNLS-G and CNLS-A
algorithms in solving the CQR and CER problems.222In this paper, we slightly
adapt the CNLS-G algorithm that is design for the usual CNLS problem (Lee2013,
*Lee2013) to solve the CQR problem, and we use our developed pyStoNED package
(https://github.com/ds2010/pyStoNED) to run the CNLS-G algorithm.  Since there
is no systematic difference in quantile (see the MC study in Section 3), we
consider only the effect of the number of observations and the number of input
variables. We repeat the simulation 10 times to average the metrics of
interest (i.e., running time and number of constraints). Furthermore, we use
the MST approach to form the initial reduced constraints.

As expected, the CNLS-A algorithm requires less time to solve the CQR and CER
problems than the CNLS-G algorithm does. In contrast to the running time of
CNLS-A, CNLS-G increases rapidly as the input dimension ($d$) increases,
especially when the high-dimensional data are estimated with a large sample
size ($n>300$). For example, according to the upper-right panel in Fig. A1, we
find that the CNLS-G algorithm takes 25 times and 124 times longer than the
CNLS-A algorithm in solving the CQR problem with $n=100$ and $n=500$,
respectively. Compared to the CQR approach, CER runs faster in both algorithms
due to the stronger convexity in the objective function.

Figure A1: Performance comparison of the CNLS-A and CNLS-G algorithms to solve
the CQR and CER problems with $\tau=0.9$.

Regarding the average number of constraints, similar to the average running
time, the CNLS-A algorithm can more efficiently find the violated constraints
than can the CNLS-G algorithm in all the scenarios. CNLS-G includes an
incremental increase in the constraints in the optimization as the number of
input dimensions or the number of observations increases. By contrast, the
constraints of the CNLS-A algorithm undergo a very small change. Furthermore,
the CER approach uses less constraints than does the CQR approach.

###  B Additional tables and figures

Table B1: Summary of statistics for the input and output variables

Variable |  | Unit | Mean | Median | Min. | Max. | Std.Dev.
---|---|---|---|---|---|---|---
Economic: |  |  |  |  |  |  |
$y$ (GDP) |  | 1012 international $ | 1.69 | 0.50 | 0.02 | 20.59 | 3.44
Social: |  |  |  |  |  |  |
$x_{1}$ |  | % | 0.44 | 0.29 | 0.00 | 2.18 | 0.44
$x_{2}$ |  | years | 77.53 | 79.20 | 66.20 | 84.20 | 5.16
$x_{3}$ |  | ladder score, 0-10 | 6.70 | 6.88 | 5.19 | 7.86 | 0.70
$x_{4}$ (Employment) |  | 109 persons | 16.38 | 4.63 | 0.14 | 153.55 | 27.40
$x_{5}$ |  | 1-100 | 35.21 | 33.74 | 25.59 | 57.83 | 6.94
Environmental: |  |  |  |  |  |  |
$x_{6}$ |  |  $\mu$g/m3 | 13.79 | 12.03 | 5.40 | 44.31 | 7.27
$x_{7}$ |  | kg/day/capita | 1.91 | 1.89 | 0.88 | 4.54 | 0.65
$x_{8}$ |  | kg/capita | 18.19 | 19.80 | 6.50 | 28.50 | 5.17
$x_{9}$ |  | kg/capita | 42.22 | 25.15 | 1.70 | 344.94 | 50.90
$x_{10}$ |  | kg/capita | 46.35 | 42.27 | 25.19 | 139.80 | 19.84
$x_{11}$ |  | tCO2/capita | 7.96 | 7.36 | 3.44 | 18.70 | 3.65
$x_{12}$ |  | EUR/tCO2 | 18.38 | 12.47 | 0.01 | 66.95 | 16.12

  * •

Notes: 1) To calculate the absolute values of the population and GDP, we also
collect two additional indicators: GDP and PPP (constant 2017 international
$).

Figure B1: Prediction errors of the $\mathcal{L}_{1}$-CQR and
$\mathcal{L}_{0}$-CQR approaches with $n=500$ and $k_{\text{true}}=2$. Figure
B2: Prediction errors of the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
approaches with $n=500$ and $k_{\text{true}}=4$. Figure B3: Accuracies of the
$\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR approaches with $n=100$ and
$k_{\text{true}}=2$. Figure B4: Accuracies of the $\mathcal{L}_{1}$-CQR and
$\mathcal{L}_{0}$-CQR approaches with $n=100$ and $k_{\text{true}}=4$. Figure
B5: Accuracies of the $\mathcal{L}_{1}$-CQR and $\mathcal{L}_{0}$-CQR
approaches with $n=500$ and $k_{\text{true}}=2$.


generality that $\min X\notin\Omega$. Consider $R\in{\cal R}$ such that for
each $i\in A$, $\rho(R_{i})=\min X$ and for each $i\in D$, $\delta(R_{i})=\max
X$. By construction, $\min X\,P_{i}\,x$ for each $x\in\Omega$ and each $i\in
N$. Then, $\min X$ Pareto dominates $f(R)$ because $f(R)\in\Omega$.

  * •

Suppose $\\{\min X,\max X\\}\subset\Omega$. Let $x\not\in\Omega$ and let
$y\in\mathbf{int}(\Omega)$ be the closest alternative to $x$ in
$\mathbf{int}(\Omega)$. Hence, $(x,y)\cap\Omega=\emptyset$ whenever $x<y$ and
$(y,x)\cap\Omega=\emptyset$ whenever $x>y$. Observe that $x$ and $y$ always
exist because $\Omega\notin\\{X,\\{\min X,\max X\\}\\}$. Consider
$R^{\prime}\in{\cal R}$ such that for each $i\in A$, $[\rho(R^{\prime}_{i})=x$
and $p(R^{\prime}_{i})=y]$, and for each $i\in D$, $\delta(R^{\prime}_{i})=y$.
Then, $\\{j\in A:p(R^{\prime}_{j})\leq^{*}y\\}=A$ and for each $\alpha<^{*}y$,
$\\{j\in A:p(R^{\prime}_{j})\leq^{*}\alpha\\}=\emptyset$. Therefore, by
Theorem 1, $\omega(p(R^{\prime}))=y$ and $f(R^{\prime})=y$. A contradiction is
reached by observing that $x\,P^{\prime}_{i}\,y$ for each $i\in N$.

Finally, we prove $(iii)$. Consider any SP rule $f$ such that $\Omega=\\{\min
X,\max X\\}\neq X$. We divide the proof into three steps:

Step 1: We show that if $A=\emptyset$, then $f$ is PE.

Since $A=\emptyset$, $N=D\neq\emptyset$. Then, for each $R\in{\cal R}$,
$\omega(p(R))=(\min X,\max X)$. Observe that for each agent, either $\min X$
or $\max X$ is the most preferred alternative of $X$. Consider first
$R\in{\cal R}$ such that all agents have the same most preferred alternative
of $X$. Assume, without loss of generality, that this most preferred
alternative is $\min X$. Then, $L_{(\min X,\max X)}(R)=N$ and, by Theorem 1,
$g_{(\min X,\max X)}(R)=\min X$ and $f(R)=\min X$. Consider now
$R^{\prime}\in{\cal R}$ such that not all agents have the same most preferred
alternative. Then, there is a non-empty set of agents $S\subset D$ whose most
preferred alternative is $\min X$, while for the remaining agents, $N\setminus
S$, $\max X$ is the most preferred alternative. Since $\Omega=\\{\min X,\max
X\\}$, we have $f(R^{\prime})\in\\{\min X,\max X\\}$. If, on the one hand,
$f(R^{\prime})=\min X$, then $f(R^{\prime})\,P^{\prime}_{i}\,x$ for each $i\in
S$ and each $x\in X\setminus\\{\min X\\}$. If, on the other hand,
$f(R^{\prime})=\max X$, then $f(R^{\prime})\,P^{\prime}_{i}\,x$ for each $i\in
N\setminus S$ and each $x\in X\setminus\\{\max X\\}$. Hence, $f$ is PE.

Step 2: We show that if $A\neq\emptyset$ and $r_{\omega}=\\{(\min X,\max
X)\\}$, then $f$ is PE.

Since $\max\Omega=\max X$ exists and $\max\Omega\notin r_{\omega}$, we know,
by Theorem 1 (exactly by condition $(iii)$ of Definition 1) that
$\emptyset\in{\cal L}(\min X,\max X)$. Then, $\emptyset$ is a minimal
coalition of ${\cal L}(\min X,\max X)$. By Theorem 1 (exactly by condition
$(ii)$ of Definition 3), there is $C\in W(g_{(\min X,\max X)})$ such that
$C\cap A=\emptyset$. Thus, $C\subseteq D$ and $C\neq\emptyset$. Note that for
each agent in $D$, either $\min X$ or $\max X$ is the most preferred
alternative of $X$. Observe also that since $r_{\omega}=\\{(\min X,\max
X)\\}$, for each $R\in{\cal R}$, $\omega(p(R))=(\min X,\max X)$. Consider
first $R^{\prime}\in{\cal R}$ such that $C\subseteq L_{(\min X,\max
X)}(R^{\prime})$. Then, by Theorem 1, $g_{(\min X,\max X)}(R^{\prime})=\min X$
and $f(R^{\prime})=\min X$. Since $C\subseteq L_{(\min X,\max X)}(R^{\prime})$
and $C\subseteq D$, we have that $\min X\,P^{\prime}_{i}\,x$ for each $i\in C$
and each $x\in X\setminus\\{\min X\\}$. Consider now $R^{\prime\prime}\in{\cal
R}$ such that $C\not\subseteq L_{(\min X,\max X)}(R^{\prime\prime})$. Then, by
Theorem 1, $g_{(\min X,\max X)}(R^{\prime\prime})=\max X$ and
$f(R^{\prime\prime})=\max X$. Since $C\not\subseteq L_{(\min X,\max
X)}(R^{\prime\prime})$ and $C\subseteq D$, we have that $\max
X\,P^{\prime\prime}_{i}\,x$ for each $i\in C\setminus L_{(\min X,\max
X)}(R^{\prime\prime})$ and each $x\in X\setminus\\{\max X\\}$. Hence, $f$ is
PE.

Step 3: We show that if $A\neq\emptyset$ and $r_{\omega}\neq\\{(\min X,\max
X)\\}$, then $f$ is not PE.

Since $\Omega=\\{\min X,\max X\\}$ and $r_{\omega}\neq\\{(\min X,\max X)\\}$,
we have that $\min X\in r_{\omega}$ and/or $\max X\in r_{\omega}$. Suppose
without loss of generality that $\min X\in r_{\omega}$. Let $x\not\in\Omega$
and consider $R\in{\cal R}$ such that for each $i\in A$, $[\rho(R_{i})=x$ and
$p(R_{i})=\min X]$, and for each $i\in D$, $\delta(R_{i})=\min X$. Since $\min
X\in r_{\omega}$ and $\\{j\in A:p(R_{j})\leq\min X\\}=A$, it follows from
Theorem 1 that $\omega(p(R))=\min X$ and $f(R)=\min X$. However,
$x\,P_{i}\,\min X$ for each $i\in N$. Thus, $f$ is not PE.$\hfill\square$

## References

  * Achuthankutty and Roy (2018) Achuthankutty, G. and S. Roy (2018): “Dictatorship on top-circular domains,” _Theory and Decision_ , 85, 479–493.
  * Alcalde-Unzu and Vorsatz (2018) Alcalde-Unzu, J. and M. Vorsatz (2018): “Strategy-proof location of public facilities,” _Games and Economic Behavior_ , 112, 21–48.
  * Barberà (2011) Barberà, S. (2011): “Strategyproof social choice,” in _Handbook of social choice and welfare_ , Elsevier, vol. 2, 731–831.
  * Barberà et al. (2010) Barberà, S., D. Berga, and B. Moreno (2010): “Individual versus group strategy-proofness: When do they coincide?” _Journal of Economic Theory_ , 145, 1648–1674.
  * Barberà et al. (2012) ——— (2012): “Domains, ranges and strategy-proofness: the case of single-dipped preferences,” _Social Choice and Welfare_ , 39, 335–352.
  * Barberà and Jackson (1994) Barberà, S. and M. Jackson (1994): “A characterization of strategy-proof social choice functions for economies with pure public goods,” _Social Choice and Welfare_ , 11, 241–252.
  * Berga (1998) Berga, D. (1998): “Strategy-proofness and single-plateaued preferences,” _Mathematical Social Sciences_ , 35, 105–120.
  * Berga and Serizawa (2000) Berga, D. and S. Serizawa (2000): “Maximal domain for strategy-proof rules with one public good,” _Journal of Economic Theory_ , 90, 39–61.
  * Black (1948a) Black, D. (1948a): “The decisions of a committee using a special majority,” _Econometrica_ , 245–261.
  * Black (1948b) ——— (1948b): “On the rationale of group decision-making,” _Journal of Political Economy_ , 56, 23–34.
  * Bossert and Peters (2014) Bossert, W. and H. Peters (2014): “Single-basined choice,” _Journal of Mathematical Economics_ , 52, 162–168.
  * Feigenbaum and Sethuraman (2015) Feigenbaum, I. and J. Sethuraman (2015): “Strategyproof mechanisms for one-dimensional hybrid and obnoxious facility location models,” in _Workshops at the twenty-ninth AAAI conference on artificial intelligence_.
  * Gibbard (1973) Gibbard, A. (1973): “Manipulation of voting schemes: a general result,” _Econometrica: journal of the Econometric Society_ , 587–601.
  * Manjunath (2014) Manjunath, V. (2014): “Efficient and strategy-proof social choice when preferences are single-dipped,” _International Journal of Game Theory_ , 43, 579–597.
  * Moulin (1980) Moulin, H. (1980): “On strategy-proofness and single peakedness,” _Public Choice_ , 35, 437–455.
  * Moulin (1984) ——— (1984): “Generalized Condorcet-winners for single peaked and single-plateau preferences,” _Social Choice and Welfare_ , 1, 127–147.
  * Rodríguez-Álvarez (2017) Rodríguez-Álvarez, C. (2017): “On single-peakedness and strategy-proofness: ties between adjacent alternatives,” _Economics Bulletin_ , 37, 1966–1974.
  * Satterthwaite (1975) Satterthwaite, M. A. (1975): “Strategy-proofness and Arrow’s conditions: Existence and correspondence theorems for voting procedures and social welfare functions,” _Journal of Economic Theory_ , 10, 187–217.
  * Thomson (2022) Thomson, W. (2022): “Where should your daughter go to college? An axiomatic analysis,” _Social Choice and Welfare_ , 1–18.